text:
function model = annotateModel(model, annotateRxns, annotateMets,modelID,modelName,modelAnnotation)
% This function annotates a model with VMH reaction and metabolite
% identifiers.
%
% function model = annotateModel(model, annotateRxns,annotateMets)
%
% INPUT
% model         Model structure
% annotateRxns  default: 1
% annotateMets  default: 1
%
% Optional:
% modelID       ID of model
% modelName     model Name
%
% OUTPUT
% model         Updated model structure
%
% Ines Thiele October 2019

if ~exist('annotateRxns','var')
    annotateRxns = 1;
end

if ~exist('annotateMets','var')
    annotateMets = 1;
end

if exist('modelID','var')
    model.modelID = modelID;
end

if exist('modelID','var')
    model.modelName = modelName;
end
if exist('modelAnnotation','var')
    model.modelAnnotation = modelAnnotation;
end
%% read input files
if annotateRxns == 1
    % load reactions from VMH
    VMH_rxns = tdfread('VMH_Reactions.tsv');
end

if annotateMets == 1
    % load metabolites from VMH
    VMH_mets = tdfread('VMH_Metabolites.tsv');
end
%% reactions
if annotateRxns == 1
    % redefine fields that should be updated
    %model.rxnNames = {};
    %model.rxnNames(1:length(model.rxns),1) = {''};
    model.rxnECNumbers = {};
    model.rxnECNumbers(1:length(model.rxns),1) = {''};
    %model.subSystems = {};
    %model.subSystems(1:length(model.rxns),1) = {''};
    model.rxnMetaNetXID = {};
    model.rxnMetaNetXID(1:length(model.rxns),1) = {''};
    model.rxnKEGGID = {};
    model.rxnKEGGID(1:length(model.rxns),1) = {''};
    model.rxnSBOTerms = {};
    model.rxnSBOTerms(1:length(model.rxns),1) = {''};
    model.rxnSEEDID = {};
    model.rxnSEEDID(1:length(model.rxns),1) = {''};

    VMH_rxns2.abbreviation = strcat('_',VMH_rxns.abbreviation);
    [a,remDash] = strtok(model.rxns,'_');
    [remDashBrak,rem2] = strtok(remDash,'['); % remove [
    [remBrak] = strtok(model.rxns,'['); % remove [
    [remDBrak,rem2] = strtok(model.rxns,'_['); % remove [
    for i = 1 : length(VMH_rxns.abbreviation)
        m =[];
        % check for organ in front and [
        tokVMH = strtok(VMH_rxns2.abbreviation(i,:),'['); % remove [
        tokVMH = strtok(tokVMH,' '); % remove space
        m = strmatch(tokVMH,remDash,'exact');
        if isempty(m)
            % no organ in front but [
            tokVMH = strtok(VMH_rxns.abbreviation(i,:),'['); % remove [
            m = strmatch(tokVMH,remBrak,'exact');
            if isempty(m)
                m = strmatch(tokVMH,remDBrak,'exact');
            end
        end
        if ~isempty(m)
            model.rxnECNumbers(m,1) = cellstr(VMH_rxns.ecnumber(i,:));
            model.rxnNames(m,1) = cellstr(VMH_rxns.description(i,:));
            model.subSystems(m,1) = cellstr(regexprep(VMH_rxns.subsystem(i,:),'"',''));
            model.rxnMetaNetXID(m,1) = cellstr(VMH_rxns.metanetx(i,:));
            model.rxnKEGGID(m,1) = cellstr(VMH_rxns.keggId(i,:));
            model.rxnSEEDID(m,1) = cellstr(VMH_rxns.seed(i,:));
            % annotate reaction SBO terms:
            if  ~isempty(strfind(VMH_rxns.abbreviation(i,:),'biomass')) ||  ~isempty(strfind(VMH_rxns.abbreviation(i,:),'Biomass'))
                % biomass reaction
                model.rxnSBOTerms(m,1) = cellstr('SBO:0000629');
                % Exchange reaction
            elseif ~isempty(strmatch('EX_',VMH_rxns.abbreviation(i,:))) ||  ~isempty(strmatch('Ex_',VMH_rxns.abbreviation(i,:)))
                model.rxnSBOTerms(m,1) = cellstr('SBO:0000627');
            elseif ~isempty(strmatch('DM_',VMH_rxns.abbreviation(i,:)))
                model.rxnSBOTerms(m,1) = cellstr('SBO:0000628');
            elseif  ~isempty(strmatch('sink_',VMH_rxns.abbreviation(i,:)))
                model.rxnSBOTerms(m,1) = cellstr('SBO:0000632');
            elseif  ~isempty(strmatch('Transport',regexprep(VMH_rxns.subsystem(i,:),'"','')))
                model.rxnSBOTerms(m,1) = cellstr('SBO:0000185');
            elseif  ~isempty(strmatch('Exchange',VMH_rxns.subsystem(i,:)))
                model.rxnSBOTerms(m,1) = cellstr('SBO:0000628');
            else % metabolicx reaction
                model.rxnSBOTerms(m,1) = cellstr('SBO:0000176');
            end
        end

    end
end
% check for biomass reaction separately again in case that the named one
% was not in the vmh file
BM1 = strmatch('biomass',model.rxns);
BM2 = strmatch('Biomass',model.rxns);
BM3 = strmatch('bio',model.rxns);
BM = unique([BM1;BM2;BM3]);
for i = 1 : length(BM)
    model.rxnSBOTerms(BM(i),1) = cellstr('SBO:0000629');
end

Sink = strmatch('sink_',model.rxns);
for i = 1 : length(Sink)
    model.rxnSBOTerms(Sink(i),1) = cellstr('SBO:0000632');
end

% annotate missing transport reactions based on subSystem with SBO term
for i = 1:length(model.rxns)
    if ~isempty(strmatch('Transport',model.subSystems(i)))
        model.rxnSBOTerms(i,1) = cellstr('SBO:0000185');
    elseif ~isempty(find(~cellfun(@isempty,strfind(model.rxns,'EX_'))))
        model.rxnSBOTerms(i,1) = cellstr('SBO:0000628');
    elseif ~isempty(find(~cellfun(@isempty,strfind(model.rxns,'DM_'))))
        model.rxnSBOTerms(i,1) = cellstr('SBO:0000627');
    elseif  ~isempty(strmatch('Demand',model.subSystems(i)))
        model.rxnSBOTerms(i,1) = cellstr('SBO:0000185');
    end
end

%% Metabolites
if annotateMets == 1
    model.metNames= {};
    model.metNames(1:length(model.mets),1) = {''};
    model.metFormulas= {};
    model.metFormulas(1:length(model.mets),1) = {''};
    model.metCharges = [];
    model.metCharges(1:length(model.mets),1) = 0;
    model.metChEBIID={};
    model.metChEBIID(1:length(model.mets),1) = {''};
    model.metHMDBID= {};
    model.metHMDBID(1:length(model.mets),1) = {''};
    model.metInChIString= {};
    model.metInChIString(1:length(model.mets),1) = {''};
    model.metKEGGID= {};
    model.metKEGGID(1:length(model.mets),1) = {''};
    model.metSmiles= {};
    model.metSmiles(1:length(model.mets),1) = {''};
    model.metMetaNetXID= {};
    model.metMetaNetXID(1:length(model.mets),1) = {''};
    model.metPubChemID= {};
    model.metPubChemID(1:length(model.mets),1) = {''};
    model.metBiGGID= {};
    model.metBiGGID(1:length(model.mets),1) = {''};
    model.metBioCycID= {};
    model.metBioCycID(1:length(model.mets),1) = {''};
    model.metSEEDID= {};
    model.metSEEDID(1:length(model.mets),1) = {''};
    model.metSBOTerms= {};
    model.metSBOTerms(1:length(model.mets),1) = {''};
    model.metChemSpider= {};
    model.metChemSpider(1:length(model.mets),1) = {''};
    model.metInchiKey= {};
    model.metInchiKey(1:length(model.mets),1) = {''};

    % note that there are a few metabolites that are not present
    [a,rem] = strtok(model.mets,'_');
    [rem,rem2] = strtok(a,'[');

    % add _ to Recon3D.mets
    VMH_mets2 = VMH_mets;
    % VMH_mets2.abbreviation = strcat('_',VMH_mets.abbreviation);
    % match both metabolite lists
    [LIA,LOCB] = ismember(rem,VMH_mets2.abbreviation);
    for i = 1 : length(LOCB)
        if LOCB(i) >0 && ~strcmp(a(i),'slack')
            model.metCharges(i,1) = str2num(VMH_mets2.charge(LOCB(i),:));
            model.metFormulas(i,1) = cellstr(VMH_mets2.chargedFormula(LOCB(i),:));
            model.metNames(i,1) =  cellstr(VMH_mets2.fullName(LOCB(i),:));
            if ~isempty(regexprep(VMH_mets2.cheBlId(LOCB(i),:),' ',''))% remove space in string
                model.metChEBIID(i,1) =  cellstr(strcat('CHEBI:',num2str(VMH_mets2.cheBlId(LOCB(i),:))));
            end
            model.metHMDBID(i,1) =  cellstr(VMH_mets2.hmdb(LOCB(i),:));
            model.metInChIString(i,1) =  cellstr(regexprep(VMH_mets2.inchiString(LOCB(i),:),'"',''));
            model.metKEGGID(i,1) =  cellstr(VMH_mets2.keggId(LOCB(i),:));
            model.metSmiles(i,1) =  cellstr(VMH_mets2.smile(LOCB(i),:));
            model.metPubChemID(i,1) =  cellstr(num2str(VMH_mets2.pubChemId(LOCB(i),:)));
            model.metMetaNetXID(i,1) =  cellstr(VMH_mets2.metanetx(LOCB(i),:));
            model.metBiGGID(i,1) =  cellstr(VMH_mets2.biggId(LOCB(i),:));
            model.metBioCycID(i,1) =  cellstr(VMH_mets2.biocyc(LOCB(i),:));
            model.metSEEDID(i,1) =  cellstr(VMH_mets2.seed(LOCB(i),:));
            model.metSBOTerms(i,1) =  cellstr('SBO:0000247'); % simple molecule
            if ~isempty(VMH_mets2.chemspider(LOCB(i),:)) %&& ~isnan(VMH_mets2.chemspider(LOCB(i),:))
                model.metChemSpider(i,1) =  cellstr(num2str(VMH_mets2.chemspider(LOCB(i),:)));
            end
            model.metInchiKey(i,1) =  cellstr(VMH_mets2.inchiKey(LOCB(i),:));
        end
    end
    % catch also those metabolites that have no organ_ in front of it
    [rem,rem2] = strtok(model.mets,'[');

    [LIA,LOCBa] = ismember(rem,VMH_mets.abbreviation);
    for i = 1 : length(LOCBa)
        if LOCBa(i) >0 && ~strcmp(a(i),'slack')
            model.metCharges(i,1) = str2num(VMH_mets2.charge(LOCBa(i),:));
            model.metFormulas(i,1) = cellstr(VMH_mets2.chargedFormula(LOCBa(i),:));
            model.metNames(i,1) =  cellstr(VMH_mets2.fullName(LOCBa(i),:));
            if ~isempty(regexprep(VMH_mets2.cheBlId(LOCBa(i),:),' ',''))% remove space in string
                model.metChEBIID(i,1) =  cellstr(strcat('CHEBI:',num2str(VMH_mets2.cheBlId(LOCBa(i),:))));
            end
            model.metHMDBID(i,1) =  cellstr(VMH_mets2.hmdb(LOCBa(i),:));
            model.metInChIString(i,1) =  cellstr(regexprep(VMH_mets2.inchiString(LOCBa(i),:),'"',''));
            model.metKEGGID(i,1) =  cellstr(VMH_mets2.keggId(LOCBa(i),:));
            model.metSmiles(i,1) =  cellstr(VMH_mets2.smile(LOCBa(i),:));
            model.metPubChemID(i,1) =  cellstr(num2str(VMH_mets2.pubChemId(LOCBa(i),:)));
            model.metMetaNetXID(i,1) =  cellstr(VMH_mets2.metanetx(LOCBa(i),:));
            model.metBiGGID(i,1) =  cellstr(VMH_mets2.biggId(LOCBa(i),:));
            model.metBioCycID(i,1) =  cellstr(VMH_mets2.biocyc(LOCBa(i),:));
            model.metSEEDID(i,1) =  cellstr(VMH_mets2.seed(LOCBa(i),:));
            model.metSBOTerms(i,1) =  cellstr('SBO:0000247'); % simple molecule
            if ~isempty(VMH_mets2.chemspider(LOCBa(i),:))% && ~isnan(VMH_mets2.chemspider(LOCBa(i),:))
                model.metChemSpider(i,1) =  cellstr(num2str(VMH_mets2.chemspider(LOCBa(i),:)));
            end
            model.metInchiKey(i,1) =  cellstr(VMH_mets2.inchiKey(LOCBa(i),:));
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function spqr_demo
%SPQR_DEMO short demo of SuiteSparseQR
%
% Example:
%   spqr_demo
%
% See also SPQR, SPQR_SOLVE, SPQR_QMULT, SPQR_MAKE, SPQR_INSTALL.

%   Copyright 2008, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

more on
help spqr
help spqr_solve
help spqr_qmult
help spqr_demo
more off

input ('Hit enter to start the SuiteSparseQR demo: ', 's') ;
fprintf ('\nTesting SuiteSparseQR functions ... please wait ...\n') ;

load west0479 ;
A = west0479 ;
rand ('state', 0) ;     %#ok
m = size (A,1) ;

figure (1)
clf

maxerr = 0 ;

for acomplex = 0:1

    if (acomplex)
        A = A + 1i * sprand (A) ;
    end
    anorm = norm (A,1) ;

    R1 = spqr (A) ;
    err = norm (R1'*R1 - A'*A, 1) / anorm^2 ;
    maxerr = max (maxerr, err) ;

    [Q,R,E] = spqr (A) ;
    err = norm (Q*R-A*E, 1) / anorm ;
    maxerr = max (maxerr, err) ;

    [H,R,E] = spqr (A, struct ('Q', 'Householder')) ;
    err = norm (spqr_qmult (H,R,1) - A*E, 1) / anorm ;
    maxerr = max (maxerr, err) ;

    if (acomplex)
        subplot (2,5,7)  ; spy (R1)  ; title ('R, no permutation (complex)') ;
        subplot (2,5,8)  ; spy (R)   ; title ('R with colamd (complex)') ;
        subplot (2,5,9)  ; spy (Q)   ; title ('Q with colamd (complex)') ;
        subplot (2,5,10) ; spy (H.H) ; title ('H with colamd (complex)') ;
    else
        subplot (2,5,1)  ; spy (A)   ; title ('A') ;
        subplot (2,5,2)  ; spy (R1)  ; title ('R, no permutation') ;
        subplot (2,5,3)  ; spy (R)   ; title ('R with colamd') ;
        subplot (2,5,4)  ; spy (Q)   ; title ('Q with colamd') ;
        subplot (2,5,5)  ; spy (H.H) ; title ('H with colamd') ;
    end
    drawnow

    % test spqr_solve: real/complex, sparse/full,
    % single/multiple right-hand-sides
    for bcomplex = 0:1
        for bsparse = 0:1
            for nrhs = 0:10
                if (bsparse)
                    b = sprand (m, nrhs, 0.1) ;
                else
                    b = rand (m, nrhs) ;
                end
                if (bcomplex)
                    b = b + 1i * sprand (b) ;
                end
                x = spqr_solve (A,b) ;
                err = norm (A*x-b,1) / max (anorm * norm(x,1) + norm (b,1), 1) ;
                maxerr = max (maxerr, err) ;
            end
        end
    end
end

fprintf ('\nQR maximum error: %g\n\n', maxerr) ;
if (maxerr > 1e-12)
    error ('One or more tests failed; error is high!') ;
end

% ------------------------------------------------------------------------------
% compare spqr_solve
% ------------------------------------------------------------------------------

fprintf ('Compare performance with MATLAB on a dense least-squares problem:\n');
A = rand (2000,1000) ;
b = rand (2000,1) ;
S = sparse (A) ;
fprintf ('\nA = rand (2000,1000) ;\nb = rand (2000,1) ;\nS = sparse (A) ;\n') ;

fprintf ('tic, x = spqr_solve(S,b) ; toc  ') ;
tic
x = spqr_solve (S, b) ;
t1 = toc ;
r1 = norm (A*x-b,1) ;
fprintf ('%% time %8.3f residual %8.3e\n', t1,r1) ;

fprintf ('tic, x = A\\b ; toc              ') ;
tic
x = A\b ;
t2 = toc ;
r2 = norm (A*x-b,1) ;
fprintf ('%% time %8.3f residual %8.3e\n', t2,r2) ;

fprintf ('tic, x = S\\b ; toc              ') ;
tic
x = S\b ;
t3 = toc ;
r3 = norm (A*x-b,1) ;
fprintf ('%% time %8.3f residual %8.3e\n', t3,r3) ;

% ------------------------------------------------------------------------------
% compare spqr with 1000-by-2000 system
% ------------------------------------------------------------------------------

fprintf ('\nA = rand (1000,2000) ;\nS = sparse (A) ;\n') ;
A = rand (1000,2000) ;
S = sparse (A) ;

fprintf ('tic, R = spqr(S) ; toc          ') ;
tic
R = spqr (S) ;
t1 = toc ;
r1 = norm (R'*R-A'*A,1) / norm(A,1)^2 ;
fprintf ('%% time %8.3f error %8.3e\n', t1,r1) ;
clear R

fprintf ('tic, R = qr(A)   ; toc          ') ;
tic
R = qr (A) ;
t2 = toc ;
R = triu (R) ;
r2 = norm (R'*R-A'*A,1) / norm(A,1)^2 ;
fprintf ('%% time %8.3f error %8.3e\n', t2,r2) ;
clear R

fprintf ('tic, R = qr(S)   ; toc          ') ;
tic
R = qr (S) ;
t3 = toc ;
r3 = norm (R'*R-A'*A,1) / norm(A,1)^2 ;
fprintf ('%% time %8.3f error %8.3e\n', t3,r3) ;
clear R

% ------------------------------------------------------------------------------
% compare spqr with 100-by-20000 system
% ------------------------------------------------------------------------------

A = rand (100,20000) ;
S = sparse (A) ;
fprintf ('\nA = rand (100,20000) ;\nS = sparse (A) ;\n') ;

fprintf ('tic, R = spqr(S) ; toc          ') ;
tic
R = spqr (S) ;                                                              %#ok
t1 = toc ;
fprintf ('%% time %8.3f\n', t1) ;
clear R

fprintf ('tic, R = qr(A)   ; toc          ') ;
tic
R = qr (A) ;                                                                %#ok
t2 = toc ;
fprintf ('%% time %8.3f\n', t2) ;
clear R

% skip the old MATLAB QR ... it's way to slow ...
%   fprintf ('tic, R = qr(S)   ; toc          ') ;
%   try
%       tic
%       R = qr (S) ;                                                        %#ok
%       t3 = toc ;
%       fprintf ('%% time %8.3f\n', t3) ;
%   catch                                                                   %#ok
%       fprintf ('%% MATLAB sparse qr failed ...\n') ;
%   end
%   clear R

fprintf ('All spqr tests passed\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
      function result = eq(A, B)
         %eq  equal for ChannelTags
         %expect that both parts are ChannelTags
         %
         % wildcards (simple '*' for fields) are automatically matched, so
         % be careful when comparing wildcards!
         %
         % empty fields are treated as intentionally empty.
         %
         % to ignore wildcards altogether, compare the fields individually!
         %
         % Wildcards must be in a scalar.
         %
         % See also matching

         %either A or B is scholar, or they are same size
         %eq is called once for each A (tested Nov16,2015)
         assert(all(size(A)==size(B)) || numel(A)==1 || numel(B)==1,...
            'Size of A and B must be equal, unless either is a scalar');

         stamatches = comparefield({A.station}, {B.station});
         chamatches = comparefield({A.channel}, {B.channel});
         locmatches = comparefield({A.location}, {B.location});
         netmatches = comparefield({A.network}, {B.network});
         result = stamatches & chamatches & locmatches & netmatches;
         if numel(A) == 1
            result = reshape(result,size(B));
         else
            result = reshape(result,size(A));
         end
      end

      function matchIdx = comparefield(A,B)
         % any field
         matchIdx = strcmp(A,B) | strcmp(A,'*') | strcmp(B,'*');
      end

author: geoscience-community-codes

repo: GISMO


text:
function legendre_function_q_test ( )

%*****************************************************************************80
%
%% LEGENDRE_FUNCTION_Q_TEST tests LEGENDRE_FUNCTION_Q.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_FUNCTION_Q_TEST:\n' );
  fprintf ( 1, '  LEGENDRE_FUNCTION_Q evaluates the Legendre Q function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N      X        Exact F       Q(N)(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx ] = legendre_function_q_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = legendre_function_q ( n, x );

    fprintf ( 1, '  %6d  %6f  %12f  %12f\n', n, x, fx, fx2(n+1) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function obj = windsorize(obj, varargin)
% Windsorize an fMRI data object to madlimit Median Absolute Deviations.
% Default = 5 MADs.
% Works across rows and columns.
% Registers this step in history.
%
% :Usage:
% ::
%
%     obj = windsorize(obj, [madlimit])

% ..
%    Calculate and display descriptives
% ..

madlimit = 5;
if ~isempty(varargin), madlimit = varargin{1}; end

dattmp = obj.dat(:);
med = median(dattmp);
mabsd = mad(dattmp);
climits = [med - madlimit * mabsd med + madlimit * mabsd];
wh_out = obj.dat < climits(1) | obj.dat > climits(2);
nout = sum(wh_out(:));
percout = 100 * nout ./ prod(size(obj.dat));

nout_by_case = sum(wh_out, 2);

sep = sprintf('____________________________________________\n');
fprintf('%sfmri_data structure .dat field (data)\n%s', sep, sep);
fprintf('Median: %3.3f\nMean: %3.3f\n', med, mean(dattmp));
fprintf('MAD: %3.3f\nSTD: %3.3f\n', mabsd, std(dattmp));
fprintf('Max:%3.3f\nMin:%3.3f\n', max(dattmp), min(dattmp));
fprintf('Control limits: %3.3f to %3.3f\n', climits);
fprintf('Outliers: %3.0f values, %3.2f%% of all values\n', nout, percout);
fprintf('Outliers by case (row): Max = %3.0f, Median = %3.0f, Min = %3.0f\n', max(nout_by_case), median(nout_by_case), min(nout_by_case));

datadj = obj.dat;
datadj(obj.dat < climits(1)) = climits(1);
datadj(obj.dat > climits(2)) = climits(2);

obj.dat = datadj;
obj.history{end+1} = sprintf('dat windsorized to %3.1f MADs', madlimit);

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function test_suite = test_bayesoptimization2
    %   Run specific demo, save values and compare the results to the expected.
    %   Works for both xUnit Test Framework package by Steve Eddins and for
    %   the built-in Unit Testing Framework (as of Matlab version 2013b).
    %
    %   See also
    %     TEST_ALL, DEMO_BAYEESOPTIMIZATION
    %
    %  Copyright (c) 2015 Jarno Vanhatalo
    %  Copyright (c) 2017 Eero Siivola

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

    % Check which package to use
    if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
    else
    initTestSuite_path = '';
    end
    if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
    else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
    end

end

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.

  % Suppress an expected warning
  Sw = warning('off','MATLAB:nearlySingularMatrix');
  rundemo(getName())
  warning(Sw); % Return the original warning visibility status
end

function testValues(testCase)
  % Test predictive mean and variance
  verifyVarsEqual(testCase, getName(), ...
    {'x'}, ...
    'RelTolElement', 0.01, 'RelTolRange', 0.02)
end

function testEI(testCase)
  % Test predictive mean and variance
  verifyVarsEqual(testCase, getName(), ...
    {'EI'}, ...
    'RelTolElement', 0.01, 'RelTolRange', 0.02)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

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


text:
function show3DModel(faces,vertices,showNormal)
%
% this function draws a 3D CG model with the normal of their faces
% the normal is computed by left-hand-rule
% http://stackoverflow.com/questions/1516296/find-the-normal-angle-of-the-face-of-a-triangle-in-3d-given-the-co-ordinated-of
% i.e. vertices are ordered clockwise with respect to its outward normal
%
% input:
% faces is Nx3 matrix for vertex index
% vertices is Kx3 matrix
%
% demo:
% load model;
% show3DModel(faces,vertices);

X = reshape(vertices(faces(:),1),size(faces))';
Y = reshape(vertices(faces(:),2),size(faces))';
Z = reshape(vertices(faces(:),3),size(faces))';

fill3(X,Y,Z,[160,160,160]/255);

cX = mean(X(1:3,:));
cY = mean(Y(1:3,:));
cZ = mean(Z(1:3,:));

p21= [X(2,:) - X(1,:); Y(2,:) - Y(1,:); Z(2,:) - Z(1,:)];
p31= [X(3,:) - X(1,:); Y(3,:) - Y(1,:); Z(3,:) - Z(1,:)];

nXYZ = cross(p21,p31);
lXYZ = sqrt(sum(nXYZ.^2,1));
nXYZ = nXYZ ./ repmat(lXYZ,3,1);

hold on;
if ~exist('showNormal','var') || showNormal
    quiver3(cX,cY,cZ,nXYZ(1,:),nXYZ(2,:),nXYZ(3,:));
end
axis equal

author: zhirongw


text:
function HRVparams = InitializeHRVparams(project_name)
%
%   settings = InitializeHRVparams('project_name')
%
%   OVERVIEW:
%       This file stores settings and should be configured before
%       each use of the PhysioNet Cardiovascular Signal Toolbox:
%       1.  Project Specific Input/Output Data type and Folders
%       2.  How much does the user trust the data
%       3.  Global Settings (Window Size) for signal segmentation
%       4.  Quality Threshold Settings
%       5.  Debug Settings
%       6.  SQI Settings
%       7.  Preprocess Settings
%       8. AF Detection Settings
%       9. Time Domain Analysis Settings
%       10. Frequency Domain Analysis Settings
%       11. SDANN and SDNNI Analysis Settings
%       12. PRSA Analysis Settings
%       13. Peak Detection Settings
%       14. Entropy  and Multiscale Entropy - MSE - Settings
%       15. Detrended Fluctuation Analysis - DFA - Settings
%       16. Poincar� plot - Settings
%       17. Heart Rate Turbulence (HRT) - Settings
%       18. Output Settings
%       19. Time of Process and Filename to Save Data
%
%   INPUT:
%       project_name = a string with the name of the project - this
%       will determine the naming convention of file folders
%
%   OUTPUT:
%       HRVparams - struct of various settings for the hrv_toolbox analysis
%
%   DEPENDENCIES & LIBRARIES:
%       PhysioNet Cardiovascular Signal Toolbox
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%
%   REFERENCE:
%   Vest et al. "An Open Source Benchmarked HRV Toolbox for Cardiovascular
%   Waveform and Interval Analysis" Physiological Measurement (In Press), 2018.
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Adriana N. Vest
%       Giulia Da Poian
%       Dependent scripts written by various authors
%       (see functions for details)
%	COPYRIGHT (C) 2018
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%%

if isempty(project_name)
    project_name = 'none';
end

% Set up input options for a specific project
switch project_name
    % Define new project name and parameters
    case 'MyProjectName'                   % Update with your project name
        HRVparams.Fs = NaN;                % Spacify sampling frequency
        HRVparams.readdata = '';           % (Optional) Specify name for data input folder
        HRVparams.writedata = '';          % (Optional) Specify name for data output folder
        HRVparams.datatype = '';           % (Optional) Spacify Data type of input
        HRVparams.ext = '';                % (Optional) Spacify file extension of input (e.g., 'mat','qrs')

    % Existing demo projects
    case 'MVanalysis'                      % Morphological variability analysis ECG
        HRVparams.Fs = 1000;
        HRVparams.readdata = '';
        HRVparams.writedata = '';
        HRVparams.datatype = '';
        HRVparams.ext = '';
    case 'demo_NSR'      % Parameters for demo using MIT nsr data
        HRVparams.readdata = strcat('TestData', filesep, 'Physionet_nsr2db');
        HRVparams.writedata = strcat('OutputData', filesep, 'ResultsNSR');
        HRVparams.ext = 'ecg';
        HRVparams.Fs = 125;
    case 'demo_RRgen'
        HRVparams.writedata = strcat('OutputData', filesep, 'ResultsRRgen');
        HRVparams.Fs = 125;
        HRVparams.demo.length = 5*60*60;   % Length of demo RR intervals in seconds
        HRVparams.demo.pe = 0.0003;        % Probability of ectopy ~ 1 per hr
        HRVparams.demo.pn = 0.0048;        % Probability of noise ~ 16 per hr
        HRVparams.demo.seed = 1;           % Seed for RRGEN
    case 'demoICU'      % Parameters for ICU demo
        HRVparams.readdata = strcat('TestData');
        HRVparams.writedata = strcat('OutputData', filesep, 'ResultsICU');
        HRVparams.ext = 'mat';
        HRVparams.Fs = 128;
    case 'demoAF'      % Parameters for AF demo
        HRVparams.readdata = strcat('TestData');
        HRVparams.writedata = strcat('OutputData', filesep, 'ResultsAFData');
        HRVparams.ext = 'mat';
        HRVparams.Fs = 128;
    case 'demoHRT'      % Parameters for HRT demo
        HRVparams.readdata = strcat('TestData', filesep, 'Physionet_nsr2db');
        HRVparams.writedata = strcat('OutputData', filesep, 'ResultsHRT');
        HRVparams.Fs = 128;
    case 'demoPVC'      % Parameters for PVC demo
        HRVparams.readdata = strcat('TestData', filesep, 'mitdb-Arrhythmia');
        HRVparams.writedata = strcat('OutputData', filesep, 'ResultsPVC');
        HRVparams.Fs = 360;
    otherwise                  % Default
        HRVparams.Fs = NaN;                          % Spacify sampling frequency
        HRVparams.writedata = 'HRV_Output';          % (Optional) Specify name for data output folder

end

if  isempty(HRVparams.writedata)
    % Default data OUTPUT folder name based on project name
    HRVparams.writedata = strcat(project_name,'_Results');
    fprintf('New OUTPUT folder: "%s"\n', HRVparams.writedata)
    mkdir(HRVparams.writedata);          % Create output folder and
elseif ~exist([pwd filesep HRVparams.writedata], 'dir')
    fprintf('New OUTPUT folder: "%s"\n',HRVparams.writedata)
    mkdir(HRVparams.writedata);          % Create output folder and
end
addpath(genpath(HRVparams.writedata));   % Add folder to search path

%% 2. How much does the user trust the data:
% This setting determines how stringently filtered the data is prior to
% calculating HRV metrics. Raw ECG or Pulse data that has been labeled
% using a peak detector (like jqrs) would require the most stringent
% filtering, whereas RR interval data that has been reviewed by a human
% technician would require the least amount of filtering.

%   - qrs detection no beats labeled - most stringent filters
%   - automatic beat detection beat typed - moderately stringent filtered
%   - hand corrected - least filtered, maybe no filter

% EXPLAIN THE THRESHOLDS BETTER
% ADD THIS TO A DEMO

HRVparams.data_confidence_level = 1;
% 1 - raw data with automatic beat detection
% 2 - raw data with automatic beat detection, but beat typed (ie N, SV,etc)
% 3 - technician reviewed data

% * ^^^^ NOT YET IN USE ^^^^  *

%% 3. Global Settings (Window Size)

HRVparams.windowlength = 300;	      % Default: 300, seconds
HRVparams.increment = 30;             % Default: 30, seconds increment
HRVparams.numsegs = 5;                % Default: 5, number of segments to collect with lowest HR
HRVparams.RejectionThreshold = .20;   % Default: 0.2, amount (%) of data that can be rejected before a
                                      % window is considered too low quality for analysis
HRVparams.MissingDataThreshold = .15; % Default: 0.15, maximum percentage of data allowable to be missing
                                      % from a window .15 = 15%
%% 5. Debug Settings

HRVparams.rawsig = 0;           % Load raw signal if it is available for debugging
HRVparams.debug = 0;

%% 6. SQI Analysis Settings
HRVparams.sqi.LowQualityThreshold = 0.9; % Default: 0.9, Threshold for which SQI represents good data
HRVparams.sqi.windowlength = 10;         % Default: 10, seconds, length of the comparison window
HRVparams.sqi.increment = 1;             % Default: 1, seconds
HRVparams.sqi.TimeThreshold = 0.1;       % Default: 0.1, seconds
HRVparams.sqi.margin = 2;                % Default: 2, seconds, Margin time not include in comparison

%% 7. Preprocess Settings

HRVparams.preprocess.figures = 0;                   % Figures on = 1, Figures off = 0
HRVparams.preprocess.gaplimit = 2;                  % Default: 2, seconds; maximum believable gap in rr intervals
HRVparams.preprocess.per_limit = 0.2;               % Default: 0.2, Percent limit of change from one interval to the next
HRVparams.preprocess.forward_gap = 3;	            % Default: 3, Maximum tolerable gap at beginning of timeseries in seconds
HRVparams.preprocess.method_outliers = 'rem';       % Default: 'rem', Method of dealing with outliers
                                                    % 'cub' = replace outlier points with cubic spline method
                                                    % 'rem' = remove outlier points
                                                    % 'pchip' = replace with pchip method
HRVparams.preprocess.lowerphysiolim = 60/160;       % Default: 60/160
HRVparams.preprocess.upperphysiolim = 60/30;        % Default: 60/30
HRVparams.preprocess.method_unphysio = 'rem';       % Default: 'rem', Method of dealing with unphysiologically low beats
                                                    % 'cub' = replace outlier points with cubic spline method
                                                    % 'rem' = remove outlier points
                                                    % 'pchip' = replace with pchip method

% The following settings do not yet have any functional effect on
% the output of preprocess.m:
HRVparams.preprocess.threshold1 = 0.9 ;	        % Default: 0.9, Threshold for which SQI represents good data
HRVparams.preprocess.minlength = 30;            % Default: 30, The minimum length of a good data segment in seconds

%% 8. AF Detection Settings and PVC detection

HRVparams.af.on = 1;              % Default: 1, AF Detection On or Off
HRVparams.af.windowlength = 30;   % Default: 30, Set to include ~30 beats in each window
HRVparams.af.increment = 30;      % Default: 30, No overlap necessary in AF feat calc

HRVparams.PVC.qrsth = 0.1;        % Default: 0.1, threshold for qrs detection

%% 9. Time Domain Analysis Settings

HRVparams.timedomain.on = 1;             % Default: 1, Time Domain Analysis 1=On or 0=Off
HRVparams.timedomain.dataoutput = 0;     % 1 = Print results to .txt file
                                         % Anything else = utputs to return variables only
                                         % returned variables
HRVparams.timedomain.alpha = 50   ;      % Default: 50 ,In msec
HRVparams.timedomain.win_tol = .15;      % Default: .15, Maximum percentage of data allowable
                                         % to be missing from a window

%% 10. Frequency Domain Analysis Settings

HRVparams.freq.on = 1;              % Default: 1, Frequency Domain Analysis 1=On or 0=Off

ULF = [0 .0033];                    % Requires window > 300 s
VLF = [0.0033 .04];                 % Requires at least 300 s window
LF = [.04 .15];                     % Requires at least 25 s window
HF = [0.15 0.4];                    % Requires at least 7 s window

HRVparams.freq.limits = [ULF; VLF; LF; HF];
HRVparams.freq.zero_mean = 1;               % Default: 1, Option for subtracting the mean from the input data
HRVparams.freq.method = 'lomb';             % Default: 'lomb'
                                            % Options: 'lomb', 'burg', 'fft', 'welch'
HRVparams.freq.plot_on = 0;

% The following settings are for debugging spectral analysis methods
HRVparams.freq.debug_sine = 0;              % Default: 0, Adds sine wave to tachogram for debugging
HRVparams.freq.debug_freq = 0.15;           % Default: 0.15
HRVparams.freq.debug_weight = .03;          % Default: 0.03

% Lomb:
HRVparams.freq.normalize_lomb = 0;	        % Default: 0
                                            % 1 = Normalizes Lomb Periodogram,
                                            % 0 = Doesn't normalize

% Burg: (not recommended)
HRVparams.freq.burg_poles = 15;    % Default: 15, Number of coefficients
                                   % for spectral estimation using the Burg
                                   % method (not recommended)

% The following settings are only used when the user specifies spectral
% estimation methods that use resampling : 'welch','fft', 'burg'
HRVparams.freq.resampling_freq = 7;             % Default: 7, Hz
HRVparams.freq.resample_interp_method = 'cub';  % Default: 'cub'
                                                % 'cub' = cublic spline method
                                                % 'lin' = linear spline method
HRVparams.freq.resampled_burg_poles = 100;      % Default: 100

%% 11. SDANN and SDNNI Analysis Settings
HRVparams.sd.on = 1;                        % Default: 1, SD analysis 1=On or 0=Off
HRVparams.sd.segmentlength = 300;           % Default: 300, windows length in seconds

%% 12. PRSA Analysis Settings

HRVparams.prsa.on = 1;             % Default: 1, PRSA Analysis 1=On or 0=Off
HRVparams.prsa.win_length = 30;    % Default: 30, The length of the PRSA signal
                                   % before and after the anchor points
                                   % (the resulting PRSA has length 2*L)
HRVparams.prsa.thresh_per = 20;    % Default: 20%, Percent difference that one beat can
                                   % differ from the next in the prsa code
HRVparams.prsa.plot_results = 0;   % Default: 0
HRVparams.prsa.scale = 2;          % Default: 2, scale parameter for wavelet analysis (to compute AC and DC)
HRVparams.prsa.min_anch = 20;       % Default: 20, minimum number of anchors point required to create the "average signal"

%% 13. Peak Detection Settings

% The following settings are for jqrs.m

HRVparams.PeakDetect.REF_PERIOD = 0.250;   % Default: 0.25 (should be 0.15 for FECG), refractory period in sec between two R-peaks
HRVparams.PeakDetect.THRES = .6;           % Default: 0.6, Energy threshold of the detector
HRVparams.PeakDetect.fid_vec = [];         % Default: [], If some subsegments should not be used for finding the optimal
                                           % threshold of the P&T then input the indices of the corresponding points here
HRVparams.PeakDetect.SIGN_FORCE = [];      % Default: [], Force sign of peaks (positive value/negative value)
HRVparams.PeakDetect.debug = 0;            % Default: 0
HRVparams.PeakDetect.ecgType = 'MECG';     % Default : MECG, options (adult MECG) or featl ECG (fECG)
HRVparams.PeakDetect.windows = 15;         % Befautl: 15,(in seconds) size of the window onto which to perform QRS detection

%% 14. Entropy Settings
% Multiscale Entropy
HRVparams.MSE.on = 1;                      % Default: 1, MSE Analysis 1=On or 0=Off
HRVparams.MSE.windowlength = [];           % Default: [], windows size in hours, default perform MSE on the entire signal
HRVparams.MSE.increment = [];              % Default: [], window increment in hours
HRVparams.MSE.RadiusOfSimilarity = 0.15;   % Default: 0.15, Radius of similarity (% of std)
HRVparams.MSE.patternLength = 2;           % Default: 2, pattern length
HRVparams.MSE.maxCoarseGrainings = 20;     % Default: 20, Maximum number of coarse-grainings
HRVparams.MSE.method = 'fir';           % Default 'fir', method use to generate the coarse-grain
                                           %                time series  options 'fir', 'butter'
HRVparams.MSE.moment = 'mean';             % Default: 'mean', moment used to coarse-grain the time series
                                           %           options 'mean', 'variance'
HRVparams.MSE.constant_r = 1;              % Default: 1, 1 use r as function of std of original time
                                           %             series, 0 compute r as function of std at each scale

% SampEn an ApEn
HRVparams.Entropy.on = 1;                     % Default: 1, MSE Analysis 1=On or 0=Off
HRVparams.Entropy.RadiusOfSimilarity = 0.15;  % Default: 0.15, Radius of similarity (% of std)
HRVparams.Entropy.patternLength = 2;          % Default: 2, pattern length

%% 15. DFA Settings

HRVparams.DFA.on = 1;               % Default: 1, DFA Analysis 1=On or 0=Off
HRVparams.DFA.windowlength = [];    % Default [], windows size in hours, default perform DFA on the entair signal
HRVparams.DFA.increment = [];       % Default: [], window increment in hours
HRVparams.DFA.minBoxSize = 4 ;      % Default: 4, Smallest box width
HRVparams.DFA.maxBoxSize = [];      % Largest box width (default in DFA code: signal length/4)
HRVparams.DFA.midBoxSize = 16;      % Medium time scale box width (default in DFA code: 16)

%% 16. Poincar� plot

HRVparams.poincare.on = 1;     % Default: 1, Poincare Analysis 1=On or 0=Off

%% 17. Heart Rate Turbulence (HRT) - Settings

HRVparams.HRT.on = 1;                        % Default: 1, HRT Analysis 1=On or 0=Off
HRVparams.HRT.BeatsBefore = 2;               % Default: 2, # of beats before PVC
HRVparams.HRT.BeatsAfter = 16;               % Default: 16, # of beats after PVC and CP
HRVparams.HRT.GraphOn = 0;                   % Default: 0, do not plot
HRVparams.HRT.windowlength = 24;             % Default 24h, windows size in hours
HRVparams.HRT.increment = 24;                % Default 24h, sliding window increment in hours
HRVparams.HRT.filterMethod = 'mean5before';  % Default mean5before, HRT filtering option

%% 18. Output Settings

HRVparams.gen_figs = 0;             % Generate figures
HRVparams.save_figs = 0;            % Save generated figures
if HRVparams.save_figs == 1
    HRVparams.gen_figs = 1;
end

% Format settings for HRV Outputs
HRVparams.output.format = 'csv';        % 'csv' - creates csv file for output
                                        % 'mat' - creates .mat file for output
HRVparams.output.separate = 1;          % Default : 1 = separate files for each subject
                                        %           0 = all results in one file
HRVparams.output.num_win = [];          % Specify number of lowest hr windows returned
                                        % leave blank if all windows should be returned
                                        % Format settings for annotations generated
HRVparams.output.ann_format = 'binary'; % 'binary'  = binary annotation file generated
                                        % 'csv'     = ASCII CSV file generated

%% 19. Filename to Save Data
HRVparams.time = datestr(now, 'yyyymmdd');              % Setup time for filename of output
HRVparams.filename = [HRVparams.time '_' project_name];

%% Export Parameter as Latex Table
% Note that if you change the order of the parameters or add parameters
% this might not work

ExportHRVparams(HRVparams);

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
%tstoolbox/mex/baker
%   Generate time-series from the iterated Baker map .
%
%   Syntax:
%
%     * x = baker(length, [eta l1 l2 x0 y0])
%
%   Input arguments:
%
%     * length - number of samples to generate
%     * [eta l1 l2 x0 y0] - vector of parameters and initial conditions
%
%   Output arguments:
%
%     * x - time series
%
%   Example:
%
%x = baker(2000, [0.6 0.25 0.4 rand(1,1) rand(1,1)]);
%plot(x(1:end-1,2), x(2:end,2), '.')
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function data=createdatamatpb(data,E,Fs,win)
%
% Helper function to create an event triggered matrix from a single
% channel of data.
% Usage: data=createdatamatpb(data,E,Fs,win)
% Inputs:
% data   (input time series as a single vector) - required
% E      (events to use as triggers) - required
% Fs     (sampling frequency of data) - required
% win    (window around triggers to use data matrix -[winl winr]) - required
%          e.g [1 1] uses a window starting 1 sec before E and
%              ending 1 sec after E if E is in secs
% Note that E, Fs, and win must have consistent units
%
% Outputs:
% data      (event triggered data)
%
if nargin < 4
   error('Need all input arguments');
end
NE=length(E);
nwinl=round(win(1)*Fs);
nwinr=round(win(2)*Fs);
nE=floor(E*Fs)+1;
datatmp=[];
for n=1:NE;
    indx=nE(n)-nwinl:nE(n)+nwinr-1;
    datatmp=[datatmp data(indx)];
end
data=datatmp;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Load far talk, close talk, and frame label.
function [Data, para, vocab] = LoadParallelWavLabel_Reverb(para, step, dataset, datatype, distance)
nCh = para.topology.useChannel;
wavlist = []; wavlistClean = [];
% note that training data do not contain real recordings.
% real dev and eval data do not have clean version
for type_i = 1:length(datatype)
    if strcmpi(datatype, 'real') || strcmpi(dataset, 'train')
        nRoom = 1;    else; nRoom = 3; end
    for dist_i = 1:length(distance)
        for roomID = 1:nRoom
            [tasklist, taskfile] = LoadTaskFile_Reverb([para.local.wavroot '/taskFiles'], dataset, datatype{type_i}, distance{dist_i}, nCh, roomID);
            wavlist = [wavlist Tasklist2wavlist(para.local.wavroot, tasklist, dataset, datatype{type_i})];

            if strcmpi(dataset, 'train')
                wavlistClean = findFiles([para.local.wsjcam0root '/data_wav/primary_microphone/si_tr'], para.local.wsjcam0ext);
            elseif ~strcmpi(datatype, 'real')
                [tasklistClean, taskfileClean] = LoadTaskFile_Reverb([para.local.wavroot '/taskFiles'], dataset, datatype{type_i}, 'cln', nCh, roomID);
                wavlistClean = [wavlistClean Tasklist2wavlist(para.local.wavroot, tasklistClean, dataset, datatype{type_i})];
            end
        end
        if strcmpi(dataset, 'train'); break; end    % For training data, far and near distances are in one list, so we can stop here.
    end
end
for i=1:length(wavlistClean)   % build an index of clean files so we can find them quickly by utterance ID
    [~,curr_uttID] = fileparts(wavlistClean{i});
    clean_struct.(['U_' curr_uttID]) = wavlistClean{i};
end

% Currently we don't support loading of frame label yet.
% switch lower(dataset)
%     case {'eval_simu'}
%         ali_file = [para.local.aliDir '_' dataset '/ali.txt'];
%     case 'tr05'
%         ali_file = [para.local.aliDir '/ali.txt'];
% end
% [ali, vocab] = LoadKaldiFrameLabel(ali_file);
vocab = [];

wavlist = wavlist(:,step:step:end);
wavlistClean = wavlistClean(:,step:step:end);

wavreader.name = 'wavfile';
wavreader.array = 1;
wavreader.multiArrayFiles = 1;

fs = 16000;
frame_size = fs*0.025;
frame_shift = fs*0.01;
if para.local.useFileName==0
    seglen = para.local.seglen;
    segshift = para.local.segshift;
end

% Load data file list
wav_noisy = {};  label = {};  wav_clean = {};
nUtt = size(wavlist,2);
for si = 1:nUtt
    words = ExtractWordsFromString_v2(wavlist{1,si}, '/');
    curr_uttID = words{end}(1:end-4);
    PrintProgress(si, nUtt, 100, curr_uttID);

    words2 = ExtractWordsFromString_v2(curr_uttID, '_');
    clean_uttID = words2{1};

    [wav] = InputReader(wavlist(:,si), wavreader);
    wav = StoreWavInt16(wav);
    if strcmpi(datatype, 'real')
        wav_c = wav; % note that for real data, we don't have clean version. so just use the noisy version. This line will only be run at test time.
    elseif strcmpi(dataset, 'train')
        wavfileClean = clean_struct.(['U_' clean_uttID]);
        wav_c = audioread(wavfileClean)';
    else
        [wav_c] = InputReader(wavlistClean(:,si), wavreader);
    end
    wav_c = StoreWavInt16(wav_c);

    % synchronize the length of label and wav
    nFr_feat = enframe_decide_frame_number(size(wav,2), frame_size, frame_shift);
    nFr_feat = min(nFr_feat, enframe_decide_frame_number(size(wav_c,2), frame_size, frame_shift));
    nFr_label = 10000;
    nFr = min([nFr_feat nFr_label]);
    requiredLen = DecideWavLen4XFrames(nFr, frame_size, frame_shift);
    wav(:,requiredLen+1:end) = [];
    wav_c(:,requiredLen+1:end) = [];

    if para.local.useFileName
        for i=1:nCh
            wavfileArray{i} = sprintf('%s 0 %2.3f', wavlist{i,si}, size(wav,2)/fs);
        end
        wav_noisy{end+1} = wavfileArray;
        if strcmpi(datatype, 'real')
            wav_clean{end+1} = wavfileArray{1};
        elseif strcmpi(dataset, 'train')
            wav_clean{end+1} = sprintf('%s 0 %2.3f', clean_struct.(['U_' clean_uttID]), size(wav,2)/fs);
        else
            wav_clean{end+1} = sprintf('%s 0 %2.3f', wavlistClean{1,si}, size(wav,2)/fs);
        end
    else
        if 1    % we use small segments rather than whole sentences as training examples
            wav_c_seg = DivideSent2Segments(wav_c, (seglen-1)*frame_shift+frame_size, segshift*frame_shift, 1);
            wav_seg = DivideSent2Segments(wav, (seglen-1)*frame_shift+frame_size, segshift*frame_shift, 1);
            wav_noisy = [wav_noisy; wav_seg];
            wav_clean = [wav_clean; wav_c_seg];
        else
            wav_noisy{end+1} = wav;
            wav_clean{end+1} = wav_c;
        end
    end
    if para.local.loadLabel
        for j = length(label)+1:length(wav_noisy)
            label{j} = [];
        end
    end
end

Data(1).data = wav_noisy;
Data(2).data = wav_clean;

para.IO.inputFeature = [1 1];
para.IO.DataSyncSet{1} = [];
para.IO.frame_rate = [16000 16000];
para.IO.isTensor = [1 1];
if para.local.useFileName
    para.IO.inputFeature([1 2]) = 0;
    wavreader.precision = 'int16';
    if isfield(para.IO, 'fileReader')
        para.IO = rmfield(para.IO, 'fileReader');
    end
    para.IO.fileReader(1) = wavreader;
    para.IO.fileReader(2) = wavreader;
    para.IO.fileReader(2).array = 0;
    para.IO.fileReader(2).multiArrayFiles = 0;
end
if para.local.loadLabel
    Data(3).data = label;
    para.IO.inputFeature(3) = 1;
    para.IO.frame_rate(3) = 100;
    para.IO.isTensor(3) = 1;
    para.IO.fileReader(3).name = '';
end
end

author: singaxiong

repo: SignalGraph


text:
function C=corrdelays(A, tdelay)
% C=corrdelays(A, tdelays)
% Computes correlation matrices of MxN matrix A (rows are time series) for
% time delayes specified in the 1xT vector TDELAY

[m,n]=size(A);
t=length(tdelay);
C=zeros(m,m,t);
for indext = 1:t
    Ashift = (circshift(A',tdelay(indext)))';
    C(:,:,indext)=corr(A',Ashift');
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function pass = test_feval( )
% Test spherefun feval.

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

% See the random number generator.
rng(9);

f = @(x,y,z) sin(x+ y.*z) + 1;
f = redefine_function_handle(f);
g = spherefun(f);
lambda = rand;
theta = rand;
pass(1) = abs( g(theta, lambda) - f(theta, lambda) ) < tol;

% feval at vectors:
lambda = rand(10, 1);
theta = rand(10, 1);
pass(2) = norm( g(theta, lambda) - f(theta, lambda), inf ) < tol;

% feval at row vectors:
lambda = rand(1, 10);
theta = rand(1, 10);
pass(3) = norm(g(theta, lambda) - f(theta, lambda) ) < tol;

% feval at vectors:
lambda = rand(2, 10);
theta = rand(2, 10);
pass(4) = norm(g(theta, lambda) - f(theta, lambda), inf ) < tol;

% feval at meshgrid:
[lambda, theta] = meshgrid(rand(3, 1));
pass(5) = norm(g(theta, lambda) - f(theta, lambda), inf) < tol;

% feval using Cartesian coordinates
f = @(x,y,z) sin(x+ y.*z) + 1;
x = rand(1, 3);
x = x/sqrt(sum(x.^2));
y = x(1, 2);
z = x(1, 3);
x = x(1, 1);
pass(6) = norm(g(x, y, z) - f(x, y, z), inf) < tol;

% feval using Cartesian coordinates at a vector of points
x = rand(10, 3);
nrm = sqrt(sum(x.^2,2));
y = x(:, 2)./nrm;
z = x(:, 3)./nrm;
x = x(:, 1)./nrm;
pass(7) = norm(g(x, y, z) - f(x, y, z), inf) < tol;

% feval using Cartesian coordinates at a vector of points
x = rand(10, 3);
y = rand(10, 3);
z = rand(10, 3);
nrm = sqrt(x.^2 + y.^2 + z.^2);
y = y./nrm;
z = z./nrm;
x = x./nrm;
pass(8) = norm(g(x, y, z) - f(x, y, z), inf) < tol;

% Check for the appropriate error flag.
try
    y = g(1,1,1);
    pass(9) = false;
catch ME
    pass(9) = strcmp(ME.identifier,'CHEBFUN:SPHEREFUN:FEVAL:pointsNotOnSphere');
end

end

function f = redefine_function_handle(f)
% nargin(f) = 2, then we are already on the sphere, if nargin(f) = 3,
% then do change of variables:

if ( nargin(f) == 3 )
    % Wrap f so it can be evaluated in spherical coordinates
    f = @(lam, th) spherefun.sphf2cartf(f, lam, th, 0);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function scales = svargplvmShowScales(model, printPlots, modalities, dims, thresh, normaliseScales, customLeg, logScales)

% SVARGPLVMSHOWSCALES Show the scales of a svargplvmModel graphically
%
% SEEALSO : svargplvmFindSharedDims
%
% COPYRIGHT: Andreas C. Damianou, 2012, 2013
%
% VARGPLVM

if nargin < 2 || isempty(printPlots)
    printPlots = 2;
end

if nargin < 3 || isempty(modalities)
    modalities = 1:length(model.comp);
else
    printPlots = true;
end

if nargin < 4 || isempty(dims)
    dims = 1:model.q;
end

if nargin < 5 || isempty(thresh)
    thresh = -Inf;
end

if nargin < 6 || isempty(normaliseScales)
    normaliseScales = false;
end

if nargin < 7
    customLeg = {}; % Custom legend
end

if nargin < 8 || isempty(logScales)
    logScales = false;
end

totalModels = length(modalities);

if printPlots
    for i=1:length(model.comp)
        if strcmp(model.comp{i}.kern.type, 'rbfardjit')
            scales{i} = model.comp{i}.kern.inputScales(dims);
        else
            scales{i} = model.comp{i}.kern.comp{1}.inputScales(dims);
        end
    end
    if logScales
        origScales = scales;
        for k=1:totalModels
            scales{modalities(k)} = scales{modalities(k)}./max(scales{modalities(k)});
            scales{modalities(k)} = log(scales{modalities(k)});%logN(scales{modalities(k)},1.01);
            a=0; b=1; % Scale data between a and b
            X=scales{modalities(k)};
            X = (X-min(min(X))).*(b-a)./(max(max(X))-min(min(X))) + a;
            scales{modalities(k)}=X;
        end
    end
    if (totalModels == 2 || normaliseScales) && ~logScales
        origScales = scales;
        for k = 1:totalModels
            scales{modalities(k)}=scales{modalities(k)}./max(scales{modalities(k)});
        end
    end
    for k=1:totalModels
       scales{modalities(k)}(scales{modalities(k)}<thresh) = 0;
    end
    if totalModels == 2
        maxScales1 = max(scales{modalities(1)});
        maxScales2 = max(scales{modalities(2)});
        %scales{modalities(1)} = scales{modalities(1)}./maxScales1;
        %scales{modalities(2)} = scales{modalities(2)}./maxScales2;
        x=1:size(scales{modalities(1)},2);
        K=0.5;
        bar1=bar(x, scales{modalities(1)}, 'FaceColor', 'b', 'EdgeColor', 'b');
        set(bar1,'BarWidth',K);
        hold on;
        bar2=bar(x, scales{modalities(2)}, 'FaceColor', 'r', 'EdgeColor', 'r');
        set(bar2,'BarWidth',K/2);
        hold off;
        if isempty(customLeg)
            legend(['scales1 (x ' num2str(maxScales1) ')'],['scales2 (x ' num2str(maxScales2) ')'])
        else
            legend(customLeg);
        end
        scales = origScales;
    else
        % Too many models. Plot one at a time
        %if totalModels > 10 && exist('hsvargplvmShowScales','file')
             %mm.layer{1} = model;
             %mm.H=1;
             %mm.layer{1}.M = totalModels;
             %hsvargplvmShowScales(mm);
        %else
            for i=1:totalModels
                if totalModels > 10
                    pause
                    %figure;
                elseif totalModels > 1
                    subplot(1,totalModels, i) ;
                end
                bar(scales{modalities(i)}); title(['Modality ' num2str(modalities(i))])
            end
        %end
    end
else
    for i=1:length(model.comp)
        if strcmp(model.comp{i}.kern.type, 'rbfardjit')
            fprintf('# Scales of model %d: ', i)
            fprintf('%.4f  ',model.comp{i}.kern.inputScales(dims));
            fprintf('\n');
            scales{i} = model.comp{i}.kern.inputScales(dims);
        else
            fprintf('# Scales of model %d: ', i)
            fprintf('%.4f  ',model.comp{i}.kern.comp{1}.inputScales(dims));
            fprintf('\n');
            scales{i} = model.comp{i}.kern.comp{1}.inputScales(dims);
        end
        scales{i}(scales{i}<thresh) = 0;
    end
end

function y=logN(x,N)
y = log(x)/log(N);

author: SheffieldML

repo: vargplvm


text:
%Protein processing II process test case
%
% Author: Jared Jacobs, jmjacobs@stanford.edu
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 8/9/2010
classdef ProteinProcessingII_Test < edu.stanford.covert.cell.sim.ProcessTestCase
    methods
        function this = ProteinProcessingII_Test(methodName)
            this = this@edu.stanford.covert.cell.sim.ProcessTestCase(methodName);
        end

        function testOneMonomerRequiringNoProcessing(this)
            m = this.process;
            m.lipoproteinMonomerIndexs = [];
            m.secretedMonomerIndexs = [];
            m.unprocessedMonomerIndexs = 1;
            m.substrates(:) = 0;
            m.enzymes(:) = 0;
            m.unprocessedMonomers = 1;
            m.processedMonomers = 0;
            m.signalSequenceMonomers = 0;

            m.evolveState();
            assertEqual(0, m.unprocessedMonomers);
            assertEqual(1, m.processedMonomers);
            assertEqual(0, m.signalSequenceMonomers);
            assertEqual(zeros(size(m.substrates)), m.substrates);
        end

        function testOneSecretedMonomer(this)
            m = this.process;
            m.lipoproteinMonomerIndexs = [];
            m.secretedMonomerIndexs = 1;
            m.unprocessedMonomerIndexs = [];
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_water) = 1;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_signalPeptidase) = 1;
            m.unprocessedMonomers = 1;
            m.processedMonomers = 0;
            m.signalSequenceMonomers = 0;

            m.evolveState();
            assertEqual(0, m.unprocessedMonomers);
            assertEqual(1, m.processedMonomers);
            assertEqual(1, m.signalSequenceMonomers);
            assertEqual(zeros(size(m.substrates)), m.substrates);
            assertEqual(1, m.enzymes(m.enzymeIndexs_signalPeptidase));
        end

        function testOneLipoprotein(this)
            m = this.process;
            m.lipoproteinMonomerIndexs = 1;
            m.secretedMonomerIndexs = [];
            m.unprocessedMonomerIndexs = [];
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_water) = 1;
            m.substrates(m.substrateIndexs_PG160) = 1;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_signalPeptidase) = 1;
            m.enzymes(m.enzymeIndexs_diacylglycerylTransferase) = 99;
            m.unprocessedMonomers = 1;
            m.processedMonomers = 0;
            m.signalSequenceMonomers = 0;

            m.evolveState();
            assertEqual(0, m.unprocessedMonomers);
            assertEqual(1, m.processedMonomers);
            assertEqual(1, m.signalSequenceMonomers);
            assertEqual(0, m.substrates(m.substrateIndexs_water));
            assertEqual(0, m.substrates(m.substrateIndexs_PG160));
            assertEqual(1, m.substrates(m.substrateIndexs_hydrogen));
            assertEqual(1, m.substrates(m.substrateIndexs_SNGLYP));
            assertEqual(1, m.enzymes(m.enzymeIndexs_signalPeptidase));
            assertEqual(99, m.enzymes(m.enzymeIndexs_diacylglycerylTransferase));
        end

        function testNoProcessingWithoutWater(this)
            m = this.process;
            m.lipoproteinMonomerIndexs = 1;
            m.secretedMonomerIndexs = 2;
            m.unprocessedMonomerIndexs = [];
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_PG160) = 1e3;
            m.enzymes(:) = 1e3;
            m.unprocessedMonomers = [1; 1];
            m.processedMonomers = [0; 0];
            m.signalSequenceMonomers = [0; 0];

            m.evolveState();
            assertEqual([1 1], m.unprocessedMonomers');
            assertEqual([0 0], m.processedMonomers');
            assertEqual([0 0], m.signalSequenceMonomers');
            assertEqual(0, m.substrates(m.substrateIndexs_water));
            assertEqual(1e3, m.substrates(m.substrateIndexs_PG160));
        end

        function testNoProcessingWithoutSignalPeptidase(this)
            m = this.process;
            m.lipoproteinMonomerIndexs = 1;
            m.secretedMonomerIndexs = 2;
            m.unprocessedMonomerIndexs = [];
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.substrates(m.substrateIndexs_PG160) = 1e3;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_diacylglycerylTransferase) = 1e3;
            m.unprocessedMonomers = [1;1];
            m.processedMonomers = [0;0];
            m.signalSequenceMonomers = [0;0];

            m.evolveState();
            assertEqual([1 1], m.unprocessedMonomers');
            assertEqual([0 0], m.processedMonomers');
            assertEqual([0 0], m.signalSequenceMonomers');
            assertEqual(1e6, m.substrates(m.substrateIndexs_water));
            assertEqual(1e3, m.substrates(m.substrateIndexs_PG160));
        end

        function testNoLipoproteinProcessingWithoutDiacylglycerylTransferase(this)
            m = this.process;
            m.lipoproteinMonomerIndexs = 1;
            m.secretedMonomerIndexs = [];
            m.unprocessedMonomerIndexs = [];
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.substrates(m.substrateIndexs_PG160) = 1e3;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_signalPeptidase) = 1e3;
            m.unprocessedMonomers = 1;
            m.processedMonomers = 0;
            m.signalSequenceMonomers = 0;

            m.evolveState();
            assertEqual(1, m.unprocessedMonomers);
            assertEqual(0, m.processedMonomers);
            assertEqual(0, m.signalSequenceMonomers);
            assertEqual(1e6, m.substrates(m.substrateIndexs_water));
            assertEqual(1e3, m.substrates(m.substrateIndexs_PG160));
        end

        % Verifies that signal peptidase processes roughly as many monomers as
        % its rate will allow when it's the limiting factor, and that the
        % monomers chosen for processing are chosen without egregious bias.
        function testLimitedSignalPeptidase_secretedProteinsOnly(this)
            m = this.process;
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_signalPeptidase) = 3;
            m.unprocessedMonomers(:) = 0;
            m.unprocessedMonomers(m.secretedMonomerIndexs) = 10;
            m.processedMonomers(:) = 0;
            m.signalSequenceMonomers(:) = 0;

            m.evolveState();
            n = m.enzymes(m.enzymeIndexs_signalPeptidase) * ...
                m.lipoproteinSignalPeptidaseSpecificRate * m.stepSizeSec;
            i = m.secretedMonomerIndexs;
            assertVectorsAlmostEqual(...
                n, sum(m.processedMonomers(i)), 'relative', 0.10);
            assertTrue(10 > max(m.processedMonomers(i)));
            assertTrue(0 < min(m.processedMonomers(i)));
            assertEqual(...
                10 * ones(size(i)), ...
                m.processedMonomers(i) + m.unprocessedMonomers(i));
            assertEqual(m.processedMonomers, m.signalSequenceMonomers);
        end

        % Verifies that signal peptidase processes roughly as many monomers as
        % its rate will allow when it's the limiting factor and there is a mix
        % of lipoproteins and secreted proteins, and that the monomers chosen
        % for processing are chosen without egregious bias.
        function testLimitedSignalPeptidase_lipoproteinsAndSecretedProteins(this)
            m = this.process;
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.substrates(m.substrateIndexs_PG160) = 1e3;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_signalPeptidase) = 20;
            m.enzymes(m.enzymeIndexs_diacylglycerylTransferase) = 1e4;
            m.unprocessedMonomers(:) = 0;
            m.unprocessedMonomers(m.secretedMonomerIndexs) = 10;
            m.unprocessedMonomers(m.lipoproteinMonomerIndexs) = 10;
            m.processedMonomers(:) = 0;
            m.signalSequenceMonomers(:) = 0;

            m.evolveState();
            n = m.enzymes(m.enzymeIndexs_signalPeptidase) * ...
                m.lipoproteinSignalPeptidaseSpecificRate * m.stepSizeSec;
            i = [m.secretedMonomerIndexs;m.lipoproteinMonomerIndexs];
            assertVectorsAlmostEqual(...
                n, sum(m.processedMonomers(i)), 'relative', 0.05);
            assertTrue(10 > max(m.processedMonomers(i)));
            assertTrue(0 < min(m.processedMonomers(i)));
            assertEqual(...
                10 * ones(size(i)), ...
                m.processedMonomers(i) + m.unprocessedMonomers(i));
            assertEqual(m.processedMonomers, m.signalSequenceMonomers);
        end

        % Verifies that diacylglyceryl transferase processes roughly as many
        % lipoproteins as its rate will allow when it's the limiting factor.
        function testLimitedDiacylglycerylTransferase(this)
            m = this.process;
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.substrates(m.substrateIndexs_PG160) = 1e3;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_signalPeptidase) = 1e6;
            m.enzymes(m.enzymeIndexs_diacylglycerylTransferase) = 1e3;
            m.unprocessedMonomers(:) = 0;
            m.unprocessedMonomers(m.lipoproteinMonomerIndexs) = 10;
            m.processedMonomers(:) = 0;
            m.signalSequenceMonomers(:) = 0;

            m.evolveState();
            n = m.enzymes(m.enzymeIndexs_diacylglycerylTransferase) * ...
                m.lipoproteinDiacylglycerylTransferaseSpecificRate * m.stepSizeSec;
            i = m.lipoproteinMonomerIndexs;
            assertElementsAlmostEqual(n, sum(m.processedMonomers(i)), 'absolute', 3);
            assertTrue(10 > max(m.processedMonomers(i)));
            assertTrue(0 < min(m.processedMonomers(i)));
        end

        function testLimitedPG160(this)
            m = this.process;
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.substrates(m.substrateIndexs_PG160) = 100;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_signalPeptidase) = 1e5;
            m.enzymes(m.enzymeIndexs_diacylglycerylTransferase) = 1e5;
            m.unprocessedMonomers(:) = 0;
            m.unprocessedMonomers(m.lipoproteinMonomerIndexs) = 10;
            m.processedMonomers(:) = 0;
            m.signalSequenceMonomers(:) = 0;
            assertTrue(sum(m.unprocessedMonomers(m.lipoproteinMonomerIndexs)) > m.substrates(m.substrateIndexs_PG160));

            m.evolveState();
            i = m.lipoproteinMonomerIndexs;
            assertIn(nnz(m.unprocessedMonomers(m.lipoproteinMonomerIndexs)), [1 Inf]);
            assertEqual(100, sum(m.processedMonomers(i)));
        end

        function testLotsOfEverything(this)
            m = this.process;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.substrates(m.substrateIndexs_PG160) = 1e4;
            m.enzymes(m.enzymeIndexs_signalPeptidase) = 1e3;
            m.enzymes(m.enzymeIndexs_diacylglycerylTransferase) = 1e2;
            i = [m.secretedMonomerIndexs;m.lipoproteinMonomerIndexs];
            m.unprocessedMonomers(i) = randi(100, size(i));
            m.processedMonomers(:) = 0;
            m.signalSequenceMonomers(:) = 0;

            m.evolveState();
            bounds = sort(...
                [m.enzymes(m.enzymeIndexs_signalPeptidase) * ...
                    m.lipoproteinSignalPeptidaseSpecificRate * ...
                    m.stepSizeSec;
                 m.enzymes(m.enzymeIndexs_diacylglycerylTransferase) * ...
                    m.lipoproteinDiacylglycerylTransferaseSpecificRate * ...
                    m.stepSizeSec]);
            assertTrue(0.99 * bounds(1) < sum(m.processedMonomers(i)));
            assertTrue(1.01 * bounds(2) > sum(m.processedMonomers(i)));
        end

        function testGeneEssentiality(this)
            m = this.process;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.substrates(m.substrateIndexs_PG160) = 1e3;
            m.enzymes(m.enzymeIndexs_signalPeptidase) = 1e3;
            m.enzymes(m.enzymeIndexs_diacylglycerylTransferase) = 1e3;
            m.unprocessedMonomers(:) = 1;
            m.processedMonomers(:) = 0;
            this.helpTestGeneEssentiality({
                'MG_086';     %prolipoprotein diacylglyceryl transferase
                'MG_210'},... %prolipoprotein signal peptidase, signal peptidase II
                @(m, i) any(i.processedMonomers(m.unprocessedMonomerIndexs) < ...
                            m.processedMonomers(m.unprocessedMonomerIndexs)) && ...
                        any(i.processedMonomers(m.lipoproteinMonomerIndexs) < ...
                            m.processedMonomers(m.lipoproteinMonomerIndexs)) && ...
                        any(i.processedMonomers(m.secretedMonomerIndexs) < ...
                            m.processedMonomers(m.secretedMonomerIndexs)));
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell