text:
function [text,AIC,BIC]=egarch_display(parameters,ll,vcv,data,p,o,q,errorType)
% Display parameters, tstats, pvals, log-likelihood and AIC/BIC
% from estimates of a EGARCH(P,O,Q) produced using egarch
%
% USAGE:
%   [TEXT] = egarch_display(PARAMETERS,LL,VCV,DATA,P,O,Q)
%   [TEXT,AIC,BIC] = egarch_display(PARAMETERS,LL,VCV,DATA,P,O,Q,ERRORTYPE)
%
% INPUTS:
%   PARAMETERS    - A 1+p+o+q column vector of parameters with
%                   [omega alpha(1) ... alpha(p) gamma(1) ... gamma(o) beta(1)
%                       ... beta(q) [nu lambda]]'
%   LL            - The log likelihood at the optimum
%   VCV           - Non-robust standard errors (inverse Hessian)
%   DATA          - A column of mean zero data
%   P             - Positive, scalar integer representing the number of
%                   symmetric innovations
%   O             - Non-negative scalar integer representing the number
%                   of asymmetric innovations (0 for symmetric processes)
%   Q             - Non-negative, scalar integer representing the number
%                   of lags of conditional variance (0 for ARCH)
%   ERRORTYPE     - [OPTIONAL] The error distribution used, valid types are:
%                     'NORMAL'    - Gaussian Innovations [DEFAULT]
%                     'STUDENTST' - T distributed errors
%                     'GED'       - Generalized Error Distribution
%                     'SKEWT'     - Skewed T distribution
%
% OUTPUTS:
%   TEXT          - Character matrix with the formatted parameters of the model
%   AIC           - Aikake Information Criteria computed from the LL
%   BIC           - Schwartz/Bayesian Information Criteria computed from the LL
%
%  See also EGARCH

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 9/1/2005

%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETER CHECKING
%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 5
        o=0;
        q=0;
        errorType='NORMAL';
    case 6
        q=0;
        errorType='NORMAL';
    case 7
        errorType='NORMAL';
    case 8
    otherwise
        error('5 to 8 inputs required')
end
if nargin==7
elseif nargin==8
    errorType='NORMAL';
end
if isempty(errorType)
    errorType='NORMAL';
end
% parameters, N by 1, real
if any(~isreal(parameters)) || size(parameters,2)~=1
    error('PARAMETERS must be a column vector.')
end
% LL
if ~isscalar(ll) || ~isreal(ll)
    error('LL must be a scalar.')
end
% VCV
if size(vcv,2)~=size(vcv,1) || any(min(eig(vcv))<=0) || size(vcv,1)~=length(parameters)
    error('VCV must be a square positive definite matrix compatible with PARAMETERS.')
end
% data
if any(~isreal(data)) || size(data,2)~=1
    error('DATA must be a T by 1 column vector.')
end
% p
if ~isscalar(p) || p<1 || floor(p)~=p
    error('P must be a postitive scalar integer')
end
% o
if ~isscalar(o) || o<0 || floor(o)~=o
    error('O must be a non-negative scalar integer')
end
% q
if ~isscalar(q) || q<0 || floor(q)~=q
    error('Q must be a non-negative scalar integer')
end
% errorType
if ~ischar(errorType)
    errorType=[];
end
switch errorType
    case 'NORMAL'
        errorType=1;
        extraP=0;
    case 'STUDENTST'
        errorType=2;
        extraP=1;
    case 'GED'
        errorType=3;
        extraP=1;
    case 'SKEWT'
        errorType=4;
        extraP=2;
    otherwise
        error('ERRORTYPE is not one of the supported distributions')
end
if length(parameters)~=(1+p+o+q+extraP)
    error('Size of PARAMETERS is not compatible with input P, O, Q')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETER CHECKING
%%%%%%%%%%%%%%%%%%%%%%%%%%

model_name = ['EGARCH(' num2str(p) ',' num2str(o) ',' num2str(q) ')'];

parameters_text=num2str(parameters,'%4.4f');
stderr=sqrt(diag(vcv));
stderr_text=num2str(stderr,'%4.4f');
tstats=parameters./stderr;
tstats_text=num2str(tstats,'%4.4f');
pvals=2-2*normcdf(abs(tstats));
pvals_text=num2str(pvals,'%4.4f');

T=length(data);
AIC = -ll/T+2*length(parameters)/T;
BIC = -ll/T+log(T)*length(parameters)/T;

% Format the output
text=[];
text{1,1}= ' ';
text{2,1}= ' ';
text{3,1}=repmat('-',1,50);
text{4,1}=model_name;
text{5,1}=repmat('-',1,50);
text{6,1}=' ';
text{7,1}=['Loglikelihood: ' sprintf('%1.2f',ll)];
text{8,1}=['AIC: ' sprintf('%1.4f',AIC)];
text{9,1}=['BIC: ' sprintf('%1.4f',BIC)];
text{10,1}=' ';

for i=1:size(text,1);
    disp(text{i,1})
end

% Format the parameter table, need to right align everything
K=size(parameters_text,1);
for i=1:K
    N=size(parameters_text,2);
    if any(parameters_text(i,:)==' ')
        numSpace=sum(parameters_text(i,:)==' ');
        parameters_text(i,:)=[repmat(' ',1,numSpace) parameters_text(i,1:N-numSpace)];
    end
    N=size(stderr_text,2);
    if any(stderr_text(i,:)==' ')
        numSpace=sum(stderr_text(i,:)==' ');
        stderr_text(i,:)=[repmat(' ',1,numSpace) stderr_text(i,1:N-numSpace)];
    end
    N=size(tstats_text,2);
    if any(tstats_text(i,:)==' ')
        numSpace=sum(tstats_text(i,:)==' ');
        tstats_text(i,:)=[repmat(' ',1,numSpace) tstats_text(i,1:N-numSpace)];
    end
    N=size(pvals_text,2);
    if any(pvals_text(i,:)==' ')
        numSpace=sum(pvals_text(i,:)==' ');
        pvals_text(i,:)=[repmat(' ',1,numSpace) pvals_text(i,1:N-numSpace)];
    end
end
% Append column labels
labels={' Parameters','   Std. Err.','     T-stat','      P-val'};
cols = {parameters_text,stderr_text,tstats_text,pvals_text};
for i=1:length(labels);
    text1=labels{i};
    text2=cols{i};
    maxcols=max(size(text1,2),size(text2,2));
    if size(text1,2)<maxcols
        labels{i} = [repmat(' ',1,maxcols-size(text1,2)) text1 ];
    end
    if size(text2,2)<maxcols
        cols{i} = [repmat(' ',size(text2,1),maxcols-size(text2,2)) text2 ];
    end
end

% Finally variable names
variable_names=cell(length(parameters),1);
variable_names{1}='omega';
index=2;
for i=1:p
    variable_names{index}=['alpha(' num2str(i) ')'];
    index=index+1;
end
for i=1:o
    variable_names{index}=['gamma(' num2str(i) ')'];
    index=index+1;
end
for i=1:q
    variable_names{index}=['beta(' num2str(i) ')'];
    index=index+1;
end

maxVarNameLength=max(cellfun('length',variable_names));
for i=1:length(variable_names)
    variable_names{i} = [repmat(' ',1,maxVarNameLength-length(variable_names{i})) variable_names{i}];
end
variable_names=cell2mat(variable_names(:));

% Output the parameters
outmat=strvcat(' ',variable_names);
for i=1:4
    outmat=[outmat repmat('  ',size(outmat,1),1) strvcat(labels{i},cols{i})];
end

text2=[];
for i=1:length(text)
    text2=strvcat(text2,text{i});
end
text=strvcat(text2,outmat);

disp(outmat)

author: bashtage

repo: mfe-toolbox


text:
% Find cycles of length 4 in a graph;
% Note: Valid for an undirected graph only
%
% INPUTS: adjacency matrix
% OUTPUTS: number of 4-cycles, for which no edges repeat in the cycle

% Other routines used: numEdges.m, numConnTriples.m, cycles3.m
% GB: last updated January 28, 2016

function l4 = cycles4(adj)

	l4 = trace(adj^4) - 2*numEdges(adj) - 4*numConnTriples(adj) - 8*cycles3(adj);

	l4 = l4/8;

author: aeolianine

repo: octave-networks-toolbox


text:
function [n,J] = dyadlength(x)
% dyadlength -- Find length and dyadic length of array
%  Usage
%    [n,J] = dyadlength(x)
%  Inputs
%    x    array of length n = 2^J (hopefully)
%  Outputs
%    n    length(x)
%    J    least power of two greater than n
%
%  Side Effects
%    A warning is issued if n is not a power of 2.
%
%  See Also
%    quadlength, dyad, dyad2ix
%
  n = length(x) ;
  J = ceil(log(n)/log(2));
  if 2^J ~= n ,
      disp('Warning in dyadlength: n != 2^J')
  end

%
%  Part of Wavelab Version 850
%  Built Tue Jan  3 13:20:40 EST 2006
%  This is Copyrighted Material
%  For Copying permissions see COPYING.m
%  Comments? e-mail wavelab@stat.stanford.edu

author: BUNPC


text:
function sol=polynsolve(polyset,ord,varnames,tol)
% POLYNSOLVE - attempt exact solution of multivariate polynomial system
%
% usage: sol=polynsolve(polyset)
%
% INPUTS: polyset, ord, varnames, tol
%  polyset is a cell array of polynomials in string or coefficient form
%    that is acceptable input for petschel.groebner.m
%  ord (optional) is the preferred ordering
%  varnames (optional) is the list of variable names if not {'x1','x2',...}
%  tol (optional) is the default zero tolerance
%    (see petschel.groebner.m for further details on ord, varnames, tol)
%
% OUTPUTS: sol
%  sol is an array containing the solutions to {polyset{:}=0}:
%    sol(i,j) is the value of xj in the i'th solution
%    If there are infinitely many solutions, sol(i,j)=NaN
%
% ALGORITHM:
%  Uses petschel.groebner bases (lex order).  If any one-variable polynomials
%  result, solve them and substitute back into the equations.
%  If "1" results, the system is not solvable.  If any multivariate polys
%  remain, those variables have an infinite number of solutions.
%
% KNOWN BUGS
%  See petschel.groebner.m for details
%
% SEE ALSO:
%  petschel.groebner, petschel.poly2str, petschel.str2poly

% Author: Ben Petschel 23/6/2009
%
% Change history:
%  23/6/2009 - first release
%  30/3/2010 - change poly representation from n-dim to rectangular array

if nargin<4

    tol = 0;

    if nargin<3

        varnames = {};

    end

end

if (nargin<2) || isempty(ord)

    ord = 'lex';

end

if (numel(polyset)>0) && ischar(polyset{1})

    polyset = petschel.str2poly(polyset,varnames);

end

gbasis=petschel.groebner(polyset,ord,varnames,tol);

% if petschel.groebner is all linear terms, it is a solution, otherwise find
% 1-variable polynomials and solve them, substituting the solutions into
% the other equations
sol = [];

i=1;

keepgoing = true;

while keepgoing && (i<=numel(gbasis))

    % search for 1-variable polynomials
    d = size(gbasis{i},2)-1;

    if numel(sol)<d
        % sol will always be a row vector, until possibly the final step
        % make sure sol has enough possible variables
        sol = [sol,nan(1,d-numel(sol))];

    end

    [tf,n,Q]=ispoly1(gbasis{i});

    if tf
        % attempt solution
        if length(Q)==1
            % have equation 1==0, so no solution
            sol = [];

            keepgoing = false;

        elseif length(Q)==2
            % Q=[a;b] so have equation a*xn+b=0 (a~=0) and solution is x=-b/a

            sol(n) = -Q(2)/Q(1);

        else
            % have polynomial in xn, so solve and substitute solution
            r = roots(Q);

            gbasisrecur = gbasis;

            sol = [];

            for j=1:length(r)
                % replace equation i with linear term (xn-r(j)==0)

                gbasisrecur{i} = linearterm(r(j),n);

                sol = [sol; petschel.polynsolve(gbasisrecur,ord,varnames,tol)];

                keepgoing = false;

            end

        end

    end

    i = i+1;

end

end % main function polynsolve(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tf,n,Q]=ispoly1(P)
% returns tf=true if P is a 1-variable polynomial
% n is the number of the variable the polynomial is in
% Q is a column vector of the polynomial coefficients

if isempty(P) || all(P(:)==0)

    tf=true;

    n=1;

    Q=0;

elseif size(P,2)==1

    tf=true;

    n=1;

    Q=P;

else

    ind = find(any(P(:,2:end)>0,1));

    if numel(ind)==1
        % is a polynomial in 1 variable if exponents of only 1 variable are >0
        tf=true;

        n=ind;

        Q=[]; % collect coefficients of the polynomial

        Q(P(:,ind+1)+1)=P(:,1);

        Q=Q(end:-1:1); % reverse order to be consistent with ROOTS

    else

        tf=false;

        n=0;

        Q=[];

    end

end

end % helper function ispoly1(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=linearterm(a,n)
% returns coefficient array of a linear term (xn-a)

Q=zeros(2,n+1);

Q(1,1)=-a;

Q(2,1)=1;

Q(2,end)=1;

end % helper function linearterm(...)

author: jmaih


text:
function fx = p18_fun ( n, x )

%*****************************************************************************80
%
%% P18_FUN evaluates the integrand for problem 18.
%
%  Integral:
%
%    Integral ( 0 <= x < +oo ) x^2 * exp ( - x / 2^beta ) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Piessens, Elise de Doncker-Kapenga,
%    Christian Ueberhuber, David Kahaner,
%    QUADPACK: A Subroutine Package for Automatic Integration,
%    Springer, 1983, page 84.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the point at which the integrand
%    is to be evaluated.
%
%    Output, real FX(N), the value of the integrand at X.
%
  beta = 1.0;

  fx(1:n) = x(1:n).^2 .* exp ( - x(1:n) / 2^beta );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% =======================================================================
%     Improved adaptive complex diffusion despeckling filter (NCDF)
% =======================================================================
% DESCRIPTION: Demo for despeckling filter
%
% _________________________________________________________________________
% REFERENCES:
% [1] Rui Bernardes, Cristina Maduro, Pedro Serranho, Adérito Araújo,
%     Sílvia Barbeiro, and José Cunha-Vaz,
%     "Improved adaptive complex diffusion despeckling filter,"
%     Opt. Express 18, 24048-24059 (2010).
%
% Work developed under the research project supported by Fundação para a
% Ciência e a Tecnologia (FCT): PTDC/SAU-BEB/103151/2008 and program
% COMPETE (FCOMP-01-0124-FEDER-010930)
%
% http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-23-24048
% http://www.aibili.pt
%
%
% DEPENDENCIES:
% needs the matlab image processing toolbox
% _________________________________________________________________________
%
%
%
% AIBILI - Association for Innovation and Biomedical Research on Light and
% Image

% Demo
clear all;
close all;

% diffusion time in seconds
TMAX         = .80;

% display parameters
minIntensity = 0;
maxIntensity = 220;
iRoi         = [250  750     250    500];

% read an image
Img_noisy    = imread('img.jpg');

% Apply filter to reduce the speckle noise
[Img_filtered, nIter, dTT] = twodncdf(Img_noisy, TMAX);

% DISPLAY results
figure(1),
imagesc(Img_noisy),    title('original image');
caxis([minIntensity maxIntensity]), axis off, colormap(gray)
rect=[iRoi(3) iRoi(1) iRoi(4)-iRoi(3) iRoi(2)-iRoi(1)];
rectangle('Position',rect,'EdgeColor','y','LineStyle','--','Linewidth',2)

figure(2),
imagesc(Img_noisy),    title('original image');
caxis([minIntensity maxIntensity]), axis off, colormap(hsv)
axis(iRoi([3 4 1 2]))

figure(3),
imagesc(Img_filtered), title('filtered image');
caxis([minIntensity maxIntensity]), axis off, colormap(hsv)
axis(iRoi([3 4 1 2]))

author: Sable

repo: mcbench-benchmarks


text:
function P=ar1spectrum(ar1,period)
%%  AR1 power spectrum.
%
% power=ar1spectrum(ar1,period)
%
%
% (c) Aslak Grinsted 2002-2014
%

% -------------------------------------------------------------------------
%The MIT License (MIT)
%
%Copyright (c) 2014 Aslak Grinsted
%
%Permission is hereby granted, free of charge, to any person obtaining a copy
%of this software and associated documentation files (the "Software"), to deal
%in the Software without restriction, including without limitation the rights
%to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
%copies of the Software, and to permit persons to whom the Software is
%furnished to do so, subject to the following conditions:
%
%The above copyright notice and this permission notice shall be included in
%all copies or substantial portions of the Software.
%
%THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
%IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
%AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
%LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
%OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
%THE SOFTWARE.
%---------------------------------------------------------------------------

freq=1./period;
P=(1-ar1.^2)./(abs(1-ar1.*exp(-2*pi*i*freq))).^2; %http://www.madsci.org/posts/archives/may97/864012045.Eg.r.html
%fixed typo in numerical recipes

author: grinsted

repo: wavelet-coherence


text:
classdef CoreStack < handle
    % CoreStack. A class for storing values on a stack
    %
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        Stack
    end

    methods
        function obj = CoreStack(stack_items)
            if nargin > 0
                obj.Stack = CoreContainerUtilities.ConvertToSet(stack_items);
            else
                obj.Stack = cell.empty;
            end
        end

        function Push(obj, stack_items)
            % Places one or more items on the stack. stack_items can be a single
            % object, a set or an array of objects
            obj.Stack = [obj.Stack CoreContainerUtilities.ConvertToSet(stack_items)];
        end

        function item = Pop(obj)
            % Retrieves the item from the top of the stack
            item = obj.Stack{end};
            obj.Stack(end) = [];
        end

        function is_empty = IsEmpty(obj)
            % Returns true if there are no items in the stack
            is_empty = isempty(obj.Stack);
        end

        function all_values = GetAndClear(obj)
            % Clears the stack and returns any remaining values in a set
            all_values = obj.Stack;
            obj.Stack = [];
        end

        function field_values = GetField(obj, field_name)
            % Returns the values of the field field_name from every item in the
            % stack
            field_values = CoreContainerUtilities.GetFieldValuesFromSet(obj.Stack, field_name);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function M = obliquecomplexfactory(n, m, transposed)
% Returns a manifold struct defining complex matrices w/ unit-norm columns.
%
% function M = obliquecomplexfactory(n, m)
% function M = obliquecomplexfactory(n, m, transposed)
%
% Oblique manifold: deals with complex matrices of size n x m such that
% each column has unit 2-norm, i.e., is a point on the unit sphere in C^n.
% The geometry is a product geometry of m unit spheres in C^n. For the
% metric, C^n is treated as R^(2n), so that the real part and imaginary
% parts are treated separately as 2n real coordinates. As such, the complex
% oblique manifold is a Riemannian submanifold of (R^2)^(n x m), with the
% usual metric <u, v> = real(u'*v).
%
% If transposed is set to true (it is false by default), then the matrices
% are transposed: a point Y on the manifold is a matrix of size m x n and
% each row has unit 2-norm. It is the same geometry, just a different
% representation.
%
% In transposed form, a point Y is such that Y*Y' is a Hermitian, positive
% semidefinite matrix of size m and of rank at most n, such that all the
% diagonal entries are equal to 1.
%
% Note: obliquecomplexfactory(1, n, true) is equivalent to (but potentially
% slower than) complexcirclefactory(n).
%
% See also: spherecomplexfactory complexcirclefactory obliquefactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Sep. 3, 2014.
% Contributors:
% Change log:
%
%   Oct. 21, 2016 (NB)
%       Formatted for inclusion in Manopt release.
%
%   July 20, 2017 (NB)
%       Distance function is now accurate for close-by points. See notes
%       inside the spherefactory file for details. Also improves distances
%       computation as part of the log function.

    if ~exist('transposed', 'var') || isempty(transposed)
        transposed = false;
    end

    if transposed
        trnsp = @(X) X.';
    else
        trnsp = @(X) X;
    end

    M.name = @() sprintf('Complex oblique manifold COB(%d, %d)', n, m);

    M.dim = @() (2*n-1)*m;

    M.inner = @(x, d1, d2) real(d1(:)'*d2(:));

    M.norm = @(x, d) norm(d(:));

    M.dist = @(x, y) norm(real(2*asin(.5*sqrt(sum(trnsp(abs(x - y).^2), 1)))));

    M.typicaldist = @() pi*sqrt(m);

    M.proj = @(X, U) trnsp(projection(trnsp(X), trnsp(U)));

    M.tangent = M.proj;

    % For Riemannian submanifolds, converting a Euclidean gradient into a
    % Riemannian gradient amounts to an orthogonal projection.
    M.egrad2rgrad = M.proj;

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, U)
        X = trnsp(X);
        egrad = trnsp(egrad);
        ehess = trnsp(ehess);
        U = trnsp(U);

        PXehess = projection(X, ehess);
        inners = sum(real(conj(X).*egrad), 1);
        rhess = PXehess - bsxfun(@times, U, inners);

        rhess = trnsp(rhess);
    end

    M.exp = @exponential;
    % Exponential on the complex oblique manifold
    function y = exponential(x, d, t)
        x = trnsp(x);
        d = trnsp(d);

        if nargin == 2
            % t = 1;
            td = d;
        else
            td = t*d;
        end

        nrm_td = sqrt(sum(real(td).^2 + imag(td).^2, 1));

        y = bsxfun(@times, x, cos(nrm_td)) + ...
            bsxfun(@times, td, sinxoverx(nrm_td));

        y = trnsp(y);
    end

    M.log = @logarithm;
    function v = logarithm(x1, x2)
        x1 = trnsp(x1);
        x2 = trnsp(x2);

        v = projection(x1, x2 - x1);
        dists = real(2*asin(.5*sqrt(sum(trnsp(abs(x2 - x1).^2), 1))));
        norms = sqrt(sum(real(v).^2 + imag(v).^2, 1));
        factors = dists./norms;
        % For very close points, dists is almost equal to norms, but
        % because they are both almost zero, the division above can return
        % NaN's. To avoid that, we force those ratios to 1.
        factors(dists <= 1e-10) = 1;
        v = bsxfun(@times, v, factors);

        v = trnsp(v);
    end

    M.retr = @retraction;
    % Retraction on the oblique manifold
    function y = retraction(x, d, t)
        x = trnsp(x);
        d = trnsp(d);

        if nargin < 3
            td = d;
        else
            td = t*d;
        end

        y = normalize_columns(x + td);

        y = trnsp(y);
    end

    M.hash = @(x) ['z' hashmd5([real(x(:)) ; imag(x(:))])];

    M.rand = @() trnsp(random(n, m));

    M.randvec = @(x) trnsp(randomvec(n, m, trnsp(x)));

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) trnsp(zeros(n, m));

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

    M.pairmean = @pairmean;
    function y = pairmean(x1, x2)
        y = trnsp(x1+x2);
        y = normalize_columns(y);
        y = trnsp(y);
    end

    % vec returns a vector representation of an input tangent vector which
    % is represented as a matrix. mat returns the original matrix
    % representation of the input vector representation of a tangent
    % vector. vec and mat are thus inverse of each other. They are
    % furthermore isometries between a subspace of R^2nm and the tangent
    % space at x.
    vect = @(X) X(:);
    M.vec = @(x, u_mat) [vect(real(trnsp(u_mat))) ; ...
                         vect(imag(trnsp(u_mat)))];
    M.mat = @(x, u_vec)    trnsp(reshape(u_vec(1:(n*m)),     [n, m])) + ...
                        1i*trnsp(reshape(u_vec((n*m+1):end), [n, m]));
    M.vecmatareisometries = @() true;

end

% Given a matrix X, returns the same matrix but with each column scaled so
% that they have unit 2-norm.
function X = normalize_columns(X)
    norms = sqrt(sum(real(X).^2 + imag(X).^2, 1));
    X = bsxfun(@times, X, 1./norms);
end

% Orthogonal projection of the ambient vector H onto the tangent space at X
function PXH = projection(X, H)

    % Compute the inner product between each vector H(:, i) with its root
    % point X(:, i), that is, real(X(:, i)' * H(:, i)).
    % Returns a row vector.
    inners = real(sum(conj(X).*H, 1));

    % Subtract from H the components of the H(:, i)'s that are parallel to
    % the root points X(:, i).
    PXH = H - bsxfun(@times, X, inners);

end

% Uniform random sampling on the sphere.
function x = random(n, m)

    x = normalize_columns(randn(n, m) + 1i*randn(n, m));

end

% Random normalized tangent vector at x.
function d = randomvec(n, m, x)

    d = randn(n, m) + 1i*randn(n, m);
    d = projection(x, d);
    d = d / norm(d(:));

end

author: NicolasBoumal

repo: manopt


text:
function [val, subCoords] = rmCoordsGet(viewType, model, param, coords)
% rmCoordsGet - wrapper for rmGet to limit values to certain coordinates
%
% val = rmCoordsGet(viewType, model, param, coords)
%

% 2008/02 SOD: split of from rmPredictedTSeries.
% 2008/06 RAS: changed view argument to viewType: this will prevent the
% rmPlotGUI code from having to carry around a large view structure,
% substantially reducing its memory footprint.
% 2009/01 RAS: allows the coords variable for gray views to be a 3xN list
% of coordinates, in addition to a vector of gray node indices.
% 2018/02 RKL: returns the coords as the 2nd output. because they will be a
% subset of the coords that are inputed
if ~exist('viewType','var') || isempty(viewType),   error('Need view type'); end;
if ~exist('model','var')    || isempty(model),      error('Need model');     end;
if ~exist('param','var')    || isempty(param),      error('Need param');     end;
if ~exist('coords','var'),                          error('Need coords');    end;
if isempty(coords), val = []; return; end
% allow the model to be a cell array -- just take the first entry
if iscell(model), model = model{1};  end

% if we input a view instead of a viewType, sort this out
if isstruct(viewType), vw = viewType; viewType = vw.viewType; end

tmp = rmGet(model, param);

switch lower(viewType),
    case 'inplane'

        % normally inplane data will be 3D, but if we have only one slice
        % it will be 2d. we need to check.
        if ~exist('vw', 'var') || isempty(vw), vw = getSelectedInplane; end
        dims = length(viewGet(vw, 'anat size'));

        % If there is one value per voxel, then we expect a matrix of
        % parameter values with the same dimensionality as our inplane
        % anatomy (usually 3d, unless we have only 1 slice in which case it
        % is 2d). So we get the values of this matrix indexed by the 3D
        % coords. So for coord (x, y, z), we would like
        %   val = tmp(x, y, z);
        % In the case of a 2D anat, we would like
        %   val = tmp(x,y);
        if length(size(tmp)) <= dims
            val = zeros(size(coords, 2), 1);
            for n = 1:length(val),
                if dims == 3,
                    val(n) = tmp(coords(1,n), coords(2,n), coords(3,n));
                elseif dims == 2
                    val(n) = tmp(coords(1,n), coords(2,n));
                end
            end;
        end

        % If there is more than one value per voxel, then we expect a
        % matrix of parameter values that has one more dimension than our inplane anatomy.
        % In this case we get all the values of this matrix whose first
        % 3 dimensionpars are indexed by the 3D coords. So for coord (x, y, z)
        % we would like
        %   val = tmp(x, y, z, :);
        % If our inplane happens to be 2d (because we have only slice),
        % then we would like
        %   val = tmp(x, y, :);

        if length(size(tmp)) > dims
            val = zeros(size(coords, 2), size(tmp, dims));
            for n = 1:size(val,1),
                 if dims == 3,
                     val(n, :) = tmp(coords(1,n), coords(2,n), coords(3,n), :);
                 elseif dims == 2
                     val(n, :) = tmp(coords(1,n), coords(2,n), :);
                 end
            end;
        end

    case 'gray'
        % allow 3xN gray coords to be specified, as well as gray node
        % indices:
        if size(coords, 1)==3
            % 3xN coords specification: remap into indices
            vw = getSelectedGray;
            allCoords = viewGet(vw, 'coords');
            [subCoords coords] = intersectCols(allCoords, coords); %#ok<ASGLU>
        end

        if numel(size(tmp)) == 2,
            val = tmp(coords);
        else
            if length(coords)==1
                % the squeeze command will mis-orient fields like beta, by
                % permuting across multiple dimensions. We want it to
                % voxels x predictors.
                val = permute(tmp(1,coords,:), [2 3 1]);
            else
                val = squeeze(tmp(1,coords,:));
            end
        end

    otherwise, error('Invalid view type.')
end;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function grad = B_splice_single_sentence2(future_grad, context, mask)
[dim,nFr,nSeg] = size(future_grad);
dim = dim/context;
half_ctx = (context-1)/2;

future_grad2 = reshape(future_grad, dim, context, nFr, nSeg);
future_grad2(:,:,end+1,:) = 0;

for i=-half_ctx:half_ctx
    all_idx{i+half_ctx+1} = min(nFr,max(1, (1:nFr) + i));
end
all_idx2 = cell2mat(all_idx');

future_grad2 = all_idx2;
future_grad2(:,end+1) = 0;

grad_idx = zeros(nFr, context);
grad_idx1 = [];
grad_idxend = [];
for i=-half_ctx:half_ctx
    offset = (nFr+1)*(i+half_ctx);
    % curr_future_grad = future_grad( (i+half_ctx)*dim+1 : (i+half_ctx+1)*dim, :,:);
    if i<0
        grad_idx(1:nFr+i, i+half_ctx+1) = (-i+1:nFr) + offset;
        %grad(:,1:end+i,:) = grad(:,1:end+i,:) + curr_future_grad(:, -i+1:end,:);
        grad_idx1 = [grad_idx1 (1:-i)+offset];
    elseif i>0
        grad_idx(i+1:nFr, i+half_ctx+1) = (1:nFr-i) + offset;

        %grad(:,i+1:end,:) = grad(:,i+1:end,:) + curr_future_grad(:, 1:end-i,:);
        %tmp2 = [tmp2 curr_future_grad(:, end-i+1:end,:)];
        grad_idxend = [grad_idxend (nFr-i+1:nFr)+offset];
    else
        grad_idx(1:nFr, i+half_ctx+1) = (1:nFr) + offset;
        %grad = grad + curr_future_grad;
    end
end
grad_idx(grad_idx==0) = nFr*context+1;

% future_grad2 = permute(future_grad2, [1 4 2 3]);
future_grad2 = reshape(future_grad2, dim*nSeg, context*(nFr+1));
grad = future_grad2(:,grad_idx);
grad = reshape(grad, dim, context, nFr);
grad = squeeze(sum(grad,2));
grad = permute(grad, [1 3 2]);

grad = future_grad2(grad_idx');

% grad(:,1,:) = grad(:,1,:) + sum(tmp1,2);
% grad(:,end,:) = grad(:,end,:) + sum(tmp2,2);

end

author: singaxiong

repo: SignalGraph


text:
classdef NicoletFile < handle
  % NICOLETFILE  Reading Nicolet .e files.
  %
  %   GETTING DATA
  %   You can load data from the .e file using the GETDATA method. The
  %   inputs to the method are the object, the segment of the data file
  %   that you want to load data from, the min, and max index you want to
  %   retrieve, and a vector of channels that you want to retrieve.
  %
  %   Example:
  %     OUT = GETDATA(OBJ, 1, [1 1000], 1:10) will return the first 1000
  %     values on the first 10 channels of the first segment of the file.
  %
  %   GETTING Nr OF SAMPLES
  %     Use the GETNRSAMPLES method to find the number of samples per
  %     channel in each data segment.
  %
  %   WARNING!
  %   The .e format allows for changes in the TimeSeries map during the
  %   recording. This results in multiple TSINFO structures. Depending on
  %   where these structures are located in the .e file, the appropriate
  %   TSINFO structure should be used. However, there seems to be a bug in
  %   the Nicolet .e file writer which sometimes renders the TSINFO structures
  %   unreadable on disk (verify with hex-edit).

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Copyright 2013 Trustees of the University of Pennsylvania
  %
  % Licensed under the Apache License, Version 2.0 (the "License");
  % you may not use this file except in compliance with the License.
  % You may obtain a copy of the License at
  %
  % http://www.apache.org/licenses/LICENSE-2.0
  %
  % Unless required by applicable law or agreed to in writing, software
  % distributed under the License is distributed on an "AS IS" BASIS,
  % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  % See the License for the specific language governing permissions and
  % limitations under the License.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % Joost Wagenaar, Jan 2015
  % Cristian Donos, Dec 2015
  % Jan Brogger, Jun 2016
  % Callum Stewart, Apr 2018

  properties
    fileName
    patientInfo
    segments
    eventMarkers
  end

  properties (Hidden)
    sections
    index
    sigInfo
    tsInfos
    chInfo
    notchFreq
    montage
    Qi
    Qii
    allIndexIDs
    useTSinfoIdx = 1
  end

  properties (Constant, Hidden)
   LABELSIZE = 32
   TSLABELSIZE = 64
   UNITSIZE = 16
   ITEMNAMESIZE  = 64
  end

  methods
    function obj = NicoletFile(filename)
      h = fopen(filename,'r','ieee-le');

      [folder, name, ext] = fileparts(filename);
      assert(strcmp(ext,'.e'), 'File extention must be .e');
      if isempty(folder)
        filename = fullfile(pwd,filename);
      end

      obj.fileName = filename;
      % Get init
      misc1 = fread(h,5, 'uint32'); %#ok<NASGU>
      unknown = fread(h,1,'uint32'); %#ok<NASGU>
      indexIdx = fread(h,1,'uint32');

      % Get TAGS structure and Channel IDS
      fseek(h, 172,'bof');
      nrTags = fread(h,1, 'uint32');
      Tags = struct();
      for i = 1:nrTags
        Tags(i).tag = deblank(cast(fread(h, 40, 'uint16'),'char')');
        Tags(i).index = fread(h,1,'uint32');
        switch Tags(i).tag
          case 'ExtraDataTags'
            Tags(i).IDStr = 'ExtraDataTags';
          case 'SegmentStream'
            Tags(i).IDStr = 'SegmentStream';
          case 'DataStream'
            Tags(i).IDStr = 'DataStream';
          case 'InfoChangeStream'
            Tags(i).IDStr = 'InfoChangeStream';
          case 'InfoGuids'
            Tags(i).IDStr = 'InfoGuids';
          case '{A271CCCB-515D-4590-B6A1-DC170C8D6EE2}'
            Tags(i).IDStr = 'TSGUID';
          case '{8A19AA48-BEA0-40D5-B89F-667FC578D635}'
            Tags(i).IDStr = 'DERIVATIONGUID';
          case '{F824D60C-995E-4D94-9578-893C755ECB99}'
            Tags(i).IDStr = 'FILTERGUID';
          case '{02950361-35BB-4A22-9F0B-C78AAA5DB094}'
            Tags(i).IDStr = 'DISPLAYGUID';
          case '{8E94EF21-70F5-11D3-8F72-00105A9AFD56}'
            Tags(i).IDStr = 'FILEINFOGUID';
          case '{E4138BC0-7733-11D3-8685-0050044DAAB1}'
            Tags(i).IDStr = 'SRINFOGUID';
          case '{C728E565-E5A0-4419-93D2-F6CFC69F3B8F}'
            Tags(i).IDStr = 'EVENTTYPEINFOGUID';
          case '{D01B34A0-9DBD-11D3-93D3-00500400C148}'
            Tags(i).IDStr = 'AUDIOINFOGUID';
          case '{BF7C95EF-6C3B-4E70-9E11-779BFFF58EA7}'
            Tags(i).IDStr = 'CHANNELGUID';
          case '{2DEB82A1-D15F-4770-A4A4-CF03815F52DE}'
            Tags(i).IDStr = 'INPUTGUID';
          case '{5B036022-2EDC-465F-86EC-C0A4AB1A7A91}'
            Tags(i).IDStr = 'INPUTSETTINGSGUID';
          case '{99A636F2-51F7-4B9D-9569-C7D45058431A}'
            Tags(i).IDStr = 'PHOTICGUID';
          case '{55C5E044-5541-4594-9E35-5B3004EF7647}'
            Tags(i).IDStr = 'ERRORGUID';
          case '{223A3CA0-B5AC-43FB-B0A8-74CF8752BDBE}'
            Tags(i).IDStr = 'VIDEOGUID';
          case '{0623B545-38BE-4939-B9D0-55F5E241278D}'
            Tags(i).IDStr = 'DETECTIONPARAMSGUID';
          case '{CE06297D-D9D6-4E4B-8EAC-305EA1243EAB}'
            Tags(i).IDStr = 'PAGEGUID';
          case '{782B34E8-8E51-4BB9-9701-3227BB882A23}'
            Tags(i).IDStr = 'ACCINFOGUID';
          case '{3A6E8546-D144-4B55-A2C7-40DF579ED11E}'
            Tags(i).IDStr = 'RECCTRLGUID';
          case '{D046F2B0-5130-41B1-ABD7-38C12B32FAC3}'
            Tags(i).IDStr = 'GUID TRENDINFOGUID';
          case '{CBEBA8E6-1CDA-4509-B6C2-6AC2EA7DB8F8}'
            Tags(i).IDStr = 'HWINFOGUID';
          case '{E11C4CBA-0753-4655-A1E9-2B2309D1545B}'
            Tags(i).IDStr = 'VIDEOSYNCGUID';
          case '{B9344241-7AC1-42B5-BE9B-B7AFA16CBFA5}'
            Tags(i).IDStr = 'SLEEPSCOREINFOGUID';
          case '{15B41C32-0294-440E-ADFF-DD8B61C8B5AE}'
            Tags(i).IDStr = 'FOURIERSETTINGSGUID';
          case '{024FA81F-6A83-43C8-8C82-241A5501F0A1}'
            Tags(i).IDStr = 'SPECTRUMGUID';
          case '{8032E68A-EA3E-42E8-893E-6E93C59ED515}'
            Tags(i).IDStr = 'SIGNALINFOGUID';
          case '{30950D98-C39C-4352-AF3E-CB17D5B93DED}'
            Tags(i).IDStr = 'SENSORINFOGUID';
          case '{F5D39CD3-A340-4172-A1A3-78B2CDBCCB9F}'
            Tags(i).IDStr = 'DERIVEDSIGNALINFOGUID';
          case '{969FBB89-EE8E-4501-AD40-FB5A448BC4F9}'
            Tags(i).IDStr = 'ARTIFACTINFOGUID';
          case '{02948284-17EC-4538-A7FA-8E18BD65E167}'
            Tags(i).IDStr = 'STUDYINFOGUID';
          case '{D0B3FD0B-49D9-4BF0-8929-296DE5A55910}'
            Tags(i).IDStr = 'PATIENTINFOGUID';
          case '{7842FEF5-A686-459D-8196-769FC0AD99B3}'
            Tags(i).IDStr = 'DOCUMENTINFOGUID';
          case '{BCDAEE87-2496-4DF4-B07C-8B4E31E3C495}'
            Tags(i).IDStr = 'USERSINFOGUID';
          case '{B799F680-72A4-11D3-93D3-00500400C148}'
            Tags(i).IDStr = 'EVENTGUID';
          case '{AF2B3281-7FCE-11D2-B2DE-00104B6FC652}'
            Tags(i).IDStr = 'SHORTSAMPLESGUID';
          case '{89A091B3-972E-4DA2-9266-261B186302A9}'
            Tags(i).IDStr = 'DELAYLINESAMPLESGUID';
          case '{291E2381-B3B4-44D1-BB77-8CF5C24420D7}'
            Tags(i).IDStr = 'GENERALSAMPLESGUID';
          case '{5F11C628-FCCC-4FDD-B429-5EC94CB3AFEB}'
            Tags(i).IDStr = 'FILTERSAMPLESGUID';
          case '{728087F8-73E1-44D1-8882-C770976478A2}'
            Tags(i).IDStr = 'DATEXDATAGUID';
          case '{35F356D9-0F1C-4DFE-8286-D3DB3346FD75}'
            Tags(i).IDStr = 'TESTINFOGUID';

          otherwise
            if isstrprop(Tags(i).tag, 'digit')
              Tags(i).IDStr = num2str(Tags(i).tag);
            else
              Tags(i).IDStr = 'UNKNOWN';
            end
        end

      end

      obj.sections = Tags;

      %% QI index
      fseek(h, 172208,'bof');
      obj.Qi=struct();
      obj.Qi.nrEntries = fread(h,1,'uint32');
      obj.Qi.misc1 = fread(h,1,'uint32');
      obj.Qi.indexIdx = fread(h,1,'uint32');
      obj.Qi.misc3 = fread(h,1,'uint32');
      obj.Qi.LQi = fread(h,1,'uint64')';
      obj.Qi.firstIdx = fread(h,nrTags,'uint64');

      % Don't know what this index is for... Not required to get data and
      % can be huge...

%       fseek(h, 188664,'bof');
%       Qindex  = struct();
%       for i = 1:obj.Qi.LQi
% %         Qindex(i).ftel = ftell(h);
%         Qindex(i).index = fread(h,2,'uint16')';  %4
%         Qindex(i).misc1 = fread(h,1,'uint32');   %8
%         Qindex(i).indexIdx = fread(h,1,'uint32'); %12
%         Qindex(i).misc2 = fread(h,3,'uint32')'; %24
%         Qindex(i).sectionIdx = fread(h,1,'uint32');%28
%         Qindex(i).misc3 = fread(h,1,'uint32'); %32
%         Qindex(i).offset = fread(h,1,'uint64'); % 40
%         Qindex(i).blockL = fread(h,1,'uint32');%44
%         Qindex(i).dataL = fread(h,1,'uint32')';%48
%       end
%       obj.Qi.index = Qindex;

      %% Get Main Index:
      % Index consists of multiple blocks, after each block is the pointer
      % to the next block. Total number of entries is in obj.Qi.nrEntries

      Index = struct();
      curIdx = 0;
      nextIndexPointer = indexIdx;
      fprintf('Parsing index ');
      curIdx2 = 1;
      while curIdx < obj.Qi.nrEntries
        if mod(curIdx2,20)
          fprintf('.');
        else
          fprintf('\n.');
        end

        fseek(h, nextIndexPointer, 'bof');
        nrIdx = fread(h,1, 'uint64');
        Index(curIdx + nrIdx).sectionIdx = 0;   % Preallocate next set of indices
        var = fread(h,3*nrIdx, 'uint64');
        for i = 1: nrIdx
%           Index(curIdx + i).sectionIdx = fread(h,1, 'uint64');
%           Index(curIdx + i).offset = fread(h,1, 'uint64');
%           Index(curIdx + i).blockL = fread(h,1, 'uint32');
%           Index(curIdx + i).sectionL = fread(h,1, 'uint32');
          Index(curIdx + i).sectionIdx = var(3*(i-1)+1);
          Index(curIdx + i).offset = var(3*(i-1)+2);
          Index(curIdx + i).blockL = mod(var(3*(i-1)+3),2^32);
          Index(curIdx + i).sectionL = round(var(3*(i-1)+3)/2^32);
        end
        nextIndexPointer = fread(h,1, 'uint64');
        curIdx = curIdx + i;
        curIdx2=curIdx2+1;

      end
      fprintf('done\n');
      obj.index = Index;
      obj.allIndexIDs = [obj.index.sectionIdx];

    %---READ DYNAMIC PACKETS---%
    dynamicPackets = struct();
    indexIdx = Tags(find(strcmp({Tags.IDStr},'InfoChangeStream'),1)).index;
    offset = Index(indexIdx).offset;
    nrDynamicPackets = Index(indexIdx).sectionL / 48;
    fseek(h, offset, 'bof');

    %Read first only the dynamic packets structure without actual data
    for i = 1: nrDynamicPackets
        dynamicPackets(i).offset = offset+i*48;
        guidmixed = fread(h,16, 'uint8')';
        guidnonmixed = [guidmixed(04), guidmixed(03), guidmixed(02), guidmixed(01), ...
                        guidmixed(06), guidmixed(05), guidmixed(08), guidmixed(07), ...
                        guidmixed(09), guidmixed(10), guidmixed(11), guidmixed(12), ...
                        guidmixed(13), guidmixed(14), guidmixed(15), guidmixed(16)];
        dynamicPackets(i).guid = num2str(guidnonmixed, '%02X');
        dynamicPackets(i).guidAsStr = sprintf('{%02X%02X%02X%02X-%02X%02X-%02X%02X-%02X%02X-%02X%02X%02X%02X%02X%02X}', guidnonmixed);
        dynamicPackets(i).date = datenum(1899,12,31) + fread(h,1,'double');
        dynamicPackets(i).datefrac = fread(h,1,'double');
        dynamicPackets(i).internalOffsetStart = fread(h,1, 'uint64')';
        dynamicPackets(i).packetSize = fread(h,1, 'uint64')';
        dynamicPackets(i).data = zeros(0, 1,'uint8');

        switch dynamicPackets(i).guid
            case 'BF7C95EF6C3B4E709E11779BFFF58EA7'
               dynamicPackets(i).IDStr = 'CHANNELGUID';
            case '8A19AA48BEA040D5B89F667FC578D635'
                dynamicPackets(i).IDStr = 'DERIVATIONGUID';
            case 'F824D60C995E4D949578893C755ECB99'
                dynamicPackets(i).IDStr = 'FILTERGUID';
            case '0295036135BB4A229F0BC78AAA5DB094'
                dynamicPackets(i).IDStr = 'DISPLAYGUID';
            case '782B34E88E514BB997013227BB882A23'
                dynamicPackets(i).IDStr = 'ACCINFOGUID';
            case 'A271CCCB515D4590B6A1DC170C8D6EE2'
                dynamicPackets(i).IDStr = 'TSGUID';
            case 'D01B34A09DBD11D393D300500400C148'
                dynamicPackets(i).IDStr = 'AUDIOINFOGUID';
            otherwise
                dynamicPackets(i).IDStr = 'UNKNOWN';
        end
    end

    %Then read the actual data from the pointers above
    for i = 1: nrDynamicPackets
        %Look up the GUID of this dynamic packet in the Tags
        % to find the section index

        infoIdx = Tags(find(strcmp({Tags.tag},dynamicPackets(i).guidAsStr),1)).index;

        %Matching index segments
        indexInstances = Index([Index.sectionIdx] == infoIdx);

        %Then, treat all these sections as one contiguous memory block
        % and grab this packet across these instances

        internalOffset = 0;
        remainingDataToRead = dynamicPackets(i).packetSize;
        %disp(['Target packet ' dynamicPackets(i).IDStr ' : ' num2str(dynamicPackets(i).internalOffsetStart) ' to ' num2str(dynamicPackets(i).internalOffsetStart+dynamicPackets(i).packetSize) ' target read length ' num2str(remainingDataToRead)]);
        currentTargetStart = dynamicPackets(i).internalOffsetStart;
        for j = 1: size(indexInstances,2)
            currentInstance = indexInstances(j);

            %hitInThisSegment = '';
            if (internalOffset <= currentTargetStart) && (internalOffset+currentInstance.sectionL) >= currentTargetStart

                startAt = currentTargetStart;
                stopAt =  min(startAt+remainingDataToRead, internalOffset+currentInstance.sectionL);
                readLength = stopAt-startAt;

                filePosStart = currentInstance.offset+startAt-internalOffset;
                fseek(h,filePosStart, 'bof');
                dataPart = fread(h,readLength,'uint8=>uint8');
                dynamicPackets(i).data = cat(1, dynamicPackets(i).data, dataPart);

                %hitInThisSegment = ['HIT at  ' num2str(startAt) ' to ' num2str(stopAt)];
                %if (readLength < remainingDataToRead)
                %    hitInThisSegment = [hitInThisSegment ' (partial ' num2str(readLength) ' )'];
                %else
                %    hitInThisSegment = [hitInThisSegment ' (finished - this segment contributed ' num2str(readLength) ' )'];
                %end
                %hitInThisSegment = [hitInThisSegment ' abs file pos ' num2str(filePosStart) ' - ' num2str(filePosStart+readLength)];

                remainingDataToRead = remainingDataToRead-readLength;
                currentTargetStart = currentTargetStart + readLength;

            end
            %disp(['    Index ' num2str(j) ' Offset: ' num2str(internalOffset) ' to ' num2str(internalOffset+currentInstance.sectionL) ' ' num2str(hitInThisSegment)]);

            internalOffset = internalOffset + currentInstance.sectionL;
        end
    end

      %% Get PatientGUID
      info = struct();

      infoProps = { 'patientID', 'firstName','middleName','lastName',...
        'altID','mothersMaidenName','DOB','DOD','street','sexID','phone',...
        'notes','dominance','siteID','suffix','prefix','degree','apartment',...
        'city','state','country','language','height','weight','race','religion',...
        'maritalStatus'};

      infoIdxStruct = Tags(find(strcmp({Tags.IDStr},'PATIENTINFOGUID'),1));
      if ~isempty(infoIdxStruct)
        infoIdx = infoIdxStruct.index;
        indexInstance = Index(find([Index.sectionIdx]==infoIdx,1));
        fseek(h, indexInstance.offset,'bof');
        guid = fread(h, 16, 'uint8'); %#ok<NASGU>
        lSection = fread(h, 1, 'uint64'); %#ok<NASGU>
%        reserved = fread(h, 3, 'uint16'); %#ok<NASGU>
        nrValues = fread(h,1,'uint64');
        nrBstr = fread(h,1,'uint64');

        for i = 1:nrValues
          id = fread(h,1,'uint64');
          switch id
            case {7,8}
              unix_time = (fread(h,1, 'double')*(3600*24)) - 2209161600;% 2208988800; %8
              obj.segments(i).dateStr = datestr(unix_time/86400 + datenum(1970,1,1));
              value = datevec( obj.segments(i).dateStr );
              value = value([3 2 1]);
            case {23,24}
              value = fread(h,1,'double');
              otherwise
              value = 0;
          end
          info.(infoProps{id}) = value;
        end

        strSetup = fread(h,nrBstr*2,'uint64');

        for i=1:2:(nrBstr*2)
          id  = strSetup(i);
          value = deblank(cast(fread(h, strSetup(i+1) + 1, 'uint16'),'char')');
          info.(infoProps{id}) = value;
        end
      end

      obj.patientInfo = info;

      %% Get INFOGUID
      infoIdx = Tags(find(strcmp({Tags.IDStr},'InfoGuids'),1)).index;
      indexInstance = Index(find([Index.sectionIdx]==infoIdx,1));
      fseek(h, indexInstance.offset,'bof');

      % Ignoring, is list of GUIDS in file.

      %% Get SignalInfo (SIGNALINFOGUID): One per file
      SIG_struct = struct();
      sensorIdx = Tags(find(strcmp({Tags.IDStr},'SIGNALINFOGUID'),1)).index;
      indexInstance = Index(find([Index.sectionIdx]==sensorIdx,1));
      fseek(h, indexInstance.offset,'bof');
      SIG_struct.guid = fread(h, 16, 'uint8');
      SIG_struct.name = fread(h, obj.ITEMNAMESIZE, '*char');
      unkown = fread(h, 152, '*char');         %#ok<NASGU>
      fseek(h, 512, 'cof');
      nrIdx = fread(h,1, 'uint16');  %783
      misc1 = fread(h,3, 'uint16'); %#ok<NASGU>
      obj.sigInfo = struct();
      for i = 1: nrIdx
        obj.sigInfo(i).sensorName = deblank(cast(fread(h, obj.LABELSIZE, 'uint16'),'char')');
        obj.sigInfo(i).transducer = deblank(cast(fread(h, obj.UNITSIZE, 'uint16'),'char')');
        obj.sigInfo(i).guid = fread(h, 16, '*uint8');
        obj.sigInfo(i).bBiPolar = logical(fread(h, 1 ,'uint32'));
        obj.sigInfo(i).bAC = logical(fread(h, 1 ,'uint32'));
        obj.sigInfo(i).bHighFilter = logical(fread(h, 1 ,'uint32'));
        obj.sigInfo(i).color =  fread(h, 1 ,'uint32');
        reserved = fread(h, 256, '*char'); %#ok<NASGU>
      end

      %% Get CHANNELINFO (CHANNELGUID)
      CH_struct = struct();
      sensorIdx = Tags(find(strcmp({Tags.IDStr},'CHANNELGUID'),1)).index;
      indexInstance = Index(find([Index.sectionIdx]==sensorIdx,1));
      fseek(h, indexInstance.offset,'bof');
            CH_struct.guid = fread(h, 16, 'uint8');
      CH_struct.name = fread(h, obj.ITEMNAMESIZE, '*char');
      fseek(h, 152, 'cof');
      CH_struct.reserved = fread(h, 16, 'uint8');
      CH_struct.deviceID = fread(h, 16, 'uint8');
      fseek(h, 488, 'cof');

      nrIdx = fread(h,2, 'uint32');  %783
      obj.chInfo = struct();
      for i = 1: nrIdx(2)
        obj.chInfo(i).sensor = deblank(cast(fread(h, obj.LABELSIZE, 'uint16'),'char')');
        obj.chInfo(i).samplingRate = fread(h,1,'double');
        obj.chInfo(i).bOn = logical(fread(h, 1 ,'uint32'));
        obj.chInfo(i).lInputID = fread(h, 1 ,'uint32');
        obj.chInfo(i).lInputSettingID = fread(h,1,'uint32');
        obj.chInfo(i).reserved = fread(h,4,'char');
        fseek(h, 128, 'cof');
      end

      curIdx = 0;
      for i = 1: length(obj.chInfo)
        if obj.chInfo(i).bOn
         obj.chInfo(i).indexID = curIdx;
         curIdx = curIdx+1;
        else
          obj.chInfo(i).indexID = -1;
        end
      end

      %% Get TS info (TSGUID):(One per segment, last used if no new for segment)
      %% To simplify things, we only read the first TSINFO.
	  tsPackets = dynamicPackets(strcmp({dynamicPackets.IDStr},'TSGUID'));

      if isempty(tsPackets)
          warning(['No TSINFO found']);
      else
          obj.tsInfos = {};
          for j = 1:length(tsPackets)

          tsPacket = tsPackets(j);
          elems = typecast(tsPacket.data(753:756),'uint32');
          alloc = typecast(tsPacket.data(757:760),'uint32');

          offset = 761;
          tsInfo = struct();
          for i = 1:elems
              internalOffset = 0;
              tsInfo(i).label = deblank(char(typecast(tsPacket.data(offset:(offset+obj.TSLABELSIZE-1))','uint16')));
              internalOffset = internalOffset + obj.TSLABELSIZE*2;
              tsInfo(i).activeSensor = deblank(char(typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+obj.LABELSIZE))','uint16')));
              internalOffset = internalOffset + obj.TSLABELSIZE;
              tsInfo(i).refSensor = deblank(char(typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','uint16')));
              internalOffset = internalOffset + 8;
              internalOffset = internalOffset + 56;
              tsInfo(i).dLowCut = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
              internalOffset = internalOffset + 8;
              tsInfo(i).dHighCut = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
              internalOffset = internalOffset + 8;
              tsInfo(i).dSamplingRate = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
              internalOffset = internalOffset + 8;
              tsInfo(i).dResolution = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
              internalOffset = internalOffset + 8;
              tsInfo(i).bMark = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+2))','uint16');
              internalOffset = internalOffset + 2;
              tsInfo(i).bNotch = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+2))','uint16');
              internalOffset = internalOffset + 2;
              tsInfo(i).dEegOffset = typecast(tsPacket.data(offset+internalOffset:(offset+internalOffset-1+8))','double');
              offset = offset + 552;
              %disp([num2str(i) ' : ' TSInfo(i).label ' : ' TSInfo(i).activeSensor ' : ' TSInfo(i).refSensor ' : ' num2str(TSInfo(i).samplingRate)]);

          end
          obj.tsInfos{j} = tsInfo;
          end
      end

      % -- -- --

      %% Get Segment Start Times
      segmentIdx = Tags(find(strcmp({Tags.IDStr}, 'SegmentStream'),1)).index;
      indexIdx = find([Index.sectionIdx] == segmentIdx, 1);
      segmentInstance = Index(indexIdx);

      nrSegments = segmentInstance.sectionL/152;
      fseek(h, segmentInstance.offset,'bof');
      obj.segments = struct();
      for i = 1: nrSegments
        dateOLE = fread(h,1, 'double');
        obj.segments(i).dateOLE = dateOLE;
        unix_time = (dateOLE*(3600*24)) - 2209161600;% 2208988800; %8
        obj.segments(i).dateStr = datestr(unix_time/86400 + datenum(1970,1,1));
        datev = datevec( obj.segments(i).dateStr );
        obj.segments(i).startDate = datev(1:3);
        obj.segments(i).startTime = datev(4:6);
        fseek(h, 8 , 'cof'); %16
        obj.segments(i).duration = fread(h,1, 'double');%24
        fseek(h, 128 , 'cof'); %152
      end

      % Get nrValues per segment and channel
      for iSeg = 1:length(obj.segments)
        if iSeg > length(obj.tsInfos)
            obj.tsInfos{iSeg} = obj.tsInfos{iSeg-1}
        end

        % Add Channel Names to segments
        obj.segments(iSeg).chName = {obj.tsInfos{iSeg}.label};
        obj.segments(iSeg).refName = {obj.tsInfos{iSeg}.refSensor};
        obj.segments(iSeg).samplingRate = [obj.tsInfos{iSeg}.dSamplingRate];
        obj.segments(iSeg).scale = [obj.tsInfos{iSeg}.dResolution];

      end

      %% Get events  - Andrei Barborica, Dec 2015
      % Find sequence of events, that are stored in the section tagged 'Events'
      idxSection = find(strcmp('Events',{Tags.tag}));
      indexIdx = find([obj.index.sectionIdx] == obj.sections(idxSection).index);
      offset = obj.index(indexIdx).offset;

      ePktLen = 272;    % Event packet length, see EVENTPACKET definition
      eMrkLen = 240;    % Event marker length, see EVENTMARKER definition
      evtPktGUID = hex2dec({'80', 'F6', '99', 'B7', 'A4', '72', 'D3', '11', '93', 'D3', '00', '50', '04', '00', 'C1', '48'}); % GUID for event packet header
      HCEVENT_ANNOTATION = '{A5A95612-A7F8-11CF-831A-0800091B5BDA}';
      HCEVENT_SEIZURE    =  '{A5A95646-A7F8-11CF-831A-0800091B5BDA}';
      HCEVENT_FORMATCHANGE      =  '{08784382-C765-11D3-90CE-00104B6F4F70}';
      HCEVENT_PHOTIC            =  '{6FF394DA-D1B8-46DA-B78F-866C67CF02AF}';
      HCEVENT_POSTHYPERVENT     =  '{481DFC97-013C-4BC5-A203-871B0375A519}';
      HCEVENT_REVIEWPROGRESS    =  '{725798BF-CD1C-4909-B793-6C7864C27AB7}';
      HCEVENT_EXAMSTART         =  '{96315D79-5C24-4A65-B334-E31A95088D55}';
      HCEVENT_HYPERVENTILATION  =  '{A5A95608-A7F8-11CF-831A-0800091B5BDA}';
      HCEVENT_IMPEDANCE         =  '{A5A95617-A7F8-11CF-831A-0800091B5BDA}';
      HCEVENT_AMPLIFIERDISCONNECT = '{A71A6DB5-4150-48BF-B462-1C40521EBD6F}';
      HCEVENT_AMPLIFIERRECONNECT = '{6387C7C8-6F98-4886-9AF4-FA750ED300DE}';
      HCEVENT_PAUSED = '{71EECE80-EBC4-41C7-BF26-E56911426FB4}';

      DAYSECS = 86400.0;  % From nrvdate.h

      fseek(h,offset,'bof');
      pktGUID = fread(h,16,'uint8');
      pktLen  = fread(h,1,'uint64');
      obj.eventMarkers = struct();
      i = 0;    % Event counter
      while (pktGUID == evtPktGUID)
          i = i + 1;
          % Please refer to EVENTMARKER structure in the Nervus file documentation
          fseek(h,8,'cof'); % Skip eventID, not used
          evtDate           = fread(h,1,'double');
          evtDateFraction   = fread(h,1,'double');
          obj.eventMarkers(i).dateOLE = evtDate;
          obj.eventMarkers(i).dateFraction = evtDateFraction;
          evtPOSIXTime = evtDate*DAYSECS + evtDateFraction - 2209161600;% 2208988800; %8
          obj.eventMarkers(i).dateStr = datestr(evtPOSIXTime/DAYSECS + datenum(1970,1,1),'dd-mmmm-yyyy HH:MM:SS.FFF'); % Save fractions of seconds, as well
          obj.eventMarkers(i).duration  = fread(h,1,'double');
          fseek(h,48,'cof');
          evtUser                       = fread(h,12,'uint16');
          obj.eventMarkers(i).user      = deblank(char(evtUser).');
          evtTextLen                    = fread(h,1,'uint64');
          evtGUID                       = fread(h,16,'uint8');
          obj.eventMarkers(i).GUID      = sprintf('{%.2X%.2X%.2X%.2X-%.2X%.2X-%.2X%.2X-%.2X%.2X-%.2X%.2X%.2X%.2X%.2X%.2X}',evtGUID([4 3 2 1 6 5 8 7 9:16]));
          fseek(h,16,'cof');    % Skip Reserved4 array
          evtLabel                      = fread(h,32,'uint16'); % LABELSIZE = 32;
          evtLabel                      = deblank(char(evtLabel).');    % Not used
          eventMarkers(i).label         = evtLabel;

          %disp(sprintf('Offset: %s, TypeGUID:%s, User:%s, Label:%s',dec2hex(offset),evtGUID,evtUser,evtLabel));

          % Only a subset of all event types are dealt with
          switch obj.eventMarkers(i).GUID
              case HCEVENT_SEIZURE
                  obj.eventMarkers(i).IDStr = 'Seizure';
                  %disp(' Seizure event');
              case HCEVENT_ANNOTATION
                  obj.eventMarkers(i).IDStr = 'Annotation';
                  fseek(h,32,'cof');    % Skip Reserved5 array
                  evtAnnotation = fread(h,evtTextLen,'uint16');
                  obj.eventMarkers(i).annotation = deblank(char(evtAnnotation).');
                  %disp(sprintf(' Annotation:%s',evtAnnotation));
              case HCEVENT_FORMATCHANGE
                  obj.eventMarkers(i).IDStr = 'Format change';
              case HCEVENT_PHOTIC
                  obj.eventMarkers(i).IDStr = 'Photic';
              case HCEVENT_POSTHYPERVENT
                  obj.eventMarkers(i).IDStr = 'Posthyperventilation';
              case HCEVENT_REVIEWPROGRESS
                  obj.eventMarkers(i).IDStr = 'Review progress';
              case HCEVENT_EXAMSTART
                  obj.eventMarkers(i).IDStr = 'Exam start';
              case HCEVENT_HYPERVENTILATION
                  obj.eventMarkers(i).IDStr = 'Hyperventilation';
              case HCEVENT_IMPEDANCE
                  obj.eventMarkers(i).IDStr = 'Impedance';
              case HCEVENT_AMPLIFIERDISCONNECT
                  obj.eventMarkers(i).IDStr = 'Amplifier Disconnect';
              case HCEVENT_AMPLIFIERRECONNECT
                  obj.eventMarkers(i).IDStr = 'Amplifier Reconnect';
              case HCEVENT_PAUSED
                  obj.eventMarkers(i).IDStr = 'Recording Paused';
              otherwise
                  obj.eventMarkers(i).IDStr = 'UNKNOWN';
          end

          % Next packet
          offset = offset + pktLen;
          fseek(h,offset,'bof');
          pktGUID = fread(h,16,'uint8');
          pktLen  = fread(h,1,'uint64');
      end

      %% Get montage  - Andrei Barborica, Dec 2015
      % Derivation (montage)
      mtgIdx  = Tags(find(strcmp({Tags.IDStr},'DERIVATIONGUID'),1)).index;
      indexIdx      = find([obj.index.sectionIdx]==mtgIdx,1);
      fseek(h,obj.index(indexIdx(1)).offset + 40,'bof');    % Beginning of current montage name
      mtgName       = deblank(char(fread(h,32,'uint16')).');
      fseek(h,640,'cof');                             % Number of traces in the montage
      numDerivations = fread(h,1,'uint32');
      numDerivations2 = fread(h,1,'uint32');

      obj.montage = struct();
      for i = 1:numDerivations
          obj.montage(i).derivationName = deblank(char(fread(h,64,'uint16')).');
          obj.montage(i).signalName1    = deblank(char(fread(h,32,'uint16')).');
          obj.montage(i).signalName2    = deblank(char(fread(h,32,'uint16')).');
          fseek(h,264,'cof');         % Skip additional info
      end

      % Display properties
      dispIdx = Tags(find(strcmp({Tags.IDStr},'DISPLAYGUID'),1)).index;
      indexIdx  = find([obj.index.sectionIdx]==dispIdx,1);
      fseek(h,obj.index(indexIdx(1)).offset + 40,'bof');    % Beginning of current montage name
      displayName          = deblank(char(fread(h,32,'uint16')).');
      fseek(h,640,'cof');                             % Number of traces in the montage
      numTraces = fread(h,1,'uint32');
      numTraces2 = fread(h,1,'uint32');

      if (numTraces == numDerivations)
          for i = 1:numTraces
              fseek(h,32,'cof');
              obj.montage(i).color = fread(h,1,'uint32'); % Use typecast(uint32(montage(i).color),'uint8') to convert to RGB array
              fseek(h,136-4,'cof');
          end
      else
          disp('Could not match montage derivations with display color table');
      end

      % Close File
      fclose(h);

    end

    function out = getNrSamples(obj, segment)
      % GETNRSAMPLES  Returns the number of samples per channel in segment.
      %
      %   OUT = GETNRSAMPLES(OBJ, SEGMENT) returns a 1xn array of values
      %   indicating the number of samples for each of the channels in the
      %   associated SEGMENT, where SEGMENT is the index of the
      %   OBJ.segments array.

      assert(length(obj.segments)>= segment, ...
        'Incorrect SEGMENT argument; must be integer representing segment index.');

      out = obj.segments(segment).samplingRate .* obj.segments(segment).duration;

    end

    function cSumSegs = getCSumSegs(obj, chI)
      % GETCSUMSEGS  Returns the cumulative sum of a channels segments
      %
      %   CSUMSEGS = GETCSUMSEGS(OBJ, CHI) returns a 1xn array of the
      %   cumulative sum of the given channel, where chI is the channel
      %   index. Different TsInfos can have different numbers of channels
      %   and different sampling rates.
      samplingRates = zeros(1, length(obj.tsInfos));
      for ts_i = 1:length(obj.tsInfos)
          if length(obj.tsInfos{ts_i}) < chI
              continue
          else
              samplingRates(ts_i) = obj.tsInfos{ts_i}(chI).dSamplingRate;
          end
      end
      cSumSegs = [0 cumsum(samplingRates.*[obj.segments.duration])];
  end

    function out = getdata(obj, segment, range, chIdx)
      % GETDATA  Returns data from Nicolet file.
      %
      %   OUT = GETDATA(OBJ, SEGMENT, RANGE, CHIDX) returns data in an nxm array of
      %   doubles where n is the number of datapoints and m is the number
      %   of channels. RANGE is a 1x2 array with the [StartIndex EndIndex]
      %   and CHIDX is a vector of channel indeces.

      % Assert range is 1x2 vector
      if isempty([obj.tsInfos{segment}.label])
          warning('Segment %d has an empty tsInfo. Skipping', segment)
          out = 0;
          return
      end

      assert(length(range) == 2, 'Range is [firstIndex lastIndex]');
      assert(length(segment) == 1, 'Segment must be single value.');

      % Reopen .e file.
      h = fopen(obj.fileName,'r','ieee-le');

      % Find sectionID for channels
      lChIdx = length(chIdx);
      sectionIdx = zeros(lChIdx,1);
      for i = 1:lChIdx
        tmp = find(strcmp(num2str(chIdx(i)-1),{obj.sections.tag}),1);
        sectionIdx(i) = obj.sections(tmp).index;
      end

      % Iterate over all requested channels and populate array.
      out = zeros(range(2) - range(1) + 1, lChIdx);
      for i = 1 : lChIdx
        % Get cumulative sum segments.
        cSumSegments = obj.getCSumSegs(chIdx(i));

        % Get sampling rate for current channel
        curSF = obj.segments(segment).samplingRate(chIdx(i));
        mult = obj.segments(segment).scale(chIdx(i));

        % Find all sections
        allSectionIdx = obj.allIndexIDs == sectionIdx(i);
        allSections = find(allSectionIdx);

        % Find relevant sections
        sectionLengths = [obj.index(allSections).sectionL]./2;
        cSectionLengths = [0 cumsum(sectionLengths)];

        skipValues = cSumSegments(segment);
        firstSectionForSegment = find(cSectionLengths > skipValues, 1) - 1 ;
        lastSectionForSegment = firstSectionForSegment + ...
          find(cSectionLengths > curSF*obj.segments(segment).duration,1) - 2 ;

        if isempty(lastSectionForSegment)
          lastSectionForSegment = length(cSectionLengths);
        end

        offsetSectionLengths = cSectionLengths - cSectionLengths(firstSectionForSegment);

        firstSection = find(offsetSectionLengths < range(1) ,1,'last');
        lastSection = find(offsetSectionLengths >= range(2),1)-1;

        if isempty(lastSection)
          lastSection = length(offsetSectionLengths);
        end

        if lastSection > lastSectionForSegment
          error('Index out of range for current section: %i > %i, on channel: %i', ...
            range(2), cSectionLengths(lastSectionForSegment+1), chIdx(i));
        end

        useSections = allSections(firstSection: lastSection) ;
        useSectionL = sectionLengths(firstSection: lastSection) ;

        % First Partial Segment
        curIdx = 1;
        curSec = obj.index(useSections(1));
        fseek(h, curSec.offset,'bof');

        firstOffset = range(1) - offsetSectionLengths(firstSection);
        lastOffset = min([range(2) useSectionL(1)]);
        lsec = lastOffset-firstOffset + 1;

        fseek(h, (firstOffset-1) * 2,'cof');
        out(1 : lsec,i) = fread(h, lsec, 'int16') * mult;
        curIdx = curIdx +  lsec;

        if length(useSections) > 1
          % Full Segments
          for j = 2: (length(useSections)-1)
            curSec = obj.index(useSections(j));
            fseek(h, curSec.offset,'bof');

            out(curIdx : (curIdx + useSectionL(j) - 1),i) = ...
              fread(h, useSectionL(j), 'int16') * mult;
            curIdx = curIdx +  useSectionL(j);
          end

          % Final Partial Segment
          curSec = obj.index(useSections(end));
          fseek(h, curSec.offset,'bof');
          out(curIdx : end,i) = fread(h, length(out)-curIdx + 1, 'int16') * mult;
        end

      end

      % Close the .e file.
      fclose(h);

    end

     function out = getdataQ(obj, segment, range, chIdx)
      % GETDATAQ  Returns data from Nicolet file. This is a "QUICK" version of getdata,
      % that uses more memory but operates faster on large datasets by reading
      % a single block of data from disk that contains all data of interest.
      %
      %   OUT = GETDATAQ(OBJ, SEGMENT, RANGE, CHIDX) returns data in an nxm array of
      %   doubles where n is the number of datapoints and m is the number
      %   of channels. RANGE is a 1x2 array with the [StartIndex EndIndex]
      %   and CHIDX is a vector of channel indeces.
      %
      % Andrei Barborica, Dec 2015
      %
      if isempty([obj.tsInfos{segment}.label])
          warning('Segment %d has an empty tsInfo. Skipping', segment)
          out = 0;
          return
      end

      % Assert range is 1x2 vector
      assert(length(range) == 2, 'Range is [firstIndex lastIndex]');
      assert(length(segment) == 1, 'Segment must be single value.');

      % Reopen .e file.
      h = fopen(obj.fileName,'r','ieee-le');

      % Find sectionID for channels
      lChIdx = length(chIdx);
      sectionIdx = zeros(lChIdx,1);
      for i = 1:lChIdx
        tmp = find(strcmp(num2str(chIdx(i)-1),{obj.sections.tag}),1);
        sectionIdx(i) = obj.sections(tmp).index;
      end

      usedIndexEntries = zeros(size([obj.index.offset]));

      % Iterate over all requested channels and populate array.
      out = zeros(range(2) - range(1) + 1, lChIdx);
      for i = 1 : lChIdx
        % Get cumulative sum segments.
        cSumSegments = obj.getCSumSegs(chIdx(i));

        % Get sampling rate for current channel
        curSF = obj.segments(segment).samplingRate(chIdx(i));
        mult = obj.segments(segment).scale(chIdx(i));

        % Find all sections
        allSectionIdx = obj.allIndexIDs == sectionIdx(i);
        allSections = find(allSectionIdx);

        % Find relevant sections
        sectionLengths = [obj.index(allSections).sectionL]./2;
        cSectionLengths = [0 cumsum(sectionLengths)];

        skipValues = cSumSegments(segment);
        firstSectionForSegment = find(cSectionLengths > skipValues, 1) - 1 ;
        lastSectionForSegment = firstSectionForSegment + ...
          find(cSectionLengths > curSF*obj.segments(segment).duration,1) - 2 ;

        if isempty(lastSectionForSegment)
          lastSectionForSegment = length(cSectionLengths);
        end

        offsetSectionLengths = cSectionLengths - cSectionLengths(firstSectionForSegment);

        firstSection = find(offsetSectionLengths < range(1) ,1,'last');
        lastSection = find(offsetSectionLengths >= range(2),1)-1;

        if isempty(lastSection)
          lastSection = length(offsetSectionLengths);
        end

        if lastSection > lastSectionForSegment
          error('Index out of range for current section: %i > %i, on channel: %i', ...
            range(2), cSectionLengths(lastSectionForSegment+1), chIdx(i));
        end

        useSections = allSections(firstSection: lastSection) ;
        useSectionL = sectionLengths(firstSection: lastSection) ;

        % First Partial Segment
        usedIndexEntries(useSections(1)) = 1;

        if length(useSections) > 1
          % Full Segments
          for j = 2: (length(useSections)-1)
            usedIndexEntries(useSections(j)) = 1;
          end

          % Final Partial Segment
          usedIndexEntries(useSections(end)) = 1;
        end

      end

      % Read a big chunk of the file, containing data of interest.
      ix = find(usedIndexEntries);
      fseek(h, obj.index(ix(1)).offset,'bof');
      dsize =  obj.index(ix(end)).offset - obj.index(ix(1)).offset + obj.index(ix(end)).sectionL;
      tmp = fread(h,dsize/2,'int16').';

      % Close the .e file.
      fclose(h);

      baseOffset = obj.index(ix(1)).offset;

      % Extract specified channels
      for i = 1 : lChIdx
        % Get cumulative sum segments.
        cSumSegments = obj.getCSumSegs(chIdx(i));

        % Get sampling rate for current channel
        curSF = obj.segments(segment).samplingRate(chIdx(i));
        mult = obj.segments(segment).scale(chIdx(i));

        % Find all sections
        allSectionIdx = obj.allIndexIDs == sectionIdx(i);
        allSections = find(allSectionIdx);

        % Find relevant sections
        sectionLengths = [obj.index(allSections).sectionL]./2;
        cSectionLengths = [0 cumsum(sectionLengths)];

        skipValues = cSumSegments(segment);
        firstSectionForSegment = find(cSectionLengths > skipValues, 1) - 1 ;
        lastSectionForSegment = firstSectionForSegment + ...
          find(cSectionLengths > curSF*obj.segments(segment).duration,1) - 2 ;

        if isempty(lastSectionForSegment)
          lastSectionForSegment = length(cSectionLengths);
        end

        offsetSectionLengths = cSectionLengths - cSectionLengths(firstSectionForSegment);

        firstSection = find(offsetSectionLengths < range(1) ,1,'last');
        lastSection = find(offsetSectionLengths >= range(2),1)-1;

        if isempty(lastSection)
          lastSection = length(offsetSectionLengths);
        end

        if lastSection > lastSectionForSegment
          error('Index out of range for current section: %i > %i, on channel: %i', ...
            range(2), cSectionLengths(lastSectionForSegment+1), chIdx(i));
        end

        useSections = allSections(firstSection: lastSection) ;
        useSectionL = sectionLengths(firstSection: lastSection) ;

        % First Partial Segment
        curIdx = 1;
        curSec = obj.index(useSections(1));
        %fseek(h, curSec.offset,'bof');

        firstOffset = range(1) - offsetSectionLengths(firstSection);
        lastOffset = min([range(2) useSectionL(1)]);
        lsec = lastOffset-firstOffset + 1;

        out(1 : lsec,i) = tmp( (curSec.offset - baseOffset)/2 + (firstOffset-1) + (1:lsec) ) * mult;
        curIdx = curIdx +  lsec;

        if length(useSections) > 1
          % Full Segments
          for j = 2: (length(useSections)-1)
            curSec = obj.index(useSections(j));
            out(curIdx : (curIdx + useSectionL(j) - 1),i) = ...
                tmp( (curSec.offset - baseOffset)/2 + (1:useSectionL(j)) ) * mult;
            curIdx = curIdx +  useSectionL(j);
          end

          % Final Partial Segment
          curSec = obj.index(useSections(end));
          out(curIdx : end,i) = tmp( (curSec.offset - baseOffset)/2 + (1:(length(out)-curIdx + 1)) ) * mult; % length(out) ??????
        end

      end
    end

    function labels = getlabels(obj,str)
        % Returns annotations containing specified string
        %
        % Cristian Donos, Dec 2015
        %
        labels=[]; counter = 1;
        for i = 1:length(obj.eventMarkers)
            if strfind(lower(obj.eventMarkers(i).annotation),lower(str))
            labels{counter,1} = obj.eventMarkers(i).annotation;  % annotation string
            labels{counter,2} = i;  % annotation index in obj.eventMarkers
            % identify segment
            time_vector = [];
            for j = 1:length(obj.segments)
                time_vector = [time_vector etime(datevec(obj.eventMarkers(i).dateStr),datevec(obj.segments(j).dateStr))];
            end
            labels{counter,3}= find(time_vector==min(time_vector(time_vector>0)));  % annotation part of this segment
            labels{counter,4}= min(time_vector(time_vector>0));  % annotation offset in seconds, relative to its segment start
            counter = counter+1;
            end
        end
    end
  end

end

author: brainstorm-tools


text:
function y = sphere_characteristic ( m, n, x )

%*****************************************************************************80
%
%% SPHERE_CHARACTERISTIC evaluates the characteristic function of a sphere.
%
%  Discussion:
%
%    The dimension of the sphere is arbitrary, and set by the input value M.
%
%    The sphere is assumed to have center at the origin and radius 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points to check.
%
%    Input, real X(M,N), the coordinates of points to be checked.
%
%    Output, real Y(N,1), is 1 if the point is inside, 0 otherwise.
%
  y = sqrt ( sum ( x.^2, 1 ) );

  y = y';

  y = ( y < 1.0 );

  return
end


text:
function datetickzoom(varargin)
%DATETICKZOOM Date formatted tick labels, automatically updated when zoomed or panned.
%   Arguments are completely identical to does of DATETICK. The argument
%   DATEFORM is reset once zoomed or panned.
%
%   See also datetick, datestr, datenum

if nargin==2 && isstruct(varargin{2}) && isfield(varargin{2},'Axes') && isscalar(varargin{2}.Axes)
    datetickdata = getappdata(varargin{2}.Axes,'datetickdata');
    if isstruct(datetickdata) && isfield(datetickdata,'axh') && datetickdata.axh==varargin{2}.Axes
        axh = datetickdata.axh;
        ax = datetickdata.ax;
        dateform = datetickdata.dateform;
        keep_ticks = datetickdata.keep_ticks;
        if keep_ticks
            set(axh,[ax,'TickMode'],'auto')
            if ~isempty(dateform)
                datetick(axh,ax,dateform,'keeplimits','keepticks')
            else
                datetick(axh,ax,'keeplimits','keepticks')
            end
        else
            if ~isempty(dateform)
                datetick(axh,ax,dateform,'keeplimits')
            else
                datetick(axh,ax,'keeplimits')
            end
        end
    end
else
    [axh,ax,dateform,keep_ticks] = parseinputs(varargin);
    datetickdata = [];
    datetickdata.axh = axh;
    datetickdata.ax = ax;
    datetickdata.dateform = dateform;
    datetickdata.keep_ticks = keep_ticks;

    setappdata(axh,'datetickdata',datetickdata);
    set(zoom(axh),'ActionPostCallback',@datetickzoom)
    set(pan(get(axh,'parent')),'ActionPostCallback',@datetickzoom)
    datetick(varargin{:})
end

function [axh,ax,dateform,keep_ticks] = parseinputs(v)
%Parse Inputs

% Defaults;
nin = length(v);
dateform = [];
keep_ticks = 0;

% check to see if an axes was specified
if nin > 0 & ishandle(v{1}) & isequal(get(v{1},'type'),'axes') %#ok ishandle return is not scalar
    % use the axes passed in
    axh = v{1};
    v(1)=[];
    nin=nin-1;
else
    % use gca
    axh = gca;
end
% Look for 'keepticks'
for i=nin:-1:max(1,nin-1),
   if strcmpi(v{i},'keepticks'),
      keep_ticks = 1;
      v(i) = [];
      nin = nin-1;
   end
end

if nin==0,
   ax = 'x';
else
   ax = v{1};
end
if nin > 1
     % The dateform (Date Format) value should be a scalar or string constant
     % check this out
     dateform = v{2};
     if (isnumeric(dateform) && length(dateform) ~= 1) && ~ischar(dateform)
         error('The Date Format value should be a scalar or string');
     end
end


text:
function mode = setColormap(mode,cmapFn)
%
%  mode = setColormap(mode,cmapFn)
%
% Author:  My guess, Dave Heeger
% Purpose:
%   If you call a color map function, say hsv,
%   this routine builds a color map entry for that mode using the hsv function.
%

numGrays  = mode.numGrays;
numColors = mode.numColors;

% mode.cmap=cmapFn(numGrays,numColors);
eval(['mode.cmap=',cmapFn,'(numGrays,numColors);']);
mode.name = cmapFn;

return


text:
% ebcdic2ascii : Converts EBCDIC formatted text to ASCII formatted text
%
% CALL : ascii=ebcdic2ascii(ebcdic);
%
% ebcdic : Array on unsigned integers
% ascii  : Array on unsigned integers
%
% (C) 2002-2004, Thomas Mejer Hansen, tmh@gfy.ku.dk/thomas@cultpenguin.com
%

%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
function ascii=ebcdic2ascii(ebcdic)

  ebcdic_arr=[
    '00';'01';'02';'03';'04';'05';'06';'07';'08';'09';
    '0A';'0B';'0C';'0D';'0E';'0F';'10';'11';'12';'13';
    '14';'15';'16';'17';'18';'19';'1A';'1B';'1C';'1D';
    '1E';'1F';'20';'21';'22';'23';'24';'25';'26';'27';
    '28';'29';'2A';'2B';'2C';'2D';'2E';'2F';'30';'31';
    '32';'33';'34';'35';'36';'37';'38';'39';'3A';'3B';
    '3C';'3D';'3E';'3F';'40';'4A';'4B';'4C';'4D';'4E';
    '4F';'50';'5A';'5B';'5C';'5D';'5E';'5F';'60';'61';
    '6A';'6B';'6C';'6D';'6E';'6F';'79';'7A';'7B';'7C';
    '7D';'7E';'7F';'81';'82';'83';'84';'85';'86';'87';
    '88';'89';'91';'92';'93';'94';'95';'96';'97';'98';
    '99';'A1';'A2';'A3';'A4';'A5';'A6';'A7';'A8';'A9';
    'C0';'C1';'C2';'C3';'C4';'C5';'C6';'C7';'C8';'C9';
    'D0';'D1';'D2';'D3';'D4';'D5';'D6';'D7';'D8';'D9';
    'E0';'E2';'E3';'E4';'E5';'E6';'E7';'E8';'E9';'F0';
    'F1';'F2';'F3';'F4';'F5';'F6';'F7';'F8';'F9';'FF';
    '00';'01';'02';'03';'37';'2D';'2E';'2F';'2F';'16';
    '05';'25';'0B';'0C';'0D';'10';'11';'12';'13';'3C';
    '3D';'32';'26';'18';'3F';'27';'1C';'1D';'1D';'1E';
    '1F';'07';'40';'5A';'7F';'7B';'5B';'6C';'50';'7D';
    '4D';'5D';'5C';'4E';'6B';'60';'4B';'61';'F0';'F1';
    'F2';'F3';'F4';'F5';'F6';'F7';'F8';'F9';'7A';'5E';
    '4C';'7E';'6E';'6F';'7C';'00';'E0';'00';'00';'6D';
    '79';'C0';'4F';'D0';'A1'];

ascii_arr=[
    '00';'01';'02';'03';'9C';'09';'86';'7F';'97';'8D';
    '8E';'0B';'0C';'0D';'0E';'0F';'10';'11';'12';'13';
    '9D';'85';'08';'87';'18';'19';'92';'8F';'1C';'1D';
    '1E';'1F';'80';'81';'82';'83';'84';'0A';'17';'1B';
    '88';'89';'8A';'8B';'8C';'05';'06';'07';'90';'91';
    '16';'93';'94';'95';'96';'04';'98';'99';'9A';'9B';
    '14';'15';'9E';'1A';'20';'A2';'2E';'3C';'28';'2B';
    '7C';'26';'21';'24';'2A';'29';'3B';'AC';'2D';'2F';
    'A6';'2C';'25';'5F';'3E';'3F';'60';'3A';'23';'40';
    '27';'3D';'22';'61';'62';'63';'64';'65';'66';'67';
    '68';'69';'6A';'6B';'6C';'6D';'6E';'6F';'70';'71';
    '72';'7E';'73';'74';'75';'76';'77';'78';'79';'7A';
    '7B';'41';'42';'43';'44';'45';'46';'47';'48';'49';
    '7D';'4A';'4B';'4C';'4D';'4E';'4F';'50';'51';'52';
    '5C';'53';'54';'55';'56';'57';'58';'59';'5A';'30';
    '31';'32';'33';'34';'35';'36';'37';'38';'39';'9F';
    '00';'01';'02';'03';'04';'05';'06';'07';'07';'08';
    '09';'0A';'0B';'0C';'0D';'10';'11';'12';'13';'14';
    '15';'16';'17';'18';'1A';'1B';'1C';'1D';'1D';'1E';
    '1F';'7F';'20';'21';'22';'23';'24';'25';'26';'27';
    '28';'29';'2A';'2B';'2C';'2D';'2E';'2F';'30';'31';
    '32';'33';'34';'35';'36';'37';'38';'39';'3A';'3B';
    '3C';'3D';'3E';'3F';'40';'5B';'5C';'5D';'5E';'5F';
    '60';'7B';'7C';'7D';'7E'
];

ascii_dec_array=hex2dec(ascii_arr);

ebcdic_dec_array=hex2dec(ebcdic_arr);

for i=1:length(ebcdic)
  m=find(ebcdic_dec_array==ebcdic(i));
  if length(m)>1,
    m=m(1);
  end
  if length(m)==0,
    ascii(i)=0;
  else
    ascii(i)=(ascii_dec_array(m));
  end
end

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function [metabolite_structure,IDsAdded] = getCas2CTD(metabolite_structure)
%
%
% The input file was obtained from http://ctdbase.org/reports/CTD_chemicals.csv.gz
% 1st col: ctd id, 3rd col cas

% INPUT
% metabolite_structure  metabolite structure
%
% OUTPUT
% metabolite_structure  Updated metabolite structure
%
%
% Ines Thiele, 2020-2021

[NUM,TXT,RAW]=xlsread('CTD_chemicals.xlsx');
annotationSource = 'CTD_id mapping from CasRegistry';
annotationType = 'automatic';

Mets = fieldnames(metabolite_structure);
for i = 1 : length(Mets)
    Cas{i,1} = Mets{i};
    Cas{i,2} = num2str(metabolite_structure.(Mets{i}).casRegistry);
end

% listing starts at row 28
b =1;
for i = 28 : size(RAW,1)
    casRAW{b,1} = num2str(RAW{i,3});
    casRAWList{b,1} = RAW{i,2};
    b = b+1;
end
a = 1;
IDsAdded = '';

for i = 1 : size(Cas,1)

    % no keggId exists
    if ~isempty(metabolite_structure.(Cas{i,1}).casRegistry) && isempty(find(isnan(metabolite_structure.(Cas{i,1}).casRegistry),1))
        if isempty(metabolite_structure.(Cas{i,1}).ctd) || ~isempty(find(isnan(metabolite_structure.(Cas{i,1}).ctd),1))
            match = strmatch(Cas{i,2},casRAW,'exact');
            if ~isempty(match)
                metabolite_structure.(Cas{i,1}).ctd = casRAWList{match};
                metabolite_structure.(Cas{i,1}).ctd_source =  [annotationSource,':',annotationType,':',datestr(now)];
                IDsAdded{a,1} = Cas{i,1};
                IDsAdded{a,2} = 'ctd';
                IDsAdded{a,3} =  metabolite_structure.(Cas{i}).ctd;
                a = a + 1;
            end
        end
    end
end
[metabolite_structure] = cleanUpMetabolite_structure(metabolite_structure);

author: opencobra

repo: cobratoolbox


text:
function [atm] = mbar2atm(mbar)
% Convert pressure from millibars to atmospheres.
% Chad Greene 2012
atm = mbar*0.000986923;


text:
function r = reval(zz, zj, fj, wj)
%REVAL   Evaluate rational function in barycentric form.
%   R = REVAL(ZZ, ZJ, FJ, WJ) returns R (vector of floats), the values of the
%   barycentric rational function with support points ZJ, function values FJ,
%   and barycentric weights WJ evaluated at the points ZZ.
%
% See also AAA, MINIMAX, PRZ.

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

zv = zz(:);                             % vectorize zz if necessary
CC = 1./bsxfun(@minus, zv, zj.');       % Cauchy matrix
r = (CC*(wj.*fj))./(CC*wj);             % vector of values

% Deal with input inf: r(inf) = lim r(zz) = sum(w.*f) / sum(w):
r(isinf(zv)) = sum(wj.*fj)./sum(wj);

% Deal with NaN:
ii = find(isnan(r));
for jj = 1:length(ii)
    if ( isnan(zv(ii(jj))) || ~any(zv(ii(jj)) == zj) )
        % r(NaN) = NaN is fine.
        % The second case may happen if r(zv(ii)) = 0/0 at some point.
    else
        % Clean up values NaN = inf/inf at support points.
        % Find the corresponding node and set entry to correct value:
        r(ii(jj)) = fj(zv(ii(jj)) == zj);
    end
end

% Reshape to input format:
r = reshape(r, size(zz));

end % End of REVAL().

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test_issue2169

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_apply_montage

%%

data0 = {};
data0.label = {
  '1'
  '2'
  '3'
  '4'
  '5'
  '6'
  '7'
  '8'
  '9'
  };
data0.trial{1} = [
  1 1 1 1 1 1 1 1 1 1
  2 2 2 2 2 2 2 2 2 2
  3 3 3 3 3 3 3 3 3 3
  4 4 4 4 4 4 4 4 4 4
  5 5 5 5 5 5 5 5 5 5
  6 6 6 6 6 6 6 6 6 6
  7 7 7 7 7 7 7 7 7 7
  8 8 8 8 8 8 8 8 8 8
  9 9 9 9 9 9 9 9 9 9
  ];
data0.time{1} = 1:10;

data1 = data0;
data1.trial{1}(1,10) = nan;

data2 = data0;
data2.trial{1}(2,10) = nan;

data3 = data0;
data3.trial{1}(3,10) = nan;

%%

cfg = [];
cfg.montage.labelold = {'1', '2', '3'};
cfg.montage.labelnew = {'1', '2', '1-2'};
cfg.montage.tra = [
  1  0 0
  0  1 0
  1 -1 0
  ];
reref0 = ft_preprocessing(cfg, data0);
reref1 = ft_preprocessing(cfg, data1);
reref2 = ft_preprocessing(cfg, data2);
reref3 = ft_preprocessing(cfg, data3);

% check that the nans don't spread

assert(~isnan(reref0.trial{1}(1,10)))
assert(~isnan(reref0.trial{1}(2,10)))
assert(~isnan(reref0.trial{1}(3,10)))

assert( isnan(reref1.trial{1}(1,10)))
assert(~isnan(reref1.trial{1}(2,10)))
assert( isnan(reref1.trial{1}(3,10)))

assert(~isnan(reref2.trial{1}(1,10)))
assert( isnan(reref2.trial{1}(2,10)))
assert( isnan(reref2.trial{1}(3,10)))

assert(~isnan(reref3.trial{1}(1,10)))
assert(~isnan(reref3.trial{1}(2,10)))
assert(~isnan(reref3.trial{1}(3,10)))

%%
% repeat the same thing as above, but now for single precision data

data0.trial{1} = single(data0.trial{1});
data1.trial{1} = single(data1.trial{1});
data2.trial{1} = single(data2.trial{1});
data3.trial{1} = single(data3.trial{1});

cfg = [];
cfg.montage.labelold = {'1', '2', '3'};
cfg.montage.labelnew = {'1', '2', '1-2'};
cfg.montage.tra = [
  1  0 0
  0  1 0
  1 -1 0
  ];
reref0 = ft_preprocessing(cfg, data0);
reref1 = ft_preprocessing(cfg, data1);
reref2 = ft_preprocessing(cfg, data2);
reref3 = ft_preprocessing(cfg, data3);

% check that the nans don't spread

assert(~isnan(reref0.trial{1}(1,10)))
assert(~isnan(reref0.trial{1}(2,10)))
assert(~isnan(reref0.trial{1}(3,10)))

assert( isnan(reref1.trial{1}(1,10)))
assert(~isnan(reref1.trial{1}(2,10)))
assert( isnan(reref1.trial{1}(3,10)))

assert(~isnan(reref2.trial{1}(1,10)))
assert( isnan(reref2.trial{1}(2,10)))
assert( isnan(reref2.trial{1}(3,10)))

assert(~isnan(reref3.trial{1}(1,10)))
assert(~isnan(reref3.trial{1}(2,10)))
assert(~isnan(reref3.trial{1}(3,10)))

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function result = substr(s, i, j)
%SUBSTR Extracts a substring using s(i:j) when last > 0, or else s(i:end+j).

    if j > 0
        result = s(i:j);
    else
        result = s(i:end+j);
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function obj = loadobj(s)
% LOADOBJ  Initialize a DagNN object from a structure.
%   OBJ = LOADOBJ(S) initializes a DagNN objet from the structure
%   S. It is the opposite of S = OBJ.SAVEOBJ().
%   If S is a string, initializes the DagNN object with data
%   from a mat-file S. Otherwise, if S is an instance of `dagnn.DagNN`,
%   returns S.

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if ischar(s) s = load(s); end
if isstruct(s)
  assert(isfield(s, 'layers'), 'Invalid model.');
  if ~isstruct(s.layers)
    warning('The model appears to be `simplenn` model. Using `fromSimpleNN` instead.');
    obj = dagnn.DagNN.fromSimpleNN(s);
    return;
  end
  obj = dagnn.DagNN() ;
  for l = 1:numel(s.layers)
    constr = str2func(s.layers(l).type) ;
    block = constr() ;
    block.load(struct(s.layers(l).block)) ;
    obj.addLayer(...
      s.layers(l).name, ...
      block, ...
      s.layers(l).inputs, ...
      s.layers(l).outputs, ...
      s.layers(l).params,...
      'skipRebuild', true) ;
  end
  obj.rebuild();
  if isfield(s, 'params')
    for f = setdiff(fieldnames(s.params)','name')
      f = char(f) ;
      for i = 1:numel(s.params)
        p = obj.getParamIndex(s.params(i).name) ;
        obj.params(p).(f) = s.params(i).(f) ;
      end
    end
  end
  if isfield(s, 'vars')
    for f = setdiff(fieldnames(s.vars)','name')
      f = char(f) ;
      for i = 1:numel(s.vars)
        p = obj.getVarIndex(s.vars(i).name) ;
        obj.vars(p).(f) = s.vars(i).(f) ;
      end
    end
  end
  for f = setdiff(fieldnames(s)', {'vars','params','layers'})
    f = char(f) ;
    obj.(f) = s.(f) ;
  end
elseif isa(s, 'dagnn.DagNN')
  obj = s ;
else
  error('Unknown data type %s for `loadobj`.', class(s));
end

repo: LapSRN


text:
function result = isNaNgeneral(x)
% Determine if x contains a NaN value
%
% function result = isNaNgeneral(x)
%
% Returns a logical value which indicates whether or not the input x
% contains a NaN value. The input x can be defined recursively by arrays,
% structs and cells.
%
% See also:

% This file is part of Manopt: www.manopt.org.
% Original author: Xiaowen Jiang, Aug. 31, 2021.
% Contributors: Nicolas Boumal
% Change log:

    if ~isstruct(x) && ~iscell(x) && ~isnumeric(x)
        up = MException('manopt:isNaNgeneral', ...
                    'isNaNgeneral should only accept structs, cells or arrays.');
        throw(up);
    end

    % recursively find NaN for each part of x
    if isstruct(x)
        result = isNaN_struct(x);
        if result > 0
            result = true;
        end
    elseif iscell(x)
        result = isNaN_cell(x);
        if result > 0
            result = true;
        end
    else
        result = any(isnan(x(:)));
    end

    % when x is a struct
    function result = isNaN_struct(x)
        elems = fieldnames(x);
        nelems = numel(elems);
        result = false;
        for ii = 1:nelems
            if isstruct(x.(elems{ii}))
                result = result + isNaN_struct(x.(elems{ii}));
            elseif iscell(x.(elems{ii}))
                result = result + isNaN_cell(x.(elems{ii}));
            else
                result = result + any(isnan(x.(elems{ii})(:)));
            end
        end
    end

    % when x is a cell
    function result = isNaN_cell(x)
        ncell = length(x);
        result = false;
        for ii = 1:ncell
            if isstruct(x{ii})
                result = result + isNaN_struct(x{ii});
            elseif iscell(x{ii})
                result = result + isNaN_cell(x{ii});
            else
                result = result + any(isnan(x{ii}(:)));
            end
        end
    end
end


text:
function c = times(a,b)
% implements e .* c and c .* e
%
% Description
% If |e| is a matrix of embeddings and |c| is a matrix of coefficients
% then |e .* c| is again a matrix of embeddings defined by
%
% $$ [\mathrm{e .* c}]_{j\ell} = mathrm{e}_{j\ell} \mathrm{c}_{j \ell}$$
%
% Syntax
%   out = e .* c
%   out = c .* e
%   out = e1 .* e2
%
% Input
%  e, e1, e2 - @embedding
%  c - double
%
% Output
%  out- @embedding
%

if isa(a,'embedding')

  if isa(b,'embedding')

    for i = 1:length(a.u), a.u{i} = a.u{i} .* b.u{i}; end

  else
    for i = 1:length(a.u), a.u{i} = a.u{i} .* b; end
  end
  c = a;

else

  for i = 1:length(b.u), b.u{i} = a .* b.u{i}; end
  c = b;

end

end

author: mtex-toolbox

repo: mtex

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


text:
function fem = genNedbasFEM3D(mesh,type)

%% Usage: Form Conforming P1 FEM degrees of freedom structure
%
% INPUTS:
% mesh --- a struct data contains very rich mesh information.
%
% OUTPUTS:
% fem --- a struct data contains the following fields:
%         fem.p: (x,y,z) coordinate of each vertex w.r.t a global DoF
%         fem.t: indices of global DoF in each element
%         fem.type: type of FEM
%         fem.ldof: number of local DoF on each element
%         fem.bas: three dimension matrix structure contains
%                  basis function on each element

% Last Modified: 07/02/2020 by Xu Zhang

%% 1. Form p,t
p = mesh.p; t = mesh.t;

%% 2. Form basis in vector form
bas = bas3DP1(p,t);
if type == 1
    bas = bas3D_ned1(bas);
end
fem = struct('p',p, 't',t, 'type','P1', 'ldof',6, 'bas',bas);

repo: ifem


text:
function tests = test_ft_preproc_padding

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

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;
dat     = randn(nchan, nsample) + 1;

result = [];
result{end+1} = ft_preproc_padding(dat, 'zero',      10);
result{end+1} = ft_preproc_padding(dat, 'mean',      10);
result{end+1} = ft_preproc_padding(dat, 'localmean', 10);
result{end+1} = ft_preproc_padding(dat, 'edge',      10);
result{end+1} = ft_preproc_padding(dat, 'mirror',    10);
result{end+1} = ft_preproc_padding(dat, 'nan',       10);
result{end+1} = ft_preproc_padding(dat, 'remove',    10);
result{end+1} = ft_preproc_padding(dat, 'zero',      20);
result{end+1} = ft_preproc_padding(dat, 'mean',      20);
result{end+1} = ft_preproc_padding(dat, 'localmean', 20);
result{end+1} = ft_preproc_padding(dat, 'edge',      20);
result{end+1} = ft_preproc_padding(dat, 'mirror',    20);
result{end+1} = ft_preproc_padding(dat, 'nan',       20);
result{end+1} = ft_preproc_padding(dat, 'remove',    20);
result{end+1} = ft_preproc_padding(dat, 'zero',      10, 20);
result{end+1} = ft_preproc_padding(dat, 'mean',      10, 20);
result{end+1} = ft_preproc_padding(dat, 'localmean', 10, 20);
result{end+1} = ft_preproc_padding(dat, 'edge',      10, 20);
result{end+1} = ft_preproc_padding(dat, 'mirror',    10, 20);
result{end+1} = ft_preproc_padding(dat, 'nan',       10, 20);
result{end+1} = ft_preproc_padding(dat, 'remove',    10, 20);

for i=1:numel(result)
  % the number of samples should be larger, or smaller upon 'remove'
  assert(size(result{i},2)~=nsample);
end

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end