text:
function image_rgb(M)
% Show a matrix of integers as a color image.
% This is like imagesc, except we know what the mapping is from integer to color.
% If entries of M contain integers in {1,2,3}, we map
% this to red/green/blue

cmap = [1 0 0; % red
	0 1 0; % green
	0 0 1; % blue
	127/255 1 212/255]; % aquamarine
image(M)
set(gcf,'colormap', cmap);

if 1
  % make dummy handles, one per object type, for the legend
  str = {};
  for i=1:size(cmap,1)
    dummy_handle(i) = line([0 0.1], [0 0.1]);
    set(dummy_handle(i), 'color', cmap(i,:));
    set(dummy_handle(i), 'linewidth', 2);
    str{i} = num2str(i);
  end
  legend(dummy_handle, str, -1);
end

if 0
[nrows ncols] = size(M);
img = zeros(nrows, ncols, 3);
for r=1:nrows
  for c=1:ncols
    q = M(r,c);
    img(r,c,q) = 1;
  end
end
image(img)
end

author: fieldtrip

repo: fieldtrip


text:
function kern = svargplvmInitDynKernel(kern, globalOpt, optionsDyn, indexInComp)

% SVARGPLVMINITDYNKERNEL
% VARGPLVM

X = optionsDyn.t;

fprintf('# Initialising kernel %s', kern.type);
if nargin > 3
    fprintf(' in place %d of the comp structure.\n', indexInComp);
else
    fprintf('\n');
end

% This is a recursive function. The recursion ending condition is for a
% kernel to not be a compound.
if isfield(kern, 'comp')
    for i=1:length(kern.comp)
        if ~isfield(kern.comp{i}, 'index') || isempty(kern.comp{i}.index)
            curInds = 1:kern.comp{i}.inputDimension;
        else
            curInds = kern.comp{i}.index;
        end
        curOptionsDyn = optionsDyn;
        curOptionsDyn.t = X(:, curInds);
        kern.comp{i} = svargplvmInitDynKernel(kern.comp{i}, globalOpt, curOptionsDyn, i);
    end
else
    % This is not a compound kernel, so we can initialise this signel
    % element.
    kernelType = kern.type;
   % if strcmp(kernelType, 'white')
   %     kern.variance = 1e-2; % Usual values: 1e-1, 1e-3
   % end

    if strcmp(kernelType, 'whitefixed')
        kern.variance = globalOpt.fixedwhiteVar;
        fprintf(1,'# fixedwhite variance: %d\n',globalOpt.fixedwhiteVar);
    end

    if strcmp(kernelType, 'bias')
        kern.variance = 0.1;
        fprintf(1,'# Bias variance: %d\n',kern.variance);
    end

    if strcmp(kernelType, 'rbfperiodic') || strcmp(kernelType, 'rbfperiodic2')
        kern.period = globalOpt.periodicPeriod;
        fprintf(1,'# periodic period %d\n',globalOpt.periodicPeriod);
         dst = dist2(optionsDyn.t, optionsDyn.t);
         lb = min(dst(:));
         ub = max(dst(:));
         inv_width = 2 / (ub+lb);
         inv_width = globalOpt.inverseWidthMult / (ub+lb); %% NEW
         kern.inverseWidth = inv_width;
         if strcmp(kernelType, 'rbfperiodic2')
             kern.factor = 2*pi/kern.period;
         end
    end

    % The following is related to the expected number of
    % zero-crossings.(larger inv.width numerator, rougher func)
    if strcmp(kernelType,'rbf') || strcmp(kernelType,'matern32') || strcmp(kernelType,'matern52')
        % NEW
         dst = dist2(optionsDyn.t, optionsDyn.t);
         lb = min(dst(:));
         ub = max(dst(:));
         %inv_width = 2 / (ub+lb);
         inv_width = globalOpt.inverseWidthMult / (ub+lb); %% NEW
        if strcmp(kernelType,'rbf')
            kern.inverseWidth = inv_width;
        elseif strcmp(kernelType, 'matern32') || strcmp(kernelType, 'matern52')
            kern.lengthScale = 1/inv_width;
        end

         %kern.comp{1}.inverseWidth = inv_width;
        %kern.inverseWidth = optionsDyn.inverseWidth./(((max(max(optionsDyn.t))-min(min(optionsDyn.t)))).^2);

        if ~globalOpt.mappingInitialisation && (exist('indexInComp') && indexInComp == 1)
            kern.variance = 1;
        end
    end
