text:
function model = mapMetsToAltModel(refModel,mapModel,metFields,mapMethod)
%mapMetsToAltModel  Map metabolites from one model to another.
%
% USAGE:
%
%   model = mapMetsToAltModel(refModel,mapModel,metFields,mapMethod);
%
%
% INPUTS:
%
%   refModel   Model structure containing metabolites that are to be mapped
%              to mapModel metabolite identifiers.
%
%   mapModel   Model structure containing metabolites to which refModel
%              will be mapped.
%
%   metFields  One or more metabolite-related fields that will be compared
%              between the two models to map metabolites from one model to
%              the other. If "metNames" is included in metFields,
%              metabolite names will be processed in three steps:
%               1) Map refModel metNames to mapModel metNames and
%                  metNamesAlt
%               2) Map refModel metNames and metNamesAlt to mapModel
%                  metNames and metNamesAlt
%               3) Make specific modifications to met names based on known
%                  differences (unique to situations where Recon3D and HMR
%                  are the input models)
%
%   mapMethod  (Optional) Specify the method used to handle metabolites in
%              refModel that map to multiple IDs in mapModel across
%              multiple metFields.
%                 'all'  (Default) For each metabolite, met IDs will be
%                        returned for all mets in mapModel that were mapped
%                        along any of the specified metFields.
%               'order'  The order in which metFields are listed will
%                        determine their priority in mapping met IDs (where
%                        first is most important/confident, and last is the
%                        least). Once a met is mapped via one of the fields
%                        listed in metFields, there will be no further
%                        attempts to map that metabolite via all subsequent
%                        fields in metFields.
%               'score'  Each met will be mapped along all fields listed in
%                        metFields. The mapped met ID(s) with the highest
%                        score for each metabolite will be kept, whereas
%                        all other mapped IDs with lower scores will be
%                        removed as potential matches. The score is
%                        calculated based on the number of fields in
%                        metFields by which the refModel met is mapped to
%                        the mapModel met. If metFields includes a second
%                        column of field weights, then the score will be
%                        weighted by these values, where mets mapped along
%                        fields with greater weights will receive higher
%                        scores than those mapped along fields with lower
%                        weights.
%
% OUTPUTS:
%
%   model     Model structure with added field "metAltModelID", which
%             contains the mapModel metIDs that were matched to the
%             refModel metabolites.
%

% handle input arguments
if nargin < 4
    mapMethod = 'all';
end

if ischar(metFields)
    % convert char to cell
    metFields = {metFields};
    metFieldWeights = 1;  % weight is irrelevant if only one metField is provided
elseif isvector(metFields)
    % ensure that metFields is a column vector (probably unnecessary)
    metFields = metFields(:);
    metFieldWeights = ones(size(metFields));  % weight all metFields equally
else
    % if metFields contains a second column, it is treated as the vector of
    % weights corresponding to the list of met fields, and is extracted
    metFieldWeights = cell2mat(metFields(:,2));
    metFields = metFields(:,1);
end

if strcmpi(mapMethod,'order')
    % if mapMethod is "order", then weight fields based on their order of appearance in metFields
    metFieldWeights = (length(metFields):-1:1)';
end

% verify that all metFields are present in both models
if any(~isfield(refModel,metFields)) || any(~isfield(mapModel,metFields))
    error('One or more of the specified metFields is not present in one or both models.');
end

% remove metNamesAlt if present in metFields - this field is automatically
% included when metNames is specified
if ismember(metFields,'metNamesAlt')
    fprintf('Note: the "metNamesAlt" entry is ignored, because it is included in the "metNames" processing step.\n');
    metFields(ismember(metFields,'metNamesAlt')) = [];
    if ~ismember('metNames',metFields)
        % warn user if they include metNamesAlt, but not metNames, in metFields
        fprintf('WARNING: comparison by "metNames" has not been specified - mets will not be compared by names.\n');
    end
end

% Add empty metNamesAlt field to either of the models if they are missing
% the field. This just makes things easier later on.
if ~isfield(refModel,'metNamesAlt')
    refModel.metNamesAlt = repmat({''},size(refModel.mets));
end
if ~isfield(mapModel,'metNamesAlt')
    mapModel.metNamesAlt = repmat({''},size(mapModel.mets));
end

% initialize output model
model = refModel;

% strip compartments from mapModel met IDs to obtain compartment-free
% met-to-name (or met-to-ID) pairs
if strcmp(mapModel.mets{1}(end),']')
    % compartment information is formatted in brackets at end of ID
    map_mets = regexprep(mapModel.mets,'\[.\]$','');
elseif length(unique(mapModel.mets)) == length(mapModel.mets)
    % If the met IDs are not all unique, assume the compartment has already
    % been removed; otherwise, assume that the last character of the ID is
    % the compartment abbreviation, and remove it.
    map_mets = regexprep(mapModel.mets,'.$','');
end

% obtain unique list of compartment-free met IDs, and only keep field
% entries associated with this unique set of metabolites
[~,uniq_met_ind] = unique(map_mets);
map_mets = map_mets(uniq_met_ind);
mapModel.metNamesAlt = mapModel.metNamesAlt(uniq_met_ind,:);
for i = 1:length(metFields)
    % also remove rows correponding to non-unique mets from metFields
    mapModel.(metFields{i}) = mapModel.(metFields{i})(uniq_met_ind,:);
end

% initialize variables
matchScores = [];
name_stage = 1;
f = 1;

% map metabolites along each field in metFields
while f <= length(metFields)

    fprintf('Mapping metabolites via %s... ',metFields{f});

    % extract IDs from model field
    switch metFields{f}

        case 'mets'

            % If comparing the "mets" field, remove the compartment
            if strcmp(refModel.mets{1}(end),']')
                ref_mets = regexprep(refModel.mets,'\[.\]$','');
            elseif length(unique(refModel.mets)) == length(refModel.mets)
                ref_mets = regexprep(refModel.mets,'.$','');
            end
            ref_ids = lower(ref_mets);
            map_ids = lower(map_mets);
            is_name = false;
            ids = repmat({''},size(refModel.mets));  % initialize matches

        case 'metNames'

            if name_stage == 1

                % METNAMES STAGE 1: Map mets via metNames field
                fprintf('(stage 1) ');
                is_name = true;
                ref_ids = refModel.metNames;
                map_ids = [mapModel.metNames,mapModel.metNamesAlt];
                ids = repmat({''},size(refModel.mets));  % initialize matches

            elseif name_stage == 2

                % METNAMES STAGE 2: Map remaining mets via metNamesAlt field
                fprintf('(stage 2) ');
                is_name = true;
                ref_ids = [refModel.metNames,refModel.metNamesAlt];
                map_ids = [mapModel.metNames,mapModel.metNamesAlt];

            elseif name_stage == 3

                % METNAMES STAGE 3: Make manual changes to metNames, and map.
                % These changes are specific to the mapping between HMR2 and
                % Recon3D, and are based on observations in naming differences.
                fprintf('(stage 3) ');
                is_name = true;
                ref_ids = [refModel.metNames,refModel.metNamesAlt];
                map_ids = [mapModel.metNames,mapModel.metNamesAlt];

                % make name adjustments to Recon3D model if present
                if isfield(refModel,'modelID') && strcmpi(refModel.modelID,'Recon3D')
                    ref_ids = applyMetNameChanges(ref_ids);
                elseif isfield(mapModel,'modelID') && strcmpi(mapModel.modelID,'Recon3D')
                    map_ids = applyMetNameChanges(map_ids);
                end

            end

        otherwise

            % comparing by field other than "mets" or "metNames"
            is_name = false;  % not a met name
            ref_ids = refModel.(metFields{f});
            map_ids = mapModel.(metFields{f});
            ids = repmat({''},size(refModel.mets));  % initialize matches

    end

    % if map_ids contains multiple columns, flatten into a single
    % column, and repeat entries of mapModel mets to maintain alignment
    non_empty = ~cellfun(@isempty,map_ids);
    if size(map_ids,2) > 1
        repMets = arrayfun(@(i) repmat(map_mets(i),sum(non_empty(i,:),2),1),(1:length(map_mets))','UniformOutput',false);
        repMets = vertcat(repMets{:});
        map_ids = map_ids';
        map_ids = map_ids(non_empty');
    else
        map_ids = map_ids(non_empty);
        repMets = map_mets(non_empty);
    end

    % assemble mapping array
    met2id = [repMets,map_ids];

    % map metabolites
    ignore_mets = ~cellfun(@isempty,ids);  % for the name-matching process, ignore mets that have been mapped in previous name-matching stages
    ids(~ignore_mets) = matchIDs(ref_ids(~ignore_mets,:),met2id,is_name);
    fprintf('Done.\n');

    if strcmpi(metFields{f},'metNames') && name_stage < 3
        % if processing metNames, move to next stage
        name_stage = name_stage + 1;
    else

        % get indices of mapped mets in each of the respective models
        ref_met_inds = arrayfun(@(i) repmat(i,numel(ids{i}),1),(1:length(ids))','UniformOutput',false);
        ref_met_inds = vertcat(ref_met_inds{:});
        [~,map_met_inds] = ismember(vertcat(ids{:}),map_mets);

        % append matches to match score matrix
        matchScores(:,:,f) = accumarray(unique([ref_met_inds,map_met_inds],'rows'),metFieldWeights(f),[length(refModel.mets),length(map_mets)]);
        f = f + 1;  % proceed to next metField

    end

end

if strcmpi(mapMethod,'order')
    % only use the match from the top-scoring (earliest-listed) metField
    matchScores(matchScores < max(matchScores,[],3)) = 0;
end

% determine the total score for each match by adding scores among all fields
matchScores = sum(matchScores,3);

if ~strcmpi(mapMethod,'all')
    % For each metabolite, only keep the mapped IDs with the highest score.
    % If there is a tie, keep all IDs that are tied.
    matchScores(matchScores < max(matchScores,[],2)) = 0;
end

% retrieve mapModel met IDs and add to output model structure
metAltModelID = arrayfun(@(i) map_mets(matchScores(i,:) > 0),(1:length(refModel.mets))','UniformOutput',false);
model.metAltModelID = flattenCell(metAltModelID,true);

end

function chgNames = applyMetNameChanges(names)

% met name changes specific to Recon3D-HMR2 metabolite mapping
chgNames = regexprep(names,'Coenzyme A','CoA');
chgNames = regexprep(chgNames,'\(.*-Density Lipoprotein\)','');
chgNames = regexprep(chgNames,'13-Cis-Retinoyl Glucuronide','13-Cis-Retinoyl-beta-D-Glucuronide');
chgNames = regexprep(chgNames,"Cytidine-5'-Diphosphate",'CDP');
chgNames = regexprep(chgNames,'Human Liver Homolog','');
chgNames = regexprep(chgNames,'Glutathionyl-Leuc4','glutathionyl-leukotriene C4');
chgNames = regexprep(chgNames,'Omega-Cooh-Tetranor-Leukotriene E3','Omega-Cooh-Tetranor-LTE3');

end

function metID = matchIDs(ids,met2id,is_name)

if is_name
    % make all met names lowercase, and remove special characters (non-word or
    % digit characters, e.g., dashes, parentheses, spaces, etc.)
    met2id(:,2) = lower(regexprep(met2id(:,2),'[^a-zA-Z0-9]',''));
    ids = lower(regexprep(ids,'[^a-zA-Z0-9]',''));
end

% compress refModel met IDs into single column of nested entries
if size(ids,2) > 1
    ids = nestCell(ids,true);
end

% find empty met ID indices to ignore
ignore_ind = cellfun(@isempty,ids);

% map refModel met IDs to mapModel met IDs
metID = repmat({''},size(ids,1),1);
metID(~ignore_ind) = cellfun(@(x) unique(met2id(ismember(met2id(:,2),x),1)),ids(~ignore_ind),'UniformOutput',false);

end

author: SysBioChalmers

repo: Human-GEM


text:
function [netSecretionFluxes, netUptakeFluxes, Y, modelStats, summary, statistics, modelsOK] = mgPipe(modPath, abunFilePath, computeProfiles, resPath, dietFilePath, infoFilePath, biomasses, hostPath, hostBiomassRxn, hostBiomassRxnFlux, figForm, numWorkers, rDiet, pDiet, lowerBMBound, upperBMBound, includeHumanMets, adaptMedium, pruneModels)
% mgPipe is a MATLAB based pipeline to integrate microbial abundances
% (coming from metagenomic data) with constraint based modeling, creating
% individuals' personalized models.
% The pipeline is divided in 3 parts:
% [PART 1] Analysis of individuals' specific microbes abundances are computed.
% [PART 2]: 1 Constructing a global metabolic model (setup) containing all the
% microbes listed in the study. 2 Building individuals' specific models
% integrating abundance data retrieved from metagenomics. For each organism,
% reactions are coupled to the objective function.
% [PART 3] Simulations under different diet regimes.
%
% USAGE:
%       [netSecretionFluxes, netUptakeFluxes, Y, modelStats, summary, statistics, modelsOK] = mgPipe(modPath, abunFilePath, computeProfiles, resPath, dietFilePath, infoFilePath, biomasses, hostPath, hostBiomassRxn, hostBiomassRxnFlux, figForm, numWorkers, rDiet, pDiet, lowerBMBound, upperBMBound, includeHumanMets, adaptMedium, pruneModels)
%
% INPUTS:
%    modPath:                char with path of directory where models are stored
%    abunFilePath:           char with path and name of file from which to retrieve abundance information
%    computeProfiles:        boolean defining whether flux variability analysis to
%                            compute the metabolic profiles should be performed.
%    resPath:                char with path of directory where results are saved
%    dietFilePath:           char with path of directory where the diet is saved
%    infoFilePath:           char with path to stratification criteria if available
%    biomasses:              Cell array containing names of biomass objective functions
%                            of models to join. Needs to be the same length as
%                            the length of models in the abundance file.
%    hostPath:               char with path to host model, e.g., Recon3D (default: empty)
%    hostBiomassRxn:         char with name of biomass reaction in host (default: empty)
%    hostBiomassRxnFlux:     double with the desired flux through the host
%                            biomass reaction (default: zero)
%    figForm:                format to use for saving figures
%    numWorkers:             integer indicating the number of cores to use for parallelization
%    rDiet:                  boolean indicating if to enable also rich diet simulations (default: 'false')
%    pDiet:                  boolean indicating if to enable also personalized diet simulations (default: 'false')
%    lowerBMBound:           lower bound on community biomass (default=0.4)
%    upperBMBound:           upper bound on community biomass (default=1)
%    includeHumanMets:       boolean indicating if human-derived metabolites
%                            present in the gut should be provided to the models (default: true)
%    adaptMedium:            boolean indicating if the medium should be adapted through the
%                            adaptVMHDietToAGORA function or used as is (default=true)
%    pruneModels:            boolean indicating whether reactions that do not carry flux on the
%                            input diet should be removed from the microbe models.
%                            Recommended for large datasets (default: false)
%
% OUTPUTS:
%    init:                   status of initialization
%    netSecretionFluxes:     Net secretion fluxes by microbiome community models
%    netUptakeFluxes:        Net uptake fluxes by microbiome community models
%    Y:                      Classical multidimensional scaling
%    modelStats:             Reaction and metabolite numbers for each model
%    summary:                Table with average, median, minimal, and maximal
%                            reactions and metabolites
%    statistics:             If info file with stratification is provided, will
%                            determine if there is a significant difference.
%    modelsOK:               Boolean indicating if the created microbiome models
%                            passed verifyModel. If true, all models passed.
%
% AUTHORS:
%   - Federico Baldini, 2017-2018
%   - Almut Heinken, 07/20: converted to function
%   - Almut Heinken, 01/21: adapted inputs

%% PIPELINE: [PART 1]
% If PART1 was already
% computed: if the associated file is already present in the results folder its
% execution is skipped else its execution starts

[mapP]=detectOutput(resPath,'mapInfo.mat');

if ~isempty(mapP)
    s= 'mapping file found: loading from resPath and skipping [PART1] analysis';
    disp(s)
    load(strcat(resPath,'mapInfo.mat'))
end

[mapP]=detectOutput(resPath,'mapInfo.mat');
if isempty(mapP)
    % check if AGORA models are in path
    if ~exist('modPath', 'var') || ~exist(modPath, 'dir')
        error('modPath is not defined. Please set the path of the model directory.');
    else
        if ~exist(modPath, 'dir')
            error(['modPath (' modPath ') does not exist.']);
        end
    end

    abundance = table2cell(readtable(abunFilePath));

    % The number of microbeNames, their names, the number of samples and their identifiers
    % are automatically detected from the input file.
    [sampNames,microbeNames,exMets]=getIndividualSizeName(abunFilePath,modPath);

    % remove rows of organisms that are not present in any sample
    if contains(version,'(R202') % for Matlab R2020a and newer
        microbeNames(sum(cell2mat(abundance(:,2:end)),2)<0.0000001,:)=[];
        abundance(sum(cell2mat(abundance(:,2:end)),2)<0.0000001,:)=[];
    else
        microbeNames(sum(str2double(abundance(:,2:end)),2)<0.0000001,:)=[];
        abundance(sum(str2double(abundance(:,2:end)),2)<0.0000001,:)=[];
    end

    % Extracellular spaces simulating the lumen are built and stored for
    % each microbe.
    [activeExMets,couplingMatrix]=buildModelStorage(microbeNames, modPath, dietFilePath, adaptMedium, includeHumanMets, numWorkers, pruneModels, biomasses);

    % Computing reaction abundance and reaction presence
    [ReactionAbundance,ReactionPresence] = fastCalculateReactionAbundance(abunFilePath, modPath, {}, numWorkers);
    writetable(cell2table(ReactionAbundance'),[resPath filesep 'ReactionAbundance.csv'], 'WriteVariableNames', false);
    writetable(cell2table(ReactionPresence'),[resPath filesep 'ReactionPresence.csv'], 'WriteVariableNames', false);

    % Computing subsystem abundance
    subsystemAbundance = calculateSubsystemAbundance([resPath filesep 'ReactionAbundance.csv']);
    writetable(cell2table(subsystemAbundance),[resPath filesep 'SubsystemAbundance.csv'], 'WriteVariableNames', false);

    % plot subsystem abundance
    data=cell2mat(subsystemAbundance(2:end,2:end));
    xlabels=subsystemAbundance(1,2:end);
    ylabels=subsystemAbundance(2:end,1);
    % remove low abundance subsystems
    ylabels(find(sum(data,2)<1),:)=[];
    data(find(sum(data,2)<1),:)=[];
    figure;
    imagesc(data)
    colormap('hot')
    colorbar
    if length(xlabels)<30
        set(gca,'xtick',1:length(xlabels));
        xticklabels(xlabels);
        xtickangle(90)
    end

    if length(xlabels)<30
        set(gca,'xtick',1:length(xlabels));
        xticklabels(xlabels);
        xtickangle(90)
    end
    set(gca,'ytick',1:length(ylabels));
    yticklabels(ylabels);
    ax=gca;

    if length(ylabels)<50
        ax.YAxis.FontSize = 8;
    else
        ax.YAxis.FontSize = 6;
    end

    set(gca,'TickLabelInterpreter', 'none');
    title('Relative reaction abundances summarized by subsystem')
    print(strcat(resPath, 'Subsystem_abundances'), figForm)

    % save mapping info
    save([resPath filesep 'mapInfo.mat'], 'mapP', 'exMets', 'activeExMets', 'sampNames', 'microbeNames', 'couplingMatrix', 'abundance','-v7.3')
end

%end of trigger for Autoload
%% PIPELINE: [PART 2.1]

% Now we detect from the content of the results folder If PART2 was already
% computed: if the associated file is already present in the results folder its
% execution is skipped else its execution starts

% If desired, a model of the host (e.g., Recon3D) can also be joined with
% the microbiome models.
if ~isempty(hostPath)
    % workaround for models that give an error in readCbModel
    try
        host = readCbModel(hostPath);
    catch
        warning('Host model could not be read through readCbModel. Consider running verifyModel.')
        modelStruct=load(hostPath);
        getfn=fieldnames(modelStruct);
        host=modelStruct.(getfn{1});
    end
else
    host = {};
end

%end of trigger for Autoload

% A  model joining all the reconstructions contained in the study
% will be created in this section. This model will be later used, integrating
% abundances coming from the metagenomic sequencing, to derive the different microbiota
% models. The result of this section will be automatically saved in the results
% folder.

%% PIPELINE: [PART 2]
% Now we will create the different microbiota models integrating the given abundances.
% Coupling constraints and personalized "cumulative biomass" objective functions
% are also added. Models that are already existent will not be recreated, and
% new microbiota models will be saved in the results folder.

% set parallel pool
if numWorkers > 1
    poolobj = gcp('nocreate');
    if isempty(poolobj)
        parpool(numWorkers)
    end
end

% create a separate setup model for each sample
% define what counts as zero abundance
tol=0.0000001;

clear('microbeNames','activeExMets','abundance')

if length(sampNames)>50
    steps=50;
else
    steps=length(sampNames);
end
% proceed in batches for improved effiency
cnt=1;
modelsWithErrors={};

for j=1:steps:length(sampNames)
    if length(sampNames)-j>=steps-1
        endPnt=steps-1;
    else
        endPnt=length(sampNames)-j;
    end
    getErrors={};

    parfor i=j:j+endPnt
        % Each personalized model will be created separately.
        % get the list of models for each sample and remove the ones not in
        % this sample

        % check if model already exists
        if ~isempty(host)
            mId = strcat('host_microbiota_model_samp_', sampNames{i,1}, '.mat');
        else
            mId = strcat('microbiota_model_samp_', sampNames{i,1}, '.mat');
        end
        if ~isfile(mId)
            mappingData=load([resPath filesep 'mapInfo.mat'])
            microbeNamesSample = mappingData.microbeNames;
            couplingMatrixSample = mappingData.couplingMatrix;
            abunRed=mappingData.abundance(:,i+1);
            abunRed=[mappingData.abundance(:,1),abunRed];

            microbeNamesSample(cell2mat(abunRed(:,2)) < tol,:)=[];
            couplingMatrixSample(cell2mat(abunRed(:,2)) < tol,:)=[];
            abunRed(cell2mat(abunRed(:,2)) < tol,:)=[];
            setupModel = fastSetupCreator(exMets, microbeNamesSample, host);

            % create personalized models for the batch
            createdModel=createPersonalizedModel(abunRed,resPath,setupModel,sampNames(i,1),microbeNamesSample,couplingMatrixSample,host,hostBiomassRxn);
            results=verifyModel(createdModel{1});
            getErrors{i} = results;
        end
    end
    for i=j:j+endPnt
        if length(getErrors) >= i
            if ~isempty(getErrors{i})
                modelsWithErrors{cnt,1} = sampNames{i,1};
                cnt=cnt+1;
            end
        end
    end
end

if isempty(modelsWithErrors)
    modelsOK = true;
else
    modelsOK = false;
end

%% PIPELINE: [PART 3]
%
% In this phase, for each microbiota model, a diet, in the form of set constraints
% to the exchanges reactions of the diet compartment, is integrated. Flux Variability
% analysis for all the exchange reactions of the diet and fecal compartment is
% also computed and saved in a file called "simRes".

load([resPath filesep 'mapInfo.mat'])
[exchanges, netProduction, netUptake, presolve, infeasModels] = microbiotaModelSimulator(resPath, exMets, sampNames, dietFilePath, hostPath, hostBiomassRxn, hostBiomassRxnFlux, numWorkers, rDiet, pDiet, computeProfiles, lowerBMBound, upperBMBound, includeHumanMets, adaptMedium);
% Finally, NMPCs (net maximal production capability) are computed in a metabolite
% resolved manner and saved in a comma delimited file in the results folder. NMPCs
% indicate the maximal production of each metabolite and are computing summing
% the maximal secretion flux with the maximal uptake flux. Similarity of metabolic
% profiles (using the different NMPCs as features) between individuals are also
% evaluated with classical multidimensional scaling.

if computeProfiles
    [netSecretionFluxes, netUptakeFluxes, Y] = mgSimResCollect(resPath, sampNames, exchanges, rDiet, pDiet, infoFilePath, netProduction, netUptake, figForm);
else
    netSecretionFluxes={};
    netUptakeFluxes={};
    Y=[];
    delete('simRes.mat','intRes.mat')
end

% get stats on microbiome models-number of reactions and metabolites
for i=1:length(sampNames)
    modelNames{i}=['microbiota_model_samp_' sampNames{i}];
end

close all

if ~isempty(infoFilePath)
    [modelStats,summary,statistics]=retrieveModelStats(resPath, modelNames, abunFilePath, numWorkers, infoFilePath);
else
    [modelStats,summary,statistics]=retrieveModelStats(resPath, modelNames, abunFilePath, numWorkers);
end
writetable(cell2table(modelStats),[resPath filesep 'ModelStatistics.csv'], 'WriteVariableNames', false);
writetable(cell2table(summary),[resPath filesep 'ModelStatsSummary.csv'], 'WriteVariableNames', false);
if ~isempty(statistics)
    writetable(cell2table(statistics),[resPath filesep 'ModelStatsStratification.csv'], 'WriteVariableNames', false);
end

end

author: opencobra

repo: cobratoolbox


text:
function [logbf,t] = spm_bms_ttest (y,prior)
% Log Bayes Factor against null for one sample t-test
% FORMAT [logbf,t] = spm_bms_ttest (y,prior)
%
% y         [N x 1] data vector
% prior     'jzs' (default), 'unit'
%
% logbf     Log Bayes Factor
% t         t-statistic
%
% Default prior is 'jzs'. The different priors are described in
% [1] Rouder et al. (2009) Psych Bull Rev, 16(2),225-237.
%
% These tests consider the quantity d = mean / standard deviation.
% The unit information prior, 'unit', uses a zero mean unit variance
% Gaussian prior over d (the prior variance of d, sig_d^2, is unity).
%
% The (Jeffreys-Zellner-Snow) JZS prior, 'jzs', uses a Cauchy over d
% and an inverse chi^2 over sig_d^2.
%_______________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_bms_ttest.m 6038 2014-06-04 15:22:42Z will $

try
    pr=prior;
catch
    pr='jzs';
end

N=length(y);
m=mean(y);
s=std(y);
sem=s/sqrt(N);
t=m/sem;
v=N-1;

switch pr,
    case 'unit',
        rs=1;
    case 'jzs',
        logbf = jzs_logbf(t,v,N);
        return
    otherwise
        disp('Unknown prior in bayes_ttest');
        return
end

% For Unit Information prior, 'unit':
%
% The minimum value of logbf obtains when t=0
% This minimum value is logbf=-0.5*log(1+N*r)
% For r=1 logbf=-3 can be gotten with N=405 samples.

% See note 5 at end of [1]
t2=t^2;
nr=1+N*rs^2;
log_num=-0.5*N*log(1+t2/v);
log_denom1=-0.5*log(nr);
log_denom2=-0.5*N*log(1+t2/(nr*v));

logbf=log_num-log_denom1-log_denom2;
logbf=-logbf;

%-------------------------------------------------------
function [logbf] = jzs_logbf(t,v,N)
% JZS Log Bayes factor
%
% See equation (1) in
% Rouder et al. (2009) Psych Bull Rev, 16(2),225-237.

t2=t^2;
num=(1+t2/v)^(-0.5*(v+1));

gmax=1000;
denom=quad(@(g)jzs_integrand(g,t2,N,v),0,gmax);

logbf=log(num/denom);
logbf=-logbf;

%--------------------------------------------------------
function [j] = jzs_integrand(g,t2,N,v)

ng=1+N*g;
t1=ng.^(-0.5);
t2=(1+t2./(ng*v)).^(0.5*(v+1));
t3=(2*pi)^(-0.5)*(g.^(-1.5)).*exp(-1./(2*g));
j=t1.*t2.*t3;

author: spm


text:
clear all;close all;clc

M_dir = GetFishDirectories();
save_masterdir = GetNestedDataDir();

range_fish = GetFishRange();

for i_fish = range_fish,
    disp(['i_fish = ' num2str(i_fish)]);

    data_dir = M_dir{i_fish};

    %% init
    save_dir = fullfile(save_masterdir,['subject_' num2str(i_fish)]);
    load(fullfile(save_dir,'CoreInfo.mat'));
    load(fullfile(save_dir,'OptionalInfo.mat'));

    %% get anatomy projections
    % x-y view
    im = max(anat_stack,[],3);
    out=imNormalize99(im);
    anat_yx = repmat(out,[1 1 3]);

    % y-z view
    im = squeeze(max(anat_stack,[],2));
    out=imNormalize99(im);
    anat_yz = repmat(out,[1 1 3]);

    % x-z view
    im = squeeze(max(anat_stack,[],1));
    out=imNormalize99(im);
    out = flipud(out'); %%%% empirically necessary...
    anat_zx = repmat(out,[1 1 3]);

    %% [Next 2 cells] delete anatomically out-of-bound cells by hand
    % only execute once (manual choice 1):
    cHolder_Anat = []; % collect cIX of all out-of-bound cells

    %% (optional) step 1: set limit to only plot ventral layer, helps to find outliers within that layer
    % draw cells first
    % (choose start and stop index to draw; low numbers ~ ventral)
    I_start = 1;
    I_stop = 1*round(numcell_full/10);

    % plot those cells
    cIX = I_start:I_stop;
    gIX = round(cIX/1000)+1;
    figure('Position',[100 0 1300 900]);
    numK = round(cIX(end)/1000)+1;
    [~, dim_totimage] = BasicDrawCellsOnAnatProj(CellXYZ,cIX,gIX,numK,anat_yx,anat_yz);%,anat_zx,'hsv','full');

    % Here manually draw polygons around cells to discard.
    % Double click to connect last vertex to first vertex, then double click again within polygon to fix.
    % then click again to start drawing the next one
    % and when finished, break loop by pressing any key

    for i = 1:100, % NOMINAL LOOP, break manually
        h_poly_yx = impoly;
        wait(h_poly_yx); % double click to finalize position!
        % update finalized polygon in bright color
        setColor(h_poly_yx,[0 1 1]);

        A = sub2ind(size(anat_yx),CellXYZ(I_start:I_stop,1),CellXYZ(I_start:I_stop,2));
        MaskArray = createMask(h_poly_yx);
        B = find(MaskArray); % find indices of pixels within ROI
        cIX2 = find(ismember(A,B));
        cHolder_Anat = union(cHolder_Anat,cIX2);
        w = waitforbuttonpress;
        if w == 1, % press any key to break
            break;
        end
    end

    %% step 2: plot all cells
    I_start = 1; % needs to be 1 for correct indexing
    I_stop = numcell_full;

    % ...here until the end of the cell is the exact dupliate of the last cell...
    cIX = I_start:I_stop;
    gIX = round(cIX/1000)+1;
    figure('Position',[100 0 1300 900]);
    numK = round(cIX(end)/1000)+1;
    [tot_image, dim_totimage] = BasicDrawCellsOnAnatProj(CellXYZ,cIX,gIX,numK,anat_yx,anat_yz);%,anat_zx,'hsv','full');

    % Here manually draw polygons around cells to discard.
    % Double click to connect last vertex to first vertex, then double click again within polygon to fix.
    % then click again to start drawing the next one
    % and when finished, break loop by pressing any key

    for i = 1:100, % NOMINAL LOOP, break manually
        h_poly_yx = impoly;
        wait(h_poly_yx); % double click to finalize position!
        % update finalized polygon in bright color
        setColor(h_poly_yx,[0 1 1]);

        A = sub2ind(size(anat_yx),CellXYZ(I_start:I_stop,1),CellXYZ(I_start:I_stop,2));
        MaskArray = createMask(h_poly_yx);
        B = find(MaskArray); % find indices of pixels within ROI
        cIX2 = find(ismember(A,B));
        cHolder_Anat = union(cHolder_Anat,cIX2);
        w = waitforbuttonpress;
        if w == 1, % press any key to break
            break;
        end
    end

    %% find extra outliers on borders of image (can't always get with polygon)
    IX = find(CellXYZ(:,1)<15 | CellXYZ(:,1)> size(anat_stack,1)-15 ...
        | CellXYZ(:,2)<15 | CellXYZ(:,2)> size(anat_stack,2)-15);
    % for fish 9, eliminating bad plans: | CellXYZ(:,3)<4 | CellXYZ(:,3)> size(anat_stack,3)-3);
    cHolder_Anat = union(cHolder_Anat,IX);

    %% test Plot: all antomy outliers
    IX_inval_anat = cHolder_Anat; % rename
    cIX = IX_inval_anat;
    gIX = (1:length(cIX))';
    figure;
    BasicDrawCellsOnAnatProj(CellXYZ,cIX,gIX,numK,anat_yx,anat_yz);
    pause(2);

    %% test Plot: all remaining cells
    I_v_Holder = ones(1,numcell_full);
    I_v_Holder(IX_inval_anat) = 0;
    cIX = find(I_v_Holder);
    gIX = (1:length(cIX))';
    figure;
    BasicDrawCellsOnAnatProj(CellXYZ,cIX,gIX,numK,anat_yx,anat_yz);
    pause(2);

    %% ...and save
    temp = fullfile(data_dir,['Fish' num2str(i_fish) '_extrainfo_anat.mat']);
    save(temp,'IX_inval_anat');

    %%

    % save_masterdir = GetNestedDataDir();
    % save_dir = fullfile(save_masterdir,['subject_' num2str(i_fish)]);

    IX_inval = IX_inval_anat;
    filename = fullfile(save_dir,'OptionalInfo.mat');
    save(filename,'IX_inval','-append');

    filename = fullfile(save_dir,'AdditionalInfo.mat');
    save(filename,'IX_inval_anat','-append');
    disp('anat outliers saved');

end

author: xiuyechen

repo: FishExplorer


text:
function r = times(a,b)
%TIMES        Taylor multiplication  a .* b
%

% written  05/21/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

  K1 = getappdata(0,'INTLAB_TAYLOR_ORDER') + 1;

  if ~isa(a,'taylor')           % non-taylor times taylor
    r = b;
    if isa(a,'intval')
      r.t = intval(r.t);
    end
    m = prod(size(a));
    if m==1                     % non-taylor scalar .* taylor
      r.t = a*r.t;
    else                        % non-taylor array .* taylor
      n = prod(b.size);
      if n==1                   % non-taylor array .* taylor scalar
        r.size = size(a);
        r.t = repmat(a(:).',K1,1) .* repmat(b.t,1,m);
      else                      % non-taylor array .* taylor array
        if ~isequal(size(a),b.size)
          error('Taylor multiplication : dimensions not compatible')
        end
        r.t = repmat(a(:).',K1,1) .* b.t;
      end
    end
  elseif ~isa(b,'taylor')       % taylor times non-taylor
    r = a;
    if isa(b,'intval')
      r.t = intval(r.t);
    end
    m = prod(size(b));
    if m==1                     % taylor scalar .* non-taylor
      r.t = r.t*b;
    else                        % taylor array .* non-taylor
      n = prod(a.size);
      if n==1                   % taylor array .* non-taylor scalar
        r.t = repmat(a.t,1,m) .* repmat(b(:).',K1,1);
      else                      % non-taylor array .* taylor array
        if ~isequal(size(b),a.size)
          error('Taylor multiplication : dimensions not compatible')
        end
        r.t = a.t .* repmat(b(:).',K1,1);
      end
    end
  else                          % both factors taylor
    m = prod(a.size);
    n = prod(b.size);
    if m==1                     % taylor scalar .* taylor
      r = b;
      if isa(a,'intval')
        r.t = intval(r.t);
      end
      for j=1:K1
        r.t(j,:) = sum(repmat(a.t(1:j),1,n).*b.t(j:-1:1,:),1);
      end
    else                        % taylor array .* taylor
      r = a;
      if isa(b.t,'intval')
        r.t = intval(r.t);
      end
      if n==1                   % taylor array .* taylor scalar
        for j=1:K1
          r.t(j,:) = sum(a.t(1:j,:).*repmat(b.t(j:-1:1),1,m),1);
        end
      else                      % taylor array .* taylor array
        if ~isequal(a.size,b.size)
          error('Taylor multiplication : dimensions not compatible')
        end
        for j=1:K1
          r.t(j,:) = sum(a.t(1:j,:).*b.t(j:-1:1,:),1);
        end
      end
    end
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function varargout = ft_transform_headmodel(varargin)

% This function is a backward compatibility wrapper for existing MATLAB scripts
% that call a function that is not part of the FieldTrip toolbox any more.
%
% Please update your code to make it future-proof.

oldname = mfilename;
newname = 'ft_transform_geometry';

ft_warning('%s is only a backward compatibility wrapper, which will soon be removed. Please call %s instead.', upper(oldname), upper(newname));

[varargout{1:nargout}] = feval(newname, varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        Sample Data and Metch Registration Sessions        %
%                                                           %
%           by Qianqian Fang <q.fang at neu.edu>            %
%                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% first of all, make sure you've already add the metch root
% folder to your matlab/octave search path list

% load the sample data, where
%   no: the node coordinates of a surface mesh
%   el: the surface triangles
%   pt: the point cloud to be registered

disp('(*)First load the mesh and point cloud. Hit Enter to continue...');
pause;

load sampledata
if(exist('OCTAVE_VERSION')~=0)
	trimesh(el,no(:,1),no(:,2),no(:,3));
else
	cla;
	trisurf(el,no(:,1),no(:,2),no(:,3));
end
title('Metch toolbox demonstration');
axis equal;
hold on;

% use metch functions to perform the registration

pnum=size(pt,1);

% define a number of point pairs to initialize the registration
disp('(*)Create 4 mapping pairs to initialize the mapping. Hit Enter to continue...');
pause;

% select 4 land-marks on the point cloud (specified by their indicies)
ptidx=[4 107 1 190];
ptselected=pt(ptidx,:);

% find the corresponding land-marks on the mesh
meshidx=[3173 1715 156 1740];
meshselected=no(meshidx,:);

% calculate the affine mapping using these point pairs
[A0,b0]=affinemap(ptselected,meshselected)

disp('(*)Display the updated points. Hit Enter to continue...');
pause;

% a rough registration from the selected point pairs
points_after_initmap=(A0*pt'+repmat(b0(:),1,pnum))';
plot3(points_after_initmap(:,1),points_after_initmap(:,2),points_after_initmap(:,3),'r.');

disp('(*)Optimize the mapping matrix to fit the surface. Hit Enter to continue...');
pause;

% set pmask: if pmask(i) is -1, it is a free nodes to be optimized
%            if pmask(i) is 0, it is fixed
%            if pmask(i) is a positive number, it is the index of
%               the mesh node to map to

pmask=-1*ones(pnum,1);
pmask(ptidx)=meshidx;

% perform mesh registration with Gauss-Newton method using A0/b0
% as initial guess
[A,b,newpos]=regpt2surf(no,el,pt,pmask,A0,b0,ones(12,1),10);
A
b

disp('(*)Display the optimized point cloud. Hit Enter to continue...');
pause;

% update point cloud with the optimized mapping
points_after_optimize=(A*pt'+repmat(b(:),1,pnum))';

plot3(points_after_optimize(:,1),points_after_optimize(:,2),points_after_optimize(:,3),'g+');

disp('(*)Project the point cloud on the surface. Hit Enter to continue...');
pause;

% project the optimized point cloud onto the surface, and make
% sure the comformity

nv=nodesurfnorm(no,el);
[d2surf,cn]=dist2surf(no,nv,points_after_optimize);
[points_after_proj eid weights]=proj2mesh(no,el,points_after_optimize,nv,cn);

disp('(*)Display the final results. Hit Enter to continue...');
pause;

plot3(points_after_proj(:,1),points_after_proj(:,2),points_after_proj(:,3),'c*');

legend('surface mesh','points after initial map','points after optimized map',...
      'points after projection');

author: fangq


text:
function acq = acqimiqr_vbmc(Xs,vp,gp,optimState,fmu,fs2,fbar,vtot)
%VBMC_ACQIMIQR Integrated median interquantile range acquisition function.

u = 0.6745; % norminv(0.75)

if isempty(Xs)
    % Return acquisition function info struct
    acq.importance_sampling = true;
    acq.importance_sampling_vp = false;
    acq.log_flag = true;
    return;
elseif ischar(Xs)
    switch lower(Xs)
        case 'islogf1'
        % Importance sampling log base proposal (shared part)
            acq = fmu;
        case 'islogf2'
        % Importance sampling log base proposal (added part)
        % (Full log base proposal is fixed + added)
            fs = sqrt(fs2);
            acq = u*fs + log1p(-exp(-2*u*fs));
        case 'islogf'
        % Importance sampling log base proposal distribution
            fs = sqrt(fs2);
            acq = fmu + u*fs + log1p(-exp(-2*u*fs));
    end
    return;
end

% Different importance sampling inputs for different GP hyperparameters?
multipleinputs_flag = size(optimState.ActiveImportanceSampling.Xa,3) > 1;

% Xs is in *transformed* coordinates

[Nx,D] = size(Xs);
Ns = size(fmu,2);
Na = size(optimState.ActiveImportanceSampling.Xa,1);

% Estimate observation noise at test points from nearest neighbor
[~,pos] = min(sq_dist(bsxfun(@rdivide,Xs,optimState.gplengthscale),gp.X_rescaled),[],2);
sn2 = gp.sn2new(pos);
% sn2 = min(sn2,1e4);
ys2 = fs2 + sn2;    % Predictive variance at test points

if multipleinputs_flag
    Xa = zeros(Na,D);
else
    Xa = optimState.ActiveImportanceSampling.Xa;
end
acq = zeros(Nx,Ns);

%% Compute integrated acquisition function via importance sampling

for s = 1:Ns
    hyp = gp.post(s).hyp;
    L = gp.post(s).L;
    Lchol = gp.post(s).Lchol;
    sn2_eff = 1/gp.post(s).sW(1)^2;

    if multipleinputs_flag
        Xa(:,:) = optimState.ActiveImportanceSampling.Xa(:,:,s);
    end

    % Compute cross-kernel matrix Ks_mat
    if gp.covfun(1) == 1    % Hard-coded SE-ard for speed
        ell = exp(hyp(1:D))';
        sf2 = exp(2*hyp(D+1));
        Ks_mat = sq_dist(gp.X*diag(1./ell),Xs*diag(1./ell));
        Ks_mat = sf2 * exp(-Ks_mat/2);

        Ka_mat = sq_dist(Xa*diag(1./ell),Xs*diag(1./ell));
        Ka_mat = sf2 * exp(-Ka_mat/2);

        %Kax_mat = sq_dist(Xa*diag(1./ell),gp.X*diag(1./ell));
        %Kax_mat = sf2 * exp(-Kax_mat/2);
        Kax_mat(:,:) = optimState.ActiveImportanceSampling.Kax_mat(:,:,s);
    else
        error('Other covariance functions not supported yet.');
    end

    if Lchol
        C = Ka_mat' - Ks_mat'*(L\(L'\Kax_mat'))/sn2_eff;
    else
        C = Ka_mat' + Ks_mat'*(L*Kax_mat');
    end

    tau2 = bsxfun(@rdivide,C.^2,ys2(:,s));
    s_pred = sqrt(max(bsxfun(@minus,optimState.ActiveImportanceSampling.fs2a(:,s)',tau2),0));

    lnw = optimState.ActiveImportanceSampling.lnw(s,:);

    zz = bsxfun(@plus,lnw,u*s_pred + log1p(-exp(-2*u*s_pred)));
    lnmax = max(zz,[],2);
    acq(:,s) = log(sum(exp(bsxfun(@minus,zz,lnmax)),2)) + lnmax;
end

if Ns > 1
    M = max(acq,[],2);
    acq = M + log(sum(exp(bsxfun(@minus,acq,M)),2)/Ns);
end

end

%SQ_DIST Compute matrix of all pairwise squared distances between two sets
% of vectors, stored in the columns of the two matrices, a (of size n-by-D)
% and b (of size m-by-D).
function C = sq_dist(a,b)

n = size(a,1);
m = size(b,1);
mu = (m/(n+m))*mean(b,1) + (n/(n+m))*mean(a,1);
a = bsxfun(@minus,a,mu); b = bsxfun(@minus,b,mu);
C = bsxfun(@plus,sum(a.*a,2),bsxfun(@minus,sum(b.*b,2)',2*a*b'));
C = max(C,0);

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
% SENSE_recon_demo.m
% - demonstrates use of SENSE_recon.m, which implements the SENSE
% reconstruction
% - figure shows how Tikhonov regularization mitigates the increase in error
% as noise variance increases
% - applies complex AWGN to kspace
%
% 2012-06-08, Mai Le

% number of coils
nc = 4;
% factor of undersampling
np = 2;
% dimension to undersample
reduced_dim = 1;
regularizer1 = 'none';
regularizer2 = 'tikhonov';
figs = 0;

image = double(imread('./data/mri/brainweb_t1.jpg','jpg'));

dims = size(image);
% introduce planar phase
[xx,yy] = ndgrid(-dims(1)/2:dims(1)/2-1,-dims(2)/2:dims(2)/2-1);
ph_max = pi/5;
ph = ph_max*xx/(dims(1)/2)+ph_max*yy/(dims(2)/2);
image = image.*exp(1i*ph);

% generate sensitivity maps
smap = mri_sensemap_sim('nx',dims(1),'ny',dims(2),'ncoil',nc);

im_rep = repmat(image,[1 1 nc]);

% apply sensitivity maps
mapped_im = im_rep.*smap;

% throw away k-space data
im_fft = [];
reduced_fft = [];
for ii = 1:nc
    im_fft(:,:,ii) = fft2(mapped_im(:,:,ii));
    if (reduced_dim == 1)
        reduced_fft(:,:,ii) = im_fft(1:np:end,:,ii);
    else
        reduced_fft(:,:,ii) = im_fft(:,1:np:end,ii);
    end
end
%% perform SENSE reconstruction for varying levels of noise
sigmas = 0:25:150;
for ii = 1:length(sigmas)
    % introduce complex AWGN
    sigma = sigmas(ii);
    noisy_reduced_fft = reduced_fft + sigma*randn(size(reduced_fft)) + i*sigma*randn(size(reduced_fft));

    % SENSE reconstruction
    recon_im1(:,:,ii) = SENSE_recon(smap, noisy_reduced_fft, reduced_dim, np, regularizer1, figs);
    recon_im2(:,:,ii) = SENSE_recon(smap, noisy_reduced_fft, reduced_dim, np, regularizer2, figs);

    error1(ii) = sqrt((sum(sum((real(image)-real(recon_im1(:,:,ii))).^2))+sum(sum((imag(image)-imag(recon_im1(:,:,ii))).^2)))/prod(dims));
    error2(ii) = sqrt((sum(sum((real(image)-real(recon_im2(:,:,ii))).^2))+sum(sum((imag(image)-imag(recon_im2(:,:,ii))).^2)))/prod(dims)) ;
end
%% plot error as a function of noise variance
figure; plot(sigmas.^2, error1,'k');
hold on; plot(sigmas.^2, error2, 'r'); legend('no regularization', 'tikhonov');
xlabel('variance of complex AWGN in kspace');
ylabel('RMS error in reconstructed images');
title('effect of Tikhonov regularization on noisy SENSE reconstruction')
%% plot reconstructed image at sigma(ii)
ii = 4; % change as desired
figure; subplot(331); imshow(real(image),[]); colorbar;
title('real of original image');
subplot(332); imshow(imag(image),[]); colorbar;
title(['imag of original image, noise variance: ' num2str(sigmas(ii).^2)]);
subplot(334); imshow(real(recon_im1(:,:,ii)),[]); colorbar;
title('real of recon image, no reg');
subplot(335); imshow(imag(recon_im1(:,:,ii)),[]); colorbar;
title('imag of recon image, no reg');
subplot(336); imshow(abs(real(image)-real(recon_im1(:,:,ii)))+abs(imag(image)-imag(recon_im1(:,:,ii))),[]); colorbar;
title('abs(differences), no reg');
subplot(337); imshow(real(recon_im2(:,:,ii)),[]); colorbar;
title('real of recon image, Tikhonov');
subplot(338); imshow(imag(recon_im2(:,:,ii)),[]); colorbar;
title('imag of recon image, Tikhonov');
subplot(339); imshow(abs(real(image)-real(recon_im2(:,:,ii)))+abs(imag(image)-imag(recon_im2(:,:,ii))),[]); colorbar;
title('abs(differences), Tikhonov');

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function gt = Kernel_TVHP( times, landmarks, sigma, type )

dt = repmat(landmarks(:), [1, length(times(:))])...
    - repmat(times(:)', [length(landmarks(:)), 1]);

switch type
    case 'gauss'
        gt = exp(-(dt.^2)./(2*sigma^2));

    case 'exp'
        gt = exp(sigma.*dt);
        gt(gt>1) = 0;

end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
% This script file load a data set using fscanf
% The default reads Northern California Hypoellipse Format
%

report_this_filefun(mfilename('fullpath'));
disp('Please make sure the has 88 characters for each line')
disp('and all blanks have been substituted by zeros')

% Lets make sure the file is closed...
safe_fclose(fid);

% reset paramteres
a = []; b = []; n = 0;

if inda == 1
    % initial selection option
    tmin   = 1910.0001;
    tmax   = 2010;
    lonmin = -180.0;
    lonmax =  180.0;
    latmin =  -90.00;
    latmax =  90.0 ;
    Mmin   = -4.0;
    Mmax   = 10.;
    mindep = -10;
    maxdep = 700;

    % call the pre-selection window
    call = 'myloadncal2000'; % this callback is needed so presel know where to return to
    presel
    return
end

% open the file and read 10000 lines at a time
[file1,path1] = uigetfile([ '*.dat'],' Earthquake Datafile');
if length(file1) >1
    fid = fopen([path1 file1],'r') ;
else
    disp('Data import canceled'); return
end

while  ferror(fid) == ''
    n = n+1;
    % this is the new Y2K compliant fornat
    % vari name   yr mo da hr mi se lat   la  lon    lo de ma1     ma he   hz
    % variabl #   1  2  3  4  5  6  7      8  9      10 11 12      13 14   15
    % position    2  4  6  8  10 14 16 17  21 24 25  29 34 36 67   69 84   88
    l = fscanf(fid,'%4d%2d%2d%2d%2d%4d%2d%*1c%4d%3d%*1c%4d%5f%2d%*32c%2d%*15c%4d',...
        [14 10000]) ;
    %if ferror(fid) ~= '' ; break; end

    b = [ -l(9,:)-l(10,:)/6000 ; l(7,:)+l(8,:)/6000 ; l(1,:);l(2,:);l(3,:);
        l(13,:)/10;l(11,:)/100;l(4,:);l(5,:); l(14,:)/100;l(12,:)/10];
    b = b';
    l =  b(:,6) >= Mmin & b(:,1) >= lonmin & b(:,1) <= lonmax & ...
        b(:,2) >= latmin & b(:,2) <= latmax & b(:,3) <= tmax  & ...
        b(:,3) >= tmin  ;
    a = [a ; b(l,:)];

    disp([ num2str(n*10000) ' earthquakes scanned, ' num2str(length(a)) ' EQ found'])
    if max(b(:,3)) >  tmax ; break; end

end
ferror(fid)
fclose(fid);

% Convert the third column into time in decimals
if length(a(1,:))== 7
    a.Date = decyear(a(:,3:5));
elseif length(a(1,:))>=9       %if catalog includes hr and minutes
    a.Date = decyear(a(:,[3:5 8 9]));
end

dep1 = 0.3*(max(a.Depth)-min(a.Depth))+min(a.Depth);
dep2 = 0.6*(max(a.Depth)-min(a.Depth))+min(a.Depth);
dep3 = max(a.Depth);

stri1 = [file1];
tim1 = minti;
tim2 = maxti;
minma2 = minma;
maxma2 = maxma;
minde = min(a.Depth);
maxde = max(a.Depth);
rad = 50.;
ic = 0;
ya0 = 0.;
xa0 = 0.;
iwl3 = 1.;
step = 3;
t1p(1) = 80.;
t2p(1) = 85.;
t3p(1) = 90.;
t4p(1) = 93.;
tresh = 10;

% save the data
[file1,path1] = uiputfile(fullfile(hodi, 'eq_data', '*.mat'), 'Save Earthquake Datafile');
sapa2 = ['save ' path1 file1 ' a'];
if length(file1) > 1; eval(sapa2);end

% call the map window
update(mainmap())

author: CelsoReyes


text:
% Copyright (C) 2012 Rik Wehbring
% Copyright (C) 1995-2012 Kurt Hornik
%
% This file is part of Octave.
%
% Octave 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.
%
% Octave 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 Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

% -*- texinfo -*-
% @deftypefn {Function File} {} finv (@var{x}, @var{m}, @var{n})
% For each element of @var{x}, compute the quantile (the inverse of
% the CDF) at @var{x} of the F distribution with @var{m} and @var{n}
% degrees of freedom.
% @end deftypefn

% Author: KH <Kurt.Hornik@wu-wien.ac.at>
% Description: Quantile function of the F distribution

function inv = finv (x, m, n)

if (nargin ~= 3)
  print_usage ();
end

if (~isscalar (m) || ~isscalar (n))
  [retval, x, m, n] = common_size (x, m, n);
  if (retval > 0)
    error ('finv: X, M, and N must be of common size or scalars');
  end
end

if ~(isreal (x) || isreal (m) || isreal (n))
  error ('finv: X, M, and N must not be complex');
end

if (isa (x, 'single') || isa (m, 'single') || isa (n, 'single'))
  inv = NaN (size (x), 'single');
else
  inv = NaN (size (x));
end

k = (x == 1) & (m > 0) & (m < Inf) & (n > 0) & (n < Inf);
inv(k) = Inf;

k = (x >= 0) & (x < 1) & (m > 0) & (m < Inf) & (n > 0) & (n < Inf);
if (isscalar (m) && isscalar (n))
  inv(k) = ((1 ./ betainv (1 - x(k), n/2, m/2) - 1) * n / m);
else
  inv(k) = ((1 ./ betainv (1 - x(k), n(k)/2, m(k)/2) - 1)...
    .* n(k) ./ m(k));
end

end


text:
function prtUtilGraphExplore(self)

locs = self.plotLocations;
graph = self.graph;
names = self.nodeNames;

if isempty(locs)
    self = self.optimizePlotLocations;
    locs = self.plotLocations;
end

hF = figure;
hA = axes;
axis(hA,[-.5 .5 -.5 .5]);

nLines = 500;
lineRndStdDev = 0;

self.plot;

hz = zoom;
hp = pan;
set(hp,'ActionPostCallback',@(e,v)mypostcallback);
set(hz,'ActionPostCallback',@(e,v)mypostcallback);
set(hp,'ActionPreCallback',@(e,v)myprecallback);
set(hz,'ActionPreCallback',@(e,v)myprecallback);
hText = [];
hLines = [];

    function myprecallback(e,v)
        try
            delete(hText);
            hText = [];
            delete(hLines);
            hLines = [];
        catch ME
            disp(ME)
        end
    end

    function mypostcallback(e,v)
        % mypostcallback(e,v)
        v = axis;
        in = locs(:,1) > v(1) & locs(:,1) < v(2) & locs(:,2) > v(3) & locs(:,2) < v(4);

        cLocs = locs(in,:);
        cGraph = graph(in,in);
        cD = diag(cGraph);
        cNames = names(in);

        %Choose nLines lines to plot, at random... Should we use degree
        %here? Note, all this is done weirdly to make plotting super fast.
        [linkI,linkJ] = find(cGraph);
        plotLocs = [];
        indices = randperm(length(linkI));
        for i = indices(1:min(length(indices),nLines));
            plotLocs = cat(1,plotLocs,cLocs([linkI(i),linkJ(i)],:),[nan nan]);
        end

        % If there are any connections to plot, plot them:
        if ~isempty(plotLocs)
            plotLocs = plotLocs + randn(size(plotLocs))*lineRndStdDev;
            hold on
            hLines = plot(plotLocs(:,1),plotLocs(:,2),'c');
            uistack(hLines,'bottom');
            set(hLines,'hittest','off');
            hold off;

        end

        % Add the text
        [~,inds] = sort(cD,'descend');
        for i = 1:min([20,length(inds)]);
            hText(i) = text(cLocs(inds(i),1),cLocs(inds(i),2),cNames{inds(i)});
            set(hText(i),'fontsize',10,'fontweight','bold');
            set(hText,'hittest','off');
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function example;

I = load_images('house');

figure('Name','Input sequence');
subplot(2,2,1); imshow(I(:,:,:,1));
subplot(2,2,2); imshow(I(:,:,:,2));
subplot(2,2,3); imshow(I(:,:,:,3));
subplot(2,2,4); imshow(I(:,:,:,4));

R = exposure_fusion(I,[1 1 1]);
figure('Name','Result');
imshow(R);

author: mahmoudnafifi


text:
function value = monomial_value ( dim_num, point_num, x, expon )

%*****************************************************************************80
%
%% MONOMIAL_VALUE evaluates a monomial.
%
%  Discussion:
%
%    This routine evaluates a monomial of the form
%
%      product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)
%
%    where the exponents are nonnegative integers.  Note that
%    if the combination 0^0 is encountered, it should be treated
%    as 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer POINT_NUM, the number of points at which the
%    monomial is to be evaluated.
%
%    Input, real X(DIM_NUM,POINT_NUM), the point coordinates.
%
%    Input, integer EXPON(DIM_NUM), the exponents.
%
%    Output, real VALUE(POINT_NUM), the value of the monomial.
%
  value(1:point_num) = 1.0;

  for dim = 1 : dim_num
    if ( 0 ~= expon(dim) )
      value(1:point_num) = value(1:point_num) .* x(dim,1:point_num).^expon(dim);
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [controlFlux, objFlux] = robustnessAnalysis(model, controlRxn, nPoints, plotResFlag, objRxn, objType)
% Performs robustness analysis for a reaction of
% interest and an objective of interest
%
% USAGE:
%
%    [controlFlux, objFlux] = robustnessAnalysis(model, controlRxn, nPoints, plotResFlag, objRxn, objType)
%
% INPUTS:
%    model:         COBRA model structure
%    controlRxn:    Reaction of interest whose value is to be controlled
%
% OPTIONAL INPUTS:
%    nPoints:        Number of points to show on plot (Default = 20)
%    plotResFlag:    Plot results (Default true)
%    objRxn:         Objective reaction to be maximized
%                    (Default = whatever is defined in model)
%    objType:        Maximize ('max') or minimize ('min') objective
%                    (Default = 'max')
%
% OUTPUTS:
%    controlFlux:    Flux values within the range of the maximum and minimum for
%                    reaction of interest
%    objFlux:        Optimal values of objective reaction at each control
%                    reaction flux value
%
% .. Author: - Monica Mo and Markus Herrgard 8/17/06

if (nargin < 3)
    nPoints = 20;
end
if (nargin < 4)
    plotResFlag = true;
end
if (nargin > 4)
    baseModel = changeObjective(model,objRxn);
else
    baseModel = model;
end
if (nargin <6)
    objType = 'max';
end

if (findRxnIDs(model,controlRxn))
    tmpModel = changeObjective(model,controlRxn);
    tmpModel.osenseStr = 'min';
    solMin = optimizeWBModel(tmpModel);
    tmpModel.osenseStr = 'max';
    solMax = optimizeWBModel(tmpModel);
else
    error('Control reaction does not exist!');
end

objFlux = [];
controlFlux = linspace(solMin.f,solMax.f,nPoints)';

showprogress(0,'Robustness analysis in progress ...');
for i=1:length(controlFlux)
    showprogress(i/length(controlFlux));
    modelControlled = changeRxnBounds(baseModel,controlRxn,controlFlux(i),'b');
    tmpModel.osenseStr = objType;
    solControlled = optimizeWBModel(modelControlled);
    objFlux(i) = solControlled.f;
end

objFlux = objFlux';

if (plotResFlag)
    plot(controlFlux,objFlux)
    xlabel(strrep(controlRxn,'_','-'));
    ylabel('Objective');
    axis tight;
end


text:
function  sort_dcm_series(root_dir)
% Sorts dcm files into folders by copying them

curr_dir=pwd;

% addpath(genpath(fileparts(fileparts(curr_dir))));

if(nargin < 1)
    error('Give a DICOM directory');
end

chdir(root_dir);
is_dicom_dir=0;
if exist('dicom','dir')
    dicomstr = 'dicom';
elseif(exist('DICOM','dir'))
    dicomstr = 'DICOM';
else
    dicomstr = '';
end
if(~exist([root_dir,filesep,dicomstr,'_ALL'],'dir'))
    copyfile(dicomstr, [dicomstr,'_ALL']);
end
if(~isempty(dicomstr))
    chdir(dicomstr);
    is_dicom_dir=1;
    %mkdir('../dicom_sorted');
    mkdir('../../dicom_sorted');
end
if(~exist('sort_dcm_series','file')) % sort_dcm_series not on path
    dicomMfilePath = fileparts(mfilename('fullpath'));
    addpath(genpath(dicomMfilePath));
end
files=get_file_list('*.IMA');
%files=get_file_list('*');

for f=1:length(files)
    dcm=dicominfo(files{f});
    dcm_path=[ regexprep(dcm.SeriesDescription,' ','_') '_' sprintf('%04d',dcm.SeriesNumber)];
    dcm_path=regexprep(dcm_path,'[*<>/&%$=+]','');
    if is_dicom_dir
        %dcm_path=[ '../dicom_sorted/' dcm_path];
        dcm_path=[ '../../dicom_sorted/' dcm_path];
    end
    if ~exist(dcm_path,'dir')
        mkdir(dcm_path);
        disp(dcm_path);
    end %if
    if exist(dcm_path,'dir')
        mv(files{f},dcm_path);
    end

end %f

chdir(root_dir);
mv('../dicom_sorted','.');
if(length(dir('DICOM')) == 2); rmdir DICOM; end
%status = rmdir('folderName')

chdir(curr_dir);

end

author: thomaskuestner


text:
function value = year_length_hebrew ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_HEBREW returns the number of days in a Hebrew year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of
%    days in the year.
%
  jed = new_year_to_jed_hebrew ( y );

  y2 = y + 1;
  jed2 = new_year_to_jed_hebrew ( y2 );

  value = round ( jed2 - jed );

  return
end


text:
report_this_filefun(mfilename('fullpath'));

[lat,lon] = meshgrat(tmap,tmapleg);

ren = interp2(X2,Y2,Z,lon,lat);

l = ren < 0.01;
ren(l) = 0.01;
mi = min(min(ren));

l = isnan(ren);
ren(l) = mi-2;

figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[33.5 35 ],'MapLonLimit',[-117.3 -116])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',8);
tightmap
view([0 90])
% hl = lightangle(45,25);
hl = camlight
lighting phong
set(gca,'projection','perspective');

load worldlo
%h = displaym(POline); set(h(1),'color',[0.9 0.9 0.9],'Linewidth',1.7)
% h2 = displaym(PPpoint);
%   h = displaym(PPtext); trimcart(h);

plotm(faults(:,2), faults(:,1),'k','Linewidth',1.4);
% plotm(mainfault(:,2), mainfault(:,1),'m','Linewidth',3);

pl = plotm(a.Latitude,a.Longitude,'ok');
set(pl,'LineWidth',0.3,'MarkerSize',2,...
    'MarkerFaceColor','k','MarkerEdgeColor','k')
pl = plotm(main(:,2),main(:,1),'hw');
set(pl,'LineWidth',1,'MarkerSize',20,...
    'MarkerFaceColor','w','MarkerEdgeColor','k')

zdatam(handlem('allline'),10000) % keep line on surface
%zdatam(handlem('alltext'),10000) % keep line on surface
caxis([0.008 0.09])
j = jet;
%j = j(64:-1:1,:);
j = [ [ 0.85 0.9 0.9] ; j];

colormap(j); brighten(0.3);

axis off; set(gcf,'color','w')

setm(gca,'ffacecolor','w')
setm(gca,'fedgecolor','k','flinewidth',3);

setm(gca,'mlabellocation',0.5)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',0.5)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','k','Fontweight','bold','FontSize',14,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.7 0.15 0.01 0.3],'TickDir','out','Ycolor','k','Xcolor','k',...
    'Fontweight','bold','FontSize',14);
set(gcf,'Inverthardcopy','off');


text:
function [pathway, stats, sPtrNew] = mtrReadCaminoPathwayFromStream(rawPaths, sPtr, mode)

% Default mode is the overriden contrack way
if notDefined('mode'), mode = 1; end

if mode==0
    % Actual camino format
    numPoints = rawPaths(sPtr);
    seedID = rawPaths(sPtr+1);
    stats = seedID;
    pathway = [];
    %pathway = reshape(rawPaths(sPtr+2:sPtr+1+3*numPoints),3,numPoints);
    sPtrNew = sPtr+3*numPoints+2;
else
    % Lightweight Query format
    numStats = rawPaths(sPtr);
    numPoints = rawPaths(sPtr+1);
    stats = rawPaths(sPtr+2:sPtr+2+numStats-1);
    pathway = reshape(rawPaths(sPtr+numStats+2:sPtr+numStats+1+3*numPoints),3,numPoints);
    sPtrNew = sPtr+3*numPoints+numStats+2;
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function peaks = PeakDetection(x,ff,varargin)
%
% peaks = PeakDetection(x,f,flag),
% Overview: R-peak detector based on max search
%
% inputs:
% x: vector of input data
% ff: approximate ECG beat-rate in Hertz, normalized by the sampling frequency
% flag: search for positive (flag=1) or negative (flag=0) peaks. By default
% the maximum absolute value of the signal, determines the peak sign.
%
% output:
% peaks: vector of R-peak impulse train
%
% Notes:
% - The R-peaks are found from a peak search in windows of length N; where
% N corresponds to the R-peak period calculated from the given f. R-peaks
% with periods smaller than N/2 or greater than N are not detected.
% - The signal baseline wander is recommended to be removed before the
% R-peak detection
%
%
% Open Source ECG Toolbox, version 1.0, November 2006
% Released under the GNU General Public License
% Copyright (C) 2006  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- GIPSA-Lab, INPG, Grenoble, France
% reza.sameni@gmail.com

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

N = length(x);
peaks = zeros(1,N);

th = .5;
rng = floor(th/ff);

if(nargin==3)
    flag = varargin{1};
else
    flag = abs(max(x))>abs(min(x));
end

if(flag)
    for j = 1 : N
        %         index = max(j-rng,1):min(j+rng,N);
        if(j>rng && j<N-rng)
            index = j-rng:j+rng;
        elseif(j>rng)
            index = N-2*rng:N;
        else
            index = 1:2*rng;
        end

        if(max(x(index))==x(j))
            peaks(j) = 1;
        end
    end
else
    for j = 1 : N
        %         index = max(j-rng,1):min(j+rng,N);
        if(j>rng && j<N-rng)
            index = j-rng:j+rng;
        elseif(j>rng)
            index = N-2*rng:N;
        else
            index = 1:2*rng;
        end

        if(min(x(index))==x(j))
            peaks(j) = 1;
        end
    end
end

% remove fake peaks
I = find(peaks);
d = diff(I);
% z = find(d<rng);
peaks(I(d<rng))=0;

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox