text:
function model=getKEGGModelForOrganism(organismID,fastaFile,dataDir,...
    outDir,keepSpontaneous,keepUndefinedStoich,keepIncomplete,...
    keepGeneral,cutOff,minScoreRatioKO,minScoreRatioG,maxPhylDist,...
    nSequences,seqIdentity,globalModel)
% getKEGGModelForOrganism
%   Reconstructs a genome-scale metabolic model based on protein homology
%   to the orthologies in KEGG. If the target species is not available in
%   KEGG, the user must select a closely related species. It is also
%   possible to circumvent protein homology search (see fastaFile parameter
%   for more details)
%
%   Input:
%   organismID          three or four letter abbreviation of the organism
%                       (as used in KEGG). If not available, use a closely
%                       related species. This is used for determing the
%                       phylogenetic distance. Use 'eukaryotes' or
%                       'prokaryotes' to get a model for the whole domain.
%                       Only applicable if fastaFile is empty, i.e. no
%                       homology search should be performed
%   fastaFile           a FASTA file that contains the protein sequences of
%                       the organism for which to reconstruct a model (opt,
%                       if no FASTA file is supplied then a model is
%                       reconstructed based only on the organism
%                       abbreviation. This option ignores all settings
%                       except for keepSpontaneous, keepUndefinedStoich,
%                       keepIncomplete and keepGeneral)
%   dataDir             directory for which to retrieve the input data.
%                       Should contain a combination of these sub-folders:
%                       -dataDir\keggdb
%                           The KEGG database files used in 1a (see below)
%                       -dataDir\fasta
%                           The multi-FASTA files generated in 1b (see
%                           below)
%                       -dataDir\aligned
%                           The aligned FASTA files as generated in 2a (see
%                           below)
%                       -dataDir\hmms
%                           The hidden Markov models as generated in 2b or
%                           downloaded from BioMet Toolbox (see below)
%                       The final directory in dataDir should be styled as
%                       proXXX_keggYY or eukXXX_keggYY, indicating whether
%                       the HMMs were trained on pro- or eukaryotic
%                       sequences, using a sequence similarity threshold of
%                       XXX %, fitting the KEGG version YY. E.g.
%                       euk90_kegg105. (opt, see note about fastaFile. Note
%                       that in order to rebuild the KEGG model from a
%                       database dump, as opposed to using the version
%                       supplied with RAVEN, you would still need to supply
%                       this)
%   outDir              directory to save the results from the quering of
%                       the hidden Markov models. The output is specific
%                       for the input sequences and the settings used. It
%                       is stored in this manner so that the function can
%                       continue if interrupted or if it should run in
%                       parallel. Be careful not to leave output files from
%                       different organisms or runs with different settings
%                       in the same folder. They will not be overwritten
%                       (opt, default is a temporary dir where all *.out
%                       files are deleted before and after doing the
%                       reconstruction)
%   keepSpontaneous     include reactions labeled as "spontaneous". (opt,
%                       default true)
%   keepUndefinedStoich	include reactions in the form n A <=> n+1 A. These
%                       will be dealt with as two separate metabolites
%                       (opt, default true)
%   keepIncomplete      include reactions which have been labelled as
%                       "incomplete", "erroneous" or "unclear" (opt,
%                       default true)
%   keepGeneral         include reactions which have been labelled as
%                       "general reaction". These are reactions on the form
%                       "an aldehyde <=> an alcohol", and are therefore
%                       unsuited for modelling purposes. Note that not all
%                       reactions have this type of annotation, and the
%                       script will therefore not be able to remove all
%                       such reactions (opt, default false)
%   cutOff              significance score from HMMer needed to assign
%                       genes to a KO (opt, default 10^-50)
%   minScoreRatioG      a gene is only assigned to KOs for which the score
%                       is >=log(score)/log(best score) for that gene. This
%                       is to prevent that a gene which clearly belongs to
%                       one KO is assigned also to KOs with much lower
%                       scores (opt, default 0.8 (lower is less strict))
%   minScoreRatioKO     ignore genes in a KO if their score is
%                       <log(score)/log(best score in KO). This is to
%                       "prune" KOs which have many genes and where some are
%                       clearly a better fit (opt, default 0.3 (lower is
%                       less strict))
%   maxPhylDist         -1: only use sequences from the same domain
%                       (Prokaryota, Eukaryota)
%                       other (positive) value: only use sequences for
%                       organisms where the phylogenetic distance is at the
%                       most this large (as calculated in getPhylDist)
%                       (opt, default Inf, which means that all sequences
%                       will be used)
%   nSequences          for each KO, use up to this many sequences from the
%                       most closely related species. This is mainly to
%                       speed up the alignment process for KOs with very
%                       many genes. This subsampling is performed before
%                       running CD-HIT (opt, default inf)
%   seqIdentity         sequence identity threshold in CD-HIT, referred as
%                       "global sequence identity" in CD-HIT User's Guide.
%                       If -1 is provided, CD-HIT is skipped (opt, default 0.9)
%   globalModel         structure containing both model and KOModel
%                       structures as generated by getModelFromKEGG. These
%                       will otherwise be loaded by via getModelFromKEGG.
%                       Providing globalKEGGmodel can speed up model
%                       generation if getKEGGModelForOrganism is run
%                       multiple times for different strains. Example:
%                       [globalModel.model,globalModel.KOModel] = getModelFromKEGG;
%                       (opt, default empty, global model is loaded by
%                       getModelFromKEGG)
%
%   Output:
%   model               the reconstructed model
%
%   PLEASE READ THIS: The input to this function can be confusing, because
%   it is intended to be run in parallel on a cluster or in multiple
%   sessions. It therefore saves a lot of intermediate results to storage.
%   This also serves the purpose of not having to do redundant
%   calculations. This, however, comes with the disadvantage of somewhat
%   trickier handling. This is what this function does:
%
%   1a. Loads files from a local KEGG FTP dump and constructs a general
%       RAVEN model representing the metabolic network. The functions
%       getRxnsFromKEGG, getGenesFromKEGG, getMetsFromKEGG summarise the
%       data into 'keggRxns.mat', 'keggGenes.mat' and 'keggMets.mat' files,
%       which are later merged into 'keggModel.mat' by getModelFromKEGG
%       function. The function getPhylDist generates 'keggPhylDist.mat'
%       file. KEGG FTP access requires a license.
%   1b. Generates protein FASTA files from the KEGG FTP dump (see 1a). One
%       multi-FASTA file for each KO in KEGG is generated.
%
%   The Step 1 has to be re-done every time KEGG updates their database (or
%   rather when the updates are large enough to warrant re-running this
%   part). Many users would probably never use this feature.
%
%   2a. Filters KO-specific protein sets. This is done by using the
%       settings "maxPhylDist" and "nSequences" to control which sequences
%       should be used for constructing Hidden Markov models (HMMs), and
%       later for matching your sequences to.
%       The most common alternatives here would be to use sequences from
%       only eukaryotes, only prokaryotes or all sequences in KEGG, but you
%       could also play around with the parameters to use e.g. only fungal
%       sequences.
%   2b. KO-specific protein FASTA files are re-organised into
%       non-redundant protein sets with CD-HIT. The user can only set
%       seqIdentity parameter, which corresponds to '-c' parameter in
%       CD-HIT, described as "sequence identity threshold". CD-HIT suggsted
%       sequence identity specific word_length (-n) parameters are used.
%   2c. Does a multi sequence alignment for multi-FASTA files obtained in
%       Step 2b for future use. MAFFT software with automatic selection of
%       alignment algorithm is used in this step ('--auto').
%   2d. Trains hidden Markov models using HMMer for each of the aligned
%       KO-specific FASTA files obtained in Step 2c. This is performed with
%       'hmmbuild' using the default settings.
%
%   Step 2 may be reasonable to be re-done if the user wants to tweak the
%   settings in proteins filtering, clustering, multi sequence alignment or
%   HMMs training steps. However, it requires to have KO-specific protein
%   FASTA files obtained in Step 1a. As such files are not provided in
%   RAVEN and BioMet ToolBox, the user can only generate these files from
%   KEGG FTP dump files, so KEGG FTP license is needed.
%
%   3a. Queries the HMMs with sequences for the organism you are making a
%       model for. This step uses both the output from step 1a and from 2d.
%       This is done with 'hmmsearch' function under default settings. The
%       significance threshold value set in 'cutOff' parameter is used
%       later when parsing '*.out' files to filter out KO hits with higher
%       value than 'cutOff' value. The results with passable E values are
%       summarised into KO-gene occurence matrix with E values in
%       intersections as 'koGeneMat'. The parameters 'minScoreRatioG' and
%       'minScoreRatioKO' are then applied to 'prune' KO-gene associations
%       (see the function descriptions above for more details). The
%       intersection values for these 'prunable' associations are converted
%       to zeroes.
%   3b. Constructs a model based on the pre-processed KO-gene association
%       matrix (koGeneMat). As the full KEGG model already has reaction-KO
%       relationships, KOs are converted into the query genes. The final
%       draft model contains only these reactions, which are associated
%       with KOs from koGeneMat. The reactions without the genes may also
%       be included, if the user set keepSpontaneous as 'true'.
%
%   The Step 3 is specific to the organism for which the model is
%   reconstructed.
%
%   In principle the function looks at which output that is already available
%   and runs only the parts that are required for step 3. This means
%   that (see the definition of the parameters for details):
%   -1a is only performed if there are no KEGG model files in the
%   RAVEN\external\kegg directory
%   -1b is only performed if not all required HMMs OR aligned FASTA files
%   OR multi-FASTA files exist in the defined dataDir. This means that this
%   step is skipped if the HMMs are downloaded from BioMet Toolbox instead
%   (see above). If not all files exist it will try to find
%   the KEGG database files in dataDir.
%   -2a is only performed if not all required HMMs OR aligned FASTA files
%   files exist in the defined dataDir. This means that this step is skipped
%   if the HMMs are downloaded from BioMet Toolbox instead (see above).
%   -2b is only performed if not all required HMMs exist in the defined
%   dataDir. This means that this step is skipped if the FASTA files or
%   HMMs are downloaded from BioMet Toolbox instead (see above).
%   -3a is performed for the required HMMs for which no corresponding .out
%   file exists in outDir. This is just a way to enable the function to be
%   run in parallel or to resume if interrupted.
%   -3b is always performed.
%
%   These steps are specific to the organism for which you are
%   reconstructing the model.
%
%   Regarding the whole pipeline, the function checks the output that is
%   already available and runs only the parts that are required for step 3.
%   This means that (see the definition of the parameters for details):
%   -1a is only performed if there are no KEGG model files in the
%   RAVEN\external\kegg directory.
%   -1b is only performed if any of required KOs do not have HMMs, aligned
%   FASTA files, clustered FASTA files and raw FASTA files in the defined
%   dataDir. This means that this step is skipped if the HMMs are
%   downloaded from BioMet Toolbox instead (see above). If not all files
%   exist it will try to find the KEGG database files in dataDir.
%   -2ab are only performed if any of required KOs do not have HMMs,
%   aligned FASTA files and clustered FASTA files in the defined dataDir.
%   This means that this step is skipped if the HMMs are downloaded from
%   BioMet Toolbox instead (see above).
%   -2c is only performed if any of required KOs do not have HMMs and
%   aligned FASTA files in the defined dataDir. This means that this step
%   is skipped if the HMMs are downloaded from BioMet Toolbox instead (see
%   above).
%   -2d is only performed if any of required KOs do not have HMMs exist in
%   the defined dataDir. This means that this step is skipped if the FASTA
%   files or HMMs are downloaded from BioMet Toolbox instead (see above).
%   -3a is performed for the required HMMs for which no corresponding .out
%   file exists in outDir. This is just a way to enable the function to be
%   run in parallel or to resume if interrupted.
%   -3b is always performed.
%
%   NOTE: it is also possible to obtain draft model from KEGG without
%   providing protein FASTA file for the target organism. In such case the
%   organism three-four letter abbreviation set as 'organismID' must exist
%   in the local KEGG database. In such case, the program just fetches all
%   the reactions, which are associated with given 'organismID'.
%
%   Usage: model=getKEGGModelForOrganism(organismID,fastaFile,dataDir,...
%    outDir,keepSpontaneous,keepUndefinedStoich,keepIncomplete,...
%    keepGeneral,cutOff,minScoreRatioKO,minScoreRatioG,maxPhylDist,...
%    nSequences,seqIdentity)

if nargin<2 || isempty(fastaFile)
    fastaFile=[];
else
    fastaFile=char(fastaFile);
end
if nargin<3
    dataDir=[];
else
    dataDir=char(dataDir);
end
if nargin<4 || isempty(outDir)
    outDir=tempdir;
    %Delete all *.out files if any exist
    delete(fullfile(outDir,'*.out'));
else
    outDir=char(outDir);
end
if nargin<5
    keepSpontaneous=true;
end
if nargin<6
    keepUndefinedStoich=true;
end
if nargin<7
    keepIncomplete=true;
end
if nargin<8
    keepGeneral=false;
end
if nargin<9
    cutOff=10^-50;
end
if nargin<10
    minScoreRatioKO=0.3;
end
if nargin<11
    minScoreRatioG=0.8;
end
if nargin<12
    maxPhylDist=inf;
    %Include all sequences for each reaction
end
if nargin<13
    nSequences=inf;
    %Include all sequences for each reaction
end
if nargin<14
    seqIdentity=0.9;
end

if isempty(fastaFile)
    fprintf(['\n*** The model reconstruction from KEGG based on the annotation available for KEGG Species ' organismID ' ***\n\n']);
else
    fprintf('\n*** The model reconstruction from KEGG based on the protein homology search against KEGG Orthology specific HMMs ***\n\n');
    %Check if query fasta exists
    fastaFile=checkFileExistence(fastaFile,2); %Copy file to temp dir
end

%Run the external binaries multi-threaded to use all logical cores assigned
%to MATLAB
cores = evalc('feature(''numcores'')');
cores = strsplit(cores, 'MATLAB was assigned: ');
cores = regexp(cores{2},'^\d*','match');
cores = cores{1};

%Get the directory for RAVEN Toolbox.
ravenPath=findRAVENroot();

%Checking if dataDir is consistent. It must point to pre-trained HMMs set,
%compatible with the the current RAVEN version. The user may have the
%required zip file already in working directory or have it extracted. If
%the zip file and directory is not here, it is downloaded from the cloud
if ~isempty(dataDir)
    hmmOptions={'euk90_kegg105','prok90_kegg105'};
    if ~endsWith(dataDir,hmmOptions) %Check if dataDir ends with any of the hmmOptions.
                                     %If not, then check whether the required folders exist anyway.
        if ~isfile(fullfile(dataDir,'keggdb','genes.pep')) && ...
                ~isfolder(fullfile(dataDir,'fasta')) && ...
                ~isfolder(fullfile(dataDir,'aligned')) && ...
                ~isfolder(fullfile(dataDir,'hmms'))
            error(['Pre-trained HMMs set is not recognised. If you want download RAVEN provided sets, it should match any of the following: ' strjoin(hmmOptions,' or ')])
        end
    else
        if isfolder(dataDir) && isfile(fullfile(dataDir,'hmms','K00844.hmm'))
            fprintf(['NOTE: Found ' dataDir ' directory with pre-trained HMMs, it will therefore be used during reconstruction\n']);
        elseif ~isfolder(dataDir) && isfile([dataDir,'.zip'])
            fprintf('Extracting the HMMs archive file... ');
            unzip([dataDir,'.zip']);
            fprintf('COMPLETE\n');
        else
            hmmIndex=strcmp(dataDir,hmmOptions);
            if ~any(hmmIndex)
                error(['Pre-trained HMMs are only provided with proteins clustered at 90% sequence identity (i.e. prok90_kegg105 and euk90_kegg105). ' ...
                    'Use either of these datasets, or otherwise download the relevant sequence data from KEGG to train HMMs with your desired sequence identity'])
            else
                fprintf('Downloading the HMMs archive file... ');
                try
                    websave([dataDir,'.zip'],['https://github.com/SysBioChalmers/RAVEN/releases/download/v2.8.0/',hmmOptions{hmmIndex},'.zip']);
                catch ME
                    if strcmp(ME.identifier,'MATLAB:webservices:HTTP404StatusCodeError')
                        error('Failed to download the HMMs archive file, the server returned a 404 error, try again later. If the problem persists please report it on the RAVEN GitHub Issues page: https://github.com/SysBioChalmers/RAVEN/issues')
                    end
                end
            end

            fprintf('COMPLETE\n');
            fprintf('Extracting the HMMs archive file... ');
            unzip([dataDir,'.zip']);
            fprintf('COMPLETE\n');
        end
        %Check if HMMs are extracted
        if ~isfile(fullfile(dataDir,'hmms','K00844.hmm'))
            error(['The HMM files seem improperly extracted and not found in ',dataDir,'/hmms. Please remove ',dataDir,' folder and rerun getKEGGModelForOrganism']);
        end
    end
end

%Check if the fasta-file contains '/' or'\'. If not then it's probably just
%a file name. Expand to full path.
if any(fastaFile)
    if ~any(strfind(fastaFile,'\')) && ~any(strfind(fastaFile,'/'))
        fastaFile=which(fastaFile);
    end
    %Create the required sub-folders in dataDir if they dont exist
    if ~isfolder(fullfile(dataDir,'keggdb'))
        mkdir(dataDir,'keggdb');
    end
    if ~isfolder(fullfile(dataDir,'fasta'))
        mkdir(dataDir,'fasta');
    end
    if ~isfolder(fullfile(dataDir,'aligned'))
        mkdir(dataDir,'aligned');
    end
    if ~isfolder(fullfile(dataDir,'hmms'))
        mkdir(dataDir,'hmms');
    end
    if ~isfolder(outDir)
        mkdir(outDir);
    end
end

%First generate the full global KEGG model. Can be provided as input.
%Otherwise, getModelFromKEGG is run. The dataDir must not be supplied as
%there is also an internal RAVEN version available
if nargin==15
    model=globalModel.model;
    KOModel=globalModel.KOModel;
elseif any(dataDir)
    [model, KOModel]=getModelFromKEGG(fullfile(dataDir,'keggdb'),keepSpontaneous,keepUndefinedStoich,keepIncomplete,keepGeneral);
else
    [model, KOModel]=getModelFromKEGG([],keepSpontaneous,keepUndefinedStoich,keepIncomplete,keepGeneral);
end
model.id=organismID;
model.c=zeros(numel(model.rxns),1);

%If no FASTA file is supplied, then just remove all genes which are not for
%the given organism ID
if isempty(fastaFile)
    %Check if organismID can be found in KEGG species list or is
    %set to "eukaryotes" or "prokaryotes"
    phylDistsFull=getPhylDist(fullfile(dataDir,'keggdb'),true);
    if ~ismember(organismID,[phylDistsFull.ids 'eukaryotes' 'prokaryotes'])
        error('Provided organismID is incorrect. Only species abbreviations from KEGG Species List or "eukaryotes"/"prokaryotes" are allowed.');
    end

    fprintf(['Pruning the model from non-' organismID ' genes... ']);
    if ismember(organismID,{'eukaryotes','prokaryotes'})
        phylDists=getPhylDist(fullfile(dataDir,'keggdb'),maxPhylDist==-1);
        if strcmp(organismID,'eukaryotes')
            proxyid='hsa';
            %Use H. sapiens here
        else
            proxyid='eco';
            %Use E. coli here
        end
        [~, phylDistId]=ismember(proxyid,phylDists.ids);
        idsToKeep=phylDists.ids(~isinf(phylDists.distMat(phylDistId,:)));
        taxIDs=cellfun(@(x) x{1},cellfun(@(x) strsplit(x,':'),model.genes,'UniformOutput',false),'UniformOutput',false);
        I=ismember(upper(taxIDs),upper(idsToKeep));
    else
        %KEGG organism IDs may have three or four letters
        organismID=strcat(organismID,':');
        %Add colon for accurate matching
        if length(organismID)==4
            I=cellfun(@(x) strcmpi(x(1:4),organismID),model.genes);
        elseif length(organismID)==5
            I=cellfun(@(x) strcmpi(x(1:5),organismID),model.genes);
        end
    end
    %Remove those genes
    model.genes=model.genes(I);
    model.rxnGeneMat=model.rxnGeneMat(:,I);
    fprintf('COMPLETE\n');
end

%First remove all reactions without genes
if keepSpontaneous==true
    fprintf('Removing non-spontaneous reactions without GPR rules... ');
    load(fullfile(ravenPath,'external','kegg','keggRxns.mat'),'isSpontaneous');
    I=~any(model.rxnGeneMat,2)&~ismember(model.rxns,isSpontaneous);
    spontRxnsWithGenes=model.rxns(any(model.rxnGeneMat,2)&~ismember(model.rxns,isSpontaneous));
else
    fprintf('Removing reactions without GPR rules... ');
    I=~any(model.rxnGeneMat,2);
end
model=removeReactions(model,I,true);
fprintf('COMPLETE\n');

%Clean gene names
fprintf('Fixing gene names in the model... ');
%Get rid of the prefix organism id
model.genes=regexprep(model.genes,'^\w+?:','');
fprintf('COMPLETE\n');

%If no FASTA file is supplied, then we are done here
if isempty(fastaFile)
    %Create grRules
    fprintf('Constructing GPR associations and annotations for the model... ');
    model.grRules=cell(numel(model.rxns),1);
    model.grRules(:)={''};
    %Add the gene associations as 'or'
    for i=1:numel(model.rxns)
        %Find the involved genes
        I=find(model.rxnGeneMat(i,:));
        if any(I)
            model.grRules{i}=['(' model.genes{I(1)}];
            for j=2:numel(I)
                model.grRules{i}=[model.grRules{i} ' or ' model.genes{I(j)}];
            end
            model.grRules{i}=[model.grRules{i} ')'];
        end
    end
    %Fix grRules and reconstruct rxnGeneMat
    [grRules,rxnGeneMat] = standardizeGrRules(model); %Give detailed output
    model.grRules = grRules;
    model.rxnGeneMat = rxnGeneMat;
    %Add geneMiriams, assuming that it follows the syntax
    %kegg.genes/organismID:geneName
    model.geneMiriams='';
    for i=1:numel(model.genes)
        model.geneMiriams{i,1}.name{1,1}='kegg.genes';
        model.geneMiriams{i,1}.value{1,1}=strcat(lower(organismID),model.genes{i,1});
    end
    %Add the description to the reactions
    for i=1:numel(model.rxns)
        if ~isempty(model.rxnNotes{i})
            model.rxnNotes(i)=strcat('Included by getKEGGModelForOrganism (without HMMs).',model.rxnNotes(i));
            model.rxnNotes(i)=strrep(model.rxnNotes(i),'.','. ');
        else
            model.rxnNotes(i)={'Included by getKEGGModelForOrganism (without HMMs)'};
        end
    end
    fprintf('COMPLETE\n\n');
    fprintf('*** Model reconstruction complete ***\n');
    return;
end

%Create a phylogenetic distance structure
phylDistStruct=getPhylDist(fullfile(dataDir,'keggdb'),maxPhylDist==-1);
[~, phylDistId]=ismember(model.id,phylDistStruct.ids);

%Calculate the real maximal distance now. An abitary large number of 1000
%is used for the "all in kingdom" or "all sequences" options. This is a bit
%inconvenient way to do it, but it is to make it fit with some older code
if isinf(maxPhylDist) || maxPhylDist==-1
    maxPhylDist=1000;
end

%Get the KO ids for which files have been generated. Maybe not the neatest
%way..
fastaFiles=listFiles(fullfile(dataDir,'fasta','*.fa'));
alignedFiles=listFiles(fullfile(dataDir,'aligned','*.fa'));
alignedWorking=listFiles(fullfile(dataDir,'aligned','*.faw'));
hmmFiles=listFiles(fullfile(dataDir,'hmms','*.hmm'));
outFiles=listFiles(fullfile(outDir,'*.out'));

%Check if multi-FASTA files should be generated. This should only be
%performed if there are IDs in the KOModel structure that haven't been
%parsed yet
missingFASTA=setdiff(KOModel.rxns,[fastaFiles;alignedFiles;hmmFiles;outFiles]);

if ~isempty(missingFASTA)
    if ~isfile(fullfile(dataDir,'keggdb','genes.pep'))
        EM=['The file ''genes.pep'' cannot be located at ' strrep(dataDir,'\','/') '/ and should be downloaded from the KEGG FTP.\n'];
        dispEM(EM);
    end
    %Only construct models for KOs which don't have files already
    fastaModel=removeReactions(KOModel,setdiff(KOModel.rxns,missingFASTA),true,true);
    %Permute the order of the KOs in the model so that constructMultiFasta
    %can be run on several processors at once
    fastaModel=permuteModel(fastaModel,randperm(RandStream.create('mrg32k3a','Seed',cputime()),numel(fastaModel.rxns)),'rxns');
    constructMultiFasta(fastaModel,fullfile(dataDir,'keggdb','genes.pep'),fullfile(dataDir,'fasta'));
else
    fprintf('Generating the KEGG Orthology specific multi-FASTA files... COMPLETE\n');
end

if isunix
    if ismac
        binEnd='.mac';
    else
        binEnd='';
    end
elseif ispc
    binEnd='';
else
    EM='Unknown OS, exiting.';
    disp(EM);
    return
end

%Check if alignment of FASTA files should be performed
missingAligned=setdiff(KOModel.rxns,[alignedFiles;hmmFiles;alignedWorking;outFiles]);
if ~isempty(missingAligned)
    if seqIdentity==-1
        fprintf('Performing the multiple alignment for KEGG Orthology specific protein sets...   0%% complete');
    else
        fprintf('Performing clustering and multiple alignment for KEGG Orthology specific protein sets...   0%% complete');
    end
    missingAligned=missingAligned(randperm(RandStream.create('mrg32k3a','Seed',cputime()),numel(missingAligned)));
    tmpFile=tempname;
    %On Windows, paths need to be translated to Unix before parsing it to WSL
    if ispc
        wslPath.tmpFile=getWSLpath(tmpFile);
        %mafft has problems writing to terminal (/dev/stderr) when running
        %on WSL via MATLAB, instead write and read progress file
        mafftOutput = tempname;
        wslPath.mafftOutput=getWSLpath(mafftOutput);
        wslPath.mafft=getWSLpath(fullfile(ravenPath,'software','mafft','mafft-linux64','mafft.bat'));
        wslPath.cdhit=getWSLpath(fullfile(ravenPath,'software','cd-hit','cd-hit'));
    end

    for i=1:numel(missingAligned)
        %This is checked here because it could be that it is created by a
        %parallel process. The faw-files are saved as temporary files to
        %kept track of which files are being worked on
        if ~isfile(fullfile(dataDir,'aligned',[missingAligned{i} '.faw'])) &&...
                ~isfile(fullfile(dataDir,'aligned',[missingAligned{i} '.fa']))
            %Check that the multi-FASTA file exists. It should do so since
            %we are saving empty files as well. Print a warning and
            %continue if not
            if ~isfile(fullfile(dataDir,'fasta',[missingAligned{i} '.fa']))
                EM=['WARNING: The multi-FASTA file for ' missingAligned{i} ' does not exist'];
                dispEM(EM,false);
                continue;
            end

            %If the multi-FASTA file is empty then save an empty aligned
            %file and continue
            s=dir(fullfile(dataDir,'fasta',[missingAligned{i} '.fa']));
            if s.bytes<=0
                fid=fopen(fullfile(dataDir,'aligned',[missingAligned{i} '.fa']),'w');
                fclose(fid);
                continue;
            end

            %Create an empty file to prevent other threads to start to work
            %on the same alignment
            fid=fopen(fullfile(dataDir,'aligned',[missingAligned{i} '.faw']),'w');
            fclose(fid);

            %First load the FASTA file, then select up to nSequences
            %sequences of the most closely related species, apply any
            %constraints from maxPhylDist, and save it as a temporary file,
            %and create the model from that

            fastaStruct=fastaread(fullfile(dataDir,'fasta',[missingAligned{i} '.fa']));
            phylDist=inf(numel(fastaStruct),1);
            for j=1:numel(fastaStruct)
                %Get the organism abbreviation
                index=strfind(fastaStruct(j).Header,':');
                if any(index)
                    abbrev=fastaStruct(j).Header(1:index(1)-1);
                    [~, index]=ismember(abbrev,phylDistStruct.ids);
                    if any(index)
                        phylDist(j)=phylDistStruct.distMat(index(1),phylDistId);
                    end
                end
            end

            %Inf means that it should not be included
            phylDist(phylDist>maxPhylDist)=[];

            %Sort based on phylDist
            [~, order]=sort(phylDist);

            %Save the first nSequences hits to a temporary FASTA file
            if nSequences<=numel(fastaStruct)
                fastaStruct=fastaStruct(order(1:nSequences));
            else
                fastaStruct=fastaStruct(order);
            end

            %Do the clustering and alignment if there are more than one
            %sequences, otherwise just save the sequence (or an empty file)
            if numel(fastaStruct)>1
                 if seqIdentity~=-1
                    cdhitInpCustom=tempname;
                    fastawrite(cdhitInpCustom,fastaStruct);
                    if seqIdentity<=1 && seqIdentity>0.7
                        nparam='5';
                    elseif seqIdentity>0.6
                        nparam='4';
                    elseif seqIdentity>0.5
                        nparam='3';
                    elseif seqIdentity>0.4
                        nparam='2';
                    else
                        EM='The provided seqIdentity must be between 0 and 1\n';
                        dispEM(EM);
                    end
                    if ispc
                        wslPath.cdhitInpCustom=getWSLpath(cdhitInpCustom);
                        [status, output]=system(['wsl "' wslPath.cdhit '" -T "' num2str(cores) '" -i "' wslPath.cdhitInpCustom '" -o "' wslPath.tmpFile '" -c "' num2str(seqIdentity) '" -n ' nparam ' -M 2000']);
                    elseif ismac || isunix
                        [status, output]=system(['"' fullfile(ravenPath,'software','cd-hit',['cd-hit' binEnd]) '" -T "' num2str(cores) '" -i "' cdhitInpCustom '" -o "' tmpFile '" -c "' num2str(seqIdentity) '" -n ' nparam ' -M 2000']);
                    end
                    if status~=0
                        EM=['Error when performing clustering of ' missingAligned{i} ':\n' output];
                        dispEM(EM);
                    end
                    %Remove the old tempfile
                    if exist(cdhitInpCustom, 'file')
                        delete([cdhitInpCustom '*']);
                    end
                else
                    %This means that CD-HIT should be skipped since
                    %seqIdentity is equal to -1
                    fastawrite(tmpFile,fastaStruct);
                end
                %Do the alignment for this file
                if ismac
                    [status, output]=system(['"' fullfile(ravenPath,'software','mafft','mafft-mac','mafft.bat') '" --auto --anysymbol --thread "' num2str(cores) '" "' tmpFile '" > "' fullfile(dataDir,'aligned',[missingAligned{i} '.faw']) '"']);
                elseif isunix
                    [status, output]=system(['"' fullfile(ravenPath,'software','mafft','mafft-linux64','mafft.bat') '" --auto --anysymbol --thread "' num2str(cores) '" "' tmpFile '" > "' fullfile(dataDir,'aligned',[missingAligned{i} '.faw']) '"']);
                elseif ispc
                    wslPath.fawFile=getWSLpath(fullfile(dataDir,'aligned',[missingAligned{i} '.faw']));
                    [status, ~]=system(['wsl "' wslPath.mafft '" --auto --anysymbol --progress "' wslPath.mafftOutput '" --thread "' num2str(cores) '" --out "' wslPath.fawFile '" "' wslPath.tmpFile '"']);
                    output=fileread(mafftOutput);
                    delete(mafftOutput);
                end
                if status~=0
                    %It could be that alignment failed because only one
                    %sequence was left after clustering. If that is the
                    %case, then the clustered file is just copied as 'faw'
                    %file
                    if any(regexp(output,'Only 1 sequence found'))
                        movefile(tmpFile,fullfile(dataDir,'aligned',[missingAligned{i} '.faw']),'f');
                    else
                        EM=['Error when performing alignment of ' missingAligned{i} ':\n' output];
                        dispEM(EM);
                    end
                end
                %Remove the old tempfile
                if exist(tmpFile, 'file')
                    delete([tmpFile '*']);
                end
            else
                %If there is only one sequence then it's not possible to do
                %a multiple alignment. Just print the sequence instead. An
                %empty file was written previously so that doesn't have to
                %be dealt with
                if numel(fastaStruct)==1
                    warnState = warning; %Save the current warning state
                    warning('off','Bioinfo:fastawrite:AppendToFile');
                    fastawrite(fullfile(dataDir,'aligned',[missingAligned{i} '.faw']),fastaStruct);
                    warning(warnState) %Reset warning state to previous settings
                end
            end
            %Move the temporary file to the real one
            movefile(fullfile(dataDir,'aligned',[missingAligned{i} '.faw']),fullfile(dataDir,'aligned',[missingAligned{i} '.fa']),'f');

            %Print the progress every 25 files
            if rem(i-1,25) == 0
                progress=num2str(floor(100*numel(listFiles(fullfile(dataDir,'aligned','*.fa')))/numel(KOModel.rxns)));
                progress=pad(progress,3,'left');
                fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress);
            end
        end
    end
    fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n');
else
    if seqIdentity==-1
        fprintf('Performing the multiple alignment for KEGG Orthology specific protein sets... COMPLETE\n');
    else
        fprintf('Performing clustering and multiple alignment for KEGG Orthology specific protein sets... COMPLETE\n');
    end
end

%Check if training of Hidden Markov models should be performed
missingHMMs=setdiff(KOModel.rxns,[hmmFiles;outFiles]);
if ~isempty(missingHMMs)
    fprintf('Training the KEGG Orthology specific HMMs...   0%% complete');
    missingHMMs=missingHMMs(randperm(RandStream.create('mrg32k3a','Seed',cputime()),numel(missingHMMs)));
    %Train models for all missing KOs
    for i=1:numel(missingHMMs)
        %This is checked here because it could be that it is created by a
        %parallel process
        if ~isfile(fullfile(dataDir,'hmms',[missingHMMs{i} '.hmm'])) && ~isfile(fullfile(dataDir,'hmms',[missingHMMs{i} '.hmw']))
            %Check that the aligned FASTA file exists. It could be that it
            %is still being worked on by some other instance of the program
            %(the .faw file should then exist). This should not happen on a
            %single computer. It doesn't throw an error, because it should
            %finalize the ones it can
            if ~isfile(fullfile(dataDir,'aligned',[missingHMMs{i} '.fa']))
                EM=['The aligned FASTA file for ' missingHMMs{i} ' does not exist'];
                dispEM(EM,false);
                continue;
            end

            %If the multi-FASTA file is empty then save an empty aligned
            %file and continue
            s=dir(fullfile(dataDir,'aligned',[missingHMMs{i} '.fa']));
            if s.bytes<=0
                fid=fopen(fullfile(dataDir,'hmms',[missingHMMs{i} '.hmm']),'w');
                fclose(fid);
                continue;
            end
            %Create a temporary file to indicate that it is working on the
            %KO. This is because hmmbuild cannot overwrite existing files
            fid=fopen(fullfile(dataDir,'hmms',[missingHMMs{i} '.hmw']),'w');
            fclose(fid);

            %Create HMM
            [status, output]=system(['"' fullfile(ravenPath,'software','hmmer',['hmmbuild' binEnd]) '" --cpu "' num2str(cores) '" "' fullfile(dataDir,'hmms',[missingHMMs{i} '.hmm']) '" "' fullfile(dataDir,'aligned',[missingHMMs{i} '.fa']) '"']);
            if status~=0
                EM=['Error when training HMM for ' missingHMMs{i} ':\n' output];
                dispEM(EM);
            end

            %Delete the temporary file
            delete(fullfile(dataDir,'hmms',[missingHMMs{i} '.hmw']));

            %Print the progress every 25 files
            if rem(i-1,25) == 0
                progress=num2str(floor(100*numel(listFiles(fullfile(dataDir,'hmms','*.hmm')))/numel(KOModel.rxns)));
                progress=pad(progress,3,'left');
                fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress);
            end
        end
    end
    fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n');
else
    fprintf('Training the KEGG Orthology specific HMMs... COMPLETE\n');
end

%Check which new .out files that should be generated. Check if training of
%Hidden Markov models should be performed
missingOUT=setdiff(KOModel.rxns,outFiles);
if ~isempty(missingOUT)
    fprintf('Querying the user-specified FASTA file against the KEGG Orthology specific HMMs...   0%% complete');
    missingOUT=missingOUT(randperm(RandStream.create('mrg32k3a','Seed',cputime()),numel(missingOUT)));
    for i=1:numel(missingOUT)
        %This is checked here because it could be that it is created by a
        %parallel process
        if ~isfile(fullfile(outDir,[missingOUT{i} '.out']))
            %Check that the HMM file exists. It should do so since %we are
            %saving empty files as well. Print a warning and continue if
            %not
            if ~isfile(fullfile(dataDir,'hmms',[missingOUT{i} '.hmm']))
                EM=['The HMM file for ' missingOUT{i} ' does not exist'];
                dispEM(EM,false);
                continue;
            end

            %Save an empty file to prevent several threads working on the
            %same file
            fid=fopen(fullfile(outDir,[missingOUT{i} '.out']),'w');
            fclose(fid);

            %If the HMM file is empty then save an out file and continue
            s=dir(fullfile(dataDir,'hmms',[missingOUT{i} '.hmm']));
            if s.bytes<=0
                continue;
            end

            %Check each gene in the input file against this model
            [status, output]=system(['"' fullfile(ravenPath,'software','hmmer',['hmmsearch' binEnd]) '" --cpu "' num2str(cores) '" "' fullfile(dataDir,'hmms',[missingOUT{i} '.hmm']) '" "' fastaFile '"']);
            if status~=0
                EM=['Error when querying HMM for ' missingOUT{i} ':\n' output];
                dispEM(EM);
            end

            %Save the output to a file
            fid=fopen(fullfile(outDir,[missingOUT{i} '.out']),'w');
            fwrite(fid,output);
            fclose(fid);

            %Print the progress every 25 files
            if rem(i-1,25) == 0
                progress=num2str(floor(100*numel(listFiles(fullfile(outDir,'*.out')))/numel(KOModel.rxns)));
                progress=pad(progress,3,'left');
                fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress);
            end
        end
    end
    fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n');
else
    fprintf('Querying the user-specified FASTA file against the KEGG Orthology specific HMMs... COMPLETE\n');
end

%***Begin retrieving the output and putting together the resulting model

fprintf('Parsing the HMM search results... ');
%Retrieve matched genes from the HMMs
koGeneMat=zeros(numel(KOModel.rxns),3000); %Make room for 3000 genes
genes=cell(3000,1);
%Store the best score for a gene in a hash list (since it will be searching
%many times)
hTable = java.util.Hashtable;

geneCounter=0;
for i=1:numel(KOModel.rxns)
    if exist(fullfile(outDir,[KOModel.rxns{i} '.out']), 'file')
        fid=fopen(fullfile(outDir,[KOModel.rxns{i} '.out']),'r');
        beginMatches=false;
        while 1
            %Get the next line
            tline = fgetl(fid);

            %Abort at end of file
            if ~ischar(tline)
                break;
            end

            if and(beginMatches,strcmp(tline,'  ------ inclusion threshold ------'))
                break;
            end

            if beginMatches==false
                %This is how the listing of matches begins
                if any(strfind(tline,'E-value '))
                    %Read one more line that is only padding
                    tline = fgetl(fid);
                    beginMatches=true;
                end
            else
                %If matches should be read
                if ~strcmp(tline,'   [No hits detected that satisfy reporting thresholds]') && ~isempty(tline)
                    elements=regexp(tline,' ','split');
                    elements=elements(cellfun(@any,elements));

                    %Check if the match is below the treshhold
                    score=str2double(elements{1});
                    gene=elements{9};
                    if score<=cutOff
                        %If the score is exactly 0, change it to a very
                        %small value to avoid NaN
                        if score==0
                            score=10^-250;
                        end
                        %Check if the gene is added already and, is so, get
                        %the best score for it
                        I=hTable.get(gene);
                        if any(I)
                            koGeneMat(i,I)=score;
                        else
                            geneCounter=geneCounter+1;
                            %The gene was not present yet so add it
                            hTable.put(gene,geneCounter);
                            genes{geneCounter}=gene;
                            koGeneMat(i,geneCounter)=score;
                        end
                    end
                else
                    break;
                end
            end
        end
        fclose(fid);
    end
end
fprintf('COMPLETE\n');

fprintf('Removing gene, KEGG Orthology associations below minScoreRatioKO, minScoreRatioG... ');
koGeneMat=koGeneMat(:,1:geneCounter);

%Remove the genes for each KO that are below minScoreRatioKO.
for i=1:size(koGeneMat,1)
    J=find(koGeneMat(i,:));
    if any(J)
        koGeneMat(i,J(log(koGeneMat(i,J))/log(min(koGeneMat(i,J)))<minScoreRatioKO))=0;
    end
end

%Remove the KOs for each gene that are below minScoreRatioG
for i=1:size(koGeneMat,2)
    J=find(koGeneMat(:,i));
    if any(J)
        koGeneMat(J(log(koGeneMat(J,i))/log(min(koGeneMat(J,i)))<minScoreRatioG),i)=0;
    end
end
fprintf('COMPLETE\n');

fprintf('Adding gene annotations to the model... ');
%Create the new model
model.genes=genes(1:geneCounter);
model.grRules=cell(numel(model.rxns),1);
model.grRules(:)={''};
model.rxnGeneMat=sparse(numel(model.rxns),numel(model.genes));

%Loop through the reactions and add the corresponding genes
for i=1:numel(model.rxns)
    if isstruct(model.rxnMiriams{i})
        %Get all KOs
        I=find(strcmpi(model.rxnMiriams{i}.name,'kegg.orthology'));
        KOs=model.rxnMiriams{i}.value(I);
        %Find the KOs and the corresponding genes
        J=ismember(KOModel.rxns,KOs);
        [~, K]=find(koGeneMat(J,:));

        if any(K)
            model.rxnGeneMat(i,K)=1;
            %Also delete KOs for which no genes were found. If no genes at
            %all were matched to the reaction it will be deleted later
            L=sum(koGeneMat(J,:),2)==0;
            model.rxnMiriams{i}.value(I(L))=[];
            model.rxnMiriams{i}.name(I(L))=[];
        end
    end
end
fprintf('COMPLETE\n');

%Find and delete all reactions without genes. This also removes genes that
%are not used (which could happen because minScoreRatioG and
%minScoreRatioKO). If keepSpontaneous==true, the spontaneous reactions
%without genes are kept in the model. Spontaneous reactions with original
%gene associations are treated in the same way, like the rest of the
%reactions - if gene associations were removed during HMM search, such
%reactions are deleted from the model
if keepSpontaneous==true
    %Not the most comprise way to delete reactions without genes, but this
    %makes the code easier to understand. Firstly the non-spontaneous
    %reactions without genes are removed. After that, the second deletion
    %step removes spontaneous reactions, which had gene associations before
    %HMM search, but no longer have after it
    fprintf('Removing non-spontaneous reactions which after HMM search no longer have GPR rules... ');
    I=~any(model.rxnGeneMat,2)&~ismember(model.rxns,isSpontaneous);
    model=removeReactions(model,I,true,true);
    I=~any(model.rxnGeneMat,2)&ismember(model.rxns,spontRxnsWithGenes);
    model=removeReactions(model,I,true,true);
else
    %Just simply check for any new reactions without genes and remove
    %it
    fprintf('Removing reactions which after HMM search no longer have GPR rules... ');
    I=~any(model.rxnGeneMat,2);
    model=removeReactions(model,I,true,true);
end
fprintf('COMPLETE\n');

fprintf('Constructing GPR rules and finalizing the model... ');
%Add the gene associations as 'or'
for i=1:numel(model.rxns)
    %Find the involved genes
    I=find(model.rxnGeneMat(i,:));
    if any(I)
        model.grRules{i}=['(' model.genes{I(1)}];
        for j=2:numel(I)
            model.grRules{i}=[model.grRules{i} ' or ' model.genes{I(j)}];
        end
        model.grRules{i}=[model.grRules{i} ')'];
    end
end

%Fix grRules and reconstruct rxnGeneMat
[grRules,rxnGeneMat] = standardizeGrRules(model,false); %Give detailed output
model.grRules = grRules;
model.rxnGeneMat = rxnGeneMat;

%Add the description to the reactions
for i=1:numel(model.rxns)
    if ~isempty(model.rxnNotes{i})
        model.rxnNotes(i)=strcat('Included by getKEGGModelForOrganism (using HMMs).',model.rxnNotes(i));
        model.rxnNotes(i)=strrep(model.rxnNotes(i),'.','. ');
    else
        model.rxnNotes(i)={'Included by getKEGGModelForOrganism (using HMMs)'};
    end
end
%Remove the temp fasta file
delete(fastaFile)
fprintf('COMPLETE\n\n*** Model reconstruction complete ***\n');
end

function files=listFiles(directory)
%Supporter function to list the files in a directory and return them as a
%cell array
temp=dir(directory);
files=cell(numel(temp),1);
for i=1:numel(temp)
    files{i}=temp(i,1).name;
end
files=strrep(files,'.fa','');
files=strrep(files,'.hmm','');
files=strrep(files,'.out','');
files=strrep(files,'.faw','');
end

author: SysBioChalmers

repo: RAVEN


text:
% eeg_eegrej() - reject porition of continuous data in an EEGLAB
%                dataset
%
% Usage:
%   >> EEGOUT = eeg_eegrej( EEGIN, regions );
%
% Inputs:
%   INEEG      - input dataset
%   regions    - array of regions to suppress. number x [beg end]  of
%                regions. 'beg' and 'end' are expressed in term of points
%                in the input dataset. Size of the array is
%                number x 2 of regions.
%
% Outputs:
%   INEEG      - output dataset with updated data, events latencies and
%                additional boundary events.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 8 August 2002
%
% See also: eeglab(), eegplot(), pop_rejepoch()

% Copyright (C) 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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 [EEG, com] = eeg_eegrej( EEG, regions);

com = '';
if nargin < 2
	help eeg_eegrej;
	return;
end;
if nargin<3
    probadded = [];
end
if isempty(regions)
	return;
end;

% For AMICA probabilities...Temporarily add model probabilities as channels
%-----------------------------------------------------
if isfield(EEG.etc, 'amica') && ~isempty(EEG.etc.amica) && isfield(EEG.etc.amica, 'v_smooth') && ~isempty(EEG.etc.amica.v_smooth) && ~isfield(EEG.etc.amica,'prob_added')
    if isfield(EEG.etc.amica, 'num_models') && ~isempty(EEG.etc.amica.num_models)
        if size(EEG.data,2) == size(EEG.etc.amica.v_smooth,2) && size(EEG.data,3) == size(EEG.etc.amica.v_smooth,3) && size(EEG.etc.amica.v_smooth,1) == EEG.etc.amica.num_models

            EEG = eeg_formatamica(EEG);
            %-------------------------------------------

            [EEG com] = eeg_eegrej(EEG,regions);

            %-------------------------------------------

            EEG = eeg_reformatamica(EEG);
            EEG = eeg_checkamica(EEG);
            return;
        else
            disp('AMICA probabilities not compatible with size of data, probabilities cannot be epoched')
            disp('Load AMICA components before extracting epochs')
            disp('Resuming rejection...')
        end
    end

end
% ------------------------------------------------------

if isfield(EEG.event, 'latency'),
   	 tmpevent = EEG.event;
     tmpalllatencies = [ tmpevent.latency ];
else tmpalllatencies = [];
end;

% handle regions from eegplot
% ---------------------------
if size(regions,2) > 2, regions = regions(:, 3:4); end;

[EEG.data EEG.xmax tmpalllatencies boundevents] = eegrej( EEG.data, ...
												  regions, EEG.xmax-EEG.xmin, tmpalllatencies);
if isfield(EEG,'srcpot') && ~isempty(EEG.srcpot) && ~isscalar(EEG.srcpot)
    [EEG.srcpot EEG.xmax garbage boundevents] = eegrej( EEG.srcpot, regions, EEG.xmax-EEG.xmin, []);
end
oldEEGpnts = EEG.pnts;
EEG.pnts   = size(EEG.data,2);
EEG.xmax   = EEG.xmax+EEG.xmin;

% add boundary events
% -------------------
if ~isempty(boundevents) % boundevent latencies will be recomputed in the function below
    [ EEG.event ] = eeg_insertbound(EEG.event, oldEEGpnts, regions);
    EEG = eeg_checkset(EEG, 'eventconsistency');
end;

com = sprintf('%s = eeg_eegrej( %s, %s);', inputname(1), inputname(1), vararg2str({ regions }));
return;

author: PatternRecognition

repo: OpenBMI


text:
function plotindividualw(w,peval,savethis,localize)
% plotindividualw(w,peval,savethis, localize)
if ~exist('localize','var')
    localize = 0;
end

[x_mu, y_mu sig] = localizew(w,peval);
wr=reshape(w, peval.nx, peval.ny, peval.ncomp);

for ii=1:peval.ncomp
    dipshow(wr(:,:,ii));
    if localize
        hold on
        scatter(x_mu(ii)-1, y_mu(ii)-1,'r')
        drawpixels([peval.nx, peval.ny],':r',-0.5)
        hold off
    end
    if savethis
        SaveImageFULL([peval.res_dir '/w' num2str(ii)],'e');
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ l, p, u ] = r8mat_lu ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_LU computes the LU factorization of an R8MAT.
%
%  Discussion:
%
%    The routine is given an M by N matrix A, and produces
%
%      L, an M by M unit lower triangular matrix,
%      U, an M by N upper triangular matrix, and
%      P, an M by M permutation matrix P,
%
%    so that
%
%      A = P' * L * U.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the M by N matrix to be factored.
%
%    Output, real L(M,M), the M by M unit lower triangular factor.
%
%    Output, real P(M,M), the M by M permutation matrix.
%
%    Output, real U(M,N), the M by N upper triangular factor.
%

%  Initialize:
%
%    U:=A
%    L:=Identity
%    P:=Identity
%
  u(1:m,1:n) = a(1:m,1:n);

  l = r8mat_identity ( m );

  p(1:m,1:m) = l(1:m,1:m);
%
%  On step J, find the pivot row, IPIV, and the pivot value PIVOT.
%
  for j = 1 : min ( m - 1, n )

    pivot = 0.0;
    ipiv = 0;

    for i = j : m

      if ( pivot < abs ( u(i,j) ) )
        pivot = abs ( u(i,j) );
        ipiv = i;
      end

    end
%
%  Unless IPIV is zero, swap rows J and IPIV.
%
    if ( ipiv ~= 0 )

      u = r8row_swap ( m, n, u, j, ipiv );

      l = r8row_swap ( m, m, l, j, ipiv );

      p = r8row_swap ( m, m, p, j, ipiv );
%
%  Zero out the entries in column J, from row J+1 to M.
%
      for i = j+1 : m

        if ( u(i,j) ~= 0.0 )

          l(i,j) = u(i,j) / u(j,j);

          u(i,j) = 0.0;

          u(i,j+1:n) = u(i,j+1:n) - l(i,j) * u(j,j+1:n);

        end

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x1,y1,x2,y2] = retrieveCoordinates(index,scale)
% compute coordinates [x1,y1,x2,y2] from an index which represent the a
% window in the scale^4 space of all the windows

image_area = scale*scale;
index1 = mod(index, image_area);
index2 = floor(index/image_area);

x1 = mod(index1, scale) + 1;
y1 = floor(index1/scale) + 1;

x2 = mod(index2, scale) + 1;
y2 = floor(index2/scale) + 1;

return

author: Cloud-CV

repo: object-proposals


text:
%CELLWRITE Write a cell array to a comma separated value file.
%   CELLWRITE(FILENAME, C) writes cell array C into FILENAME as comma
%   separated values.
%
%   NOTE: This function is not completely compatible with CSVWRITE.
%   Offsets are not supported and 0 values are not omitted.
%
%   See also CSVWRITE, CSVREAD, DLMREAD, DLMWRITE, WK1READ, WK1WRITE.
function cellwrite(filename, cellarray)
% The cell array is traversed, the contents of each cell are converted
% to a string, and a CSV file is written using low level fprintf
% statements.
[rows, cols] = size(cellarray);
fid = fopen(filename, 'w');
for i_row = 1:rows
    file_line = '';
    for i_col = 1:cols
        contents = cellarray{i_row,i_col};
        if isnumeric(contents)
            contents = num2str(contents);
        elseif isempty(contents)
            contents = '';
        end
        if i_col<cols
            file_line = [file_line, contents, ','];
        else
            file_line = [file_line, contents];
        end
    end
    count = fprintf(fid,'%s\n',file_line);
end
st = fclose(fid);
if st == -1
    error('Bad file write')
end

author: Sable

repo: mcbench-benchmarks


text:
function gc = getGraphClass(obj)
if isHG2
    gc = lower(obj.Type);
else
    gc = lower(handle(obj).classhandle.name);
end
end

author: plotly


text:
function [ a, rcond, z, info ] = dpoco ( a, lda, n )

%*****************************************************************************80
%
%% DPOCO factors a real symmetric positive definite matrix and estimates its condition.
%
%  Discussion:
%
%    If RCOND is not needed, DPOFA is slightly faster.
%
%    To solve A*X = B, follow DPOCO by DPOSL.
%
%    To compute inverse(A)*C, follow DPOCO by DPOSL.
%
%    To compute determinant(A), follow DPOCO by DPODI.
%
%    To compute inverse(A), follow DPOCO by DPODI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the symmetric matrix to be factored.  Only the
%    diagonal and upper triangle are used.
%
%    Input, integer LDA, the leading dimension of the array A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(LDA,N), an upper triangular matrix R so that A = R'*R
%    where R' is the transpose.  The strict lower triangle is unaltered.
%    If INFO /= 0, the factorization is not complete.
%
%    Output, real RCOND, an estimate of the reciprocal
%    condition of A.  For the system A*X = B, relative perturbations in
%    A and B of size EPSILON may cause relative perturbations in X of
%    size EPSILON/RCOND.  If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate underflows.
%
%    Output, real Z(N), a work vector whose contents are usually
%    unimportant.  If A is close to a singular matrix, then Z is an
%    approximate null vector in the sense that
%      norm(A*Z) = RCOND * norm(A) * norm(Z).
%    If INFO /= 0, Z is unchanged.
%
%    Output, integer INFO, error flag.
%    0, for normal return.
%    K, signals an error condition.  The leading minor of order K is not
%    positive definite.
%

%
%  Find norm of A using only upper half.
%
  for j = 1 : n
    z(j) = dasum ( j, a(1:j,j), 1 );
    for i = 1 : j-1
      z(i) = z(i) + abs ( a(i,j) );
    end
  end

  anorm = max ( z(1:n) );
%
%  Factor.
%
  [ a, info ] = dpofa ( a, lda, n );

  if ( info ~= 0 )
    return
  end
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where R'*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve R' * W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  for k = 1 : n

    if ( z(k) ~= 0.0 )
      ek = - abs ( ek ) * r8_sign ( z(k) );
    end

    if ( a(k,k) < abs ( ek - z(k) ) )
      s = a(k,k) / abs ( ek - z(k) );
      z(1:n) = s * z(1:n);
      ek = s * ek;
    end

    wk = ek - z(k);
    wkm = -ek - z(k);
    s = abs ( wk );
    sm = abs ( wkm );
    wk = wk / a(k,k);
    wkm = wkm / a(k,k);

    if ( k + 1 <= n )

      for j = k+1 : n
        sm = sm + abs ( z(j) + wkm * a(k,j) );
        z(j) = z(j) + wk * a(k,j);
        s = s + abs ( z(j) );
      end

      if ( s < sm )
        t = wkm - wk;
        wk = wkm;
        z(k+1:n) = z(k+1:n) + t * a(k,k+1:n);
      end

    end

    z(k) = wk;

  end

  z(1:n) = z(1:n) / dasum ( n, z(1:n), 1 );