end

author: SheffieldML

repo: vargplvm


text:
function d = det3x3(x)

% DET3X3 computes determinant of matrix x, using explicit analytic definition
% if size(x) = [3 3 K M]

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

siz = size(x);
if all(siz(1:2)==3)
  d = x(1,1,:,:).*x(2,2,:,:).*x(3,3,:,:) - ...
      x(1,1,:,:).*x(2,3,:,:).*x(3,2,:,:) - ...
      x(1,2,:,:).*x(2,1,:,:).*x(3,3,:,:) + ...
      x(1,2,:,:).*x(2,3,:,:).*x(3,1,:,:) + ...
      x(1,3,:,:).*x(2,1,:,:).*x(3,2,:,:) - ...
      x(1,3,:,:).*x(2,2,:,:).*x(3,1,:,:);
else
  ft_error('not implemented');
  % write for loop for the higher dimensions, using normal inv
end

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Xtensor = tucker2multiarray(X)
% Converts a 3d Tucker form tensor to a multiarray.
%
% function Xtensor = tucker2multiarray(X)
%
% X has fields U1, U2, U3, and G.
%
% The matrices U1 (n1-by-r1), U2 (n2-by-r2) and U3 (n3-by-r3) are
% orthogonal matrices.
% G (r1-by-r2-by-r3) is a multidimensional array.
%
% See also: fixedrankfactory_tucker_preconditioned

% This file is part of Manopt: www.manopt.org.
% Original authors: Hiroyuki Kasai and Bamdev Mishra, June 05, 2015.
% Contributors:
% Change log:

    U1 = X.U1;
    U2 = X.U2;
    U3 = X.U3;
    G = X.G;

    % Tensor size
    n1 = size(U1, 1);
    n2 = size(U2, 1);
    n3 = size(U3, 1);

    % Core size
    [r1, r2, r3] = size(G);

    % Multplication by U1
    G1 = reshape(G, r1, r2*r3);
    GU1 = reshape(U1*G1, n1, r2, r3);

    % Further multplication by U2
    G2 = reshape(permute(GU1, [2 1 3]), r2, n1*r3);
    GU1U2 = permute(reshape(U2*G2, n2, n1, r3), [2 1 3]);

    % Further multplication by U3
    G3 = reshape(permute(GU1U2, [3 1 2]), r3, n1*n2);
    GU1U2U3 = permute(reshape(U3*G3, n3, n1, n2), [2 3 1]);

    Xtensor = GU1U2U3;% Full tensor

end

author: NicolasBoumal

repo: manopt


text:
function test_ft_plot_patch

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_plot_patch

hdat = [1:10 10:-1:1];
vdat = rand(1,10);
vdat = [vdat vdat(end:-1:1)+1];
ft_plot_patch(hdat, vdat)


text:
function autoContour40SUVMax(structNum)
% autoContour40SUVMax
%
% This function creates a contour at 40% threshold level of Max scan value.
% structNum is the structure index in planC.
%
% Created DK 11/21/06
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS

indexS = planC{end};

scanSet = getStructureAssociatedScan(structNum);

[struct3D slices] = getStruct3DMaskedData(structNum);

siz = size(struct3D);

newStructNum = length(planC{indexS.structures})+ 1;
C = max(struct3D(:));

for i = 1:length(slices)
    slice = struct3D(:,:,i);

    [c I] = max(slice(:));

    if c > (C*4/10)
        [x y]=ind2sub([siz(1) siz(2)],I);

        BW = roicolor(slice,C*4/10,C);

        L = bwlabel(BW, 4);

        region = L(x,y);

        ROI = L == region;

        if isempty(find(ROI))
            continue
        end

        [contour, sliceValues] = maskToPoly(ROI, slices(i), scanSet,planC);

        if(length(contour.segments) > 1)
            longestDist = 0;
            longestSeg =  [];

            for j = 1:length(contour.segments)
                segmentV = contour.segments(j).points(:,1:2);
                curveLength = 0;

                for k = 1:size(segmentV,1) - 1
                    curveLength = curveLength + sepsq(segmentV(k,:)', segmentV(k+1,:)');
                end

                if curveLength > longestDist
                    longestDist = curveLength;
                    longestSeg = j;
                end

            end
            tmp = contour.segments(longestSeg).points;

        else
            if isempty(contour.segments.points)
                continue;
            else
                tmp = contour.segments.points;
            end

        end

        planC{indexS.structures}(newStructNum).contour(slices(i)).segments.points = tmp;
    end

end

for l = slices(end)+1 : length(planC{indexS.scan}(scanSet).scanInfo)
    planC{indexS.structures}(newStructNum).contour(l).segments.points = [];
end
stateS.structsChanged = 1;

planC{indexS.structures}(newStructNum).strUID = createUID('structure');
planC{indexS.structures}(newStructNum).assocScanUID = planC{indexS.structures}(structNum).assocScanUID;
planC{indexS.structures}(newStructNum).structureName = [planC{indexS.structures}(structNum).structureName '_40_SUVMax'];

planC = getRasterSegs(planC, newStructNum);

planC = setUniformizedData(planC);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% funSimLogNormProbCov returns SINR-based k-coverage probability
% under log-normal shadowing based on repeated simulations of
% of model outlined in [1]
%
% simPCovk=funSimLogNormProbCov(tValues,betaConst,K,lambda,sigma,W,diskRadius,simNumb,k)
% simPCovk is the k-coverage probability
% tValues are the SINR threshold values. tValues can be a vector
% betaConst is the pathloss exponent.
% K = path-loss constant
% lambda = density of base station/nodes of cellular network
% sigma = standard deviation of log-normal shadowing (not in dB)
% W = noise constant
% diskRadius = radius of simulation region (Warning: if too small, edge
% effects will not be sufficiently included and disagreement with analytic
% results may occurr)
% simNumb = number of simulations
% k = coverage number (often set to one)
% All input values (except tValues) are scalars
%
% Author: H.P. Keeler, Inria Paris/ENS, 2013
%
% References
% [1] H.P. Keeler, B. Błaszczyszyn and M. Karray,
% 'SINR-based k-coverage probability in cellular networks with arbitrary
% shadowing', accepted at ISIT, 2013

function simPCovk=funSimLogNormProbCov(tValues,betaConst,K,lambda,sigma,W,diskRadius,simNumb,k)

if nargin==8
    k=1;
end

tNumb=length(tValues);

%%% Simulation Section %%%
%(uniformly) randomly places nodes on a disk of radius diskRadius
diskArea=pi*diskRadius^2;
coveredNumbk=zeros(size(tValues));
ESTwoBeta=exp(sigma^2*(2-betaConst)/betaConst^2);
%rescale lambda - see foot note 5 in [1]
lambdaSim=lambda*ESTwoBeta;
for i=1:simNumb
    randNumb=poissrnd(lambdaSim*diskArea);
    %shadowing distribution can be constant if lambda is rescaled - see [1]
    shadowRand=ones(randNumb,1);
    %random distances from the typical node
    rRand=diskRadius*sqrt(rand(randNumb,1)); %uniform in cartesion, not polar coordinates

    signalRand=shadowRand.*(K*rRand).^(-betaConst);
    interferTotal=sum(signalRand); %total inteference in network
    SINR=signalRand./((interferTotal-signalRand)+W); %calculate SINR for each node in the network

    for j=1:tNumb
        T=tValues(j);
        %counts how many nodes are exactly k or more connected/covered
        if sum(SINR>=T)>=k
            coveredNumbk(j)=coveredNumbk(j)+1;
        end
    end

end

simPCovk=coveredNumbk/simNumb;

author: Sable

repo: mcbench-benchmarks


text:
function [ row, col, a ] = r8sp_dif2 ( m, n, nz_num )

%*****************************************************************************80
%
%% R8SP_DIF2 returns the DIF2 matrix in R8SP format.
%
%  Example:
%
%    N = 5
%
%    2 -1  .  .  .
%   -1  2 -1  .  .
%    . -1  2 -1  .
%    .  . -1  2 -1
%    .  .  . -1  2
%
%  Properties:
%
%    A is banded, with bandwidth 3.
%
%    A is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is a special case of the TRIS or tridiagonal scalar matrix.
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A is positive definite.
%
%    A is an M matrix.
%
%    A is weakly diagonally dominant, but not strictly diagonally dominant.
%
%    A has an LU factorization A = L * U, without pivoting.
%
%      The matrix L is lower bidiagonal with subdiagonal elements:
%
%        L(I+1,I) = -I/(I+1)
%
%      The matrix U is upper bidiagonal, with diagonal elements
%
%        U(I,I) = (I+1)/I
%
%      and superdiagonal elements which are all -1.
%
%    A has a Cholesky factorization A = L * L', with L lower bidiagonal.
%
%      L(I,I) =    sqrt ( (I+1) / I )
%      L(I,I-1) = -sqrt ( (I-1) / I )
%
%    The eigenvalues are
%
%      LAMBDA(I) = 2 + 2 * COS(I*PI/(N+1))
%                = 4 SIN^2(I*PI/(2*N+2))
%
%    The corresponding eigenvector X(I) has entries
%
%       X(I)(J) = sqrt(2/(N+1)) * sin ( I*J*PI/(N+1) ).
%
%    Simple linear systems:
%
%      x = (1,1,1,...,1,1),   A*x=(1,0,0,...,0,1)
%
%      x = (1,2,3,...,n-1,n), A*x=(0,0,0,...,0,n+1)
%
%    det ( A ) = N + 1.
%
%    The value of the determinant can be seen by induction,
%    and expanding the determinant across the first row:
%
%      det ( A(N) ) = 2 * det ( A(N-1) ) - (-1) * (-1) * det ( A(N-2) )
%                = 2 * N - (N-1)
%                = N + 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2000
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969,
%    ISBN: 0882756494,
%    LC: QA263.68
%
%    Morris Newman, John Todd,
%    Example A8,
%    The evaluation of matrix inversion programs,
%    Journal of the Society for Industrial and Applied Mathematics,
%    Volume 6, Number 4, pages 466-476, 1958.
%
%    John Todd,
%    Basic Numerical Mathematics,
%    Volume 2: Numerical Algebra,
%    Birkhauser, 1980,
%    ISBN: 0817608117,
%    LC: QA297.T58.
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer NZ_NUM, the number of nonzero elements in
%    the matrix.
%
%    Output, integer ROW(NZ_NUM), COL(NZ_NUM), the row and
%    column indices of the nonzero elements.
%
%    Output, real A(NZ_NUM), the nonzero elements of the matrix.
%
  row = zeros(nz_num,1);
  col = zeros(nz_num,1);
  a = zeros(nz_num,1);

  k = 0;
  for i = 1 : m

    if ( 0 < i - 1 )
      k = k + 1;
      row(k) = i;
      col(k) = i - 1;
      a(k) = -1.0;
    end

    k = k + 1;
    row(k) = i;
    col(k) = i;
    a(k) = 2.0;

    if ( i < n )
      k = k + 1;
      row(k) = i;
      col(k) = i + 1;
      a(k) = -1.0;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [hCBodySurf, hCBodySurfXForm] = ma_initOrbPlot(hFig, orbitDispAxes, bodyInfo)
    %ma_initOrbPlot Summary of this function goes here
    %   Detailed explanation goes here

    hold(orbitDispAxes,'on');
    set(orbitDispAxes,'XTickLabel',[]);
    set(orbitDispAxes,'YTickLabel',[]);
    set(orbitDispAxes,'ZTickLabel',[]);
    grid(orbitDispAxes,'on');

    axis(orbitDispAxes, 'equal');

    axis(orbitDispAxes, 'tight');

    if(~isempty(bodyInfo))
        dRad = bodyInfo.radius;
        [X,Y,Z] = sphere(50);

        if(not(isempty(bodyInfo.surftexturefile)))
            try
                I = bodyInfo.getSurfaceTexture();
                if(not(isempty(I)) && not(any(any(any(isnan(I))))))
                    hCBodySurf = surf(orbitDispAxes, dRad*X,dRad*Y,dRad*Z, 'CData',I, 'FaceColor','texturemap', 'BackFaceLighting','lit', 'FaceLighting','gouraud', 'EdgeLighting','gouraud', 'LineStyle','none');
                else
                    hCBodySurf = createUntexturedSphere(bodyInfo, orbitDispAxes, dRad, X, Y, Z);
                end

            catch ME
                if(exist('hCBodySurf','var'))
                    delete(hCBodySurf); %#ok<NODEF>
                    warning('Could not create textured celestial body "%s".  Error message: \n\n%s', bodyInfo.name, ME.message);
                end

                hCBodySurf = createUntexturedSphere(bodyInfo, orbitDispAxes, dRad, X, Y, Z);
            end

        else
            hCBodySurf = createUntexturedSphere(bodyInfo, orbitDispAxes, dRad, X, Y, Z);
        end

        hCBodySurfXForm = hgtransform('Parent', orbitDispAxes);
        set(hCBodySurf,'Parent',hCBodySurfXForm);
        material(hCBodySurf,'dull');
    else
        hCBodySurf = [];
    end

    mColor = colorFromColorMap(bodyInfo.bodycolor);
    plot3(orbitDispAxes, 0, 0, 0,'Marker','o','MarkerEdgeColor',mColor,'MarkerFaceColor',mColor,'MarkerSize',3);

    hold(orbitDispAxes, 'off');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = cart2sph2(varargin)
%CART2SPH2 Convert cartesian coordinates to spherical coordinates.
%
%   [THETA PHI RHO] = cart2sph2([X Y Z])
%   [THETA PHI RHO] = cart2sph2(X, Y, Z)
%
%   The following convention is used:
%   THETA is the colatitude, in radians, 0 for north pole, +pi for south
%   pole, pi/2 for points with z=0.
%   PHI is the azimuth, in radians, defined as matlab cart2sph: angle from
%   Ox axis, counted counter-clockwise.
%   RHO is the distance of the point to the origin.
%   Discussion on choice for convention can be found at:
%   http://www.physics.oregonstate.edu/bridge/papers/spherical.pdf
%
%   Example:
%   cart2sph2([1 0 0])  returns [pi/2 0 1];
%   cart2sph2([1 1 0])  returns [pi/2 pi/4 sqrt(2)];
%   cart2sph2([0 0 1])  returns [0 0 1];
%
%   See also
%     angles3d, sph2cart2, cart2sph, cart2sph2d
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-02-18
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

% parse input angles based on input argument number
if length(varargin) == 1
    xyz = varargin{1};
elseif length(varargin) == 3
    xyz = [varargin{1} varargin{2} varargin{3}];
end

% ensure z coordinate is specified
if size(xyz, 2) == 2
    xyz(:,3) = 1;
end

% convert to spherical coordinates
[p, t, r] = cart2sph(xyz(:,1), xyz(:,2), xyz(:,3));

% format output arguments
if nargout == 1 || nargout == 0
    varargout{1} = [pi/2-t p r];
elseif nargout==2
    varargout{1} = pi/2-t;
    varargout{2} = p;
else
    varargout{1} = pi/2-t;
    varargout{2} = p;
    varargout{3} = r;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [stats] = est_checkMVARStability(varargin)
%
% Test the stability of a fitted VAR model. See [1-2] for mathematical
% details on testing VAR stability. A stable VAR process is also a
% stationary VAR process [2].
%
% Inputs:
%
%   EEG:        EEGLAB data structure
%   MODEL:      SIFT MODEL structure
%   typeproc:   reserved for future use. Use 0
%
% Optional:
%
%   <Name,Value> pairs containing model fitting parameters. See
%   est_fitMVAR(). Generally, these should be left unspecified.
%
% Outputs:
%
%   stats
%       .stability:  [numwindows x 1] vector of results of stability tests. 1
%                    indicates stable VAR process for that window, 0 indicates
%                    an unstable VAR process.
%
%       .lambda:     [numwindows x nchs*morder] matrix of eigenvalues of VAR
%                    process. All eigenvalues should be < 1 for stable VAR process
%
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%     Theoretical Handbook and User Manual. Chapters 3,6.
%     Available at: http://www.sccn.ucsd.edu/wiki/Sift
% [2] Lutkepohl, H. (2007) New Introduction to Time Series Analysis.
%     Springer.
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

g = arg_define([0 2],varargin, ...
        arg_norep({'EEG','ALLEEG'},mandatory,[],'EEGLAB dataset'), ...
        arg_norep({'MODEL','Model'},mandatory,[],'MVAR MODEL object'), ...
        arg_nogui({'winStartIdx','WindowStartIndices'},[],[],'Starting indices for windows. This is a vector of sample points (start of windows) at which to estimate windowed VAR model. Default is empty (use all windows)','cat','Options'), ...
        arg({'prctWinToSample','WindowSamplePercent'},100,[1 100],'Percent of windows to sample','cat','Options'), ...
        arg({'verb','VerbosityLevel'},2,{int32(0) int32(1) int32(2)},'Verbosity level. 0 = no output, 1 = text, 2 = graphical') ...
        );

% commit EEG and MODEL variables to workspace
[data g] = hlp_splitstruct(g,{'EEG','MODEL'});
arg_toworkspace(data);
clear data;

morder = MODEL.morder;

% window size in points
% winLenPnts = floor(MODEL.winlen*EEG.srate);

if isempty(g.winStartIdx)
    % starting point of each window (points)
    g.winStartIdx  = round(MODEL.winStartTimes*EEG.srate)+1;
end

if g.prctWinToSample<100
    % randomly select percentage of windows to work with
    randwin = randperm(length(g.winStartIdx));
    randwin = sort(randwin(1:ceil(length(g.winStartIdx)*g.prctWinToSample/100)));
    g.winStartIdx = g.winStartIdx(randwin);
    g.winArrayIndex = randwin;
end

% get the array indices of the windows we are working with
g.winArrayIndex = getindex(MODEL.winStartTimes,(g.winStartIdx-1)/EEG.srate);

% initialize waitbar
if g.verb==2
    waitbarTitle = sprintf('Checking stability %s...', ...
        fastif(isempty(EEG.condition),'',['for ' EEG.condition]));

    multiWaitbar(waitbarTitle,'Reset');
    multiWaitbar(waitbarTitle,'ResetCancel',true);
    multiWaitbar(waitbarTitle, ...
                 'Color', hlp_getNextUniqueColor, ...
                 'CanCancel','on', ...
                 'CancelFcn',@(a,b) disp('[Cancel requested. Please wait...]'));
end

numWins = length(g.winStartIdx);

stats.stability = zeros(1,numWins);
[nchs Mp] = size(MODEL.AR{1});
stats.lambda = zeros(numWins,Mp);
%lambda = [];
I = eye(nchs*morder-nchs,nchs*morder-nchs);
O = zeros(nchs*morder-nchs,nchs);
for t=1:numWins
    % get the array index of the window we are working with
    winArrIdx = g.winArrayIndex(t);

    % rewrite VAR[p] process as VAR[1]
    A = [MODEL.AR{winArrIdx} ; [I O]];
    stats.lambda(t,:) = log(abs(eig(A)));
    stats.stability(t) = all(stats.lambda(t,:)<0);

    if g.verb==2
        % update waitbar
        drawnow;
        cancel = multiWaitbar(waitbarTitle,t/numWins);
        if cancel && hlp_confirmWaitbarCancel(waitbarTitle)
            stats = [];
            return;
        end
    end

end

stats.winStartIdx = g.winStartIdx;
stats.winStartTimes = MODEL.winStartTimes(g.winArrayIndex);
stats.winArrayIndex = g.winArrayIndex;

% clean up
if g.verb==2
    multiWaitbar(waitbarTitle,'Close');
end

repo: REST