text:
function [ExcelID,WorkbookID,SheetNum,shapeID,imageHandleList]=excelReport(view,comment,ExcelID,WorkbookID,SheetNum,imSize,startPos)
%
% [ExcelID,WorkbookID,SheetNum,shapeID,imageHandleList]= ...
%                 excelReport(comment,ExcelID,WorkbookID,SheetNum,imSize,startPos)
%
% PURPOSE: Automates the generation of an Excel spreadsheet with
% Screenshots of the current datatype / scan
% Assumes that mrmesh has been invoked to visualize the data

mrGlobals;

meshTypes={'Folded','Relaxed'};
if (~exist('imSize'))
imSize=[200 200];
end
if (~exist('startPos'))
startPos=[20 20];
end

% Script to generate screenshots
% Saves in current dir:

viewList={'back','left','right','bottom','top'};
viewVectors={[pi -pi/2 0],[pi 0 0],[0 0 pi],[pi/2 -pi/2 0],[-pi/2 -pi/2 0]};
view=getSelectedVolume;

idList = viewGet(view,'allwindowids');
thisID=view.meshNum3d;
% Set up an excel worksheet to accept the data
% If no handle has been passed, we initiate a new instance of excel

if (~exist('ExcelID','var'))
    ExcelID = actxserver('Excel.Application');
end

if (isempty(ExcelID))
    ExcelID = actxserver('Excel.Application');
end
disp(ExcelID)

set(ExcelID, 'Visible', 1);

if (~exist('WorkbookID','var'))
         Workbooks = ExcelID.Workbooks;
         WorkbookID = invoke(Workbooks, 'Add');
end
if (isempty(WorkbookID))
         Workbooks = ExcelID.Workbooks;
         WorkbookID = invoke(Workbooks, 'Add');
end

if (~exist('SheetNum','var'))
    SheetNum=1;

end

         Sheets = ExcelID.ActiveWorkBook.Sheets;
         SheetID = get(Sheets, 'Item', SheetNum);

        invoke(SheetID, 'Activate');

ActivesheetID = ExcelID.Activesheet;
shapeID=WorkbookID.ActiveSheet.Shapes

for thisView=1:length(viewList);
    cam.actor=0;
    cam.rotation=rotationMatrix3d(viewVectors{thisView})
    mrMesh('localhost',thisID,'set',cam)
    date_string=datestr(now,30);
    dt=dataTYPES(view.curDataType).name;
    scanNum=int2str(getCurScan(view));
    %mt=meshTypes{v.meshNum3d};
    filename=fullfile(pwd,[dt,'_',scanNum,'_',date_string,'_',viewList{thisView},'.bmp'])
    c.filename=fullfile(pwd,'test.bmp');

    pause(2);
    [a,b,result]=mrMesh('localhost',thisID,'screenshot',c)
    pause(5);
    disp(thisView)

     a=invoke(shapeID,'AddPicture',c.filename,1,1,startPos(1),startPos(2),imSize(1),imSize(2));
%    set(a,'AlternativeText',[comment,' : ',filename]);

    startPos(1)=startPos(1)+imSize(1);
    imageHandleList(thisView)=a;

end

% We can also copy and paste figures from matlab to the excel spreadsheet.

% To copy a figure from matlab try : print -dmeta


text:
%> @brief Return a set of indices separating contiguous segments of logical
%> values.
%> @author Eric Cousineau <eacousineau@gmail.com>
function [C] = contiguous_indices(B)

% Following: http://stackoverflow.com/q/2212201/170413
B = logical(B(:)');
count = length(B);

if all(B)
    % Special case for all values being high
    C = {1:count};
    return;
end

dB = diff(B);
starts = find(dB == 1) + 1;
ends = find(dB == -1);

% If we start high, then we need to insert an initial index
if length(ends) >= 1 && (isempty(starts) || ends(1) < starts(1))
    starts = [1, starts];
end

% The count for ends and starts may differ if there is a lone point at the
% end
if length(starts) > length(ends)
    if isempty(ends) || starts(end) > ends(end)
        % Means we ended at transition going from low to high
        ends(end + 1) = count;
    else
        % Augment ends with last starting point
        ends(end + 1) = starts(end);
    end
end

nranges = length(ends);
C = cell(1, nranges);
for i = 1:nranges
    C{i} = starts(i):ends(i);
end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [r,a] =  training(a,d)

  if isa(d,'group') & strcmp(d.group,'one_for_each') & length(d.child)==length(a.child)

        %% pass only nth data object into nth algorithm

    r=[];
    for i=1:length(a.child)
      [res,a.child{i}]=train(a.child{i},d.child{i}); r{i}=res;
    end

  else  %% standard method

    r=[];
    for i=1:length(a.child)
      [res,a.child{i}]=train(a.child{i},d); r{i}=res;
    end
  end

  if length(r)==1 r=r{1}; else r=group(r); end;
  if isa(r,'group') r.group=a.group; end;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function batchTestAllReconstructionFunctions(modelFolder,testResultsFolder,inputDataFolder,reconVersion, numWorkers)
% Part of the DEMETER pipeline. This function performs all quality
% control/quality assurance tests on a batch of reconstructions and saves
% the results for each reconstruction in the input folder.
%
% USAGE:
%
%   batchTestAllReconstructionFunctions(modelFolder,testResultsFolder,inputDataFolder,reconVersion, numWorkers)
%
% INPUTS
% modelFolder           Folder with COBRA models (draft or refined
%                       reconstructions) to analyze
% testResultsFolder     Folder where the test results should be saved
% inputDataFolder       Folder with experimental data and database files
%                       to load
% reconVersion          Name of the refined reconstruction resource
% numWorkers            Number of workers in parallel pool
%
% .. Author:
%   - Almut Heinken, 09/2020

% set a solver if not done yet
global CBT_LP_SOLVER
solver = CBT_LP_SOLVER;
if isempty(solver)
    initCobraToolbox(false); %Don't update the toolbox automatically
end

% initialize parallel pool
if numWorkers>0 && ~isempty(ver('parallel'))
    % with parallelization
    poolobj = gcp('nocreate');
    if isempty(poolobj)
        parpool(numWorkers)
    end
end
environment = getEnvironment();

% Runs through all tests
dInfo = dir(modelFolder);
modelList={dInfo.name};
modelList=modelList';
modelList(~contains(modelList(:,1),'.mat'),:)=[];

fields = {
    'Mass_imbalanced'
    'Charge_imbalanced'
    'Mets_without_formulas'
    'Leaking_metabolites'
    'ATP_from_O2'
    'Blocked_reactions'
    'RefinedReactionsCarryingFlux'
    'BlockedRefinedReactions'
    'Incorrect_Gene_Rules'
    'Incorrect_Compartments'
    'Carbon_sources_TruePositives'
    'Carbon_sources_FalseNegatives'
    'Fermentation_products_TruePositives'
    'Fermentation_products_FalseNegatives'
    'growsOnDefinedMedium'
    'growthOnKnownCarbonSources'
    'Biomass_precursor_biosynthesis_TruePositives'
    'Biomass_precursor_biosynthesis_FalseNegatives'
    'Metabolite_uptake_TruePositives'
    'Metabolite_uptake_FalseNegatives'
    'Secretion_products_TruePositives'
    'Secretion_products_FalseNegatives'
    'Bile_acid_biosynthesis_TruePositives'
    'Bile_acid_biosynthesis_FalseNegatives'
    'Drug_metabolism_TruePositives'
    'Drug_metabolism_FalseNegatives'
    'PutrefactionPathways_TruePositives'
    'PutrefactionPathways_FalseNegatives'
    };

%% load the results from existing test suite run and restart from there
Results=struct;

alreadyAnalyzedStrains={};

for i=1:length(fields)
    if isfile([testResultsFolder filesep fields{i} '_' reconVersion '.txt'])
        savedResults = readInputTableForPipeline([testResultsFolder filesep fields{i} '_' reconVersion '.txt']);
        Results.(fields{i}) = savedResults;
        if size(savedResults,2)>0
            alreadyAnalyzedStrains = Results.(fields{i})(:,1);
        end
    else
        Results.(fields{i})={};
    end
end

% propagate strains to empty fields
for i=1:length(fields)
    if isempty(Results.(fields{i}))
        Results.(fields{i})=alreadyAnalyzedStrains;
    end
end

% remove already analyzed reconstructions
if size(Results.(fields{1}),1)>0
    [C,IA]=intersect(strrep(modelList(:,1),'.mat',''),Results.(fields{1})(:,1));
    modelList(IA,:)=[];
end

% define the intervals in which the testing and regular saving will be
% performed
if length(modelList)>5000
    steps=1000;
elseif length(modelList)>1000
    steps=500;
elseif length(modelList)>200
    steps=200;
else
    steps=25;
end

%% Start the test suite
for i =1:steps:length(modelList)
    tmpData={};
    if length(modelList)>steps-1 && (length(modelList)-1)>=steps-1
        endPnt=steps-1;
    else
        endPnt=length(modelList)-i;
    end

    modelsToLoad={};
    for j=i:i+endPnt
        modelsToLoad{j}={};
        if j <= length(modelList)
            modelsToLoad{j}=[modelFolder filesep modelList{j}];
        end
    end
    parfor j=i:i+endPnt
%for j=i:i+endPnt
        if j <= length(modelList)
            restoreEnvironment(environment);
            changeCobraSolver(solver, 'LP');
            % prevent creation of log files
            changeCobraSolverParams('LP', 'logFile', 0);

            try
                model=readCbModel(modelsToLoad{j},'modelName','model');
            catch
                modelsToLoad{j}
                model=load(modelsToLoad{j});
                fieldnames(model)
                model=model.model;
            end

            microbeID=strrep(modelList{j},'.mat','');
            testResults = runTestsOnModel(model, microbeID, inputDataFolder);
            tmpData{j}=testResults;
        end
    end
    for k=1:length(fields)
        for j=i:i+endPnt
            if j <= length(modelList)
                res=tmpData{j};

                if strcmp(fields{k},'Blocked_reactions') && length(modelList) + length(alreadyAnalyzedStrains) > 10000
                    % do not save for very large-scale resources-file would be
                    % enormous
                    Results.(fields{k}){size(Results.(fields{k}),1)+1,1} = res.(fields{k}){1,1};
                    Results.(fields{k}){size(Results.(fields{k}),1),2} = size(res.(fields{k}),2)-1;
                else
                    Results.(fields{k})(size(Results.(fields{k}),1)+1,1:size(res.(fields{k}),2)) = res.(fields{k})(1,1:end);
                end
            end
        end
    end

    %% print the results regularly to avoid having to repeat simulations
    % only if there were any findings that need to be reported
    for j=1:length(fields)
        if size(Results.(fields{j}),2)>1 || any(contains(fields{j},{'TruePositives','FalseNegatives','growsOnDefinedMedium'}))
            table2print=cell2table(Results.(fields{j}));
            writetable(table2print,[testResultsFolder filesep fields{j} '_' reconVersion],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
        end
    end
end

end

author: opencobra

repo: cobratoolbox


text:
function mat2file(a,val,varargin)
% Function for writing to file_array objects
% FORMAT mat2file(a,val,ind1,ind2,ind3,...)
% a      - file_array object
% val    - values to write
% indx   - indices for dimension x (int32)
%
% This function is normally called by file_array/subsasgn.
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: mat2file.m 7147 2017-08-03 14:07:01Z spm $

%-This is merely the help file for the compiled routine
error('mat2file.c not compiled - see Makefile');

author: spm


text:
function [V2Dr,V2Ds] = GradVandermonde2D(N,r,s)

% function [V2Dr,V2Ds] = GradVandermonde2D(N,r,s)
% Purpose : Initialize the gradient of the modal basis (i,j) at (r,s) at order N

V2Dr = zeros(length(r),(N+1)*(N+2)/2); V2Ds = zeros(length(r),(N+1)*(N+2)/2);

% find tensor-product coordinates
[a,b] = rstoab(r,s);

% Initialize matrices
sk = 1;
for i=0:N
  for j=0:N-i
    [V2Dr(:,sk),V2Ds(:,sk)] = GradSimplex2DP(a,b,i,j);
    sk = sk+1;
  end
end
return;

repo: Aero-matlab


text:
function tan = p00_tan ( problem, option, nvar, x )

%*****************************************************************************80
%
%% P00_TAN determines a tangent vector at X.
%
%  Discussion:
%
%    If X is a solution of F(Y) = 0, then the vector TAN
%    is tangent to the curve of solutions at X.
%
%    If X is not a solution of F(Y) = 0, then the vector TAN
%    is tangent to the curve F(Y) = F(X) at X.
%
%    The vector will have unit euclidean norm.
%
%    The sign of TAN will be chosen so that the determinant
%    of F'(X) augmented with a final row equal to TAN will be positive.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the evaluation point.
%
%    Output, real TAN(NVAR), a tangent vector at X.
%

%
%  Compute the jacobian.
%
  jac = p00_jac ( problem, option, nvar, x );
%
%  Compute the QR factorization of JAC'.
%
  [ q, r ] = qr ( jac' );

  tan(1:nvar,1) = q(1:nvar,nvar);
%
%  Choose the sign of TAN by the determinant condition.
%
  jac(nvar,1:nvar) = tan(1:nvar,1);

  d = det ( jac );

  if ( d < 0.0 )
    tan(1:nvar) = - tan(1:nvar);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function buckling_spring_theta ( l_num, theta_num )

%*****************************************************************************80
%
%% BUCKLING_SPRING_THETA plots LAMBDA(L,THETA) and MU(L,THETA) lines for fixed L.
%
%  Discussion:
%
%    This routine fixes L and plots LAMBDA and MU as functions of THETA.
%
%  Usage:
%
%    buckling_spring_theta ( l_num, theta_num )
%
%    * l_num is the number of L values along each THETA line;
%    * theta_num is the number of THETA lines to draw.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer L_NUM, the number of values of L,
%    which controls the smoothness of each line.
%    A value of 101 is plenty.
%
%    Input, integer THETA_NUM, the number of values of THETA,
%    which sets the number of lines to be drawn.
%    A value of 101 is plenty.
%
  l_lo = 0.25;
  l_hi = 1.75;

  theta_lo = - 3 * pi / 8;
  theta_hi = + 3 * pi / 8;

  for j = 1 : l_num

    if ( l_num == 1 )
      L = 0.5 * ( l_lo + l_hi );
    else
      L = ( ( l_num - j     ) * l_lo ...
          + (         j - 1 ) * l_hi ) ...
          / ( l_num     - 1 );
    end

    for i = 1 : theta_num

      if ( theta_num == 1 )
        theta = 0.5 * ( theta_lo + theta_hi );
      else
        theta = ( ( theta_num - i     ) * theta_lo ...
                + (             i - 1 ) * theta_hi ) ...
                / ( theta_num     - 1 );
      end

      lambda(i,j) = ( 1 - L ) * cos ( theta ) + theta * sin ( theta ) / 4 / L;
      mu(i,j) = - theta * cos ( theta ) / 2 / L + 2 * ( 1 - L ) * sin ( theta );

    end

  end

  plot ( lambda, mu )
%
%  Limit the plot range.
%  Suppress this command to see the whole plot.
%
  axis ( [ 0.1 0.9 -0.07 0.07 ] );
%
%  Label the plot.
%
  xlabel ( '\lambda(L,\theta)' );
  ylabel ( '\mu(L,\theta)' );
  title ( 'Each line is a fixed value of L' );

  return
end


text:
function om = opf_setup(mpc, mpopt)
%OPF  Constructs an OPF model object from a MATPOWER case struct.
%   OM = OPF_SETUP(MPC, MPOPT)
%
%   Assumes that MPC is a MATPOWER case struct with internal indexing,
%   all equipment in-service, etc.
%
%   See also OPF, EXT2INT, OPF_EXECUTE.

%   MATPOWER
%   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% options
dc  = strcmp(upper(mpopt.model), 'DC');
alg = upper(mpopt.opf.ac.solver);
use_vg = mpopt.opf.use_vg;
vcart = ~dc && mpopt.opf.v_cartesian;

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

%% define flag to indicate whether we are tied to legacy formulation
%% implemented by legacy MINOS and PDIPM solvers (e.g. with hard-coded
%% costs and constrain
if strcmp(alg, 'MINOPF') || strcmp(alg, 'PDIPM') || ...
        strcmp(alg, 'TRALM') || strcmp(alg, 'SDPOPF')
    legacy_formulation = 1;
    if vcart
        error('Option ''opf.v_cartesian'' = 1 is not compatible with ''opf.solver.ac''=''%s''.', alg);
    end
    if mpopt.opf.current_balance
        error('Option ''opf.current_balance'' = 1 is not compatible with ''opf.solver.ac''=''%s''.', alg);
    end
else
    legacy_formulation = 0;
end
if ~dc && ( ~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ...
                    ~isequal(mpopt.exp.sys_wide_zip_loads.pw, [1 0 0]) || ...
            ~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ...
                    ~isequal(mpopt.exp.sys_wide_zip_loads.qw, [1 0 0]) )
    if vcart
        warning('Voltage dependent loads are not supported with option ''opf.v_cartesian'' = 1. Reverting to constant power load model.');
    end
    if mpopt.opf.current_balance
        warning('Voltage dependent loads are not supported with option ''opf.current_balance'' = 1. Reverting to constant power load model.');
    end
end

%% data dimensions
nb   = size(mpc.bus, 1);    %% number of buses
nl   = size(mpc.branch, 1); %% number of branches
ng   = size(mpc.gen, 1);    %% number of dispatchable injections
nnle = 0;                   %% number of nonlinear user-defined equality cons
nnli = 0;                   %% number of nonlinear user-defined inequality cons
if isfield(mpc, 'A')
  nlin = size(mpc.A, 1);    %% number of linear user constraints
else
  nlin = 0;
end
if isfield(mpc, 'N')
  nw = size(mpc.N, 1);      %% number of general cost vars, w
else
  nw = 0;
end

if dc
  %% ignore reactive costs for DC
  mpc.gencost = pqcost(mpc.gencost, ng);

  %% reduce A and/or N from AC dimensions to DC dimensions, if needed
  if nlin || nw
    acc = [nb+(1:nb) 2*nb+ng+(1:ng)];   %% Vm and Qg columns
    if nlin && size(mpc.A, 2) >= 2*nb + 2*ng
      %% make sure there aren't any constraints on Vm or Qg
      if any(any(mpc.A(:, acc)))
        error('opf_setup: attempting to solve DC OPF with user constraints on Vm or Qg');
      end
      mpc.A(:, acc) = [];               %% delete Vm and Qg columns
    end
    if nw && size(mpc.N, 2) >= 2*nb + 2*ng
      %% make sure there aren't any costs on Vm or Qg
      if any(any(mpc.N(:, acc)))
        [ii, jj] = find(mpc.N(:, acc));
        ii = unique(ii);    %% indices of w with potential non-zero cost terms from Vm or Qg
        if any(mpc.Cw(ii)) || (isfield(mpc, 'H') && ~isempty(mpc.H) && ...
                any(any(mpc.H(:, ii))))
          error('opf_setup: attempting to solve DC OPF with user costs on Vm or Qg');
        end
      end
      mpc.N(:, acc) = [];               %% delete Vm and Qg columns
    end
  end
else    %% AC
  if use_vg     %% adjust bus voltage limits based on generator Vg setpoint
    %% gen connection matrix, element i, j is 1 if, generator j at bus i is ON
    Cg = sparse(mpc.gen(:, GEN_BUS), (1:ng)', mpc.gen(:, GEN_STATUS) > 0, nb, ng);
    Vbg = Cg * sparse(1:ng, 1:ng, mpc.gen(:, VG), ng, ng);
    Vmax = max(Vbg, [], 2); %% zero for non-gen buses, else max Vg of gens @ bus
    ib = find(Vmax);                %% buses with online gens
    Vmin = max(2*Cg - Vbg, [], 2);  %% same as Vmax, except min Vg of gens @ bus
    Vmin(ib) = 2 - Vmin(ib);

    if use_vg == 1      %% use Vg setpoint directly
        mpc.bus(ib, VMAX) = Vmax(ib);   %% max set by max Vg @ bus
        mpc.bus(ib, VMIN) = Vmin(ib);   %% min set by min Vg @ bus
        mpc.bus(ib, VM) = mpc.bus(ib, VMAX);
    elseif use_vg > 0 && use_vg < 1     %% fractional value
        %% use weighted avg between original Vmin/Vmax limits and Vg
        mpc.bus(ib, VMAX) = (1-use_vg) * mpc.bus(ib, VMAX) + use_vg * Vmax(ib);
        mpc.bus(ib, VMIN) = (1-use_vg) * mpc.bus(ib, VMIN) + use_vg * Vmin(ib);
    else
        error('opf_setup: option ''opf.use_vg'' (= %g) cannot be negative or greater than 1', use_vg);
    end
  end
  if isfield(mpc, 'user_constraints')
    if isfield(mpc.user_constraints, 'nle')
      for k = 1:length(mpc.user_constraints.nle)
        nnle = nnle + mpc.user_constraints.nle{k}{2};
      end
    end
    if isfield(mpc.user_constraints, 'nli')
      for k = 1:length(mpc.user_constraints.nli)
        nnli = nnli + mpc.user_constraints.nli{k}{2};
      end
    end
  end
end

%% convert single-block piecewise-linear costs into linear polynomial cost
pwl1 = find(mpc.gencost(:, MODEL) == PW_LINEAR & mpc.gencost(:, NCOST) == 2);
% p1 = [];
if ~isempty(pwl1)
  x0 = mpc.gencost(pwl1, COST);
  y0 = mpc.gencost(pwl1, COST+1);
  x1 = mpc.gencost(pwl1, COST+2);
  y1 = mpc.gencost(pwl1, COST+3);
  m = (y1 - y0) ./ (x1 - x0);
  b = y0 - m .* x0;
  mpc.gencost(pwl1, MODEL) = POLYNOMIAL;
  mpc.gencost(pwl1, NCOST) = 2;
  mpc.gencost(pwl1, COST:COST+1) = [m b];
end

%% create (read-only) copies of individual fields for convenience
[baseMVA, bus, gen, branch, gencost, Au, lbu, ubu, mpopt, ...
    N, fparm, H, Cw, z0, zl, zu, userfcn] = opf_args(mpc, mpopt);

%% warn if there is more than one reference bus
refs = find(bus(:, BUS_TYPE) == REF);
if length(refs) > 1 && mpopt.verbose > 0
  errstr = ['\nopf_setup: Warning: Multiple reference buses.\n', ...
              '           For a system with islands, a reference bus in each island\n', ...
              '           may help convergence, but in a fully connected system such\n', ...
              '           a situation is probably not reasonable.\n\n' ];
  fprintf(errstr);
end

%% set up initial variables and bounds
Va   = bus(:, VA) * (pi/180);
Vau = Inf(nb, 1);       %% voltage angle limits
Val = -Vau;
Vau(refs) = Va(refs);   %% voltage angle reference constraints
Val(refs) = Va(refs);
Pg   = gen(:, PG) / baseMVA;
Pmin = gen(:, PMIN) / baseMVA;
Pmax = gen(:, PMAX) / baseMVA;
if ~dc
  Vm   = bus(:, VM);
  Qg   = gen(:, QG) / baseMVA;
  Qmin = gen(:, QMIN) / baseMVA;
  Qmax = gen(:, QMAX) / baseMVA;
  if vcart
    V = Vm .* exp(1j*Va);
    Vr = real(V);
    Vi = imag(V);
  end
end

%% find/prepare polynomial generator costs
cpg = [];
cqg = [];
%% reduce order where highest order has zero coefficient
gc = mpc.gencost;
while 1
    ii = find(  gc(:, MODEL) == POLYNOMIAL & ...
                gc(:, NCOST) > 2 & ...
                gc(:, COST) == 0    );
    if isempty(ii), break; end  %% exit loop if there are none
    gc(ii, NCOST) = gc(ii, NCOST) - 1;          %% reduce order
    gc(ii, COST:end-1) = gc(ii, COST+1:end);    %% shift coefficients
    gc(ii, end) = 0;
end
[pcost qcost] = pqcost(gc, ng);
ip0 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) == 1);   %% constant
ip1 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) == 2);   %% linear
ip2 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) == 3);   %% quadratic
ip3 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) > 3);    %% cubic or greater
if ~isempty(ip2) || ~isempty(ip1) || ~isempty(ip0)
    kpg = zeros(ng, 1);
    cpg = zeros(ng, 1);
    if ~isempty(ip2)
        Qpg = zeros(ng, 1);
        Qpg(ip2) = 2 * pcost(ip2, COST) * baseMVA^2;
        cpg(ip2) = cpg(ip2) + pcost(ip2, COST+1) * baseMVA;
        kpg(ip2) = kpg(ip2) + pcost(ip2, COST+2);
    else
        Qpg = [];   %% no quadratic terms
    end
    if ~isempty(ip1)
        cpg(ip1) = cpg(ip1) + pcost(ip1, COST) * baseMVA;
        kpg(ip1) = kpg(ip1) + pcost(ip1, COST+1);
    end
    if ~isempty(ip0)
        kpg(ip0) = kpg(ip0) + pcost(ip0, COST);
    end
end
if ~isempty(qcost)
    iq0 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) == 1);   %% constant
    iq1 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) == 2);   %% linear
    iq2 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) == 3);   %% quadratic
    iq3 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) > 3);    %% cubic or greater
    if ~isempty(iq2) || ~isempty(iq1) || ~isempty(iq0)
        kqg = zeros(ng, 1);
        cqg = zeros(ng, 1);
        if ~isempty(iq2)
            Qqg = zeros(ng, 1);
            Qqg(iq2) = 2 * qcost(iq2, COST) * baseMVA^2;
            cqg(iq2) = cqg(iq2) + qcost(iq2, COST+1) * baseMVA;
            kqg(iq2) = kqg(iq2) + qcost(iq2, COST+2);
        else
            Qqg = [];   %% no quadratic terms
        end
        if ~isempty(iq1)
            cqg(iq1) = cqg(iq1) + qcost(iq1, COST) * baseMVA;
            kqg(iq1) = kqg(iq1) + qcost(iq1, COST+1);
        end
        if ~isempty(iq0)
            kqg(iq0) = kqg(iq0) + qcost(iq0, COST);
        end
    end
end

%% branch voltage angle difference limits
[Aang, lang, uang, iang]  = makeAang(baseMVA, branch, nb, mpopt);
nang = length(iang);

if dc               %% DC model
  %% check generator costs
  if ~isempty(ip3)
    error('opf_setup: DC OPF cannot handle polynomial costs with higher than quadratic order.');
  end

  %% more problem dimensions
  nv    = 0;            %% number of voltage magnitude vars
  nq    = 0;            %% number of Qg vars
  q1    = [];           %% index of 1st Qg column in Ay

  %% power mismatch constraints
  [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
  neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng);   %% Pbus w.r.t. Pg
  Amis = [B neg_Cg];
  bmis = -(bus(:, PD) + bus(:, GS)) / baseMVA - Pbusinj;

  %% branch flow constraints
  il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
  nl2 = length(il);         %% number of constrained lines
  if nl2
    upf = branch(il, RATE_A) / baseMVA - Pfinj(il);
    upt = branch(il, RATE_A) / baseMVA + Pfinj(il);
  else
    upf = [];
    upt = [];
  end

  user_vars = {'Va', 'Pg'};
  ycon_vars = {'Pg', 'y'};
else                %% AC model
  %% more problem dimensions
  nv    = nb;           %% number of voltage magnitude vars
  nq    = ng;           %% number of Qg vars
  q1    = 1+ng;         %% index of 1st Qg column in Ay

  %% find branches with flow limits
  il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
  nl2 = length(il);         %% number of constrained lines

  %% build admittance matrices
  [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);

  %% dispatchable load, constant power factor constraints
  [Avl, lvl, uvl]  = makeAvl(mpc);

  %% generator PQ capability curve constraints
  [Apqh, ubpqh, Apql, ubpql, Apqdata] = makeApq(baseMVA, gen);

  if vcart
      user_vars = {'Vr', 'Vi', 'Pg', 'Qg'};
      nodal_balance_vars = {'Vr', 'Vi', 'Pg', 'Qg'};
      flow_lim_vars = {'Vr', 'Vi'};
  else
      user_vars = {'Va', 'Vm', 'Pg', 'Qg'};
      nodal_balance_vars = {'Va', 'Vm', 'Pg', 'Qg'};
      flow_lim_vars = {'Va', 'Vm'};
  end
  ycon_vars = {'Pg', 'Qg', 'y'};

  %% nonlinear constraint functions
  if mpopt.opf.current_balance
    mis_cons = {'rImis', 'iImis'};
    fcn_mis = @(x)opf_current_balance_fcn(x, mpc, Ybus, mpopt);
    hess_mis = @(x, lam)opf_current_balance_hess(x, lam, mpc, Ybus, mpopt);
  else
    mis_cons = {'Pmis', 'Qmis'};
    fcn_mis = @(x)opf_power_balance_fcn(x, mpc, Ybus, mpopt);
    hess_mis = @(x, lam)opf_power_balance_hess(x, lam, mpc, Ybus, mpopt);
  end
  fcn_flow = @(x)opf_branch_flow_fcn(x, mpc, Yf(il, :), Yt(il, :), il, mpopt);
  hess_flow = @(x, lam)opf_branch_flow_hess(x, lam, mpc, Yf(il, :), Yt(il, :), il, mpopt);
  if vcart
    fcn_vref = @(x)opf_vref_fcn(x, mpc, refs, mpopt);
    hess_vref = @(x, lam)opf_vref_hess(x, lam, mpc, refs, mpopt);
    veq = find(mpc.bus(:, VMIN) == mpc.bus(:, VMAX));
    viq = find(mpc.bus(:, VMIN) ~= mpc.bus(:, VMAX));
    nveq = length(veq);
    nvlims = length(viq);
    if nveq
      fcn_veq = @(x)opf_veq_fcn(x, mpc, veq, mpopt);
      hess_veq = @(x, lam)opf_veq_hess(x, lam, mpc, veq, mpopt);
    end
    fcn_vlim = @(x)opf_vlim_fcn(x, mpc, viq, mpopt);
    hess_vlim = @(x, lam)opf_vlim_hess(x, lam, mpc, viq, mpopt);
    fcn_ang = @(x)opf_branch_ang_fcn(x, Aang, lang, uang);
    hess_ang = @(x, lam)opf_branch_ang_hess(x, lam, Aang, lang, uang);
  end

  %% nonlinear cost functions
  if ~isempty(ip3)
    cost_Pg = @(x)opf_gen_cost_fcn(x, baseMVA, pcost, ip3);
  end
  if ~isempty(qcost) && ~isempty(iq3)
    cost_Qg = @(x)opf_gen_cost_fcn(x, baseMVA, qcost, iq3);
  end
end

%% basin constraints for piece-wise linear gen cost variables
if (strcmp(alg, 'PDIPM') && mpopt.pdipm.step_control) || strcmp(alg, 'TRALM')
  %% SC-PDIPM or TRALM, no CCV cost vars
  ny = 0;
  Ay = sparse(0, ng+nq);
  by =[];
else
  ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
  ny = size(ipwl, 1);   %% number of piece-wise linear cost vars
  [Ay, by] = makeAy(baseMVA, ng, gencost, 1, q1, 1+ng+nq);
end
if any(gencost(:, MODEL) ~= POLYNOMIAL & gencost(:, MODEL) ~= PW_LINEAR)
    error('opf_setup: some generator cost rows have invalid MODEL value');
end

%% more problem dimensions
nx    = nb+nv + ng+nq;  %% number of standard OPF control variables
if nlin
  nz = size(mpc.A, 2) - nx; %% number of user z variables
  if nz < 0
    error('opf_setup: user supplied A matrix must have at least %d columns.', nx);
  end
else
  nz = 0;               %% number of user z variables
  if nw                 %% still need to check number of columns of N
    if size(mpc.N, 2) ~= nx;
      error('opf_setup: user supplied N matrix must have %d columns.', nx);
    end
  end
end

%% construct OPF model object
om = opf_model(mpc);
if ~isempty(pwl1)
  om.userdata.pwl1 = pwl1;
end
om.userdata.iang = iang;
if dc
  %% user data
  om.userdata.Bf = Bf;
  om.userdata.Pfinj = Pfinj;

  %% optimization variables
  om.add_var('Va', nb, Va, Val, Vau);
  om.add_var('Pg', ng, Pg, Pmin, Pmax);

  %% linear constraints
  om.add_lin_constraint('Pmis', Amis, bmis, bmis, {'Va', 'Pg'});    %% nb
  om.add_lin_constraint('Pf',  Bf(il,:), -upt, upf, {'Va'});        %% nl2
  om.add_lin_constraint('ang', Aang, lang, uang, {'Va'});           %% nang

  %% quadratic generator costs
  if ~isempty(cpg)
    om.add_quad_cost('polPg', Qpg, cpg, kpg, {'Pg'});
  end
else
  %% user data
  om.userdata.Apqdata = Apqdata;

  %% optimization variables
  if vcart
      Vclim = 1.1 * bus(:, VMAX);
      om.add_var('Vr', nb, Vr, -Vclim, Vclim);
      om.add_var('Vi', nb, Vi, -Vclim, Vclim);
  else
      om.add_var('Va', nb, Va, Val, Vau);
      om.add_var('Vm', nb, Vm, bus(:, VMIN), bus(:, VMAX));
  end
  om.add_var('Pg', ng, Pg, Pmin, Pmax);
  om.add_var('Qg', ng, Qg, Qmin, Qmax);

  %% nonlinear constraints
  om.add_nln_constraint(mis_cons, [nb;nb], 1, fcn_mis, hess_mis, nodal_balance_vars);
  if legacy_formulation
    om.add_nln_constraint({'Sf', 'St'}, [nl;nl], 0, fcn_flow, hess_flow, flow_lim_vars);
  else
    om.add_nln_constraint({'Sf', 'St'}, [nl2;nl2], 0, fcn_flow, hess_flow, flow_lim_vars);
  end
  if vcart
    om.userdata.veq = veq;  %% buses with voltage magnitude equality constraints
    om.userdata.viq = viq;  %% buses with voltage magnitude limits
    om.add_nln_constraint('Vref', length(refs), 1, fcn_vref, hess_vref, {'Vr', 'Vi'});
    if nveq
      om.add_nln_constraint('Veq', nveq, 1, fcn_veq, hess_veq, {'Vr', 'Vi'});
    end
    om.add_nln_constraint({'Vmin', 'Vmax'}, [nvlims;nvlims], 0, fcn_vlim, hess_vlim, {'Vr', 'Vi'});
    om.add_nln_constraint({'angL', 'angU'}, [nang;nang], 0, fcn_ang, hess_ang, {'Vr', 'Vi'});
  end

  %% linear constraints
  om.add_lin_constraint('PQh', Apqh, [], ubpqh, {'Pg', 'Qg'});      %% npqh
  om.add_lin_constraint('PQl', Apql, [], ubpql, {'Pg', 'Qg'});      %% npql
  om.add_lin_constraint('vl',  Avl, lvl, uvl,   {'Pg', 'Qg'});      %% nvl
  if ~vcart
    om.add_lin_constraint('ang', Aang, lang, uang, {'Va'});         %% nang
  end

  %% polynomial generator costs
  if ~legacy_formulation
    %% quadratic/linear generator costs
    if ~isempty(cpg)
      om.add_quad_cost('polPg', Qpg, cpg, kpg, {'Pg'});
    end
    if ~isempty(cqg)
      om.add_quad_cost('polQg', Qqg, cqg, kqg, {'Qg'});
    end

    %% higher order polynomial generator costs
    if ~isempty(ip3)
      om.add_nln_cost('polPg', 1, cost_Pg, {'Pg'});
    end
    if ~isempty(qcost) && ~isempty(iq3)
      om.add_nln_cost('polQg', 1, cost_Qg, {'Qg'});
    end
  end
end

%% y vars, constraints for piece-wise linear gen costs
if ny > 0
  om.add_var('y', ny);
  om.add_lin_constraint('ycon', Ay, [], by, ycon_vars);             %% ncony
  if dc || ~legacy_formulation
    om.add_quad_cost('pwl', [], ones(ny, 1), 0, {'y'});
  end
end

%% add user vars, constraints and costs (as specified via A, ..., N, ...)
if nz > 0
  om.add_var('z', nz, z0, zl, zu);
  user_vars{end+1} = 'z';
end
if nlin
  om.add_lin_constraint('usr', mpc.A, lbu, ubu, user_vars);         %% nlin
end
if nnle
  for k = 1:length(mpc.user_constraints.nle)
    nlc = mpc.user_constraints.nle{k};
    fcn  = eval(['@(x)' nlc{3} '(x, nlc{6}{:})']);
    hess = eval(['@(x, lam)' nlc{4} '(x, lam, nlc{6}{:})']);
    om.add_nln_constraint(nlc{1:2}, 1, fcn, hess, nlc{5});
  end
end
if nnli
  for k = 1:length(mpc.user_constraints.nli)
    nlc = mpc.user_constraints.nli{k};
    fcn  = eval(['@(x)' nlc{3} '(x, nlc{6}{:})']);
    hess = eval(['@(x, lam)' nlc{4} '(x, lam, nlc{6}{:})']);
    om.add_nln_constraint(nlc{1:2}, 0, fcn, hess, nlc{5});
  end
end
if nw
  user_cost.N = mpc.N;
  user_cost.Cw = Cw;
  if ~isempty(fparm)
    user_cost.dd = fparm(:, 1);
    user_cost.rh = fparm(:, 2);
    user_cost.kk = fparm(:, 3);
    user_cost.mm = fparm(:, 4);
  end
  if ~isempty(H)
    user_cost.H = H;
  end
  om.add_legacy_cost('usr', user_cost, user_vars);
end

%% execute userfcn callbacks for 'formulation' stage
om = run_userfcn(userfcn, 'formulation', om, mpopt);

%% implement legacy user costs using quadratic or general non-linear costs
cp = om.params_legacy_cost();   %% construct/fetch the parameters
[N, H, Cw, rh, mm] = deal(cp.N, cp.H, cp.Cw, cp.rh, cp.mm);
[nw, nx] = size(N);
if nw
    if any(cp.dd ~= 1) || any(cp.kk)    %% not simple quadratic form
        if dc                           %% (includes "dead zone" or
            if any(cp.dd ~= 1)          %%  quadratic "penalty")
                error('opf_setup: DC OPF can only handle legacy user-defined costs with d = 1');
            end
            if any(cp.kk)
                error('opf_setup: DC OPF can only handle legacy user-defined costs with no "dead zone", i.e. k = 0');
            end
        elseif ~legacy_formulation
            %% use general nonlinear cost to implement legacy user cost
            legacy_cost_fcn = @(x)opf_legacy_user_cost_fcn(x, cp);
            om.add_nln_cost('usr', 1, legacy_cost_fcn);
        end
    else                                %% simple quadratic form
        %% use a quadratic cost to implement legacy user cost
        if dc || ~legacy_formulation
            %% f = 1/2 * w'*H*w + Cw'*w, where w = diag(mm)*(N*x - rh)
            %% Let: MN = diag(mm)*N
            %%      MR = M * rh
            %%      HMR  = H  * MR;
            %%      HtMR = H' * MR;
            %%  =>   w = MN*x - MR
            %% f = 1/2 * (MN*x - MR)'*H*(MN*x - MR) + Cw'*(MN*x - MR)
            %%   = 1/2 * x'*MN'*H*MN*x +
            %%          (Cw'*MN - 1/2 * MR'*(H+H')*MN)*x +
            %%          1/2 * MR'*H*MR - Cw'*MR
            %%   = 1/2 * x'*Q*w + c'*x + k

            M    = sparse(1:nw, 1:nw, mm, nw, nw);
            MN   = M * N;
            MR   = M * rh;
            HMR  = H  * MR;
            HtMR = H' * MR;
            Q = MN' * H * MN;
            c = full(MN' * (Cw - 1/2*(HMR+HtMR)));
            k = (1/2 * HtMR - Cw)' * MR;
            om.add_quad_cost('usr', Q, c, k);
        end
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function confs = loadOptions(inputFile, confs)

fin = fopen(inputFile, 'rt');

for k = 1:length(confs.vars)
    if confs.args(k) < 10
        [numOut, varargout] = Get_Vals(fin, confs.vars{k}, confs.args(k));
        if numOut > 0
            strVar = sprintf('confs.%s=str2num(varargout);', confs.vars{k});
            eval(strVar);
        end
    else
        [numOut, varargout] = Get_Vals(fin, confs.vars{k}, 1);
        if numOut > 0
            strVar = sprintf('confs.%s=strtrim(varargout);', confs.vars{k});
            eval(strVar);
        end
    end
end

fclose(fin);

return;

function [numOut, vals] = Get_Vals(fInput, Keyword, NumVals)

%global userdata;

numOut = 0;
vals = '';

fseek(fInput, 0, 'bof');
delim = ' ,';

while 1
    tline = fgetl(fInput);
    if ~ischar(tline),   break,   end;
    [key, remain] = strtok(tline, delim);
    if strcmpi(strtrim(key), strtrim(Keyword))
        while numOut < NumVals
           if ~isempty(remain)
               numOut = numOut+1;
               vals = remain;
           end
        end
        break;
    end
end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% Lists m-files.
ls -lrt *.mat

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function Idenoised = nsf(I)
% noise supression function(NSF-Eq.(6)in the paper) is used to denoise the enhanced image.
warning off;
I = double(I);
h = fspecial('gaussian');
I = imfilter(I,h);
[m,n] = size(I);
[count,x] = hist(I(:),100);
[d,ind] = max(count);
t = 2;
x1 = x(1:ind+t);
count1 = count(1:ind+t);
coeffs = create_Fit(x1,count1);

for i = 1:m
    for j = 1:n
        [d,ind] = min(abs(I(i,j)-x));
        g(i,j) = shrinkage(I(i,j),ind,count,coeffs);
    end
end

Idenoised = I.* g;
end

function g = shrinkage(x,k,count,coeffs)
%             (1-w)*p(x|noise)
% g = 1 - -----------------------------
%          w*p(x|edge)+(1-w)*p(x|noise)

a1 = coeffs(1);
b1 = coeffs(2);
c1 = coeffs(3);

if x < b1
    g = 0.05;
else
    p1 = a1*exp(-((x-b1)^2/c1^2));

    p = count(k);
    g = abs(1 - (p1/p));
end
end

author: Sable

repo: mcbench-benchmarks


text:
function fcstr = directCollocation(obj, name, x, dx)
    % Return the SymFunction object of the direct collocation constraint
    % given the state (x) and derivatives (dx).
    %
    % Parameters:
    % name: the name suffix of the function @type char
    % x: the state SymVariable @type SymVariable
    % dx: the derivative of states @type SymVariable

    T  = [SymVariable('t0');SymVariable('tf')];
    Ts = T(2) - T(1);
    N = SymVariable('nNode');
    numState = length(x);
    nNode = obj.NumNode;
    switch obj.Options.CollocationScheme
        case 'HermiteSimpson'

            xn = SymVariable('xn',[numState,1]);
            dxn = SymVariable('dxn',[numState,1]);
            xm = SymVariable('xm',[numState,1]);
            dxm = SymVariable('dxm',[numState,1]);

            int_x = [xn - x - (2.*Ts./(N-1)).*(dxn + 4.*dxm + dx)./6;
                xm - (x+xn)./2 - 2.*Ts.*(dx-dxn)./(8*(N-1))];

            if isnan(obj.Options.ConstantTimeHorizon)
                fcstr = SymFunction(['hs_int_' name],int_x,{T,x,dx,xm,dxm,xn,dxn},{N});
            else
                fcstr = SymFunction(['hs_int_' name],int_x,{x,dx,xm,dxm,xn,dxn},{T,N});
            end

        case 'Trapezoidal'
            xn = SymVariable('xn',[numState,1]);
            dxn = SymVariable('dxn',[numState,1]);

            int_x = xn - x - (Ts./(N-1)).*(dxn + dx)./2;

            if isnan(obj.Options.ConstantTimeHorizon)
                fcstr = SymFunction(['tr_int_' name],int_x,{T,x,dx,xn,dxn},{N});
            else
                fcstr = SymFunction(['tr_int_' name],int_x,{x,dx,xn,dxn},{T,N});
            end

        case 'PseudoSpectral'
            t = sym('t');
            p = legendreP(nNode-1,t);
            dp = jacobian(p,t);

            roots = vpasolve(dp*(1-t)*(1+t)==0);

            D_LGL = zeros(nNode);
            for i=1:nNode
                for j=1:nNode
                    if i==j
                        if j== 1
                            D_LGL(i,j) = - (nNode-1)*(nNode)/4;
                        elseif j==nNode
                            D_LGL(i,j) = (nNode-1)*(nNode)/4;
                        else
                            D_LGL(i,j) = 0;
                        end
                    else
                        D_LGL(i,j) = subs(p,t,roots(i))/(subs(p,t,roots(j))*(roots(i) - roots(j)));
                    end
                end
            end

            KD_LGL = kron(D_LGL, eye(numState));

            xn = cell(1,nNode);
            dxn = cell(1,nNode);

            for i=1:nNode
                xn{i} = SymVariable(['x',num2str(i)],[numState,1]);
                dxn{i} = SymVariable(['dx',num2str(i)],[numState,1]);
            end

            X = transpose(flatten([xn{:}]));
            dX = transpose(flatten([dxn{:}]));
            int_x = dX.*(Ts./2) - KD_LGL*X;

            dep = [xn;dxn];
            if isnan(obj.Options.ConstantTimeHorizon)
                fcstr = SymFunction(['ps_int_' name '_' obj.Name],int_x,[{T},dep(:)']);
            else
                fcstr = SymFunction(['ps_int_' name '_' obj.Name],int_x,dep(:)',{T});
            end
    end

end


text:
function value = pyramid_unit_monomial_3d ( alpha, beta, gamma )

%*****************************************************************************80
%
%% PYRAMID_UNIT_MONOMIAL_3D: monomial integral in a unit pyramid in 3D.
%
%  Discussion:
%
%    This routine returns the integral of X^ALPHA Y^BETA Z^GAMMA over
%    the unit pyramid.
%
%    The unit pyramid is defined as:
%
%    - ( 1 - Z ) <= X <= 1 - Z
%    - ( 1 - Z ) <= Y <= 1 - Z
%              0 <= Z <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer ALPHA, BETA, GAMMA, the exponents of
%    X, Y and Z in the monomial.
%
%    Output, real PYRAMID_UNIT_MONOMIAL_3D, the volume of the pyramid.
%
  value = 0.0;

  if ( mod ( alpha, 2 ) == 0 & mod ( beta, 2 ) == 0 )

    i_hi = 2 + alpha + beta;

    for i = 0 : i_hi
      value = value + r8_mop ( i ) * r8_choose ( i_hi, i ) / ( i + gamma + 1 );
    end

    value = value * 2.0 / ( alpha + 1 ) * 2.0 / ( beta + 1 );

  end

  return
end


text:
function N = spm_ecat2nifti(fname,opts)
% Import ECAT 7 images from CTI PET scanners
% FORMAT N = spm_ecat2nifti(fname)
% fname    - name of ECAT file
% opts     - options structure
%
% N        - NIfTI object (written in current directory)
%__________________________________________________________________________
% Copyright (C) 2005-2015 Wellcome Trust Centre for Neuroimaging

% John Ashburner & Roger Gunn
% $Id: spm_ecat2nifti.m 6315 2015-01-23 17:09:06Z guillaume $

if nargin==1
    opts = struct('ext',spm_file_ext);
else
    if opts.ext(1) ~= '.', opts.ext = ['.' opts.ext]; end
end

fp = fopen(fname,'r','ieee-be');
if fp == -1,
    error(['Can''t open "' fname '".']);
end

mh = ECAT7_mheader(fp);
if ~strcmp(mh.MAGIC_NUMBER,'MATRIX70v') &&...
        ~strcmp(mh.MAGIC_NUMBER,'MATRIX71v') &&...
        ~strcmp(mh.MAGIC_NUMBER,'MATRIX72v')
    fclose(fp);
    error(['"' fname '" does not appear to be ECAT 7 format.']);
end

if mh.FILE_TYPE ~= 7
    fclose(fp);
    error(['"' fname '" does not appear to be an image file.']);
end

list    = s7_matlist(fp);
matches = find((list(:,4) == 1) | (list(:,4) == 2));
llist   = list(matches,:);

for i=1:size(llist,1)
    sh(i) = ECAT7_sheader(fp,llist(i,2));
end
fclose(fp);

for i=1:size(llist,1)
    dim      = [sh(i).X_DIMENSION sh(i).Y_DIMENSION sh(i).Z_DIMENSION];
    dtype    = [4 1];
    off      = 512*llist(i,2);
    scale    = sh(i).SCALE_FACTOR*mh.ECAT_CALIBRATION_FACTOR;
    inter    = 0;
    dati     = file_array(fname,dim,dtype,off,scale,inter);

    dircos   = diag([-1 -1 -1]);
    step     = ([sh(i).X_PIXEL_SIZE sh(i).Y_PIXEL_SIZE sh(i).Z_PIXEL_SIZE]*10);
    start    = -(dim(1:3)'/2).*step';
    mat      = [[dircos*diag(step) dircos*start] ; [0 0 0 1]];

    matnum   = sprintf('%.8x',list(i,1));
    nam      = spm_file(fname,'basename');
    fnameo   = fullfile(pwd,[nam '_' matnum opts.ext]);
    dato     = file_array(fnameo,dim,[4 spm_platform('bigend')],0,scale,inter);

    N        = nifti;
    N.dat    = dato;
    N.mat    = mat;
    N.mat0   = mat;
    N.mat_intent  = 'aligned';
    N.mat0_intent = 'scanner';
    N.descrip     = sh(i).ANNOTATION;
    N.timing = struct('toffset',sh(i).FRAME_START_TIME/1000,'tspace',sh(i).FRAME_DURATION/1000);
    create(N);
    for j=1:dim(3)
        N.dat(:,:,j) = dati(:,:,j);
    end
    N.extras     = struct('mh',mh,'sh',sh(i));
end

%==========================================================================
% function list = s7_matlist(fp)
%==========================================================================
function list = s7_matlist(fp)
%S7_MATLIST List the available matrixes in an ECAT 7 file.
%   LIST = S7_MATLIST(FP) lists the available matrixes
%   in the file specified by FP.
%
%   Columns in LIST:
%       1 - Matrix identifier.
%       2 - Matrix subheader record number
%       3 - Last record number of matrix data block.
%       4 - Matrix status:
%           1 - exists - rw
%           2 - exists - ro
%           3 - matrix deleted
%

% I believe fp should be opened with:
%     fp = fopen(filename,'r','ieee-be');

fseek(fp,512,'bof');
block = fread(fp,128,'int');
if size(block,1) ~= 128
    list = [];
    return;
end
block = reshape(block,4,32);
list  = [];
while block(2,1) ~= 2
    if block(1,1)+block(4,1) ~= 31
        list = []; return;
    end
    list = [list block(:,2:32)];

    fseek(fp,512*(block(2,1)-1),'bof');
    block = fread(fp,128,'int');
    if size(block,1) ~= 128, list = []; return; end
    block = reshape(block,4,32);
end
list = [list block(:,2:(block(4,1)+1))];
list = list';

%==========================================================================
% function SHEADER=ECAT7_sheader(fid,record)
%==========================================================================
function SHEADER=ECAT7_sheader(fid,record)
%
% Sub header read routine for ECAT 7 image files
%
% Roger Gunn, 260298
off                     = (record-1)*512;
status                  = fseek(fid, off,'bof');
data_type               = fread(fid,1,'uint16',0);
num_dimensions          = fread(fid,1,'uint16',0);
x_dimension             = fread(fid,1,'uint16',0);
y_dimension             = fread(fid,1,'uint16',0);
z_dimension             = fread(fid,1,'uint16',0);
x_offset                = fread(fid,1,'float32',0);
y_offset                = fread(fid,1,'float32',0);
z_offset                = fread(fid,1,'float32',0);
recon_zoom              = fread(fid,1,'float32',0);
scale_factor            = fread(fid,1,'float32',0);
image_min               = fread(fid,1,'int16',0);
image_max               = fread(fid,1,'int16',0);
x_pixel_size            = fread(fid,1,'float32',0);
y_pixel_size            = fread(fid,1,'float32',0);
z_pixel_size            = fread(fid,1,'float32',0);
frame_duration          = fread(fid,1,'uint32',0);
frame_start_time        = fread(fid,1,'uint32',0);
filter_code             = fread(fid,1,'uint16',0);
x_resolution            = fread(fid,1,'float32',0);
y_resolution            = fread(fid,1,'float32',0);
z_resolution            = fread(fid,1,'float32',0);
num_r_elements          = fread(fid,1,'float32',0);
num_angles              = fread(fid,1,'float32',0);
z_rotation_angle        = fread(fid,1,'float32',0);
decay_corr_fctr         = fread(fid,1,'float32',0);
corrections_applied     = fread(fid,1,'uint32',0);
gate_duration           = fread(fid,1,'uint32',0);
r_wave_offset           = fread(fid,1,'uint32',0);
num_accepted_beats      = fread(fid,1,'uint32',0);
filter_cutoff_frequency = fread(fid,1,'float32',0);
filter_resolution       = fread(fid,1,'float32',0);
filter_ramp_slope       = fread(fid,1,'float32',0);
filter_order            = fread(fid,1,'uint16',0);
filter_scatter_fraction = fread(fid,1,'float32',0);
filter_scatter_slope    = fread(fid,1,'float32',0);
annotation              = fread(fid,40,'char',0);
mt_1_1                  = fread(fid,1,'float32',0);
mt_1_2                  = fread(fid,1,'float32',0);
mt_1_3                  = fread(fid,1,'float32',0);
mt_2_1                  = fread(fid,1,'float32',0);
mt_2_2                  = fread(fid,1,'float32',0);
mt_2_3                  = fread(fid,1,'float32',0);
mt_3_1                  = fread(fid,1,'float32',0);
mt_3_2                  = fread(fid,1,'float32',0);
mt_3_3                  = fread(fid,1,'float32',0);
rfilter_cutoff          = fread(fid,1,'float32',0);
rfilter_resolution      = fread(fid,1,'float32',0);
rfilter_code            = fread(fid,1,'uint16',0);
rfilter_order           = fread(fid,1,'uint16',0);
zfilter_cutoff          = fread(fid,1,'float32',0);
zfilter_resolution      = fread(fid,1,'float32',0);
zfilter_code            = fread(fid,1,'uint16',0);
zfilter_order           = fread(fid,1,'uint16',0);
mt_4_1                  = fread(fid,1,'float32',0);
mt_4_2                  = fread(fid,1,'float32',0);
mt_4_3                  = fread(fid,1,'float32',0);
scatter_type            = fread(fid,1,'uint16',0);
recon_type              = fread(fid,1,'uint16',0);
recon_views             = fread(fid,1,'uint16',0);
fill                    = fread(fid,1,'uint16',0);
annotation              = deblank(char(annotation.*(annotation>0))');

SHEADER = struct('DATA_TYPE', data_type, ...
    'NUM_DIMENSIONS', num_dimensions, ...
    'X_DIMENSION', x_dimension, ...
    'Y_DIMENSION', y_dimension, ...
    'Z_DIMENSION', z_dimension, ...
    'X_OFFSET', x_offset, ...
    'Y_OFFSET', y_offset, ...
    'Z_OFFSET', z_offset, ...
    'RECON_ZOOM', recon_zoom, ...
    'SCALE_FACTOR', scale_factor, ...
    'IMAGE_MIN', image_min, ...
    'IMAGE_MAX', image_max, ...
    'X_PIXEL_SIZE', x_pixel_size, ...
    'Y_PIXEL_SIZE', y_pixel_size, ...
    'Z_PIXEL_SIZE', z_pixel_size, ...
    'FRAME_DURATION', frame_duration, ...
    'FRAME_START_TIME', frame_start_time, ...
    'FILTER_CODE', filter_code, ...
    'X_RESOLUTION', x_resolution, ...
    'Y_RESOLUTION', y_resolution, ...
    'Z_RESOLUTION', z_resolution, ...
    'NUM_R_ELEMENTS', num_r_elements, ...
    'NUM_ANGLES', num_angles, ...
    'Z_ROTATION_ANGLE', z_rotation_angle, ...
    'DECAY_CORR_FCTR', decay_corr_fctr, ...
    'CORRECTIONS_APPLIED', corrections_applied, ...
    'GATE_DURATION', gate_duration, ...
    'R_WAVE_OFFSET', r_wave_offset, ...
    'NUM_ACCEPTED_BEATS', num_accepted_beats, ...
    'FILTER_CUTOFF_FREQUENCY', filter_cutoff_frequency, ...
    'FILTER_RESOLUTION', filter_resolution, ...
    'FILTER_RAMP_SLOPE', filter_ramp_slope, ...
    'FILTER_ORDER', filter_order, ...
    'FILTER_SCATTER_CORRECTION', filter_scatter_fraction, ...
    'FILTER_SCATTER_SLOPE', filter_scatter_slope, ...
    'ANNOTATION', annotation, ...
    'MT_1_1', mt_1_1, ...
    'MT_1_2', mt_1_2, ...
    'MT_1_3', mt_1_3, ...
    'MT_2_1', mt_2_1, ...
    'MT_2_2', mt_2_2, ...
    'MT_2_3', mt_2_3, ...
    'MT_3_1', mt_3_1, ...
    'MT_3_2', mt_3_2, ...
    'MT_3_3', mt_3_3, ...
    'RFILTER_CUTOFF', rfilter_cutoff, ...
    'RFILTER_RESOLUTION', rfilter_resolution, ...
    'RFILTER_CODE', rfilter_code, ...
    'RFILTER_ORDER', rfilter_order, ...
    'ZFILTER_CUTOFF', zfilter_cutoff, ...
    'ZFILTER_RESOLUTION', zfilter_resolution, ...
    'ZFILTER_CODE', zfilter_code, ...
    'ZFILTER_ORDER', zfilter_order, ...
    'MT_4_1', mt_4_1, ...
    'MT_4_2', mt_4_2, ...
    'MT_4_3', mt_4_3, ...
    'SCATTER_TYPE', scatter_type, ...
    'RECON_TYPE', recon_type, ...
    'RECON_VIEWS', recon_views, ...
    'FILL', fill);

%==========================================================================
% function [MHEADER]=ECAT7_mheader(fid)
%==========================================================================
function [MHEADER]=ECAT7_mheader(fid)
%
% Main header read routine for ECAT 7 image files
%
% Roger Gunn, 260298

status                  = fseek(fid, 0,'bof');
magic_number            = fread(fid,14,'char',0);
original_file_name      = fread(fid,32,'char',0);
sw_version              = fread(fid,1,'uint16',0);
system_type             = fread(fid,1,'uint16',0);
file_type               = fread(fid,1,'uint16',0);
serial_number           = fread(fid,10,'char',0);
scan_start_time         = fread(fid,1,'uint32',0);
isotope_name            = fread(fid,8,'char',0);
isotope_halflife        = fread(fid,1,'float32',0);
radiopharmaceutical     = fread(fid,32,'char',0);
gantry_tilt             = fread(fid,1,'float32',0);
gantry_rotation         = fread(fid,1,'float32',0);
bed_elevation           = fread(fid,1,'float32',0);
intrinsic_tilt           = fread(fid,1,'float32',0);
wobble_speed             = fread(fid,1,'uint16',0);
transm_source_type       = fread(fid,1,'uint16',0);
distance_scanned         = fread(fid,1,'float32',0);
transaxial_fov           = fread(fid,1,'float32',0);
angular_compression      = fread(fid,1,'uint16',0);
coin_samp_mode           = fread(fid,1,'uint16',0);
axial_samp_mode          = fread(fid,1,'uint16',0);
ecat_calibration_factor  = fread(fid,1,'float32',0);
calibration_units        = fread(fid,1,'uint16',0);
calibration_units_type   = fread(fid,1,'uint16',0);
compression_code         = fread(fid,1,'uint16',0);
study_type               = fread(fid,12,'char',0);
patient_id               = fread(fid,16,'char',0);
patient_name             = fread(fid,32,'char',0);
patient_sex              = fread(fid,1,'char',0);
patient_dexterity        = fread(fid,1,'char',0);
patient_age              = fread(fid,1,'float32',0);
patient_height           = fread(fid,1,'float32',0);
patient_weight           = fread(fid,1,'float32',0);
patient_birth_date       = fread(fid,1,'uint32',0);
physician_name           = fread(fid,32,'char',0);
operator_name            = fread(fid,32,'char',0);
study_description        = fread(fid,32,'char',0);
acquisition_type         = fread(fid,1,'uint16',0);
patient_orientation      = fread(fid,1,'uint16',0);
facility_name            = fread(fid,20,'char',0);
num_planes               = fread(fid,1,'uint16',0);
num_frames               = fread(fid,1,'uint16',0);
num_gates                = fread(fid,1,'uint16',0);
num_bed_pos              = fread(fid,1,'uint16',0);
init_bed_position        = fread(fid,1,'float32',0);
bed_position = zeros(15,1);
for bed=1:15
    tmp = fread(fid,1,'float32',0);
    if ~isempty(tmp), bed_position(bed) = tmp; end
end
plane_separation         = fread(fid,1,'float32',0);
lwr_sctr_thres           = fread(fid,1,'uint16',0);
lwr_true_thres           = fread(fid,1,'uint16',0);
upr_true_thres           = fread(fid,1,'uint16',0);
user_process_code        = fread(fid,10,'char',0);
acquisition_mode         = fread(fid,1,'uint16',0);
bin_size                 = fread(fid,1,'float32',0);
branching_fraction       = fread(fid,1,'float32',0);
dose_start_time          = fread(fid,1,'uint32',0);
dosage                   = fread(fid,1,'float32',0);
well_counter_corr_factor = fread(fid,1,'float32',0);
data_units               = fread(fid,32,'char',0);
septa_state              = fread(fid,1,'uint16',0);
fill                     = fread(fid,1,'uint16',0);

magic_number             = deblank(char(magic_number.*(magic_number>32))');
original_file_name       = deblank(char(original_file_name.*(original_file_name>0))');
serial_number            = deblank(char(serial_number.*(serial_number>0))');
isotope_name             = deblank(char(isotope_name.*(isotope_name>0))');
radiopharmaceutical      = deblank(char(radiopharmaceutical.*(radiopharmaceutical>0))');
study_type               = deblank(char(study_type.*(study_type>0))');
patient_id               = deblank(char(patient_id.*(patient_id>0))');
patient_name             = deblank(char(patient_name.*(patient_name>0))');
patient_sex              = deblank(char(patient_sex.*(patient_sex>0))');
patient_dexterity        = deblank(char(patient_dexterity.*(patient_dexterity>0))');
physician_name           = deblank(char(physician_name.*(physician_name>0))');
operator_name            = deblank(char(operator_name.*(operator_name>0))');
study_description        = deblank(char(study_description.*(study_description>0))');
facility_name            = deblank(char(facility_name.*(facility_name>0))');
user_process_code        = deblank(char(user_process_code.*(user_process_code>0))');
data_units               = deblank(char(data_units.*(data_units>0))');

MHEADER = struct('MAGIC_NUMBER', magic_number, ...
    'ORIGINAL_FILE_NAME', original_file_name, ...
    'SW_VERSION', sw_version, ...
    'SYSTEM_TYPE', system_type, ...
    'FILE_TYPE', file_type, ...
    'SERIAL_NUMBER', serial_number, ...
    'SCAN_START_TIME', scan_start_time, ...
    'ISOTOPE_NAME', isotope_name, ...
    'ISOTOPE_HALFLIFE', isotope_halflife, ...
    'RADIOPHARMACEUTICAL', radiopharmaceutical, ...
    'GANTRY_TILT', gantry_tilt, ...
    'GANTRY_ROTATION', gantry_rotation, ...
    'BED_ELEVATION', bed_elevation, ...
    'INTRINSIC_TILT', intrinsic_tilt, ...
    'WOBBLE_SPEED', wobble_speed, ...
    'TRANSM_SOURCE_TYPE', transm_source_type, ...
    'DISTANCE_SCANNED', distance_scanned, ...
    'TRANSAXIAL_FOV', transaxial_fov, ...
    'ANGULAR_COMPRESSION', angular_compression, ...
    'COIN_SAMP_MODE', coin_samp_mode, ...
    'AXIAL_SAMP_MODE', axial_samp_mode, ...
    'ECAT_CALIBRATION_FACTOR', ecat_calibration_factor, ...
    'CALIBRATION_UNITS', calibration_units, ...
    'CALIBRATION_UNITS_TYPE', calibration_units_type, ...
    'COMPRESSION_CODE', compression_code, ...
    'STUDY_TYPE', study_type, ...
    'PATIENT_ID', patient_id, ...
    'PATIENT_NAME', patient_name, ...
    'PATIENT_SEX', patient_sex, ...
    'PATIENT_DEXTERITY', patient_dexterity, ...
    'PATIENT_AGE', patient_age, ...
    'PATIENT_HEIGHT', patient_height, ...
    'PATIENT_WEIGHT', patient_weight, ...
    'PATIENT_BIRTH_DATE', patient_birth_date, ...
    'PHYSICIAN_NAME', physician_name, ...
    'OPERATOR_NAME', operator_name, ...
    'STUDY_DESCRIPTION', study_description, ...
    'ACQUISITION_TYPE', acquisition_type, ...
    'PATIENT_ORIENTATION', patient_orientation, ...
    'FACILITY_NAME', facility_name, ...
    'NUM_PLANES', num_planes, ...
    'NUM_FRAMES', num_frames, ...
    'NUM_GATES', num_gates, ...
    'NUM_BED_POS', num_bed_pos, ...
    'INIT_BED_POSITION', init_bed_position, ...
    'BED_POSITION', bed_position, ...
    'PLANE_SEPARATION', plane_separation, ...
    'LWR_SCTR_THRES', lwr_sctr_thres, ...
    'LWR_TRUE_THRES', lwr_true_thres, ...
    'UPR_TRUE_THRES', upr_true_thres, ...
    'USER_PROCESS_CODE', user_process_code, ...
    'ACQUISITION_MODE', acquisition_mode, ...
    'BIN_SIZE', bin_size, ...
    'BRANCHING_FRACTION', branching_fraction, ...
    'DOSE_START_TIME', dose_start_time, ...
    'DOSAGE', dosage, ...
    'WELL_COUNTER_CORR_FACTOR', well_counter_corr_factor, ...
    'DATA_UNITS', data_units, ...
    'SEPTA_STATE', septa_state, ...
    'FILL', fill);


text:
function value = bmi_metric ( w_kg, h_m )

%*****************************************************************************80
%
%% BMI_METRIC computes the body mass index given metric measurements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real W_KG, the body weight in kilograms.
%
%    Input, real H_M, the body height in meters.
%
%    Output, real VALUE, the body mass index.
%
  value = w_kg / ( h_m * h_m );

  return
end


text:
function [varargout]=constrainedDelaunayTetGen(V,C)

%% Create TetGen input structure

if ~isempty(C)
    inputStruct.Faces=C; %Add face constraints if not empty
    inputStruct.stringOpt='-pQY';
else
    inputStruct.stringOpt='-Q';
end
inputStruct.Nodes=V;
inputStruct.holePoints=[];
inputStruct.faceBoundaryMarker=ones(size(C,1),1); %Face boundary markers
inputStruct.regionPoints=[]; %region points

%% Run TetGen
[meshOutput]=runTetGen(inputStruct); %Run tetGen

TR = triangulation(meshOutput.elements,meshOutput.nodes);

switch nargout
    case 1
        varargout{1}=TR;
    case 2
        varargout{1}=TR;
        varargout{2}=meshOutput;
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