%
%  Solve R * Y = W.
%
  for k = n : -1 : 1

    if ( a(k,k) < abs ( z(k) ) )
      s = a(k,k) / abs ( z(k) );
      z(1:n) = s * z(1:n);
    end

    z(k) = z(k) / a(k,k);
    t = -z(k);
    z(1:k-1) = daxpy ( k-1, t, a(1:k-1,k)', 1, z(1:k-1), 1 );

  end

  z(1:n) = z(1:n) / dasum ( n, z(1:n), 1 );
  ynorm = 1.0;
%
%  Solve R' * V = Y.
%
  for k = 1 : n

    z(k) = z(k) - ddot ( k-1, a(1:k-1,k), 1, z(1:k-1), 1 );

    if ( a(k,k) < abs ( z(k) ) )
      s = a(k,k) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    z(k) = z(k) / a(k,k);

  end

  s = 1.0 / dasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;
%
%  Solve R * Z = V.
%
  for k = n : -1 : 1

    if ( a(k,k) < abs ( z(k) ) )
      s = a(k,k) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    z(k) = z(k) / a(k,k);
    t = -z(k);
    z(1:k-1) = daxpy ( k-1, t, a(1:k-1,k)', 1, z(1:k-1), 1 );

  end
%
%  Make ZNORM = 1.0.
%
  s = 1.0 / dasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end


text:
function x = cvx_c2r( x, dim, cleanup_eps )

%
% Determine expansion dimension
%

sx = size( x );
if nargin < 2,
    dim = [ find( sx > 1 ), 1 ];
    dim = dim( 1 );
elseif ~isnumeric( dim ) || dim <= 0 || dim ~= floor( dim ),
    error( 'Second argument must be a dimension.' );
end
sx = [ sx, ones( 1, dim - length( sx ) ) ];
nd = length( sx );

%
% Perform the sparse case differently
%

if isnumeric( x ) && issparse( x ) && dim <= 2,
    [ rr, cc, vv ] = find( x );
    vr = real( vv );
    vi = imag( vv );
    if nargin > 2,
        ndxs = find( vr & vi );
        if ~isempty( ndxs ),
            temp = abs( vr(ndxs) ./ vi(ndxs) );
            vr(ndxs(temp<=cleanup_eps)) = 0;
            vi(ndxs(temp>=1.0./cleanup_eps)) = 0;
        end
    end
    if dim == 1,
        rr = 2 * rr; rr = [ rr - 1 ; rr ];
        cc = [ cc ; cc ];
        sx( 1 ) = 2 * sx( 1 );
    else
        cc = 2 * cc; cc = [ cc - 1 ; cc ];
        rr = [ rr ; rr ];
        sx( 2 ) = 2 * sx( 2 );
    end
    x = sparse( rr, cc, [ vr ; vi ], sx( 1 ), sx( 2 ) );
    return
end

%
% Permute if necessary
%

perm = [];
if any( sx( 1 : dim - 1 ) ~= 1 ),
    perm = [ dim, 1 : dim - 1, dim + 1 : nd ];
    x = permute( x, perm );
    sx = sx( perm );
    dim = 1;
end

%
% Perform expansion and possibly cleanup
%

x = x( : ).';
sx( dim ) = 2 * sx( dim );
xr = real( x );
xi = imag( x );
if nargin > 2,
    ndxs = find( xr & xi );
    if ~isempty( ndxs ),
        temp = abs( xr(ndxs) ./ xi(ndxs) );
        xr(ndxs(temp<=cleanup_eps)) = 0;
        xi(ndxs(temp>=1.0./cleanup_eps)) = 0;
    end
end
x = reshape( [ xr ; xi ], sx );

%
% Reverse permute if necessary
%

if ~isempty( perm ),
    x = ipermute( x, perm );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function area = annulus_sector_area_2d ( r1, r2, theta1, theta2 )

%*****************************************************************************80
%
%% ANNULUS_SECTOR_AREA_2D computes the area of an annular sector in 2D.
%
%  Discussion:
%
%    An annular sector with center (XC,YC), inner radius R1 and
%    outer radius R2, and angles THETA1, THETA2, is the set of points
%    (X,Y) so that
%
%      R1**2 <= (X-XC)**2 + (Y-YC)**2 <= R2**2
%
%    and
%
%      THETA1 <= THETA ( X - XC, Y - YC ) <= THETA2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, R2, the inner and outer radii.
%
%    Input, real THETA1, THETA2, the angles.
%
%    Output, real AREA, the area.
%
  area = 0.5 * ( theta2 - theta1 ) * ( r2 + r1 ) * ( r2 - r1 );

  return
end


text:
function [bndinfo2, err] = updateBoundaryInfo(bndinfo, result, im)

if isstruct(result) && isfield(result, 'regions')
    regions = result.regions;

    empty = false(numel(regions), 1);
    splab = zeros(bndinfo.nseg, 1);
    for k = 1:numel(regions)
        splab(regions{k}) = k;
        empty(k) = isempty(regions{k});
    end
    regions(empty) = [];
else
    splab = result;
    regions = cell(max(splab), 1);
    empty = false(numel(regions), 1);
    for k = 1:numel(regions)
        regions{k} = find(splab==k);
        empty(k) = isempty(regions{k});
    end
    regions(empty) = [];
end

wseg = splab(bndinfo.wseg);

% XXX only reading image until seg2fragments gets fixed
%im = imread(['/usr1/projects/dhoiem/iccv07/iccvGroundTruth/images/' bndinfo.imname]);
%im = im2double(im);

[edges, juncts, neighbors, wseg] = seg2fragments(wseg, im, 1);
bndinfo2 = processBoundaryInfo(wseg, edges, neighbors);

if isfield(bndinfo, 'imname')
    bndinfo2.imname = bndinfo.imname;
end

%bndinfo2.imname = bndinfo.imname;

stats = regionprops(bndinfo.wseg, 'Area');
area = cat(1, stats(:).Area);
bndinfo2.spArea = area;

if isfield(bndinfo, 'type')
    bndinfo2.type = bndinfo.type;
    bndinfo2.names = bndinfo.names;
    [bndinfo2.labels, err] = iccvTransferLabels(bndinfo.labels, regions, area);
    bndinfo2 = processGtBoundaryLabels(bndinfo2);

else
    err = nan;
end

bndinfo2 = orderfields(bndinfo2);


text:
function [u,v,t]=rotro2pl(r)
%ROTRO2PL find the plane and rotation angle of a rotation matrix [u,v,t]=r
% Inputs:
%
%     R(n,n)   Rotation matrix
%
% Outputs:
%
%     U(n,1) and V(n,1) are orthonormal vectors defining a plane in n-dimensional space
%     T is the rotation angle in radians from U towards V with 0<=T<=pi. If T
%       is omitted it U and V will be separated by T instead of being orthogonal

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

n=size(r,1);
[q,e]=schur(r);
[m,i]=max(abs(e(2:n+1:n^2)));
z=e(i+1,i)<0; % =1 if negative
uv=q(:,i+z:1-2*z:i+1-z);
u=uv(:,1);
% the following code selects unique values of u and v
% v=uv(:,2);
% f=u.*v;         % maximize inner product of u.^2 and v.^2
% g=(v+u).*(v-u);
% t=atan2(sum(f.*g),sum(g.^2/4-f.^2))/4;
% c=cos(t);
% s=sin(t);
% uv=uv*[c s; -s c];
% a=sum(uv)<0;
% uv=uv*[1-a(1)-a(2) a(2)-a(1); a(1)-a(2) 1-a(1)-a(2)];
% u=uv(:,1);
if nargout>2
    v=uv(:,2);
    t=atan2(abs(e(i+1,i)),e(i,i));
else
    [s,c]=atan2sc(abs(e(i+1,i)),e(i,i));
    v=uv*[c;s];
end

author: decouples


text:
function r = acot(a)
%ACOT         Taylor inverse cotangent  acot(a)
%
%Thanks to George Corliss for providing the Taylor expansion
%

% written  06/03/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                   % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = getappdata(0,'INTLAB_TAYLOR_ORDER');

  r = a;
  ct = a.t;
  N = size(a.t,2);
  r.t(1,:) = acot(a.t(1,:));
  ct1 = 1+a.t(1,:).^2;     % 1+a^2
  r.t(2,:) = -a.t(2,:) ./ ct1 ;
  for j=2:K
    ct(j,:) = sum( a.t(1:j,:).*a.t(j:-1:1,:) , 1 );
    r.t(j+1,:) = ( j*a.t(j+1,:) + sum( repmat((1:j-1)',1,N).*r.t(2:j,:).*ct(j:-1:2,:) , 1 ) ) ./ ( (-j)*ct1 );
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function g = gammaFun(j, k)
%GAMMAFUN   Get a function handle to a gamma function.
%   G = GAMMAFUN(J, K) returns a function handle to the gamma function (J, K).
%
% See also EXPINTEG/GAMMAEVAL.

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

% Trivial case:
if ( j == 0 )
    g = @(z) (exp(k*z) - 1)./z;

% Compute them recursively using the recurrence formula:
else
    g = @(z) 0*z;
    for m = 1:j
        gm = expinteg.gammaFun(j-m, k);
        g = @(z) g(z) + (-1)^(m-1)/m*gm(z);
    end
    if ( j <= k )
        g = @(z) (g(z) - nchoosek(k, j))./z;
    else
        g = @(z) g(z)./z;
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function cMap=kwJet(varargin)

switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

cMap=[0 0 0; gjet(4); 1 1 1];

[cMap]=resampleColormap(cMap,n);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function anorm = update_gbound(anorm,alpha,beta,j)
%UPDATE_GBOUND   Update Gerscgorin estimate of 2-norm
%  ANORM = UPDATE_GBOUND(ANORM,ALPHA,BETA,J) updates the Gersgorin bound
%  for the tridiagonal in the Lanczos process after the J'th step.
%  Applies Gerscgorins circles to T_K'*T_k instead of T_k itself
%  since this gives a tighter bound.

if j==1 % Apply Gerscgorin circles to T_k'*T_k to estimate || A ||_2
  i=j;
  % scale to avoid overflow
  scale = max(abs(alpha(i)),abs(beta(i+1)));
  alpha(i) = alpha(i)/scale;
  beta(i) = beta(i)/scale;
  anorm = 1.01*scale*sqrt(alpha(i)^2+beta(i+1)^2 + abs(alpha(i)*beta(i+1)));
elseif j==2
  i=1;
  % scale to avoid overflow
  scale = max(max(abs(alpha(1:2)),max(abs(beta(2:3)))));
  alpha(1:2) = alpha(1:2)/scale;
  beta(2:3) = beta(2:3)/scale;

  anorm = max(anorm, scale*sqrt(alpha(i)^2+beta(i+1)^2 + ...
      abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
      abs(beta(i+1)*beta(i+2))));
  i=2;
  anorm = max(anorm,scale*sqrt(abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1))) );
elseif j==3
  % scale to avoid overflow
  scale = max(max(abs(alpha(1:3)),max(abs(beta(2:4)))));
  alpha(1:3) = alpha(1:3)/scale;
  beta(2:4) = beta(2:4)/scale;
  i=2;
  anorm = max(anorm,scale*sqrt(abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
      abs(beta(i+1)*beta(i+2))) );
  i=3;
  anorm = max(anorm,scale*sqrt(abs(beta(i)*beta(i-1)) + ...
      abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1))) );
else
  % scale to avoid overflow
  %  scale = max(max(abs(alpha(j-2:j)),max(abs(beta(j-2:j+1)))));
  %  alpha(j-2:j) = alpha(j-2:j)/scale;
  %  beta(j-2:j+1) = beta(j-2:j+1)/scale;

  % Avoid scaling, which is slow. At j>3 the estimate is usually quite good
  % so just make sure that anorm is not made infinite by overflow.
  i = j-1;
  anorm1 = sqrt(abs(beta(i)*beta(i-1)) + ...
      abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
      abs(beta(i+1)*beta(i+2)));
  if isfinite(anorm1)
    anorm = max(anorm,anorm1);
  end
  i = j;
  anorm1 = sqrt(abs(beta(i)*beta(i-1)) + ...
      abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
      beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
      abs(alpha(i)*beta(i+1)));
  if isfinite(anorm1)
    anorm = max(anorm,anorm1);
  end
end

author: andrewssobral

repo: lrslibrary


text:
function x = cs_reach(L,b)                                                  %#ok
%CS_REACH non-recursive reach (interface to CSparse cs_reach)
% find nonzero pattern of x=L\sparse(b).  L must be sparse, real, and lower
% triangular.  b must be a real sparse vector.
%
% Example:
%   x = cs_reach(L,b)
% See also: cs_demo

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

error ('cs_reach mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function mt = SupIntensity_HP(t, History, para, options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute the super bound of intensity function of Hawkes processes
%
% Parameters of Hawkes processes
% para.mu: base exogenous intensity
% para.A: coefficients of impact function
% para.kernel: 'exp', 'gauss'
% para.w: bandwith of kernel
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(History)
    mt = sum(para.mu);
else
    Time = History(1, :);
    index = Time<=t;
    Time = Time(index);
    Event = History(2, index);

    MT = sum(para.mu)*ones(1, options.M);
    for m=1:options.M
        t_current = t+(m-1)*options.tstep/options.M;

        basis = Kernel(t_current-Time(:), para);
        A = para.A(Event, :, :);

        for c = 1:size(para.A, 3);
            MT(m) = MT(m) + sum(sum(basis.*A(:,:,c)));
        end
    end
    mt = max(MT);
end

mt = mt.*(mt>0);
end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
% vgg_mrdivs  Solves equation system Y*diag(s) = A*X with unkowns A, s.
%
% A = vgg_mrdivs(X,Y) solves (overdetermined) equation system Y*diag(s) = A*X
% by linear method (DLT algorithm).
% Parameters:
%   X ... double (N,K)
%   Y ... double (M,K)
%   A ... double (M,N)
%   s ... double (1,K)
%
% Preconditioning of the points not included in the function. Use vgg_conditioner_*.
%
% Typical usage:
%   1. Estimating an image homography from K pairs of corresponding points.
%      If 3-by-K matrices x and y are the points in homogeneous coordinates, the 3-by-3 homography
%      matrix is obtained as H = vgg_mrdivs(x,y).
%
%   2. Estimating 3-by-4 camera projection matrix P from corresponding pairs of image and scene points.
%      For image points x (3xK matrix) and scene points X (4xK matrix) do P = vgg_mrdivs(X,x).

% (c) werner@robots.ox.ac.uk

% Algorithm:
%
% For n-th point pair X(:,n) and Y(:,n) we have
%  s*X(:,n) = A*Y(:,n)
% We eliminate s what results in MY*(MY-1)/2 linear homogenous equations
% for elements of A. We solve this system by svd or eig.

function A = vgg_mrdivs(X,Y)

[MX,N] = size(X);
[MY,NY] = size(Y);
if N ~= NY, error('Matrices A, B must have equal number of columns.'); end

 % construct the measurement matrix
W = zeros(MX*MY,MY*(MY-1)/2*N);
k = 1;
for i = 1:MY
  for j = 1:i-1
    W([[1:MX]+MX*(j-1) [1:MX]+MX*(i-1)],[1:N]+N*(k-1)) =  [(ones(MX,1)*Y(i,:)).*X; -(ones(MX,1)*Y(j,:)).*X];
    k = k+1;
  end
end

% solve the system || A'(:)' * W || ---> min
[dummy,s,A] = svd(W',0);
A = reshape(A(:,end),MX,MY)';

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function [uniqName,scanUID] = genScanUniqName(planC, scanNum)
% previously from register_scans.m
if nargin < 2
    scanNum = 1;
end

randPart = num2str(floor(rand*1000));

if iscell(planC)
    indexS = planC{end};
    scanUID = planC{indexS.scan}(scanNum).scanUID;
    uniqName = [scanUID randPart];
else
    [~,f,ext] = fileparts(planC);
    if strcmp(ext,'.gz')
        [~,f,~] = fileparts(f);
    end
    orgRoot = '1.3.6.1.4.1.9590.100.1.2';
    scanUID = char(javaMethod('createUID','org.dcm4che3.util.UIDUtils',...
        orgRoot));
    %scanUID = dicomuid;
    uniqName = [f '_' scanUID];
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function y = han(x,frac)
    % apply a hamming window to signal x of length x*(1+frac).
    % equivalent to zero-padding, windowing and then truncating.

    if nargin < 2 || isempty(frac)
        frac = 0.5;
    end

    if isinf(frac);
        y = x;
    else
        sx = size(x);
        h1 = hamming(sx(1) + 2*round(sx(1)*frac));
        h2 = hamming(sx(2) + 2*round(sx(2)*frac));
        %y = x.*(h1(round(sx(1)*frac)+(1:sx(1)))*h2(round(sx(2)*frac)+(1:sx(2)))');
        y = bsxfun(@times, bsxfun(@times,h2(round(sx(2)*frac)+(1:sx(2)))',x), h1(round(sx(1)*frac)+(1:sx(1))));
        if length(sx) == 3
            temp = hamming(sx(3) + 2*round(sx(3)*frac));
            h3(1,1,1:sx(3)) = temp(round(sx(3)*frac)+(1:sx(3)));
            y = bsxfun(@times, y, h3);
        end
    end

author: flatironinstitute

repo: NoRMCorre


text:
function[num,a,b]=blocknum(x,delta)
%BLOCKNUM  Numbers the contiguous blocks of an array.
%
%   Suppose X is a column vector which contains blocks of identical
%   values, e.g. X=[0 0 0 1 1 3 2 2 2]';
%
%   N=BLOCKNUM(X) gives each contiguous block of identical values a
%   unique number, in order beginning with 1.  Each elements of N
%   specifies the number of the block to which the corresponding
%   element of X belongs.
%
%   In the above example, N=[1 1 1 2 2 3 4 4 4]';
%
%   [N,A,B]=BLOCKNUM(X) also returns arrays A and B which are indices
%   into the first and last point, respectively, of each block.  In
%   the above example, A=[1 4 6 7]' and B=[3 5 6 9]';
%
%   [...]=BLOCKNUM(X,D) defines the junction between two blocks as
%   locations where ABS(DIFF(X))>D.  Thus D=1 is a 'rate of change'
%   definition, and X=[1 2 3 5 6 10 16 17 18]'; will yield the same
%   result for N as in the previous example.
%
%   See also BLOCKLEN.
%
%   Usage: num=blocknum(x);
%          [num,a,b]=blocknum(x);
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2014 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(x,'--t')
  blocknum_test;return;
end

if nargin==1
  delta=0;
end

dx=diff(x);
index=find(abs(dx)>delta);
num=zeros(size(x));
if ~isempty(index)
  num(index+1)=1;
end
num(1)=1;
num=cumsum(num);

if nargout>=2
   a=[1;find(diff(num)~=0)+1];
end
if nargout>=3
  b=[find(diff(num)~=0);length(num)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[]=blocknum_test
x= [0 0 0 1 1 3 2 2 2]';
y= [1 1 1 2 2 3 4 4 4]';

reporttest('BLOCKNUM with D==0',all(y==blocknum(x)))

x=[1 2 3 5 6 10 16 17 18]';
reporttest('BLOCKNUM with D==1',all(y==blocknum(x,1)))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [ a x e ] = sdr_ialm(A, B, y, beta, gamma, epsilon, maxIter)

%-------------------------------------------------------------------------------------
% Dec 2014
%
% by Jian Lai, jlai1@ntu.edu.sg
% solve the equation (10) Sparse and Dense Hybrid Representation in our PAMI paper
% min ||a||_1 + \gamma ||x||_2^2 + \beta ||e||_1 s.t. y = Aa + Bx + e.

% input
% A: class-specific dictionary, each column is a atom
% B: nonclass-specific dictionary, each column is a atom
% beta: the balacing parameter of e
% gamma: the balancing parameter of x
% epsilon: the tolerance of the convergent checking
% maxIter: max iteration limit

% output
% a: the coefficients of A
% x: the coefficinets of B
% e: the error
%-------------------------------------------------------------------------------------

% check input
if nargin < 5
    error('Too few arguments') ;
end

if nargin < 6
    epsilon = 1e-6;
end

if nargin < 7
    maxIter = 1000;
end

% initialization
a=zeros(size(A,2),1);
x=zeros(size(B,2),1);
e=zeros(size(y));
phi=zeros(size(y));
[m, n]=size(A);
mu=1;
mu_bar = 1e+6;
rho=1.5;
iter = 0;
converged = false;
BTB=B'*B;

while ~converged

    iter = iter + 1;

	% optimize e
    Resi_e = y - A*a - B*x + phi/mu;
    e = max(Resi_e - beta/mu, 0)+min(Resi_e + beta/mu, 0);

    %optimize x
    Resi_x = B'*( phi + mu * (y-A*a-e) );
    x=(2*gamma*eye(n)+mu*BTB)\Resi_x;

    %optmize a
    Resi_a=y - B * x - e + phi/mu;
    a=SolveHomotopy(A, Resi_a, 'lambda', 1/mu, 'tolerance', 1e-5, 'stoppingcriterion', 3);

    %update \phi and \mu
    Z = y - A*a - B*x - e ;
    phi = phi + mu*Z;
    mu = min(mu*rho, mu_bar);

    %check convergence
    if norm(Z) < epsilon
        converged = true;
    end

    if ~converged && iter >= maxIter
        disp('Maximum iterations reached') ;
        converged = 1 ;
    end

end

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function esportspect_Callback(hObject, eventdata, handles)
% Plotting Defferntial Expression

[~,detectionName] = fileparts(handles.current_detection_file);
spectname = [detectionName ' Call ' num2str(handles.data.currentcall) '.png'];
[FileName,PathName] = uiputfile(spectname,'Save Spectrogram');

% Cancel if cancelled
if isnumeric(FileName)
    return
end

% Get the spectrogram from the display
I = get(handles.spect,'CData');
Ydata = get(handles.spect,'Ydata');
Ylim = get(handles.focusWindow,'Ylim');

% Set limits
Ymax = find(Ydata>=Ylim(2),1);
Ymin = find(Ydata>=Ylim(1),1);
I = flipud(I(Ymin:Ymax,:));

Clim = get(handles.focusWindow,'Clim');
I = mat2gray(I,Clim);
fullFileName = fullfile(PathName, FileName); % Add Figure Path
cmap = handles.data.cmap;
I2 = gray2ind(I,256);
RGB3 = ind2rgb(I2, inferno(255));
imwrite(RGB3,fullFileName,'png','BitDepth',8); % Re-change it to colored one

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function EEG = hlp_importConnData(varargin)
% Convert/Import connectivity matrices and (optionally) source locations
% into SIFT data structures and store in an EEGLAB data structure
% This can be used to visualize, in SIFT, connectivity data computed
% using another toolbox, such as SPM, GCCA, or Fieldtrip.

% Author: Tim Mullen, 2012, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

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

g= arg_define(varargin,...
    arg({'connDataStruct','ConnectivityDataStructure'},[],mandatory,'Connectivity data structure. This is a structure where each field ConnDataStruct.(<MethodName>) = [vars_to x vars_from x freqs x times] contains a 2D, 3D, or 4D connectivity matrix for a given MethodName. All matrices must be of equal size and first two dimensions must be equal size (num_vars_to==num_vars_from). If this describes an asymmetric (directed) graph, edges in the graph point from columns (vars_from) to rows (vars_to). If there is no frequency information, the 3rd dimension is a singleton (matrix is still 4-D), if there is no time information then the matrix is 3-dimensional (matlab always squeezes out trailing singleton dimensions). If there is neither frequency nor time information, matrix is 2D.'), ...
    arg({'winCenterTimes','WinCenterTimes'},[],[],'Time-stamps (sec). This is a vector of time-stamps for corresponding time-varying connectivity (winCenterTimes=[], if no time-varying connectivity)','shape','row'), ...
    arg({'erWinCenterTimes','EventRelatedWinCenterTimes'},[],[],'Event-related time-stamps (sec). This is a vector of "event-related" time-stamps for event-related analysis. Time t=0 is an event to which trials have been locked, t>0 is post-event and t<0 is pre-event. If empty (default), value is same as winCenterTimes.','shape','row'), ...
    arg({'freqs','Frequencies'},0,[],'Frequencies (Hz). This is a vector of frequencies for corresponding frequency-varying connectivity (freqs=0, if no frequency information)','shape','row'), ...
    arg({'winlen','SlidingWindowLength'},0,[],'Sliding Window Length (sec). This is the length of the sliding window (in sec) if you used a sliding-window analysis to get time-varying connectivity. Default is 0 (no sliding window used)'), ...
    arg({'variableNames','VariableNames','curComponentNames'},{},[],'Variable IDs. This is a cell array of ''user-friendly'' names for each variable -- e.g. anatomical names for each source, or channel names for each channel','type','cellstr','shape','row'), ...
    arg('MNI_XYZ',[],[],'MNI Coords (X,Y,Z) for each variable. An EEGLAB dipfit structure can also be provided. Otherwise, this is an optional matrix of size [num_vars x 3] where each column contains, respectively, X,Y,Z MNI coordinates for the variables (e.g. channels or sources). Subject-specific coordinates can also be provided, but then a subject-specific head model will need to be provided when rendering a brainmovie. This goes into the .posxyz field of a new EEGLAB dipfit structure','shape','matrix'), ...
    arg('EEG',[],[],'EEGLAB data structure. Optional EEGLAB data struct. SIFT data objects will be stored in here. If ommitted, an empty EEGLAB dataset will be created.') ...
    );

% Input checks and conversions
% -------------------------------------------------------------------------
fnames = fieldnames(g.connDataStruct);
if isempty(g.connDataStruct)
    error('Input argument ''ConnectivityDataStructure'' cannot be empty');
end

% get connectivity matrix dimensions
% and check that all connectivity matrices are of equal size
[numVars numVars numFreqs numTimes] = size(g.connDataStruct.(fnames{1}));
for k=2:length(fnames)
    if        ndims(g.connDataStruct.(fnames{k})) ~= ndims(g.connDataStruct.(fnames{1})) ...
       || any( size(g.connDataStruct.(fnames{k})) ~=  size(g.connDataStruct.(fnames{1})) )
        error('All connectivity matrices in ConnectivityDataStructure must be of equal size');
    end
end

% check time-stamp vectors
if isempty(g.winCenterTimes)
    g.winCenterTimes = 0;
end
if isempty(g.erWinCenterTimes)
    g.erWinCenterTimes = g.winCenterTimes;
end
if length(g.winCenterTimes)~=numTimes
    error('The number of time-stamps (length of ''WinCenterTimes'') does not agree with the corresponding dimension (dim 4) of the connectivity matrices');
end
if length(g.erWinCenterTimes)~=length(g.winCenterTimes)
    error('The number of event-related time-stamps (length of ''EventRelatedWinCenterTimes'') must equal the length of WinCenterTimes.');
end

% check frequency vector
if length(g.freqs)~=numFreqs
    error('The number of frequencies (length of ''Frequencies'') does not agree with the corresponding dimension (dim 3) of the connectivity matrices');
end

% check variable names
variableInds = 1:numVars;
if isempty(g.variableNames)
    g.variableNames = cellstr(num2str(variableInds'))';
elseif length(g.variableNames) ~= length(variableInds)
    error('The length of VariableNames must equal the number of variables');
end

% check MNI coords matrix
if ~isempty(g.MNI_XYZ) && ~isstruct(g.MNI_XYZ)
    if size(g.MNI_XYZ,1)~=numVars || size(g.MNI_XYZ,2) ~=3
        error('''MNI_XYZ'' must be either a dipfit structure or a matrix of size [num_vars x 3]');
    end
end

% Make Conn object
% -------------------------------------------------------------------------
Conn = struct(...
    'winCenterTimes'        , g.winCenterTimes     , ... % vector of time-stamps for corresponding time-varying connectivity (g.winCenterTimes=0, if no time-varying connectivity)
    'erWinCenterTimes'      , g.erWinCenterTimes   , ... % vector of "event-related" time-stamps for event-related analysis. Time t=0 is an event to which trials have been locked, t>0 is post-event and t<0 is pre-event. Otherwise, can be the same as g.winCenterTimes
    'freqs'                 , g.freqs                ... % vector of frequencies (Hz) for corresponding frequency-varying connectivity (g.freqs=0, if no frequency information)
    );
Conn = catstruct(Conn,g.connDataStruct);

% Make (abbreviated) model object
% -------------------------------------------------------------------------
MODEL = struct(...
    'winlen', 0 ...     % this is the only field you need from MODEL. It contains the length of the sliding window (in sec) if you used a sliding-window analysis to get time-varying connectivity. Otherwise, just set to 0
    );

% Make (abbreviated) SIFT data structure.
% -------------------------------------------------------------------------
CAT = hlp_sift_emptyset(...
    'curComps'          ,   variableInds    , ... % numbers for variables (e.g. channels or sources) here
    'curComponentNames' ,   g.variableNames   , ... % 'user-friendly' name for each variable -- e.g. anatomical names
    'nbchan'            ,   numVars         , ... % number of variables (e.g. numVars above)
    'Conn'              ,   Conn            , ... % Connectivity structure (see above)
    'MODEL'             ,   MODEL             ... % Model structure (see above)
    );

% Make EEGLAB data structure
% -------------------------------------------------------------------------
if isempty(g.EEG)
    % if g.EEG dataset not provided, create an empty one
    g.EEG = eeg_emptyset;
end
g.EEG.CAT = CAT; % add the the SIFT ('CAT') structure

% Make dipfit structure with source/channel coordinates
% -------------------------------------------------------------------------

% These coordinates are typically (although not necessarily) in MNI space.
% If you are working with channels and have a locations file, you can also
% import the channel locations using EEGLAB routines (c.f. pop_chanedit or
% pop_readlocs). You can then convert the channel locations to a dipfit
% structure using SIFT routine hlp_chanlocs2dipfit(). This will allow you
% to render a 'channel-level' BrainMovie
%
% Below is an example of a super-basic dipfit structure which contains MNI
% source coordinates [X Y Z] copied from X, Y, Z coordinate vectors
%
% See http://sccn.ucsd.edu/wiki/A08:_DIPFIT#DIPFIT_structure_and_functions
% for details on the fields in dipfit structure.
%
% NOTE: SIFT by default uses the freesurfer segmentations of the standard
%       (Colins) MNI brain for surface visualization. If you use your own
%       head model, you will need to supply this surface file via the
%       'Layers' option of vis_causalBrainMovie3D().
%       Also, if using default MNI brain and the network looks rotated
%       90 deg, try changing dipfit.coordformat to 'mni'

if ~isempty(g.MNI_XYZ)

    if isstruct(g.MNI_XYZ)
        % just copy the provided dipfit struct
        g.EEG.dipfit = g.MNI_XYZ;
    else
        % build a new dipfit struct
        dipfit.hdmfile  = '';
        dipfit.mrifile  = fullfile(fileparts(utl_whichfile('eeglab')),'plugins','dipfit2.2','standard_BEM','standard_mri.mat');
        dipfit.chanfile = '';
        dipfit.chansel  = 1:numVars;
        dipfit.coordformat = 'MNI';
        dipfit.coord_transform = [0 0 0 0 0 0 1 1 1];

        for i=1:numVars
            dipfit.model(i).posxyz      = g.MNI_XYZ(i,:);
            dipfit.model(i).momxyz      = [0 0 0];
            dipfit.model(i).rv          = 0;
            dipfit.model(i).select      = 1;
            dipfit.model(i).diffmap     = [];
            dipfit.model(i).sourcepot   = [];
            dipfit.model(i).datapot     = [];
        end
        g.EEG.dipfit = dipfit;
    end
end

EEG = g.EEG;

repo: REST


function cleanFiles(cleanBefore)
% clean output files in ./tex using make
%FIXME: this file appears to be unused (but it is useful)
%FIXME: adapt this file to take the output directory into account
    if cleanBefore && exist(fullfile('tex','Makefile'),'file')
        fprintf(1, 'Cleaning output files...\n');
        cwd = pwd;
        try
            cd('tex');
            [exitCode, output] = system('make distclean');
            fprintf(1,'%s\n', output);
            assert(exitCode==0, 'Exit code 0 means correct execution');
        catch
            % This might happen when make is not present
            fprintf(2, '\tNot completed succesfully\n\n');
        end
        cd(cwd);
    end
end


text:
function c=batmask()
%BATMASK  Load a Gabor multiplier symbol for the 'bat' test signal
%   Usage:  c=batmask;
%
%   `batmask` loads a Gabor multiplier with a 0/1 symbol that masks out
%   the main contents of the 'bat' signal. The symbol fits a Gabor
%   multiplier with lattice given by $a=10$ and $M=40$.
%
%   The mask was created manually using a image processing program. The
%   mask is symmetric, such that the result will be real valued if the
%   multiplier is applied to a real valued signal using a real valued
%   window.
%
%   See also:  bat

%   AUTHOR : Peter L. Søndergaard
%   TESTING: TEST_BATMASK
%   REFERENCE: OK

if nargin>0
  error('This function does not take input arguments.')
end;

f=mfilename('fullpath');

c=load('-ascii',[f,'.asc']);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function  tau = dtauint(z, pp, om, omb, omt, h)
% calculates the integrant for integration of the optical length dtau over the redshift z
% dtau in Mpc^-1 is extracted from interpolation in a table which has been splined.
% pp is the spline of the interpolation for tau values as f of the redshift in this table.
% futher arguments in dtau_tabel are :matter fraction (om), baryon fraction (omb and rel Hubble constant (h)
% ha is the full Hubble factor = sqrt(rho/rho_crit)
%
% D Vangheluwe 31 mrt 2005
% modified for curved space with extra argument omt

global GL_cmb_ka1
ka1 = GL_cmb_ka1;
ha = (1e5 * h) * sqrt(om * (1 + z) .^3 + (ka1/h^2) * (1 + z) .^4 + (omt - om - ka1/h^2) + (1 - omt) * (1 + z) .^2);
tau = ppval(pp, z) .* (1 + z) .^ 2 ./ha;


text:
close all; clear all;

%% Setting of the problem
global s
s = 0.3;
maxIt = 4;
pde = fonedata;
pde.s = s;
pde.L = 1;
option.solver = 'direct';
% option.solver = 'amg';
option.gNquadorder = 4;
[node,elem] = squaremesh([-1,1,-1,1],2);
% [node,elem] = delmesh(node,elem,'x>0 & y<0');

%% Finite element approximation
N = zeros(maxIt,1);
energy = zeros(maxIt,1);
for k = 1:maxIt
    [node,elem] = uniformrefine(node,elem);
    tic;
    [uh,eqn] = fracLapP1P1(node,elem,pde,option);
    Nv = size(node,1);
    figure(1); showresult(node,elem,uh(1:Nv));
    pause(0.1)
    N(k) = length(uh);
    energy(k) = (uh'*eqn.A*uh)/2 - sum(eqn.b(1:Nv).*uh(1:Nv));
    fprintf('#Dof %d   ',N(k));
    toc;
end
[node,elem] = uniformrefine(node,elem);
tic;
[uh,eqn,info] = fracLapP1P1(node,elem,pde,option);
toc;
display(length(uh));
fineEnergy = (uh'*eqn.A*uh)/2 - sum(eqn.b.*uh);

%% Plot convergence rates
energyError = sqrt(2*(energy(1:k)-fineEnergy));
figure(2)
showrate(N(1:k),energyError(1:k),1,'r-*','Energy Error');

%% Display error
colname = {' #Dof   sqrt(2*(E(uh)-E(u)))'};
disptable(colname,N,[],energyError,'%0.5e')

repo: ifem


text:
function [tri, x] = delaunay_sphere(N, xi)
% DELAUNAY_SPHERE  Triangulation of a sphere or a surface topologically
% equivalent to a sphere so that it can be used as a mesh
%
% [TRI, X] = delaunay_sphere(N)
%
%   Number of samples of the sphere (in latitude and longitude).
%
%   TRI is a 3-column matrix with a triangulation of the sphere. Each
%   element in TRI is an index to a row in X. Each row represents the three
%   vertices of a triangle on the sphere.
%
%   X is a 3-column matrix with the Euclidean coordinates of the
%   triangulation vertices.
%
%   The surface can be plotted running
%
%     trisurf(tri, x(:, 1), x(:, 2), x(:, 3));
%
% ... = delaunay_sphere(N, XI)
%
%   XI is an (N+1,N+1,3)-array with the coordinates of the vertices of a
%   surface topologically equivalent to a sphere. By default, XI is the
%   concatenation of the output of sphere(N), and the surface is a unit
%   sphere.
%
%   The surface can be plotted running
%
%     surf(xi(:, :, 1), xi(:, :, 2), xi(:, :, 3))
%
% See also: scimat_tri_to_raster.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 2);
nargoutchk(0, 2);

if (nargin > 1 ...
        && (size(xi, 1) ~= N+1 || size(xi, 2) ~= N+1 || size(xi, 3) ~= 3))
    error('XI must be a (N+1, N+1, 3) array')
end

% sphere surface
[x, y, z] = sphere(N);

% remove repeated points
[x, idx] = unique([x(:) y(:) z(:)], 'rows');

% triangulate the sphere and extract the surface
tri = DelaunayTri(x);
tri = freeBoundary(tri);

% if an input point configuration was provided, we use it to replace the
% unit sphere points coordinates
if (nargin > 1)
    xi = reshape(xi, (N+1)*(N+1), 3);
    x = xi(idx, :);
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text: %% SO3FunHandle

author: mtex-toolbox

repo: mtex

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