text:
%IMPEXP_TRIANGLE Import/export grid data in Triangle format.
%
%   [ DATA ] = IMPEXP_TRIANGLE( FILE_NAME, MODE, DATA, HMAX, HMAX_E, FID_LOG ) Import
%   or export of Triangle grid data in .node/.ele/.edge/.poly format. FILE_NAME is a string
%   specifying the (root) file name to process. MODE can either be a string indicating
%   import or export. For export, DATA can be either a whole fea struct or just the geom
%   struct which is used to compute the poly information used by Triangle. HMAX and HMAX_E
%   specifies grid sizes for subdomains and edges when exporting. A DATA grid struct is
%   output when importng. FID_LOG is an optional log file for message output.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function Y = clfKnnFwd( clf, X )
% Apply a k-nearest neighbor classifier to X.
%
% USAGE
%  Y = clfKnnFwd( clf, X )
%
% INPUTS
%  clf     - trained model
%  X       - nxp data array
%
% OUTPUTS
%  Y       - nx1 vector of labels predicted according to the clf
%
% EXAMPLE
%
% See also CLFKNN, CLFKNNTRAIN
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2008 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

if( ~strcmp(clf.type,'knn')); error( ['incorrect type: ' clf.type] ); end
if( size(X,2)~= clf.p ); error( 'Incorrect data dimension' ); end

metric = clf.metric;
Xtrain = clf.Xtrain;
Ytrain = clf.Ytrain;
k = clf.k;

% get nearest neighbors for each X point
D = pdist2( X, Xtrain, metric );
Y = clfKnnDist( D, Ytrain, k );

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function model = gpdisimCreate(numGenes, numProteins, times, geneVals, ...
			       geneVars, options, annotation)

% GPDISIMCREATE Create a GPDISIM model.
% The GPSIM model is a model for estimating the protein
% concentration in a small gene network where several genes are
% governed by one protein. The model is based on Gaussian processes
% and simple linear differential equations of the form
%
% dx(t)/dt = B + Cf(t) - Dx(t)
%
% where x(t) is a given genes concentration and f(t) is the protein
% concentration.
%
% FORMAT
% DESC creates a model for single input motifs with Gaussian
% processes.
% ARG numGenes : number of genes to be modelled in the system.
% ARG numProteins : number of proteins to be modelled in the
% system.
% ARG times : the time points where the data is to be modelled.
% ARG geneVals : the values of each gene at the different time points.
% ARG geneVars : the varuabces of each gene at the different time points.
% ARG options : options structure, the default options can be
% generated using gpsimOptions.
% ARG annotation : annotation for the data (gene names, etc.) that
% is stored with the model. (Optional)
% RETURN model : model structure containing default
% parameterisation.
%
% SEEALSO : modelCreate, gpsimOptions
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007

% SHEFFIELDML

if any(size(geneVars)~=size(geneVals))
  error('The gene variances have a different size matrix to the gene values.');
end

if(numGenes ~= (size(geneVals, 2) - 1))
  error('The number of genes given does not match the dimension of the gene values given.')
end

if(size(times, 1) ~= size(geneVals, 1))
  error('The number of time points given does not match the number of gene values given')
end

model.type = 'gpdisim';

kernType1{1} = 'multi';
kernType2{1} = 'multi';
kernType1{2} = 'rbf';
for i = 1:numGenes
  kernType1{i+2} = 'disim';
end
tieParam = {'di_decay', 'inverse width', 'di_variance', 'rbf(_| . )variance'};

model.y = geneVals(:);
model.yvar = geneVars(:);

model.includeNoise = options.includeNoise;

if model.includeNoise
  model.yvar = zeros(size(geneVars(:)));
else
  model.yvar = geneVars(:);
end

% Check if we have a noise term.
if model.includeNoise
  % Create a new multi kernel to contain the noise term.
  kernType2{1} = 'multi';

  % Set the new multi kernel to just contain 'white' kernels.
  for i = 1:numGenes+1
    kernType2{i+1} = 'white';
  end
  if isfield(options, 'singleNoise') & options.singleNoise
    tieParam{5} = 'white \d+ variance';
  end

  % Now create model with a 'cmpnd' (compound) kernel build from two
  % multi-kernels. The first multi-kernel is the sim-sim one the next
  % multi-kernel is the white-white one.
  model.kern = kernCreate(times, {'cmpnd', kernType1, kernType2});
  simMultiKernName = 'model.kern.comp{1}';
else
  model.kern = kernCreate(times, kernType1);
  simMultiKernName = 'model.kern';
end
simMultiKern = eval(simMultiKernName);

% This is if we need to place priors on parameters ...
if isfield(options, 'addPriors') && options.addPriors,
  for i = 1:length(simMultiKern.numBlocks)
    % Priors on the sim kernels.
    eval([simMultiKernName '.comp{i}.priors = priorCreate(''gamma'');']);
    eval([simMultiKernName '.comp{i}.priors.a = 1;']);
    eval([simMultiKernName '.comp{i}.priors.b = 1;']);
    %model.kern.comp{i}.priors = priorCreate('gamma');
    %model.kern.comp{i}.priors.a = 1;
    %model.kern.comp{i}.priors.b = 1;
    if i == 1
      % For first kernel place prior on inverse width.
      % model.kern.comp{i}.priors.index = [1 2];
      eval([simMultiKernName '.comp{i}.priors.index = [1 2];']);
    elseif i == 2
      %model.kern.comp{i}.priors.index = [1 3 4 5];
      eval([simMultiKernName '.comp{i}.priors.index = [1 3 4 5];']);
    else
      % For other kernels don't place prior on inverse width --- as
      % they are all tied together and it will be counted multiple
      % times.
      %model.kern.comp{i}.priors.index = [4 5];
      eval([simMultiKernName '.comp{i}.priors.index = [4 5];']);
    end
  end

  % Prior on the b values.
  model.bprior = priorCreate('gamma');
  model.bprior.a = 1;
  model.bprior.b = 1;
end

model.kern = modelTieParam(model.kern, tieParam);
if model.includeNoise,
  for i = 1:numGenes+1,
    model.kern.comp{2}.comp{i}.variance = 1e-2;
  end
end

% The decays and sensitivities are actually stored in the kernel.
% We'll put them here as well for convenience.
model.delta = 10;
model.sigma = 1;
for i = 2:simMultiKern.numBlocks
  eval([simMultiKernName '.comp{i}.di_decay = model.delta;']);
  eval([simMultiKernName '.comp{i}.di_variance = model.sigma^2;']);
  model.D(i-1) = simMultiKern.comp{i}.decay;
  model.S(i-1) = sqrt(simMultiKern.comp{i}.variance);
end

rand('seed',0);
model.numParams = numGenes + model.kern.nParams;
model.numGenes = numGenes;
model.mu = mean(geneVals(:, 2:end));
% model.B = model.D.*model.mu;
model.B = model.D.*geneVals(1, 2:end);
model.m = model.y;
model.t = times;

model.optimiser = options.optimiser;

if isfield(options, 'fix')
  model.fix = options.fix;
end

% The basal transcriptions rates must be postitive.
model.bTransform = optimiDefaultConstraint('positive');

if nargin > 6,
  model.annotation = annotation;
end

model.options = options;

% This forces kernel compute.
params = gpdisimExtractParam(model);
model = gpdisimExpandParam(model, params);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [dir,x0,y0] = boundarydir(x,y,orderout)
%BOUNDARYDIR Determine the direction of a sequence of planar points.
%   DIR = BOUNDARYDIR(X,Y) determines the direction of travel of a
%   closed, nonintersecting sequence of planar points with coordinates
%   contained in column vectors X and Y. Values of DIR are 'cw'
%   (clockwise) and 'ccw' (counterclockwise). The direction of travel is
%   with respect to the image coordinate system defined in Chapter 2 of
%   the book.
%
%   [DIR,X0,Y0] = BOUNDARYDIR(X,Y,ORDEROUT) determines the direction DIR
%   of the input sequence, and also outputs the sequence with its
%   direction of travel as specified in ORDEROUT. Valid values of this
%   parameter as 'cw' and 'ccw'. The coordinates of the output sequence
%   are column vectors X0 and Y0.
%
%   The input sequence is assumed to be nonintersecting, and it cannot
%   have duplicate points, with the exception of the first and last
%   points possibly being the same, a condition often resulting from
%   boundary-following functions, such as bwboundaries.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Preliminaries.
% Make sure coordinates are column vectors.
x = x(:);
y = y(:);

% If the first and last points are the same, delete the last point.
% The point will be restored later.
restore = false;
if x(1) == x(end) && y(1) == y(end)
   x = x(1:end-1);
   y = y(1:end-1);
   restore = true;
end
% Check for duplicate points.
if length([x y]) ~= length(unique([x y],'rows'))
   error('No duplicate points except first and last are allowed.')
end

% The topmost, leftmost point in the sequence is always a convex
% vertex.
x0 = x;
y0 = y;
cx = find(x0 == min(x0));
cy = find(y0 == min(y0(cx)));
x1 = x0(cx(1));
y1 = y0(cy(1));
% Scroll data so that the first point in the sequence is (x1, y1),
% the guaranteed convex point.
I = find(x0 == x1 & y0 == y1);
x0 = circshift(x0, [-(I - 1), 0]);
y0 = circshift(y0, [-(I - 1), 0]);

% Form the matrix needed to check for travel direction. Only three
% points are needed: (x1, y1), the point before it, and the point
% after it.
A = [x0(end) y0(end) 1; x0(1) y0(1) 1; x0(2) y0(2) 1];
dir = 'cw';
if det(A) > 0
   dir = 'ccw';
end

% Prepare outputs.
if nargin == 3
   x0 = x; % Reuse x0 and y0.
   y0 = y;
   if ~strcmp(dir,orderout)
      x0(2:end) = flipud(x0(2:end)); % Reverse order of travel.
      y0(2:end) = flipud(y0(2:end));
   end
   if restore
      x0(end + 1) = x0(1);
      y0(end + 1) = y0(1);
   end
end

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
function [err, errK] = getL2error3CR(node,elem,uvecexact,uh,quadOrder)
%% GETL2ERROR3CR L2 norm of the approximation error for vectorial CR element
%    need Matlab ver > 2016b for the vectorization to work
%
%    Reference: Error analysis of a decoupled finite element method for quad-curl problems
%               https://arxiv.org/abs/2102.03396
%
%    see alos: getL2error3
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('quadOrder','var'); quadOrder = 3; end
elem2face = dof3face(elem);

%%
NT = size(elem,1);
err = zeros(NT,1);
[lambda,weight] = quadpts3(quadOrder);
phi = 1-3*lambda;

nQuad = size(lambda,1);
for p = 1:nQuad
    uhp = uh(elem2face(:,1),:)*phi(p,1) + ...
        uh(elem2face(:,2),:)*phi(p,2) + ...
        uh(elem2face(:,3),:)*phi(p,3) + ...
        uh(elem2face(:,4),:)*phi(p,4);
    % quadrature points in the x-y coordinate
    pxyz = lambda(p,1)*node(elem(:,1),:) ...
         + lambda(p,2)*node(elem(:,2),:) ...
         + lambda(p,3)*node(elem(:,3),:) ...
         + lambda(p,4)*node(elem(:,4),:);
    err = err + weight(p)*sum((uvecexact(pxyz) - uhp).^2,2);
end

%%
d12 = node(elem(:,2),:)-node(elem(:,1),:);
d13 = node(elem(:,3),:)-node(elem(:,1),:);
d14 = node(elem(:,4),:)-node(elem(:,1),:);
volume = abs(dot(mycross(d12,d13,2),d14,2)/6);
err(isnan(err)) = 0; % singular point is excluded
err = volume.*err;
if nargout > 1; errK = sqrt(err); end
err = sqrt(sum(err));

end

repo: ifem


text:
function  [ id, type, rep, field, symm ] = mm_header_read ( input_unit )

%*****************************************************************************80
%
%% MM_HEADER_READ reads the header line from a Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the input unit identifier.
%
%    Output, character ( len = 14 ) ID, the Matrix Market identifier.
%    This value must be '%%MatrixMarket'.
%
%    Output, character ( len = 6 ) TYPE, the Matrix Market type.
%    This value must be 'matrix'.
%
%    Output, character ( len = 10 ) REP, the Matrix Market 'representation'
%    indicator.  Possible values include:
%    'coordinate'   (for sparse data)
%    'array'        (for dense data)
%    'elemental'    (to be added)
%
%    Output, character ( len = 7 ) FIELD, the Matrix Market 'field'.
%    Possible values include:
%    'real'
%    'double'
%    'complex'
%    'integer'
%    'pattern'
%
%    Output, character ( len = 19 ) SYMM, the Matrix Market symmetry.
%    Possible values include:
%    'symmetric'
%    'hermitian'
%    'skew-symmetric'
%    'general'
%
  header = fgets ( input_unit );

  if ( header == -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MM_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  An empty header line was read.\n' );
    error ( 'MM_HEADER_READ - Fatal error!' )
  end

  [ id,    header ] = strtok ( header );
  [ type,  header ] = strtok ( header );
  [ rep,   header ] = strtok ( header );
  [ field, header ] = strtok ( header );
  [ symm,  header ] = strtok ( header );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%Transcription
%
% @wholeCellModelID Process_Transcription
% @name             Transcription
% @description
%   Biology
%   ===============
%   Transcription is the first step in the synthesis of functional gene
%   products where RNA polymerase and several accessory enzymes translate
%   transcription units, or regions of the DNA containing 1 or more genes, into
%   RNA molecules. Following transcription the RNA molecules follow one of
%   two pathways:
%   - mRNAs are used as templates for translation (see translation process)
%   - r/s/tRNAs which are transcribed as molecules containing multiple genes are
%     cleavaged into their individual genes (see RNA processing process), are
%     modified at several bases (see RNA modification process) to improve their
%     stability and enhance their catalytic activity, and finally act as
%     ribozymes (r/sRNAs) or as adaptors between mRNAs and the amino acids they
%     code for (tRNAs).
%
%   Transcription begins with the recruitment of RNA polymerase to a promoter
%   with the help of the sigma initiation factor and possiblity transcription
%   factors. Next elongation factors are recruited, RNA begins to be
%   polymerized, and sigma factor is released. Finally the RNA polymerase
%   reaches  a terminator at the the end of the transcription unit, and
%   with the help of termination factors releases the polymerized RNA and
%   dissociates from the DNA. Termination in E. coli occurs via either a
%   Rho-dependent (50%) or Rho-independent mechanism (50%). Rho-dependent
%   termination is catalyzed by the hexameric ATP-dependent helicase Rho.
%   Rho is not essential in B. subtilis [PUB_0234]. Rho-independent
%   termination occur via the intrinsic properties RNA which disrupt RNA
%   polymerase-DNA binding. Terminator hairpins are not predicted in M.
%   genitalium. In E. coli termination is incompetition with
%   antitermination. However antitermination has not been reported in any
%   mcyoplasma [PUB_0182].
%
%   As soon as the RNA begins to polymerized, even prior to termination, the
%   mRNA transcripts may be bound by ribosomes and polymerized. For simplicity,
%   our model doesn't represent this phenomenon.
%
%   Knowledge Base
%   ===============
%   The transcription unit structure was compiled from several sources:
%   - Primary reports of cotranscribed genes [PUB_0176, PUB_0182, PUB_0186,
%     PUB_0188, PUB_0188, PUB_0244, PUB_0247, PUB_0248, PUB_0249, PUB_0251]
%   - OperonDB database of cotranscribed genes [PUB_0250]
%   - Conservation of gene order across multiple species
%   - Related function of adjacent genes
%   - Expression levels as measured by microarrays of adjacent genes [PUB_0569]
%   - Strandedness of adjacent genes
%   - Weiner, Hermann, and Browning computational model of promoters and
%     transcription unit start sites [PUB_0411]
%
%   The transcription unit structure is organized in the knowledge base, and is
%   loaded into this process by the initializeConstants method.
%
%   The knowledge base also contains the measured expression and half-lives of
%   many transcripts. These values are loaded by initializeConstants and fit by
%   simulation.fitConstants to be consistent with other processes.
%
%   Representation
%   ===============
%   substrates, enzymes, boundEnzymes, and RNAs represent the counts of
%   metabolites, free transcription enzymes, transcription enzymes bound to RNAs
%   and RNA polymerases, and nascent RNAs.
%
%   rnaPolymerases.states represents the current state / pseudostate (actively
%   transcribing, specifically bound, non-specifically bound, free,
%   non-existent) of each RNA polymerase, where each state is indicated by the
%   enumeration:
%   - rnaPolymeraseActivelyTranscribingValue
%   - rnaPolymeraseSpecificallyBoundValue
%   - rnaPolymeraseNonSpecificallyBoundValue
%   - rnaPolymeraseFreeValue
%   - rnaPolymeraseNotExistValue (state exists as a way to account for memory
%     allocated for future RNA polymerases)
%   For actively transcribing polymerases rnaPolymerases.states also represents
%   the position of the polymerase along the transcription unit.
%
%   That is entries of rnaPolymerases.states with the following values
%   corresponding to these states:
%     >= RNAPolymerases.activelyTranscribingValue: RNA polymerases position on genome actively transcribing
%     == RNAPolymerases.specificallyBoundValue:    RNA polymerase specifically bound
%     == RNAPolymerases.nonSpecificallyBoundValue: RNA polymerase non-specifically bound
%     == RNAPolymerases.freeValue:                 RNA polymerase free
%     == RNAPolymerases.notExistValue:             RNA polymerase doesn't exist
%
%   transcripts.boundTranscriptionUnits represents the particular transcription
%   unit to which each actively transcribing and specifically bound polymerase
%   is bound.
%
%   rnaPolymeraseStateExpectations represents the expected occupancies of the
%   RNA polymerase states.
%
%   transcriptionUnitBindingProbabilities represents the relative affinity of
%   RNA polymerases for the promoters of each transcription unit.
%   transcriptionFactorBindingProbFoldChange represents the fold change affect of
%   transcription factors on the relative affinities of the RNA polymerse for
%   the promoters. RNA polymerases are assigned to
%   transcription units weighted by the product of
%   transcriptionUnitBindingProbabilities and
%   transcriptionFactorBindingProbFoldChange.
%
%   Initialization
%   ===============
%   All RNAs are initialized to the mature state. This is implemented by the
%   simulation class initializeState method.
%
%   RNA polymerases are initialized to their steady state:
%   - Each RNA polymerase is randomly assigned (with replacement) to one of the
%     actively transcribing, specifically bound, non-specifically bound, or free
%     states weighted by the expected occupancy of each state
%     (rnaPolymeraseStateExpectations)
%   - Actively transcribing and specifically bound polymerases randomly assigned
%     to transcription units weighted by their transcription rates
%     (transcriptionUnitBindingProbabilities)
%   - Each transcription unit to which an actively transcribing polymerase has
%     been assigned is divided into 1 segment for each polymerase
%   - Actively transcribing polymerases randomly assigned to positions within
%     the assigned segment of their assigned transcription unit (positions near
%     the segment border are not allowed to prevent polymerases from being too
%     close to each other) with uniform probably.
%
%   Simulation
%   ===============
%   Evolves the state of RNA polymerase using a markov chain model with four
%   states:
%   - actively translating
%   - specifically bound
%   - non-specifically bound
%   - free
%
%   Transition probabilities are designed to maintain the occupancy of each
%   state within a narrow window around their expected values. Transition
%   probability are determined by four logistic control functions. These can
%   be tuned with the constants
%   - rnaPolymeraseStateExpectations
%
%   RNA polymerase are created in the free state.
%
%   Actively transcribing state:
%   1. Release sigma factor if after first second of elongation
%   2. Elongate transcript according to nucleic acid limits (substrates)
%      if elongation factors are available
%   3. If transcription complete and termination factor available
%      - release transcript
%      - transition RNA polymerase to free state
%      - increment gene expression
%      Otherwise remain in active state
%
%   Specifically bound State:
%   - Can transition to active, specifically bound, non-specifically bound,
%     or free states
%   - Transition into state only if a free sigma factor is available
%     1. Decrement number of free sigma factors
%     2. Pick a transcription unit (tu) to bind to according to
%     Expression transcription unit i~prob(ribosome releases tu i|ribosome active)
%                      =prob(ribosome within RNA polymerase elongation rate bases of length of tu i|ribosome active)
%                      =prob(ribosome within RNA polymerase elongation rate bases of length of tu i|ribosome active, bound to tu i)*prob(ribosome bound to tu i|ribosome active)
%                      =[(RNA polymerase transcription rate)/(length of tu i)] * [(length of tu i)*prob(binding tu i|binding)]
%     prob(binding tu i | binding)~expression tu i
%
%   Non-specifically bound state:
%   - Can transition to specifically bound, non-specifically bound, or free
%     states
%
%   Free state:
%   - Can transition to specifically bound, non-specifically bound, or free
%     states
%
%   Algorithm
%   +++++++++++++++
%   1. Randomly transition RNA polymerases among activlely transcribing,
%      specifically bound, non-specifically, bound, and free states weighted by
%      state transition probabilities. Update rnaPolymerases.states.
%   2. Randomly assign RNA polymerases entering the specifically bound to
%      specific transcription units weighted by the product of
%      transcriptionUnitBindingProbabilities and
%      transcriptionFactorBindingProbFoldChange. Update
%      transcripts.boundTranscriptionUnits.
%   3. Assign RNA polymerase entering the actively transcribing state sigma
%      factors. Update enzymes and boundEnzymes.
%   4. Simulate RNA polymerization by actively transcribing RNA polymerases with
%      the aid of elongation factors. Allocate available nucleic acids among the
%      actively transcribing RNA polymerases. Release sigma factors from RNA
%      polymerases that started at the beginning of the transcription unit and
%      progressed. Update rnaPolymerases.states. Update substrates. Update enzymes
%      and boundEnzymes.
%   5. If termination factors are available dissociate RNA polymerases which
%      have reached the terminus of the transcription they're bound to, and
%      release RNAs. Update rnaPolymerases.states and
%      transcripts.boundTranscriptionUnits. Increment RNAs.
%
%   References
%   ===========
%   1. McClure, W. R. 1985. Mechanism and control of transcription
%      initiation in prokaryotes. Annu. Rev. Biochem. 54:171-204.
%      [PUB_0775]
%   2. Ciampi MS (2006). Rho-dependent terminators and transcription
%      termination. Microbiology. 152(9):2515-28 [PUB_0233].
%   3. Nudler E, Gottesman ME (2002). Transcription termination and
%      anti-termination in E. coli. Genes Cells. 7(8):755-68. [PUB_0662]
%   4. Washio T, Sasayama J, Tomita M (1998). Analysis of complete genomes
%      suggests that many prokaryotes do not rely on hairpin formation in
%      transcription termination. Nucleic Acids Res. 26(23):5456-63
%      [PUB_0234]
%   5. Peterson JD, Umayam LA, Dickinson T, Hickey EK, White O (2001). The
%      Comprehensive Microbial Resource. Nucleic Acids Res. 29(1):123-5.
%      [PUB_0182]
%   6. Shepherd N, Dennis P, Bremer H (2001). Cytoplasmic RNA Polymerase in
%      Escherichia coli. J Bacteriol. 183(8): 2527-34. [PUB_0784]
%   7. Klumpp S, Hwa T (2008). Growth-rate-dependent partitioning of RNA
%      polymerases in bacteria. Proc Natl Acad Sci U S A. 105(21):
%      20245-50. [PUB_0785]
%   8. Grigorova IL, Phleger NJ, Mutalik VK, Gross CA (2006). Insights into
%      transcriptional regulation and sigma competition from an equilibrium
%      model of RNA polymerase binding to DNA. Proc Natl Acad Sci U S A.
%      103(14): 5332-7. [PUB_0786]
%
% Author: Markus Covert, mcovert@stanford.edu
% Author: Jayodita Sanghvi, jayodita@stanford.edu
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 1/4/2010
%
%TODO: require Mg2+ or Mn2+ as cofactor for transcription
classdef Transcription < edu.stanford.covert.cell.sim.Process & edu.stanford.covert.cell.sim.ChromosomeProcessAspect
    %property annotations
    properties (Constant)
        optionNames__              = {}; %names of option properties
        fixedConstantNames__       = {   %names of fixed constant properties
            'rnaPolymeraseElongationRate';
            'transcriptionUnitBaseCounts';
            'transcriptionUnitBindingProbabilities';
            };
        fittedConstantNames__      = {   %names of fitted constant properties
            'transcriptionUnitBindingProbabilities'
            };
        localStateNames__          = {   %names of simulation state properties redundant with timecourses in this or other processes or the simulation
            'RNAs'};
    end

    %IDs, names, and local indices
    properties
        stimuliWholeCellModelIDs = {}; %whole cell model IDs of stimuli

        substrateWholeCellModelIDs = { %whole cell model IDs of substrates
            'ATP'; 'CTP'; 'GTP'; 'UTP';
            'AMP'; 'CMP'; 'GMP'; 'UMP';
            'ADP'; 'PPI'; 'H2O'; 'H'};
        substrateIndexs_ntp         = (1:4)'; %indices within substrates of NTPs
        substrateIndexs_nmp         = (5:8)'; %indices within substrates of NMPs
        substrateIndexs_atp         = 1;      %index within substrates of ATP
        substrateIndexs_adp         = 9;      %index within substrates of ADP
        substrateIndexs_diphosphate = 10;     %index within substrates of diphosphate
        substrateIndexs_water       = 11;     %index within substrates of water
        substrateIndexs_hydrogen    = 12;     %index within substrates of hydrogen

        enzymeWholeCellModelIDs = {       %enzyme whole cell model ids
            'MG_249_MONOMER';             %RNA polymerase sigma factor RpoD
            'MG_282_MONOMER';             %transcription elongation factor GreA
            'MG_141_MONOMER';             %transcription termination factor NusA
            'MG_027_MONOMER';             %transcription termination/antitermination protein NusB
            'RNA_POLYMERASE';             %DNA-directed RNA polymerase
            'RNA_POLYMERASE_HOLOENZYME'}; %DNA-directed RNA polymerase holoenzyme
        enzymeIndexs_transcriptionFactors    = (1:4)'; %indices within enzymes of transcription factors
        enzymeIndexs_sigmaFactor             = 1;      %index within enzymes of sigma factor RpoD
        enzymeIndexs_elongationFactor        = 2;      %index within enzymes of elongation factor GreA
        enzymeIndexs_terminationFactor       = 3;      %index within enzymes of termination factor NusA
        enzymeIndexs_antiterminationFactor   = 4;      %index within enzymes of termination/antitermination protein NusB
        enzymeIndexs_rnaPolymerase           = 5;      %index within enzymes of DNA-directed RNA polymerase
        enzymeIndexs_rnaPolymeraseHoloenzyme = 6;      %index within enzymes of DNA-directed RNA polymerase holoenzyme

        complexIndexs_DnaA_ATP                         %indices within protein complexes of DnaA-ATP
        transcriptionUnitIndexs_DnaAR12345Boxes        %indices of transcription units containing the functional DnaA boxes R1-5
    end

    %fixed biological constants
    properties
        rnaPolymeraseElongationRate                   %RNA polymerase elongation rate (50 nucleotides per second per RNA polymerase) [PUB_0562, PUB_0563]
        transcriptionUnitBaseCounts                   %transcription unit base counts
        transcriptionUnitBindingProbabilities         %transcription unit binding probabilities
        stateTransitionProbabilities                  %transition probabilities among RNA polymerase states
    end

    %global state, stored locally for convenience
    properties
        RNAs                                          %copy number of transcription units
    end

    %global state (referenced locally for convenience)
    properties
        rnaPolymerases   %RNA Polymerase state class
        transcripts      %New Transcripts state class
    end

    %constructor
    methods
        function this = Transcription(wholeCellModelID, name)
            this = this@edu.stanford.covert.cell.sim.Process(wholeCellModelID, name);
        end
    end

    %communication between process/simulation
    methods
        %set references to state objects
        function storeObjectReferences(this, simulation)
            this.storeObjectReferences@edu.stanford.covert.cell.sim.Process(simulation);
            this.storeObjectReferences@edu.stanford.covert.cell.sim.ChromosomeProcessAspect(simulation);

            this.rnaPolymerases = simulation.state('RNAPolymerase');
            this.transcripts = simulation.state('Transcript');

            this.states = [this.states; {this.rnaPolymerases; this.transcripts}];
        end

        %initialize constants
        function initializeConstants(this, knowledgeBase, simulation, varargin)
            %call super class method
            this.initializeConstants@edu.stanford.covert.cell.sim.Process(knowledgeBase, simulation, varargin{:});
            this.initializeConstants@edu.stanford.covert.cell.sim.ChromosomeProcessAspect(knowledgeBase, simulation, varargin{:});

            this.transcriptionUnitBaseCounts = reshape([knowledgeBase.transcriptionUnits.baseCount],[],length(knowledgeBase.transcriptionUnits))';
            this.transcriptionUnitBaseCounts = this.transcriptionUnitBaseCounts(:, this.substrateMetaboliteGlobalIndexs);

            this.transcriptionUnitBindingProbabilities  = zeros(size(this.transcripts.transcriptionUnitLengths));

            this.complexIndexs_DnaA_ATP = find(ismember({knowledgeBase.proteinComplexs.wholeCellModelID}', {
                'MG_469_1MER_ATP'      %DnaA-ATP 1mer
                'MG_469_2MER_1ATP_ADP' %DnaA 2mer-(1)ATP-(1)ADP
                'MG_469_2MER_ATP'      %DnaA-ATP 2mer
                'MG_469_3MER_2ATP_ADP' %DnaA 3mer-(2)ATP-(1)ADP
                'MG_469_3MER_ATP'      %DnaA-ATP 3mer
                'MG_469_4MER_3ATP_ADP' %DnaA 4mer-(3)ATP-(1)ADP
                'MG_469_4MER_ATP'      %DnaA-ATP 4mer
                'MG_469_5MER_4ATP_ADP' %DnaA 5mer-(4)ATP-(1)ADP
                'MG_469_5MER_ATP'      %DnaA-ATP 5mer
                'MG_469_6MER_5ATP_ADP' %DnaA 6mer-(5)ATP-(1)ADP
                'MG_469_6MER_ATP'      %DnaA-ATP 6mer
                'MG_469_7MER_6ATP_ADP' %DnaA 7mer-(6)ATP-(1)ADP
                'MG_469_7MER_ATP'}));  %DnaA-ATP 7mer;
            [~, idxs] = ismember({
                'Functional box R1'
                'Functional box R2'
                'Functional box R3'
                'Functional box R4'
                'Functional box R5'}, {knowledgeBase.genomeFeatures.name});
            tus = [knowledgeBase.genomeFeatures(idxs).transcriptionUnits];
            this.transcriptionUnitIndexs_DnaAR12345Boxes = unique([tus.idx])';
        end

        %retrieve state from simulation
        function copyFromState(this)
            this.copyFromState@edu.stanford.covert.cell.sim.Process();

            this.RNAs = this.rna.counts(this.rna.nascentIndexs, this.compartment.cytosolIndexs, :);
        end

        %send state to simulation
        function copyToState(this)
            this.copyToState@edu.stanford.covert.cell.sim.Process();

            this.rna.counts(this.rna.nascentIndexs, this.compartment.cytosolIndexs, :) = this.RNAs;
        end
    end

    %allocate memory for state
    methods
        function allocateMemoryForState(this, numTimePoints)
            this.allocateMemoryForState@edu.stanford.covert.cell.sim.Process(numTimePoints);

            numTranscriptionUnits = length(this.transcripts.transcriptionUnitLengths);

            this.RNAs = zeros(numTranscriptionUnits, 1, numTimePoints);
        end
    end

    %model
    methods
        %Calculate
        %- contribution to FBA objective
        %- minimum expression consistent with cell cycle length
        function [bmProd, byProd, minEnzExp, maxEnzExp] = calcResourceRequirements_LifeCycle(this, ~, states)
            %% initialize
            bmProd = zeros(size(this.substrateWholeCellModelIDs));
            byProd = zeros(size(this.substrateWholeCellModelIDs));
            minEnzExp = zeros(size(this.enzymeWholeCellModelIDs));
            maxEnzExp = Inf(size(this.enzymeWholeCellModelIDs));

            %% substrate and byproducts
            %Contributions, by process, to RNA metabolism:
            %                 NTPs + H2O --(transcription)--> nascent RNA + PPi + H
            %          nascent RNA + H2O --(processing)-----> processed RNA + NMPs
            %processed RNA + metabolites --(modification)---> modified RNA + metabolites
            %               modified RNA --(decay)----------> NMPs + modified NMPs
            %               NMPs + 2 ATP --(charging)-------> NTPs + 2 ADP
            %        modified NMPs + H2O --(catabolism)-----> modified bases + Pi
            %
            %Here we only account for the contributions of transcription.
            invMat = edu.stanford.covert.util.ComputationUtil.invertCompositionMatrix(this.rna.nascentRNAMatureRNAComposition);

            transcribedNTPs = this.transcriptionUnitBaseCounts(:, this.substrateIndexs_nmp)' * invMat * states.rnaProductions;
            bmProd(this.substrateIndexs_ntp)         = bmProd(this.substrateIndexs_ntp)         + transcribedNTPs;
            byProd(this.substrateIndexs_diphosphate) = byProd(this.substrateIndexs_diphosphate) + sum(transcribedNTPs);
            bmProd(this.substrateIndexs_water)       = bmProd(this.substrateIndexs_water)       + sum(states.rnaProductions);
            byProd(this.substrateIndexs_hydrogen)    = byProd(this.substrateIndexs_hydrogen)    + sum(states.rnaProductions);

            %% enzymes
            %RNA polymerase
            fractionInitiating = this.rnaPolymeraseElongationRate * sum(invMat * states.rnaProductions0) / ... %fraction of active RNA polymerases that are initiating
                (this.transcripts.transcriptionUnitLengths' * invMat * states.rnaProductions0);
            minEnzExp(this.enzymeIndexs_rnaPolymerase) = 4.5867 * ...
                (this.transcripts.transcriptionUnitLengths' * invMat * states.rnaProductions0) / ...
                this.rnaPolymeraseElongationRate / ...
                this.rnaPolymerases.stateExpectations(this.rnaPolymerases.activelyTranscribingIndex) / ...
                (1 - fractionInitiating);

            %sigma, elongation, termination factors
            minEnzExp(this.enzymeIndexs_sigmaFactor) = minEnzExp(this.enzymeIndexs_rnaPolymerase) * ...
                (this.rnaPolymerases.stateExpectations(this.rnaPolymerases.activelyTranscribingIndex) * fractionInitiating + ...
                this.rnaPolymerases.stateExpectations(this.rnaPolymerases.specificallyBoundIndex));
            minEnzExp(this.enzymeIndexs_elongationFactor)  = 2;
            minEnzExp(this.enzymeIndexs_terminationFactor) = 2;
        end

        %Initialize RNA polymerase, bound transcription units, transcription factor states
        %Assumptions:
        %- Metabolic cost of the transcription of these nucleic acids is negligible
        %- Probability that so many RNA polymerases bind a given transcription unit
        %  that they can't pack along the transcription unit without steric
        %  interactions is negligible. Thus we don't try to handle this case
        %  separately
        function initializeState(this)
            rnaPols = this.rnaPolymerases;
            trnspts = this.transcripts;

            %% states
            trnsptStrnds = 2 - trnspts.transcriptionUnitDirections;
            trnsptLens = trnspts.transcriptionUnitLengths;
            trnsptDirs = 2 * trnspts.transcriptionUnitDirections - 1;
            trnspt5Coords = trnspts.transcriptionUnitFivePrimeCoordinates;
            trnsptStarts = ...
                trnspt5Coords - (1 - trnspts.transcriptionUnitDirections) .* (trnspts.transcriptionUnitLengths-1);
            trnsptStarts(trnsptDirs == 1) = ...
                trnsptStarts(trnsptDirs == 1) - ...
                this.enzymeDNAFootprints5Prime(this.enzymeIndexs_rnaPolymerase);
            trnsptStarts(trnsptDirs == -1) = ...
                trnsptStarts(trnsptDirs == -1) + ...
                this.enzymeDNAFootprints5Prime(this.enzymeIndexs_rnaPolymerase);

            %% enzymes
            freTFs = this.enzymes(this.enzymeIndexs_transcriptionFactors);
            bndTFs = this.boundEnzymes(this.enzymeIndexs_transcriptionFactors);
            frePls = this.enzymes(this.enzymeIndexs_rnaPolymerase);
            bndPls = this.boundEnzymes(this.enzymeIndexs_rnaPolymerase);
            freHls = this.enzymes(this.enzymeIndexs_rnaPolymeraseHoloenzyme);
            bndHls = this.boundEnzymes(this.enzymeIndexs_rnaPolymeraseHoloenzyme);
            pls = frePls + bndPls + freHls + bndHls;

            %% allocate
            rnaPols.states = repmat(rnaPols.notExistValue, [2 * pls, 1, 1]);
            rnaPols.positionStrands = zeros([2 * pls, 2, 1]);
            trnspts.boundTranscriptionUnits = repmat(trnspts.nullTranscriptValue, [2 * pls, 1, 1]);
            trnspts.boundTranscriptProgress = repmat(trnspts.nullTranscriptValue, [2 * pls, 1, 1]);
            trnspts.boundTranscriptChromosome = repmat(trnspts.nullTranscriptValue, [2 * pls, 1, 1]);

            %% initialize
            pState = rnaPols.stateExpectations;
            probs = this.computeRNAPolymeraseTUBindingProbabilities();
            pBinds = trnspts.transcriptionUnitLengths .* probs(:, 1);
            pSpBinds = pBinds;

            for i = 1:pls
                %check whether conditions necessary for actively transcribing state are met
                if frePls == 0 || ~any(pBinds) || freTFs(this.enzymeIndexs_elongationFactor) == 0
                    pState(rnaPols.activelyTranscribingIndex) = 0;
                end

                %check whether conditions necessary for specifically bound state are met
                if frePls == 0 || ~any(pSpBinds) || freTFs(this.enzymeIndexs_sigmaFactor) == 0
                    pState(rnaPols.specificallyBoundIndex) = 0;
                end

                %check whether conditions necessary for non-specifically bound state are met
                if frePls == 0
                    pState(rnaPols.nonSpecificallyBoundIndex) = 0;
                end

                if ~any(pState)
                    pState(rnaPols.freeIndex) = 1;
                end

                %randomly select state
                state = this.randStream.randsample(4, 1, true, pState);

                %Actively transcribing, specifically bound: randomly select
                %transcription unit
                switch state
                    case rnaPols.activelyTranscribingIndex
                        tmp_pBinds = pBinds;
                        posStrnds = zeros(0, 2);
                        while any(tmp_pBinds)
                            iTU = this.randStream.randsample(numel(pBinds), 1, true, tmp_pBinds);
                            tmp_pBinds(iTU) = 0;
                            [~, posStrnds] = this.bindProteinToChromosomeStochastically(...
                                this.enzymeIndexs_rnaPolymerase, 1, [trnsptStarts(iTU)+1 trnsptStrnds(iTU)], trnsptLens(iTU)-2);
                            if ~isempty(posStrnds)
                                break;
                            end
                        end
                        if isempty(posStrnds)
                            throw(MException('Transcription:error', 'unable to bind chromosome'));
                        end

                        posStrnds( isodd(posStrnds(:, 2)), 1) = posStrnds( isodd(posStrnds(:, 2)), 1) + this.enzymeDNAFootprints5Prime(this.enzymeIndexs_rnaPolymerase);
                        posStrnds(iseven(posStrnds(:, 2)), 1) = posStrnds(iseven(posStrnds(:, 2)), 1) + this.enzymeDNAFootprints3Prime(this.enzymeIndexs_rnaPolymerase);
                        if trnspts.transcriptionUnitDirections(iTU)
                            state = posStrnds(:, 1) - trnspt5Coords(iTU) + 1;
                        else
                            state = trnspt5Coords(iTU) - posStrnds(:, 1) + 1;
                        end

                        rnaPols.states(i) = state;
                        rnaPols.positionStrands(i, :) = posStrnds;
                        trnspts.boundTranscriptProgress(i) = state;
                        trnspts.boundTranscriptionUnits(i) = iTU;
                        trnspts.boundTranscriptChromosome(i) = 1; %initialize to first chromosome

                        frePls = frePls - 1;
                        bndPls = bndPls + 1;

                        %freTFs(this.enzymeIndexs_elongationFactor) = freTFs(this.enzymeIndexs_elongationFactor) - 1;
                        %bndTFs(this.enzymeIndexs_elongationFactor) = bndTFs(this.enzymeIndexs_elongationFactor) + 1;
                    case rnaPols.specificallyBoundIndex
                        while any(pSpBinds)
                            iTU = this.randStream.randsample(numel(pSpBinds), 1, true, pSpBinds);
                            pSpBinds(iTU) = 0;
                            [~, ~, posStrnds] = this.bindProteinToChromosome([trnspt5Coords(iTU)-trnsptDirs(iTU) trnsptStrnds(iTU)], this.enzymeIndexs_rnaPolymeraseHoloenzyme, 1);
                            if ~isempty(posStrnds)
                                break;
                            end
                        end
                        if isempty(posStrnds)
                            throw(MException('Transcription:error', 'unable to bind chromosome'));
                        end

                        rnaPols.states(i) = rnaPols.specificallyBoundValue;
                        rnaPols.positionStrands(i, :) = posStrnds;
                        trnspts.boundTranscriptProgress(i) = trnspts.nullTranscriptValue;
                        trnspts.boundTranscriptionUnits(i) = iTU;
                        trnspts.boundTranscriptChromosome(i) = 1; %initialize to first chromosome

                        bndHls = bndHls + 1;
                        frePls = frePls - 1;
                        freTFs(this.enzymeIndexs_sigmaFactor) = freTFs(this.enzymeIndexs_sigmaFactor) - 1;
                    case rnaPols.nonSpecificallyBoundIndex
                        [~, posStrnds] = this.bindProteinToChromosomeStochastically(this.enzymeIndexs_rnaPolymerase, 1);
                        if isempty(posStrnds)
                            throw(MException('Transcription:error', 'unable to bind chromosome'));
                        end
                        posStrnds( isodd(posStrnds(:, 2)), 1) = posStrnds( isodd(posStrnds(:, 2)), 1) + this.enzymeDNAFootprints5Prime(this.enzymeIndexs_rnaPolymerase);
                        posStrnds(iseven(posStrnds(:, 2)), 1) = posStrnds(iseven(posStrnds(:, 2)), 1) + this.enzymeDNAFootprints3Prime(this.enzymeIndexs_rnaPolymerase);

                        rnaPols.states(i) = rnaPols.nonSpecificallyBoundValue;
                        rnaPols.positionStrands(i, :) = posStrnds;
                        trnspts.boundTranscriptionUnits(i) = trnspts.nullTranscriptValue;
                        trnspts.boundTranscriptProgress(i) = trnspts.nullTranscriptValue;
                        trnspts.boundTranscriptChromosome(i) = trnspts.nullTranscriptValue;

                        frePls = frePls - 1;
                        bndPls = bndPls + 1;
                    case rnaPols.freeIndex
                        rnaPols.states(i) = rnaPols.freeValue;
                        rnaPols.positionStrands(i, :) = 0;
                        trnspts.boundTranscriptionUnits(i) = trnspts.nullTranscriptValue;
                        trnspts.boundTranscriptProgress(i) = trnspts.nullTranscriptValue;
                        trnspts.boundTranscriptChromosome(i) = trnspts.nullTranscriptValue;
                end
            end

            %% store states of enzymes
            this.enzymes(this.enzymeIndexs_transcriptionFactors)      = freTFs;
            this.boundEnzymes(this.enzymeIndexs_transcriptionFactors) = bndTFs;
            this.enzymes(this.enzymeIndexs_rnaPolymerase)             = frePls;
            this.boundEnzymes(this.enzymeIndexs_rnaPolymerase)        = bndPls;
            this.enzymes(this.enzymeIndexs_rnaPolymeraseHoloenzyme)   = freHls;
            this.boundEnzymes(this.enzymeIndexs_rnaPolymeraseHoloenzyme) = bndHls;

            %% decrement counts of RNAs
            this.copyToState();

            comp = this.compartment;
            rna = this.rna;

            matureRnaWt = max(0, sum(rna.dryWeight) - trnspts.dryWeight);
            initRnaCnts = rna.counts;
            while sum(rna.dryWeight) > matureRnaWt
                idx = this.randStream.randsample(size(rna.counts, 1), 1, false, rna.counts(:, comp.cytosolIndexs));
                rna.counts(idx, comp.cytosolIndexs) = rna.counts(idx, comp.cytosolIndexs) - 1;
            end
            rna.counts(:, comp.cytosolIndexs) = ...
                + rna.counts(:, comp.cytosolIndexs) ...
                + rna.updateExternalState(-(initRnaCnts(:, comp.cytosolIndexs) - rna.counts(:, comp.cytosolIndexs)), true);

            this.copyFromState();
        end

        %resource requirements
        function result = calcResourceRequirements_Current(this)
            result = zeros(size(this.substrates));
            result(this.substrateIndexs_ntp) = 2 * ...
                sum(this.metabolite.nmpComposition, 2) * this.rnaPolymerases.nActive * ...
                this.rnaPolymeraseElongationRate;
            result(this.substrateIndexs_water) = this.rnaPolymerases.nActive;
        end

        %simulation
        function evolveState(this)
            %% define and allocate variables

            %states
            c = this.chromosome;
            rnaPols = this.rnaPolymerases;
            trnspts = this.transcripts;

            %numbers of enzymes
            nFreTfs  = this.enzymes(this.enzymeIndexs_transcriptionFactors);
            nBndTFs  = this.boundEnzymes(this.enzymeIndexs_transcriptionFactors);
            nFrePols = this.enzymes(this.enzymeIndexs_rnaPolymerase);
            nBndPols = this.boundEnzymes(this.enzymeIndexs_rnaPolymerase);
            nFreHols = this.enzymes(this.enzymeIndexs_rnaPolymeraseHoloenzyme);
            nBndHols = this.boundEnzymes(this.enzymeIndexs_rnaPolymeraseHoloenzyme);
            nTotPols = nFrePols + nBndPols + nFreHols + nBndHols;

            %properties of RNA polymerases
            pStTrns = this.stateTransitionProbabilities; %probabilities of state transitions

            %transcription unit properties
            tuDirs = 2 * trnspts.transcriptionUnitDirections - 1;
            tuStrnds = c.transcriptionUnitStrands;
            tuLens = trnspts.transcriptionUnitLengths;
            tu5Coords = trnspts.transcriptionUnitFivePrimeCoordinates;
            tuSeqs = trnspts.transcriptionUnitSequences;  %sequences

            %% Update states of RNA polymerases

            %allocate space to store states of new RNA polymerases
            if nTotPols > size(rnaPols.states, 1)
                rnaPols.states = [
                    rnaPols.states;
                    rnaPols.notExistValue(ones(2 * nTotPols - size(rnaPols.states, 1), 1), 1)];
                rnaPols.positionStrands = [
                    rnaPols.positionStrands;
                    zeros(2 * nTotPols - size(rnaPols.positionStrands, 1), 2)];
                trnspts.boundTranscriptionUnits = [
                    trnspts.boundTranscriptionUnits;
                    trnspts.nullTranscriptValue(ones(2 * nTotPols - size(trnspts.boundTranscriptChromosome, 1), 1), 1)];
                trnspts.boundTranscriptProgress = [
                    trnspts.boundTranscriptProgress;
                    trnspts.nullTranscriptValue(ones(2 * nTotPols - size(trnspts.boundTranscriptChromosome, 1), 1), 1)];
                trnspts.boundTranscriptChromosome = [
                    trnspts.boundTranscriptChromosome;
                    trnspts.nullTranscriptValue(ones(2 * nTotPols - size(trnspts.boundTranscriptChromosome, 1), 1), 1)];
            end

            %indices of RNA polymerases in various states
            actPols = find(rnaPols.states >= rnaPols.activelyTranscribingValue);
            sbPols = find(rnaPols.states == rnaPols.specificallyBoundValue);
            freeNsbPols = find(...
                rnaPols.states == rnaPols.nonSpecificallyBoundValue | ...
                rnaPols.states == rnaPols.freeValue);
            sbPols = sbPols(this.randStream.randperm(numel(sbPols)));
            freeNsbPols = freeNsbPols(this.randStream.randperm(numel(freeNsbPols)));

            %number of new RNA polymerases
            nNewPols = nTotPols - (numel(actPols) + numel(sbPols) + numel(freeNsbPols));

            %dissociate any free RNA polymerase holoenzymes (created say by
            %their release from chromosome by another protein) into RNA
            %polymerase core and sigma factor
            nFrePols = nFrePols + nFreHols;
            nFreTfs(this.enzymeIndexs_sigmaFactor) = nFreTfs(this.enzymeIndexs_sigmaFactor) + nFreHols;
            nFreHols = 0;

            %initiating specifically bound polymerases
            nInitMax = this.randStream.stochasticRound(numel(sbPols) * pStTrns(rnaPols.activelyTranscribingIndex, rnaPols.specificallyBoundIndex)); %#ok<*PROP>
            if nInitMax > 0
                releasedProteins = this.releaseProteinFromSites(rnaPols.positionStrands(sbPols(1:nInitMax), :), false, this.enzymeIndexs_rnaPolymeraseHoloenzyme, true, true);
                if ~isequal(releasedProteins(this.enzymeIndexs_rnaPolymeraseHoloenzyme), nInitMax)
                    throw(MException('Transcription:error', 'Unable to release protein'));
                end

                %try to initiate RNA polymerases
                positions = tu5Coords(trnspts.boundTranscriptionUnits(sbPols(1:nInitMax)));
                tfs = this.bindProteinToChromosome([positions rnaPols.positionStrands(sbPols(1:nInitMax), 2)], ...
                    this.enzymeIndexs_rnaPolymeraseHoloenzyme);
                rnaPols.states(sbPols(tfs)) = rnaPols.activelyTranscribingValue;
                trnspts.boundTranscriptProgress(sbPols(tfs)) = trnspts.newTranscriptValue;
                rnaPols.positionStrands(sbPols(tfs), 1) = positions(tfs, 1);

                %rebind RNA polymerases that couldn't move forward
                if ~all(this.bindProteinToChromosome(rnaPols.positionStrands(sbPols(~tfs), :), ...
                        this.enzymeIndexs_rnaPolymeraseHoloenzyme))
                    throw(MException('Transcription:error', 'Unable to unbind protein'));
                end
            end

            %specifically bound polymerases becoming non-specifically bound / free
            nUnsb = max(0, numel(sbPols) - nInitMax - this.randStream.stochasticRound(numel(sbPols) * ...
                pStTrns(rnaPols.specificallyBoundIndex, rnaPols.specificallyBoundIndex)));
            if nUnsb > 0
                idxs = sbPols(end-nUnsb+1:end);
                posStrnds = rnaPols.positionStrands(idxs, :);
                releasedProteins = this.releaseProteinFromSites(posStrnds, false, this.enzymeIndexs_rnaPolymeraseHoloenzyme, true, true);
                if ~isequal(releasedProteins(this.enzymeIndexs_rnaPolymeraseHoloenzyme), nUnsb)
                    throw(MException('Transcription:error', 'Unable to unbind protein'));
                end

                rnaPols.states(idxs) = rnaPols.freeValue;
                rnaPols.positionStrands(idxs, :) = 0;
                trnspts.boundTranscriptionUnits(idxs) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptProgress(idxs) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptChromosome(idxs) = trnspts.nullTranscriptValue;

                nBndHols = nBndHols - nUnsb;
                nFrePols = nFrePols + nUnsb;
                nFreTfs(this.enzymeIndexs_sigmaFactor) = nFreTfs(this.enzymeIndexs_sigmaFactor) + nUnsb;
            end

            %free polymerases/non-specifically bound becoming specifically bound
            nFreeNsbPols = this.randStream.stochasticRound(min([
                numel(freeNsbPols) * pStTrns(rnaPols.specificallyBoundIndex, rnaPols.freeIndex)
                nFreTfs(this.enzymeIndexs_sigmaFactor)]));
            if nFreeNsbPols > 0
                %unbind from non-specifically bound site
                idxs = freeNsbPols(1:nFreeNsbPols);
                nNsbPols = sum(rnaPols.states(idxs) == rnaPols.nonSpecificallyBoundValue);
                releasedProteins = this.releaseProteinFromSites(rnaPols.positionStrands(idxs, :), false, this.enzymeIndexs_rnaPolymerase, true, true);
                if ~isequal(releasedProteins(this.enzymeIndexs_rnaPolymerase), nNsbPols)
                    throw(MException('Transcription:error', 'Unable to release protein'));
                end

                rnaPols.states(idxs) = rnaPols.freeValue;
                rnaPols.positionStrands(idxs, :) = 0;
                trnspts.boundTranscriptionUnits(idxs) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptProgress(idxs) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptChromosome(idxs) = trnspts.nullTranscriptValue;

                nFrePols = nFrePols + nNsbPols;
                nBndPols = nBndPols - nNsbPols;

                %bind to new specifically bound site
                pBnds = this.computeRNAPolymeraseTUBindingProbabilities();

                if nnz(c.polymerizedRegions) == 2
                    [~, iTU, posStrnds, nFreeNsbPols] = this.bindProteinToChromosome( ...
                        [tu5Coords-tuDirs tuStrnds], ...
                        this.enzymeIndexs_rnaPolymeraseHoloenzyme, nFreeNsbPols, pBnds(:, 1), ...
                        true, true, 1, false, []);
                    iChr = 1;
                else
                    [~, iTU, posStrnds, nFreeNsbPols] = this.bindProteinToChromosome( ...
                        [tu5Coords-tuDirs tuStrnds
                        tu5Coords-tuDirs tuStrnds + 2], ...
                        this.enzymeIndexs_rnaPolymeraseHoloenzyme, nFreeNsbPols, pBnds(:), ...
                        true, true, 1, false, []);
                    iTU = mod(iTU - 1, numel(tuLens)) + 1;
                    iChr = ceil(posStrnds(:, 2) / 2);
                end

                idxs = freeNsbPols(1:nFreeNsbPols);
                rnaPols.states(idxs) = rnaPols.specificallyBoundValue;
                rnaPols.positionStrands(idxs, :) = posStrnds;
                trnspts.boundTranscriptionUnits(idxs) = iTU;
                trnspts.boundTranscriptProgress(idxs) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptChromosome(idxs) = iChr;

                nBndHols = nBndHols + nFreeNsbPols;
                nFrePols = nFrePols - nFreeNsbPols;
                nFreTfs(this.enzymeIndexs_sigmaFactor) = nFreTfs(this.enzymeIndexs_sigmaFactor) - nFreeNsbPols;
            end

            %free polymerases/non-specifically bound becoming free/non-specifically bound
            nsbPols = find(rnaPols.states == rnaPols.nonSpecificallyBoundValue);
            if ~isempty(nsbPols)
                %set non-specifically bound to free
                releasedProteins = this.releaseProteinFromSites(rnaPols.positionStrands(nsbPols, :), false, this.enzymeIndexs_rnaPolymerase, true, true);
                if ~isequal(releasedProteins(this.enzymeIndexs_rnaPolymerase), numel(nsbPols))
                    throw(MException('Transcription:error', 'Unable to release protein'));
                end

                rnaPols.states(nsbPols) = rnaPols.freeValue;
                rnaPols.positionStrands(nsbPols, :) = 0;
                trnspts.boundTranscriptionUnits(nsbPols) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptProgress(nsbPols) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptChromosome(nsbPols) = trnspts.nullTranscriptValue;

                nFrePols = nFrePols + numel(nsbPols);
                nBndPols = nBndPols - numel(nsbPols);
            end

            freePols = find(rnaPols.states == rnaPols.freeValue);
            nNsbPols = this.randStream.stochasticRound(min([...
                nFrePols
                numel(freePols) * pStTrns(rnaPols.nonSpecificallyBoundIndex, rnaPols.freeIndex) / (pStTrns(rnaPols.nonSpecificallyBoundIndex, rnaPols.freeIndex) + pStTrns(rnaPols.freeIndex, rnaPols.freeIndex))]));
            if nNsbPols > 0
                %set some to non-specifically bound state
                [nNsbPols posStrnds] = this.bindProteinToChromosomeStochastically(this.enzymeIndexs_rnaPolymerase, nNsbPols);
                posStrnds( isodd(posStrnds(:, 2)), 1) = posStrnds( isodd(posStrnds(:, 2)), 1) + this.enzymeDNAFootprints5Prime(this.enzymeIndexs_rnaPolymerase);
                posStrnds(iseven(posStrnds(:, 2)), 1) = posStrnds(iseven(posStrnds(:, 2)), 1) + this.enzymeDNAFootprints3Prime(this.enzymeIndexs_rnaPolymerase);

                rnaPols.states(freePols(1:nNsbPols)) = rnaPols.nonSpecificallyBoundValue;
                rnaPols.positionStrands(freePols(1:nNsbPols), :) = posStrnds;
                trnspts.boundTranscriptionUnits(freePols(1:nNsbPols)) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptProgress(freePols(1:nNsbPols)) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptChromosome(freePols(1:nNsbPols)) = trnspts.nullTranscriptValue;

                nFrePols = nFrePols - nNsbPols;
                nBndPols = nBndPols + nNsbPols;
            end

            %set newly created RNA polymerases to free state
            if nNewPols > 0
                newPols = find(rnaPols.states == rnaPols.notExistValue, nNewPols, 'first');
                rnaPols.states(newPols) = rnaPols.freeValue;
                rnaPols.positionStrands(newPols, :) = 0;
                trnspts.boundTranscriptionUnits(newPols) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptProgress(newPols) = trnspts.nullTranscriptValue;
                trnspts.boundTranscriptChromosome(newPols) = trnspts.nullTranscriptValue;
            end

            %% Transcribe active sequences
            %compute progress of transcription and cost
            usedNTPs = zeros(4, 1);
            if ~isempty(actPols) && nFreTfs(this.enzymeIndexs_elongationFactor) && this.rnaPolymeraseElongationRate > 0
                posStrnds = rnaPols.positionStrands(actPols, :);
                iTUs = trnspts.boundTranscriptionUnits(actPols);

                %temporarily release RNA polymerase from old position strands
                releasedProteins = ...
                    this.releaseProteinFromSites(posStrnds, false, this.enzymeIndexs_rnaPolymerase, true, true) + ...
                    this.releaseProteinFromSites(posStrnds, false, this.enzymeIndexs_rnaPolymeraseHoloenzyme, true, true);
                if ~isequal(releasedProteins(this.enzymeIndexs_rnaPolymerase) + releasedProteins(this.enzymeIndexs_rnaPolymeraseHoloenzyme), numel(actPols))
                    throw(MException('Transcription:error', 'Unable to unbind protein'));
                end

                %extract active sequences (part of sequence that should be transcribed)
                actSeqs = cell(size(actPols)); %active sequences (those can be transcribed)
                for idx = 1:numel(actPols)
                    i = actPols(idx);
                    left = rnaPols.states(i);
                    actSeqs{idx} = tuSeqs{trnspts.boundTranscriptionUnits(i)}(left:end);
                end
                actSeqs = char(actSeqs);
                actSeqs = actSeqs(:, 1:min(end, this.rnaPolymeraseElongationRate));

                %elongation limits by RNA polymerase kinetic rate and transcription unit extent
                elngMax = min(this.rnaPolymeraseElongationRate, ...
                    tuLens(trnspts.boundTranscriptionUnits(actPols)) - rnaPols.states(actPols) + 1);

                %elongation limits by RNA-polymerase / RNA-polymerase interactions
                for strnd = 1 : 4
                    idxs = find(posStrnds(:, 2) == strnd);
                    if numel(idxs) <= 1
                        continue;
                    end
                    [~, order] = sort(posStrnds(idxs, 1));

                    if isodd(strnd)
                        elngMax(idxs(order)) = min(...
                            elngMax(idxs(order)), ...
                            diff([posStrnds(idxs(order), 1); posStrnds(idxs(order(1)), 1) + c.sequenceLen]) - ...
                            this.enzymeDNAFootprints(this.enzymeIndexs_rnaPolymerase));
                    else
                        elngMax(idxs(order)) = min(...
                            elngMax(idxs(order)), ...
                            diff([posStrnds(idxs(order(end)), 1) - c.sequenceLen; posStrnds(idxs(order), 1)]) - ...
                            this.enzymeDNAFootprints(this.enzymeIndexs_rnaPolymerase));
                    end
                end

                %elongation limits by RNA-polymerase DNA interactions (Note: independent of other RNA polymerases)
                [~, ~, ~, tmp_elngMax] = c.isRegionAccessible(...
                    posStrnds, (elngMax + 1) .* tuDirs(iTUs), ...
                    [], this.enzymeGlobalIndexs(this.enzymeIndexs_rnaPolymerase), true, [], true, false);
                if ~all(tmp_elngMax)
                    throw(MException('Transcription:error', 'Unable to bind proteins'));
                end
                elngMax = min(abs(tmp_elngMax) - 1, elngMax);

                %elongation limits by NTPs
                for i = 1:size(actSeqs, 1)
                    actSeqs(i, elngMax(i) + 1:end) = ' ';
                end
                [elngProg, this.substrates(this.substrateIndexs_ntp), usedNTPs] = ...
                    edu.stanford.covert.cell.sim.util.polymerize(...
                    actSeqs, this.substrates(this.substrateIndexs_ntp), 'ACGU', ' ', 0, 0, this.randStream);

                %release sigma factor after initiation
                nInits = sum(...
                    rnaPols.states(actPols) == rnaPols.activelyTranscribingValue & ...
                    elngProg > 0);
                nBndHols = nBndHols - nInits;
                nFreTfs(this.enzymeIndexs_sigmaFactor) = nFreTfs(this.enzymeIndexs_sigmaFactor) + nInits;
                nBndPols = nBndPols + nInits;

                %rebind RNA polymerase to new positions
                rnaPols.states(actPols) = rnaPols.states(actPols) + elngProg;
                rnaPols.positionStrands(actPols, 1) = posStrnds(:, 1) + elngProg .* tuDirs(iTUs);
                trnspts.boundTranscriptProgress(actPols) = rnaPols.states(actPols);

                initPols = rnaPols.states(actPols) == rnaPols.activelyTranscribingValue;
                if ...
                        ~all(this.bindProteinToChromosome(rnaPols.positionStrands(actPols(~initPols), :), this.enzymeIndexs_rnaPolymerase)) || ...
                        ~all(this.bindProteinToChromosome(rnaPols.positionStrands(actPols( initPols), :), this.enzymeIndexs_rnaPolymeraseHoloenzyme))
                    throw(MException('Transcription:error', 'Unable to bind proteins'));
                end
            end

            %release polymerases which have completed transcription, and increment
            %expression of the corresponding transcription unit
            if ~isempty(actPols) && nFreTfs(this.enzymeIndexs_elongationFactor) && nFreTfs(this.enzymeIndexs_terminationFactor)
                trmPols = actPols(rnaPols.states(actPols) > ...
                    tuLens(trnspts.boundTranscriptionUnits(actPols)));
                nTrmPols = numel(trmPols);
                if nTrmPols > this.substrates(this.substrateIndexs_water)
                   order = this.randStream.randperm(nTrmPols);
                   nTrmPols = this.substrates(this.substrateIndexs_water);
                   trmPols = trmPols(order(1:nTrmPols));
                end

                if nTrmPols > 0
                    if isscalar(trmPols)
                        this.RNAs(trnspts.boundTranscriptionUnits(trmPols)) = ...
                            this.RNAs(trnspts.boundTranscriptionUnits(trmPols)) + 1;
                    else
                        this.RNAs = this.RNAs + ...
                            histc(trnspts.boundTranscriptionUnits(trmPols, 1), 1:numel(tuLens));
                    end

                    %release proteins from chromsome
                    releasedProteins = this.releaseProteinFromSites(rnaPols.positionStrands(trmPols, :), false, this.enzymeIndexs_rnaPolymerase, true, true);
                    if ~isequal(releasedProteins(this.enzymeIndexs_rnaPolymerase), numel(trmPols))
                        throw(MException('Transcription:error', 'Unable to unbind protein'));
                    end

                    nFrePols = nFrePols + nTrmPols;
                    nBndPols = nBndPols - nTrmPols;

                    this.substrates(this.substrateIndexs_water)    = this.substrates(this.substrateIndexs_water)    - nTrmPols;
                    this.substrates(this.substrateIndexs_hydrogen) = this.substrates(this.substrateIndexs_hydrogen) + nTrmPols;

                    rnaPols.states(trmPols) = rnaPols.freeValue;
                    rnaPols.positionStrands(trmPols, :) = 0;
                    trnspts.boundTranscriptionUnits(trmPols) = trnspts.nullTranscriptValue;
                    trnspts.boundTranscriptProgress(trmPols) = trnspts.nullTranscriptValue;
                    trnspts.boundTranscriptChromosome(trmPols) = trnspts.nullTranscriptValue;
                end
            end

            %% diphosphate released by polymerization
            this.substrates(this.substrateIndexs_diphosphate) = ...
                this.substrates(this.substrateIndexs_diphosphate) + ...
                sum(usedNTPs);

            %% store enzyme state
            this.enzymes(this.enzymeIndexs_transcriptionFactors) = nFreTfs;
            this.boundEnzymes(this.enzymeIndexs_transcriptionFactors) = nBndTFs;
            this.enzymes(this.enzymeIndexs_rnaPolymerase) = nFrePols;
            this.boundEnzymes(this.enzymeIndexs_rnaPolymerase) = nBndPols;
            this.enzymes(this.enzymeIndexs_rnaPolymeraseHoloenzyme) = nFreHols;
            this.boundEnzymes(this.enzymeIndexs_rnaPolymeraseHoloenzyme) = nBndHols;
        end
    end

    %model helper functions
    methods
        function pBinds = computeRNAPolymeraseTUBindingProbabilities(this, protectDnaABoxes)
            c = this.chromosome;
            r = this.rnaPolymerases;

            %relative probability RNA polymerase binds each transcription
            %unit
            pBinds = this.transcriptionUnitBindingProbabilities(:, ones(1, 2)) .* ...
                r.transcriptionFactorBindingProbFoldChange .* ...
                r.supercoilingBindingProbFoldChange;

            %if functional DnaA boxes R1-5 are occupied by DnaA-ATP,
            %don't allow RNA polymerase to bind since it would knockoff
            %DnaA-ATP
            if nargin == 1 || protectDnaABoxes
                [posStrnds, complexs] = find(...
                    c.complexBoundSites(c.transcriptionUnitStartCoordinates(this.transcriptionUnitIndexs_DnaAR12345Boxes) + ...
                    (1:c.transcriptionUnitLengths(this.transcriptionUnitIndexs_DnaAR12345Boxes))' - 1, :));
                if any(ismembc(complexs, this.complexIndexs_DnaA_ATP) & posStrnds(:, 2) <= 2)
                    pBinds(this.transcriptionUnitIndexs_DnaAR12345Boxes, 1) = 0;
                end
                if any(ismembc(complexs, this.complexIndexs_DnaA_ATP) & posStrnds(:, 2) > 2)
                    pBinds(this.transcriptionUnitIndexs_DnaAR12345Boxes, 2) = 0;
                end
            end
        end

        function pStProbs = calcStateTransitionProbabilities(this, nPols, ntpProdRate, tuBindingProbs, method)
            import edu.stanford.covert.util.ComputationUtil;
            import edu.stanford.covert.util.ConstantUtil;

            rna = this.rna;
            rnaPols = this.rnaPolymerases;

            pStProbs = zeros(4, 4);
            pStProbs(rnaPols.nonSpecificallyBoundIndex, rnaPols.activelyTranscribingIndex) = ...
                ntpProdRate / (nPols * rnaPols.stateExpectations(rnaPols.activelyTranscribingIndex)) / ...
                (tuBindingProbs' * rna.lengths(rna.nascentIndexs));
            if nargin < 5 || isequal(method, 'handTuned')
                pStProbs(rnaPols.activelyTranscribingIndex, rnaPols.specificallyBoundIndex) = 1;
                pStProbs(rnaPols.nonSpecificallyBoundIndex, rnaPols.specificallyBoundIndex) = 0;
                pStProbs(rnaPols.specificallyBoundIndex, rnaPols.nonSpecificallyBoundIndex) = 0.10;
            else
                pStProbs(rnaPols.activelyTranscribingIndex, rnaPols.specificallyBoundIndex) = ...
                    pStProbs(rnaPols.nonSpecificallyBoundIndex, rnaPols.activelyTranscribingIndex) * ...
                    rnaPols.stateExpectations(rnaPols.activelyTranscribingIndex) / ...
                    rnaPols.stateExpectations(rnaPols.specificallyBoundIndex);
                pStProbs(rnaPols.specificallyBoundIndex, rnaPols.nonSpecificallyBoundIndex) = ...
                    pStProbs(rnaPols.nonSpecificallyBoundIndex, rnaPols.activelyTranscribingIndex) * ...
                    rnaPols.stateExpectations(rnaPols.activelyTranscribingIndex) / ...
                    rnaPols.stateExpectations(rnaPols.nonSpecificallyBoundIndex);
            end

            pStProbs(rnaPols.activelyTranscribingIndex, rnaPols.activelyTranscribingIndex) = ...
                1 - pStProbs(rnaPols.nonSpecificallyBoundIndex, rnaPols.activelyTranscribingIndex);
            pStProbs(rnaPols.specificallyBoundIndex, rnaPols.specificallyBoundIndex) = ...
                + 1 ...
                - pStProbs(rnaPols.activelyTranscribingIndex, rnaPols.specificallyBoundIndex)...
                - pStProbs(rnaPols.nonSpecificallyBoundIndex, rnaPols.specificallyBoundIndex);
            pStProbs(rnaPols.nonSpecificallyBoundIndex, rnaPols.nonSpecificallyBoundIndex) = ...
                1 - pStProbs(rnaPols.specificallyBoundIndex, rnaPols.nonSpecificallyBoundIndex);
            pStProbs(:, rnaPols.freeIndex) = pStProbs(:, rnaPols.nonSpecificallyBoundIndex);
        end
    end

    %get methods of dependent local state
    methods
        function value = getDryWeight(this)
            if size(this.RNAs, 3) == 1
                value = this.getDryWeight@edu.stanford.covert.cell.sim.Process() + ...
                    this.rna.molecularWeights(this.rna.nascentIndexs)' * this.RNAs ...
                    / edu.stanford.covert.util.ConstantUtil.nAvogadro;
            else
                value = this.getDryWeight@edu.stanford.covert.cell.sim.Process() + ...
                    permute(this.rna.molecularWeights(this.rna.nascentIndexs)' * permute(this.RNAs,[1 3 2]),[1 3 2]) ...
                    / edu.stanford.covert.util.ConstantUtil.nAvogadro;
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
% Derive_acrobot.m
%
% This script derives the equations of motion for the simple acrobot robot:
% a double pendulum with a motor at the elbow joint.
%

syms q1 q2 dq1 dq2 ddq1 ddq2 'real'   % states
syms u 'real' % actuation
syms m1 m2 g l1 l2 'real' % physical parameters

%%%% Unit vectors:
i = sym([1;0]);
j = sym([0;1]);

e1 = cos(q1)*(-j) + sin(q1)*(i);    % shoulder -> elbow
e2 = cos(q2)*(-j) + sin(q2)*(i);    % elbow -> wrist

%%%% State vectors:
z = [q1;q2;dq1;dq2];
dz = [dq1;dq2;ddq1;ddq2];

%%%% Kinematics:
p1 = l1*e1;
p2 = p1 + l2*e2;

dp1 = jacobian(p1,z)*dz;  %Chain rule to get velocity of elbow joint
dp2 = jacobian(p2,z)*dz;

ddp1 = jacobian(dp1,z)*dz;
ddp2 = jacobian(dp2,z)*dz;

%%%% Define a function for doing '2d' cross product: dot(a x b, k)
cross2d = @(a,b)(a(1)*b(2) - a(2)*b(1));

%%%% Angular momentum balance of system about shoulder joint:
sumTorques1 = cross2d(p1,-m1*g*j) + cross2d(p2,-m2*g*j);
sumInertial1 = cross2d(p1,m1*ddp1) + cross2d(p2,m2*ddp2);
eqn1 = sumTorques1-sumInertial1;

%%%% Angular momentum balance of outer link about elbow joint:
sumTorques2 = cross2d(p2-p1,-m2*g*j) + u;
sumInertial2 = cross2d(p2-p1,m2*ddp2);
eqn2 = sumTorques2-sumInertial2;

%%%% Write out dynamics in matrix form:   MM*ddq = ff
ddq = [ddq1;ddq2];
eqns = [eqn1;eqn2];
[MM,ff] = equationsToMatrix(eqns,ddq);

%%%% Generate an optimized matlab function for dynamics:
matlabFunction(MM(1,1),MM(1,2),MM(2,1),MM(2,2),ff(1),ff(2),G,B,...
    'file','autoGen_acrobotDynamics.m',...
    'vars',{q1,q2,dq1,dq2,u,m1,m2,g,l1,l2},...
    'outputs',{'M11','M12','M21','M22','f1','f2'});

%%%% Compute the energy of the system:
U = m1*g*dot(p1,j) + m2*g*dot(p2,j);   %Potential Energy
T = 0.5*m1*dot(dp1,dp1) + 0.5*m2*dot(dp2,dp2);   %Kinetic Energy

%%%% Generate an optimized matlab function for energy:
matlabFunction(U,T,...
    'file','autoGen_acrobotEnergy.m',...
    'vars',{q1,q2,dq1,dq2,m1,m2,g,l1,l2},...
    'outputs',{'U','T'});

%%%% Generate a function for computing the kinematics:
matlabFunction(p1,p2,dp1,dp2,...
    'file','autoGen_acrobotKinematics.m',...
    'vars',{q1,q2,dq1,dq2,l1,l2},...
    'outputs',{'p1','p2','dp1','dp2'});

author: MatthewPeterKelly

repo: OptimTraj


text:
function [a,b,alpha,p,chiopt,Cab,Calphap]=...
    wtls_line(xin,yin,uxin,uyin,guck);
% function [a,b,alpha,p,chiopt,Cab,Calphap]=...
%   wtls_line(xin,yin,uxin,uyin,guck);
%
% weighted total least squares (wtls) fit of a straigth line
% to a set of points with uncertainties in both coordinates
%
% input:        xin     abscissa vector
%               yin     ordinate vector
%               uxin    (standard) uncertainties of xin, same size as xin
%               uyin    (standard) uncertainties of yin, same size as yin
%               guck    flag, if >0, chisquare is plotted over whole
%                       alpha-interval, optimum alpha is shown as red +
%
% output:       a, b    usual straight line parameters
%                       y=a*x+b
%               alpha,p more stable parametrisation
%                       y*cos(alpha)-x*sin(alpha)-p=0
%                       alpha: slope angle in radians
%                       p: distance of straight line from (0,0)
%                       conversion: a=tan(alpha),b=p/cos(alpha)
%               chiopt  minimum chisquare found
%               Cab     covariances, [var(a),var(b),cov(a,b)]
%               Calphap covariances, [var(alpha),var(p),cov(alpha,p)]
%
% algorithm:    M. Krystek & M. Anton
%               Physikalisch-Technische Bundesanstalt Braunschweig, Germany
%               Meas. Sci. Tech. 18 (2007), pp3438-3442
%
% tested for Matlab 6 and Matlab 7
% testdata: script file pearson_york_testdata.m
%
% 2007-03-08

tol=1e-6; %"tolerance" parameter of fnimbnd, see there

global ux uy x y
% inputs
if nargin<5,guck=0;end
% force column vectors
x=xin(:);
y=yin(:);
ux=uxin(:);
uy=uyin(:);
% "initial guess":
p0=polyfit(x,y,1);
alpha0=atan(p0(1));
% one-dimensional search, use p=p^
[alphaopt,chiopt,exitflag] = ...
    fminbnd(@chialpha,alpha0-pi/2,alpha0+pi/2,optimset('TolX',tol));
% get optimum p from alphaopt
alpha=alphaopt;
uk2=ux.^2*sin(alpha)^2+uy.^2*cos(alpha)^2;
u2=1./mean(1./uk2);
w=u2./uk2;
xbar=mean(w.*x);
ybar=mean(w.*y);
p=ybar*cos(alpha)-xbar*sin(alpha);
% convert to a, b parameters y=a*x+b
a=sin(alpha)/cos(alpha);
b=p/cos(alpha);
% --- uncertainty calculation, covariance matrix = 2*inv(Hessian(chi2)) ---
n=length(x);
vk=y.*cos(alpha)-x*sin(alpha)-p;
vka=-y*sin(alpha)-x*cos(alpha);
vkaa=-vk-p;
fk=vk.*vk;
fka=2*vk.*vka;
fkaa=2*(vka.^2+vk.*vkaa);
gk=uk2;
gka=2*sin(alpha)*cos(alpha)*(ux.^2-uy.^2);
gkaa=2*(ux.^2-uy.^2)*(cos(alpha)^2-sin(alpha)^2);
Hpp=2*n/u2;
Halphap=-2*sum((vka.*gk-gka.*vk)./gk.^2);
Halphaalpha=...
    sum(fkaa./gk-2*fka.*gka./gk.^2+2*gka.^2.*fk./gk.^3-gkaa.*fk./gk.^2);
NN=2/(Hpp*Halphaalpha-Halphap^2);
var_p=NN*Halphaalpha;
var_alpha=NN*Hpp;
cov_alphap=-NN*Halphap;
Calphap=[var_alpha,var_p,cov_alphap];
% ------ convert to a & b covariance matrix, following DIN 1319 (4)
var_a=var_alpha/cos(alpha)^4;
var_b=(var_alpha*p*p*sin(alpha)^2+var_p*cos(alpha)^2+...
    2*cov_alphap*p*sin(alpha)*cos(alpha))/cos(alpha)^4;
cov_ab=(var_alpha*p*sin(alpha)+cov_alphap*cos(alpha))/cos(alpha)^4;
Cab=[var_a,var_b,cov_ab];
% ------ end of uncertainty calculation -----------------------------------
%
%-------------------- plotting section ------------------------------------
if guck~=0
    alphatest=linspace(alpha-pi/2,alpha+pi/2,1000);
    for k=1:length(alphatest),chitest(k)=chialpha(alphatest(k));end
    plot(alphatest,chitest,alphaopt,chiopt,'+r')
    xlabel('alpha')
    ylabel('\chi^2')
    grid on
    disp([mfilename,' :: paused, hit any key to continue'])
    pause
end
% ------------------ end of plotting section ------------------------------
return
% ------------------ subfunction ------------------------------------------
function chi=chialpha(alpha)
    global ux uy x y
    uk2=ux.^2*sin(alpha)^2+uy.^2*cos(alpha)^2;
    u2=1./mean(1./uk2);
    w=u2./uk2;
    xbar=mean(w.*x);
    ybar=mean(w.*y);
    p=ybar*cos(alpha)-xbar*sin(alpha);
    chi=sum((y*cos(alpha)-x*sin(alpha)-p).^2./uk2);
return
% ----------- end of subfunction ------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function obj = addpoints(obj, xyz, varargin)
% Plots points on fmridisplay objects (e.g., montages of slices)
%
% :Usage:
% ::
%
%     newax = addpoints(obj, xyz, varargin)
%
% Registers handles with the object (referred to as obj)
%
% - enter xyz as n x 3 list of coordinates in mm to plot (world space)
% - Points or text labels or both
% - Flexible slice spacing, colors, marker sizes/styles, axis layout (one row/standard square)
% - axial, saggital, or coronal orientation handled automatically
% - Multiple different sets of points can be plotted in different colors/text labels
%
% :Optional Inputs:
%
% Takes all inputs of plot_points_on_slice.  See help for additional
% documentation of options.
%
%   **{'text', 'textcodes'}:**
%        cell array of text values corresponding to points
%
%   **{'condf' 'colorcond'}:**
%        vector of integers to define color conditions
%
%   **'close_enough':**
%        mm within which to plot; defined automatically based on slice distance if not entered
%
%   **'color':**
%        string, 'b', or vector, [1 0 0], to define colors; cell if condf is used, e.g., {'b' 'g'}
%
%   **{'marker', 'MarkerStyle'}:**
%        e.g., 'o', 'v', 's'
%
%   **{'MarkerSize', 'markersize'}:**
%
%   **{'MarkerFaceColor', 'markerfacecolor'}:**
%        see color above
%
% :Examples:
%
% Plot points (i.e., coordinate locations) for xyz coords:
% ::
%
%    o2 = addpoints(o2, DB.xyz, 'MarkerFaceColor', 'b', 'Marker', 'o', 'MarkerSize', 4);
%    o2 = addpoints(o2, DB.xyz, 'text', DB.textcodes, 'condf', DB.condf, 'color', {'b' 'g'});
%    o2 = removepoints(o2);

wh_montage = 1:length(obj.montage); % select which montages; default = all

whm = strcmp(varargin, 'wh_montages') | strcmp(varargin, 'wh_montage') | strcmp(varargin, 'which_montages') | strcmp(varargin, 'which montages');
if any(whm)
    whm = find(whm);
    wh_montage = varargin{whm(1) + 1};
end

% Montages
% -------------------------------------------------------------------------

for i = wh_montage

    if ~isfield(obj.montage{i}, 'plotted_point_handles') || isempty(obj.montage{i}.plotted_point_handles)

        obj.montage{i}.plotted_point_handles = [];

    end

    % Plot it
    pointhan = plotpoints(xyz, obj.montage{i}, varargin{:});

    obj.montage{i}.plotted_point_handles = [obj.montage{i}.plotted_point_handles pointhan];

end

end  % main function

function pointhan = plotpoints(xyz, montagestruct, varargin)

% fixed fields

myview = montagestruct.orientation;
slicemm = montagestruct.slice_mm_coords;
axhan = montagestruct.axis_handles;

% optional inputs

textcodes = [];
texthandles = [];
condf = [];

close_enough = min(diff(slicemm)) ./ 2;
if isempty(close_enough), close_enough = 8; end

color = 'k';
marker = 'o';
markersize = 12;
markerfacecolor = 'k';

% ------------------------------------------------------
% parse inputs
% ------------------------------------------------------

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            % reserved keywords

            case {'text', 'textcodes'} % do not pass on to slice plot...
                textcodes = varargin{i + 1};
                varargin{i+1} = [];
                varargin{i} = [];

            case {'condf' 'colorcond'}, condf = varargin{i + 1};

            case 'close_enough', close_enough = varargin{i + 1};

            case {'Color','color'} % do not pass on...
                color = varargin{i+1};
                varargin{i+1} = [];
                varargin{i} = [];

            case {'marker', 'Marker', 'MarkerStyle'}, marker = varargin{i+1}; varargin{i+1} = [];

            case {'MarkerSize', 'markersize'}, markersize = varargin{i+1}; varargin{i+1} = [];

            case {'MarkerFaceColor', 'markerfacecolor'}, markerfacecolor = varargin{i+1}; varargin{i+1} = [];

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% SETUP
% -----------------------------------------------

fprintf('Plotting points within %3.2f mm of slices\n', close_enough);

if isempty(condf)
    condf = ones(size(xyz, 1), 1);

    if iscell(color)
        error('Color should not be a cell array.');
    end
    color = {color};

else
    if ~iscell(color)
        error('When entering condf, enter colors cell with same number of entries.');
    end
end

u = unique(condf);
n = length(u);

% Do the work for each slice
% -----------------------------------------------

pointhan = [];

for i = 1:length(axhan)
    axes(axhan(i));

    for j = 1:n % For each color code

        % Select coordinates (2D) for plot
        slicexyz = xyz;

        switch myview
            case 'axial'
                whcol = 3;

            case {'sagg', 'sagittal', 'saggital'}

                whcol = 1;

            case {'cor', 'coronal'}

                whcol = 2;

            otherwise
                error('Unknown slice orientation.')
        end

        wh_to_plot = condf == u(j) & abs(slicexyz(:, whcol) - slicemm(i)) <= close_enough;

        slicexyz = slicexyz(wh_to_plot, :);

        slicexyz(:, whcol) = [];

        if ~isempty(slicexyz)

            if isempty(textcodes)

                pointhan(end + 1) = plot(slicexyz(:, 1), slicexyz(:, 2), '.', 'color', color{j}, 'Marker', marker, 'MarkerSize', markersize, 'MarkerFaceColor', markerfacecolor);

            else
                my_text = textcodes(wh_to_plot);

                pointhan = [pointhan plottext(slicexyz(:, 1), slicexyz(:, 2), my_text, color{j}, myview)];

            end

        end

    end % condition color codes

end % slices

end % function

function texthandles = plottext(x, y, textcodes, color, orientation)

% adjust for font placement
switch orientation
    case 'sagittal'
        xshift = 6;
    case 'axial'
        xshift = -6;

end

texthandles = [];
for j = 1:length(textcodes)
    % text labels

    if ischar(color)

        texthandles(j) = text(x(j) + xshift, y(j), textcodes{j},'Color',color,'FontSize',12,'FontWeight','bold');

    else
        texthandles(j) = text(x(j)  + xshift, y(j), textcodes{j},'Color',color,'FontSize',12,'FontWeight','bold');

    end

end

end % plottext

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [C,phi,S12,S1,S2,f,zerosp,confC,phistd,Cerr]=coherencysegcpb(data1,data2,win,params,segave,fscorr)
% Multi-taper coherency,cross-spectrum and individual spectra with segmenting
%   computed by segmenting two univariate time series into chunks - continuous and binned point process
%
% Usage:
% [C,phi,S12,S1,S2,f,zerosp,confC,phistd,Cerr]=coherencysegcpb(data1,data2,win,params,segave,fscorr)
% Input:
% Note units have to be consistent. See chronux.m for more information.
%       data1 (column vector, continuous data) -- required
%       data2 (column vector, binned point process data) -- required
%       win   (length of segments) - required
%       params: structure with fields tapers, pad, Fs, fpass, err
%       - optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
%                                   [0 p] or 0 - no error bars) - optional. Default 0.
%       segave (average over segments for 1, don't average for 0)
%       fscorr   (finite size corrections, 0 (don't use finite size corrections) or
%                1 (use finite size corrections) - optional
%                (available only for spikes). Defaults 0.
% Output:
%       C (magnitude of coherency - frequencies x segments if segave=0; dimension frequencies if segave=1)
%       phi (phase of coherency - frequencies x segments if segave=0; dimension frequencies if segave=1)
%       S12 (cross spectrum -  frequencies x segments if segave=0; dimension frequencies if segave=1)
%       S1 (spectrum 1 - frequencies x segments if segave=0; dimension frequencies if segave=1)
%       S2 (spectrum 2 - frequencies x segments if segave=0; dimension frequencies if segave=1)
%       f (frequencies)
%       zerosp (1 for trials where no spikes were found, 0 otherwise)
%       confC (confidence level for C at 1-p %) - only for err(1)>=1
%       phistd - theoretical/jackknife (depending on err(1)=1/err(1)=2) standard deviation for phi
%                Note that phi + 2 phistd and phi - 2 phistd will give 95% confidence
%                bands for phi - only for err(1)>=1
%       Cerr  (Jackknife error bars for C - use only for Jackknife - err(1)=2)

if nargin < 3; error('Need data1 and data2 and size of segment'); end;
if nargin < 4; params=[]; end;
[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear tapers pad fpass trialave
if nargin < 5 || isempty(segave); segave=1;end;
if nargin < 6 || isempty(fscorr); fscorr=0; end;
if nargout > 7 && err(1)==0;
%   Errors computed only if err(1) is non-zero. Need to change params and run again.
    error('When errors are desired, err(1) has to be non-zero.');
end;
if nargout > 9 && err(1)~=2;
    error('Cerr computed only for Jackknife. Correct inputs and run again');
end;

N=check_consistency(data1,data2);
dt=1/Fs; % sampling interval
T=N*dt; % length of data in seconds
E=0:win:T-win; % fictitious event triggers
win=[0 win]; % use window length to define left and right limits of windows around triggers
data1=createdatamatc(data1,E,Fs,win); % segmented data 1
data2=createdatamatpb(data2,E,Fs,win); % segmented data 2
params.trialave=segave;
if nargout==7;
  [C,phi,S12,S1,S2,f,zerosp]=coherencycpb(data1,data2,params,fscorr); % compute coherency for segmented data
elseif nargout==9;
  [C,phi,S12,S1,S2,f,zerosp,confC,phistd]=coherencycpb(data1,data2,params,fscorr); % compute coherency for segmented data
elseif nargout==10;
  [C,phi,S12,S1,S2,f,zerosp,confC,phistd,Cerr]=coherencycpb(data1,data2,params,fscorr); % compute coherency for segmented data
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function test_external_images

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY hsv2rgb

filelist = {'hsv2rgb'};

% NOTE: hsv2rgb was part of standard matlab in older matlab versions

[ftver, ftpath] = ft_version;
restoredefaultpath

% ensure that the the 'compat' (i.e. external/signal) is used
global ft_default
ft_default.toolbox.signal = 'compat';
addpath(ftpath);
ft_defaults;

for k = 1:numel(filelist)
  assert(exist(filelist{k}, 'file')==2);

  fprintf('testing the functionality of %s\n', filelist{k});
  switch filelist{k}
    case 'hsv2rgb'
      assert(isequal(hsv2rgb([0 1 1]),[1 0 0]));
      assert(isequal(hsv2rgb([1 0 1]),[1 1 1]));
    otherwise
      ft_error('function %s is not part of the official external/images directory', filelist{k});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compare the external/images output with the matlab version
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tmp = rand(10,10,3);
for k = 1:numel(filelist)
  try
    % this only works for the functions that create a window
    funhandle = str2func(filelist{k});
    output1.(filelist{k}) = funhandle(tmp);
  end
end

restoredefaultpath

% ensure that the the matlab version is used
ft_default.toolbox.imagesc = 'matlab';
addpath(ftpath);
ft_defaults;

for k = 1:numel(filelist)
  try
    funhandle = str2func(filelist{k});
    output2.(filelist{k}) = funhandle(tmp);
  end
end

fn1 = fieldnames(output1);
fn2 = fieldnames(output2);
assert(isequal(sort(fn1), sort(fn2)));
for k = 1:numel(fn1)
  assert(isalmostequal(output1.(fn1{k}), output2.(fn1{k}), 'abstol', 10*eps));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [mu,varOut,numVals] = knuthMeanVar(val,newCalcFlag)
% function [mu,varOut] = knuthMeanVar(val,newCalcFlag)
%
% Knuth method for mean and standard deviation (like ITK)
%
% Input: value to add to the existing mean and variance calculation.
%
% Output: Mean and variance (var(xV,1))
%
% Example:
% newCalcFlag = true;
% knuthMeanVar([],true);
% valV = [1:10,100:105];
% for i = 1:length(valV)
%     [mu,varOut,numVals] = knuthMeanVar(valV(i));
% end
%
% To start a new calculation:
% newCalcFlag = true;
% knuthMeanVar([],newCalcFlag);
% valV = [1:5];
% for i = 1:length(valV)
%     [mu,varOut,numVals] = knuthMeanVar(valV(i));
% end
%
% APA, 6/24/2019

if exist('newCalcFlag','var') && newCalcFlag
    clear persistent n
    clear persistent muPrev
    clear persistent varN
end

persistent n;
persistent muPrev;
persistent varN;

if isempty(val)
    numVals = n;
    mu = muPrev;
    varOut = varN/n;
    return;
end
if isempty(muPrev)
    muPrev = 0;
    varN = 0;
    n = 0;
end
n = n + 1;
mu = muPrev + (val - muPrev)/n;
varN = varN + (val-mu)*(val-muPrev);
muPrev = mu;
varOut = varN/n;
numVals = n;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR