text:
function [TruePositives, FalseNegatives] = testDrugMetabolism(model, microbeID, biomassReaction, database)
% Performs an FVA and reports those drug metabolites (exchange reactions)
% that can be taken up and/or secreted by the model and should be secreted according to
% data (true positives) and those bile acid metabolites that cannot be secreted by
% the model but should be secreted according to in vitro data (false
% negatives).
%
% INPUT
% model             COBRA model structure
% microbeID         Microbe ID in carbon source data file
% biomassReaction   Biomass objective functions (low flux through BOF
%                   required in analysis)
% database          Structure containing rBioNet reaction and metabolite
%                   database
%
% OUTPUT
% TruePositives     Cell array of strings listing all drug metabolites
%                   (exchange reactions) that can be taken up and/or secreted by the model
%                   and in in vitro data.
% FalseNegatives    Cell array of strings listing all drug metabolites
%                   (exchange reactions) that cannot be taken up and/or secreted by the model
%                   but should be secreted according to in vitro data.
%
% .. Author:
%       Almut Heinken, April 2020
%                      March 2022 - changed code to string-matching to make
%                      it more robust

global CBT_LP_SOLVER
if isempty(CBT_LP_SOLVER)
    initCobraToolbox
end

% read drug metabolism table
dataTable = readInputTableForPipeline('drugTable.txt');

corrRxns = {'eALP','EX_r788(e)','EX_r406(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cALP','EX_r788(e)','EX_r406(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cAzNAD','EX_bzd(e)','EX_abz_ala_b(e)','EX_5asa(e)','EX_olsa(e)','EX_neopront(e)','EX_sanilamide(e)','EX_ssz(e)','EX_sulfp(e)','EX_tab(e)','EX_pront(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cAzNADP','EX_bzd(e)','EX_abz_ala_b(e)','EX_5asa(e)','EX_olsa(e)','EX_neopront(e)','EX_sanilamide(e)','EX_ssz(e)','EX_sulfp(e)','EX_tab(e)','EX_pront(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'eBGal','EX_gal(e)','EX_fru(e)','EX_lactl(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cBGal','EX_lactl(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cBRV','EX_brv(e)','EX_srv(e)','EX_bvu(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cCda','EX_dfduri(e)','EX_dfdcytd(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'eCda','EX_dfduri(e)','EX_dfdcytd(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cCgr','EX_digoxin(e)','EX_digitoxin(e)','EX_dihydro_digitoxin(e)','EX_dihydro_digoxin(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cCodA','EX_5fura(e)','EX_fcsn(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'eCodA','EX_5fura(e)','EX_fcsn(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cDpdNAD','EX_dh5fura(e)','EX_5fura(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cDpdNADP','EX_dh5fura(e)','EX_5fura(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'eDpdNADP','EX_dh5fura(e)','EX_5fura(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cHpd','EX_4hphac(e)','EX_pcresol(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'TdcA','EX_34dhphe(e)','EX_dopa(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'Dadh','EX_mtym(e)','EX_dopa(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cNAT','EX_5asa(e)','EX_isnzd(e)','EX_acisnzd(e)','EX_ac5asa(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'eNAT','EX_5asa(e)','EX_isnzd(e)','EX_acisnzd(e)','EX_ac5asa(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'eNit','EX_nzp(e)','EX_chlphncl(e)','EX_nchlphncl(e)','EX_7a_czp(e)','EX_anzp(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'cNit','EX_nzp(e)','EX_chlphncl(e)','EX_nchlphncl(e)','EX_7a_czp(e)','EX_anzp(e)','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','','';'eUidA','EX_1hibup_S(e)','EX_1hibupglu_S(e)','EX_1hmdgluc(e)','EX_1ohmdz(e)','EX_2hatvacid(e)','EX_2hatvacidgluc(e)','EX_2hatvlac(e)','EX_2hatvlacgluc(e)','EX_2hibup_S(e)','EX_2hibupglu_S(e)','EX_2oh_cbz(e)','EX_2oh_cbz_glc(e)','EX_2oh_mtz(e)','EX_2oh_mtz_glc(e)','EX_3meacmp(e)','EX_3hibup_S(e)','EX_3hibupglu_S(e)','EX_3oh_cbz(e)','EX_3oh_cbz_glc(e)','EX_3oh_dlor(e)','EX_3oh_dlor_glc(e)','EX_3oh_mdea(e)','EX_3oh_mdea_glc(e)','EX_3oh_mxn(e)','EX_3oh_mxn_glc(e)','EX_4dh_tpno(e)','EX_4dh_tpno_1glc(e)','EX_4hmdgluc(e)','EX_4oh_dcf(e)','EX_4oh_dcf_glc(e)','EX_4oh_kp(e)','EX_4oh_kp_glc(e)','EX_4oh_levole(e)','EX_4oh_levole_glc(e)','EX_4oh_meth(e)','EX_4oh_meth_glc(e)','EX_4oh_propl(e)','EX_4oh_propl_glc(e)','EX_4oh_trz(e)','EX_4oh_trz_glc(e)','EX_4oh_vcz(e)','EX_4oh_vcz_glc(e)','EX_4ohmdz(e)','EX_5oh_sulfp(e)','EX_5oh_sulfp_glc(e)','EX_5ohfvs(e)','EX_5ohfvsglu(e)','EX_6bhglz(e)','EX_6bhglzglc(e)','EX_6ohfvs(e)','EX_6ohfvsglu(e)','EX_7bhglz(e)','EX_7bhglzglc(e)','EX_7oh_efv(e)','EX_7oh_efv_glc(e)','EX_814dioh_efv(e)','EX_814dioh_efv_glc(e)','EX_8oh_efv(e)','EX_8oh_efv_glc(e)','EX_ac_amn_b_gly(e)','EX_ac_amn_b_gly_glc(e)','EX_acmp(e)','EX_acmp_glc(e)','EX_acmpglu(e)','EX_alpz_4oh(e)','EX_alpz_4oh_glc(e)','EX_alpz_aoh(e)','EX_alpz_aoh_glc(e)','EX_am14(e)','EX_am14_glc(e)','EX_am1ccs(e)','EX_am1cglc(e)','EX_am5(e)','EX_am5_glc(e)','EX_am6(e)','EX_am6_glc(e)','EX_amio(e)','EX_amio_c(e)','EX_amio_c_glc(e)','EX_amio_glc(e)','EX_amn_b_gly(e)','EX_amn_b_gly_glc(e)','EX_amntd_m6(e)','EX_amntd_m6_glc(e)','EX_atvacid(e)','EX_atvacylgluc(e)','EX_atvethgluc(e)','EX_atvlac(e)','EX_atvlacgluc(e)','EX_bhpm(e)','EX_bhpm_glc(e)','EX_bilr_355(e)','EX_bilr_M10(e)','EX_bilr_M12(e)','EX_bilr_M14(e)','EX_bilr_M15(e)','EX_bilr_M16(e)','EX_bilr_M7(e)','EX_bsn(e)','EX_bsn_glc(e)','EX_bz(e)','EX_bz_glc(e)','EX_caribup_s(e)','EX_caribupglu_S(e)','EX_cbz(e)','EX_cbz_glc(e)','EX_cd6168(e)','EX_cd6168_glc(e)','EX_chlphncl(e)','EX_chlphncl_glc(e)','EX_clcxb(e)','EX_clcxb_c(e)','EX_clcxb_c_glc(e)','EX_clcxb_glc(e)','EX_clobi_c(e)','EX_clobi_glc(e)','EX_crvsm1(e)','EX_crvsm23(e)','EX_cvm1gluc(e)','EX_cvm23gluc(e)','EX_daa(e)','EX_daa_glc(e)','EX_dcf(e)','EX_dcf_glc(e)','EX_ddea(e)','EX_ddea_glc(e)','EX_des_astzl(e)','EX_des_astzl_glc(e)','EX_digoxin(e)','EX_digoxin_glc(e)','EX_dlb(e)','EX_dlb_glc(e)','EX_dnpz_5des(e)','EX_dnpz_6des(e)','EX_dnpz_m11(e)','EX_dnpz_m12(e)','EX_dnpz_m13(e)','EX_dnpz_m14(e)','EX_dnpz_m9(e)','EX_doh_etr(e)','EX_doh_etr_glc(e)','EX_doh_vcz(e)','EX_doh_vcz_glc(e)','EX_dxo(e)','EX_dxo_glc(e)','EX_efv(e)','EX_efv_glc(e)','EX_eltr(e)','EX_eltr_glc(e)','EX_eltr_m3(e)','EX_eltr_m4(e)','EX_eztmb(e)','EX_eztmb_glc(e)','EX_fvs(e)','EX_fvsgluc(e)','EX_fvstet(e)','EX_fvstetglu(e)','EX_glc3meacp(e)','EX_gltmn(e)','EX_gltmn_glc(e)','EX_gmfl(e)','EX_gmfl_glc(e)','EX_gmfl_mI(e)','EX_gmfl_mI_glc(e)','EX_gmfl_mII(e)','EX_gmfl_mII_glc(e)','EX_gmfl_mIII(e)','EX_gmfl_mIII_glc(e)','EX_gtacmp(e)','EX_hst(e)','EX_hst_3_glc(e)','EX_hst_3_s(e)','EX_hst_37_diglc(e)','EX_hst_3glc_7s(e)','EX_hst_7_glc(e)','EX_hst_7_s(e)','EX_hst_7glc_3s(e)','EX_ibup_S(e)','EX_ibupgluc(e)','EX_imn(e)','EX_imn_glc(e)','EX_inv(e)','EX_inv_m1(e)','EX_isosorbide_5mn(e)','EX_isosorbide_5mn_glc(e)','EX_kprofen(e)','EX_kprofen_glc(e)','EX_lst4exp(e)','EX_lstn(e)','EX_lstn1gluc(e)','EX_lstnm4(e)','EX_lstnm7(e)','EX_mdz(e)','EX_mdz_glc(e)','EX_mdzglc(e)','EX_miso(e)','EX_miso_glc(e)','EX_mrphn(e)','EX_mrphn_3glc(e)','EX_mrphn_6glc(e)','EX_mtz(e)','EX_mtz_glc(e)','EX_N_oh_phtn(e)','EX_N_oh_phtn_glc(e)','EX_nsldp_m5(e)','EX_nsldp_m5_glc(e)','EX_nverp(e)','EX_nverp_glc(e)','EX_odsm_egltmn(e)','EX_odsm_egltmn_glc(e)','EX_odsm_gltmn(e)','EX_odsm_gltmn_glc(e)','EX_oh_etr(e)','EX_oh_etr_glc(e)','EX_oh_pbl(e)','EX_oh_pbl_glc(e)','EX_phppa(e)','EX_phppa_glc(e)','EX_phtn_glc(e)','EX_prob(e)','EX_prob_glc(e)','EX_pront(e)','EX_pront_glc(e)','EX_propl(e)','EX_propl_glc(e)','EX_prx_mI(e)','EX_prx_mI_glc(e)','EX_prx_mII(e)','EX_prx_mII_glc(e)','EX_ptvst(e)','EX_ptvstgluc(e)','EX_pvs(e)','EX_pvsgluc(e)','EX_R_6oh_warf(e)','EX_R_6oh_warf_glc(e)','EX_R_7oh_warf(e)','EX_R_7oh_warf_glc(e)','EX_R_8oh_warf(e)','EX_R_8oh_warf_glc(e)','EX_r406(e)','EX_r406_glc(e)','EX_r529(e)','EX_r529_glc(e)','EX_rep(e)','EX_rep_glc(e)','EX_rpn_104557(e)','EX_rpn_104557_cb_glc(e)','EX_rpn_96990(e)','EX_rpn_96990_glc(e)','EX_rpn_oh(e)','EX_rpn_oh_glc(e)','EX_rsv(e)','EX_rsvgluc(e)','EX_S_4oh_warf(e)','EX_S_4oh_warf_glc(e)','EX_S_6oh_warf(e)','EX_S_6oh_warf_glc(e)','EX_sb_611855(e)','EX_sb_y(e)','EX_sch_488128(e)','EX_sch_57871(e)','EX_sch_57871_glc(e)','EX_sfnd_1689(e)','EX_sfnd_1689_glc(e)','EX_sftz(e)','EX_sftz_glc(e)','EX_simvgluc(e)','EX_smap(e)','EX_smap_glc(e)','EX_smvacid(e)','EX_sn38(e)','EX_sn38g(e)','EX_spz(e)','EX_spz_glc(e)','EX_spz_sfn(e)','EX_spz_sfn_glc(e)','EX_stg(e)','EX_stg_m3(e)','EX_stg_m4(e)','EX_tat(e)','EX_tgz(e)','EX_tgz_glc(e)','EX_thsacmp(e)','EX_tlf_a(e)','EX_tlf_a_1a(e)','EX_tlf_a_1b(e)','EX_tlf_a_1x(e)','EX_tlf_a_2(e)','EX_tlf_a_2a(e)','EX_tlf_a_3(e)','EX_tlf_a_4(e)','EX_tlf_a_m1(e)','EX_tlf_a_m2(e)','EX_tlf_a_m3(e)','EX_tlf_a_m4(e)','EX_tlf_a_m4a(e)','EX_tlf_a_m5(e)','EX_tlf_a_m5a(e)','EX_tlf_a_m5b(e)','EX_tlf_a_m6(e)','EX_tlf_a_m9(e)','EX_tlms(e)','EX_tlms_glc(e)','EX_tmacmp(e)','EX_tolcp(e)','EX_tolcp_ac(e)','EX_tolcp_ac_glc(e)','EX_tolcp_am(e)','EX_tolcp_am_glc(e)','EX_tolcp_glc(e)','EX_tpno_1glc_4g(e)','EX_tpno_4g(e)','EX_tpno_4glc(e)','EX_tpnoh(e)','EX_tsacmgluc(e)';'cUidA','EX_1hibup_S(e)','EX_1hibupglu_S(e)','EX_1hmdgluc(e)','EX_1ohmdz(e)','EX_2hatvacid(e)','EX_2hatvacidgluc(e)','EX_2hatvlac(e)','EX_2hatvlacgluc(e)','EX_2hibup_S(e)','EX_2hibupglu_S(e)','EX_2oh_cbz(e)','EX_2oh_cbz_glc(e)','EX_2oh_mtz(e)','EX_2oh_mtz_glc(e)','EX_3meacmp(e)','EX_3hibup_S(e)','EX_3hibupglu_S(e)','EX_3oh_cbz(e)','EX_3oh_cbz_glc(e)','EX_3oh_dlor(e)','EX_3oh_dlor_glc(e)','EX_3oh_mdea(e)','EX_3oh_mdea_glc(e)','EX_3oh_mxn(e)','EX_3oh_mxn_glc(e)','EX_4dh_tpno(e)','EX_4dh_tpno_1glc(e)','EX_4hmdgluc(e)','EX_4oh_dcf(e)','EX_4oh_dcf_glc(e)','EX_4oh_kp(e)','EX_4oh_kp_glc(e)','EX_4oh_levole(e)','EX_4oh_levole_glc(e)','EX_4oh_meth(e)','EX_4oh_meth_glc(e)','EX_4oh_propl(e)','EX_4oh_propl_glc(e)','EX_4oh_trz(e)','EX_4oh_trz_glc(e)','EX_4oh_vcz(e)','EX_4oh_vcz_glc(e)','EX_4ohmdz(e)','EX_5oh_sulfp(e)','EX_5oh_sulfp_glc(e)','EX_5ohfvs(e)','EX_5ohfvsglu(e)','EX_6bhglz(e)','EX_6bhglzglc(e)','EX_6ohfvs(e)','EX_6ohfvsglu(e)','EX_7bhglz(e)','EX_7bhglzglc(e)','EX_7oh_efv(e)','EX_7oh_efv_glc(e)','EX_814dioh_efv(e)','EX_814dioh_efv_glc(e)','EX_8oh_efv(e)','EX_8oh_efv_glc(e)','EX_ac_amn_b_gly(e)','EX_ac_amn_b_gly_glc(e)','EX_acmp(e)','EX_acmp_glc(e)','EX_acmpglu(e)','EX_alpz_4oh(e)','EX_alpz_4oh_glc(e)','EX_alpz_aoh(e)','EX_alpz_aoh_glc(e)','EX_am14(e)','EX_am14_glc(e)','EX_am1ccs(e)','EX_am1cglc(e)','EX_am5(e)','EX_am5_glc(e)','EX_am6(e)','EX_am6_glc(e)','EX_amio(e)','EX_amio_c(e)','EX_amio_c_glc(e)','EX_amio_glc(e)','EX_amn_b_gly(e)','EX_amn_b_gly_glc(e)','EX_amntd_m6(e)','EX_amntd_m6_glc(e)','EX_atvacid(e)','EX_atvacylgluc(e)','EX_atvethgluc(e)','EX_atvlac(e)','EX_atvlacgluc(e)','EX_bhpm(e)','EX_bhpm_glc(e)','EX_bilr_355(e)','EX_bilr_M10(e)','EX_bilr_M12(e)','EX_bilr_M14(e)','EX_bilr_M15(e)','EX_bilr_M16(e)','EX_bilr_M7(e)','EX_bsn(e)','EX_bsn_glc(e)','EX_bz(e)','EX_bz_glc(e)','EX_caribup_s(e)','EX_caribupglu_S(e)','EX_cbz(e)','EX_cbz_glc(e)','EX_cd6168(e)','EX_cd6168_glc(e)','EX_chlphncl(e)','EX_chlphncl_glc(e)','EX_clcxb(e)','EX_clcxb_c(e)','EX_clcxb_c_glc(e)','EX_clcxb_glc(e)','EX_clobi_c(e)','EX_clobi_glc(e)','EX_crvsm1(e)','EX_crvsm23(e)','EX_cvm1gluc(e)','EX_cvm23gluc(e)','EX_daa(e)','EX_daa_glc(e)','EX_dcf(e)','EX_dcf_glc(e)','EX_ddea(e)','EX_ddea_glc(e)','EX_des_astzl(e)','EX_des_astzl_glc(e)','EX_digoxin(e)','EX_digoxin_glc(e)','EX_dlb(e)','EX_dlb_glc(e)','EX_dnpz_5des(e)','EX_dnpz_6des(e)','EX_dnpz_m11(e)','EX_dnpz_m12(e)','EX_dnpz_m13(e)','EX_dnpz_m14(e)','EX_dnpz_m9(e)','EX_doh_etr(e)','EX_doh_etr_glc(e)','EX_doh_vcz(e)','EX_doh_vcz_glc(e)','EX_dxo(e)','EX_dxo_glc(e)','EX_efv(e)','EX_efv_glc(e)','EX_eltr(e)','EX_eltr_glc(e)','EX_eltr_m3(e)','EX_eltr_m4(e)','EX_eztmb(e)','EX_eztmb_glc(e)','EX_fvs(e)','EX_fvsgluc(e)','EX_fvstet(e)','EX_fvstetglu(e)','EX_glc3meacp(e)','EX_gltmn(e)','EX_gltmn_glc(e)','EX_gmfl(e)','EX_gmfl_glc(e)','EX_gmfl_mI(e)','EX_gmfl_mI_glc(e)','EX_gmfl_mII(e)','EX_gmfl_mII_glc(e)','EX_gmfl_mIII(e)','EX_gmfl_mIII_glc(e)','EX_gtacmp(e)','EX_hst(e)','EX_hst_3_glc(e)','EX_hst_3_s(e)','EX_hst_37_diglc(e)','EX_hst_3glc_7s(e)','EX_hst_7_glc(e)','EX_hst_7_s(e)','EX_hst_7glc_3s(e)','EX_ibup_S(e)','EX_ibupgluc(e)','EX_imn(e)','EX_imn_glc(e)','EX_inv(e)','EX_inv_m1(e)','EX_isosorbide_5mn(e)','EX_isosorbide_5mn_glc(e)','EX_kprofen(e)','EX_kprofen_glc(e)','EX_lst4exp(e)','EX_lstn(e)','EX_lstn1gluc(e)','EX_lstnm4(e)','EX_lstnm7(e)','EX_mdz(e)','EX_mdz_glc(e)','EX_mdzglc(e)','EX_miso(e)','EX_miso_glc(e)','EX_mrphn(e)','EX_mrphn_3glc(e)','EX_mrphn_6glc(e)','EX_mtz(e)','EX_mtz_glc(e)','EX_N_oh_phtn(e)','EX_N_oh_phtn_glc(e)','EX_nsldp_m5(e)','EX_nsldp_m5_glc(e)','EX_nverp(e)','EX_nverp_glc(e)','EX_odsm_egltmn(e)','EX_odsm_egltmn_glc(e)','EX_odsm_gltmn(e)','EX_odsm_gltmn_glc(e)','EX_oh_etr(e)','EX_oh_etr_glc(e)','EX_oh_pbl(e)','EX_oh_pbl_glc(e)','EX_phppa(e)','EX_phppa_glc(e)','EX_phtn_glc(e)','EX_prob(e)','EX_prob_glc(e)','EX_pront(e)','EX_pront_glc(e)','EX_propl(e)','EX_propl_glc(e)','EX_prx_mI(e)','EX_prx_mI_glc(e)','EX_prx_mII(e)','EX_prx_mII_glc(e)','EX_ptvst(e)','EX_ptvstgluc(e)','EX_pvs(e)','EX_pvsgluc(e)','EX_R_6oh_warf(e)','EX_R_6oh_warf_glc(e)','EX_R_7oh_warf(e)','EX_R_7oh_warf_glc(e)','EX_R_8oh_warf(e)','EX_R_8oh_warf_glc(e)','EX_r406(e)','EX_r406_glc(e)','EX_r529(e)','EX_r529_glc(e)','EX_rep(e)','EX_rep_glc(e)','EX_rpn_104557(e)','EX_rpn_104557_cb_glc(e)','EX_rpn_96990(e)','EX_rpn_96990_glc(e)','EX_rpn_oh(e)','EX_rpn_oh_glc(e)','EX_rsv(e)','EX_rsvgluc(e)','EX_S_4oh_warf(e)','EX_S_4oh_warf_glc(e)','EX_S_6oh_warf(e)','EX_S_6oh_warf_glc(e)','EX_sb_611855(e)','EX_sb_y(e)','EX_sch_488128(e)','EX_sch_57871(e)','EX_sch_57871_glc(e)','EX_sfnd_1689(e)','EX_sfnd_1689_glc(e)','EX_sftz(e)','EX_sftz_glc(e)','EX_simvgluc(e)','EX_smap(e)','EX_smap_glc(e)','EX_smvacid(e)','EX_sn38(e)','EX_sn38g(e)','EX_spz(e)','EX_spz_glc(e)','EX_spz_sfn(e)','EX_spz_sfn_glc(e)','EX_stg(e)','EX_stg_m3(e)','EX_stg_m4(e)','EX_tat(e)','EX_tgz(e)','EX_tgz_glc(e)','EX_thsacmp(e)','EX_tlf_a(e)','EX_tlf_a_1a(e)','EX_tlf_a_1b(e)','EX_tlf_a_1x(e)','EX_tlf_a_2(e)','EX_tlf_a_2a(e)','EX_tlf_a_3(e)','EX_tlf_a_4(e)','EX_tlf_a_m1(e)','EX_tlf_a_m2(e)','EX_tlf_a_m3(e)','EX_tlf_a_m4(e)','EX_tlf_a_m4a(e)','EX_tlf_a_m5(e)','EX_tlf_a_m5a(e)','EX_tlf_a_m5b(e)','EX_tlf_a_m6(e)','EX_tlf_a_m9(e)','EX_tlms(e)','EX_tlms_glc(e)','EX_tmacmp(e)','EX_tolcp(e)','EX_tolcp_ac(e)','EX_tolcp_ac_glc(e)','EX_tolcp_am(e)','EX_tolcp_am_glc(e)','EX_tolcp_glc(e)','EX_tpno_1glc_4g(e)','EX_tpno_4g(e)','EX_tpno_4glc(e)','EX_tpnoh(e)','EX_tsacmgluc(e)'};

TruePositives = {};
FalseNegatives = {};
% find microbe index in drug metabolism table
mInd = find(strcmp(dataTable(:,1), microbeID));
if isempty(mInd)
    warning(['Microbe "', microbeID, '" not found in drug data file.'])
else
    % perform FVA
    % set BOF
    if ~any(ismember(model.rxns, biomassReaction)) || nargin < 3
        error(['Biomass reaction "', biomassReaction, '" not found in model.'])
    end
    model = changeObjective(model, biomassReaction);
    % set a low lower bound for biomass
    %     model = changeRxnBounds(model, biomassReaction, 1e-3, 'l');
    % list exchange reactions
    exchanges = model.rxns(strncmp('EX_', model.rxns, 3));
    % open all exchanges
    model = changeRxnBounds(model, exchanges, -1000, 'l');
    model = changeRxnBounds(model, exchanges, 1000, 'u');

    % get the reactions to test
    rxns = {};
    for i=2:size(dataTable,2)
        if contains(version,'(R202') % for Matlab R2020a and newer
            if dataTable{mInd,i}==1
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
                rxns = union(rxns,corrRxns(findCorrRxns,2:end));
            end
        else
            if strcmp(dataTable{mInd,i},'1')
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
                rxns = union(rxns,corrRxns(findCorrRxns,2:end));
            end
        end
    end

    % flux variability analysis on reactions of interest
    rxns = unique(rxns);
    rxns = rxns(~cellfun('isempty', rxns));
    rxnsInModel=intersect(rxns,model.rxns);
    if ~isempty(rxnsInModel)
        currentDir=pwd;
        try
            [minFlux, maxFlux, ~, ~] = fastFVA(model, 0, 'max', 'ibm_cplex', ...
                rxnsInModel, 'S');
        catch
            warning('fastFVA could not run, so fluxVariability is instead used. Consider installing fastFVA for shorter computation times.');
            cd(currentDir)
            [minFlux, maxFlux] = fluxVariability(model, 0, 'max', rxnsInModel);
        end

        % active flux
        flux = rxnsInModel(minFlux < -1e-6);
        flux = union(flux,rxnsInModel(maxFlux > 1e-6));
    else
        flux = {};
    end
    % check all exchanges corresponding to drug uptake/secretion
    % in this case, all of them should be carrying flux
    for i=2:size(dataTable,2)
        findCorrRxns = [];
        if contains(version,'(R202') % for Matlab R2020a and newer
            if dataTable{mInd,i}==1
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
            end
        else
            if strcmp(dataTable{mInd,i},'1')
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
            end
        end
        if ~isempty(findCorrRxns)
            allEx = corrRxns(findCorrRxns,2:end);
            allEx = allEx(~cellfun(@isempty, allEx));

            TruePositives = union(TruePositives, intersect(allEx, flux));
            FalseNegatives = union(FalseNegatives, setdiff(allEx, flux));
        end
    end
end

% replace reaction IDs with metabolite names
if ~isempty(TruePositives)
    TruePositives = TruePositives(~cellfun(@isempty, TruePositives));
    TruePositives=strrep(TruePositives,'EX_','');
    TruePositives=strrep(TruePositives,'(e)','');

    for i=1:length(TruePositives)
        TruePositives{i}=database.metabolites{find(strcmp(database.metabolites(:,1),TruePositives{i})),2};
    end
end

% warn about false negatives
if ~isempty(FalseNegatives)
    FalseNegatives = FalseNegatives(~cellfun(@isempty, FalseNegatives));
    FalseNegatives=strrep(FalseNegatives,'EX_','');
    FalseNegatives=strrep(FalseNegatives,'(e)','');
    for i = 1:length(FalseNegatives)
        FalseNegatives{i}=database.metabolites{find(strcmp(database.metabolites(:,1),FalseNegatives{i})),2};
    end
end

end

author: opencobra

repo: cobratoolbox


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: distances and SSD forces for hand data
%
%   - data                 Hand, Omega=(0,20)x(0,25), level=3:7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             MI
%   - transformation       translation2D
%==============================================================================

clear, close all, help(mfilename);

fprintf('setup data, viewer, interpolation, transformation\n');
setup2DhandData; level = 6; omega = ML{level}.omega; m = ML{level}.m;
viewImage('reset',viewPara{:},'axis','off');
imgModel('reset','imgModel','splineInter');
trafo('reset','trafo','translation2D');
R  = imgModel('coefficients',ML{level}.R,[],omega,'out',0);
xc = reshape(getCellCenteredGrid(omega,m),[],2);
yc = reshape(trafo([0.75;0],xc(:)),[],2);

fprintf('%20s : %s\n','viewImage',viewImage);
fprintf('%20s : %s\n','imgModel',imgModel);

% compute R(xc), T(yc), dT and forces
Rc = imgModel(R,omega,xc);
[Tc,dT] = imgModel(R,omega,yc);
dT = full(spdiags(dT,[0,size(dT,1)]));
F  = spdiags(Tc-Rc,0,length(Tc),length(Tc))*dT;

% compute lengthes for plots
fac = 128/norm(dT(:),'inf');
lengthdT = sqrt(sum(dT.^2,2));
normdT   = max(lengthdT);
J  = find(lengthdT>1e-2*normdT);
lengthF = sum(F.^2,2);
normF   = sqrt(max(lengthF));
K  = find(lengthF>1e-2*normF);

% plot R, T, T-R
FAIRfigure(1,'figname',mfilename); clf;
subplot(1,3,1);  viewImage(Rc,omega,m);            th(1) = title('R');
subplot(1,3,2);  viewImage(Tc,omega,m);            th(2) = title('T');
subplot(1,3,3);  viewImage(128+(Tc-Rc)/2,omega,m); th(3) = title('T-R');
set(th,'fontsize',30);
FAIRpause;

% plot \partial_j T and \nabla T
figure(2); clf;  colordef(gcf,'black');
subplot(1,3,1);  viewImage(128+fac*dT(:,1),omega,m);    th(1) = title('\partial_1T');
subplot(1,3,2);  viewImage(128+fac*dT(:,2),omega,m);    th(2) = title('\partial_2T');
subplot(1,3,3);  viewImage(128+Tc/2,omega,m); hold on;
qh = quiver(xc(J,1),xc(J,2),dT(J,1)/normdT,dT(J,2)/normdT,2);
set(qh,'color','b','linewidth',1.5);
th(3) = title('\nabla T');
set(th,'fontsize',30);
FAIRpause;

% plot \nabla T, T-R, forces
figure(3); clf;  colordef(gcf,'black');
subplot(1,3,1);  viewImage(128+Tc/2,omega,m); hold on;
qh = quiver(xc(J,1),xc(J,2),dT(J,1)/normdT,dT(J,2)/normdT,0);
set(qh,'color','b','linewidth',1.5);
th(1) = title('\nabla T');

subplot(1,3,2);  vh = viewImage(128+(Tc-Rc)/2,omega,m); th(2) = title('T-R');

subplot(1,3,3);  viewImage(128+Tc/2,omega,m); hold on;
qh = quiver(xc(K,1),xc(K,2),F(K,1)/normF,F(K,2)/normF,2);
set(qh,'color','r','linewidth',1.5);
th(3) = title('forces');
set(th','fontsize',30);

%==============================================================================


text:
function varargout = circleAsPolygon(circle, varargin)
%CIRCLEASPOLYGON Convert a circle into a series of points
%
%   Note: this function is deprecated, use "circleToPolygon" instead
%
%   P = circleAsPolygon(CIRCLE, N);
%   convert circle given as [x0 y0 r], where x0 and y0 are coordinate of
%   center, and r is the radius, into an array of  [(N+1)x2] double,
%   containing x and y values of points.
%   The polygon is closed
%
%   P = circleAsPolygon(CIRCLE);
%   uses a default value of N=64 points
%
%   Example
%   circle = circleAsPolygon([10 0 5], 16);
%   figure;
%   drawPolygon(circle);
%
%   See also:
%   circles2d, circleToPolygon
%
%
% ---------
% author : David Legland
% created the 06/04/2005.
% Copyright 2010 INRA - Cepia Software Platform.
%

%   HISTORY
%   20/04/2007: return a closed polygon with N+1 vertices, use default N=64

warning('matGeom:deprecated', ...
    'function "circleAsPolygon" is deprecated, use "circleToPolygon" instead');

% format output
if nargout <= 1
    varargout = {circleToPolygon(circle, varargin{:})};
else
    [x, y] = circleToPolygon(circle, varargin{:});
    varargout = {x, y};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function varargout = process_stdrow( varargin )
% PROCESS_STDROW: Uniformize the list of rows for a set of time-frequency files.
%
% USAGE: OutputFiles = process_stdrow('Run', sProcess, sInputs)
%        [DestRowNames, AllRowNames, iRowsSrc, iRowsDest, msgError] = process_stdrow('GetUniformRows', FileNames, Method)
%        [DestRowNames, AllRowNames, iRowsSrc, iRowsDest, msgError] = process_stdrow('GetUniformRows', LoadedFiles, Method)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2013-2016

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Uniform signal names';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Standardize';
    sProcess.Index       = 301;
    sProcess.Description = '';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'timefreq', 'matrix'};
    sProcess.OutputTypes = {'timefreq', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 2;
    % Definition of the options
    % === TARGET LIST OF ROW NAMES
    sProcess.options.method.Comment = {'Keep only the common row names
=> Remove all the others', ...
                                       'Keep all the row names
=> Fill the missing signals with zeros', ...
                                       'Use the first file in the list as a template'};
    sProcess.options.method.Type    = 'radio';
    sProcess.options.method.Value   = 1;
    % === OVERWRITE
    sProcess.options.overwrite.Comment = 'Overwrite input files';
    sProcess.options.overwrite.Type    = 'checkbox';
    sProcess.options.overwrite.Value   = 0;
    sProcess.options.overwrite.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    switch (sProcess.options.method.Value)
        case 1,    Comment = [sProcess.Comment, ' (remove extra)'];
        case 2,    Comment = [sProcess.Comment, ' (add missing)'];
        case 3,    Comment = [sProcess.Comment, ' (use first)'];
    end
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Options
    switch (sProcess.options.method.Value)
        case 1,  Method = 'common';
        case 2,  Method = 'all';
        case 3,  Method = 'first';
    end
    isOverwrite = sProcess.options.overwrite.Value;
    OutputFiles = {};

    % ===== GET TARGET LIST OF SIGNALS =====
    % Get list of common signals
    [DestRowNames, AllRowNames, iRowsSrc, iRowsDest, msgError] = GetUniformRows({sInputs.FileName}, Method);
    % Report errors and warnings
    if ~isempty(msgError)
        if isempty(DestRowNames)
            bst_report('Error', sProcess, sInputs, msgError);
            return;
        else
            bst_report('Warning', sProcess, sInputs, msgError);
        end
    end
    % Check if there are any difference in the row names
    if all(cellfun(@(c)isequal(c,DestRowNames), AllRowNames))
        bst_report('Error', sProcess, sInputs, 'All the input files have identical row names.');
        return;
    end

    % ===== PROCESS FILES =====
    % Process each input file
    for iInput = 1:length(sInputs)
        % If it's a file that was not changed: skip to next file
        if isequal(DestRowNames, AllRowNames{iInput})
            OutputFiles{iInput} = file_fullpath(sInputs(iInput).FileName);
            continue;
        end

        % List of added rows
        iAddedRows = setdiff(1:length(DestRowNames), iRowsDest{iInput});
        iRemRows   = setdiff(1:length(AllRowNames{iInput}), iRowsSrc{iInput});
        % Add a history entry
        strHistory = 'Uniform list of signals:';
        if ~isempty(iAddedRows)
            strTmp = '';
            for i = 1:length(iAddedRows)
                strTmp = [strTmp, DestRowNames{iAddedRows(i)}, ','];
            end
            strHistory = [strHistory, sprintf(' %d added (%s)', length(iAddedRows), strTmp(1:end-1))];
        end
        if ~isempty(iRemRows)
            strTmp = '';
            for i = 1:length(iRemRows)
                strTmp = [strTmp, AllRowNames{iInput}{iRemRows(i)}, ','];
            end
            strHistory = [strHistory, sprintf(' %d removed (%s)', length(iRemRows), strTmp(1:end-1))];
        end

        % Load the data file
        fileMat = load(file_fullpath(sInputs(iInput).FileName));
        newFileMat = fileMat;
        % Saved structure depends on the file type
        switch lower(sInputs(iInput).FileType)
            case 'timefreq'
                newFileMat.TF = zeros(length(DestRowNames), size(fileMat.TF,2), size(fileMat.TF,3));
                newFileMat.TF(iRowsDest{iInput},:) = fileMat.TF(iRowsSrc{iInput},:);
                newFileMat.RowNames = DestRowNames;
                % Do not keep the Std field in the output
                if isfield(newFileMat, 'Std') && ~isempty(newFileMat.Std)
                    newFileMat.Std = [];
                end
            case 'matrix'
                newFileMat.Value = zeros(length(DestRowNames), size(fileMat.Value,2));
                newFileMat.Value(iRowsDest{iInput},:) = fileMat.Value(iRowsSrc{iInput},:);
                newFileMat.Description = DestRowNames;
        end
        % Add comment
        newFileMat.Comment = [newFileMat.Comment ' | stdrow'];
        % Add a history entry
        newFileMat = bst_history('add', newFileMat, 'stdrow', strHistory);

        % Overwrite the input file
        if isOverwrite
            OutputFiles{iInput} = file_fullpath(sInputs(iInput).FileName);
            bst_save(OutputFiles{iInput}, newFileMat, 'v6');
        % Create a new file
        else
            % Output filename: add file tag
            OutputFiles{iInput} = strrep(file_fullpath(sInputs(iInput).FileName), '.mat', '_stdrow.mat');
            OutputFiles{iInput} = file_unique(OutputFiles{iInput});
            % Save file
            bst_save(OutputFiles{iInput}, newFileMat, 'v6');
            % Add file to database structure
            db_add_data(sInputs(iInput).iStudy, OutputFiles{iInput}, newFileMat);
        end
    end
end

%% ===== GET UNIFORM ROWS =====
function [DestRowNames, AllRowNames, iRowsSrc, iRowsDest, msgError] = GetUniformRows(FileNames, Method)
    AllRowNames   = cell(1,length(FileNames));
    nRows         = zeros(1,length(FileNames));
    iRowsSrc      = cell(1,length(FileNames));
    iRowsDest     = cell(1,length(FileNames));
    unionRowNames = {};
    interRowNames = {};
    DestRowNames  = [];
    msgError      = [];
    % Check all the input files
    for iInput = 1:length(FileNames)
        % Load row names
        switch (file_gettype(FileNames{iInput}))
            case 'timefreq'
                if ischar(FileNames{iInput})
                    fileMat = in_bst_timefreq(FileNames{iInput}, 0, 'DataType', 'RowNames', 'RefRowNames');
                else
                    fileMat = FileNames{iInput};
                end
                % Check file type: Cannot process source files
                if strcmpi(fileMat.DataType, 'results') || ~iscell(fileMat.RowNames)
                    msgError = 'STDROW> Cannot process source maps, or any file that does not have explicit row names.';
                    return;
                elseif (length(fileMat.RefRowNames) > 1) && ~isequal(fileMat.RowNames, fileMat.RefRowNames)
                    msgError = 'STDROW> Cannot process connectivity [AxB] results.';
                    return;
                end
                % Add row names to the list
                AllRowNames{iInput} = fileMat.RowNames;
            case 'matrix'
                if ischar(FileNames{iInput})
                    fileMat = in_bst_matrix(FileNames{iInput}, 'Description');
                else
                    fileMat = FileNames{iInput};
                end
                % Check file type
                if (size(fileMat.Description,2) > 1)
                    msgError = 'Cannot process a matrix file in which the "Description" fields has more than one column.';
                    return;
                end
                % Add row names to the list
                AllRowNames{iInput} = fileMat.Description;
            otherwise
                error('Unsupported file format.');
        end

        % Remove the @filename at the end of the row names
        for iRow = 1:length(AllRowNames{iInput})
            iAt = find(AllRowNames{iInput}{iRow} == '@', 1);
            if ~isempty(iAt) && any(AllRowNames{iInput}{iRow}(iAt+1:end) == '/')
                AllRowNames{iInput}{iRow} = strtrim(AllRowNames{iInput}{iRow}(1:iAt-1));
            end
        end

        % Keep track of row numbers
        nRows(iInput) = length(AllRowNames{iInput});
        % Union of all the row names
        unionRowNames = union(unionRowNames, AllRowNames{iInput});
        % Intersection of all the row names
        if isempty(interRowNames)
            interRowNames = AllRowNames{iInput};
        else
            interRowNames = intersect(interRowNames, AllRowNames{iInput});
        end
    end
    % Check if there are rowns left
    if isempty(interRowNames) && strcmpi(Method, 'common')
        msgError = 'No common row names in these data sets.';
        return;
    end

    % ===== COMMON ROW LIST =====
    % Get the row list that has the more/less rows
    switch (Method)
        % Only common rows
        case 'common'
            % Get the minimum number of rows
            [tmp, iRef] = min(nRows);
            % Get rows
            DestRowNames = AllRowNames{iRef};
            % Remove unecessary rows
            iRemove = find(~ismember(DestRowNames, interRowNames));
            if ~isempty(iRemove)
                DestRowNames(iRemove) = [];
            end
        % All rows
        case 'all'
            % Get the maximum number of rows
            [tmp, iRef] = max(nRows);
            % Get rows
            DestRowNames = AllRowNames{iRef};
            % Add all the other rows
            iAdd = find(~ismember(unionRowNames, DestRowNames));
            if ~isempty(iAdd)
                if (size(DestRowNames,1) == 1)
                    DestRowNames = [DestRowNames, unionRowNames{iAdd}];
                else
                    DestRowNames = [DestRowNames; unionRowNames{iAdd}];
                end
            end
        % First file
        case 'first'
            DestRowNames = AllRowNames{1};
    end

    % ===== GET INDICES FOR EACH FILE =====
    % Process each input file
    for iInput = 1:length(FileNames)
        % If it's a file that was not changed: skip to next file
        if isequal(DestRowNames, AllRowNames{iInput})
            iRowsSrc{iInput}  = 1:length(DestRowNames);
            iRowsDest{iInput} = 1:length(DestRowNames);
            continue;
        end
        % Create list of orders for rows
        for iChan = 1:length(DestRowNames)
            iTmp = find(strcmpi(DestRowNames{iChan}, AllRowNames{iInput}));
            iTmp = setdiff(iTmp, iRowsSrc{iInput});
            if (length(iTmp) > 1)
                if (length(iTmp) > 1)
                    msgError = 'Several signals with the same name, re-ordering might be inaccurate.';
                    iTmp = iTmp(1);
                end
            end
            if ~isempty(iTmp)
                iRowsDest{iInput}(end+1) = iChan;
                iRowsSrc{iInput}(end+1)  = iTmp;
            end
        end
    end
end

author: brainstorm-tools


text:
function [ a, b ] = p08_ab ( m )

%*****************************************************************************80
%
%% P08_AB evaluates the limits of the optimization region for problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Marcin Molga, Czeslaw Smutnicki,
%    Test functions for optimization needs.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Output, real A(M), B(M), the lower and upper bounds.
%
  a(1:m) = - 1.0;
  b(1:m) = + 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Copyright (C) 2016, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym airy (@var{k}, @var{x})
%% @defmethodx @@sym airy (@var{x})
%% Symbolic Airy functions of first/second kind and their derivatives.
%%
%% @var{k} can be 0, 1, 2, or 3; as in the documentation for the
%% non-symbolic Airy function, @pxref{airy}.
%%
%% Examples:
%% @example
%% @group
%% syms z
%% @c doctest: +SKIP_UNLESS(pycall_sympy__ ('return Version(spver) > Version("1.4")'))
%% Ai = airy(0, z)
%%   @result{} Ai = (sym) Ai(z)
%%
%% @c doctest: +SKIP_UNLESS(pycall_sympy__ ('return Version(spver) > Version("1.4")'))
%% Bi = airy(2, z)
%%   @result{} Bi = (sym) Bi(z)
%%
%% @c doctest: +SKIP_UNLESS(pycall_sympy__ ('return Version(spver) > Version("1.4")'))
%% Bi_prime = airy(3, z)
%%   @result{} Bi_prime = (sym) Bi'(z)
%%
%% @c doctest: +SKIP_UNLESS(pycall_sympy__ ('return Version(spver) > Version("1.4")'))
%% diff(Bi, z)
%%   @result{} (sym) Bi'(z)
%%
%% @c doctest: +SKIP_UNLESS(pycall_sympy__ ('return Version(spver) > Version("1.4")'))
%% diff(Bi, z, z)
%%   @result{} (sym) z⋅Bi(z)
%% @end group
%% @end example
%%
%% @seealso{@@sym/besselj, @@sym/bessely, @@sym/besseli, @@sym/besselk, @@sym/besselh}
%% @end defmethod

function A = airy(k, x)

  if (nargin == 2)
    % no-op
  elseif (nargin == 1)
    x = k;
    k = 0;
  else
    print_usage ();
  end

  assert(isscalar(k))

  if (logical(k == 0))
    A = elementwise_op ('airyai', sym(x));
  elseif (logical(k == 1))
    A = elementwise_op ('airyaiprime', sym(x));
  elseif (logical(k == 2))
    A = elementwise_op ('airybi', sym(x));
  elseif (logical(k == 3))
    A = elementwise_op ('airybiprime', sym(x));
  else
    error('airy: expecting K = 0, 1, 2, or 3')
  end

end

%!test
%! syms z
%! a = airy(0, z);
%! ap = airy(1, z);
%! assert (isequal (diff (a), ap))
%! assert (isequal (diff (ap), z*a))

%!test
%! syms z
%! b = airy(2, z);
%! bp = airy(3, z);
%! assert (isequal (diff (b), bp))
%! assert (isequal (diff (bp), z*b))

%!test
%! % default to k=0
%! syms z
%! a = airy(0, z);
%! a2 = airy(z);
%! assert (isequal (a, a2))

%!error airy(0, sym('x'), 2)
%!error <expecting K = 0, 1, 2, or 3> airy(4, sym('z'))
%!error <expecting K = 0, 1, 2, or 3> airy(-1, sym('z'))

%!test
%! % symbolic k
%! syms z
%! b1 = airy(2, z);
%! b2 = airy(sym(2), z);
%! assert (isequal (b1, b2))

%!test
%! % doubles, relative error
%! X = [1 2 pi; 4i 5 6+6i];
%! Xs = sym(X);
%! for k = 0:3
%!   A = double(airy(k, Xs));
%!   B = airy(k, X);
%!   assert (all (all (abs(A - B) < 500*eps*abs(A))))
%! end

%!test
%! % round-trip
%! syms x
%! for k = 0:3
%!   A = airy(k, 10);
%!   q = airy(k, x);
%!   h = function_handle(q);
%!   B = h(10);
%!   assert (abs(A-B) < 500*eps*abs(A))
%! end

repo: octsympy


text:
function vl_xmkdir(path, varargin)
% VL_XMKDIR  Create a directory recursively.
%   VL_XMKDIR(PATH) creates all directory specified by PATH if they
%   do not exist (existing directories are skipped).
%
%   The function accepts the following options:
%
%   Pretend:: [false]
%     Set to true to avoid creating any directory but print which
%     directories would be created (implies 'Verbose',true).
%
%   Verbose:: [false]
%     Set to true to print the operations performed.
%
%   See also: VL_HELP().

% Authors: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

pretend = false ;
verbose = false ;

for ai=1:2:length(varargin)
  opt = varargin{ai} ;
  arg = varargin{ai+1} ;
  switch lower(opt)
    case 'pretend'
      pretend = arg ;
    case 'verbose'
      verbose = arg ;
    otherwise
      error('Unknown option ''%s''.', opt) ;
  end
end

verbose = verbose | pretend ;

if isempty(path)
  return ;
end

[subPath, name, ext] = fileparts(path) ;
name = [name ext] ;

% subPath is equal to path if path = '/'
if ~strcmp(subPath, path)
  vl_xmkdir(subPath, varargin{:}) ;
end

% name is empty if path = '*/'
if isempty(name)
  return ;
end

if verbose
  fprintf('%s: directory ''%s'' in ''%s''', mfilename, name, subPath) ;
end

if ~exist(path, 'dir')
  if pretend && verbose
    fprintf(' does not exist, would create.\n') ;
  else
    mkdir(subPath, name) ;
    if verbose
      fprintf(' did not exist, created.\n') ;
    end
  end
else
  if verbose
    if pretend
      fprintf(' exists, would not create.\n') ;
    else
      fprintf(' exists, not creating.\n') ;
    end
  end
end

author: jianxiongxiao

repo: ProfXkit


text:
classdef prtClassRvmForContext < prtClass
    % prtClassRvmForContext  Relevance vector machine classifier for
    % context-dependent applications

    properties (SetAccess=private)
        name = 'Relevance Vector Machine'  % Relevance Vector Machine
        nameAbbreviation = 'RVM'           % RVM
        isNativeMary = false;  % False
        useForContext = true;
    end

    properties
        kernels = prtKernelDc & prtKernelDirect;  % The kernels to be used

        verboseText = true;  % Whether or not to display text during training
        verbosePlot = false;  % Whether or not to plot during training

        maxIterations = 1000;       % The maximum number of iterations
        convergenceThreshold = .01;  % Learning tolerance;

        prior = struct('a',1e-6,'b',1e-6','c',1e-6,'d',1e-6);
        weights = [];

    end

    methods

        function Obj = prtClassRvm(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end

        function Obj = set.maxIterations(Obj,val)
            if ~prtUtilIsPositiveInteger(val)
                error('prt:prtClassRvm:maxIterations','maxIterations must be a positive integer');
            end
            Obj.learningMaxIterations = val;
        end

        function Obj = set.convergenceThreshold(Obj,val)
            if ~prtUtilIsPositiveScalar(val)
                error('prt:prtClassRvm:convergenceThreshold','convergenceThreshold must be a positive scalar');
            end
            Obj.learningConvergedTolerance = val;
        end

        function Obj = set.kernels(Obj,val)
            assert(numel(val)==1 &&  isa(val,'prtKernel'),'prt:prtClassRvm:kernels','kernels must be a prtKernel');

            Obj.kernels = val;
        end

        function Obj = set.verbosePlot(Obj,val)
            assert(isscalar(val) && (islogical(val) || prtUtilIsPositiveInteger(val)),'prt:prtClassRvm:verbosePlot','verbosePlot must be a logical value or a positive integer');
            Obj.verbosePlot = val;
        end

        function Obj = set.verboseText(Obj,val)
            assert(isscalar(val) && islogical(val),'prt:prtClassRvm:verboseText','verboseText must be a logical value, but value provided is a %s',class(val));
            Obj.verboseText = val;
        end

        function varargout = plot(Obj)
            % plot - Plot output confidence of the prtClassRvm object
            %
            %   CLASS.plot plots the output confidence of the prtClassRvm
            %   object. The dimensionality of the dataset must be 3 or
            %   less, and verboseStorage must be true.

            HandleStructure = plot@prtClass(Obj);

            holdState = get(gca,'nextPlot');
            hold on;
            Obj.sparseKernels.plot;
            set(gca, 'nextPlot', holdState);

            varargout = {};
            if nargout > 0
                varargout = {HandleStructure};
            end
        end

    end

    methods (Access=protected, Hidden = true)

        function Obj = trainAction(Obj,DataSet)
            %Note: do not assume that getTargets returns a double array or
            %values "0" and "1", instead use this:
            Y = Obj.getMinusOneOneTargets(DataSet);
            Y(Y==-1) = 0;

            Obj.kernels = Obj.kernels.train(DataSet);
            GRAMM = Obj.kernels.run_OutputDoubleArray(DataSet);

            prior = Obj.prior;
            prior.A = diag(repmat(prior.a,1,size(GRAMM,2)));
            prior.S = eye(size(GRAMM,2));
            prior.m = zeros(size(GRAMM,2),1);

            if isfield(DataSet.observationInfo,'contextPosterior')
                for i = 1:DataSet.nObservations
                    pCgivenX(i,:) = DataSet.observationInfo(i).contextPosterior;
                end
            else
                pCgivenX = ones(DataSet.nObservations,1);
            end
            pCgivenXmat = repmat(pCgivenX',size(GRAMM,2),1);

            % Run VB
            iter = 0;
            finished = false;
            while ~finished
                iter = iter + 1;

                phi = GRAMM';
                N = size(phi,2);              % Number of samples
                D = size(phi,1);              % Dimensionality
                t = Y;

                if iter == 1 %% Initialization
                    a = repmat(prior.a,1,D);
                    b = repmat(prior.b,1,D);
                    S = prior.S;
                    m = prior.m;
                end

                xi = sqrt(diag(phi'*(S+m*m')*phi));
                lambda = (1./(4*xi)).*tanh(xi./2);

                A = diag(a./b);
                S = inv(A + 2*pCgivenXmat.*phi*diag(lambda)*phi');

                m = .5*S*sum(pCgivenXmat.*repmat(2*t'-1,D,1).*phi,2);

                a = repmat(prior.a,D,1) + .5;
                b = diag(diag(repmat(prior.b,1,D)) + .5*(diag(m.^2) + diag(diag(S))));

                % Calculate Negative Free Energy (Convergence Criterion)
                sig = (1 + exp(-xi)).^-1;
                logF = sum(pCgivenX.*(log(sig) + .5*(2*t-1).*(phi'*m) - .5*xi - lambda.*((diag(phi'*(S+m*m')*phi)) - xi.^2)));

                KLDalpha = 0;
                for d = 1:D
                    KLDalpha = KLDalpha + prtRvUtilGammaKld(a(d),b(d),prior.a,prior.b);
                end
                if det(S) ~= 0
                    KLDw  = .5*(-sum(log(diag(A))) -log(det(S)) + trace(A.*S)' + sum(diag(A).*m.^2) - N);
                else
                    KLDw  = .5*(-sum(log(diag(A))) -log(realmin) + trace(A.*S)' + sum(diag(A).*m.^2) - N);
                end

                NFE(iter) = logF - KLDalpha - KLDw;
                if Obj.verbosePlot
                    figure(666),plot(NFE)
                end

                if iter > 1
                    if NFE(iter) > 0
                        Lpct = 100*(NFE(iter) - NFE(iter-1))/NFE(iter-1);
                    else
                        Lpct = 100*(NFE(iter-1) - NFE(iter))/NFE(iter-1);
                    end
                else
                    Lpct = nan;
                end

                %% Progress mssages, check convergence
                if Obj.verboseText
                    fprintf(['Iteration #',num2str(iter),': Negative Free Energy = ',num2str(NFE(iter)),' (',num2str(Lpct),'%%)\n'])
                end

                if iter > 1
                    if abs(Lpct) <= Obj.convergenceThreshold;
                        finished = true;
                        Obj.weights = m;
                        if Obj.verboseText
                            fprintf('NFE Converged! Congratulation.\n')
                        end
                    elseif iter == Obj.maxIterations;
                        finished = true;
                        Obj.weights = m;
                        if Obj.verboseText
                            fprintf('Max iterations reached. Get out of here!\n')
                        end
                    elseif isinf(Lpct)
                        finished = true;
                        Obj.weights = m;
                        if Obj.verboseText
                            fprintf('Inf found...exiting\n')
                        end
                    end
                end
            end

        end

        function DataSetOut = runAction(Obj,DataSet)

            GRAMM = Obj.kernels.run_OutputDoubleArray(DataSet);
            x = DataSet.X;
            w = Obj.weights;
            for n = 1:size(x,1);
                y(n,:) = w'*GRAMM(n,:)';
            end
            Yout = (1 + exp(-y)).^(-1);

            DataSetOut = prtDataSetClass(Yout,DataSet.Y);
            if isa(DataSet,'prtDataSetClassContext')
                DataSetOut = prtDataSetClassContext(DataSetOut,prtDataSetClass);
            end
        end
    end

    methods (Access=protected, Hidden = true)

        function y = getMinusOneOneTargets(Obj, DataSet) %#ok<MANU>
            yMat = double(DataSet.getTargetsAsBinaryMatrix());
            y = nan(size(yMat,1),1);
            y(yMat(:,1) == 1) = -1;
            y(yMat(:,2) == 1) = 1;
        end

        function G = regularizeGramInnerProduct(Obj, gram)
            nBasis = size(gram,2);

            sigmaSquared = 1e-6;

            %Check to make sure the problem is well-posed.  This can be fixed either
            %with changes to kernels, or by regularization
            G = gram'*gram;
            while rcond(G) < 1e-6
                if sigmaSquared == eps && Obj.verboseText
                    %warning('prt:prtClassRvm:illConditionedG','RVM initial G matrix ill-conditioned; regularizing diagonal of G to resolve; this can be modified by changing kernel parameters\n');
                    fprintf('\n\tRegularizing Gram matrix...\n');
                end
                G = (sigmaSquared*eye(nBasis) + gram'*gram);
                sigmaSquared = sigmaSquared*2;
            end

        end

        function verboseIterationPlot(Obj,DataSet,relevantIndices)
            DsSummary = DataSet.summarize;

            [linGrid, gridSize,xx,yy] = prtPlotUtilGenerateGrid(DsSummary.lowerBounds, DsSummary.upperBounds, Obj.plotOptions); %#ok<ASGLU>

            localKernels = Obj.kernels.train(DataSet);
            cKernels = localKernels.retainKernelDimensions(relevantIndices);
            cPhiDataSet = cKernels.run(prtDataSetClass([xx(:),yy(:)]));
            cPhi = cPhiDataSet.getObservations;

            confMap = reshape(prtRvUtilNormCdf(cPhi*Obj.beta(relevantIndices)),gridSize);
            imagesc(xx(1,:),yy(:,1),confMap,[0,1])
            colormap(Obj.plotOptions.twoClassColorMapFunction());
            axis xy
            hold on
            plot(DataSet);
            cKernels.plot();
            hold off;

            set(gcf,'color',[1 1 1]);
            drawnow;
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function p = abss(p);
%ABSS         Polynomial with absolute values of coefficients
%
%   r = abss(p);
%
% Obsolete, replaced by mag (thanks to Arnold for better naming).
%

% written  08/28/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 10/18/08     S.M. Rump  obsolete
%

  p = mag(p);

repo: OpenMAS


text:
function [a] = spm_uncat(x,a)
% converts a matrix into an array
% FORMAT [a] = spm_uncat(x,a)
% x - matrix
% a - cell array
%
% i.e. a = spm_uncat(spm_cat(a),a)
%
% see also spm_vec and spm_unvec
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_uncat.m 5219 2013-01-29 17:07:07Z spm $

% fill in varargout
%--------------------------------------------------------------------------
[p,q] = size(a);
for i = 1:p
    for j = 1:q
          m(i,j) = size(a{i,j},1);
          n(i,j) = size(a{i,j},2);
    end
end
m     = max(m',[],1);
n     = max(n ,[],1);
for i = find(m)
    y     = x(1:m(i),:);
    x     = x((m(i) + 1):end,:);
    for j = 1:q
        a{i,j} = y(:,1:n(j));
        y      = y(:,(n(j) + 1):end);
    end
end

author: spm


text:
classdef SelectiveSearchSegmentation < handle
    %SELECTIVESEARCHSEGMENTATION  Selective search segmentation algorithm
    %
    % The class implements the algorithm described in [uijlings2013selective].
    %
    % ## References
    % [uijlings2013selective]:
    % > Jasper RR Uijlings, Koen EA van de Sande, Theo Gevers, and Arnold WM
    % > Smeulders. "Selective search for object recognition". International
    % > journal of computer vision, 104(2):154-171, 2013.
    %
    % See also: cv.SelectiveSearchSegmentation.process
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    %% Constructor/destructor
    methods
        function this = SelectiveSearchSegmentation()
            %SELECTIVESEARCHSEGMENTATION  Constructor
            %
            %     obj = cv.SelectiveSearchSegmentation()
            %
            % See also: cv.SelectiveSearchSegmentation.process
            %
            this.id = SelectiveSearchSegmentation_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.SelectiveSearchSegmentation
            %
            if isempty(this.id), return; end
            SelectiveSearchSegmentation_(this.id, 'delete');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.SelectiveSearchSegmentation.empty,
            %  cv.SelectiveSearchSegmentation.load
            %
            SelectiveSearchSegmentation_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.SelectiveSearchSegmentation.clear,
            %  cv.SelectiveSearchSegmentation.load
            %
            b = SelectiveSearchSegmentation_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.SelectiveSearchSegmentation.load
            %
            SelectiveSearchSegmentation_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.SelectiveSearchSegmentation.save
            %
            SelectiveSearchSegmentation_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.SelectiveSearchSegmentation.save,
            %  cv.SelectiveSearchSegmentation.load
            %
            name = SelectiveSearchSegmentation_(this.id, 'getDefaultName');
        end
    end

    %% SelectiveSearchSegmentation
    methods
        function setBaseImage(this, img)
            %SETBASEIMAGE  Set a image used by switch functions to initialize the class
            %
            %     obj.setBaseImage(img)
            %
            % ## Input
            % * __img__ The image.
            %
            % See also: cv.SelectiveSearchSegmentation.process
            %
            SelectiveSearchSegmentation_(this.id, 'setBaseImage', img);
        end

        function switchToSingleStrategy(this, varargin)
            %SWITCHTOSINGLESTRATEGY  Initialize the class with the 'Single strategy' parameters
            %
            %     obj.switchToSingleStrategy()
            %     obj.switchToSingleStrategy('OptionName',optionValue, ...)
            %
            % ## Options
            % * __K__ The k parameter for the graph segmentation. default 200
            % * __Sigma__ The sigma parameter for the graph segmentation.
            %   default 0.8
            %
            % As described in [uijlings2013selective].
            %
            % See also: cv.SelectiveSearchSegmentation.switchToSelectiveSearchFast
            %
            SelectiveSearchSegmentation_(this.id, 'switchToSingleStrategy', varargin{:});
        end

        function switchToSelectiveSearchFast(this, varargin)
            %SWITCHTOSELECTIVESEARCHFAST  Initialize the class with the 'Selective search fast' parameters
            %
            %     obj.switchToSelectiveSearchFast()
            %     obj.switchToSelectiveSearchFast('OptionName',optionValue, ...)
            %
            % ## Options
            % * __BaseK__ The k parameter for the first graph segmentation.
            %   default 150
            % * __IncK__ The increment of the k parameter for all graph
            %   segmentations. default 150
            % * __Sigma__ The sigma parameter for the graph segmentation.
            %   default 0.8
            %
            % As described in [uijlings2013selective].
            %
            % See also: cv.SelectiveSearchSegmentation.switchToSelectiveSearchQuality
            %
            SelectiveSearchSegmentation_(this.id, 'switchToSelectiveSearchFast', varargin{:});
        end

        function switchToSelectiveSearchQuality(this, varargin)
            %SWITCHTOSELECTIVESEARCHQUALITY  Initialize the class with the 'Selective search fast' parameters
            %
            %     obj.switchToSelectiveSearchQuality()
            %     obj.switchToSelectiveSearchQuality('OptionName',optionValue, ...)
            %
            % ## Options
            % * __BaseK__ The k parameter for the first graph segmentation.
            %   default 150
            % * __IncK__ The increment of the k parameter for all graph
            %   segmentations. default 150
            % * __Sigma__ The sigma parameter for the graph segmentation.
            %   default 0.8
            %
            % As described in [uijlings2013selective].
            %
            % See also: cv.SelectiveSearchSegmentation.switchToSingleStrategy
            %
            SelectiveSearchSegmentation_(this.id, 'switchToSelectiveSearchQuality', varargin{:});
        end

        function addImage(this, img)
            %ADDIMAGE  Add a new image in the list of images to process
            %
            %     obj.addImage(img)
            %
            % ## Input
            % * __img__ The image.
            %
            % See also: cv.SelectiveSearchSegmentation.clearImages
            %
            SelectiveSearchSegmentation_(this.id, 'addImage', img);
        end

        function clearImages(this)
            %CLEARIMAGES  Clear the list of images to process
            %
            %     obj.clearImages()
            %
            % See also: cv.SelectiveSearchSegmentation.addImage
            %
            SelectiveSearchSegmentation_(this.id, 'clearImages');
        end

        function addGraphSegmentation(this, varargin)
            %ADDGRAPHSEGMENTATION  Add a new graph segmentation in the list of graph segementations to process
            %
            %     obj.addGraphSegmentation()
            %     obj.addGraphSegmentation('OptionName',optionValue, ...)
            %
            % ## Options
            % * __Sigma__ The sigma parameter, used to smooth image.
            %   default 0.5
            % * __K__ The k parameter of the algorithm. default 300
            % * __MinSize__ The minimum size of segments. default 100
            %
            % See also: cv.SelectiveSearchSegmentation.clearGraphSegmentations,
            %  cv.GraphSegmentation
            %
            SelectiveSearchSegmentation_(this.id, 'addGraphSegmentation', varargin{:});
        end

        function clearGraphSegmentations(this)
            %CLEARGRAPHSEGMENTATIONS  Clear the list of graph segmentations to process
            %
            %     obj.clearGraphSegmentations()
            %
            % See also: cv.SelectiveSearchSegmentation.addGraphSegmentation,
            %  cv.GraphSegmentation
            %
            SelectiveSearchSegmentation_(this.id, 'clearGraphSegmentations');
        end

        function addStrategy(this, stype, varargin)
            %ADDSTRATEGY  Add a new strategy in the list of strategy to process
            %
            %     obj.addStrategy(stype)
            %     obj.addStrategy('Multiple', stype, stype, ...)
            %
            % ## Input
            % * __stype__ The strategy type for the selective search
            %   segmentation algorithm, one of:
            %   * __Color__ Color-based strategy.
            %   * __Size__ Size-based strategy.
            %   * __Texture__ Texture-based strategy.
            %   * __Fill__ Fill-based strategy.
            %   * __Multiple__ Regroup multiple strategies, where all
            %     sub-strategies have equal weights.
            %
            % The classes are implemented from the algorithm described in
            % [uijlings2013selective].
            %
            % See also: cv.SelectiveSearchSegmentation.clearStrategies
            %
            SelectiveSearchSegmentation_(this.id, 'addStrategy', stype, varargin{:});
        end

        function clearStrategies(this)
            %CLEARSTRATEGIES  Clear the list of strategy to process
            %
            %     obj.clearStrategies()
            %
            % See also: cv.SelectiveSearchSegmentation.addStrategy
            %
            SelectiveSearchSegmentation_(this.id, 'clearStrategies');
        end

        function rects = process(this)
            %PROCESSIMAGE  Based on all images, graph segmentations and stragies, computes all possible rects and return them
            %
            %     rects = obj.process()
            %
            % ## Output
            % * __rects__ The list of rects as a Nx4 numeric matrix
            %   `[x,y,w,h; ...]`. The first ones are more relevents than the
            %   lasts ones.
            %
            % See also: cv.SelectiveSearchSegmentation.setBaseImage
            %
            rects = SelectiveSearchSegmentation_(this.id, 'process');
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% GEOM3D Geometry 3D Toolbox
% Version 1.22 06-Jun-2018 .
%
%   Creation, transformations, algorithms and visualization of geometrical
%   3D primitives, such as points, lines, planes, polyhedra, circles and
%   spheres.
%
%   Euler Angles are defined as follow:
%   PHI is the azimut, i.e. the angle of the projection on horizontal plane
%   with the Ox axis, with value beween 0 and 180 degrees.
%   THETA is the latitude, i.e. the angle with the Oz axis, with value
%   between -90 and +90 degrees.
%   PSI is the 'roll', i.e. the rotation around the (PHI, THETA) direction,
%   with value in degrees
%   See also the 'angles3d' page.
%
%   Base format for primitives:
%   Point:      [x0 y0 z0]
%   Vector:     [dx dy dz]
%   Line:       [x0 y0 z0 dx dy dz]
%   Edge:       [x1 y1 z1 x2 y2 z2]
%   Plane:      [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
%   Sphere:     [x0 y0 z0 R]
%   Circle:     [x0 y0 z0 R PHI THETA PSI] (origin+center+normal+'roll').
%   Ellipsoid:  [x0 y0 z0 A B C PHI THETA PSI]
%   Cylinder:   [X1 Y1 Z1 X2 Y2 Z2 R]
%   Box:        [xmin xmax ymin ymax zmin zmax]. Used for clipping shapes.
%
%   Polygons are represented by N-by-3 array of points, the last point is
%   not necessarily the same as the first one. Points must be coplanar.
%
%
% 3D Points
%   points3d                    - Description of functions operating on 3D points.
%   midPoint3d                  - Middle point of two 3D points or of a 3D edge.
%   isCoplanar                  - Tests input points for coplanarity in 3-space.
%   transformPoint3d            - Transform a point with a 3D affine transform.
%   distancePoints3d            - Compute euclidean distance between pairs of 3D Points.
%   clipPoints3d                - Clip a set of points by a box or other 3d shapes.
%   drawPoint3d                 - Draw 3D point on the current axis.
%
% 3D Vectors
%   vectors3d                   - Description of functions operating on 3D vectors.
%   transformVector3d           - Transform a vector with a 3D affine transform.
%   normalizeVector3d           - Normalize a 3D vector to have norm equal to 1.
%   vectorNorm3d                - Norm of a 3D vector or of set of 3D vectors.
%   hypot3                      - Diagonal length of a cuboidal 3D box .
%   crossProduct3d              - Vector cross product faster than inbuilt MATLAB cross.
%   vectorAngle3d               - Angle between two 3D vectors.
%   isParallel3d                - Check parallelism of two 3D vectors.
%   isPerpendicular3d           - Check orthogonality of two 3D vectors.
%   drawVector3d                - Draw vector at a given position.
%
% Angles
%   angles3d                    - Conventions for manipulating angles in 3D.
%   anglePoints3d               - Compute angle between three 3D points.
%   sphericalAngle              - Compute angle between points on the sphere.
%   angleSort3d                 - Sort 3D coplanar points according to their angles in plane.
%   randomAngle3d               - Return a 3D angle uniformly distributed on unit sphere.
%
% Coordinate transforms
%   sph2cart2                   - Convert spherical coordinates to cartesian coordinates.
%   cart2sph2                   - Convert cartesian coordinates to spherical coordinates.
%   cart2sph2d                  - Convert cartesian coordinates to spherical coordinates in degrees.
%   sph2cart2d                  - Convert spherical coordinates to cartesian coordinates in degrees.
%   cart2cyl                    - Convert cartesian to cylindrical coordinates.
%   cyl2cart                    - Convert cylindrical to cartesian coordinates.
%
% 3D Lines and Edges
%   lines3d                     - Description of functions operating on 3D lines.
%   edges3d                     - Description of functions operating on 3D edges.
%   createLine3d                - Create a line with various inputs.
%   createEdge3d                - Create an edge between two 3D points, or from a 3D line.
%   fitLine3d                   - Fit a 3D line to a set of points.
%   parallelLine3d              - Create 3D line parallel to another one.
%   projPointOnLine3d           - Project a 3D point orthogonally onto a 3D line.
%   distancePointLine3d         - Euclidean distance between 3D point and line.
%   isPointOnLine3d             - Test if a 3D point belongs to a 3D line.
%   distancePointEdge3d         - Minimum distance between a 3D point and a 3D edge.
%   linePosition3d              - Return the position of a 3D point projected on a 3D line.
%   distanceLines3d             - Minimal distance between two 3D lines.
%   transformLine3d             - Transform a 3D line with a 3D affine transform.
%   reverseLine3d               - Return same 3D line but with opposite orientation.
%   midPoint3d                  - Middle point of two 3D points or of a 3D edge.
%   edgeLength3d                - Return the length of a 3D edge.
%   clipEdge3d                  - Clip a 3D edge with a cuboid box.
%   lineToEdge3d                - Convert a 3D straight line to a 3D finite edge.
%   edgeToLine3d                - Convert a 3D edge to a 3D straight line.
%   clipLine3d                  - Clip a line with a box and return an edge.
%   drawEdge3d                  - Draw 3D edge in the current axes.
%   drawLine3d                  - Draw a 3D line clipped by the current axes.
%
% Planes
%   planes3d                    - Description of functions operating on 3D planes.
%   createPlane                 - Create a plane in parametrized form.
%   fitPlane                    - Fit a 3D plane to a set of points.
%   normalizePlane              - Normalize parametric representation of a plane.
%   parallelPlane               - Parallel to a plane through a point or at a given distance.
%   reversePlane                - Return same 3D plane but with opposite orientation.
%   isPlane                     - Check if input is a plane.
%   transformPlane3d            - Transform a 3D plane with a 3D affine transform.
%   planesBisector              - Bisector plane between two other planes.
%   projPointOnPlane            - Return the orthogonal projection of a point on a plane.
%   intersectPlanes             - Return intersection line between 2 planes in space.
%   intersectThreePlanes        - Return intersection point between 3 planes in space.
%   intersectLinePlane          - Intersection point between a 3D line and a plane.
%   intersectEdgePlane          - Return intersection point between a plane and a edge.
%   distancePointPlane          - Signed distance betwen 3D point and plane.
%   projLineOnPlane             - Return the orthogonal projection of a line on a plane.
%   isBelowPlane                - Test whether a point is below or above a plane.
%   medianPlane                 - Create a plane in the middle of 2 points.
%   planeNormal                 - Compute the normal to a plane.
%   planePosition               - Compute position of a point on a plane.
%   planePoint                  - Compute 3D position of a point in a plane.
%   dihedralAngle               - Compute dihedral angle between 2 planes.
%   drawPlane3d                 - Draw a plane clipped by the current axes.
%
% 3D Polygons and curves
%   polygons3d                  - Description of functions operating on 3D polygons.
%   polygonCentroid3d           - Centroid (or center of mass) of a polygon.
%   polygonArea3d               - Area of a 3D polygon.
%   polygon3dNormalAngle        - Normal angle at a vertex of the 3D polygon.
%   intersectLinePolygon3d      - Intersection point of a 3D line and a 3D polygon.
%   intersectRayPolygon3d       - Intersection point of a 3D ray and a 3D polygon.
%   clipConvexPolygon3dHP       - Clip a convex 3D polygon with Half-space.
%   drawPolygon3d               - Draw a 3D polygon specified by a list of vertex coords.
%   drawPolyline3d              - Draw a 3D polyline specified by a list of vertex coords.
%   fillPolygon3d               - Fill a 3D polygon specified by a list of vertex coords.
%
% 3D Triangles
%   triangleArea3d              - Area of a 3D triangle.
%   distancePointTriangle3d     - Minimum distance between a 3D point and a 3D triangle.
%   intersectLineTriangle3d     - Intersection point of a 3D line and a 3D triangle.
%
% 3D circles and ellipses
%   circles3d                   - Description of functions operating on 3D circles.
%   fitCircle3d                 - Fit a 3D circle to a set of points.
%   circle3dPosition            - Return the angular position of a point on a 3D circle.
%   circle3dPoint               - Coordinates of a point on a 3D circle from its position.
%   circle3dOrigin              - Return the first point of a 3D circle.
%   drawCircle3d                - Draw a 3D circle.
%   drawCircleArc3d             - Draw a 3D circle arc.
%   drawEllipse3d               - Draw a 3D ellipse.
%
% Spheres
%   spheres                     - Description of functions operating on 3D spheres.
%   createSphere                - Create a sphere containing 4 points.
%   intersectLineSphere         - Return intersection points between a line and a sphere.
%   intersectPlaneSphere        - Return intersection circle between a plane and a sphere.
%   drawSphere                  - Draw a sphere as a mesh.
%   drawSphericalEdge           - Draw an edge on the surface of a sphere.
%   drawSphericalTriangle       - Draw a triangle on a sphere.
%   fillSphericalTriangle       - Fill a triangle on a sphere.
%   drawSphericalPolygon        - Draw a spherical polygon.
%   fillSphericalPolygon        - Fill a spherical polygon.
%   sphericalVoronoiDomain      - Compute a spherical voronoi domain.
%
% Smooth surfaces
%   equivalentEllipsoid         - Equivalent ellipsoid of a set of 3D points.
%   fitEllipse3d                - Fit an ellipse to a set of points.
%   ellipsoidSurfaceArea        - Approximated surface area of an ellipsoid.
%   oblateSurfaceArea           - Approximated surface area of an oblate ellipsoid.
%   prolateSurfaceArea          - Approximated surface area of a prolate ellipsoid.
%   cylinderSurfaceArea         - Surface area of a cylinder.
%   intersectLineCylinder       - Compute intersection points between a line and a cylinder.
%   revolutionSurface           - Create a surface of revolution from a planar curve.
%   surfaceCurvature            - Curvature on a surface from angle and principal curvatures.
%   drawEllipsoid               - Draw a 3D ellipsoid.
%   drawTorus                   - Draw a torus (3D ring).
%   drawCylinder                - Draw a cylinder.
%   drawEllipseCylinder         - Draw a cylinder with ellipse cross-section.
%   drawSurfPatch               - Draw a 3D surface patch, with 2 parametrized surfaces.
%
% Bounding boxes management
%   boxes3d                     - Description of functions operating on 3D boxes.
%   boundingBox3d               - Bounding box of a set of 3D points.
%   orientedBox3d               - Object-oriented bounding box of a set of 3D points.
%   intersectBoxes3d            - Intersection of two 3D bounding boxes.
%   mergeBoxes3d                - Merge 3D boxes, by computing their greatest extent.
%   box3dVolume                 - Volume of a 3-dimensional box.
%   randomPointInBox3d          - Generate random point(s) within a 3D box.
%   drawBox3d                   - Draw a 3D box defined by coordinate extents.
%
% Geometric transforms
%   transforms3d                - Conventions for manipulating 3D affine transforms.
%   fitAffineTransform3d        - Fit an affine transform using two point sets.
%   registerPoints3dAffine      - Fit 3D affine transform using iterative algorithm.
%   createTranslation3d         - Create the 4x4 matrix of a 3D translation.
%   createScaling3d             - Create the 4x4 matrix of a 3D scaling.
%   createRotationOx            - Create the 4x4 matrix of a 3D rotation around x-axis.
%   createRotationOy            - Create the 4x4 matrix of a 3D rotation around y-axis.
%   createRotationOz            - Create the 4x4 matrix of a 3D rotation around z-axis.
%   createBasisTransform3d      - Compute matrix for transforming a basis into another basis.
%   eulerAnglesToRotation3d     - Convert 3D Euler angles to 3D rotation matrix.
%   isTransform3d               - Check if input is a affine transformation matrix.
%   rotation3dToEulerAngles     - Extract Euler angles from a rotation matrix.
%   createRotation3dLineAngle   - Create rotation around a line by an angle theta.
%   rotation3dAxisAndAngle      - Determine axis and angle of a 3D rotation matrix.
%   createRotationVector3d      - Calculates the rotation between two vectors.
%   createRotationVectorPoint3d - Calculates the rotation between two vectors.
%   recenterTransform3d         - Change the fixed point of an affine 3D transform.
%   composeTransforms3d         - Concatenate several space transformations.
%
% Various drawing Functions
%   drawGrid3d                  - Draw a 3D grid on the current axis.
%   drawAxis3d                  - Draw a coordinate system and an origin.
%   drawAxisCube                - Draw a colored cube representing axis orientation.
%   drawCube                    - Draw a 3D centered cube, eventually rotated.
%   drawCuboid                  - Draw a 3D cuboid, eventually rotated.
%   drawPlatform                - Draw a rectangular platform with a given size.
%   drawLabels3d                - Draw text labels at specified 3D positions.
%
%
%   Credits:
%   * Several functions contributed by Sven Holcombe
%   * function isCoplanar was originally written by Brett Shoelson.
%   * Songbai Ji enhanced file intersectPlaneLine (6/23/2006).
%   * several functions contributed by oqilipo
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2005-11-07
% Homepage: http://github.com/mattools/matGeom
% http://www.pfl-cepia.inra.fr/index.php?page=geom3d
% Copyright 2005 INRA

% In development:
%   clipPolygon3dHP             - clip a 3D polygon with Half-space.
%   drawPartialPatch            - draw surface patch, with 2 parametrized surfaces.

% Deprecated:
%   vectorCross3d               - Vector cross product faster than inbuilt MATLAB cross.
%   inertiaEllipsoid            - Inertia ellipsoid of a set of 3D points.

% Others

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Checks whether the graph is directed, using the matrix transpose function.
%
% INPUTS: adjacency matrix, nxn
% OUTPUTS: boolean variable, 0 or 1
%
% Note: one-liner alternative: S=not(isSymmetric(adj));
% GB: last updated, Sep 23, 2012

function S=isDirected(adj)

S = true;
if adj==transpose(adj); S = false; end

author: aeolianine

repo: octave-networks-toolbox


text:
function res = all(a,dim)
%ANY          Like Matlab function "all" for Taylor
%
%Call
%
%   L = all(A)
%   L = all(A,dim)
%
%Same functionality as Matlab/all for intval quantity A
%

% written  05/21/09     S.M. Rump
%

  if nargin==1
    res = any(reshape(all(a.t),a.size));
  else
    res = any(reshape(all(a.t),a.size),dim);
  end