text:
function [U,B_k,V] = lanc_b(A,p,k,reorth)
%LANC_B Lanczos bidiagonalization.
%
% B_k = lanc_b(A,p,k,reorth)
% [U,B_k,V] = lanc_b(A,p,k,reorth)
%
% Performs k steps of the Lanczos bidiagonalization process with
% starting vector p, producing a lower bidiagonal matrix
%           [b_11               ]
%           [b_21 b_22          ]
%     B_k = [     b_32 .        ]
%           [          . b_kk   ]
%           [            b_k+1,k]
% such that
%     A*V = U*B_k ,
% where U and V consist of the left and right Lanczos vectors.
%
% Reorthogonalization is controlled by means of reorth:
%    reorth = 0 : no reorthogonalization,
%    reorth = 1 : reorthogonalization by means of MGS,
%    reorth = 2 : Householder-reorthogonalization.
% No reorthogonalization is assumed if reorth is not specified.

% Reference: G. H. Golub & C. F. Van Loan, "Matrix Computations",
% 3. Ed., Johns Hopkins, 1996.  Section 9.3.4.
% Referred to as "bidiag1" by Paige and Saunders.

% Per Christian Hansen, IMM, April 8, 2001.

% Initialization.
if (k<1), error('Number of steps k must be positive'), end
if (nargin < 4), reorth = 0; end
if (reorth < 0 | reorth > 2), error('Illegal reorth'), end
if (nargout==2), error('Not enough output arguments'), end
[m,n] = size(A);
B_k = sparse(k+1,k);
if (nargout>1 | reorth==1)
  U = zeros(m,k); V = zeros(n,k); UV = 1;
else
  UV = 0;
end
if (reorth==2)
  if (k>=n), error('No. of iterations must satisfy k < n'), end
  HHU = zeros(m,k); HHV = zeros(n,k);
  HHalpha = zeros(1,k); HHbeta = HHalpha;
end

% Prepare for Lanczos iteration.
v = zeros(n,1);
beta = norm(p);
if (beta==0), error('Starting vector must be nonzero'), end
if (reorth==2)
  [beta,HHbeta(1),HHU(:,1)] = gen_hh(p);
end
u = p/beta;
if (UV), U(:,1) = u; end

% Perform Lanczos bidiagonalization with/without reorthogonalization.
for i=1:k

  r = A'*u - beta*v;
  if (reorth==0)
    alpha = norm(r); v = r/alpha;
  elseif (reorth==1)
    for j=1:i-1, r = r - (V(:,j)'*r)*V(:,j); end
    alpha = norm(r); v = r/alpha;
  else
    for j=1:i-1
      r(j:n) = app_hh(r(j:n),HHalpha(j),HHV(j:n,j));
    end
    [alpha,HHalpha(i),HHV(i:n,i)] = gen_hh(r(i:n));
    v = zeros(n,1); v(i) = 1;
    for j=i:-1:1
      v(j:n) = app_hh(v(j:n),HHalpha(j),HHV(j:n,j));
    end
  end
  B_k(i,i) = alpha; if (UV), V(:,i) = v; end

  p = A*v - alpha*u;
  if (reorth==0)
    beta = norm(p); u = p/beta;
  elseif (reorth==1)
    for j=1:i, p = p - (U(:,j)'*p)*U(:,j); end
    beta = norm(p); u = p/beta;
  else
    for j=1:i
      p(j:m) = app_hh(p(j:m),HHbeta(j),HHU(j:m,j));
    end
    [beta,HHbeta(i+1),HHU(i+1:m,i+1)] = gen_hh(p(i+1:m));
    u = zeros(m,1); u(i+1) = 1;
    for j=i+1:-1:1
      u(j:m) = app_hh(u(j:m),HHbeta(j),HHU(j:m,j));
    end
  end
  B_k(i+1,i) = beta; if (UV), U(:,i+1) = u; end

end

if (nargout==1), U = B_k; end

repo: REST


text:
% Description:
%
%     Create a design matrix for group-covariate interactions.
%
% Syntax:
%
%     [ X, terms ] = mGC2X(groups, covariates, a, b)
%
% Inputs:
%
%     groups     - [ N x G ] (int)    - columns of qualitative variables
%     covariates - [ N x C ] (double) - columns of quantitative variables
%     a          - [ 1 x 1 ] (int)    - group offset index for terms
%     b          - [ 1 x 1 ] (int)    - covariate offset index for terms
%     options    - [ 1 x P ] (cell)   - see Options
%     columns    - [ T x 2 ] (int)    - see Details
%
% Outputs:
%
%     X     - [ N x M ] (double)
%     terms - [ 1 x M ] (cell)
%
% Details:
%
%     The first column of the columns input in this case is an index into the
%     columns of groups.  The second column of the columns input is an index
%     into the columns of covariates.  For example, [ 2 1 ] would refer to an
%     interaction between column 2 of groups and column 1 of covariates.
%
% Options:
%
%     'over-determined'  - use over-determined coding for the design matrix
%     'sigma-restricted' - use sigma-restricted coding for the design matrix
%     'verbose'          - display extra information to the command window
%
% Examples:
%
% Notes:
%
% Author(s):
%
%     William Gruner (williamgruner@gmail.com)
%
% References:
%
% Acknowledgements:
%
%     Many thanks to Dr. Erik Erhardt and Dr. Elena Allen of the Mind Research
%     Network (www.mrn.org) for their continued collaboration.
%
% Version:
%
%     $Author: williamgruner $
%     $Date: 2010-04-12 07:14:07 -0600 (Mon, 12 Apr 2010) $
%     $Revision: 494 $

function [ X, terms ] = mGC2X(groups, covariates, a, b, options, columns)

    if ~exist('a', 'var') || isempty(a)
        a = 0;
    end

    if ~exist('b', 'var') || isempty(b)
        b = 0;
    end

    if ~exist('options', 'var')
        options = {};
    end

    if ~exist('columns', 'var')

        columns = [];

        for i = 1 : size(groups, 2)
            for j = 1 : size(covariates, 2)
                columns(end + 1, :) = [ i j ];
            end
        end

    end

    terms = {};
    X     = [];

    [ G, g ] = mG2X(groups, 0, intersect(options, ...
        { 'over-determined' 'sigma-restricted' }));

    if ~isempty(strmatch('verbose', options, 'exact'))
        fprintf('\n')
    end

    for i = 1 : size(columns, 1)

        if ~isempty(strmatch('verbose', options, 'exact'))
            fprintf('Factor %d%d represents the interaction between column %d of groups and column %d of covariates.\n', ...
                columns(i, 1) + a, columns(i, 2) + b, columns(i, 1), columns(i, 2))
        end

        I = mFindTerms(columns(i, 1), g);
        x = G(:, I) .* repmat(covariates(:, columns(i, 2)), 1, length(I));
        X = cat(2, X, x);

        for j = 1 : length(I)
            terms{end + 1} = columns(i, :) + [ a b ];
        end

    end

author: Sable

repo: mcbench-benchmarks


text:
function graph( plain, stego, bit, datasig, window, N, L )

bits = str2num(reshape(bit, N, 1))';

figure('Name','Data Bits','NumberTitle','off')
stairs(bits,'linewidth',1);
axis([0 length(bit) -0.2 1.2]);
grid minor;

figure('Name','Direct Sequence Spread Spectrum','NumberTitle','off')
subplot(2,1,1);
plot(datasig); axis([0 N*L -2 2]); grid minor;
title('Spread spectrum with data bits');

subplot(2,1,2);
plot(window); axis([0 N*L -2 2]); grid minor;
title('Hanning windowed data signal');

figure('Name','Audio Signals','NumberTitle','off')
subplot(2,1,1);
plot(plain(1:N*L,1)); axis([0 N*L -2 2]);
title('Plain signal');

subplot(2,1,2);
plot(stego(1:N*L,1)); axis([0 N*L -2 2]);
title('Stego signal');

t = floor(N/2);
figure('Name',['Segment: ' num2str(t)],'NumberTitle','off')
subplot(2,1,1);
plot(plain(t*L+1:(t+1)*L, 1)); axis([0 L -2 2]);
title('Plain signal');

subplot(2,1,2);
plot(stego(t*L+1:(t+1)*L, 1)); axis([0 L -2 2]);
title('Stego signal');

end

author: ktekeli

repo: audio-steganography-algorithms


text:
function X = kr(U,varargin)
%KR Khatri-Rao product.
%   kr(A,B) returns the Khatri-Rao product of two matrices A and B, of
%   dimensions I-by-K and J-by-K respectively. The result is an I*J-by-K
%   matrix formed by the matching columnwise Kronecker products, i.e.
%   the k-th column of the Khatri-Rao product is defined as
%   kron(A(:,k),B(:,k)).
%
%   kr(A,B,C,...) and kr({A B C ...}) compute a string of Khatri-Rao
%   products A o B o C o ..., where o denotes the Khatri-Rao product.
%
%   See also kron.

%   Version: 21/10/10
%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)

if ~iscell(U), U = [U varargin]; end
K = size(U{1},2);
if any(cellfun('size',U,2)-K)
    error('kr:ColumnMismatch', ...
          'Input matrices must have the same number of columns.');
end
J = size(U{end},1);
X = reshape(U{end},[J 1 K]);
for n = length(U)-1:-1:1
    I = size(U{n},1);
    A = reshape(U{n},[1 I K]);
    X = reshape(bsxfun(@times,A,X),[I*J 1 K]);
    J = I*J;
end
X = reshape(X,[size(X,1) K]);

author: andrewssobral


text:
function [W, Zhat] = ica(X,Nsources,Wprev)
%ICA  Perform independent component analysis
%  [W, ZHAT] = ICA(X) performs independent component analysis on data
%  observation matrix X of physiological observations.
%
%   Inputs:
%       X               = Observation matrix: rows should be observations, cols samples.
%       Nsources        = Number of source signals (optional).
%       Wprev           = Timepoint at which to start process (default = 0 seconds).
%
%   Outputs:
%       W               = Demixing matrix.
%       Zhat            = Source signal matrix.
%
%  This approach uses Cardoso's ICA
%  algorithm to estimate sources (ZHAT) and the de-mixing matrix W, an
%  approximation to A^{-1}, the original (unknown) mixing matrix.
%
% Daniel McDuff, Ethan Blackford, Justin Estepp, December 2018
% Based on an implementation by: G D Clifford (2004) gari AT mit DOT edu
% Licensed under the MIT License and the RAIL AI License.

[nRows, nCols] = size(X);
if nRows > nCols
    fprintf('Warning - The number of rows is cannot be greater than the number of columns.\n');
    error('Please transpose input.');
end

if Nsources > min([nRows nCols])
    Nsources = min([nRows nCols]);
    fprintf('Warning - The number of soures cannot exceed number of observation channels. \n')
    fprintf('The number of sources will be reduced to the number of observation channels (%i) \n',Nsources)
end

if exist('Vprev','var')
    [Winv, Zhat] = jade(X,Nsources,Wprev);
else
    [Winv, Zhat] = jade(X,Nsources);
end
W = pinv(Winv);

author: danmcduff

repo: iphys-toolbox


text:
function fx1 = p11_fx1 ( x )

%*****************************************************************************80
%
%% P11_FX1 evaluates the derivative of the function for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX1, the first derivative of the function at X.
%
  epsilon = 0.00001;

  fx1 = - 4.0 * x.^3 * ( epsilon + 256.0 ) / ( 16.0 * x.^4 + epsilon ).^2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef GemImageDisplayParameters < CoreBaseClass
    % GemImageDisplayParameters. Parameters for visualising an image
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    events
        OpacityChanged
        ShowImageChanged
    end

    properties (SetObservable)
        Window = 1600          % The image window (in HU for CT images)
        Level = -600           % The image level (in HU for CT images)
        Opacity = 100
        ShowImage = true
        BlackIsTransparent = false
        OpaqueColour
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%  converts a basic IRIS file into a RISE one
%
%  iris2rise(irisFileName)
%  iris2rise(irisFileName,riseFileName)
%  iris2rise(irisFileName,riseFileName,stderr_name)
%

author: jmaih


text:
function r8blt_print ( n, ml, a, title )

%*****************************************************************************80
%
%% R8BLT_PRINT prints a R8BLT matrix.
%
%  Discussion:
%
%    The R8BLT storage format is appropriate for a banded lower triangular matrix.
%    The matrix is assumed to be zero below the ML-th subdiagonal.
%    The matrix is stored in an ML+1 by N array, in which the diagonal
%    appears in the first row, followed by successive subdiagonals.
%    Columns are preserved.
%
%  Example:
%
%    N = 5, ML = 2
%
%    A11   0   0   0   0
%    A21 A22   0   0   0
%    A31 A32 A33   0   0
%      0 A42 A43 A44   0
%      0   0 A53 A54 A55
%                --- ---
%                    ---
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ML, the lower bandwidth.
%
%    Input, real A(ML+1,N), the R8BLT matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8blt_print_some ( n, ml, a, 1, 1, n, n, title );

  return
end


text:
function demoDrawPoint3d
%DEMODRAWPOINT3D Demo of drawPoint3d
%
%   Example
%     demoDrawPoint3d
%
%   See also
%

% ------
% Author: oqilipo
% Created: 2018-01-12, using R2017b
% Copyright 2018

%% Initialisation
ptsx = rand(1,5);
ptsy = rand(1,5);
ptsz = rand(1,5);

pts2 = rand(5,3);

%% 1
figure('color','w');
drawPoint3d(ptsx, ptsy, ptsz)
hold on
drawPoint3d(gca,pts2)
%% 2
figure('color','w')
drawPoint3d(pts2,'g.')
hold on
drawPoint3d(gca,ptsx, ptsy, ptsz,'r')
%% 3
figure('color','w')
props.Marker='o';
props.MarkerEdgeColor='r';
props.MarkerFaceColor='g';
props.LineStyle='none';
drawPoint3d(ptsx, ptsy, ptsz,props,'MarkerSize',12)
hold on
drawPoint3d(gca,pts2,props,'MarkerSize',5)
%% 4
figure('color','w')
drawPoint3d(ptsx, ptsy, ptsz,'MarkerSize',8,props,'Marker','s')
hold on
drawPoint3d(gca, pts2,'MarkerSize',12,props,'Marker','h')
%% 5
figure('color','w')
drawPoint3d(ptsx, ptsy, ptsz,props)
hold on
drawPoint3d(gca,pts2,props)

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function pass = test_subsref(pref)
% Test Chebfun3v subsref() command.

if ( nargin == 0)
    pref = chebfunpref;
end
tol = 1000*pref.cheb3Prefs.chebfun3eps;

% Test recursive subsref:
f = chebfun3(@(x,y,z) sin(x.*y.*z));
F = [f; f; f];

G = F(1);
exactCore = G.core;

pass(1) = norm(exactCore(:) - F(1).core(:)) < tol;

% Composition of a CHEBFUN3V (2 components) with three CHEBFUN3.
f1 = chebfun3(@(x,y,z) x);
f2 = chebfun3(@(x,y,z) y);
f3 = chebfun3(@(x,y,z) z);
G = chebfun3v(@(x,y,z) x + y, @(x,y,z) y);
H = G(f1, f2, f3);
pass(2) = ( norm(H - G) < tol );

% Composition of a CHEBFUN3V (3 components) with three CHEBFUN3.
f1 = chebfun3(@(x,y,z) x);
f2 = chebfun3(@(x,y,z) y);
f3 = chebfun3(@(x,y,z) z);
G = chebfun3v(@(x,y,z) x + y, @(x,y,z) y, @(x,y,z) z);
H = G(f1, f2, f3);
pass(3) = ( norm(H - G) < tol );

% Composition of a CHEBFUN3V (2 components) with three CHEBFUN2.
f1 = chebfun2(@(x,y) x);
f2 = chebfun2(@(x,y) y);
f3 = chebfun2(@(x,y) x+y);
G = chebfun3v(@(x,y,z) x + y, @(x,y,z) y + z, [ -1, 1, -1, 1, -2, 2 ]);
H = G(f1, f2, f3);
H_true = chebfun2v(@(x,y) x + y, @(x,y) x + 2*y);
pass(4) = ( norm(H - H_true) < tol );

% Composition of a CHEBFUN3V (3 components) with three CHEBFUN2.
f1 = chebfun2(@(x,y) x);
f2 = chebfun2(@(x,y) y);
f3 = chebfun2(@(x,y) x+y);
G = chebfun3v(@(x,y,z) x + y, @(x,y,z) y + z, @(x,y,z) x + z, ...
    [ -1, 1, -1, 1, -2, 2 ]);
H = G(f1, f2, f3);
H_true = chebfun2v(@(x,y) x + y, @(x,y) x + 2*y, @(x,y) 2*x + y);
pass(5) = ( norm(H - H_true) < tol );

% Composition of a CHEBFUN3V with a CHEBFUN2V.
G = chebfun3v(@(x,y,z) 2*x, @(x,y,z) y - x, @(x,y,z) z, ...
    [ -1, 1, -1, 1, -2, 2 ]);
F = chebfun2v(@(x,y) x, @(x,y) y, @(x,y) x + y);
H = G(F);
H_true = chebfun2v(@(x,y) 2*x, @(x,y) y - x, @(x,y) x + y);
pass(6) = ( norm(H - H_true) < tol );

% Composition of a CHEBFUN3V with a CHEBFUN3V.
G = chebfun3v(@(x,y,z) x + y + z, @(x,y,z) x - y, [ -2, 2, -2, 2, 0, 2 ]);
F = chebfun3v(@(x,y,z) 2*x, @(x,y,z) x + y, @(x,y,z) z + 1);
H = G(F);
H_true = chebfun3v(@(x,y,z) 3*x + y + z + 1, @(x,y,z) x - y);
pass(7) = ( norm(H - H_true) < tol );

% Test composition with one inf by 3 CHEBFUN:
F = chebfun(@(t) [ t, t, t ]);
G = chebfun3v(@(x,y,z) x + y + z, @(x,y,z) x);
H = G(F);
H_true = chebfun(@(t) [ 3*t, t ]);
pass(8) = ( norm(H - H_true) < tol );

% Test composition with three CHEBFUNs:
f = chebfun(@(t) t);
G = chebfun3v(@(x,y,z) x + y + z, @(x,y,z) x);
H = G(f, f, f);
H_true = chebfun(@(t) [ 3*t, t ]);
pass(9) = ( norm(H - H_true) < tol );

% Test composition with a SPHEREFUNV:
f = spherefunv(@(x,y,z) x, @(x,y,z) y, @(x,y,z) z);
g = chebfun3v(@(x,y,z) x, @(x,y,z) 2*y, @(x,y,z) z);
h_true = spherefunv(@(x,y,z) x, @(x,y,z) 2*y, @(x,y,z) z);
h = g(f);
pass(10) = ( norm(h - h_true) < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function uv = compute_flow(this, init, gt)
%
%COMPUTE_FLOW   Compute flow field
%   UV = COMPUTE_FLOW(THIS[, INIT]) computes the flow field UV with
%   algorithm THIS and the optional initialization INIT.
%
%   This is a member function of the class 'classic_nl_optical_flow'.
%
% Authors: Deqing Sun, Department of Computer Science, Brown University
% Contact: dqsun@cs.brown.edu
% $Date: $
% $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  tic

  % Frame size
  sz = [size(this.images, 1), size(this.images, 2)];

  % If we have no initialization argument, initialize with all zeros
  if (nargin < 2)
    uv = zeros([sz, 2]);
  else
    uv = init;
  end

  % Preprocess input (gray) sequences
  if this.texture
      % Perform ROF structure texture decomposition
      images  = structure_texture_decomposition_rof( this.images, 1/8, 100, this.alp);

  elseif this.fc

      % Laplacian in flowfusion
      f = fspecial('gaussian', [5 5], 1.5);
      %f = fspecial('gaussian', [3 3], 1);
      images  = this.images- this.alp*imfilter(this.images, f, 'symmetric');

      % Gaussian pre-filering
      %         f = fspecial('gaussian', [3 3], 0.5);  % Li & Huttenlocher
      %         images  = imfilter(this.images, f, 'symmetric');

      % Sobel edge magnitude
      %         Dy = fspecial('sobel')/8;
      %         Dx = Dy';
      %         Ix = imfilter(this.images, Dx, 'symmetric');
      %         Iy = imfilter(this.images, Dy, 'symmetric');
      %         images = sqrt(Ix.^2 + Iy.^2);

      images  = scale_image(images, 0, 255);
  else
      images  = scale_image(this.images, 0, 255);
  end;

  if this.auto_level
      % Automatic determine pyramid level
      % this.pyramid_levels  =  1 + floor( log(min(size(images, 1), size(images,2))/16) / log(this.pyramid_spacing) );

      % largest size around 16
      N1 = 1 + floor( log(max(size(images, 1), size(images,2))/16)/...
          log(this.pyramid_spacing) );
      % smaller size shouldn't be less than 6
      N2 = 1 + floor( log(min(size(images, 1), size(images,2))/6)/...
          log(this.pyramid_spacing) );

      % satisfies at least one(N2)
      this.pyramid_levels  =  min(N1, N2);
  end;

  % Construct image pyramid, using filter setting in Bruhn et al in "Lucas/Kanade.." (IJCV2005') page 218

  % For gnc stage 1
  factor            = sqrt(2);  % sqrt(3) worse
  smooth_sigma      = sqrt(this.pyramid_spacing)/factor;   % or sqrt(3) recommended by Manuel Werlberger
  f                 = fspecial('gaussian', 2*round(1.5*smooth_sigma) +1, smooth_sigma);
  pyramid_images    = compute_image_pyramid(images, f, this.pyramid_levels, 1/this.pyramid_spacing);

  % For segmentation purpose
  org_pyramid_images = compute_image_pyramid(this.images, f, this.pyramid_levels, 1/this.pyramid_spacing);
  org_color_pyramid_images = compute_image_pyramid(this.color_images, f, this.pyramid_levels, 1/this.pyramid_spacing);

  % For gnc stage 2 to gnc_iters
  smooth_sigma      = sqrt(this.gnc_pyramid_spacing)/factor;
  f                 = fspecial('gaussian', 2*round(1.5*smooth_sigma) +1, smooth_sigma);
  gnc_pyramid_images= compute_image_pyramid(images, f, this.gnc_pyramid_levels, 1/this.gnc_pyramid_spacing);

  % For segmentation purpose
  org_gnc_pyramid_images = compute_image_pyramid(this.images, f, this.gnc_pyramid_levels, 1/this.gnc_pyramid_spacing);
  org_color_gnc_pyramid_images = compute_image_pyramid(this.color_images, f, this.gnc_pyramid_levels, 1/this.gnc_pyramid_spacing);

  for ignc = 1:this.gnc_iters

      if this.display
          disp(['GNC stage: ', num2str(ignc)])
      end

      if ignc == 1
          pyramid_levels = this.pyramid_levels;
          pyramid_spacing = this.pyramid_spacing;

      else
          pyramid_levels = this.gnc_pyramid_levels;
          pyramid_spacing = this.gnc_pyramid_spacing;

      end;

      % Iterate through all pyramid levels starting at the top
      for l = pyramid_levels:-1:1

          if this.display
              disp(['-Pyramid level: ', num2str(l)])
          end

          % Generate copy of algorithm with single pyramid level and the
          % appropriate subsampling
          small = this;

          if ignc == 1
              nsz               = [size(pyramid_images{l}, 1) size(pyramid_images{l}, 2)];
              small.images      = pyramid_images{l};
              small.max_linear  = 1;             % number of linearization performed per warping, 1 OK for quadratic formulation
              im1   = org_pyramid_images{l}(:,:,1);
              small.color_images      = org_color_pyramid_images{l};

          else
              small.images         = gnc_pyramid_images{l};
              nsz   = [size(gnc_pyramid_images{l}, 1) size(gnc_pyramid_images{l}, 2)];
              im1   = org_gnc_pyramid_images{l}(:,:,1);

              small.color_images      = org_color_gnc_pyramid_images{l};
          end;

          % Rescale the flow field
          uv        = resample_flow(uv, nsz);

          small.seg = im1;

         % Adaptively determine half window size for the area term
          small.affine_hsz      = min(4, max(2, ceil(min(nsz)/75)) );

          % Run flow method on subsampled images
          uv = compute_flow_base(small, uv);
      end

      % Update GNC parameters (linearly)
      if this.gnc_iters > 1
          new_alpha  = 1 - ignc / (this.gnc_iters-1);
          this.alpha = min(this.alpha, new_alpha);
          this.alpha = max(0, this.alpha);
      end;

      fprintf('GNC stage %d finished, %3.2f minutes passed \t \t', ignc, toc/60);

      if nargin == 3
          [aae stdae aepe] = flowAngErr(gt(:,:,1), gt(:,:,2), uv(:,:,1), uv(:,:,2), 0);        % ignore 0 boundary pixels
          fprintf('AAE %3.3f STD %3.3f average end point error %3.3f \n', aae, stdae, aepe);
      else
          fprintf('\n');
      end;

  end;

author: kristinbranson

repo: JAABA


text:
function imData=read_file(path_to_file,sframe,num2read,options,im_info)

% Reads uncompressed multipage .tiff, .hdf5, .avi or .raw files
% Usage:  my_data=read_file('path_to_data_file, start frame, num to read);

% INPUTS:
% path_to_file:     location of file to be read
% sframe:           first frame to read (optional, default: 1)
% num2read:         number of frames to read (optional, default: read the whole file)
% options:          options for reading .raw or .bin files
% im_info:          information about the file (if already present)

% OUTPUT:
% imData:           data in array format

% Written by Eftychios A. Pnevmatikakis, Simons Foundation

if nargin<2 || isempty(sframe); sframe = 1; end
if nargin<3 || isempty(num2read); num2read = Inf; end

[~,~,ext] = fileparts(path_to_file);

if strcmpi(ext,'.tiff') || strcmpi(ext,'.tif') || strcmpi(ext,'.btf')
    if ~exist('im_info','var')
        im_info = imfinfo(path_to_file);
    end
    TifLink = Tiff(path_to_file, 'r');
    num2read = min(num2read,length(im_info)-sframe+1);
    imData = zeros(im_info(1).Height,im_info(1).Width,num2read,'like',TifLink.read());
    for i=1:num2read
       TifLink.setDirectory(i+sframe-1);
       imData(:,:,i)=TifLink.read();
    end
    TifLink.close()
    %imData = loadtiff(path_to_file,sframe,num2read);
elseif strcmpi(ext,'.hdf5') || strcmpi(ext,'.h5')
%     info = hdf5info(path_to_file);
%     dims = info.GroupHierarchy.Datasets.Dims;
%     name = info.GroupHierarchy.Datasets.Name;
    info = h5info(path_to_file);
    dims = info.Datasets.Dataspace.Size;
    name = info.Datasets.Name;
    if nargin < 3
        num2read = dims(end)-sframe+1;
    end
    num2read = min(num2read,dims(end)-sframe+1);
%    imData = h5read(path_to_file,name,[ones(1,length(dims)-1),sframe],[dims(1:end-1),num2read]);
    imData = h5read(path_to_file,['/',name],[ones(1,length(dims)-1),sframe],[dims(1:end-1),num2read]);
elseif strcmpi(ext,'.avi')
    v = VideoReader(path_to_file);
    if v.Duration*v.FrameRate < sframe
        imData = [];
    else
        if nargin < 3
            num2read = v.Duration*v.FrameRate-sframe+1;
        end
        num2read = min(num2read,v.Duration*v.FrameRate-sframe+1);
        Y1 = readFrame(v);
        imData = zeros(v.Height,v.Width,num2read,'like',Y1);
        i = 1;
        if sframe == 1
            imData(:,:,i-sframe+1) = Y1;
        end
        while hasFrame(v)  && (i - sframe + 1 < num2read)
            video = readFrame(v);
            i = i + 1;
            if i >= sframe
                imData(:,:,i-sframe+1) = video;
            end
            if i - sframe + 1 >= num2read
                break;
            end
        end
    end
elseif strcmpi(ext,'.raw')
    if nargin < 4 || ~isfield(options,'d1') || isempty(options.d1);
        options.d1 = input('What is the total number of rows? \n');
    end
    if nargin < 4 || ~isfield(options,'d2') || isempty(options.d2);
        options.d2 = input('What is the total number of columns? \n');
    end
    if nargin < 4 || ~isfield(options,'bitsize') || isempty(options.bitsize);
        options.bitsize = 2;
    end
    imData = read_raw_file(path_to_file,sframe,num2read,[options.d1,options.d2],options.bitsize);
else
    error('Unknown file extension. Only .tiff, .avi and .hdf5 files are currently supported');
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function alphaMets = alphabetizeMetFormulas(metFormulas)
%alphabetizeMetFormulas  Reorder metabolite formulas alphabetically.
%
% USAGE:
%
%   alphaMets = alphabetizeMetFormulas(metFormulas);
%
% INPUTS:
%
%   metFormulas     Cell array vector of metabolite atomic formulas.
%                   NOTE: formulas containing characters other than letters
%                         or digits (A-Z,a-z,0-9), or those containing the
%                         phrase "FULLR", will NOT be reordered.
%
% OUTPUTS:
%
%   alphaMets       Cell array vector of metabolite formulas rearranged in
%                   alphabetical order of the element symbols.
%                   Ex: C17H27N5O4SR2 -> C17H27N5O4R2S
%

% Ignore metabolite formulas containing special characters, or "FULLR",
% which appears in the Recon3D model.
ignoreMets = ~cellfun(@isempty,regexp(metFormulas,'[^A-Za-z0-9]','match')) | ...
             contains(metFormulas,'FULLR');

% initialize alphabetized met formulas with original met formulas
alphaMets = metFormulas;

% reorder metabolite formulas
alphaMets(~ignoreMets) = cellfun(@(F) strjoin(sort(regexp(F,'[A-Z][a-z]*\d*','match')),''),alphaMets(~ignoreMets),'UniformOutput',false);

author: SysBioChalmers

repo: Human-GEM


text:
function y = heat2d_fun ( x )

%*****************************************************************************80
%
%% HEAT2D_FUN demonstrates MATLAB's SPMD command for parallel programming.
%
%  Discussion:
%
%    A black and white image X is input.
%
%    A copy of X is distributed across the workers.
%
%    The pixels are regarded as numeric temperatures.
%
%    Each worker applies the explicit heat equation stencil to the pixels
%    in its local part.
%
%    The client combines the results and displays them.
%
%    100 "time" steps are taken, using a CFL coefficient of 0.20.
%    (For this problem, CFL should not be more than 0.25!)
%
%    In this example, we make sure that neighboring workers can share
%    their data, and eliminate the bands that would otherwise show up
%    because of artificial boundaries.
%
%    This is a very artificial calculation, but by using an interesting
%    initial condition, we can make the point that the heat equation "blurs"
%    data, and that by treating the interfaces correctly, we can solve a
%    big problem as a connected set of small problems.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, image X, the initial black and white image.
%
%    Output, image Y, the black and white image after 100 "time" steps of blurring.
%

%
%  CFL is the Courant-Friedrichs-Loewy coefficient, which for this 2D problem
%  should not be larger than 0.25.
%
  cfl = 0.20;
%
%  Since the image is black and white, it is a 2D array.
%  Hence, it will be distributed by columns.
%
  xd = distributed ( x );
%
%  Each worker retrieves its portion of the distributed picture,
%  turn it into a numerical array, and compute the previous and
%  next worker.
%
%  We only need to do these operations once, so do them OUTSIDE the time loop.
%
  spmd

    xl = getLocalPart ( xd );

    xl = double ( xl );

    if ( labindex ~= 1 )
      previous = labindex - 1;
    else
      previous = numlabs;
    end

    if ( labindex ~= numlabs )
      next = labindex + 1;
    else
      next = 1;
    end
%
%  Do 100 steps of blurring.
%
  for i = 1 : 100

    spmd
%
%  Each worker sends its first column to the previous worker,
%  and receives corresponding data from the next worker.
%
      column = labSendReceive ( previous, next, xl(:,1) );

      if ( labindex < numlabs )
        xl = [ xl, column ];
      end
%
%  Each worker sends its last column to the next worker,
%  and receives corresponding data from the previous worker.
%
      column = labSendReceive ( next, previous, xl(:,end) );

      if ( 1 < labindex )
        xl = [ column, xl ];
      end
%
%  Now apply the stencil for the explicit heat equation
%  on the interior points.
%
%  Here, we are modeling:
%
%    dH/dt = - d/dx k ( dH/dx )
%
      [ nr, nc ] = size ( xl );

      c1 = 2:nr-1;
      c2 = 2:nc-1;

      xl(c1,c2) = ( 1 - 4 * cfl ) * xl(c1,c2) ...
        + cfl * ( xl(c1-1,c2) + xl(c1+1,c2) ...
                + xl(c1,c2-1) + xl(c1,c2+1)  );
%
%  West boundary conditions on lab #1.
%
      if ( labindex == 1 )
        xl(:,1) = xl(:,2);
      end
%
%  North and south boundary conditions occur on every lab.
%
      xl(1,:) = xl(2,:);
      xl(end,:) = xl(end-1,:);
%
%  East boundary conditions on last lab.
%
      if ( labindex == numlabs )
        xl(:,end) = xl(:,end-1);
      end
%
%  Now strip off the extra columns.
%
      if ( 1 < labindex )
        xl = xl(:,2:end);
      end

      if ( labindex < numlabs )
        xl = xl(:,1:end-1);
      end
%
%  Convert to unsigned short ints so we can make an image.
%
      yl = uint8 ( xl );

    end
%
%  We are working with a black and white image, so we can simply
%  concatenate the submatrices to get the whole object.
%
    y = [ yl{:} ];
    imshow ( y )
    pause

  end

  return
end


text:
function testFoobar
%testFoobar Passing M-file test

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

author: ricedsp


text:
% ASTROTIK by Francesco Santilli
% R2BP (Restricted Two Bodies Problem)
% M2E inverts the keplerian equation for an elliptical orbit and finds
% eccentric anomaly with a numeric iterative method.
%
% Usage: E = M2E(M,e)
%
% where: M(k) = mean anomaly [rad]
%        e = eccentricity [-] (0<e<1)
%        E(k) = eccentric anomaly [rad]

function E = M2E(M,e)

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

    K = check(M,1);
    check(e,0)

    if (e<=0 || e>=1)
        error('e must be in the range (0,1)')
    end

    E = zeros(K,1);
    for k = 1:K
        F = @(E) E-e*sin(E)-M(k);
        E1 = M(k)-e-eps;
        E2 = M(k)+e+eps;
        E(k) = regulafalsi(F,E1,E2);
    end

end


text:
function [m] = cellmsq(x, dim)

% [M] = CELLMSQ(X, DIM) computes the mean of squares, across all cells in x along
% the dimension dim.
%
% X should be an linear cell-array of matrices for which the size in at
% least one of the dimensions should be the same for all cells

nx = size(x);
if ~iscell(x) || length(nx)>2 || all(nx>1),
  error('incorrect input for cellmsq');
end

if nargin==1,
  scx1 = cellfun('size', x, 1);
  scx2 = cellfun('size', x, 2);
  if     all(scx2==scx2(1)), dim = 2; %let second dimension prevail
  elseif all(scx1==scx1(1)), dim = 1;
  else   error('no dimension to compute mean for');
  end
end

nx   = max(nx);
nsmp = cellfun('size', x, dim);
ssmp = cellfun(@sumsq,   x, repmat({dim},1,nx), 'UniformOutput', 0);
m    = sum(cell2mat(ssmp), dim)./sum(nsmp);

function [s] = sumsq(x, dim)

s = sum(x.^2, dim);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [i1_rect,i2_rect] = rect_pollefeys(F,i1,i2)

% function [i1_rect,i2_rect] = rect_pollefeys(i1,i2,F,P1,P2)
%
%	[i1_rect,i2_rect]=rect_pollefeys(i1,i2,F,P1,P2)
%
%
%Rectifies a pair of images related by fundamental F
%
%IN:
%	i1 - Matlab image
%	i2 - Matlab image
%	F - Fundamental matrix (p1'*F*p2=0). Assumes that image coordiantes
%		are 1..width where pixel centers are at integer locations.
%
%OUT:
%	fig - handle to the figure
[w,h,~] = size(i1);
line  = 0;

% C1 = null(P1);
% e2 = P2*C1;
e1 = null(F);
e2 = null(F');
e2 = e2/e2(3);
% C2 = null(P2);
% e1 = P1*C2;
e1 = e1/e1(3);
thetas = get_theta_bounds(e1,e2,F,[w,h]);
theta = thetas(1);
if 0,
    h0 = my_figure(gr1,gr2,F);
    ud = get(gcf, 'UserData');
end

while theta<thetas(2)
    line = line+1;
    l1 = l_from_theta_p(theta,e1);
    l2 = F_transfer_l(F',l1,e1);
    if 0,
        pts=get_line_points(l1,ud.sizes(1,:));
        axes(ud.ah(1));
        delete(ud.l1);
        delete(ud.l2);
        ud.l1=plot(pts(1,:), pts(2,:), [ud.color '-'], ...
            'LineWidth', ud.size, 'EraseMode','xor');
        pts=get_line_points(l2,ud.sizes(1,:));
        axes(ud.ah(2));
        ud.l2=plot(pts(1,:), pts(2,:), [ud.color '-'], ...
            'LineWidth', ud.size, 'EraseMode','xor');
    end
    [ps1{line},rr1(:,line)] = get_ps(l1,e1,[w,h]);
    [ps2{line},rr2(:,line)] = get_ps(l2,e2,[w,h]);
    dtheta(line) = min(1/rr1(1,line),1/rr2(1,line));
    theta = theta + dtheta(line);
end

%% project images onto new (r,theta) coordtinates
dim = nan(3,1);
dim(1) = size(ps1,2);
rr1_w = max(rr1(2,:))-min(rr1(1,:));
rr2_w = max(rr2(2,:))-min(rr2(1,:));
dim(2) = round(max(rr2_w,rr1_w));
dim(3) = 3;

i1_rect = im_project(i1,ps1,rr1,dim);
i2_rect = im_project(i2,ps2,rr2,dim);

function h0 = my_figure(i1,i2,F)

h0 = figure;
ah1 = axes('Parent', h0, 'Position',[0 0 .5 1]);
h1 = imshow(i1); hold on; title('Image 1');
set(h1, 'ButtonDownFcn','vgg_gui_F(''b1'');');

ah2 = axes('Parent',h0, ...
    'Position',[.5 0 .5 1], ...
    'Tag','Axes2');
h2=imshow(i2); hold on; title('Image 2');
set(h2, 'ButtonDownFcn','vgg_gui_F(''b2'');');

point=plot(-1000, -1000,'EraseMode','xor');
l1=plot([-1000, -1001], [-1000 -1000], 'r-','EraseMode','xor');
l2=plot([-1000, -1001], [-1000 -1000], 'r-','EraseMode','xor');

s1=size(i1); s2=size(i2);
t(:,:,1)=F';  t(:,:,2)=F;  F=t;

ud=struct('h0', h0, 'h',[h1 h2], 'ah', [ah1, ah2], ...
    'sizes', [s1(1:2); s2(1:2)], ...
    'color', 'k', 'size', 1, ...
    'F', F, 'l1', l1,'l2',l2 );

set(h0,'UserData',ud);

function pts=get_line_points(l,sz)
a=l(1); b=l(2);c=l(3);
h=sz(1); w=sz(2);

% This might cause 'divide by zero' warning:
ys=c/-b ;
yf=-(a*w+c)/b;
xs=c/-a;
xf=-(b*h+c)/a;

m1 = [[xs;1] [xf;h] [1;ys] [w;yf]];
w2 = [(xs<=w & xs>=1) (xf<=w & xf>=1) (ys<=h & ys>=1) (yf<=h & yf>=1)];
v = w2>0;
pts = [m1(:,v)];


text:
function model = createEmptyFields(model,fieldNames, fieldDefinitions)
% Create the specified model field with its default values. Works only for
% fields defined in the toolbox if fieldDefinitions are not supplied.
%
% USAGE:
%    model = createEmptyFields(model,fieldName, fieldDefinitions)
%
% INPUTS:
%
%    model:                 The model to add a field to
%    fieldNames:            The names of the fields to add.
%
% OPTIONAL INPUTS:
%    fieldDefinitions:      The Specifications of the field. Only necessary
%                           if a field which is not defined yet should be
%                           added.
%
% OUTPUTS:
%
%    model:                 The original model struct with the specified
%                           field added.
%
% Author:
%    Thomas Pfau Nov 2017

if ischar(fieldNames)
    fieldNames = {fieldNames};
end

if ~exist('fieldDefinitions','var') || isempty(fieldDefinitions)
    fieldDefinitions = getDefinedFieldProperties('SpecificFields',fieldNames);
end

if isempty(fieldDefinitions)
    error('Supplied Field name invalid.')
end
for field = 1:length(fieldNames)
    %Get the dimensions
    xdim = fieldDefinitions{field,2};
    ydim = fieldDefinitions{field,3};

    %And adjust if necessary
    if isnan(xdim)
        xdim = 1;
    end

    %If its listed as a char, it refers to a field (and the first dimension
    %of that field)
    if ischar(xdim)
        if isfield(model,xdim)
            xdim = size(model.(xdim),1);
        else
            xdim = 0;
        end
    end

    if isnan(ydim)
        ydim = 1;
    end

    %If its listed as a char, it refers to a field (and the first dimension
    %of that field)
    if ischar(ydim)
        if isfield(model,ydim)
            ydim = size(model.(ydim),1);
        else
            ydim = 0;
        end
    end

    %Get the field definitions
    fieldType = fieldDefinitions{field,7};
    defaultValue = fieldDefinitions{field,5};
    %define i, for checks relating to model fields which have non ''
    %defaults
    i = 1;
    switch fieldType
        case 'sparse'
            model.(fieldNames{field}) = sparse(xdim,ydim);
        case 'cell'
            model.(fieldNames{field}) =cell(xdim,ydim);
            % need to both check, whether the eval is a char AND if it is empty.
            if xdim > 0 && ischar(eval(defaultValue)) && strcmp(eval(defaultValue),'')
                model.(fieldNames{field})(:) = {''};
            else
                for i = 1:xdim
                    eval(['currentvalue = ' defaultValue ';' ]);
                    model.(fieldNames{field}){i} = currentvalue;
                end
            end
        case 'numeric'
            model.(fieldNames{field}) = repmat(defaultValue,xdim,ydim);
        case 'sparselogical'
            model.(fieldNames{field}) = sparse(repmat(logical(defaultValue),xdim,ydim));
        case 'char'
            model.(fieldNames{field}) = repmat(defaultValue,xdim,ydim);
    end
end

author: opencobra

repo: cobratoolbox


text:
function SNR = VBA_spm_getSNR(y,u,hrf)

% Computes fMRI SNR from controlled time windows
% function SNR = spm_getSNR(y,uu,hrf)
% IN:
%   - y: the nXt1 ROI time series
%   - u: the nuXt2 input (experimental control)
%   - hrf: the hrf (should be same sampling rate than y)
% OUT:
%   - SNR: the SNR (power ratio: signal vs noise)

ru = VBA_spm_resample(full(u),size(y,2)./size(u,2));
nu = size(u,1);
for i=1:nu
    tmp = conv(ru(i,:),hrf);
    tmp = tmp./std(tmp);
    X(:,i) = tmp(2:size(y,2)+1)';
end
iX = pinv(X'*X)*X';
ind = find(abs(mean(ru,1))>1e-1);
nreg = size(y,1);
for i=1:nreg
    beta = iX*y(i,:)';
    yc = X*beta;
    yp(:,i) = yc(ind);
    np(:,i) = y(i,ind)'-yc(ind);
end
SNR = sum(yp(:).^2)./sum(np(:).^2);

author: MBB-team

repo: VBA-toolbox


text:
function ppTrigger(value,raiseTime, delayTime)
% function - ppTrigger
%
% Send a value to the parallel port
%
% Synopsis:
%   ppTrigger(value)
%   ppTrigger(value, raiseTime)
%   ppTrigger(value, raiseTime, delayTime)
%
% Description:
%   The function uses ppWrite to write a value on the parallel port if
%   the global variable IO_ADDR is set. If the global variable
%   NOUZZ_UDP_SOCKET is set the values are also send via udp.
%
%   Currently the delay is not implemented for nouzz.
%
% Input:
%       value: The value you want write on the parallel port
%   raiseTime:[Optional] The time in ms how long the value should stay on
%             the parallel port. (Default: 10ms)
%   delayTime:[Optional] The time in ms how long the value should be
%              delayed. (Default: 0ms)

% Author ??? (Max Sagebaum)
%   2011/11/16 - Max Sagebaum
%                - Added documentation
%                - Added delayTime argument

global IO_ADDR
global NOUZZ_UDP_SOCKET
%global ETH_MARKERS

if ~isempty(IO_ADDR)
    if exist('delayTime','var')
      ppWrite(IO_ADDR, value, raiseTime, delayTime);
    elseif exist('raiseTime','var')
      ppWrite(IO_ADDR, value, raiseTime);
    else
      ppWrite(IO_ADDR, value);
    end
end

if ~isempty(NOUZZ_UDP_SOCKET)
    if exist('delayTime','var')
        if 0 ~= delayTime
            warning('ppTrigger: delay is not implemented for nouzz');
        end
    end
    pnet(NOUZZ_UDP_SOCKET, 'write', [uint8('M') typecast(uint16(value),'uint8')]);
    pnet(NOUZZ_UDP_SOCKET, 'writepacket');
end

%if (~isempty(ETH_MARKERS))
%    fprintf('redirecting trigger eth: %d\n', value);
 %   rawethmex('sendto', 1, sprintf('S%3d',value), 3);
%end

author: PatternRecognition

repo: OpenBMI


text:
function w = K_symmetrised(psi,q1,q2,CS,SS,varargin)
% evaluate kernel modulo symmetries
%
% Input
%  psi    - @SO3Kernel
%  q1, q2 - @quaternion(s)
%  CS, SS - crystal , specimen @symmetry
%
% Options
%  exact -
%  epsilon -
%
% Description
% K(q1,q2) = Sum(S) Sum(l) A_l Tr T_l(s1^-1 q1 s2)

% only the pur rotational part is of interest
% TODO
qCS = unique(quaternion(CS));
qSS = unique(quaternion(SS));

if check_option(varargin,'exact')
  epsilon = pi;
else
  epsilon = min(pi,get_option(varargin,'epsilon',psi.halfwidth*3.5));
end

% how to use sparse matrix representation
if isa(q1,'SO3Grid')
  lg1 = length(q1);
else
  lg1 = -length(q1);
end
if isa(q2,'SO3Grid')
  lg2 = length(q2);
else
  lg2 = -length(q2);
end

if epsilon>2*pi/CS.Laue.multiplicityZ % full matrixes

  q1 = quaternion(q1);
  q2 = quaternion(q2);
  w = zeros(length(q1),length(q2));

	for iks = 1:length(qCS)
		for ips = 1:length(qSS) % for all symmetries

			sg    = qSS(ips) * q1 * qCS(iks);  % rotate g1
      omega = abs(dot_outer(sg,q2));      % calculate full distance matrix
      w = w + psi.eval(omega);

		end
  end

elseif (lg1>0 || lg2>0) && ~check_option(varargin,'old')

  w = sparse(abs(lg1),abs(lg2));

  % sum over specimen symmetry
  if (lg1 >= lg2)              % first argument is SO3Grid
    for issq = 1:length(qSS)
      d = abs(dot_outer(q1,qSS(issq)*quaternion(q2),'epsilon',epsilon,...
        'nospecimensymmetry'));
        w = w + spfun(@psi.eval,d);
    end
  else                         % second argument is SO3Grid
    for issq = 1:length(qSS)
      d = abs(dot_outer(q2,qSS(issq)*quaternion(q1),'epsilon',epsilon,...
        'nospecimensymmetry'));
      w = w + spfun(@psi.eval,d.');
    end
  end

else

	q1 = quaternion(q1);
	q2 = quaternion(q2);

	w = sparse(length(q1),length(q2));

	for iks = 1:length(qCS)
		for ips = 1:length(qSS) % for all symmetries

      if abs(lg1) > abs(lg2)
        sg    = qSS(ips) * q2 * qCS(iks);  % rotate g1
        omega = abs(dot_outer(q1,sg));      % calculate full distance matrix
      else
        sg    = qSS(ips) * q1 * qCS(iks);  % rotate g1
        omega = abs(dot_outer(sg,q2));      % calculate full distance matrix
      end

%  z = find(omega>cos(epsilon));
%  if length(z) > length(omega)/length(CS)/10, w = full(w); end
%  w(z) = w(z) +  kk.K(omega(z));

%  if length(z) > numel(omega)/length(CS)/10, w = full(w); end

      [y,x] = find(omega>cos(epsilon));
      dummy = sparse(y,x,psi.eval(omega(sub2ind(size(w),y,x))),length(q1),length(q2));

      w = w + dummy;

		end
  end

end
%nnz(w)
w = w / length(qCS) / length(qSS);

author: mtex-toolbox

repo: mtex

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


text:
function [y,zf,u,p]=randfilt(pb,pa,ny,zi)
%RANDFILT Generate filtered gaussian noise without initial transient
%
%  Inputs: pb(1,:)  Numerator polynomial of discrete time filter
%          pa(1,:)  Denominator polynomial of discrete time filter
%          ny       Number of output samples required
%          zi       Optional initial filter state
%
% Outputs: y(ny,1)  Filtered Gaussian noise
%          zf       final filter state (can be used to extend the noise sequence)
%          u        The state covariance matrix, <zf*zf'>, is u*u'
%          p        Is the expected value of y(i)^2
%
% zf and zi are the output and optional input state as defined in filter()
% If zi is not specified, random numbers with the correct covariance will be used.
% output u*u' is the state covariance matrix for filter(). Output p is the
% mean power of the output signal y.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: randfilt.m 5885 2015-03-18 09:35:01Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% first normalize the denominator polynomial if necessary

if pa(1)~=1
    pb=pb/pa(1); pa=pa/pa(1);
end

% check to see if we must generate zi

if nargin<4 | nargout>2
    lb=length(pb);
    la=length(pa);

    k=max(la,lb)-1;
    l=la-1;
    ii=k+1-l:k;

    % form controllability matrix

    q=zeros(k,k);
    [z,q(:,1)]=filter(pb,pa,1);
    for i=2:k [z,q(:,i)]=filter(pb,pa,0,q(:,i-1)); end

    % we generate m through the step-down procedure
    s=randn(k,1);
    if l
        m=zeros(l,l);
        g=pa;
        for i=1:l
            g=(g(1)*g(1:end-1)-g(end)*g(end:-1:2))/sqrt((g(1)-g(end))*(g(1)+g(end)));
            m(i,i:l)=g;
        end
        s(ii)=triu(toeplitz(pa(1:l)))*(m\s(ii));
        if nargout>2
            u=q;
            u(:,ii)=q(:,ii)*triu(toeplitz(pa(1:l)))/m;
        end
    else
        if nargout>2
            u=q;
        end
    end
    if nargin < 4
        if k
            zi=q*s;
        else
            zi=[];
        end
    end
end
if nargout>2
    if ~numel(u)
        p=pb(1).^2;
    else
        p=u(1,:)*u(1,:)'+pb(1).^2;
    end
end
if nargin>2 && ny>0
    [y,zf]=filter(pb,pa,randn(ny,1),zi);
else
    zf=zi;
    y=[];
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
classdef ParadigmSpecRCSP < ParadigmDataflowSimplified
    % Advanced paradigm for oscillatory processes via the Spectrally weighted regularized CSP algorithm.
    %
    % This method is a combination of the features of the Spec-CSP [1] and the RCSP [7] methods.
    %
    % The Spec-CSP paradigm [1] is a more advanced variant of CSP, developed for the Berlin
    % Brain-Computer Interface (BBCI); the primary focus was motor imagery BCI, but the algorithm was
    % designed from the outset to be applicable for a wider range of applications. The implementation
    % closely follows the TR [2].
    %
    % The RCSP (Regularized Common Spatial Patterns) method [7] is a recently-developed regularized
    % extension of the CSP method, which combines multiple regularization terms (of which we here
    % implemented covariance shrinkage and Tikhonov regularization); these terms are themselves
    % derived from earlier regularized CSP methods, such as TRCSP (Tikhonov Regularized CSP).
    %
    % In addition to the aforementioned papers this method also use multi-taper spectral estimation
    % instead of the Welch method to regularize the spectral estimates, unlike the Spec-CSP method.
    %
    % The paradigm is applicable to the majority of oscillatory processes, and is the most advanced
    % spatio-spectrally adaptive method that is currently provided in the toolbox. Whenever the exact
    % frequency and location of some (conjectured) oscillatory process is not known exactly, Spec-RCSP
    % can be used, and typically gives better results than CSP with an appropriately unrestricted (e.g.,
    % broad-band) spectral filter. Several other methods exist to adapt the spectrum to a process of
    % interest, among others Common Spatio-Spectral Patterns [3], Common Sparse Spectral Spatial Pattern
    % [4], r^2-based heuristics [5], automated parameter search, and manual selection based on visual
    % inspection. Several of these methods have been shown to give approx. comparable results [2]. An
    % alternative and competitive method, especially when there complex interactions between frequency
    % bands and time periods are to be modeled is the Dual-Augmented Lagrange paradigm
    % (para_dal/para_dal_hf).
    %
    % The method iteratively optimizes spatial and spectral filters in alternation and extracts
    % log-variance features from the resulting (filtered) signal. These features are subsequently
    % processed by a (typically simple) machine learning component, by default LDA. Learning is
    % therefore significantly slower than CSP. An option is to impose custom prior knowledge on the
    % relevant data spectrum, for example by placing a focus on the alpha rhythm, without ruling out
    % other frequencies. Note that there are parameters which constrain the spectrum: one is the
    % frequency prior and the other is the spectral filter that is applied before running the alorithm;
    % both need to be adapted when the considered spectrum shall be extended (e.g. to high-gamma
    % oscillations). Other parameters which are frequently adapted are the time window of interest and
    % the learner component (e.g., logistic regression is a good alternative choice).
    %
    % Some application areas include detection of major brain rhythm modulations (e.g. theta, alpha,
    % beta, gamma), for example related to relaxation/stress, aspects of workload, emotion,
    % sensori-motor imagery, and in general cortical idle oscillations in various modalities.
    %
    % Example: Consider the goal of predicting the emotion felt by a person at a given time. A possible
    % calibration data set for this task would contain a sequence of blocks in each of which the subject
    % is one out of several possible emotions, indicated by events 'e1','e2','e3','e4' covering these
    % blocks at regular rate. The data might for example be induced via guided imagery [6].
    %
    %   calib = io_loadset('data sets/bruce/emotions.eeg')
    %   myapproach = {'SpecRCSP' 'SignalProcessing',{'EpochExtraction',[-2.5 2.5]}};
    %   [loss,model,stats] = bci_train('Data',calib,'Approach',myapproach, 'TargetMarkers',{'e1','e2','e3','e4'});
    %
    %
    % References:
    %  [1] Tomioka, R., Dornhege, G., Aihara, K., and Mueller, K.-R. "An iterative algorithm for spatio-temporal filter optimization."
    %      In Proceedings of the 3rd International Brain-Computer Interface Workshop and Training Course 2006.
    %  [2] Ryota Tomioka, Guido Dornhege, Guido Nolte, Benjamin Blankertz, Kazuyuki Aihara, and Klaus-Robert Mueller
    %      "Spectrally Weighted Common Spatial Pattern Algorithm for Single Trial EEG Classification",
    %      Mathematical Engineering Technical Reports (METR-2006-40), July 2006.
    %  [3] Steven Lemm, Benjamin Blankertz, Gabriel Curio, and Klaus-Robert M?ller.
    %      "Spatio-spectral filters for improving classification of single trial EEG."
    %      IEEE Trans Biomed Eng, 52(9):1541-1548, 2005.
    %  [4] G. Dornhege, B. Blankertz, M. Krauledat, F. Losch, G. Curio, and K.-R. M?ller,
    %      "Combined optimization of spatial and temporal filters for improving brain-computer interfacing,"
    %      IEEE Transactions on Biomedical Engineering, vol. 53, no. 11, pp. 2274?2281, 2006.
    %  [5] Benjamin Blankertz, Ryota Tomioka, Steven Lemm, Motoaki Kawanabe, and Klaus-Robert Mueller.
    %      "Optimizing spatial filters for robust EEG single-trial analysis."
    %      IEEE Signal Process Mag, 25(1):41-56, January 2008
    %  [6] Onton J & Makeig S. "Broadband high-frequency EEG dynamics during emotion imagination."
    %      Frontiers in Human Neuroscience, 2009.
    %  [7] Lotte, F., Guan, G., "Regularizing Common Spatial Patterns to Improve BCI Designs: Unified Theory and New Algorithms."
    %      IEEE Trans Biomed Eng 58, 2 (2011), 355-362.
    %
    % Name:
    %   Spectrally Weighted RCSP
    %
    %                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                           2015-11-05

    methods

        function defaults = preprocessing_defaults(self) %#ok<MANU>
            defaults = {'IIRFilter',{'Frequencies',[0.5 5],'Mode','highpass'}, 'EpochExtraction',[0.5 3.5], 'Resampling',100};
        end

        function model = feature_adapt(self,varargin) %#ok<INUSL>
            args = arg_define(varargin, ...
                arg_norep('signal'), ...
                ... % standard CSP parameters
                arg({'patterns','PatternPairs'},3,uint32([1 1 64 1000]),'Number of CSP patterns (times two).','cat','Feature Extraction'), ...
                ... % multi-taper spectral estimation
                arg({'freqwnd','FrequencyRange','SpectralPrior'},[7 30],[0 Inf],'Frequency range of interest. This is the overall frequency range within which to compute the cross spectrum.'), ...
                arg({'bandwidth','TimeBandwidth'},5,[],'Spectral smoothing. Controls the bias vs. variance of the spectral estimation. Reasonable values are 1 to 3 (1 being fairly noisy, and 3 being fairly smooth but 5x slower)'), ...
                arg({'tapers','Tapers'},[],[],'Number of tapers. Should be an integer smaller than 2*TimeBandwith; default 2*TimeBandwidth-1'), ...
                arg({'subsample_spectrum','SubsampleSpectrum'},1,[1 Inf],'Sub-sample the spectrum. This is the sub-sampling factor.'), ...
                arg({'robust_spectrum','RobustSpectrum'},false,[],'Robust spectral statistics. Can help when some trials are noisy.'), ...
                arg({'normalize_spectrum','NormalizeSpectrum'},false,[],'Normalize spectrum. This corrects for the 1/f characteristic of the spectrum, which is useful when employing the power-based spec-csp prior to low frequencies.'), ...
                ... % RCSP regularization parameters
                arg({'cov_shrinkage','CovarianceShrinkage','covariance_shrinkage'},0,[0 1],'Covariance shrinkage. This is the shrinkage parameter gamma that blends between the empirical covariance (if 0) and the identity matrix (if 1). A good search range is: search(0:0.1:0.9)'), ...
                arg({'tikhonov_reg','TikhonovRegularization','tikhonov_regularization'},0,[0 Inf],'Tikhonov regularization. This parameter regularizes the CSP objective function. A good search range is: search(10.^(-10:-1))'), ...
                arg({'robust_cov','RobustCovariance','robust_covariance'},false,[],'Robust covariance estimation. Can help when some trials are noisy.'), ...
                ... % Spec-CSP specific parameters (rarely used)
                arg({'pp','ParameterP'},0,[-1 1],'Regularization parameter p''. Can be searched over -1:0.5:1.','cat','Feature Extraction','guru',true), ...
                arg({'qp','ParameterQ'},1,[0 4],'Regularization parameter q''. Can be searched over 0:0.5:4.','cat','Feature Extraction','guru',true), ...
                arg({'steps','MaxIterations'},3,uint32([1 3 10 50]),'Number of iterations. A step is spatial optimization, followed by spectral optimization.','cat','Feature Extraction'));

            [signal,n_of,pp,qp,steps] = deal(args.signal,args.patterns,args.pp,args.qp,args.steps);
            if signal.nbchan == 1
                error('Spec-RCSP does intrinsically not support single-channel data (it is a spatial filter).'); end
            if signal.nbchan < args.patterns
                error('Spec-RCSP prefers to work on at least as many channels as you request output patterns. Please reduce the number of pattern pairs.'); end
            p = pp+qp; q = qp;  % re-parameterize the hyper-parameters p' and q' into p and q
            % number of nC=Channels, nS=Samples, nT=trials, nF=frequency bins
            [nC,nS,nT] = size(signal.data); %#ok<ASGLU>
            for c=2:-1:1
                % Extension: compute the per-class multi-taper cross-spectrum across trials
                X{c} = exp_eval_optimized(set_picktrials(signal,'rank',c));
                F{c} = 2*utl_calc_crossspec('Signal',X{c},'FrequencyRange',args.freqwnd,'TimeBandwidth',args.bandwidth, ...
                    'Tapers',args.tapers,'SubsampleSpectrum',args.subsample_spectrum, 'Measure','real',...
                    'FilteredSubsampling',true,'FeatureFilters',false);
                F{c} = permute(F{c},[2 3 1 4]);
                % perform 1/f normalization
                if args.normalize_spectrum
                    F{c} = bsxfun(@times,F{c},reshape((1:size(F{c},3))/size(F{c},3),1,1,[])); end
                % compute the cross-spectrum V as an average over trials
                if args.robust_cov && X{c}.trials > 1
                    V{c} = reshape(geometric_median(reshape(F{c},[],X{c}.trials)')',nC,nC,[]);
                else
                    V{c} = mean(F{c},4);
                end
            end
            nF = size(F{1},3);
            % 1. initialize the filter set alpha and the number of filters J
            J = 1; alpha{J}(1:nF) = 1;
            % 2. for each step
            for step=1:steps
                % 3. for each set of spectral coefficients alpha{j} (j=1,...,J)
                for j=1:J
                    % 4. calculate sensor covariance matrices for each class from alpha{j}
                    for c=2:-1:1
                        Sigma{c} = zeros(nC); %#ok<*AGROW>
                        for b=1:nF
                            Sigma{c} = Sigma{c} + alpha{j}(b)*V{c}(:,:,b); end
                        % Extension: optionally implement shrinkage towards identity
                        if args.cov_shrinkage > 0
                            Sigma{c} = (1-args.cov_shrinkage)*Sigma{c} + args.cov_shrinkage*eye(nC); end
                    end

                    if args.tikhonov_reg == 0
                        % 5. solve the generalized eigenvalue problem Eq. (2) in Spec-CSP
                        [VV,DD] = eig(Sigma{1},Sigma{1}+Sigma{2});
                        iVV = inv(VV)';
                        % and retain n_of top eigenvectors at both ends of the eigenvalue spectrum...
                        W{j} = {VV(:,1:n_of), VV(:,end-n_of+1:end)};
                        P{j} = {iVV(:,1:n_of), iVV(:,end-n_of+1:end)};
                        % as well as the top eigenvalue for each class
                        lambda(j,:) = [DD(1), DD(end)];
                    else
                        [VV,DD,iVV] = deal({});
                        % Extension: use Tikhonov regularization as in TRCSP
                        for c=2:-1:1
                            M{c} = Sigma{c}/(Sigma{3-c} + args.tikhonov_reg*eye(nC));
                            try
                                [VV{c}, DD{c}] = eig(M{c});
                                [DD{c},order] = sort(diag(DD{c}),'descend');
                                VV{c} = VV{c}(:,order);
                                iVV{c} = inv(VV{c});
                                if ~all(isfinite(iVV{c}(:)))
                                    error('Divergence.'); end
                            catch
                                fprintf('ParadigmSpecRCSP: encountered a divergence.\n');
                                % keep going, results like this will be weeded out by a parameter
                                % search
                                VV{c} = randn(size(M{c}));
                                DD{c} = randn(size(M{c}));
                                iVV{c} = randn(size(M{c}));
                            end
                        end
                        % wrap up (note that the roles of patterns and filters are confused in this
                        % formulation, and that we are using the largest EVs from each matrix)
                        W{j} = {iVV{2}(1:n_of,:)',iVV{1}(1:n_of,:)'};
                        P{j} = {VV{2}(:,1:n_of), VV{1}(:,1:n_of)};
                        % ... however, the outer loop finds the minimal first lambda, so we invert it
                        lambda(j,:) = [1/DD{2}(1), DD{1}(1)];
                    end
                end

                % 7. set W{c} from all W{j}{c} such that lambda(j,c) is minimal/maximal over j
                W = {W{argmin(lambda(:,1))}{1}, W{argmax(lambda(:,2))}{2}};
                P = {P{argmin(lambda(:,1))}{1}, P{argmax(lambda(:,2))}{2}};
                % 8. for each projection w in the concatenated [W{1},W{2}]...
                Wcat = [W{1} W{2}]; J = 2*n_of;
                Pcat = [P{1} P{2}];
                for j=1:J
                    w = Wcat(:,j);
                    % 9. calcualate (across trials within each class) mean and variance of the w-projected cross-spectrum components
                    for c=2:-1:1
                        % part of Eq. (3)
                        s{c} = zeros(size(F{c},4),nF);
                        for k=1:nF
                            for t = 1:size(s{c},1)
                                s{c}(t,k) = w'*F{c}(:,:,k,t)*w; end
                        end
                        if args.robust_spectrum
                            mu_s{c} = median(s{c});
                            var_s{c} = median(abs(bsxfun(@minus,s{c},mu_s{c})))*1.4826;
                        else
                            mu_s{c} = mean(s{c});
                            var_s{c} = var(s{c});
                        end
                    end
                    % 10. update alpha{j} according to Eqs. (4) and (5)
                    for c=2:-1:1
                        for k=1:nF
                            % Eq. (4)
                            alpha_opt{c}(k) = max(0, (mu_s{c}(k)-mu_s{3-c}(k)) / (var_s{1}(k) + var_s{2}(k)) );
                            % Eq. (5), with prior from Eq. (6)
                            alpha_tmp{c}(k) = alpha_opt{c}(k).^q * ((mu_s{1}(k) + mu_s{2}(k))/2).^p;
                        end
                    end
                    % ... as the maximum for both classes
                    alpha{j} = max(alpha_tmp{1},alpha_tmp{2});
                    % and normalize alpha{j} so that it sums to unity
                    alpha{j} = alpha{j} / sum(alpha{j});
                end
            end
            alpha = vertcat(alpha{:})';
            model = struct('W',{Wcat},'P',{Pcat},'nF',{nF},'nC',{nC},'alpha',{alpha},'chanlocs',{signal.chanlocs}, ...
                'freqwnd',{args.freqwnd},'bandwidth',{args.bandwidth},'tapers',{args.tapers}, ...
                'subsample_spectrum',{args.subsample_spectrum},'nTrials',[X{1}.trials,X{2}.trials]);
        end

        function features = feature_extract(self,signal,mdl) %#ok<*INUSL>
            spec = 2*utl_calc_crossspec('Signal',signal,'FrequencyRange',mdl.freqwnd,'TimeBandwidth',mdl.bandwidth, ...
                'Tapers',mdl.tapers,'SubsampleSpectrum',mdl.subsample_spectrum, 'Measure','real',...
                'FilteredSubsampling',true,'FeatureFilters',false);
            features = zeros(size(signal.data,3),size(mdl.W,2));
            [W, alpha] = deal(mdl.W, mdl.alpha);
            for f=1:size(W,2)
                C = permute(sum(bsxfun(@times,spec,alpha(:,f)),1),[2,3,4,1]);
                for t=1:size(signal.data,3)
                        features(t,f) = log(W(:,f)'*C(:,:,t)*W(:,f)); end
            end
        end

        function visualize_model(self,varargin) %#ok<*INUSD>
            args = arg_define([0 3],varargin, ...
                arg_norep({'myparent','Parent'},[],[],'Parent figure.'), ...
                arg_norep({'featuremodel','FeatureModel'},[],[],'Feature model. This is the part of the model that describes the feature extraction.'), ...
                arg_norep({'predictivemodel','PredictiveModel'},[],[],'Predictive model. This is the part of the model that describes the predictive mapping.'), ...
                arg({'patterns','PlotPatterns'},true,[],'Plot patterns instead of filters. Whether to plot spatial patterns (forward projections) rather than spatial filters.'), ...
                arg({'paper','PaperFigure'},false,[],'Use paper-style font sizes. Whether to generate a plot with font sizes etc. adjusted for paper.'), ...
                arg_nogui({'nosedir_override','NoseDirectionOverride'},'',{'','+X','+Y','-X','-Y'},'Override nose direction.'));
            arg_toworkspace(args);

            % no parent: create new figure
            if isempty(myparent)
                myparent = figure('Name','Common Spatial Patterns'); end
            % determine nose direction for EEGLAB graphics
            try
                nosedir = args.fmodel.signal.info.chaninfo.nosedir;
            catch
                disp_once('Nose direction for plotting not store in model; assuming +X');
                nosedir = '+X';
            end
            if ~isempty(nosedir_override)
                nosedir = nosedir_override; end
            % number of pairs, and index of pattern per subplot
            np = size(featuremodel.W,2)/2; idxp = [1:np np+(2*np:-1:np+1)]; idxf = [np+(1:np) 2*np+(2*np:-1:np+1)];
            % for each CSP pattern...
            for p=1:np*2
                subplot(4,np,idxp(p),'Parent',myparent);
                if args.patterns
                    plotdata = featuremodel.P(:,p);
                else
                    plotdata = featuremodel.W(:,p);
                end
                topoplot(plotdata,featuremodel.chanlocs,'nosedir',nosedir);
                subplot(4,np,idxf(p),'Parent',myparent);
                alpha = featuremodel.alpha(:,p);
                range = 1:max(find(alpha)); %#ok<MXFND>
                if ~isfield(featuremodel,'freqs')
                    featuremodel.freqs = linspace(featuremodel.freqwnd(1),featuremodel.freqwnd(2),length(range)); end
                pl=plot(featuremodel.freqs(range),featuremodel.alpha(range,p));
                xlim([min(featuremodel.freqs(range)) max(featuremodel.freqs(range))]);
                l1 = xlabel('Frequency in Hz');
                l2 = ylabel('Weight');
                t=title(['Spec-CSP Pattern ' num2str(p)]);
                if args.paper
                    set([gca,t,l1,l2],'FontUnits','normalized');
                    set([gca,t,l1,l2],'FontSize',0.2);
                    set(pl,'LineWidth',2);
                end
            end
            try set(gcf,'Color',[1 1 1]); end
        end

        function layout = dialog_layout_defaults(self)
            layout = {'SignalProcessing.Resampling.SamplingRate', 'SignalProcessing.FIRFilter.Frequencies', ...
                'SignalProcessing.FIRFilter.Type', 'SignalProcessing.EpochExtraction', '', ...
                'Prediction.FeatureExtraction', '', ...
                'Prediction.MachineLearning.Learner'};
        end

        function tf = needs_voting(self)
            tf = true;
        end

    end
end

repo: REST


text:
% Isoparametric Formulation Implementation
% clear memory
clear all

%% Physical Parameters
% numberElements: number of elements
numberElements = 6;

% Initial Area
A = 12.5e-4; % m^2

% Modulus of elasticity
E = 2e11; % N/m^2

% Length of bar
Ltotal = 1.5; % m

%% Exact solutions
xx = 0:0.1:1.5;
displacements_exact = (40000*xx.^3/3-45000)/(E*A);
stress_exact = (40000*xx.^2)/A;

figure(3)

% Displacements
subplot(1,2,1); hold on;
plot(xx,displacements_exact,'-r');

% Stress
subplot(1,2,2); hold on;

plot(xx,stress_exact,'-r');

test = [1 2 4 8];
for numberElements = test

% Build nodes coordinates and middle points,
delta_x = Ltotal/numberElements;
i = 0:numberElements;
x = delta_x*i';x_mid = (x(2:end)-x(1:end-1))/2 + x(1:end-1);

% Length of each element
L = x(2:end)-x(1:end-1);

%% Build Elements
% numberNodes: number of nodes
numberNodes = numberElements + 1;

% generation of coordinates and connectivities
elementNodes = zeros(numberElements,2);
elementNodes(1,:) = [1,2];
for i = 2:numberElements
    elementNodes(i,:) = elementNodes(i-1,:) + 1;
end
nodeCoordinates = x;

% for structure:
   % displacements: displacement vector
   % force : force vector
   % stiffness: stiffness matrix
force=zeros(numberNodes,1);
stiffness=zeros(numberNodes,numberNodes);
% computation of the system stiffness matrix and force vector
for e=1:numberElements;
  % elementDof: element degrees of freedom (Dof)
  elementDof=elementNodes(e,:) ;
  detJacobian=L(e)/2;
  invJacobian=1/detJacobian;
  ngp = 2;
  [w,xi]=gauss1d(ngp);
  xc=0.5*(nodeCoordinates(elementDof(1))+nodeCoordinates(elementDof(2)));
  for ip=1:ngp;
      [shape,naturalDerivatives]=shapeFunctionL2(xi(ip));
      B=naturalDerivatives*invJacobian;
      stiffness(elementDof,elementDof)=...
      stiffness(elementDof,elementDof)+ B'*B*w(ip)*detJacobian*E*A;
      force(elementDof) = force(elementDof)+...
          -80000*(xc+detJacobian*xi(ip))*shape'*detJacobian*w(ip);
  end
end

%% BCs and solution
% prescribed dofs
prescribedDof = find(nodeCoordinates == 1.5); % fixed at node x = 1.5 m
% solution
GDof=numberNodes;
displacements=solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness, ...
    numberNodes,prescribedDof,force)

%% Stress
stress = zeros(numberElements,1);
for e=1:numberElements;
    stress(elementNodes(e,:)) = E*B*displacements(elementNodes(e,:));
end

%% plot figures

switch numberElements
    case 1
        line='r*--';
    case 2
        line='g*--';
    case 4
        line='k*--';
    case 8
        line='m*--';
end

% Displacements
subplot(1,2,1); plot(nodeCoordinates,displacements,line)

% Stress
subplot(1,2,2); stairs(nodeCoordinates,stress,line)

end %end for
subplot(1,2,1); legend('Exact','NEL=1','NEL=2','NEL=4','NEL=8',2); hold off
subplot(1,2,2); legend('Exact','NEL=1','NEL=2','NEL=4','NEL=8',2); hold off

repo: Aero-matlab


text:
%%                   testinitOrthogonal.m

% This test file runs the Orthogonality Promoting initializer. The
% code builds a synthetic formulation of the Phase Retrieval problem, b =
% |Ax| and computes an estimate to x. The code finally outputs the
% correlation achieved.

% PAPER TITLE:
%              Phase retrieval with one or two diffraction patterns by
%              alternating projection with the null initializer.
% ARXIV LINK:
%              https://arxiv.org/pdf/1510.07379.pdf
% PAPER TITLE:
%              Solving Systems of Random Quadratic Equations via Truncated
%              Amplitude Flow.
% ARXIV LINK:
%              https://arxiv.org/pdf/1605.08285.pdf

% 1.) Each test script for an initializer starts out by defining the length
% of the unknown signal, n and the number of measurements, m. These
% mesurements can be complex by setting the isComplex flag to be true.

% 2.) We then build the test problem by generating random gaussian
% measurements and using b0 = abs(Ax) where A is the measurement matrix.

% 3.) We run x0 = initSpectral(A,[],b0,n,isTruncated,isScaled) which runs
% the initialier and and recovers the test vector with high correlation.

% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START----------------------------------

clc
clear
close all

% Parameters
n = 500;           % number of unknowns
m = 8*n;          % number of measurements
isComplex = true;  % use complex matrices? or just stick to real?

%  Build the test problem
xt = randn(n,1)+isComplex*randn(n,1)*1i; % true solution
A = randn(m,n)+isComplex*randn(m,n)*1i;  % matrix
b0 = abs(A*xt);                          % data

% Invoke the truncated spectral initial method
x0 = initOrthogonal(A,[],b0,n);

% Calculate the correlation between the recovered signal and the true signal
correlation = abs(x0'*xt/norm(x0)/norm(xt));

fprintf('correlation: %f\n', correlation);

author: tomgoldstein

repo: phasepack-matlab


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = sigmoid(z).*(1 - sigmoid(z));

% =============================================================

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function out = mkconstarray(class, value, size)
%MKCONSTARRAY creates a constant array of a specified numeric class.
%   A = MKCONSTARRAY(CLASS, VALUE, SIZE) creates a constant array
%   of value VALUE and of size SIZE.

%   Copyright 1993-2003 The MathWorks, Inc.
%   $Revision: 1.8.4.1 $  $Date: 2003/01/26 06:00:35 $

out = repmat(feval(class, value), size);

author: gpeyre

repo: matlab-toolboxes


text:
%% gtitle
% Below is a demonstration of the features of the |gtitle| function

%%
clear; close all; clc;

%% Syntax
% |hf=gtitle(textString,);|

%% Description
% The |gtitle| function can be used to create a title at the top of a
% figure window. The title is tied to the figure rather than any axes and
% can therefore be combined with figures with multiple subplots.

%% Examples

%%
% Settings
titleString='Title on top'; %The title string
fontSize=25; %The desired font size

%%
% Visualize something with a title
hf=cFigure;

hText=gtitle(titleString,fontSize,hf); %Create the title

surf(peaks(25));
colormap(gjet(250));
axisGeom;
camlight headlight;
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/ra_mmv_phase_transition_snr_30db_s_4.mat';
N = 256;
S = 4;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
pta.SNR = 30;

% options for CoSaMP MMV solver
solver_options.RankAwareResidual = true;
P = 2;

% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K, P, solver_options).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
function y_pred = MTrick1(TrainX,TrainY,TestX,alpha,beta,numK,numCircle)

for id = 1:length(TrainY)
    if TrainY(id) == 2 || TrainY(id)==0
        TrainY(id) = -1;
    end
end

G0 = [];
for i = 1:length(TrainY)
    if TrainY(i) == 1
        G0(i,1) = 1;
        G0(i,2) = 0;
    else
        G0(i,1) = 0;
        G0(i,2) = 1;
    end
end

TrainXY = scale_cols(TrainX,TrainY);
fprintf('......start to train logistic regression model.........\n');
w00 = zeros(size(TrainXY,1),1);
lambda = exp(linspace(-0.5,6,20));
wbest = [];
f1max = -inf;
for i = 1:length(lambda)
    w_0 = train_cg(TrainXY,w00,lambda(i));
    f1 = logProb(TrainXY,w_0);
    if f1 > f1max
        f1max = f1;
        wbest = w_0;
        se_lambda = lambda(i);
    end
end
% csvwrite(strcat('model_','test','.model'),wbest);
% wbest = load(strcat('model_','test','.model'));

ptemp = 1./(1 + exp(-wbest'*TrainX));
oriA = getResult(ptemp,TrainY);
ptemp = 1./(1 + exp(-wbest'*TestX));
DataSetX = [TrainX TestX];
% set some variables
Learn.Verbosity = 1;
Learn.Max_Iterations = 200;
Learn.heldout = .1; % for tempered EM only, percentage of held out data
Learn.Min_Likelihood_Change = 1;
Learn.Folding_Iterations = 20; % for TEM only: number of fiolding
% in iterations
Learn.TEM = 0; %tempered or not tempered

[Pw_z,Pz_d,Pd,Li,perp,eta] = pLSA(DataSetX,[],numK,Learn);
pz = Pz_d*Pd';
pw = Pw_z*pz;
A = Pw_z;
% for i = 1:size(Pw_z,1)
%     A(i,:) = A(i,:).*pz';
% end
%
% for i = 1:size(Pw_z,2)
%     for j = 1:length(A(:,i))
%         if pw(j) > 0
%             A(j,i) = A(j,i)./pw(j);
%         else
%             A(j,i) = 1/size(Pw_z,2);
%         end
%     end
% end
pwz = A;
clear A;
% csvwrite(strcat('pwz_common.pwz'),pwz);

% pwz = load(strcat('pwz_common.pwz'));

Fs = pwz;
Ft = Fs;

Gs = G0;
Gt = [];
for i = 1:size(TestX,2)
    Gt(i,1) = ptemp(i);
    Gt(i,2) = 1 - ptemp(i);
end

Xs = TrainX;
Xt = TestX;
Xs = Xs/sum(sum(Xs));
Xt = Xt/sum(sum(Xt));

b = 1/(size(Gs,1));

S = ones(size(Fs,2),size(Gs,2));
for i = 1:size(S,1)
    S(i,:) = S(i,:)/sum(S(i,:));
end

fvalue = trace(Xs'*Xs-2*Xs'*Fs*S*Gs'+Gs*S'*Fs'*Fs*S*Gs')+alpha*b*trace(Gs*Gs'-2*Gs*G0'+G0*G0')+beta*trace(Xt'*Xt-2*Xt'*Ft*S*Gt'+Gt*S'*Ft'*Ft*S*Gt');
tempf = 0;
for circleID = 1:numCircle
    tempM = (Fs*S*Gs'*Gs*S');
    tempM1 = Xs*Gs*S';
    for i = 1:size(Fs,1)
        for j = 1:size(Fs,2)
            if tempM(i,j)~=0
                Fs(i,j) = Fs(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                Fs(i,j) = 0;
            end
        end
    end
    for i = 1:size(Fs,1)
        if sum(Fs(i,:))~= 0
            Fs(i,:) = Fs(i,:)/sum(Fs(i,:));
        else
            for j = 1:size(Fs,2)
                Fs(i,j) = 1/(size(Fs,2));
            end
        end
    end
    tempM = (Gs*S'*(Fs')*Fs*S+alpha.*b*Gs);
    tempM1 = Xs'*Fs*S + alpha.*b*G0;
    for i = 1:size(Gs,1)
        for j = 1:size(Gs,2)
            if tempM(i,j)~=0
                Gs(i,j) = Gs(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                Gs(i,j) = 0;
            end
        end
    end
    for i = 1:size(Gs,1)
        if sum(Gs(i,:))~= 0
            Gs(i,:) = Gs(i,:)/sum(Gs(i,:));
        else
            for j = 1:size(Gs,2)
                Gs(i,j) = 1/(size(Gs,2));
            end
        end
    end

    tempM = (Ft*S*Gt'*Gt*S');
    tempM1 = Xt*Gt*S';
    for i = 1:size(Ft,1)
        for j = 1:size(Ft,2)
            if tempM(i,j)~=0
                Ft(i,j) = Ft(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                Ft(i,j) =0;
            end
        end
    end
    for i = 1:size(Ft,1)
        if sum(Ft(i,:))~= 0
            Ft(i,:) = Ft(i,:)/sum(Ft(i,:));
        else
            for j = 1:size(Ft,2)
                Ft(i,j) = 1/(size(Ft,2));
            end
        end
    end

    tempM = (Gt*S'*Ft'*Ft*S);
    tempM1 = Xt'*Ft*S;
    for i = 1:size(Gt,1)
        for j = 1:size(Gt,2)
            if tempM(i,j)~=0
                Gt(i,j) = Gt(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                Gt(i,j) = 0;
            end
        end
    end
    for i = 1:size(Gt,1)
        if sum(Gt(i,:))~= 0
            Gt(i,:) = Gt(i,:)/sum(Gt(i,:));
        else
            for j = 1:size(Gt,2)
                Gt(i,j) = 1/(size(Gt,2));
            end
        end
    end

    tempM = (Fs'*Fs*S*Gs'*Gs+beta*Ft'*Ft*S*Gt'*Gt);
    tempM1 = Fs'*Xs*Gs+beta*Ft'*Xt*Gt;
    for i = 1:size(S,1)
        for j = 1:size(S,2)
            if tempM(i,j)~=0
                S(i,j) = S(i,j)*(tempM1(i,j)/tempM(i,j))^(0.5);
            else
                S(i,j) = 0;
            end
        end
    end

    fvalue = trace(Xs'*Xs-2*Xs'*Fs*S*Gs'+Gs*S'*Fs'*Fs*S*Gs')+alpha*b*trace(Gs*Gs'-2*Gs*G0'+G0*G0')+beta*trace(Xt'*Xt-2*Xt'*Ft*S*Gt'+Gt*S'*Ft'*Ft*S*Gt');

    pp = [];
    for i = 1:size(TestX, 2)
        if sum(Gt(i,:))~= 0
            pp(1,i) = Gt(i,1)/sum(Gt(i,:));
        else
            pp(1,i) = 0.5;
        end
    end
%     fprintf('the %g iteration is %g,the value of objective is %g\n',circleID,getResult(pp,TestY),fvalue);

    if circleID == 1
        tempf = fvalue;
    end
    if circleID > 1
        if abs(tempf - fvalue) < 10^(-11)
            break;
        end
        tempf = fvalue;
    end
end
for i = 1:size(TestX,2)
    if pp(i) > 0.5
        y_pred(i) = 1;
    else
        y_pred(i) = 0;
    end
end
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function node_display ( v, bi, p )

%% NODE_DISPLAY displays nodes within a boundary.
%
%  Parameters:
%
%    Input, real V(V_NUM,2), the coordinates of vertices used to define the curves.
%
%    Input, integer BI(BI_NUM), a sequence of indices into V.  Each closed curve
%    is defined by giving a sequence of indices, which is terminated by
%    repeating the starting index.  Thus, BI = { 3, 1, 5, 3, 4, 2, 9, 7, 4 }
%    describes two curves: ( 3, 1, 5, 3 ) and (4, 2, 9, 7, 4 ).
%
%    Input,  real P(P_NUM,2), the coordinates of nodes.
%
  bi_num = length ( bi );

  clf
  hold on;
  next = 1;
  s = bi(1);
  t2 = s;
  draw = 1;

  while ( next < bi_num )
    t1 = t2;
    next = next + 1;
    t2 = bi(next);
    if ( draw )
      line ( [ v(t1,1), v(t2,1) ], [ v(t1,2), v(t2,2) ], 'LineWidth', 2, 'Color', 'k' );
      if ( t2 == s )
        draw = 0;
      end
    else
      s = t2;
      draw = 1;
    end

  end

  plot ( p(:,1), p(:,2), 'b.', 'MarkerSize', 15 );

  plot ( v(:,1), v(:,2), 'r.', 'MarkerSize', 15 );

  axis equal
  grid on
  hold off

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ a, info ] = spofa ( a, lda, n )

%*****************************************************************************80
%
%% SPOFA factors a real symmetric positive definite matrix.
%
%  Discussion:
%
%    SPOFA is usually called by SPOCO, but it can be called
%    directly with a saving in time if RCOND is not needed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the symmetric matrix to be  factored.  Only the
%    diagonal and upper triangle are used.
%
%    Input, integer LDA, the leading dimension of the array A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(LDA,N), an upper triangular matrix R so that A = R'*R
%    where R' is the transpose.  The strict lower triangle is unaltered.
%    If INFO /= 0, the factorization is not complete.
%
%    Output, integer INFO, error flag.
%    0, for normal return.
%    K, signals an error condition.  The leading minor of order K is not
%    positive definite.
%
  for j = 1 : n

    s = 0.0;

    for k = 1 : j-1
      t = a(k,j) - sdot ( k-1, a(1:k-1,k), 1, a(1:k-1,j), 1 );
      t = t / a(k,k);
      a(k,j) = t;
      s = s + t * t;
    end

    s = a(j,j) - s;

    if ( s <= 0.0 )
      info = j;
      return
    end

    a(j,j) = sqrt ( s );

  end

  info = 0;

  return
end


text:
function [ m, n, nst, ist, jst, ast ] = st_transpose ( m, n, nst, ist, jst, ...
  ast )

%*****************************************************************************80
%
%% ST_TRANSPOSE transposes an ST matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer M, the number of row.
%
%    Input/output, integer N, the number of columns.
%
%    Input/output, integer NST, the number of nonzeros.
%
%    Input/output, integer IST(NST), JST(NST), the row and column indices.
%
%    Input/output, real AST(NST), the nonzero values.
%
  t = m;
  m = n;
  n = t;

  t          = ist(1:nst);
  ist(1:nst) = jst(1:nst);
  jst(1:nst) = t(1:nst);

  return
end


text:
function uu = VBA_getU(u,options,dim,flag)
% reshape micro-resolution driving input u
% function uu = VBA_getU(u,options,dim,flag)
% IN:
%   - u: macro-/micro- resolution input
%   - options: the options structure (see VBA_check.m)
%   - dim: the dim structure (see VBA_check.m)
%   - flag: gives the transformation direction (from micro- to macro- or
%   back...)
% OUT:
%   - uu: reshaped micro-resolution driving input

if options.microU && ~isequal(options.decim,1)
    switch flag
        case '2macro'
            u0 = u;
            nut = dim.u*options.decim;
            uu = zeros(nut,dim.n_t);
            for t=1:dim.n_t
                uu(:,t) = VBA_vec(u0(:,(t-1)*options.decim+1:t*options.decim));
            end
        case 'back2micro'
            u0 = u;
            uu = zeros(dim.u,options.decim*dim.n_t);
            for t=1:dim.n_t
                ut = reshape(u0(:,t),dim.u,options.decim);
                uu(:,(t-1)*options.decim+1:t*options.decim) = ut;
            end
    end
else % do not change input
    uu = u;
end


text:
function [elem,HB] = uniformcoarsen3red(elem)
%% UNIFORMCOARSEN3RED uniform coarsening of red refinement
%
% [elem,HB] = uniformcoarsen3red(elem) remove grid points added by uniform
% refinement in 3-D.
%
% It is mainly used to get multilevel decomposition in multigrid methods.
% The input matrix elem stands for the fine mesh and the output one for the
% coarse mesh. The HB records the hierarchical structure of added points
% going from the coarse to the fine mesh such that HB(:,2:3) are two parent
% nodes of HB(:,1).
%
%   See also: uniformrefine3, mg, uniformcoarsen, uniformcoarsen3
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

HB = [];
NT = size(elem,1);
if mod(NT,8)==0
    NTc = NT/8; % number of triangles in the coarse grid
else
%     display('Not from red refinement');
    return
end

%% Find points
t1 = 1:NTc; t2 = t1+NTc; t3 = t2+NTc; t4 = t3+NTc;
if any(elem(t1,2)~=elem(t2,1)) || any(elem(t1,3)~=elem(t3,1)) || ...
   any(elem(t4,1)~=elem(t1,4)) || any(elem(t4,2)~=elem(t2,4))
%     display('Not from red refinement');
    return
end
p1 = elem(t1,1);
p2 = elem(t2,2);
p3 = elem(t3,3);
p4 = elem(t4,4);
p5 = elem(t1,2);
p6 = elem(t1,3);
p7 = elem(t1,4);
p8 = elem(t2,3);
p9 = elem(t2,4);
p10 = elem(t3,4);

%% Remove tetrahedron
elem(t1,:) = [p1 p2 p3 p4];
elem = elem(t1,:);

%% Record HB
HB(p5,:) = [p5 p1 p2];
HB(p6,:) = [p6 p1 p3];
HB(p7,:) = [p7 p1 p4];
HB(p8,:) = [p8 p2 p3];
HB(p9,:) = [p9 p2 p4];
HB(p10,:) = [p10 p3 p4];
Nc = max(elem(:));
HB = HB(Nc+1:end,:);


text:
function [c] = cellcov(x, y, dim, flag)

% [C] = CELLCOV(X, DIM) computes the covariance, across all cells in x along
% the dimension dim. When there are three inputs, covariance is computed between
% all cells in x and y
%
% X (and Y) should be linear cell-array(s) of matrices for which the size in at
% least one of the dimensions should be the same for all cells

if nargin<4 && iscell(y)
  flag = 1;
elseif nargin<4 && isnumeric(y)
  flag = dim;
end

if nargin<3 && iscell(y)
  scx1 = cellfun('size', x, 1);
  scx2 = cellfun('size', x, 2);
  if     all(scx2==scx2(1)), dim = 2; %let second dimension prevail
  elseif all(scx1==scx1(1)), dim = 1;
  else   error('no dimension to compute covariance for');
  end
elseif nargin<=3 && isnumeric(y)
  dim = y;
end

if isnumeric(y), y = []; end

nx = size(x);
if ~iscell(x) || length(nx)>2 || all(nx>1),
  error('incorrect input for cellcov');
end

if flag,
  mx   = cellmean(x, 2);
  x    = cellvecadd(x, -mx);
  if ~isempty(y),
    my = cellmean(y, 2);
    y  = cellvecadd(y, -my);
  end
end

nx   = max(nx);
nsmp = cellfun('size', x, dim);
if isempty(y),
  csmp = cellfun(@covc, x, repmat({dim},1,nx), 'UniformOutput', 0);
else
  csmp = cellfun(@covc, x, y, repmat({dim},1,nx), 'UniformOutput', 0);
end
nc   = size(csmp{1});
c    = sum(reshape(cell2mat(csmp), [nc(1) nc(2) nx]), 3)./sum(nsmp);

function [c] = covc(x, y, dim)

if nargin==2,
  dim = y;
  y   = x;
end

if dim==1,
  c = x'*y;
elseif dim==2,
  c = x*y';
end


text:
function x = beta_binomial_cdf_inv ( cdf, a, b, c )

%*****************************************************************************80
%
%% BETA_BINOMIAL_CDF_INV inverts the Beta Binomial CDF.
%
%  Discussion:
%
%    A simple discrete approach is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%
%    Input, real A, B, parameters of the PDF.
%    0.0D+00 < A,
%    0.0D+00 < B.
%
%    Input, integer C, a parameter of the PDF.
%    0 <= C.
%
%    Output, integer X, the smallest X whose cumulative density function
%    is greater than or equal to CDF.
%
  if ( cdf <= 0.0 )

    x = 0;

  else

    cum = 0.0;

    for y = 0 : c

      pdf = beta ( a + y, b + c - y ) / ( ( c + 1 ) ...
        * beta ( y + 1, c - y + 1 ) * beta ( a, b ) );

      cum = cum + pdf;

      if ( cdf <= cum )
        x = y;
        return
      end

    end

    x = c;

  end

  return
end


text:
function calpak_test74 ( )

%*****************************************************************************80
%
%% CALPAK_TEST74 tests YMDF_INC_JULIAN, YMDF_NEXT_JULIAN, YMDF_PREV_JULIAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2013
%
%  Author:
%
%    John Burkardt
%
  days = 10.25;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST74\n' );
  fprintf ( 1, '  For the Julian calendar:\n' );
  fprintf ( 1, '  YMDF_INC_JULIAN increments a date by days;\n' );
  fprintf ( 1, '  YMDF_NEXT_JULIAN computes the next day,\n' );
  fprintf ( 1, '  YMDF_PREV_JULIAN computes the previous day.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  YMDF date    Tomorrow	    Yesterday	   +10.25 days\n' );
  fprintf ( 1, '\n' );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed = jed_test ( i );

    if ( jed < 0 )
      break
    end

    [ y1, m1, d1, f1 ] = jed_to_ymdf_julian ( jed );
    s1 = ymdf_to_s_julian ( y1, m1, d1, f1 );

    [ y2, m2, d2, f2 ] = ymdf_next_julian ( y1, m1, d1, f1 );
    s2 = ymdf_to_s_julian ( y2, m2, d2, f2 );

    [ y3, m3, d3, f3 ] = ymdf_prev_julian ( y1, m1, d1, f1 );
    s3 = ymdf_to_s_julian ( y3, m3, d3, f3 );

    [ y4, m4, d4, f4 ] = ymdf_inc_julian ( y1, m1, d1, f1, days );
    s4 = ymdf_to_s_julian ( y4, m4, d4, f4 );

    fprintf ( 1, '  %15s  %15s  %15s  %15s\n', s1, s2, s3, s4 );

  end

  return
end


text:
classdef MimZeroIsTransparent < MimGuiPlugin
    % MimZeroIsTransparent. Gui Plugin for switching transparency of zero
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimZeroIsTransparent is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Show zero transparent'
        SelectedText = 'Show zero as black'
        ToolTip = 'Changes whether zero values in the overlay are shown as black or transparent'
        Category = 'Segmentation display'
        Visibility = 'Dataset'
        Mode = 'Segment'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'show_black.png'
        Location = 24
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.BlackIsTransparent = ~gui_app.ImagePanel.BlackIsTransparent;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ~strcmp(gui_app.ImagePanel.Mode, MimModes.View3DMode) && gui_app.ImagePanel.ShowOverlay;
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.BlackIsTransparent;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit