text:
function y = sinspace(d1, d2, n, factor)
%SINSPACE cosine spaced vector.
%   SINSPACE(X1, X2) generates a row vector of 100 sine spaced points
%   between X1 and X2.
%
%   SINSPACE(X1, X2, N) generates N points between X1 and X2.
%
%   A sine spaced vector clusters the elements toward X2:
%    X1 |      |      |      |     |    |   |  | || X2
%
%   Make n negative to cluster the elements toward X1:
%    X1 || |  |   |    |     |      |      |      | X2
%
%   For -2 < N < 2, SINSPACE returns X2.
%
%   SINSPACE(X1, X2, N, W) clusters the elements to a lesser degree as
%   dictated by W. W = 0 returns a normal sine spaced vector. W = 1 is the
%   same as LINSPACE(X1, X2, N). Experiment with W < 0 and W > 1 for
%   different clustering patterns.
%
%   Author:     Sky Sartorius
%
%   See also COSSPACE, LINSPACE, LOGSPACE.

if nargin == 2
    n = 100;
end
if nargin < 4
    factor = false;
end

if n<0
    n = floor(double(-n));
    y = d2 + (d1-d2)*sin(pi/2*(1-(0:(n-1))/(n-1)));
else
    n = floor(double(n));
    y = d1 + (d2-d1)*sin(pi/(2*(n-1))*(0:n-1));
end

if factor
    y = (1-factor)*y+factor*[d1+(0:n-2)*(d2-d1)/(n-1) d2];
end

author: Sable

repo: mcbench-benchmarks


text:
function [ u, v, p ] = uvp_spiral ( nu, rho, n, x, y, t )

%*****************************************************************************80
%
%% UVP_SPIRAL returns velocity and pressure for the spiral flow.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Maxim Olshanskii, Leo Rebholz,
%    Application of barycenter refined meshes in linear elasticity
%    and incompressible fluid dynamics,
%    ETNA: Electronic Transactions in Numerical Analysis,
%    Volume 38, pages 258-274, 2011.
%
%  Parameters:
%
%    Input, real NU, the kinematic viscosity.
%
%    Input, real RHO, the fluid density.
%
%    Input, integer N, the number of nodes.
%
%    Input, real X(N), Y(N), the coordinates of nodes.
%
%    Input, real T, the current time.
%
%    Output, real U(N), V(N), the X and Y velocity.
%
%    Output, real P(N), the pressure.
%
  u = ( 1.0 + nu * t ) * 2.0 ...
    * x.^2 .* ( x - 1.0 ).^2 ...
    .* y .* ( 2.0 * y - 1.0 ) .* ( y - 1.0 );

  v = - ( 1.0 + nu * t ) * 2.0 ...
    * x .* ( 2.0 * x - 1.0 ) .* ( x - 1.0 ) ...
    .* y.^2 .* ( y - 1.0 ).^2;

  p = rho * y;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [crop,crop_norm,crop_mean,crop_std] = corner2image(img, p, tsize)
%   (r1,c1) ***** (r3,c3)            (1,1) ***** (1,cols)
%     *             *                  *           *
%      *             *       ----->     *           *
%       *             *                  *           *
%     (r2,c2) ***** (r4,c4)              (rows,1) **** (rows,cols)
afnv_obj = corners2affine(p, tsize);
map_afnv = afnv_obj.afnv;

img_map = IMGaffine_c(double(img), map_afnv, tsize);

[crop, crop_mean, crop_std] = whitening( reshape(img_map, prod([tsize 3]), 1) ); % crop is a vector
crop_norm = norm(crop);
crop = crop/crop_norm;

author: yuxng


text:
%COMPOSERT  Combines two rotation-and-shift transformations
%
%     S = cv.composeRT(rvec1, tvec1, rvec2, tvec2)
%
% ## Input
% * __rvec1__ First rotation vector, 3x1 float vector.
% * __tvec1__ First translation vector, 3x1 float vector.
% * __rvec2__ Second rotation vector, 3x1 float vector.
% * __tvec2__ Second translation vector, 3x1 float vector.
%
% ## Output
% * __S__ A scalar struct with the following fields:
%   * __rvec3__ Rotation vector of the superposition, 3x1 float vector.
%   * __tvec3__ Translation vector of the superposition, 3x1 float vector.
%   * __dr3dr1__, __dr3dt1__, __dr3dr2__, __dr3dt2__, __dt3dr1__, __dt3dt1__,
%     __dt3dr2__, __dt3dt2__ Derivatives of `rvec3` or `tvec3` with regard to
%     `rvec1`, `tvec1`, `rvec2`, and `tvec2`, respectively. Each derivative is
%     a 3x3 float matrix.
%
% The function computes:
%
%     rvec3 = rodrigues^-1( rodorigues(rvec2) * rodrigues(rvec1) )
%     tvec3 = rodrigues(rvec2) * tvec1 + tvec2
%
% where `rodrigues` denotes a rotation vector to a rotation matrix
% transformation, and `rodrigues^-1` denotes the inverse transformation.
% See cv.Rodrigues for details.
%
% Also, the function can compute the derivatives of the output vectors with
% regards to the input vectors (see cv.matMulDeriv). The function us used
% inside cv.stereoCalibrate but can also be used in your own code where
% Levenberg-Marquardt or another gradient-based solver is used to optimize a
% function that contains a matrix multiplication.
%
% See also: cv.Rodrigues, cv.matMulDeriv, cv.stereoCalibrate
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [sol,fit,popo,fits] = gaevolve(varargin)
% GAEVOLVE  Evolution of a population
%
% function [sol,fit,popo] = gaevolve(pop,funcs,fparams,params,maxiters)
%
%  This function controls the evolution of a population via simple
% genetical crossing, mutations and selection of the units. The initial
% population is given and the functions to compute the fit value, the
% crossing and the mutation are given. The maximum number of iterations is
% defined. An hooking function can be defined to hook the single evolution
% iteration step before update.
%
%  Params
%  ------
% IN:
%  pop      = The initial population: a (column) cell of units.
%  funcs    =  A cell array containing the function handlers for the
%             following functions:
%       ffit  = The fit function: recives a unit and compute its fit value.
%       fcros =  The crossing function: recives two units and cross them
%               generating a new unit.
%       fmut  =  The mutation function: recives a unit and mutate it into a
%               new unit.
%       fhook = The hooking function called before iterations: recives the
%              actual output values; returns a boolean value that says if
%              the updating process must stop.
%  fparams  = The cell of function aditional parameters for the funcs.
%  params   =  An array containing the parameters the following of the
%             updating process:
%       psel  = The selection percentage.
%       prmut = The mutation probability.
%  maxiters = The maximum number of iterations. (def=1000)
% OUT:
%  sol      = The final solution unit.
%  fit      = The final fit.
%  popo     = The final population.
%  fits     = All the actual population fits.

% Parsing of params:
[pop,ffit,fcros,fmut,fhook,pfit,pcros,pmut,phook,psel,prmut,maxiters] = ...
    ParseParams(varargin{:});

% Init:
iters = 0;
szpop = numel(pop);
nsel = round(psel*szpop);

% Computing the initial solution:
[sol,fit,fits] = ComputeSolution(pop,ffit,pfit);
popo = pop;

% Iterating:
while iters<maxiters && ~feval(fhook,sol,fit,pop,fits,phook{:})

    % ------------   Init    ------------

    % The new population:
    popo = cell(szpop,1);

    % The actual vacant element:
    actual = 1;

    % ------------ Selection ------------

    % Selection on the fits:
    [sfits,sinds] = sort(fits);
    popo(1:nsel) = pop(sinds(szpop-nsel+1:szpop));
    actual = nsel + 1;

    % ------------ Crossing  ------------

    % Computing the probabilities:
    totFit = sum(fits);
    if totFit==0
        probs = ones([1,szpop])/szpop;
    else
        probs = fits/totFit;
    end

    % The cumulative probs:
    for ind=2:szpop
        probs(ind) = probs(ind-1) + probs(ind);
    end

    % Iterating on all the other units:
    while actual<=szpop
        % Selecting the first unit:
        u1p = rand; u1 = 1;
        while probs(u1)<u1p && u1<szpop u1 = u1+1; end

        % Selecting the second unit:
        u2 = u1;
        while u1==u2 && szpop>1
            % A random unit selection:
            u2p = rand; u2 = 1;
            while probs(u2)<u2p && u2<szpop u2 = u2+1; end
        end

        % Generating a new unit:
        popo{actual} = feval(fcros,pop{u1},pop{u2},pcros{:});

        % Next:
        actual = actual + 1;
    end

    % ------------ Mutation  ------------

    % Iterating on all the units:
    for ind=1:szpop
        % Must be mutated?
        if rand<=prmut
            % Mutating:
            popo{ind} = feval(fmut,popo{ind},pmut{:});
        end
    end

    % ------------  Saving   ------------

    % The population:
    pop = popo;

    % Computing the solution:
    [sol,fit,fits] = ComputeSolution(pop,ffit,pfit);

    % A new iteration:
    iters = iters + 1;
end

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

% Computing the solution given a population:
function [sol,fit,fits] = ComputeSolution(pop,ffit,pfit)

% Computing the fits:
nel = numel(pop);
fits = zeros(1,nel);
for ind=1:nel
    % Computing a fit:
    fits(ind) = feval(ffit,pop{ind},pfit{:});
end

% Localizing and extracting the solution:
[fit,pos] = max(fits);
sol = pop(pos);

% ------------------------ LOCAL FUNCTIONS ------------------------

% Parsing of parameters:
function [pop,ffit,fcros,fmut,fhook, ...
          pfit,pcros,pmut,phook,psel,prmut,maxiters] = ...
    ParseParams(varargin)

% Minimum number of parameters is 1 (the image)
if nargin<4 error('Too few parameters'); end
if nargin>5 error('Too many parameters'); end

% Set up variables
varnames = {'pop','funcs','fparams','params','maxiters'};
for ind=1:nargin
    eval([varnames{ind} ' = varargin{ind} ;']);
end

% Default parameters
if nargin<5 maxiters=1000; end

% Checking types:
if ~isa(funcs,'cell') || numel(funcs)>4 || numel(funcs)<3
    error('The functions array must contain 3 or 4 function_handles!');
end
if ~isa(fparams,'cell') || size(fparams,1)>4
    error('The function parameter cell must contain at least 4 rows!');
end
if ~isa(params,'double') || numel(params)~=2
    error('The double parameters array must contain psel and pmut values!');
end

% Extracting the functions:
ffit = funcs{1};
fcros = funcs{2};
fmut = funcs{3};
if numel(funcs)>3
    fhook = funcs{4};
else
    fhook = @(sol,fit,popo)(false);
end

% Extracting function params:
npar = size(fparams,2);
if npar>=1 pfit=fparams{1}; else pfit={}; end
if npar>=2 pcros=fparams{2}; else pcros={}; end
if npar>=3 pmut=fparams{3}; else pmut={}; end
if npar>=4 phook=fparams{4}; else phook={}; end

% Extracting parameters:
psel = params(1);
prmut = params(2);


text:
% STE_SWLP Stabilized weighted linear prediction using short-time-energy
% weighting
%   [A,w] = ste_swlp(s,p,m,k)
%
% Description
%   This function fits linear prediction coefficients to the analysis
%   frame using the basic form of stabilized weighted linear prediction
%   (SWLP) that weights each value of the squared prediction error by the
%   short-time energy (STE) of the previous samples.
%
% Inputs
%   s      : Speech signal frame [samples]
%   p      : Order of SWLP analysis
%   m      : Length of the STE window (default to p)
%   k      : Delay of the STE window (default to 1)
%
% Outputs
%   A      : Linear prediction inverse filter coefficients
%   w      : the STE weighting function
%
% Notes
%   SWLP generally gives smoother spectral shapes than the corresponding
%   weighted linear prediction (WLP) using the same weighting function.
%
% Example
%   A = ste_swlp(s,p) gives the linear predictive inverse filter
%   coefficients optimized using STE-SWLP
%
% References
%  [1] C. Magi, J. Pohjalainen, T. Bäckström and P. Alku, "Stabilised
%  weighted linear prediction", Speech Communication, vol. 51, no. 5, pp.
%  401411, 2009.
%  [2] J. Pohjalainen, H. Kallasjoki, K. J. Palomäki, M. Kurimo and P.
%  Alku, "Weighted linear prediction for speech analysis in noisy
%  conditions", in Proc. Interspeech, Brighton, UK, pp. 1315-1318,
%  2009.
%
% Copyright (c) 2013 Aalto University
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file 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 Lesser General Public License for more
%  details.
%
% This function is part of the Common Speech Processing Repository
% http://covarep.github.io/covarep/
%
% Octave Compatible
%
% Author
%  Jouni Pohjalainen jouni.pohjalainen@aalto.fi
%
% $Id <info set by the versioning system> $

function [A,w] = ste_swlp(s,p,m,k)

% handle the special case of all-zero frames
if all(s==0)
  s = eps*randn(size(s));
end

N = length(s);

% default value for the STE window length
if nargin<3
    m = p;
end

% default value for the STE window lag
if nargin<4
    k = 1;
end

% compute STE weighting function
w = stew(s,m,p,k)+eps;

w = w(1:(N+p));

% initialize partial weights for recursive computation (see [2], Eqs.
% (5)-(7), for the partial-weight formulation of SWLP)
Z = zeros(N+p,p+1);
Z(:,1) = sqrt(w);

% delayed and weighted versions of the signal
Y = zeros(N+p,p+1);
Y(:,1) = Z(:,1).*[s;zeros(p,1)];

% recursion for partial weights and weighting of differently lagged
% versions of the signal
for i1=1:p
    Z((i1+1):(N+p),i1+1) = max([ones(N+p-i1,1) sqrt(w((i1+1):(N+p))./w(i1:(N+p-1)))],[],2) .* Z(i1:(N+p-1),i1);
    Y(:,i1+1) = Z(:,i1+1) .* [zeros(i1,1);s;zeros(p-i1,1)];
end

% compute weighted autocorrelations
R = (Y'*Y)/N;

% solve the p predictor coefficients
A = R(2:end,2:end)\R(2:end,1);

% convert to inverse filter form
A = [1;-A]';

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

function w = stew(x,m,p,k)
% STE weighting for an order p predictor using window of length m delayed
% by k samples
w = conv([zeros(k,1);x(:)].^2,ones(m,1));
w = [w;zeros(p-m,1)];

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [tArchive,GbestRBF] = UpdateRBF(Problem,tArchive,SwarmRBF,GbestRBF)
% Update solutions in RBF-assisted swarm

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    [N,D] = size(SwarmRBF);
    D     = D - 1;
    [value,best] = min(SwarmRBF(:,1+D));

    new = Problem.Evaluation(SwarmRBF(best,1:D));
    SwarmRBF(best,D+1) = new.objs;
    if new.objs < GbestRBF(:,D+1)
        GbestRBF = SwarmRBF(best,:);
    end
    tArchive = [tArchive,new];

end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% GETPROBLEM  Script for locating, reading and convert to SeDuMi, a problem
%   with the name specified by the variable `pname'. You have to define the
%   variable 'pname' before using this script.
%
% Usage: pname = 'problemname'; getproblem
%   Problem returned in (At,b,c,K).
%
% IMPORTANT
%   1. If you have an environment variable named SDPPACK, containing
%   the directory name of SDPPACK, then this searches in SDPPACK/problems.
%   2. If you have an environment variable named SDPLIB, then it searches
%   for problems in the directory SDPLIB/. These problems should be in
%   sparse SDPA format, with extension ".dat-s".
%   3. If LIPSOL is installed as a MATLAB Toolbox, then LIPSOL's  `findprob'
%   is used for searching MPS-type models.
%   4. Files with extension ".gz" are also recognized, so that you can
%   store your problem libraries more efficiently.
%
% See also  fromsdpa, frompack, prelp, feasreal, feascpx, sedumi,
%   postprocess (LIPSOL).

%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA%

global PROBDIR

% ----------------------------------------
% status = {0 don't know yet, 1 SDPPack, 2 MPS, 3 SDPA}
% ----------------------------------------
status = 0;
PROBDIR = str2mat([]);
% ----------------------------------------
% Try to get problem in SDPA format, from SDPLIB
% ----------------------------------------
if status == 0
    SDPAPATH = getenv('SDPLIB');
    if ~isempty(SDPAPATH)
        MATNAME = [SDPAPATH '/' pname '.dat-s'];
        if exist(MATNAME,'file')
            [At,b,c,K] = fromsdpa(MATNAME); status = 3;
        elseif exist([MATNAME '.gz'],'file')
            unix(['gunzip -c ' MATNAME '.gz > /tmp/default.mat']);
            [At,b,c,K] = fromsdpa('/tmp/default.mat'); status = 3;
        end
    end
    clear SDPAPATH
end
% ----------------------------------------
% Try to get problem in SDPPACK format
% ----------------------------------------
if status == 0
    PACKPATH = getenv('SDPPACK');
    if ~isempty(PACKPATH)
        MATNAME = [PACKPATH '/problems/' pname '.mat'];
        if exist(MATNAME,'file')
            unix(['/bin/cp ' MATNAME ' /tmp/default.mat']); status = 1;
        elseif exist([MATNAME '.gz'],'file')
            unix(['gunzip -c ' MATNAME '.gz > /tmp/default.mat']); status = 1;
        end
    end
    clear PACKPATH
end
% ----------------------------------------
% Try to get problem in LIPSOL format
% ----------------------------------------
if(status == 0)
    if(exist('findprob','file') > 1)
        if (findprob([],pname) == 1)
            status = 2;
        end
    end
end
% ----------------------------------------
% Move the problem into [At,b,c,K] structure
% ----------------------------------------
if(status==0)
    error('Problem not found.  Check spelling?')
elseif(status==1)
    load /tmp/default
    % ----------------------------------------
    % In SDPPACK, 0s are sometimes used  as emptyset and vice versa.
    % ----------------------------------------
    if isempty(blk.l)
        blk.l = 0;
    end
    if ~isempty(blk.q)
        if(blk.q == 0)
            blk.q = [];
        end
    end
    if ~isempty(blk.s)
        if(blk.s == 0)
            blk.s = [];
        end
    end
    [At,c] = frompack(A,b,C,blk);
    if isempty(blk.l)
        blk.l = 0;
    end
    K = blk;
    % ----------------------------------------
    % LP problem: preprocess using LIPSOL.
    % ----------------------------------------
elseif(status==2)
    [At,b,c,nx,lbs] = prelp([]);
    At = At';
    clear K;
    K.l = length(c);
end
clear status
% clear global PROBDIR   %would disable LIPSOL.

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function test_suite = test_mdwt
initTestSuite;

function test_mdwt_1D
  x = makesig('LinChirp', 8);
  h = daubcqf(4, 'min');
  L = 2;  % For 8 values in x we would normally be L=2
  [y, L] = mdwt(x, h, L);
  y_corr = [1.1097 0.8767 0.8204 -0.5201 -0.0339 0.1001 0.2201 -0.1401];
  L_corr = 2;
assertVectorsAlmostEqual(y, y_corr, 'relative', 0.001);
assertEqual(L, L_corr);

function test_mdwt_2D
  x = [1 2 3 4; 5 6 7 8 ; 9 10 11 12; 13 14 15 16];
  h = daubcqf(4);
  y = mdwt(x, h);
  y_corr = [34.0000 -3.4641 0.0000 -2.0000; -13.8564 0.0000 0.0000 -2.0000; -0.0000 0.0000 -0.0000 -0.0000; -8.0000 -8.0000 0.0000 -0.0000];
assertVectorsAlmostEqual(y, y_corr, 'relative', 0.001);

function test_mdwt_compute_L1
  x = [1 2];
  h = daubcqf(4, 'min');
  [y, L] = mdwt(x, h);
assertEqual(L, 1);

function test_mdwt_compute_L2
  x = [1 2 3 4];
  h = daubcqf(4, 'min');
  [y, L] = mdwt(x, h);
assertEqual(L, 2);

function test_mdwt_compute_L3
  x = [1 2 3 4 5 6 7 8];
  h = daubcqf(4, 'min');
  [y, L] = mdwt(x, h);
assertEqual(L, 3);

function test_mdwt_compute_bad_L
  L = -1;
  x = [1 2 3 4 5 6 7 8 9];
  h = daubcqf(4, 'min');
  mdwtHandle = @() mdwt(x, h);
assertExceptionThrown(mdwtHandle, '');

function test_mdwt_empty_input
  mdwtHandle = @() mdwt([], [0 0 0 0]);
assertExceptionThrown(mdwtHandle, '');

author: ricedsp


text:
% clear all;
addpath('../msckf/utils');
load('../datasets/dataset3.mat')

rng(42);

%Set up appropriate structs
calibParams.c_u = cu;
calibParams.c_v = cv;
calibParams.f_u = fu;
calibParams.f_v = fv;
calibParams.b = b;

%Noise parameters
y_var = [1,1,1,1];
mu = zeros(4,1);
sigma = diag(y_var);

%Generate new landmarks
newLmNum = 500;

if exist('newLmPos', 'var')
    oldLmNum = size(newLmPos,2);
    oldLmPos = newLmPos;
else
    oldLmNum = 0;
    oldLmPos = [];
end

% xRange = linspace(min(rho_i_pj_i(1,:)), max(rho_i_pj_i(1,:)), newLmNum);
% yRange = linspace(min(rho_i_pj_i(2,:)), max(rho_i_pj_i(2,:)), newLmNum);
% zRange = linspace(min(rho_i_pj_i(3,:)), max(rho_i_pj_i(3,:)), newLmNum);
% [newLmPosX, newLmPosY, newLmPosZ ] = meshgrid(xRange, yRange, zRange);

xRange = [min(rho_i_pj_i(1,:)) - 5, max(rho_i_pj_i(1,:)) + 5];
yRange = [min(rho_i_pj_i(2,:)) - 5, max(rho_i_pj_i(2,:)) + 5];
zRange = [min(rho_i_pj_i(3,:)) - 5, max(rho_i_pj_i(3,:)) + 5];

newLmPosX = range(xRange)*rand(1,newLmNum - oldLmNum) + xRange(1);
newLmPosY = range(yRange)*rand(1,newLmNum - oldLmNum) + yRange(1);
newLmPosZ = range(zRange)*rand(1,newLmNum - oldLmNum) + zRange(1);

newLmPos = [newLmPosX(:)'; newLmPosY(:)'; newLmPosZ(:)'];
newLmPos = [oldLmPos, newLmPos];

rho_i_pj_i = newLmPos;
T_cv = [C_c_v -C_c_v*rho_v_c_v; 0 0 0 1];
y_k_j = [];

for k = 1:length(t)
    C_vi = axisAngleToRotMat(theta_vk_i(:,k));
    T_vi = [C_vi -C_vi*r_i_vk_i(:,k); 0 0 0 1];
    T_ci = T_cv*T_vi;

    %Add observations
    for lm_i = 1:size(newLmPos, 2)
        p_li_i = newLmPos(:,lm_i);
        p_lc_c = homo2cart(T_ci*cart2homo(p_li_i));
        [yMeas] = stereoCamProject(p_lc_c, calibParams);
        if all(yMeas > 0) ...
            && all(yMeas([1,3]) <= 640) && all(yMeas([2,4]) <= 480)...
            && p_lc_c(3) > 0    % in view and in front of camera
                noise = mvnrnd(mu, sigma)';
                y_k_j(:,k, lm_i) = yMeas + noise;
        else
                y_k_j(:,k, lm_i) = -1*ones(4,1);
        end
    end
end

save(['../datasets/dataset3_fresh2_',num2str(newLmNum),'lessnoisy.mat']);

repo: MSCKF


text:
function circular_arc ( xcenter, ycenter, radius, angmin, angmax )

%*****************************************************************************80
%
%% CIRCULAR_ARC draws a circular arc of a given angular size and radius.
%
%  Discussion:
%
%    It is assumed that a plot has already been begun.
%
%  Example:
%
%    r = rand ( 100, 1 );
%    t = pi * rand ( 100, 1 );
%    plot ( r .* cos ( t ), r .* sin ( t ), 'b*' );
%    axis equal
%    axis square
%    circular_arc ( 0.0, 0.0, 1.0, 0, 180 );
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XCENTER, YCENTER, the X and Y coordinates of
%    the center of the circle on which the arc lies.
%
%    Input, real RADIUS, the radius of the circle on which the arc lies.
%
%    Input, real ANGMIN, ANGMAX, the minimum and maximum
%    angles of the circle.  ANGMIN and ANGMAX are both
%    measured in degrees.  ANGMIN and ANGMAX determine the
%    portion of the circle to be drawn.  If ANGMIN=0.0 and
%    ANGMAX=90.0, for instance, an arc of 90 degrees will be drawn.
%
  color = [ 0.75, 0.75, 0.75 ];

  nval = 65;

  if ( radius == 0.0 )
    return
  end
%
%  Set up the data defining the circular arc, using NVAL equally
%  spaced points along the circumference.
%
  if ( nval == 1 )
    angle(1) = 0.5 * ( angmax + angmin );
  else
    angle(1:nval) = linspace ( angmin, angmax, nval );
  end

  xval(1:nval) = xcenter + radius * cos ( angle * pi / 180 );
  yval(1:nval) = ycenter + radius * sin ( angle * pi / 180 );

  line ( xval, yval, 'Color', color );

  return
end


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nickabattista@gmail.com
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs, non-invariant beams*)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nickabattista@gmail.com) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the spring attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info)

%springs_info: col 1: starting spring pt (by lag. discretization)
%              col 2: ending spring pt. (by lag. discretization)
%              col 3: spring stiffness
%              col 4: spring resting lengths

%RL = springs_info(:,4); % resting-length vector

freq = 0.875;
%theta=freq*current_time*2*pi;
%phase = sin(theta);

%size(springs_info)

%springs_info(1:3,:)
%pause();

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL
%springs_info(319:end,4) = phase;

springs_info(319:end,4) = abs( sin(0.875*current_time*2*pi) );

% if ( mod(current_time,0.875) < 0.875/2 )
%     springs_info(319:end,4) = abs( sin(0.875*current_time*2*pi) );
% else
%     springs_info(319:end,4) = 0;%abs( sin(0.875*current_time*2*pi) );
% end

author: nickabattista


text:
function removeFolderContents(folder)
warning('off','MATLAB:RMDIR:RemovedFromPath')
if isemptyFolder(folder)
    return
end
o = dir(folder);
for ii = 1:length(o)
    if strcmp(o(ii).name,'.')
        continue
    end
    if strcmp(o(ii).name,'..')
        continue
    end
    try
        if o(ii).isdir
            rmdir([folder, '/', o(ii).name],'s')
        else
            delete([folder, '/', o(ii).name])
        end
    catch
    end
end
warning('on','MATLAB:RMDIR:RemovedFromPath')

author: BUNPC


%%  min cost of each iteration

iter = [1;2;3;4;5];
min_cost_2000 = [min_cost_0_2000; min_cost_1_2000; min_cost_2_2000; min_cost_3_2000; min_cost_4_2000];
min_cost_3000 = [min_cost_0_3000; min_cost_1_3000; min_cost_2_3000; min_cost_3_3000; min_cost_4_3000];
min_cost_1000 = [min_cost_0_1000; min_cost_1_1000; min_cost_2_1000; min_cost_3_1000; min_cost_4_1000];
min_cost_4000 = [min_cost_0_4000; min_cost_1_4000; min_cost_2_4000; min_cost_3_4000; min_cost_4_4000];

ff = figure;
plot(iter, min_cost_4000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'm'); hold on;
plot(iter, min_cost_3000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'b'); hold on;
plot(iter, min_cost_2000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'r'); hold on;
plot(iter, min_cost_1000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'k'); hold on;

grid on;
l=legend('b=4000', 'b=3000',  'b=2000', 'b=1000');
title(l, 'Cost after Progm');
xlabel('Iteration', 'FontWeight', 'bold', 'FontSize', 14);
ylabel('Cost', 'FontWeight', 'bold', 'FontSize', 14);
set(ff, 'Renderer', 'painters');
set(gca,'Xtick', 1:1:5);
set(gca,'FontSize',11);
set(gca,'Xticklabel',{'Iter 1','Iter 2', 'Iter 3', 'Iter 4', 'Iter 5'})

%% totalNum and Kept ratio

iter = [1;2;3;4;5];
totalNum_4000 = [totalNum_0_4000; totalNum_1_4000; totalNum_2_4000; totalNum_3_4000; totalNum_4_4000 ];
totalNum_3000 = [totalNum_0_3000; totalNum_1_3000; totalNum_2_3000; totalNum_3_3000; totalNum_4_3000];
totalNum_2000 = [totalNum_0_2000; totalNum_1_2000; totalNum_2_2000; totalNum_3_2000; totalNum_4_2000 ];
totalNum_1000 = [totalNum_0_1000; totalNum_1_1000; totalNum_2_1000; totalNum_3_1000; totalNum_4_1000 ];

% ff = figure;
% plot(iter, totalNum_2000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'r'); hold on;
% plot(iter, totalNum_1000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'm'); hold on;
% plot(iter, totalNum_4000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'k'); hold on;
% plot(iter, totalNum_4000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'm'); hold on;

% grid on;
% l=legend('b=2000', 'b=3000', 'b=1000', 'b=500');
% title(l, 'Number of Points after Progm');
% xlabel('Iteration', 'FontWeight', 'bold', 'FontSize', 14);
% ylabel('Number of Points', 'FontWeight', 'bold', 'FontSize', 14);
% set(ff, 'Renderer', 'painters');
% set(gca,'Xtick', 1:1:8, 'fontsize',14);
% set(gca,'Xticklabel',{'Iter 1','Iter 2', 'Iter 3', 'Iter 4', 'Iter 5', 'Iter 6', 'Iter 7', 'Iter 8'});
% set(gca,'FontSize',12);

iter = [1;2;3;4;5];
ratio_2000 = totalNum_2000 / totalNum;
ratio_3000 = totalNum_3000 / totalNum;
ratio_1000 = totalNum_1000 / totalNum;
ratio_4000 = totalNum_4000 / totalNum;

ff = figure;
plot(iter, ratio_4000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'm'); hold on;
plot(iter, ratio_3000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'b'); hold on;
plot(iter, ratio_2000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'r'); hold on;
plot(iter, ratio_1000, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'k'); hold on;

grid on;
l=legend('b=4000', 'b=3000', 'b=2000', 'b=1000');
title(l, 'Ratio after Progm');
xlabel('Iteration', 'FontWeight', 'bold', 'FontSize', 14);
ylabel('Compression Ratio', 'FontWeight', 'bold', 'FontSize', 14);
set(ff, 'Renderer', 'painters');
set(gca,'Xtick', 1:1:5);
set(gca,'Xticklabel',{'Iter 1','Iter 2', 'Iter 3', 'Iter 4', 'Iter 5'});
set(gca,'FontSize',11);

%%  Total time cost
%
% time_sum_2000 = sum([time_sum_0, time_sum_1, time_sum_2, time_sum_3, time_sum_4, time_sum_5, time_sum_6, time_sum_7]);
% time_sum_3000 = sum([time_sum_0, time_sum_1, time_sum_2, time_sum_3, time_sum_4, time_sum_5, time_sum_6, time_sum_7]);
% time_sum_1200 = sum([time_sum_0, time_sum_1, time_sum_2, time_sum_3, time_sum_4, time_sum_5, time_sum_6, time_sum_7]);
% time_sum_1000 = sum([time_sum_0, time_sum_1, time_sum_2, time_sum_3, time_sum_4, time_sum_5, time_sum_6, time_sum_7]);
% time_sum_800 = sum([time_sum_0, time_sum_1, time_sum_2, time_sum_3, time_sum_4, time_sum_5, time_sum_6, time_sum_7]);
%
%
%
% iter = [1;2;3;4;5];
% times = [time_sum_2000; time_sum_3000; time_sum_1200; time_sum_1000; time_sum_800];
% ff = figure;
% plot(iter, times, '--bs', 'MarkerSize', 10,  'MarkerFaceColor', 'r'); hold on;
% grid on;
% l=legend('time costs');
% title(l, 'Total time under Diff b');
% xlabel('Iteration', 'FontWeight', 'bold', 'FontSize', 12);
% ylabel('Time Costs (second)', 'FontWeight', 'bold', 'FontSize', 12);
% set(ff, 'Renderer', 'painters');
% set(gca,'Xtick', 1:1:8, 'fontsize',10);
% set(gca,'Xticklabel',{'b=2000','b=3000', 'b=1200', 'b=1000', 'b=800'});


text:
function [ y, symm ] = cvx_s_symmetric( m, n, symm ) %#ok
%CVX_S_SYMMETRIC Symmetric matrices (lower triangle storage).
if m ~= n,
    error( 'Symmetric structure requires square matrices.' );
end
symm = false;
nsq = n * n;
ntr = 0.5 * ( nsq + n );
c  = 0 : n - 1;
c  = c( ones( 1, n ), : );
r  = c';
mn = min( r, c );
mx = max( r, c );
y  = mx + mn .* ( n - 0.5 * ( mn + 1 ) ) + 1;
y  = sparse( y( : ), 1 : nsq, 1, ntr, nsq );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
%DIPIM Fixed mapping running any DIPimage command with one input image
%
%	B = DIPIM(A,COMMAND,PAR1,PAR2, ...)
%	B = A*DIPIM([],COMMAND,PAR1,PAR2, ...)
%	B = A*DIPIM(COMMAND,PAR1,PAR2, ...)
%
% INPUT
%   A        Dataset or datafiles with object images
%   COMMAND  String with name of a DIPimage command
%   PAR1     Additional parameters needed for COMMAND
%
% OUTPUT
%   B        Datasetor datafile with results
%
% DESCRIPTION
% This is a general routine to run DIPimage commands on the images stored
% in a dataset or datafile. It is thereby appropriate for grey value
% operations. For logical and morphological operations DIPBIN should be
% used. See DIPIMAGE for a list of all possible commands. A few suggections:
%     stretch              - Grey-value stretching
%     hist_equalize        - Histogram equalization
%     lut                  - Look-up Table (with interpolation)
%     convolve             - General convolution filter
%     gaussf               - Gaussian filter
%     unif                 - Uniform filter
%     maxf                 - Maximum filter
%     minf                 - Minimum filter
%     medif                - Median filter
%     percf                - Percentile filter
%     varif                - Variance filter
%     gabor                - Gabor filter
%     derivative           - Derivative filters
%     dx                   - First derivative in the X-direction
%     dy                   - First derivative in the Y-direction
%     dz                   - First derivative in the Z-direction
%     laplace              - Laplace operator
%     laplace_plus_dgg     - Laplace + Dgg
%     laplace_min_dgg      - Laplace - Dgg
%     hessian              - Hessian matrix of an image
%     prewittf             - Prewitt derivative filter
%     sobelf               - Sobel derivative filter
%     threshold            - Thresholding
%     hist2image           - Backmaps a 2D histogram ROI to the images
%     minima               - Detect local minima
%     maxima               - Detect local maxima
%     watershed            - Watershed
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, DIP_IMAGE, DIPBIN

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = dipim(varargin)

  checktoolbox('diplib');
	argin = shiftargin(varargin,'char');
  argin = setdefaults(argin,[],'');
  varargin = cell(1,numel(argin)-2);
  [a,command,varargin{:}] = deal(argin{:});
	if isempty(command)
    error('No DIPimage command supplied')
  end
  if mapping_task(argin,'definition')
    b = define_mapping(argin,'fixed');
    b = setname(b,'DIPimage');
  elseif isa(a,'prdataset') % allows datafiles too
		isobjim(a);
    b = filtim(a,mfilename,{command varargin{:}});
  elseif isa(a,'double') % here we have a single image
		a = dip_image(a);
    b = runcommand(command,a,varargin{:});
  elseif isa(a,'dip_image') % here we have a single DIP image
    b = runcommand(command,a,varargin{:});
  else
    error('Illegal call')
  end

 return

  function b = runcommand(command,a,varargin)

  try
    b = feval(command,a,varargin{:});
  catch ME
    if strcmp(ME.message, 'Argument # 1: image data type not supported')
      error('Wrong routine for this DIPimage command. Try DIPBIN.')
    end
    rethrow(ME);
  end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [ pn, dist ] = circle_arc_point_near_2d ( r, center, theta1, ...
  theta2, p )

%*****************************************************************************80
%
%% CIRCLE_ARC_POINT_NEAR_2D : nearest point on a circular arc.
%
%  Discussion:
%
%    A circular arc is defined by the portion of a circle (R,C)
%    between two angles (THETA1,THETA2).
%
%    Thus, a point (X,Y) on a circular arc satisfies
%
%      ( X - C(1) ) * ( X - C(1) ) + ( Y - C(2) ) * ( Y - C(2) ) = R * R
%
%    and
%
%      Theta1 <= Theta <= Theta2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2), the center of the circle.
%
%    Input, real THETA1, THETA2, the angles defining the arc,
%    in radians.  Normally, THETA1 < THETA2.
%
%    Input, real P(2), the point to be checked.
%
%    Output, real PN(2), a point on the circular arc which is
%    nearest to the point.
%
%    Output, real DIST, the distance to the nearest point.
%
  ndim = 2;
%
%  Special case, the zero circle.
%
  if ( r == 0.0 )
    pn(1:ndim) = center(1:ndim);
    dist = sqrt ( sum ( ( p(1:ndim) - pn(1:ndim) ).^2 ) );
    return
  end
%
%  Determine the angle made by the point.
%
  theta = atan4 ( p(2) - center(2), p(1) - center(1) );
%
%  If the angle is between THETA1 and THETA2, then you can
%  simply project the point onto the arc.
%
  if ( r8_modp ( theta  - theta1,  2.0 * pi ) <= ...
       r8_modp ( theta2 - theta1,  2.0 * pi ) )

    r2 = sqrt ( sum ( ( p(1:ndim) - center(1:ndim) ).^2 ) );

    pn(1:ndim) = center(1:ndim) + ( p(1:ndim) - center(1:ndim) ) * r / r2;
%
%  Otherwise, if the angle is less than the negative of the
%  average of THETA1 and THETA2, it's on the side of the arc
%  where the endpoint associated with THETA2 is closest.
%
%  As usual, MATLAB has some ninny problem with using cosine in a vector.
%  So I have to write the assignment as two scalar assignments.
%
  elseif ( r8_modp ( theta - 0.5 * ( theta1 + theta2 ), 2.0 * pi ) <= pi )

    pn(1) = center(1) + r * cos ( theta2 );
    pn(2) = center(2) + r * sin ( theta2 );
%
%  Otherwise, the endpoint associated with THETA1 is closest.
%
  else

    pn(1) = center(1) + r * cos ( theta1 );
    pn(2) = center(2) + r * sin ( theta1 );

  end

  dist = sqrt ( sum ( ( p(1:ndim) - pn(1:ndim) ).^2 ) );

  return
end


text:
%DNA damage process test case
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 7/13/2010
classdef DNADamage_Test < edu.stanford.covert.cell.sim.ReactionProcessTestCase
    %constants
    properties (Constant=true)
        expected_essentialGenes = {};
    end

    %constructor
    methods
        function this = DNADamage_Test(name)
            this = this@edu.stanford.covert.cell.sim.ReactionProcessTestCase(name);
        end
    end

    %simple test fixture
    methods
        function loadSimpleTestFixture(this)
            import edu.stanford.covert.cell.sim.constant.ChromosomeSequence;
            import edu.stanford.covert.cell.sim.state.Chromosome;

            %process
            m = this.process;

            %% chromosome state
            c = Chromosome([], []);
            c.monomerDNAFootprints = zeros(0, 1);
            c.complexDNAFootprints = zeros(0, 1);
            c.monomerDNAFootprintBindingStrandedness = repmat(c.dnaStrandedness_dsDNA, 0, 1);
            c.complexDNAFootprintBindingStrandedness = repmat(c.dnaStrandedness_dsDNA, 0, 1);
            c.monomerDNAFootprintRegionStrandedness = repmat(c.dnaStrandedness_dsDNA, 0, 1);
            c.complexDNAFootprintRegionStrandedness = repmat(c.dnaStrandedness_dsDNA, 0, 1);
            c.reactionBoundMonomer = zeros(0, 1);
            c.reactionBoundComplex = zeros(0, 1);
            c.reactionMonomerCatalysisMatrix = zeros(0, 0);
            c.reactionComplexCatalysisMatrix = zeros(0, 0);
            c.reactionThresholds = ...
                sum(c.reactionMonomerCatalysisMatrix, 2) + ...
                sum(c.reactionComplexCatalysisMatrix, 2);
            m.states = {c};
            m.chromosome = c;

            %parameters
            c.relaxedBasesPerTurn = 10.5;
            c.equilibriumSuperhelicalDensity = -0.06;
            c.supercoiledSuperhelicalDensityTolerance = 0.10;

            %genome
            seq = reshape([
                'AGCATTTTCAGGCACGGGGATACGGTTCATGGCGCTCCCGGAGACGACAA'
                'TCCAACCCGCTTGCGCCTTATAATCATTTCTCGAGCGAAAAACAAAAGCC'
                'ATAGTGTCTCTGCGCATCACTCGCACCAATGAATTCCGCGTCTGATAGCC'
                'TGCTTTGAGTCCCGGCTAATGTACAGTCCGATCACCCAACACCGGCAAGC'
                'GTTGCGGCACTAAGCTGGCACGATATACCGGTCGGGCGCCTGGCAGCACC'
                'AACTGTCGATAGGATCAGCGCTGTATCTAGAATGTGAGCTTGGCGGTCAG'
                'CTATTCTTGAATGAATATTCTTAGGTCGAAGCCTGCATTCAGCACCGCGG'
                'CGCTGCTAATTCTTTCGTTGGTTGGGCCAGCGAATGTGCGCCCTTCCGTA'
                'CTGGTCATGGTGCGGAGAAGCAACGTAATGACCGGAACATCCATTACGAC'
                'CCTAATCGAAGCTGACAGTTACTAGCGCTAGACGAGACGTACCAGGGAAG'],...
                1,[]);
            c.sequence = ChromosomeSequence(seq);
            c.sequenceLen = length(seq);
            c.sequenceGCContent = getGCContent(c.sequence);

            %% enzymes
            m.enzymeWholeCellModelIDs = cell(0,1);
            m.enzymeNames = cell(0,1);
            m.enzymeMolecularWeights = zeros(0,1);

            m.enzymes       = zeros(0, 1);
            m.boundEnzymes  = zeros(0, 1);

            %% substrates
            %whole cell model IDs, molecular weights
            m.substrateWholeCellModelIDs = {
                'AD'; 'CSN'; 'DA'; 'DAMP'; 'DC'; 'DCMP'; 'DG'; 'DGMP'; 'DHTHY'; 'DR5P';
                'DT'; 'DTMP'; 'ENU'; 'GLC'; 'GN'; 'H'; 'H2O'; 'MNU'; 'NH3'; 'NU'; 'O4eTHY';
                'O4mTHY'; 'PTRC'; 'THY'; 'THY64THY'; 'URA'; 'UVB_radiation'; 'butylene_glycol';
                'cyclodG'; 'dRibose5P_dRibose5P'; 'dhpURA'; 'gamma_radiation'; 'glchm5URA';
                'hm5URA'; 'ho5URA'; 'hydrogen_radical'; 'hydroxyl_radical'; 'm5CSN'; 'putTHY'};
            m.substrateMolecularWeights = [
                135.1265; 111.1018; 251.2414; 329.2055; 227.2167; 305.1808; 267.2408;
                345.2049; 128.1289; 212.0942; 242.2280; 320.1921; 117.1063; 180.1554;
                151.1259;   1.0079;  18.0152; 103.0798;  17.0304;  89.0533; 154.1661;
                140.1396;  90.1668; 126.1131; 252.2262; 112.0866;        0;  90.1206;
                265.2250; 407.1811; 214.2179;        0; 304.2527; 142.1125; 128.0860;
                1.0079;  17.0073; 125.1283; 212.2483];

            %names
            m.substrateNames = m.substrateWholeCellModelIDs;

            %indices
            m.substrateIndexs_radiation             = m.substrateIndexs({'UVB_radiation';'gamma_radiation'});
            m.substrateIndexs_hydrogen              = m.substrateIndexs({'H'});
            m.substrateIndexs_hydroxylRadical       = m.substrateIndexs({'hydroxyl_radical'});
            m.substrateIndexs_water                 = m.substrateIndexs({'H2O'});
            m.substrateIndexs_DR5P                  = m.substrateIndexs({'DR5P'});
            m.substrateIndexs_dRibose5P_dRibose5P   = m.substrateIndexs({'dRibose5P_dRibose5P'});
            m.substrateIndexs_unmodifiedNucleobases = m.substrateIndexs({'AD';'CSN';'GN';'THY'});
            m.substrateIndexs_unmodifiedNucleosides = m.substrateIndexs({'DA';'DC';'DG';'DT'});
            m.substrateIndexs_unmodifiedDNMPs       = m.substrateIndexs({'DAMP';'DCMP';'DGMP';'DTMP'});
            m.substrateIndexs_modifiedNucleobases   = m.substrateIndexs({
                'DHTHY'; 'O4eTHY'; 'O4mTHY'; 'THY64THY'; 'URA'; 'dhpURA';
                'glchm5URA'; 'hm5URA'; 'ho5URA'; 'm5CSN'; 'putTHY'});
            m.substrateIndexs_modifiedNucleosides   = m.substrateIndexs({'cyclodG'});
            m.substrateIndexs_modifiedDNMPs         = m.substrateIndexs({});

            m.substrateMetaboliteLocalIndexs = (1:numel(m.substrateWholeCellModelIDs))';
            m.substrateMetaboliteGlobalIndexs = m.substrateMetaboliteLocalIndexs;

            c.metabolite.dr5pIndexs = m.substrateIndexs({'DR5P'}) ;
            c.metabolite.waterIndexs = m.substrateIndexs({'H2O'});
            c.metabolite.dnmpIndexs = m.substrateIndexs({'AD'; 'CSN'; 'GN'; 'THY'});
            c.metabolite.m6ADIndexs = numel(m.substrateWholeCellModelIDs)+1;

            c.metabolite.molecularWeights = [m.substrateMolecularWeights; 0];

            %% reactions -- 1 reaction of each type
            %whole cell model IDs
            m.reactionWholeCellModelIDs = {
                'DNADamage_THYTHY_THY64THY_UVB_radiation';
                'DNADamage_hm5U_dhpURA_butylene_glycol';
                'DNADamage_hm5U_putTHY_PTRC';
                'DNADamage_THY_O4eTHY_ENU';
                'DNADamage_hm5U_glchm5URA_GLC';
                'DNADamage_THY_O4mTHY_MNU';
                'DNADamage_THY_DHTHY_hydrogen_radical';
                'DNADamage_DG_cyclodG_UVB_radiation';
                'DNADamage_URA_ho5URA_hydroxyl_radical';
                'DNADamage_SpontaneousBaseDeamination_m5C';
                'DNADamage_SpontaneousBaseLoss_thymine';
                'DNADamage_SpontaneousAbasicSiteStrandBreak'};

            %names
            m.reactionNames = m.reactionWholeCellModelIDs;

            %types
            m.reactionTypes = {
                'UV-B photodimerization';
                'base alkylation';
                'base amination';
                'base ethylation';
                'base glucosyl transfer';
                'base methylation';
                'base reduction';
                'photooxidation';
                'radiation (gamma-ray) induced base oxidation';
                'spontaneous base deamination';
                'spontaneous base loss';
                'strand break'};

            %indices
            m.reactionIndexs_spontaneousBaseLoss        = find(strcmp(m.reactionTypes,'spontaneous base loss'));
            m.reactionIndexs_spontaneousBaseDeamination = find(strcmp(m.reactionTypes,'spontaneous base deamination'));
            m.reactionIndexs_gammaRayBaseOxidation      = find(strcmp(m.reactionTypes,'radiation (gamma-ray) induced base oxidation'));
            m.reactionIndexs_uvbDimerization            = find(strcmp(m.reactionTypes,'UV-B photodimerization'));
            m.reactionIndexs_strandBreak                = find(strcmp(m.reactionTypes,'strand break'));
            m.reactionIndexs_baseReduction              = find(strcmp(m.reactionTypes,'base reduction'));
            m.reactionIndexs_baseAmination              = find(strcmp(m.reactionTypes,'base amination'));
            m.reactionIndexs_baseGlucosylTransfer       = find(strcmp(m.reactionTypes,'base glucosyl transfer'));
            m.reactionIndexs_baseAlkylation             = find(strcmp(m.reactionTypes,'base alkylation'));
            m.reactionIndexs_baseEthylation             = find(strcmp(m.reactionTypes,'base ethylation'));
            m.reactionIndexs_baseMethylation            = find(strcmp(m.reactionTypes,'base methylation'));
            m.reactionIndexs_photooxidation             = find(strcmp(m.reactionTypes,'photooxidation'));

            %ReactionProcess inherited reaction properties
            m.reactionCatalysisMatrix = zeros(length(m.reactionWholeCellModelIDs), 0);
            m.enzymeBounds = repmat([-Inf Inf], length(m.reactionWholeCellModelIDs), 1);

            %% reaction stoichiometry matrices, reaction sequences
            %allocate memory
            m.reactionStoichiometryMatrix              = zeros(length(m.substrateWholeCellModelIDs), length(m.reactionWholeCellModelIDs));
            m.reactionRadiationStoichiometryMatrix     = zeros(length(m.substrateWholeCellModelIDs), length(m.reactionWholeCellModelIDs));
            m.reactionSmallMoleculeStoichiometryMatrix = zeros(length(m.substrateWholeCellModelIDs), length(m.reactionWholeCellModelIDs));
            m.reactionDNAStoichiometryMatrix           = zeros(length(m.substrateWholeCellModelIDs), length(m.reactionWholeCellModelIDs));
            m.reactionVulnerableMotifs                 = cell(length(m.reactionWholeCellModelIDs), 1);
            m.reactionVulnerableMotifTypes             = cell(length(m.reactionWholeCellModelIDs), 1);
            m.reactionDamageTypes                      = cell(length(m.reactionWholeCellModelIDs), 1);
            m.reactionBounds                           = zeros(length(m.reactionWholeCellModelIDs), 2);

            %DNADamage_THYTHY_THY64THY_UVB_radiation
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'THY'}), m.reactionIndexs_uvbDimerization) = -2;
            m.reactionRadiationStoichiometryMatrix(m.substrateIndexs({'UVB_radiation'}), m.reactionIndexs_uvbDimerization) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'THY64THY'}), m.reactionIndexs_uvbDimerization) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_uvbDimerization} = 'TT';
            m.reactionVulnerableMotifTypes{m.reactionIndexs_uvbDimerization} = 'sequence';
            m.reactionDamageTypes{m.reactionIndexs_uvbDimerization} = 'intrastrandCrossLinks';
            m.reactionBounds(m.reactionIndexs_uvbDimerization, 2) = 1.8300e-008;

            %DNADamage_hm5U_dhpURA_butylene_glycol
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'butylene_glycol'}), m.reactionIndexs_baseAlkylation) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'hm5URA'}), m.reactionIndexs_baseAlkylation) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'dhpURA'}), m.reactionIndexs_baseAlkylation) = 1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H2O'}), m.reactionIndexs_baseAlkylation) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_baseAlkylation}  = m.substrateIndexs({'hm5URA'});
            m.reactionVulnerableMotifTypes{m.reactionIndexs_baseAlkylation} = 'damagedBases';
            m.reactionDamageTypes{m.reactionIndexs_baseAlkylation} = 'damagedBases';
            m.reactionBounds(m.reactionIndexs_baseAlkylation, 2) = 0;

            %DNADamage_hm5U_putTHY_PTRC
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'hm5URA'}), m.reactionIndexs_baseAmination) = -1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'PTRC'}), m.reactionIndexs_baseAmination) = -1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H'}), m.reactionIndexs_baseAmination) = 2;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H2O'}), m.reactionIndexs_baseAmination) = 1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'putTHY'}), m.reactionIndexs_baseAmination) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_baseAmination}  = m.substrateIndexs({'hm5URA'});
            m.reactionVulnerableMotifTypes{m.reactionIndexs_baseAmination} = 'damagedBases';
            m.reactionDamageTypes{m.reactionIndexs_baseAmination} = 'damagedBases';
            m.reactionBounds(m.reactionIndexs_baseAmination, 2) = 0;

            %DNADamage_THY_O4eTHY_ENU
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'ENU'}), m.reactionIndexs_baseEthylation) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'THY'}), m.reactionIndexs_baseEthylation) = -1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'NU'}), m.reactionIndexs_baseEthylation) = 1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'O4eTHY'}), m.reactionIndexs_baseEthylation) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_baseEthylation}  = 'T';
            m.reactionVulnerableMotifTypes{m.reactionIndexs_baseEthylation} = 'sequence';
            m.reactionDamageTypes{m.reactionIndexs_baseEthylation} = 'damagedBases';
            m.reactionBounds(m.reactionIndexs_baseEthylation, 2) = 0;

            %DNADamage_hm5U_glchm5URA_GLC
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'GLC'}), m.reactionIndexs_baseGlucosylTransfer) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'hm5URA'}), m.reactionIndexs_baseGlucosylTransfer) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'glchm5URA'}), m.reactionIndexs_baseGlucosylTransfer) = 1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H2O'}), m.reactionIndexs_baseGlucosylTransfer) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_baseGlucosylTransfer}  = m.substrateIndexs({'hm5URA'});
            m.reactionVulnerableMotifTypes{m.reactionIndexs_baseGlucosylTransfer} = 'damagedBases';
            m.reactionDamageTypes{m.reactionIndexs_baseGlucosylTransfer} = 'damagedBases';
            m.reactionBounds(m.reactionIndexs_baseGlucosylTransfer, 2) = 0;

            %DNADamage_THY_O4mTHY_MNU
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'MNU'}), m.reactionIndexs_baseMethylation) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'THY'}), m.reactionIndexs_baseMethylation) = -1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'NU'}), m.reactionIndexs_baseMethylation) = 1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'O4mTHY'}), m.reactionIndexs_baseMethylation) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_baseMethylation}  = 'T';
            m.reactionVulnerableMotifTypes{m.reactionIndexs_baseMethylation} = 'sequence';
            m.reactionDamageTypes{m.reactionIndexs_baseMethylation} = 'damagedBases';
            m.reactionBounds(m.reactionIndexs_baseMethylation, 2) = 0;

            %DNADamage_THY_DHTHY_hydrogen_radical
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'hydrogen_radical'}), m.reactionIndexs_baseReduction) = -2;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'THY'}), m.reactionIndexs_baseReduction) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'DHTHY'}), m.reactionIndexs_baseReduction) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_baseReduction}  = 'T';
            m.reactionVulnerableMotifTypes{m.reactionIndexs_baseReduction} = 'sequence';
            m.reactionDamageTypes{m.reactionIndexs_baseReduction} = 'damagedBases';
            m.reactionBounds(m.reactionIndexs_baseReduction, 2) = 0;

            %DNADamage_DG_cyclodG_UVB_radiation
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'DG'}), m.reactionIndexs_photooxidation) = -1;
            m.reactionRadiationStoichiometryMatrix(m.substrateIndexs({'UVB_radiation'}), m.reactionIndexs_photooxidation) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'cyclodG'}), m.reactionIndexs_photooxidation) = 1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H'}), m.reactionIndexs_photooxidation) = 2;
            m.reactionVulnerableMotifs{m.reactionIndexs_photooxidation}  = 'G';
            m.reactionVulnerableMotifTypes{m.reactionIndexs_photooxidation} = 'sequence';
            m.reactionDamageTypes{m.reactionIndexs_photooxidation} = 'damagedBases';
            m.reactionBounds(m.reactionIndexs_photooxidation, 2) = 0;

            %DNADamage_URA_ho5URA_hydroxyl_radical
            m.reactionRadiationStoichiometryMatrix(m.substrateIndexs({'gamma_radiation'}), m.reactionIndexs_gammaRayBaseOxidation) = -1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H2O'}), m.reactionIndexs_gammaRayBaseOxidation) = -1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H'}), m.reactionIndexs_gammaRayBaseOxidation) = 2;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'URA'}), m.reactionIndexs_gammaRayBaseOxidation) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'ho5URA'}), m.reactionIndexs_gammaRayBaseOxidation) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_gammaRayBaseOxidation}  = m.substrateIndexs({'URA'});
            m.reactionVulnerableMotifTypes{m.reactionIndexs_gammaRayBaseOxidation} = 'damagedBases';
            m.reactionDamageTypes{m.reactionIndexs_gammaRayBaseOxidation} = 'damagedBases';
            m.reactionBounds(m.reactionIndexs_gammaRayBaseOxidation, 2) = 5.8999e-008;

            %DNADamage_SpontaneousBaseDeamination_m5C
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H2O'}), m.reactionIndexs_spontaneousBaseDeamination) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'m5CSN'}), m.reactionIndexs_spontaneousBaseDeamination) = -1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'NH3'}), m.reactionIndexs_spontaneousBaseDeamination) = 1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'THY'}), m.reactionIndexs_spontaneousBaseDeamination) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_spontaneousBaseDeamination}  = m.substrateIndexs({'m5CSN'});
            m.reactionVulnerableMotifTypes{m.reactionIndexs_spontaneousBaseDeamination} = 'damagedBases';
            m.reactionDamageTypes{m.reactionIndexs_spontaneousBaseDeamination} = 'damagedBases';
            m.reactionBounds(m.reactionIndexs_spontaneousBaseDeamination, 2) = 2.2000e-014;

            %DNADamage_SpontaneousBaseLoss_thymine
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'DTMP'}), m.reactionIndexs_spontaneousBaseLoss) = -1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H2O'}), m.reactionIndexs_spontaneousBaseLoss) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'DR5P'}), m.reactionIndexs_spontaneousBaseLoss) = 1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'THY'}), m.reactionIndexs_spontaneousBaseLoss) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_spontaneousBaseLoss}  = 'T';
            m.reactionVulnerableMotifTypes{m.reactionIndexs_spontaneousBaseLoss} = 'sequence';
            m.reactionDamageTypes{m.reactionIndexs_spontaneousBaseLoss} = 'abasicSites';
            m.reactionBounds(m.reactionIndexs_spontaneousBaseLoss, 2) = 6.0000e-012;

            %DNADamage_SpontaneousAbasicSiteStrandBreak
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'dRibose5P_dRibose5P'}), m.reactionIndexs_strandBreak) = -1;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H2O'}), m.reactionIndexs_strandBreak) = -1;
            m.reactionDNAStoichiometryMatrix(m.substrateIndexs({'DR5P'}), m.reactionIndexs_strandBreak) = 2;
            m.reactionSmallMoleculeStoichiometryMatrix(m.substrateIndexs({'H'}), m.reactionIndexs_strandBreak) = 1;
            m.reactionVulnerableMotifs{m.reactionIndexs_strandBreak}  = true;
            m.reactionVulnerableMotifTypes{m.reactionIndexs_strandBreak}  = 'abasicSites';
            m.reactionDamageTypes{m.reactionIndexs_strandBreak} = 'strandBreaks';
            m.reactionBounds(m.reactionIndexs_strandBreak, 2) = 1.9254e-006;

            %reaction DNA product
            [i, j] = find(m.reactionDNAStoichiometryMatrix > 0);
            m.reactionDNAProduct = zeros(numel(m.reactionWholeCellModelIDs), 1);
            m.reactionDNAProduct(j) = m.substrateMetaboliteGlobalIndexs(ismembc2(i, m.substrateMetaboliteLocalIndexs));

            %reaction radiation
            [i, j] = find(m.reactionRadiationStoichiometryMatrix);
            m.reactionRadiation = zeros(numel(m.reactionWholeCellModelIDs), 1);
            m.reactionRadiation(j) = i;

            %reactionStoichiometryMatrix
            m.reactionStoichiometryMatrix = ...
                m.reactionRadiationStoichiometryMatrix + ...
                m.reactionSmallMoleculeStoichiometryMatrix + ...
                m.reactionDNAStoichiometryMatrix;

            %% initial state
            m.substrates = zeros(length(m.substrateWholeCellModelIDs), 1);

            m.chromosome.initialize();
        end
    end

    %tests
    methods
        function testConstants(this)
            m = this.process;
            c = m.chromosome;

            %assert that are sorted
            assertTrue(issorted(m.substrateMetaboliteLocalIndexs));

            %assert that this process doesn't create or modify m6AD
            assertFalse(any(m.substrateMetaboliteGlobalIndexs == c.metabolite.m6ADIndexs));
            for i = 1:numel(m.reactionVulnerableMotifs)
                assertFalse(isequal(m.reactionVulnerableMotifs{i}, c.metabolite.m6ADIndexs));
            end

            %check that reaction properties were loaded correctly from KB
            substrateWholeCellModelIDs               = m.substrateWholeCellModelIDs;
            substrateMetaboliteGlobalIndexs          = m.substrateMetaboliteGlobalIndexs;
            substrateMetaboliteLocalIndexs           = m.substrateMetaboliteLocalIndexs;
            reactionWholeCellModelIDs                = m.reactionWholeCellModelIDs;
            reactionTypes                            = m.reactionTypes;
            reactionDamageTypes                      = m.reactionDamageTypes;
            reactionVulnerableMotifs                 = m.reactionVulnerableMotifs;
            reactionVulnerableMotifTypes             = m.reactionVulnerableMotifTypes;
            reactionSmallMoleculeStoichiometryMatrix = m.reactionSmallMoleculeStoichiometryMatrix;
            reactionDNAStoichiometryMatrix           = m.reactionDNAStoichiometryMatrix;
            reactionRadiationStoichiometryMatrix     = m.reactionRadiationStoichiometryMatrix;
            reactionBounds                           = m.reactionBounds;

            this.loadSimpleTestFixture();

            [~, sbIdxs1, sbIdxs2] = intersect(substrateWholeCellModelIDs, m.substrateWholeCellModelIDs);
            [~, rxIdxs1, rxIdxs2] = intersect(reactionWholeCellModelIDs, m.reactionWholeCellModelIDs);

            for i=1:length(rxIdxs1)
                motif = reactionVulnerableMotifs{rxIdxs1(i)};
                mMotif = m.reactionVulnerableMotifs{rxIdxs2(i)};
                if isnumeric(motif) || isnumeric(mMotif)
                    motif = substrateWholeCellModelIDs{substrateMetaboliteLocalIndexs(motif == substrateMetaboliteGlobalIndexs)};
                    mMotif = m.substrateWholeCellModelIDs{m.substrateMetaboliteLocalIndexs(mMotif == m.substrateMetaboliteGlobalIndexs)};
                end

                assertEqual(reactionTypes{rxIdxs1(i)},                                    m.reactionTypes{rxIdxs2(i)},                                    sprintf('Error in definition of reaction %s', m.reactionWholeCellModelIDs{rxIdxs2(i)}));
                assertEqual(reactionDamageTypes{rxIdxs1(i)},                              m.reactionDamageTypes{rxIdxs2(i)},                              sprintf('Error in definition of reaction %s', m.reactionWholeCellModelIDs{rxIdxs2(i)}));
                assertEqual(motif,                                                        mMotif,                                                         sprintf('Error in definition of reaction %s', m.reactionWholeCellModelIDs{rxIdxs2(i)}));
                assertEqual(reactionVulnerableMotifTypes{rxIdxs1(i)},                     m.reactionVulnerableMotifTypes{rxIdxs2(i)},                     sprintf('Error in definition of reaction %s', m.reactionWholeCellModelIDs{rxIdxs2(i)}));
                assertEqual(reactionSmallMoleculeStoichiometryMatrix(sbIdxs1,rxIdxs1(i)), m.reactionSmallMoleculeStoichiometryMatrix(sbIdxs2,rxIdxs2(i)), sprintf('Error in definition of reaction %s', m.reactionWholeCellModelIDs{rxIdxs2(i)}));
                assertEqual(reactionDNAStoichiometryMatrix(sbIdxs1,rxIdxs1(i)),           m.reactionDNAStoichiometryMatrix(sbIdxs2,rxIdxs2(i)),           sprintf('Error in definition of reaction %s', m.reactionWholeCellModelIDs{rxIdxs2(i)}));
                assertEqual(reactionRadiationStoichiometryMatrix(sbIdxs1,rxIdxs1(i)),     m.reactionRadiationStoichiometryMatrix(sbIdxs2,rxIdxs2(i)),     sprintf('Error in definition of reaction %s', m.reactionWholeCellModelIDs{rxIdxs2(i)}));
                assertEqual(reactionBounds(rxIdxs1(i),1),                                 m.reactionBounds(rxIdxs2(i),1),                                 sprintf('Error in definition of reaction %s', m.reactionWholeCellModelIDs{rxIdxs2(i)}));
                assertElementsAlmostEqual(reactionBounds(rxIdxs1(i),2), m.reactionBounds(rxIdxs2(i),2), 'relative','1e-4', sprintf('Error in definition of reaction %s',m.reactionWholeCellModelIDs{rxIdxs2(i)}));
            end
        end

        function testReactionStoichiometryMatrixFactorization(this)
            m = this.process;

            dMW1 = m.reactionStoichiometryMatrix'*m.substrateMolecularWeights;
            dMW2 = (...
                m.reactionSmallMoleculeStoichiometryMatrix+...
                m.reactionDNAStoichiometryMatrix+...
                m.reactionRadiationStoichiometryMatrix)'*m.substrateMolecularWeights;

            assertElementsAlmostEqual(dMW1, dMW2);
        end

        function testBaseAlkylation(this)
            this.loadSimpleTestFixture();

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_hm5U_dhpURA_butylene_glycol'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 1;

            %set DNA to undamaged
            this.clearDNADamages();
            m.chromosome.damagedBases(1,1) = m.substrateMetaboliteGlobalIndexs(m.substrateMetaboliteLocalIndexs == find(m.reactionDNAStoichiometryMatrix(:,rxIdx)<0));

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertEqual([1 1], subs);

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testBaseAmination(this)
            this.loadSimpleTestFixture();

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_hm5U_putTHY_PTRC'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 1;

            %set DNA to undamaged
            this.clearDNADamages();
            m.chromosome.damagedBases(1,1) = m.substrateMetaboliteGlobalIndexs(m.substrateMetaboliteLocalIndexs == find(m.reactionDNAStoichiometryMatrix(:,rxIdx)<0));

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertEqual([1 1], subs);

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testBaseEthylation(this)
            this.loadSimpleTestFixture();

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_THY_O4eTHY_ENU'});

            %set amount of substrates
            m.substrates = ...
                max(0, -m.reactionSmallMoleculeStoichiometryMatrix(:, rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 10 / size(m.chromosome.sequence,1);

            %set DNA to undamaged
            this.clearDNADamages();

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertTrue(~isempty(subs));
            assertEqual(m.reactionVulnerableMotifs{rxIdx}, this.reactedSequences(subs, rxIdx));

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testBaseGlucosylTransfer(this)
            this.loadSimpleTestFixture();

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_hm5U_glchm5URA_GLC'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 1;

            %set DNA to undamaged
            this.clearDNADamages();
            m.chromosome.damagedBases(1,1) = m.substrateMetaboliteGlobalIndexs(m.substrateMetaboliteLocalIndexs == find(m.reactionDNAStoichiometryMatrix(:,rxIdx)<0));

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertEqual([1 1], subs);

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testBaseMethylation(this)
            this.loadSimpleTestFixture();

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_THY_O4mTHY_MNU'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 10 / size(m.chromosome.sequence,1);

            %set DNA to undamaged
            this.clearDNADamages();

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertEqual(m.reactionVulnerableMotifs{rxIdx}, this.reactedSequences(subs, rxIdx));

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testBaseReduction(this)
            this.loadSimpleTestFixture();

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_THY_DHTHY_hydrogen_radical'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 10 / size(m.chromosome.sequence,1);

            %set DNA to undamaged
            this.clearDNADamages();

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertEqual(m.reactionVulnerableMotifs{rxIdx}, this.reactedSequences(subs, rxIdx));

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testPhotooxidation(this)
            this.loadSimpleTestFixture();

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_DG_cyclodG_UVB_radiation'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionRadiationStoichiometryMatrix(:,rxIdx)) + ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 10 / size(m.chromosome.sequence,1);

            %set DNA to undamaged
            this.clearDNADamages();

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertEqual(m.reactionVulnerableMotifs{rxIdx}, this.reactedSequences(subs, rxIdx));

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                max(0,-m.reactionRadiationStoichiometryMatrix(:,rxIdx)) + ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testGammaRayBaseOxidation(this)
            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_URA_ho5URA_hydroxyl_radical'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionRadiationStoichiometryMatrix(:,rxIdx)) + ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 1;

            %set DNA to undamaged
            this.clearDNADamages();
            m.chromosome.damagedBases(1,1) = m.substrateMetaboliteGlobalIndexs(m.substrateMetaboliteLocalIndexs == find(m.reactionDNAStoichiometryMatrix(:,rxIdx)<0));

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damagedthis.reactedSequences(subs, rxIdx)
            assertEqual([1 1], subs);

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                max(0,-m.reactionRadiationStoichiometryMatrix(:,rxIdx)) + ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testSpontaneousBaseDeamination(this)
            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_SpontaneousBaseDeamination_m5C'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 1;

            %set DNA to undamaged
            this.clearDNADamages();
            m.chromosome.damagedBases(1,1) = m.substrateMetaboliteGlobalIndexs(m.substrateMetaboliteLocalIndexs == find(m.reactionDNAStoichiometryMatrix(:,rxIdx)<0));

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertEqual([1 1], subs);

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testSpontaneousBaseLoss(this)
            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_SpontaneousBaseLoss_thymine'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 10 / size(m.chromosome.sequence,1);

            %set DNA to undamaged
            this.clearDNADamages();

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.abasicSites);

            %correct sequence damaged
            assertEqual(m.reactionVulnerableMotifs{rxIdx}, this.reactedSequences(subs, rxIdx));

            %correct damage recorded
            assertEqual(true, unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.abasicSites, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testStrandBreak(this)
            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_SpontaneousAbasicSiteStrandBreak'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 1;

            %set DNA to undamaged
            this.clearDNADamages();
            m.chromosome.abasicSites(1, 1) = true;

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.strandBreaks);

            %correct sequence damaged
            assertEqual([1 1], subs);

            %correct damage recorded
            assertEqual(true, unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.strandBreaks | m.chromosome.strandBreaks5, m.chromosome.getDamagedSites(true, true, true, true, false, false, true) ~= 0);

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testUVBDimerization(this)
            import edu.stanford.covert.util.CircularSparseMat;

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_THYTHY_THY64THY_UVB_radiation'});

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionRadiationStoichiometryMatrix(:,rxIdx)) + ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 10 / size(m.chromosome.sequence,1);

            %set DNA to undamaged
            this.clearDNADamages();

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.intrastrandCrossLinks);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertEqual(m.reactionVulnerableMotifs{rxIdx}, this.reactedSequences(subs, rxIdx));

            %correct damage recorded
            assertEqual(find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.intrastrandCrossLinks + m.chromosome.intrastrandCrossLinks5, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = size(subs,1);
            assertEqual(m.substrates, ...
                max(0,-m.reactionRadiationStoichiometryMatrix(:,rxIdx)) + ...
                nReactions * max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)));
        end

        function testMultipleReactions(this)
            this.loadSimpleTestFixture();

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_hm5U_dhpURA_butylene_glycol'; 'DNADamage_THY_O4eTHY_ENU'});

            %set amount of substrates
            m.substrates = sum(...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)),...
                2);

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = [1; 10 / size(m.chromosome.sequence,1)];

            %set DNA to undamaged
            this.clearDNADamages();
            m.chromosome.damagedBases(1,1) = m.substrateMetaboliteGlobalIndexs(m.substrateMetaboliteLocalIndexs == find(m.reactionDNAStoichiometryMatrix(:,rxIdx(1))<0));

            %evolve state
            m.evolveState();

            %find locations of damages
            [subs, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);

            %correct sequence damaged
            assertEqual([1 1], subs(vals == find(m.reactionDNAStoichiometryMatrix(:,rxIdx(1))>0),:));
            assertEqual(m.reactionVulnerableMotifs{rxIdx(2)}, this.reactedSequences(subs(vals == find(m.reactionDNAStoichiometryMatrix(:,rxIdx(2))>0),:), rxIdx(2)));

            %correct damage recorded
            [~, vals] = find(m.chromosome.damagedBases);
            [~, idxs] = ismember(vals, m.substrateMetaboliteGlobalIndexs);
            vals = m.substrateMetaboliteLocalIndexs(idxs, :);
            assertEqual(find(any(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0,2)), unique(vals));

            %no other damages occured
            assertEqual(m.chromosome.damagedBases, m.chromosome.getDamagedSites(true, true, true, true, false, false, true));

            %substrates updated correctly
            nReactions = sum(m.reactionDNAStoichiometryMatrix(vals,rxIdx)>0,1)';
            assertEqual(m.substrates, ...
                max(0,m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx)) * nReactions);
        end

        function testFairness(this)
            this.loadSimpleTestFixture();

            m = this.process;
            rxIdx = m.reactionIndexs({'DNADamage_hm5U_dhpURA_butylene_glycol'});
            dnaSubstrate = find(m.reactionDNAStoichiometryMatrix(:,rxIdx)<0);
            dnaProduct = find(m.reactionDNAStoichiometryMatrix(:,rxIdx)>0);

            dnaSubstrate = m.substrateMetaboliteGlobalIndexs(m.substrateMetaboliteLocalIndexs == dnaSubstrate);
            dnaProduct = m.substrateMetaboliteGlobalIndexs(m.substrateMetaboliteLocalIndexs == dnaProduct);

            %set amount of substrates
            m.substrates = ...
                max(0,-m.reactionSmallMoleculeStoichiometryMatrix(:,rxIdx));
            initial_substrates = m.substrates;

            %increase rate of damage
            m.reactionBounds(:,2) = 0;
            m.reactionBounds(rxIdx, 2) = 1;

            %set DNA to undamaged
            this.clearDNADamages();
            m.chromosome.damagedBases(1:2,1) = dnaSubstrate;
            damagedBases = m.chromosome.damagedBases;

            counts = zeros(2,1);
            for i=1:200
                %reset state
                m.substrates = initial_substrates;
                m.chromosome.damagedBases = damagedBases;

                %evolve state
                m.evolveState();

                %find locations of damages
                [subs, vals] = find(m.chromosome.damagedBases);
                subs = subs(vals == dnaProduct,1);

                %update counts
                counts(subs) = counts(subs) + 1;
            end

            %check positions modified fairly
            assertTrue(max(counts)-min(counts) < 0.2 * max(counts));
        end

        function testGeneEssentiality(this)
            m = this.process;
            c = m.chromosome;
            c.polymerizedRegions(1, 1:2) = size(c.sequence, 1);

            m.substrates(:) = 1e3;
            m.substrates(m.substrateIndexs_radiation) = 1;

            %increase rates of damage
            m.reactionBounds(:,2) = 10 / size(m.chromosome.sequence,1);

            damagedSites = m.chromosome.getDamagedSites(true, true, true, true, false, false, true);
            this.helpTestGeneEssentiality(...
                {},...
                @(m,i) nnz(m.chromosome.getDamagedSites(true, true, true, true, false, false, true)) > nnz(damagedSites));
        end
    end

    %helper methods
    methods
        function clearDNADamages(this)
            m = this.process;

            m.chromosome.initialize();
        end

        function sequence = reactedSequences(this, subs, rxIdx)
            m = this.process;
            positions = mod(...
                repmat(subs(:,1), 1, length(m.reactionVulnerableMotifs{rxIdx})) + ...
                repmat(2*mod(subs(:,2),2)-1, 1, length(m.reactionVulnerableMotifs{rxIdx})) .* ...
                repmat(0:length(m.reactionVulnerableMotifs{rxIdx})-1, size(subs,1), 1) ...
                -1, size(m.chromosome.sequence,1))+1;
            strands = repmat(subs(:,2), 1, length(m.reactionVulnerableMotifs{rxIdx}));

            sequence = reshape(unique(m.chromosome.sequence.subsequence(positions, strands),'rows'), size(m.reactionVulnerableMotifs{rxIdx}));
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [model, B, elapse] = AGH2_learn(A, maxbits, Anchor,s)
%   This is a function of Two Layer AGH (Anchor Graph Hashing) learning.
%
%	Usage:
%	[model, B,elapse] = AGH2_learn(A, maxbits, Anchor,s)
%
%	      A: Rows of vectors of data points. Each row is sample point
%   maxbits: Code length
%    Anchor: Anchors (landmarks), Each row is sample point
%         s: Number of nearest anchor to learn a representation
%            of the sample vector
%
%     model: Used for encoding a test sample point.
%	      B: The binary code of the input data A. Each row is sample point
%    elapse: The coding time (training time).
%
%
%
%   version 2.0 --Nov/2016
%   version 1.0 --Jan/2014
%
%   Written by  Yue Lin (linyue29@gmail.com)
%               Deng Cai (dengcai AT gmail DOT com)
%

tmp_T = tic;

options.CodingMethod = 'Gaussian';

nAnchor = 1500;
if ~exist('Anchor','var')
    [~,Anchor]=litekmeans(A,nAnchor,'MaxIter',5,'Replicates',1);
end
if ~exist('s','var')
    s = 50;
end

[B, W, Thres, sigma] = TwoLayerAGH_Train(A, Anchor, maxbits, s, 0, options);

model.W = W;
model.Thres = Thres;
model.sigma = sigma;
model.Anchor = Anchor;
model.s = s;
model.options = options;

elapse = toc(tmp_T);
end

author: ZJULearning

repo: MatlabFunc


text:
function [PatientPositionCODE] = dicomrt_getPatientPosition(study)
% dicomrt_getPatientPosition(study)
%
% Get Patient Position
% Return a number (CODE) which correspond to one of the supported cases
%
% Patient Position codes:
%
% Head First Supine (HFS) - CODE=1
% Feet First Supine (FFS) - CODE=2
% Head First Prone  (HFP) - CODE=3
% Feet First Prone  (FFP) - CODE=4
%
% Example:
%
% [B]=dicomrt_getPatientPosition(A)
%
% if patient position is HFS, dicomrt_getPatientPosition returns "1" in B.
%
% See also: dicomrt_loaddose, dicomrt_ctcreate, dicomrt_createwphantom
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

if iscell(study) == 1 % this is a rtplan warmed for MC export: Patient Position already assigned
    if size(study,2) > 2
        CODE=study{1,2}(9);
    else % this is CT data
        study=study{1,1};
    end
end

if isstruct(study) == 1
    if strcmpi(study.Modality,'RTPLAN')
        try
            PatientPosition=getfield(study,'PatientSetupSequence','Item_1','PatientPosition');
        catch
            PatientPosition=input('dicomrt_getPatientPosition: Please specify Patient Position: HFS(default),FFS,HFP,FFP: ','s');
            if isempty(PatientPosition)==1
                PatientPosition='HFS';
            end
        end
    elseif strcmpi(study.Modality,'CT')
        PatientPosition=getfield(study,'PatientPosition');
    else
        error('dicomrt_getPatientPosition: could not retrieve PatientPosition field. Exit now')
    end

    if strcmpi(PatientPosition, 'HFS')
        PatientPositionCODE = 1;
    elseif strcmpi(PatientPosition, 'FFS')
        PatientPositionCODE = 2;
    elseif strcmpi(PatientPosition, 'HFP')
        PatientPositionCODE = 3;
    elseif strcmpi(PatientPosition, 'FFP')
        PatientPositionCODE = 4;
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function epoch= getBestEpoch(recalls, recallNs, N)
    if nargin<3, N= 5; end

    % This command is generally OK
    % [~, epoch]= max(recalls(recallNs==N,:));
    % but the following does tie-breaking

    nRecalls= length(recallNs);
    nEpochs= size(recalls, 2);
    assert( nRecalls==size(recalls, 1) );

    posN= find( recallNs==N, 1);
    assert(~isempty(posN));
    posNs= posN:-1:1;
    if posN<nRecalls
        posNs= [posNs, (posN+1):nRecalls];
    end

    potential= true(1, nEpochs);
    for posN= posNs
        maxVal= max(recalls(posN, potential));
        isMax= abs(recalls(posN,:)-maxVal)<1e-6;
        potential= potential & isMax;
        if sum(potential)<=1
            break;
        end
    end

    epoch= find( potential, 1 );
    assert(~isempty(epoch));
end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
function runSimSCFDMA()

SP.FFTsize = 512;
SP.inputBlockSize = 16;
SP.CPsize = 20;
%SP.subband = 15;
SP.subband = 0;

SP.SNR = [0:2:20];
SP.numRun = 10^5;

% TS 25.104
pedAchannel = [1 10^(-9.7/20) 10^(-22.8/20)];
pedAchannel = pedAchannel/sqrt(sum(pedAchannel.^2));
vehAchannel = [1 0 10^(-1/20) 0 10^(-9/20) 10^(-10/20) 0 0 0 10^(-15/20) 0 0 0 10^(-20/20)];
vehAchannel = vehAchannel/sqrt(sum(vehAchannel.^2));
idenChannel = 1;

SP.channel = idenChannel;
%SP.channel = pedAchannel;
%SP.channel = vehAchannel;

SP.equalizerType ='ZERO';
%SP.equalizerType ='MMSE';

[SER_ifdma SER_lfdma] = scfdma(SP);

save scfdma_awgn


text:
function [ahm, AHM_rf, AHM_sf, AHM_sk, AHM_sc, AHM_u, AHM_rho] = ...
    retroProjAhmPntFromPinHoleOnRob(Rf, Sf, Sk, Sc, u, n)

% RETROPROJAHMPNTFROMPINHOLEONROB Retro-project ahm from pinhole on robot.
%
%   AHM = RETROPROJAHMPNTFROMPINHOLEONROB(RF, SF, SK, SC, U, N) gives the
%   retroprojected AHM in World Frame from an observed pixel U. RF and SF
%   are Robot and Sensor Frames, SK and SD are camera calibration and
%   distortion correction parameters. U is the pixel coordinate and N is
%   the non-observable inverse depth. AHM is a 7-vector :
%     AHM = [X Y Z U V W IDepth]'
%
%   [AHM, AHM_rf, AHM_sf, AHM_k, AHM_c, AHM_u, AHM_n] = ... returns the
%   Jacobians wrt RF.x, SF.x, SK, SC, U and N.
%
%   See also INVPINHOLEAHM, FROMFRAMEAHM.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

% Frame World -> Robot  :  Rf
% Frame Robot -> Sensor :  Sf

if(isempty(Sc))
    % AHM in Sensor Frame
    [ahms, AHMS_u, AHMS_rho, AHMS_sk] = invPinHoleAhm(u,n,Sk) ;
else
    % AHM in Sensor Frame
    [ahms, AHMS_u, AHMS_rho, AHMS_sk, AHMS_sc] = invPinHoleAhm(u,n,Sk,Sc) ;
end

[ahmr, AHMR_sf, AHMR_ahms] = fromFrameAhm(Sf,ahms);
[ahm , AHM_rf , AHM_ahmr]  = fromFrameAhm(Rf,ahmr);

AHM_ahms = AHM_ahmr*AHMR_ahms;
AHM_sk   = AHM_ahms*AHMS_sk ;
AHM_sf   = AHM_ahmr*AHMR_sf;

if(isempty(Sc))
    AHM_sc = [] ;
else
    AHM_sc = AHM_ahms*AHMS_sc ;
end

AHM_u = AHM_ahms*AHMS_u ;
AHM_rho = AHM_ahms*AHMS_rho ;

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function pdf = bernoulli_pdf ( x, a )

%*****************************************************************************80
%
%% BERNOULLI_PDF evaluates the Bernoulli PDF.
%
%  Discussion:
%
%    PDF(X)(A) = A**X * ( 1.0D+00 - A )**( X - 1 )
%
%    X = 0 or 1.
%
%    The Bernoulli PDF describes the simple case in which a single trial
%    is carried out, with two possible outcomes, called "success" and
%    "failure"; the probability of success is A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the number of successes on a single trial.
%    X = 0 or 1.
%
%    Input, real A, the probability of success on one trial.
%    0.0D+00 <= A <= 1.0.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < 0 )
    pdf = 0.0;
  elseif ( x == 0 )
    pdf = 1.0 - a;
  elseif ( x == 1 )
    pdf = a;
  else
    pdf = 0.0;
  end

  return
end