text:
function cvx_setup( varargin )

% CVX_SETUP   Sets up and tests the cvx distribution.
%    This function is to be called any time CVX is installed on a new machine,
%    to insure that the paths are set properly and the MEX files are compiled.

global cvx___
try

    cvx___ = [];
    squares = {}; %#ok
    nret = false;
    oldpath = '';
    line = '---------------------------------------------------------------------------';

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Get version and portability information %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    cvx_version( '-install', varargin{:} );
    if ~isfield( cvx___, 'loaded' ) || ~cvx___.loaded, %#ok
        error( 'CVX:Expected', 'Error detected by cvx_version' );
    end
    isoctave = cvx___.isoctave;
    mpath = cvx___.where;
    fs = cvx___.fs;
    ps = cvx___.ps;

    %%%%%%%%%%%%%%%%%%%%%%%
    % Reset the CVX paths %
    %%%%%%%%%%%%%%%%%%%%%%%

    oldpath = cvx_startup( false );
    if ~isempty( oldpath ),
        fprintf( 'Saving updated path...' );
        nret = true;
        s = warning('off'); %#ok
        stat = savepath;
        warning(s);
        if stat,
            fprintf('failed. (see below)\n');
        else
            fprintf('done.\n');
            oldpath = [];
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%
    % Search for solvers %
    %%%%%%%%%%%%%%%%%%%%%%

    try
        selected = cvx___.solvers.names{cvx___.solvers.map.default};
    catch
        selected = 'sdpt3';
    end
    cvx___.solvers = struct( 'selected', 0, 'active', 0, 'list', [], 'names', {{}}, 'map', struct( 'default', 0 ) );
    fprintf( 'Searching for solvers...' ); nret = true;
    shimpath = [ mpath, fs, 'shims', fs ];
    solvers = dir( shimpath );
    solvers = { solvers(~[solvers.isdir]).name };
    if isoctave, efilt = '\.m$'; else efilt = '\.(m|p)$'; end
    solvers = solvers( ~cellfun( @isempty, regexp( solvers, efilt ) ) );
    solvers = unique( cellfun( @(x)x(1:end-2), solvers, 'UniformOutput', false ) );
    sconfig =  struct;
    solvers = struct( ...
        'name', '', 'version', '', 'location', '', 'fullpath', '', ...
        'error', '', 'warning', '', 'path', '', 'solve', [], ...
        'settings', struct, 'sname', solvers, 'spath', shimpath, ...
        'params', struct, 'config', sconfig, 'eargs', {{}} );
    solver2 = which( 'cvx_solver_shim', '-all' );
    if ~isempty(solver2) && ~iscell(solver2),
        solver2 = { solver2 };
    end
    for k = 1 : length(solver2),
        tsolv = solver2{k};
        ndxs = find(tsolv==fs,1,'last');
        solvers(end+1).spath = tsolv(1:ndxs); %#ok
        solvers(end).sname = tsolv(ndxs+1:end-2);
    end
    cur_d = pwd;
    nsolvers = [];
    nshims = length(solvers);
    for k = 1 : nshims,
        tsolv = solvers(k);
        try
            cd(tsolv.spath);
            tsolv = feval( tsolv.sname, tsolv );
            nsolv = length(tsolv);
            if nsolv == 1,
                if ~isempty(tsolv.error), tsolv.solve = []; end
            else
               sndx = 0;
               for e = 0 : 1,
                    for j = 1 : nsolv,
                        if e ~= ( isempty( tsolv(j).error ) && ~isempty( tsolv(j).solve ) ),
                            if e, tsolv(j).solve = []; end
                            sndx = sndx + 1;
                            if sndx > 1,
                                tsolv(j).name = sprintf( '%s_%d', tsolv(j).name, sndx );
                            end
                        end
                    end
                end
            end
        catch errmsg
            errmsg = cvx_error( errmsg, 63, false, '  ' );
            if isempty( tsolv.name ),
                tsolv.name = [ tsolv.spath, tsolv.sname ];
            end
            tsolv.error = sprintf( 'unexpected error:\n%s', errmsg );
            tsolv.solve = [];
        end
        try
            nsolvers = [ nsolvers, tsolv ]; %#ok
        catch
            tsolv = solvers(k);
            tsolv.error = 'This solver shim is not compatible with CVX 3.0. Please contact the authors for an update.';
            nsolvers = [ nsolvers, tsolv ]; %#ok
        end
    end
    solvers = nsolvers;
    cd( cur_d );

    %%%%%%%%%%%%%%%%%%%%%%%%%
    % Process solver errors %
    %%%%%%%%%%%%%%%%%%%%%%%%%

    plurals = { 's', '', 's' };
    nsolv   = length(solvers);
    srej    = cellfun( @isempty, { solvers.solve } );
    sgood   = ~srej;
    fprintf( '%d shim%s found.\n', nshims, plurals{min(nshims+1,3)} );
    lens = [ 0, 0 ];
    for k = 1 : nsolv,
        lens = max( lens, [ length(solvers(k).name), length(solvers(k).version) ] );
    end
    fmt = sprintf( ' %%s  %%-%ds   %%-%ds    %%s\\n', lens(1), lens(2) );
    if any( sgood ),
        ngood = nnz(sgood);
        fprintf( '%d solver%s initialized (* = default):\n', ngood, plurals{min(ngood+1,3)} );
        stats =  { ' ', '*' };
        for k = find(sgood),
            fprintf( fmt, stats{1+strcmpi(solvers(k).name,selected)}, solvers(k).name, solvers(k).version, solvers(k).location );
            if ~isempty(solvers(k).warning),
                cvx_error( [ 'WARNING: ', solvers(k).warning ], 63, false, '        ' );
            end
        end
    else
        fprintf( [ ...
            'No valid solvers were found. This suggests a corrupt installation. Please\n', ...
            'try re-installing the files and re-running cvx_setup. If the same error\n', ...
            'occurs, please contact CVX support.\n' ] );
        error('CVX:Unexpected','No valid solvers were found.');
    end
    if any( srej ),
        t1 = srej & ( cellfun( @(x)isempty(x)||strncmp(x,'http',4), {solvers.error} ) );
        if any(t1),
            nrej = nnz(t1);
            fprintf( '%d solver%s not found:\n', nrej, plurals{min(nrej+1,3)} );
            for k = find(t1),
                err = solvers(k).error;
                if isempty(err), err = ''; end
                fprintf( fmt, ' ', solvers(k).name, err, '' ); %#ok
            end
        end
        t2 = cellfun( @(x)isequal(x,'A CVX Professional license is required.'), {solvers.error} );
        if any(t2),
            nrej = nnz(t2);
            fprintf( '%d solver%s require%s a CVX Professional license:\n', nrej, plurals{min(nrej+1,3)}, plurals{max(4-nrej,2)} );
            for k = find(t2),
                fprintf( fmt, ' ', solvers(k).name, solvers(k).version, solvers(k).location ); %#ok
            end
        end
        t3 = srej & ( ~t1 & ~t2 );
        if any(t3),
            nrej = nnz(t3);
            fprintf( '%d solver%s skipped due to errors:\n', nrej, plurals{min(nrej+1,3)} );
            for k = find(t3),
                fprintf( fmt, ' ', solvers(k).name, solvers(k).version, solvers(k).location ); %#ok
                cvx_error( solvers(k).error, 63, false, '        ' );
            end
        end
    end
    solvers = solvers(sgood);
    cvx___.solvers.list  = solvers;
    cvx___.solvers.names = { solvers.name };
    cvx___.solvers.map   = struct( 'default', 0 );
    for k = 1 : length(solvers),
        cvx___.solvers.map.(lower(solvers(k).name)) = k;
        if strcmpi( solvers(k).name, selected ),
            cvx___.solvers.map.default = k;
            cvx___.solvers.selected = k;
        end
    end
    if cvx___.solvers.selected == 0,
        cvx___.solvers.selected = 1;
        cvx___.solvers.map.default = 1;
        fprintf( [ ...
            'WARNING: The default solver %s is missing; %s has been selected as a\n', ...
            '    new default. If this was unexpected, try re-running cvx_setup.\n' ], ...
            selected, cvx___.solvers.list(cvx___.solvers.selected).name, lower(cvx___.solvers.list(cvx___.solvers.selected).name) );
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Create the global data structure %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    cvx_global;
    if isempty( cvx___.solvers.list ),
        error('CVX:Unexpected','No valid solvers were found.');
    end

    %%%%%%%%%%%%%%%%%%%%
    % Save preferences %
    %%%%%%%%%%%%%%%%%%%%

    fprintf( 'Saving updated preferences...' ); nret = true;
    try
        cvx_save_prefs( 'save' );
        fprintf( 'done.\n' ); nret = false;
    catch errmsg
        fprintf( 'unexpected error:\n' ); nret = false;
        cvx_error( errmsg, 67, true, '    ' );
        fprintf( 'Please attempt to correct this error and re-run cvx_setup. If you cannot,\n' );
        fprintf( 'you will be forced to re-run cvx_setup every time you start MATLAB.\n' );
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%
    % Test the distribution %
    %%%%%%%%%%%%%%%%%%%%%%%%%

    fprintf( 'Testing with a simple model...' ); nret = true;
    need_cc = false;
    try
        m = 16; n = 8;
        A = randn(m,n);
        b = randn(m,1);
        cvx_begin('quiet')
            variable('x(n)');
            minimize( norm(A*x-b,1) );
        cvx_end
        fprintf( 'done!\n' ); nret = false;
    catch exc
        if any(strfind(exc.message,'clear classes')),
            fprintf( 'problem (see below).\n' );
            need_cc = true;
        else
            rethrow( exc );
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Quick instructions on changing the solver %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if length( cvx___.solvers.list ) > 1,
        fprintf( '%s\n', line );
        fprintf( 'To change the default solver, type "cvx_solver <solver_name>".\n')
        fprintf( 'To save this change for future sessions, type "cvx_save_prefs".\n' );
        fprintf( 'Please consult the users'' guide for more information.\n' );
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Instruct the user to save the path %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if ~isempty( oldpath )
        need_upr = false;
        need_disclaim = true;
        user_path = which( 'startup.m' );
        if isempty( user_path ),
            user_path = userpath;
            if length(user_path) <= 1,
                need_upr = true;
                user_path = system_dependent('getuserworkfolder', 'default');
                if ~isempty( user_path ),
                    if isempty( strfind( user_path, [ fs, 'MATLAB' ] ) ),
                        user_path = [ user_path, fs, 'MATLAB' ];
                    end
                end
            elseif user_path(end) == ps,
                user_path(end) = '';
            end
            if ~isempty( user_path ),
                user_file = [ user_path, fs, 'startup.m' ];
            else
                user_file = '';
            end
        else
            user_file = user_path;
            user_path = user_path(1:end-10);
        end
        fprintf( '%s\n', line );
        fprintf('NOTE: the MATLAB path has been changed to point to the CVX distribution. To\n' );
        fprintf('use CVX without having to re-run CVX_SETUP every time MATLAB starts, you\n' );
        fprintf('will need to save this path permanently. This script attempted to do this\n' );
        if fs == '/',
            fprintf('for you, but failed---likely due to UNIX permissions restrictions.\n' );
            nextword = 'To solve the problem';
        else
            fprintf('for you, but failed, due to the Windows User Access Control (UAC) settings.\n');
            fprintf('Click here for a MATLAB document that discusses the issue.\n\n');
            fprintf('To solve this problem, please take the following steps:\n');
            fprintf('    1) Exit MATLAB.\n');
            fprintf('    2) Right click on the MATLAB icon, and select "Run as administrator."\n' );
            fprintf('    3) Re-run "cvx_setup".\n\n');
            nextword = 'Alternatively, if you do not have administrator access';
        end
        if exist( user_file, 'file' ),
            fprintf( '%s, edit the file\n    %s\nand add the following line to the end of the file:\n', nextword, user_file );
            fprintf( '    run %s%scvx_startup.m\n', mpath, fs );
        elseif ~isempty( user_path ),
            fprintf( '%s, create a new file\n    %s\ncontaining the following line:\n', nextword, user_file );
            fprintf( '    run %s%scvx_startup.m\n', mpath, fs );
        else
            fprintf( '%s, create a startup.m file containing the line:\n', nextword );
            fprintf( '    run %s%scvx_startup.m\n', mpath, fs );
            fprintf( 'Consult the MATLAB documentation for the proper location for that file.\n' );
            need_disclaim = false;
        end
        if need_upr,
            fprintf( 'Then execute the following MATLAB commands:\n    userpath reset; startup\n' );
        end
        if need_disclaim,
            fprintf( 'Please consult the MATLAB documentation for more information about the\n' );
            fprintf( 'startup.m file and its proper placement and usage.\n' );
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Warn about class conflict with previous version of CVX %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if need_cc,
        fprintf( '%s\n', line );
        fprintf('WARNING: CVX was unable to run the test model due to a conflict with the\n' );
        fprintf('previous version of CVX. If no other errors occurred, then the setup was\n' );
        fprintf('still successful; however, to use CVX, you will need to re-start MATLAB.\n' );
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Warn about signal processing toolbox conflict %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    squares = which( 'square', '-all' );
    if iscell( squares ) && length( squares ) > 1,
        squares = squares(~cellfun(@(x)any(strfind(x,[fs,'@cvx',fs])),squares));
        if length(squares) == 1 || ~strncmp(squares{1},[mpath,fs],length(mpath)+1),
            squares = {};
        end
    else
        squares = {};
    end
    if ~isempty(squares),
        fprintf( '%s\n', line );
        fprintf('WARNING: CVX includes a function\n    %s\n', squares{1} );
        fprintf('that conflicts with a function of the same name found here:\n    %s\n', squares{2} );
        fprintf('If you wish to use this second function, you will need to rename or delete\n' );
        fprintf('the CVX version, as it will likely produce different results. This will not\n' );
        fprintf('affect the use of SQUARE() within CVX models, because CVX relies on a\n' );
        fprintf('different, internal version of SQUARE() when constructing CVX models.\n')
    end

catch errmsg

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Restore the environment in the event of an error %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nret, fprintf( '\n' ); end
    switch errmsg.identifier,
        case { 'CVX:Expected', 'CVX:Licensing' },
            unexpected = false;
            if ~isempty( errmsg.message ),
                cvx_error( errmsg, 67, 'ERROR: ', '    ' );
            end
        case 'CVX:Unexpected',
            unexpected = true;
        otherwise,
            cvx_error( errmsg, 67, 'UNEXPECTED ERROR: ', '    ' );
            unexpected = true;
    end
    if ~isempty( oldpath ),
        path( oldpath );
    end
    clear global cvx___
    if unexpected,
        fprintf( 'Please report this error to support, and include entire output of\n' );
        fprintf( 'CVX_SETUP in your support request.\n' );
    else
        fprintf( 'The installation of CVX was not completed. Please correct the error\n' );
        fprintf( 'and re-run CVX_SETUP.\n' );
    end

end

fprintf( '%s\n\n', line );

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [subIndOut,L_valid]=snapSubInd(subIndIn,siz)

%Initialize outputs
L_valid=false(size(subIndIn));
subIndOut=subIndIn;

%Fixing out of range indices
for q=1:1:size(subIndIn,2)
    subInd=subIndIn(:,q); %The current subscript index set
    L_valid(:,q)=(subInd>0) & (subInd<=siz(q)); %Logic for valid indices not out of range

    subIndToFix=subInd(~L_valid(:,q)); %Indices to fix

    %Snapping out of range indices to boundary
    subIndToFix=(subIndToFix.*(subIndToFix>1))+(subIndToFix<1); %Fix smaller than 1
    subIndToFix=(subIndToFix.*(subIndToFix<=siz(q)))+siz(q).*(subIndToFix>siz(q)); %Fix larger than siz(q)

    %Storing fixed indices in output
    subIndOut_current=subInd;
    subIndOut_current(~L_valid(:,q))=subIndToFix;
    subIndOut(:,q)=subIndOut_current;
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [modelOut,nTotalAtomTransitions] = checkRXNFiles(model, RXNFileDir)
% Checks whether the set of RXN files coresponding to a model have the consistent stoichiometry and are elementally balanced
%
% INPUTS:
%    model:         Directed stoichiometric hypergraph
%                   Represented by a matlab structure with following fields:
%
%                     * .S - The `m` x `n` stoichiometric matrix for the metabolic network
%                     * .model.mets - An `m` x 1 array of metabolite identifiers. Should match
%                       metabolite identifiers in `RXNfiles`.
%                     * .model.rxns - An `n` x 1 array of reaction identifiers. Should match
%                       RXNfile names in `RXNFileDir`.
%                     * .lb -  An `n` x 1 vector of lower bounds on fluxes.
%                     * .ub - An `n` x 1 vector of upper bounds on fluxes.
%
%    RXNFileDir:    Path to directory containing `RXNfiles` with atom mappings
%                   for internal reactions in `S`. File names should
%                   correspond to reaction identifiers in input `model.rxns`.
%                   e.g. git clone https://github.com/opencobra/ctf ~/fork-ctf
%                        then RXNFileDir = ~/fork-ctf/rxns/atomMapped
%
% OUTPUT:
% metRXNBool: `m` x 1 vector, true if metabolite identified in at least one RXN file
% RXNBool: `n` x 1 boolean vector, true if RXN file exists
% RXNParsedBool: `n` x 1 boolean vector, true if RXN file could be parsed
% RXNAtomsConservedBool: `n` x 1 boolean vector, true if atoms in RXN file are conserved
% RXNStoichiometryMatchBool:  `n` x 1 boolean vector, true if RXN stoichiometry matches model.S stoichiometry
% RXNStoichiometryMatchUptoProtonsBool:   `n` x 1 boolean vector, true if RXN stoichiometry matches model.S stoichiometry when ingnoring protons.
% RXNSubstrateTransitionNumbersOrdered:   `n` x 1 boolean vector, true if RXN file with substrate transition numbers ordered 1:q.
% RXNProductTransitionNumbersOrdered:   `n` x 1 boolean vector, true if RXN file with product transition numbers ordered 1:q.
% RXNTransitionNumbersMatching:    `n` x 1 boolean vector, true if RXN file with matching numbering of atoms between substrates and products.
% RXNMatchingElementBool:     `n` x 1 boolean vector, true if RXN file with matching elements between substrates and products.
%
% .. Authors: - Ronan M. T. Fleming, 2022.

fprintf('Checking quality of RXN files...\n');

if ~exist('RXNFileDir','var')
    RXNFileDir=pwd;
end
RXNFileDir = [regexprep(RXNFileDir,'(/|\\)$',''), ''];

[nMets,nRxns]= size(model.S);
if length(unique(model.mets))~=nMets
    disp(setdiff(model.mets,unique(model.mets)))
    error('duplicate metabolites')
end
if length(unique(model.rxns))~=nRxns
    disp(setdiff(model.rxns,unique(model.rxns)))
    error('duplicate reactions')
end

[metRXNBool,RXNBool,internalRxnBool] = findRXNFiles(model,RXNFileDir);

%preallocate outputs
RXNParsedBool = NaN*ones(nRxns,1);
RXNAtomsConservedBool = NaN*ones(nRxns,1);
RXNStoichiometryMatchBool = NaN*ones(nRxns,1);
RXNStoichiometryMatchUptoProtonsBool = NaN*ones(nRxns,1);
RXNSubstrateTransitionNumbersOrdered = NaN*ones(nRxns,1);
RXNTransitionNumbersMatching = NaN*ones(nRxns,1);
RXNProductTransitionNumbersOrdered = NaN*ones(nRxns,1);
RXNMatchingElementBool = NaN*ones(nRxns,1);

%identify the protons in the atom mapped subset
pat = '[' + lettersPattern(1) + ']';
hBool = strcmp(model.mets,'h') | matches(model.mets,pat);

checkDecompartmentaliseRXN=1;
printAtomTransitionStats=0;
nTotalAtomTransitions = 0;
for i = 1:nRxns
    if RXNBool(i)

        if checkDecompartmentaliseRXN==1
            % Read atom mapping from RXNfile to test if it is decompartmentalised
            [atomMets,metEls, metNrs, atomTransitionNrs,isSubstrate,instances] = readRXNFile(model.rxns{1},RXNFileDir);

            decompartmentaliseRXN=0;
            atomMetAbbr  = atomMets{1};
            metAbbr = model.mets{1};
            if ~strcmp(atomMetAbbr(end),metAbbr(end))
                if strcmp(atomMetAbbr(end),']')
                    decompartmentaliseRXN=1;
                elseif strcmp(metAbbr(end),']')
                    for j=1:length(model.mets)
                        model.mets{j} = model.mets{j}(1:end-2);
                    end
                end
                checkDecompartmentaliseRXN=0;
            end
        end
        rxn = model.rxns{i};

        if strcmp(rxn,'AKGDm')
            %disp(rxn)
        end

        try
            %read in each RXN file
            [atomMets,metEls, metNrs, atomTransitionNrs,isSubstrate,instances] = readRXNFile(rxn,RXNFileDir);

            if decompartmentaliseRXN
                for k=1:length(atomMets)
                    atomMets{k,1}=atomMets{k,1}(1:end-3);
                end
            end
            RXNParsedBool(i)=1;
        catch ME
            RXNParsedBool(i)=0;
            fprintf('%s%s\n',rxn,' could not be parsed.');
            disp(getReport(ME))
        end

        %check there are an even number of atoms in the rxn file
        if mod(length(metEls),2)==0
            RXNAtomsConservedBool(i)=1;
        else
            RXNAtomsConservedBool(i)=0;
            rxnFormula = printRxnFormula(model, 'rxnAbbrList',model.rxns{i},'printFlag',0);
            fprintf('%s%s%s%s\n',rxn,' ',rxnFormula{1},' does not conserve atoms.');
        end

        % Check that stoichiometry in rxnfile matches the one in S
        rxnMets = unique(atomMets);
        ss = model.S(:,i);
        as = zeros(size(ss));
        for j = 1:length(rxnMets)
            rxnMet = rxnMets{j};

            if isSubstrate(strcmp(atomMets,rxnMet))
                as(strcmp(model.mets,rxnMet)) = -max(instances(strcmp(atomMets,rxnMet)));
            else

                as(strcmp(model.mets,rxnMet)) = max(instances(strcmp(atomMets,rxnMet)));
            end
        end
        if all(as == ss)
            RXNStoichiometryMatchBool(i)=1;
        else
            RXNStoichiometryMatchBool(i)=0;
            if all(as == ss  | hBool)
                fprintf('%s%s\n',rxn, ' stoichiometry matches upto protons.')
                RXNStoichiometryMatchUptoProtonsBool(i)=1;
            else
                RXNStoichiometryMatchUptoProtonsBool(i)=0;
                fprintf('%s%s\n',rxn, ' stoichiometry in model and RXNfile do not match:')
                fprintf('%s\t,', 'In model:')
                printRxnFormula(model,'rxnAbbrList',rxn)
                fprintf('%s\t,', 'In RXNfile:')
                model.S(:,ismember(model.rxns,rxn))=as;
                printRxnFormula(model,'rxnAbbrList',rxn)
                fprintf('\n');
            end
        end

        %checks specific for atom transitions
        if ~all(atomTransitionNrs==0)
            if printAtomTransitionStats == 0
                printAtomTransitionStats=1;
            end
            if all(sort(atomTransitionNrs(isSubstrate)) == (1:sum(isSubstrate))')
                RXNSubstrateTransitionNumbersOrdered(i)=1;
            else
                RXNSubstrateTransitionNumbersOrdered(i)=0;
                fprintf([rxn, '.rxn, Substrate transition numbers not ordered 1:q.\n'])
            end
            if all(all(sort(atomTransitionNrs(~isSubstrate)) == (1:sum(~isSubstrate))'))
                RXNProductTransitionNumbersOrdered(i)=1;
            else
                RXNProductTransitionNumbersOrdered(i)=0;
                fprintf([rxn, '.rxn, Product transition numbers not ordered 1:q.\n'])
            end
            if all(sort(atomTransitionNrs(isSubstrate)) == sort(atomTransitionNrs(~isSubstrate)))
                RXNTransitionNumbersMatching(i)=1;
            else
                RXNTransitionNumbersMatching(i)=0;
                fprintf([rxn, '.rxn, Substrate and product transition numbers not matching order 1:q.\n'])
            end

            nAtomTransitions = max(atomTransitionNrs);
            matchingElementBool=false(nAtomTransitions,1);
            for k=1:nAtomTransitions
                if strcmp(metEls(atomTransitionNrs==k & isSubstrate),metEls(atomTransitionNrs==k & ~isSubstrate))
                    matchingElementBool(k)=1;
                end
            end
            if all(matchingElementBool)
                RXNMatchingElementBool(i)=1;
            else
                RXNMatchingElementBool(i)=0;
                fprintf('%s%s%s%s%u%s\n',RXNfileName,' ',rxnFormula,' contains ', nnz(~matchingElementBool), ' atom transitions violating elemental conservation.');
            end
            %calculate the total number of atom transitions
            nTotalAtomTransitions  = nTotalAtomTransitions + nAtomTransitions;
        end
    end
end

if printAtomTransitionStats
    fprintf('\n%s\n','RXN file with atom mapping summary:')
else
    fprintf('\n%s\n','RXN file summary:')
end
fprintf('%u%s\n',nMets,' metabolites in model')
fprintf('%u%s\n',nnz(metRXNBool),' metabolites in atom mapped subset of model')
fprintf('%u%s\n',nRxns,' reactions in model')
fprintf('%u%s\n',nnz(internalRxnBool),' internal reactions in model')
fprintf('%u%s\n',nnz(RXNBool),' reactions in atom mapped subset of model')
fprintf('%u%s\n',nnz(RXNParsedBool==1),' RXN files parsed.')
fprintf('%u%s\n',nnz(RXNAtomsConservedBool==1),' RXN files with the number of atom conserved between substrates and products.')
fprintf('%u%s\n',nnz(RXNStoichiometryMatchBool==1),' RXN files matching stoichiometry of the model.')
fprintf('%u%s\n',nnz(internalRxnBool & RXNStoichiometryMatchUptoProtonsBool==1),' RXN files matching stoichiometry of the model, upto protons.')
fprintf('%u%s\n',nnz(internalRxnBool & RXNStoichiometryMatchBool~=1 & ~(RXNStoichiometryMatchUptoProtonsBool==1)),' RXN files not matching stoichiometry of model, even ignoring protons.')
if printAtomTransitionStats
    fprintf('%u%s\n',nnz(RXNSubstrateTransitionNumbersOrdered==1),' RXN files with substrate transition numbers ordered 1:q.')
    fprintf('%u%s\n',nnz(RXNProductTransitionNumbersOrdered==1),' RXN files with product transition numbers ordered 1:q.')
    fprintf('%u%s\n',nnz(RXNTransitionNumbersMatching==1),' RXN files with matching numbering of atoms between substrates and products.')
    fprintf('%u%s\n',nnz(RXNMatchingElementBool==1),' RXN files with matching elements between substrates and products.')
end

model.metRXNBool=metRXNBool;
model.RXNBool=RXNBool;
model.RXNParsedBool=RXNParsedBool;
model.RXNAtomsConservedBool=RXNAtomsConservedBool;
model.RXNStoichiometryMatchBool=RXNStoichiometryMatchBool;
model.RXNStoichiometryMatchUptoProtonsBool=RXNStoichiometryMatchUptoProtonsBool;
model.RXNSubstrateTransitionNumbersOrdered=RXNSubstrateTransitionNumbersOrdered;
model.RXNProductTransitionNumbersOrdered=RXNProductTransitionNumbersOrdered;
model.RXNTransitionNumbersMatching=RXNTransitionNumbersMatching;
model.RXNMatchingElementBool=RXNMatchingElementBool;
modelOut=model;

author: opencobra

repo: cobratoolbox


text:
function [curve, params, errors, T] = fit_poly_to_fragment(fragment, order)
%
%[curve, params, errors] = fit_poly_to_fragment(fragment, order)
%
% Fit a polynomial curve of specified order to an edge fragment.  A
% fragment is simply an Nx2 vector of (x,y) coordinates.
%
% Can also return fit error and the polynomial parameters.
%

% Ensure there are enough points:
if(isempty(fragment))
    error('Supplied fragment contains no points!');
end

N = size(fragment,1);

% Reduce the order if not enough points are provided to fit the requested
% order polynomial (e.g. we need N=4 points to fit a cubic, if only N=3
% points are provided, this will reduce the order to 2, thereby enabling
% successful fitting of a quadratic)
while( N < (order+1) )
    order = order-1;
end

t = linspace(0,1,N)';

fit_normals = false;
% if(nargin==2)
%     order = normal_angles;
%     fit_normals = false;
% end

T = ones(N,1);
for(i_order = 1:order)
    T = [t.^i_order T];
end

weights = ones(N, 1);
weights(1) = N/4;
weights(end) = N/4;
if(~fit_normals)
    % Fit the polynomial so that it simply tries to pass through the
    % fragment's vertex coordinates.
    params = (T .* repmat(weights, [1, order+1])) \ (fragment .* repmat(weights, [1, 2]));

%     % Constrained least squares to get the start and end points to exactly
%     % match up to the input coordinates:
%     for(i=1:2)
%         params(:,i) = lsqlin(T(2:end-1,:), fragment(2:end-1,i), ...
%             T([1 end],:), fragment([1 end],i),T([1 end],:), fragment([1 end],i));
%     end
else
    error(['Um, actually, trying to constrain the slopes this way will NOT ' ...
           'work. It is not linear because we do not know the length of the ' ...
           'normal vectors, only the direction.']);

%     params = T \ fragment;
%     L = L_helper(params, t, order);
%
%     % Fit the polynomial so that it tries to pass through the coordinates of the
%     % input fragment AND tries to have slope matching the orientations of
%     % each vertex in the fragment.
%     T_orient = [ones(N,1) zeros(N,1)];
%     for(i_order = 1:(order-1))
%         T_orient = [(i_order+1)*t.^i_order T_orient];
%     end
%
%     for(i=1:100)
%         TT = [blkdiag(T,T); blkdiag(-T_orient, T_orient)];
%
%         params = TT \ [fragment(:); L.*sin(normal_angles); L.*cos(normal_angles)];
%         params = reshape(params, [order+1 2]);
%
%         L = L_helper(params, t, order);
%     end

end

curve = T*params;
if(nargout>=3)
    errors = abs(curve - fragment);
end

return;

function L = L_helper(params, t, order)

dx = params(end-1,1);
dy = params(end-1,2);
for(i_order = 2:order)
    dx = dx + i_order*params(end-i_order,1)*t.^(i_order-1);
    dy = dy + i_order*params(end-i_order,2)*t.^(i_order-1);
end

L = sqrt(dx.^2 + dy.^2);
return;

author: Cloud-CV

repo: object-proposals


text:
% prettyMap 2.0
%--------------------------------------------------------------------------
%
% plot with projection an area of the world
%
% POSSIBLE SINTAXES:
%   prettyPolar(map);
%
%   prettyPolar(map, shape);
%   prettyPolar(map, projection);
%   prettyPolar(map, lineCol);
%
%   prettyPolar(map, phiGrid, lambdaGrid);
%   prettyPolar(map, shape, projection);
%   prettyPolar(map, shape, lineCol);
%   prettyPolar(map, projection, shape);
%   prettyPolar(map, projection, lineCol);
%
%   prettyPolar(map, phiGrid, lambdaGrid, shape);
%   prettyPolar(map, phiGrid, lambdaGrid, projection);
%   prettyPolar(map, phiGrid, lambdaGrid, lineCol);
%
%   prettyPolar(map, phiMin, phiMax, lambdaMin, lambdaMax)
%   prettyPolar(map, phiMin, phiMax, shape, projection);
%   prettyPolar(map, phiMin, phiMax, shape, lineCol);
%   prettyPolar(map, phiMin, phiMax, projection, shape);
%   prettyPolar(map, phiMin, phiMax, projection, lineCol);
%
%   prettyPolar(map, phiMin, phiMax, lambdaMin, lambdaMax, shepe);
%   prettyPolar(map, phiMin, phiMax, lambdaMin, lambdaMax, projection);
%   prettyPolar(map, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
%
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax)
%   prettyPolar(map, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
%   prettyPolar(map, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
%   prettyPolar(map, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
%   prettyPolar(map, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
%
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shepe);
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection);
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
%
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
%
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape, lineCol);
%   prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection, lineCol);
%
% EXAMPLE:
%   prettyPolar(map, phiMin, phiMax, lambdaMin, lambdaMax, 'Miller Cylindrical');
%
% INPUT:
%   map             matrix containing data of the whole world to be shown
%   phiGrid         array [degree]
%   lambdaGrid      array [degree]
%   phiMin          minimum latitude    [degree]
%   phiMax          maximum latitude    [degree]
%   lambdaMin       minimum longitude   [degree]
%   lambdaMax       maximum longitude   [degree]
%   projection      type of projection to be used "standard" is the default
%   shape           shapefile to load as coast (or country) contour
%                       - fill          fill coasts coarse
%                       - coast         only coasts coarse
%                       - 50m           1:50000000 scale country contours
%                       - 30m           1:30000000 scale country contours
%                       - 10m           1:10000000 scale country contours
%   lineCol         [1 1 1] array of RGB component to draw the contour lines
%
% DEFAULT VALUES:
%    projection = 'Lambert'
%
% AVAILABLE PROJECTION:
%    * Lambert
%      Stereographic
%      Orthographic
%      Azimuthal Equal-area
%      Azimuthal Equidistant
%      Gnomonic
%      Satellite
%      Albers Equal-Area Conic
%      Lambert Conformal Conic
%      Mercator
%    * Miller Cylindrical
%    * Equidistant Cylindrical (world map)
%      Oblique Mercator
%      Transverse Mercator
%      Sinusoidal
%      Gall-Peters
%      Hammer-Aitoff
%      Mollweide
%      Robinson
%    * UTM
%
% SEE ALSO:
%   mapPlot, mapPlot3D, quiver
%
% REQUIREMENTS:
%   M_Map: http://www.eos.ubc.ca/~rich/map.html
%   shape files with contours
%
% VERSION: 2.1
%
% CREDITS:
%   http://www.eos.ubc.ca/~rich/map.html
%
%   Andrea Gatti
%   DIIAR - Politecnico di Milano
%   2013-12-19
%
function prettyPolar(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape, lineCol)

%shape = 'coast';
%shape = 'fill';
%shape = '10m';
%shape = '30m';
%shape = '50m';

% lineCol = [0 0 0];
limitsOk = false;

% Manage opening a new figure;
tohold = false;
if length(findall(0,'Type','figure'))>=1
    if ishold
        clf;
        tohold = true;
    else
        figure;
    end
end

switch (nargin)
    case 1                                                                % prettyMap(map);
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        phiMin = 90;
        phiMax = -90;
        lambdaMin = -180;
        lambdaMax = 180;

        deltaPhi = (phiMax-phiMin)/size(map,1);
        deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

        phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
        lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
    case 2
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        if (ischar(phiGrid))
            if (sum(strcmp(phiGrid,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}]))) % prettyMap(map, shape);
                shape = phiGrid;
            else                                                          % prettyMap(map, projection);
                projection = phiGrid;
            end
        elseif (length(phiGrid) == 3)                                     % prettyMap(map, lineCol);
            lineCol = phiGrid;
        end

        phiMin = 90;
        phiMax = -90;
        lambdaMin = -180;
        lambdaMax = 180;

        deltaPhi = (phiMax-phiMin)/size(map,1);
        deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

        phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
        lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
    case 3
        shape = 'coast';
        lineCol = [0 0 0];
        if (ischar(phiGrid))
            projection = 'Miller Cylindrical';
            if (sum(strcmp(phiGrid,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = phiGrid;
                if (ischar(lambdaGrid))
                    projection = lambdaGrid;                              % prettyMap(map, shape, projection);
                else
                    lineCol = lambdaGrid;                                 % prettyMap(map, shape, lineCol);
                end
            else
                projection = phiGrid;
                if (ischar(lambdaGrid))
                    shape = lambdaGrid;                                   % prettyMap(map, projection, shape);
                else
                    lineCol = lambdaGrid;                                 % prettyMap(map, projection, lineCol);
                end
            end
            phiMin = 90;
            phiMax = -90;
            lambdaMin = -180;
            lambdaMax = 180;

            deltaPhi = (phiMax-phiMin)/size(map,1);
            deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

            phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
            lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
        else                                                              % prettyMap(map, phiGrid, lambdaGrid);
            projection = 'Miller Cylindrical';
            phiMin = max(phiGrid);
            phiMax = min(phiGrid);
            lambdaMin = min(lambdaGrid);
            lambdaMax = max(lambdaGrid);
        end
    case 4
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        if (ischar(phiMin))
            if (sum(strcmp(phiMin,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))  % prettyMap(map, phiGrid, lambdaGrid, shape);
                shape = phiMin;
            else                                                          % prettyMap(map, phiGrid, lambdaGrid, projection);
                projection = phiMin;
            end
        elseif (length(phiMin) == 3)                                      % prettyMap(map, phiGrid, lambdaGrid, lineCol);
            lineCol = phiMin;
        end

        phiMin = max(phiGrid);
        phiMax = min(phiGrid);
        lambdaMin = min(lambdaGrid);
        lambdaMax = max(lambdaGrid);
    case 5
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        if (ischar(phiMin))
            if (sum(strcmp(phiMin,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = phiMin;
                if (ischar(phiMax))
                    projection = phiMax;                                  % prettyMap(map, phiMin, phiMax, shape, projection);
                else
                    lineCol = phiMax;                                     % prettyMap(map, phiMin, phiMax, shape, lineCol);
                end
            else
                projection = phiMin;
                if (ischar(phiMax))
                    shape = phiMax;                                       % prettyMap(map, phiMin, phiMax, projection, shape);
                else
                    lineCol = phiMax;                                     % prettyMap(map, phiMin, phiMax, projection, lineCol);
                end
            end
            phiMin = max(phiGrid);
            phiMax = min(phiGrid);
            lambdaMin = min(lambdaGrid);
            lambdaMax = max(lambdaGrid);
        else                                                             %  prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax);
            limitsOk = true;
            lambdaMin = phiMin;
            lambdaMax = phiMax;
            phiMin = phiGrid;
            phiMax = lambdaGrid;

            if (phiMin < phiMax)
                tmp = phiMin;
                phiMin = phiMax;
                phiMax = tmp;
            end

            deltaPhi = (phiMax-phiMin)/size(map,1);
            deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

            phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
            lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
         end
    case 6
        shape = 'coast';
        lineCol = [0 0 0];
        limitsOk = true;
        projection = 'Lambert';
        if (ischar(lambdaMin))
            if (sum(strcmp(lambdaMin,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = lambdaMin;                                        % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, shape);
            else
                projection = lambdaMin;                                   % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, projection);
            end
        elseif (length(lambdaMin) == 3)
            lineCol = lambdaMin;                                          % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
        end

        lambdaMax = phiMax;
        lambdaMin = phiMin;
        phiMin = phiGrid;
        phiMax = lambdaGrid;

        if (phiMin < phiMax)
            tmp = phiMin;
            phiMin = phiMax;
            phiMax = tmp;
        end

        deltaPhi = (phiMax-phiMin)/size(map,1);
        deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

        phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
        lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
    case 7
        shape = 'coast';
        lineCol = [0 0 0];
        limitsOk = true;
        projection = 'Lambert';
        if (ischar(lambdaMin))
            if (sum(strcmp(lambdaMin,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = lambdaMin;
                if (ischar(lambdaMax))
                    projection = lambdaMax;                               % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
                else
                    lineCol = lambdaMax;                                  % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
                end
            else
                projection = lambdaMin;
                if (ischar(lambdaMax))
                    shape = lambdaMax;                                    % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
                else
                    lineCol = lambdaMax;                                  % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
                end
            end

            lambdaMin = phiMin;
            lambdaMax = phiMax;
            phiMin = phiGrid;
            phiMax = lambdaGrid;

            if (phiMin < phiMax)
                tmp = phiMin;
                phiMin = phiMax;
                phiMax = tmp;
            end

            deltaPhi = (phiMax-phiMin)/size(map,1);
            deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

            phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
            lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
        else
            projection = 'lambert';                                       % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax);
        end
    case 8
        shape = 'coast';
        lineCol = [0 0 0];
        limitsOk = true;
        if (ischar(projection))
            if (sum(strcmp(projection,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = projection;                                       % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape);
                projection = 'Lambert';
            else
                                                                          % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection);
            end
        elseif (length(projection) == 3)
            lineCol = projection;                                         % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
            projection = 'Lambert';
        end
    case 9
        lineCol = [0 0 0];
        limitsOk = true;
        if (ischar(projection))
            if (sum(strcmp(projection,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                tmp = shape;
                shape = projection;
                if (ischar(tmp))
                    projection = tmp;                                     % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
                else
                    lineCol = tmp;                                        % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
                    projection = 'UTM';
                end
            else
                if (ischar(shape))
                                                                          % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
                else
                    lineCol = shape;                                      % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
                    shape = 'coast';
                end
            end
        end
    case 10                                                               % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape, lineCol)
        limitsOk = true;
       if (sum(strcmp(projection,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))   % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection, lineCol)
            tmp = shape;
            shape = projection;
            projection = tmp;
       end
end

if (phiMin < phiMax)
    tmp = phiMin;
    phiMin = phiMax;
    phiMax = tmp;
end

lambdaGrid = sort(lambdaGrid);

[val idMax] = max(diff(lambdaGrid));
if sum(diff(lambdaGrid) == val) == 1
    lambdaGrid(1:idMax) = lambdaGrid(1:idMax)+360;
    if ~limitsOk
        lambdaMax = lambdaGrid(idMax);
        lambdaMin = lambdaGrid(idMax+1);
    end
end
lambdaGrid = sort(lambdaGrid);

if(lambdaMax<lambdaMin)
    lambdaMax = lambdaMax+360;
end

projection = 'Transverse Mercator';

% setup the projection
if (strcmpi(projection,[{'lambert'}]) && abs(phiMax==-phiMin))
    projection='Miller Cylindrical';
end

if (sum(strcmpi(projection,[{'lambert'},{'UTM'},{'Sinusoidal'},{'Transverse Mercator'},{'Oblique Mercator'},{'Miller Cylindrical'}])))
    m_proj(projection,'long',[lambdaMin lambdaMax],'lat',[phiMax phiMin]);
else
    m_proj(projection);
end

% Printing projection
fprintf('Using projection: %s\n', projection);

if sum(diff(lambdaGrid)<-200)
    lambdaGrid(lambdaGrid<0)=lambdaGrid(lambdaGrid<0)+360;
end
% plot the map
m_pcolor(lambdaGrid,phiGrid, map);

% set the light
shading flat;

% read shapefile
if (~strcmp(shape,'none'))
	if (~strcmp(shape,'coast')) && (~strcmp(shape,'fill'))
    	if (strcmp(shape,'10m'))
        	M=m_shaperead('countries_10m');
	    elseif (strcmp(shape,'30m'))
    	    M=m_shaperead('countries_30m');
	    else
    	    M=m_shaperead('countries_50m');
	    end
    	[xMin,yMin] = m_ll2xy(lambdaMin,phiMin);
	    [xMax,yMax] = m_ll2xy(lambdaMax,phiMax);
    	for k=1:length(M.ncst)
        	lamC = M.ncst{k}(:,1);
	        ids = lamC < lambdaMin;
   	    	lamC(ids) = lamC(ids) + 360;
        	phiC = M.ncst{k}(:,2);
        	[x,y] = m_ll2xy(lamC,phiC);
        	if sum(~isnan(x))>1
            	x(find(abs(diff(x))>=abs(xMax-xMin)*0.90)+1) = nan; % Romove lines that occupy more than th 90% of the plot
	            line(x,y,'color', lineCol);
    	    end
	    end;
    else
        if (strcmp(shape,'coast'))
        	m_coast('line','color', lineCol);
        else
            m_coast('patch',lineCol);
        end
	end
end

m_grid('box','fancy','tickdir','in');
drawnow;
if (phiMax >= 0)
    h=get(gca,'Children'); hf = findobj(h(1:end),'Tag','m_grid_fancybox1');
    delete(h(18));
    delete(hf([1:5 8]));
else
    h=get(gca,'Children');
    hf = findobj(h(1:end),'Tag','m_grid_fancybox1');
    %delete(h(18));
    delete(hf([1:4 6 7]));
end
colorbar;

if tohold
    hold on;
end

author: goGPS-Project


text:
%MAXPOOLING implements a max-pooling operation for 4D tensors shaped as:
% (rows, cols, channels, samples) and returns the pooled data with the
% corresponding indices.
%
%   [m, idx] = MaxPooling(IM, [2 2])
%
% IM can also be a 2D tensor, the missing dims are set to 1.

% AUTORIGHTS
% Copyright (C) 2011 Jonathan Masci <jonathan@idsia.ch>
%
% This file is available under the terms of the
% GNU GPLv2.

author: kyunghyuncho

repo: deepmat


text:
function i4vec_index_delete_one_test ( )

%*****************************************************************************80
%
%% I4VEC_INDEX_DELETE_ONE_TEST tests I4VEC_INDEX_DELETE_ONE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n_max = 25;
  n = 0;
  x = [];
  indx = [];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_INDEX_DELETE_ONE_TEST\n' );
  fprintf ( 1, '  I4VEC_INDEX_DELETE_ONE deletes one copies of a\n' );
  fprintf ( 1, '  particular value.\n' );

  xval = 8;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  xval = 7;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  b = 0;
  c = 20;
  seed = 123456789;

  for i = 1 : 20
    [ xval, seed ] = i4_uniform_ab ( b, c, seed );
    fprintf ( 1, '  %6d\n', xval );
    [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );
  end

  xval = 7;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  xval = 8;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)  X(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %6d  %6d  %6d\n', i, indx(i), x(i), x(indx(i)) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Call I4VEC_INDEX_DELETE_ONE to delete a value of 8:\n' );

  xval = 8;
  [ n, x, indx ] = i4vec_index_delete_one ( n, x, indx, xval );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)  X(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %6d  %6d  %6d\n', i, indx(i), x(i), x(indx(i)) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function result = cube_unit_3d ( func )

%*****************************************************************************80
%
%% CUBE_UNIT_3D approximates an integral inside the unit cube in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that:
%
%      -1 <= X <= 1,
%      -1 <= Y <= 1,
%      -1 <= Z <= 1.
%
%  Discussion:
%
%    An 8 point third degree formula is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function which evaluates F(X,Y,Z), of the form
%      function value = func ( x, y, z )
%
%    Output, real RESULT, the approximate integral of the function.
%
  s = 1.0E+00 / sqrt ( 3.0E+00 );
  w = 1.0E+00 / 8.0E+00;

  x = s;
  y = s;
  z = s;

  quad = w * ( ...
      feval ( func,  x,  y,  z ) ...
    + feval ( func,  x,  y, -z ) ...
    + feval ( func,  x, -y,  z ) ...
    + feval ( func,  x, -y, -z ) ...
    + feval ( func, -x,  y,  z ) ...
    + feval ( func, -x,  y, -z ) ...
    + feval ( func, -x, -y,  z ) ...
    + feval ( func, -x, -y, -z ) );

  volume = cube_unit_volume_nd ( 3 );
  result = quad * volume;

  return
end


text:
function z = entropy(x)
% Compute entropy H(x) of a discrete variable x.
% Written by Mo Chen (mochen80@gmail.com).
    n = numel(x);
    x = reshape(x,1,n);
    [u,~,label] = unique(x);
    p = full(mean(sparse(1:n,label,1,n,numel(u),n),1));
    z = -dot(p,log2(p+eps));
endfunction

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [diagnostics,out] = VBA_getDiagnostics(posterior,out)
% derives post-hoc diagnostics of VBA's model inversion
% function [diagnostics,out] = VBA_getDiagnostics(posterior,out)
% Post-hoc diagnostics include: goodness-of-fit metrics, Volterra kernels,
% null model evidence, posterior entropies, etc...
% IN:
%   - posterior,out: VBA's output structures
% OUT:
%   - diagnostics: a structure containing the following fields:
%       .pgx: mean of the prior predictive density
%       .pvy: variance of the prior predictive density
%       .kernels: Volterra kernels (if options.kernelSize>0)
%       .efficiency: posterior entropies
%       .DKL: prior/posterior Kullback-Leibler divergences
%       .LLH0: log-evidence of the null model
%       .MT_x: microtime hidden states time series
%       .MT_gx: microtime predicted data time series
%       .microTime: microtime sampling grid
%       .sampleInd: sub-indexing of the microtime sampling grid
%       .dy: data residuals structure (e.g., autocorrelation...)
%       .dx: states innovations structure
%       .C: posterior correlation matrix

u = out.u;
y = out.y;

% get goodness-tof-fit metrics
try; out.fit; catch; out.fit = VBA_getFit(posterior,out); end

% derive Volterra kernels
if out.dim.n_t>1 && out.options.kernelSize>0
    try
        kernels = VBA_getVolterraKernels(posterior,out);
    catch
        VBA_disp('  *** could not derive kernels!\n',out.options);
        kernels = [];
    end
else
    kernels = [];
end

% get null model (H0) evidence
[LLH0] = VBA_LMEH0(y,out.options);

% Entropies and KL divergences
gsi = find([out.options.sources.type]==0) ;
if isempty(gsi)
    efficiency.sigma = NaN;
    DKL.sigma = NaN;
else
    for iG=1:numel(gsi)
        efficiency.sigma(iG) = -out.suffStat.Ssigma(iG);
        m0 = out.options.priors.a_sigma(iG)/out.options.priors.b_sigma(iG);
        v0 = out.options.priors.a_sigma(iG)/out.options.priors.b_sigma(iG)^2;
        m = posterior.a_sigma(iG)/posterior.b_sigma(iG);
        v = posterior.a_sigma(iG)/posterior.b_sigma(iG)^2;
        DKL.sigma(iG) = VBA_KL(m0,v0,m,v,'Gamma');
    end
end
if out.dim.n > 0 % hidden states and initial conditions
    efficiency.X = -out.suffStat.SX;
    efficiency.X0 = -out.suffStat.SX0;
    if isinf(out.options.priors.a_alpha) && isequal(out.options.priors.b_alpha,0)
        efficiency.alpha = NaN;
        DKL.alpha = NaN;
    else
        efficiency.alpha = -out.suffStat.Salpha;
        m0 = out.options.priors.a_alpha./out.options.priors.b_alpha;
        v0 = out.options.priors.a_alpha./out.options.priors.b_alpha^2;
        m = posterior.a_alpha(end)./posterior.b_alpha(end);
        v = posterior.a_alpha(end)./posterior.b_alpha(end)^2;
        DKL.alpha = VBA_KL(m0,v0,m,v,'Gamma');
    end
    try
        DKL.X = 0;
        for t=1:out.dim.n_t
            IN = out.options.params2update.x{t};
            m0 = out.options.priors.muX(IN,t);
            v0 = out.options.priors.SigmaX.current{t}(IN,IN);
            m = posterior.muX(IN,t);
            v = posterior.SigmaX.current{t}(IN,IN);
            DKL.X = DKL.X + VBA_KL(m0,v0,m,v,'Normal');
        end
    catch
        DKL.X = NaN;
    end
    IN = out.options.params2update.x0;
    m0 = out.options.priors.muX0(IN);
    v0 = out.options.priors.SigmaX0(IN,IN);
    m = posterior.muX0(IN);
    v = posterior.SigmaX0(IN,IN);
    DKL.X0 = VBA_KL(m0,v0,m,v,'Normal');
else
    efficiency.X = NaN;
    efficiency.X0 = NaN;
    efficiency.alpha = NaN;
    DKL.X = NaN;
    DKL.X0 = NaN;
    DKL.alpha = NaN;
end
if out.dim.n_phi > 0 % observation parameters
    efficiency.Phi = -out.suffStat.Sphi;
    IN = out.options.params2update.phi;
    m0 = out.options.priors.muPhi(IN);
    v0 = out.options.priors.SigmaPhi(IN,IN);
    if ~out.options.OnLine
        m = posterior.muPhi(IN);
        v = posterior.SigmaPhi(IN,IN);
    else
        m = posterior.muPhi(IN,end);
        v = posterior.SigmaPhi{end}(IN,IN);
    end
    DKL.Phi = VBA_KL(m0,v0,m,v,'Normal');
else
    efficiency.Phi = NaN;
    DKL.Phi = NaN;
end
if out.dim.n_theta > 0 % evolution parameters
    efficiency.Theta = -out.suffStat.Stheta;
    IN = out.options.params2update.theta;
    m0 = out.options.priors.muTheta(IN);
    v0 = out.options.priors.SigmaTheta(IN,IN);
    if ~out.options.OnLine
        m = posterior.muTheta(IN);
        v = posterior.SigmaTheta(IN,IN);
    else
        m = posterior.muTheta(IN,end);
        v = posterior.SigmaTheta{end}(IN,IN);
    end
    DKL.Theta = VBA_KL(m0,v0,m,v,'Normal');
else
    efficiency.Theta = NaN;
    DKL.Theta = NaN;
end

% get prior predictive density
try
    [muy,Vy] = VBA_getLaplace(u,out.options.f_fname,out.options.g_fname,out.dim,out.options,0,'diag');
catch
    muy =[];
    Vy = [];
end

% get micro-time posterior hidden-states estimates
try
    [MT_x,MT_gx,microTime,sampleInd] = VBA_microTime(posterior,u,out);
catch
    MT_x = [];
    MT_gx = [];
    microTime = 1:out.dim.n_t;
    sampleInd = 1:out.dim.n_t;
end

% get residuals structure: data noise
for iS = 1:numel(out.options.sources)
    % source wise
    ySource = out.options.sources(iS).out ;
    res = out.suffStat.dy(ySource,:);
    if out.options.sources(iS).type==0
        gi = find(iS==find([out.options.sources(:).type]==0));
        iQyt = out.options.priors.iQy(:,gi);
        res = getWeightedResiduals(res,iQyt);
    end
    % remove skipped
    res(out.options.isYout(ySource,:)==1) = NaN ;
    dy(iS).dy = res(:);
    dy(iS).R = VBA_spm_autocorr(res);
    dy(iS).m = VBA_nanmean(dy(iS).dy);
    dy(iS).v = VBA_nanvar(dy(iS).dy);
    [dy(iS).ny,dy(iS).nx] = hist(dy(iS).dy,10);
    dy(iS).ny = dy(iS).ny./sum(dy(iS).ny);
    d = diff(dy(iS).nx);
    d = abs(d(1));
    dy(iS).d = d;
    spgy = sum(exp(-0.5.*(dy(iS).m-dy(iS).nx).^2./dy(iS).v));
    dy(iS).grid = dy(iS).nx(1):d*1e-2:dy(iS).nx(end);
    dy(iS).pg = exp(-0.5.*(dy(iS).m-dy(iS).grid).^2./dy(iS).v);
    dy(iS).pg = dy(iS).pg./spgy;
    if  out.options.sources(iS).type==0
        igs = sum([out.options.sources(1:iS).type]==0);
        shat = posterior.a_sigma(igs)./posterior.b_sigma(igs);
        spgy = sum(exp(-0.5.*shat.*dy(iS).nx.^2));
        dy(iS).pg2 = exp(-0.5.*shat.*dy(iS).grid.^2);
        dy(iS).pg2 = dy(iS).pg2./spgy;
    end
end

% get residuals structure: state noise
if ~isempty(out.suffStat.dx)
    wdx = getWeightedResiduals(out.suffStat.dx,out.options.priors.iQx);
    dx.dx = wdx(:);
    dx.m = mean(dx.dx);
    dx.v = var(dx.dx);
    [dx.ny,dx.nx] = hist(dx.dx,10);
    dx.ny = dx.ny./sum(dx.ny);
    d = diff(dx.nx);
    d = abs(d(1));
    dx.d = d;
    spgy = sum(exp(-0.5.*(dx.m-dx.nx).^2./dx.v));
    dx.grid = dx.nx(1):d*1e-2:dx.nx(end);
    dx.pg = exp(-0.5.*(dx.m-dx.grid).^2./dx.v);
    dx.pg = dx.pg./spgy;
    ahat = posterior.a_alpha(end)./posterior.b_alpha(end);
    spgy = sum(exp(-0.5.*ahat.*dx.nx.^2));
    dx.pg2 = exp(-0.5.*ahat.*dx.grid.^2);
    dx.pg2 = dx.pg2./spgy;
else
    dx.dx = [];
end

% get parameters posterior correlation matrix
if out.dim.n > 0 && isinf(out.options.priors.a_alpha) && isequal(out.options.priors.b_alpha,0)
    S = out.suffStat.ODE_posterior.SigmaPhi;
else
    S = NaN*zeros(out.dim.n+out.dim.n_theta+out.dim.n_phi);
    ind = 0;
    if out.dim.n_phi > 0
        if iscell(posterior.SigmaPhi) % online version
            SP = posterior.SigmaPhi{end};
        else
            SP = posterior.SigmaPhi;
        end
        S(1:out.dim.n_phi,1:out.dim.n_phi) = SP;
        ind = out.dim.n_phi;
    end
    if out.dim.n_theta > 0
        if iscell(posterior.SigmaTheta) % online version
            SP = posterior.SigmaTheta{end};
        else
            SP = posterior.SigmaTheta;
        end
        S(ind+1:ind+out.dim.n_theta,ind+1:ind+out.dim.n_theta) = SP;
        ind = ind + out.dim.n_theta;
    end
    if out.dim.n > 0 && out.options.updateX0
        if iscell(posterior.SigmaX0) % online version
            SP = posterior.SigmaX0{end};
        else
            SP = posterior.SigmaX0;
        end
        S(ind+1:ind+out.dim.n,ind+1:ind+out.dim.n) = SP;
    end
end
C = VBA_cov2corr(S);
C = C + diag(NaN.*diag(C));
tick = [0];
ltick = [];
ticklabel = cell(0,0);
if out.dim.n_phi > 0
    ltick = [ltick,tick(end)+out.dim.n_phi/2];
    tick = [tick,out.dim.n_phi];
    ticklabel{end+1} = 'phi';
end
if out.dim.n_theta > 0
    ltick = [ltick,tick(end)+out.dim.n_theta/2];
    tick = [tick,tick(end)+out.dim.n_theta];
    ticklabel{end+1} = 'theta';
end
if out.dim.n > 0 && out.options.updateX0
    ltick = [ltick,tick(end)+out.dim.n/2];
    tick = [tick,tick(end)+out.dim.n];
    ticklabel{end+1} = 'x0';
end
tick = tick +0.5;
tick = tick(2:end-1);
ltick = ltick + 0.5;

% wrap up
diagnostics.pgx = reshape(muy,out.dim.p,[]);
diagnostics.pvy = reshape(Vy,out.dim.p,[]);
diagnostics.kernels = kernels;
diagnostics.efficiency = efficiency;
diagnostics.DKL = DKL;
diagnostics.LLH0 = LLH0;
diagnostics.MT_x = MT_x;
diagnostics.MT_gx = MT_gx;
diagnostics.microTime = microTime;
diagnostics.sampleInd = sampleInd;
diagnostics.dy = dy;
diagnostics.dx = dx;
diagnostics.ltick = ltick;
diagnostics.tick = tick;
diagnostics.ticklabel = ticklabel;
diagnostics.C = C;
out.diagnostics = diagnostics;

function wdx = getWeightedResiduals(dx,iQx)
% weigths residuals according to (state/data) precision matrix
wdx = zeros(size(dx));
for t = 1:size(dx,2)
    sqrtiQ = VBA_sqrtm (iQx{t});
    wdx(:,t) = sqrtiQ*dx(:,t);
end

author: MBB-team

repo: VBA-toolbox