text:
%Translation
%
% @wholeCellModelID Process_Translation
% @name             Translation
% @description
%   Biology
%   ===============
%   Translation is the first step in the synthesis of function proteins whereby
%   the ribosome, accessory enzymes, and tRNAs transcode mRNAs and produce
%   amino acid polymers. Following translation proteins must be matured to
%   become fully functional:
%    1. Protein Processing I: N-terminal methionine deformylation and cleavage
%    2. Protein Translocation: Translocation of integral membrane proteins,
%       lipoproteins, and secreted proteins into and through the cell membrane
%    3. Protein Processing II: Diacylglyceryl addition to lipoproteins. Signal
%       sequence cleavage of lipo- and secreted proteins.
%    4. Protein Folding: Folding and prosthetic group coordination of proteins.
%    5. Protein Modification: Addition of covalently attached chemical groups to
%       proteins.
%    6. Macromolecular Complexation: Incorporation of protein monomers into
%       larger macromolecular assemblies.
%    7. Ribosome Assembly: Special case of macromolecular complexation of the
%       ribosome which requires GTPases.
%    8. Terminal Organelle Assembly: Special case of translocation which
%       sequesters several protein species in a membrane-bound bleb referred to
%       as the terminal organelle.
%    9. Protein Activation: Functional activation and inactivation of proteins
%       in response to their surrounding chemical environment.
%
%   Translation begins with the recruitment of the 30S and 50S ribosomal
%   particles and initiation factor 3 (IF3) to an mRNA molecule. Next the
%   ribosomal particles assembly into a 70S ribosome on the mRNA molecule with
%   the help of initiation factors. Third the ribosome polymerizes amino
%   acids with the help of elongation factors. Finally release factor
%   (MG_258) recognizes the stop codon UAG, release factor hydolyzes the
%   peptidyl tRNA bond and dissociates, ribosome recycling factor
%   dissociates the E-site tRna, elongation factor G release the release
%   factor and 50S ribosome, and initiation factor 3 dissociates the 30S
%   ribosome, P-site tRNA, and mRNA.
%
%   When a ribosome stalls its lastest tRNA is expelled and replaced by the
%   tRNA-like domain of a tmRNA molecule. Next the mRNA-like domain of the tmRNA
%   expels the bound mRNA. Third the ribosome resumes polymerization, now
%   using the tmRNA's mRNA-like domain as its template. This results in the
%   production of an amino acid polymer containing a C-terminal proteolysis tag.
%   Finally, the proteolysis tag will be recognized by the protein degradation
%   machinery, and the amino acid polymer will be degraded into its individual
%   component amino acids.
%
%   This process simulates protein translation by ribosomes and acessory
%   initiation, elongation, and termination factors. The process also simulates
%   the identification of stalled ribosomes by the tmRNA, the replacement of the
%   tRNA and mRNA with the tmRNA, and the synthesis of the proteolysis tag
%   encoded by the tmRNA's mRNA-like domain.
%
%   Knowledge Base
%   ===============
%   The knowledge base contains the region of the genome each mRNA species
%   corresponds to, and the location of the tmRNA mRNA-like domain on the
%   genome. This information is converted into a sequence of tRNAs which are
%   required to polymerize each protein and proteolysis tag by the knowledge
%   base class' computeTRNASequences method.
%
%   Representation
%   ===============
%   The properties substrates, enzymes, monomers, mRNAs, freeTRNAs,
%   aminoacylatedTRNAs, and aminoacylatedTMRNA represents the counts of
%   metabolites, translation enzymes, nascent peptides, mature mRNAs, free
%   tRNAs, aminoacylated tRNAs, and aminoacylated tmRNAs.
%
%   monomerTRNASequences and tmRNAProteolysisTagTRNASequence represent the
%   sequences of tRNAs required to polymerize each protein and proteolysis tag.
%   These tRNA sequences are computed from the DNA sequence of the mRNA gene /
%   tmRNA gene by the knowledge  base classs' computeTRNASequences method.
%
%   ribosomeStates represents the state/pseudostate (actively translating, free,
%   non-existent) of each ribosome. boundMRNAs indicates the mRNA
%   species each actively translating ribosome is bound to. The non-existent
%   pseudostate is used to keep track of elements in ribosomeStates
%   and boundMRNAs which have been allocated.
%   nascentMonomerLengths and proteolysisTagLengths represeent
%   the position of actively  translating ribosomes on mRNAs, and if stalled on
%   a tmRNA.
%
%   Initialization
%   ===============
%   All protein monomers are initialized to the mature state, and in their
%   correct localization. This is achieved by the simulation class
%   initializeState method.
%
%   Ribosomes are initialized to their steady state:
%   - Each ribosome is randomly assigned (without replacement) to mRNA species,
%     weighted by the current expression of the mRNAs,.
%   - Each ribosome is randomly assigned to positions within the assigned mRNA
%     with uniform probably.
%   - No ribosomes are initialized to the stalled state, since the expected
%     occupancy is negligible. No tmRNAs are initialized to the bound state.
%
%   tRNAs and tmRNAs are initialized to the aminoacyated state by
%   simulation.initializeState.
%
%   Simulation
%   ===============
%   Evolves the states of ribosomes among two states:
%   - actively translating
%   - free
%
%   Transition to from the free state to the actively translating state is
%   allowed when
%   - ribosome binding factor A is present
%   - initiation factors (IF-1, IF-2, and IF-3) are present
%   - one unit of energy (GTP) is available.
%   At this time the ribosome binds an mRNA randomly with uniform probability.
%
%   At the next iteration the initiation factors are released when the first
%   amino acid, f-methionine binds and elongation begins, assuming that
%   elongation factors (EF-tu, TS, and G) and energy (GTP) is available.
%
%   ASSUMPTION made here is that one of each is sufficient for each ribosome,
%   but need a separate set for each ribosome. So each ribosome can only
%   translate if a full set exists for it. Also, for cases in which
%   translation of a peptide finishes partway through the time step, the EFs
%   are free, but haven't added the complexity of another ribosome being able
%   to take up EFs partway through the timestep.
%
%   Finally, once all amino acids of a protein have been translated,
%   termination occurs if
%   - at least one terminator (RF-1) available
%   - at least one recycling factor available
%   - at least one elongation factor G available
%   - at least one trigger factor available
%   - energy is available.
%   The ribosome will be available to bind mRNA at the following iteration.
%
%   Ribosomes are created in the free state.
%
%   Algorithm
%   +++++++++++++++
%   1. Up to limit of elongation factors, randomly select actively translating
%      ribosomes to elongate.
%   2. Up to limit of initiation factors and energy, randomly select ribosomal
%      particles to initiatate. These ribosome will be able to start elongating
%      at the next time step. Randomly select mRNA species for each initiating
%      mRNA to bind to, weighted by the counts of each mRNA species. Update
%      ribosomeStates, boundMRNAs. Update substrates.
%   3. Allocate available amino acids and energy among actively translating
%      ribosomes. Update ribosome nascentMonomerLengths, or
%      proteolysisTagLengths for stalled ribosomes with the number of
%      polymerized bases. Update substrates.
%   4. If ribosome has reached end of (t)mRNA and termination factors
%      available, increase count of protein and dissolve (t)mRNA-ribosome
%      complex. Update substrates.
%   5. If ribosome hasn't advanced, then with a small probability transition
%      ribosome to stalled state. Expel mRNA and replace with tmRNA. Update
%      ribosomeStates. Update substrates.
%
%   References
%   ===============
%   1. Petry S, Weixlbaumer A, Ramakrishnan V (2008). The termination of
%      translation. Curr Opin Struct Biol. 18(1):70-7. [PUB_0226].
%
% Author: Markus Covert, mcovert@stanford.edu
% Author: Jayodita Sanghvi, jayodita@stanfod.edu
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 1/4/2010

%TODO: Bind specific mRNAs rather than kinds of mRNAs.
%
%TODO: Prevent ribosomes from colliding/overlapping on an mRNA. Currently this
%      is only somewhat ensured by allowing only one ribosome to bind each mRNA
%      per time step.
%
%TODO: Associate elongation factors to mRNA for just 1 amino acid, not
%      elongationRate amino acids. With this change, all mRNAs will
%      translate at each second but with a slower rate. Related to the
%      assumption above. That is we should translate total of min(energy,
%      min(EFs)*rate*this.simulationTime) amino acids, where rate is a
%      rate of EF recycling per second, that is how many times an
%      individual EF can be used per second.
classdef Translation < edu.stanford.covert.cell.sim.Process
    %property annotations
    properties (Constant)
        optionNames__              = {}; %names of option properties
        fixedConstantNames__       = {   %names of fixed constant properties
            'ribosomeElongationRate'
            'tmRNABindingProbability'
		    };
        fittedConstantNames__      = {}; %names of fitted constant properties
        localStateNames__          = {   %names of simulation state properties redundant with timecourses in this or other processes or the simulation
            'mRNAs'
            'freeTRNAs'
            'aminoacylatedTRNAs'
            'aminoacylatedTMRNA'
            'monomers'
            };
    end

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

        substrateWholeCellModelIDs = { %whole cell model IDs of substrates
            'ALA'; 'ARG'; 'ASN'; 'ASP'; 'CYS'; 'GLN'; 'GLU'; 'GLY'; 'HIS'; 'ILE';
            'LEU'; 'LYS'; 'MET'; 'PHE'; 'PRO'; 'SER'; 'THR'; 'TRP'; 'TYR'; 'VAL'; 'FMET';
            'GTP'; 'GDP'; 'PI'; 'H2O'; 'H'};
        substrateIndexs_aminoAcids  = (1:20)'; %indices within substrates of amino acids
        substrateIndexs_fMet        = 21;      %index within substrates of formylmethionine
        substrateIndexs_gtp         = 22;      %index within substrates of GTP
        substrateIndexs_gdp         = 23;      %index within substrates of GDP
        substrateIndexs_phosphate   = 24;      %index within substrates of phosphate
        substrateIndexs_water       = 25;      %index within substrates of water
        substrateIndexs_hydrogen    = 26;      %index within substrates of hydrogen

        enzymeWholeCellModelIDs = { %whole cell model ids of enzymes
            'MG_173_MONOMER';   %translation initiation factor IF-1
            'MG_142_MONOMER';   %translation initiation factor IF-2
            'MG_196_MONOMER';   %translation initiation factor IF-3
            'MG_089_DIMER';     %translation elongation factor G
            'MG_026_MONOMER';   %translation elongation factor P
            'MG_451_DIMER';     %translation elongation factor Tu
            'MG_433_DIMER';     %translation elongation factor Ts
            'MG_258_MONOMER';   %peptide chain release factor 1
            'MG_435_MONOMER';   %ribosome recycling factor
            'RIBOSOME_30S';     %30S ribosomal subunit
            'RIBOSOME_30S_IF3'; %30S ribosomal subunit - translation initiation factor IF-3 complex
            'RIBOSOME_50S';     %50S ribosomal subunit
            'RIBOSOME_70S';     %70S ribosome
            'MG_0004';          %transfer-messenger RNA, tmRNA, MCS6, 10Sa RNA
            'MG_059_MONOMER';   %SsrA-binding protein
            'MG_083_MONOMER'};  %peptidyl-tRNA hydrolase
        enzymeIndexs_translationFactors          = (1:9)';    %indices within enzymes of translation factors
        enzymeIndexs_initiationFactors           = (1:3)';    %indices within enzymes of initiation factors
        enzymeIndexs_initiationFactor1           = 1;         %index within enzymes of initiation factor 1
        enzymeIndexs_initiationFactor2           = 2;         %index within enzymes of initiation factor 2
        enzymeIndexs_initiationFactor3           = 3;         %index within enzymes of initiation factor 3
        enzymeIndexs_elongationFactors           = (4:7)';    %indices within enzymes of elongation factors
        enzymeIndexs_elongationGFactor           = 4;         %index within enzymes of elongation factor G
        enzymeIndexs_terminationFactor           = 8;         %index within enzymes of termination factor
        enzymeIndexs_recyclingFactor             = 9;         %index within enzymes of recylcing factor
        enzymeIndexs_ribosome30S                 = 10;        %index within enzymes of 30S ribosomal particle
        enzymeIndexs_ribosome30SIF3              = 11;        %index within enzymes of 30S ribosomal particle-initiation factor 3 complex
        enzymeIndexs_ribosome50S                 = 12;        %index within enzymes of 50S ribosomal particle
        enzymeIndexs_ribosome70S                 = 13;        %index within enzymes of 70S ribosome
        enzymeIndexs_abortiveTranslationEnzymes  = (14:16)';  %indices within enzymes of abortive translation enzymes
        enzymeIndexs_tmRNA                       = 14;        %index within enzymes of tmRNA
        enzymeIndexs_tmRNABindingProtein         = 15;        %index within enzymes of tmRNA binding protein
        enzymeIndexs_peptidylTRNAHydrolase       = 16;        %index within enzymes of peptidyl-tRNA hydrolase

        freeTRNAWholeCellModelIDs          %IDs of free (non-aminoacylated) tRNAs
        aminoacylatedTRNAWholeCellModelIDs %IDs of aminoacylated tRNAs
        aminoacylatedTMRNAWholeCellModelID %ID of aminoacylated tmRNA

        tmRNAGlobalIndex %index of tmRNA within simulation.matureRNAIndexs
    end

    %fixed biological constants
    properties
        ribosomeElongationRate             %ribosomse elongation rate (amino acids per second per ribosome) [PUB_0564]
        tmRNABindingProbability            %probability that tmRNA is recruited to non-elongating, active ribosome
    end

    %global state (stored locally for convenience)
    properties
        monomers                       %numbers of protein monomers
        mRNAs                          %numbers of mRNAs
        freeTRNAs                      %numbers of free (non-aminoacylated) tRNAs
        freeTMRNA                      %numbers of free tmRNA
        aminoacylatedTRNAs             %numbers of aminoacylated tRNAs
        aminoacylatedTMRNA             %numbers of aminoacylated tmRNA
        boundTMRNA                     %numbers of bound tmRNA
    end

    %global state (referenced locally for convenience)
    properties
        polypeptide                    %New polypeptides state class
        ribosome                       %Ribosome state class
    end

    %constructor
    methods
        function this = Translation(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.polypeptide = simulation.state('Polypeptide');
            this.ribosome = simulation.state('Ribosome');
            this.states = [this.states; {this.polypeptide; this.ribosome}];
        end

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

            %tmRNA
            this.tmRNAGlobalIndex = this.enzymeRNAGlobalIndexs(this.enzymeRNALocalIndexs==this.enzymeIndexs_tmRNA);

            %RNA
            this.freeTRNAWholeCellModelIDs          = this.rna.wholeCellModelIDs(this.rna.matureIndexs(this.rna.matureTRNAIndexs));
            this.aminoacylatedTRNAWholeCellModelIDs = this.rna.wholeCellModelIDs(this.rna.aminoacylatedIndexs(this.rna.matureTRNAIndexs));
            this.aminoacylatedTMRNAWholeCellModelID = this.rna.wholeCellModelIDs(this.rna.aminoacylatedIndexs(this.tmRNAGlobalIndex));
        end

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

            r = this.rna;
            m = this.monomer;

            if size(r.counts, 3) == 1
                this.mRNAs = ...
                    r.matureRNAGeneComposition(this.gene.mRNAIndexs, r.matureMRNAIndexs) * ...
                    r.counts(r.matureIndexs(r.matureMRNAIndexs), this.compartment.cytosolIndexs);
            else
                this.mRNAs = multiprod(...
                    r.matureRNAGeneComposition(this.gene.mRNAIndexs, r.matureMRNAIndexs), ...
                    r.counts(r.matureIndexs(r.matureMRNAIndexs), this.compartment.cytosolIndexs, :), [1 2], 1);
            end
            this.monomers           = m.counts(m.nascentIndexs, this.compartment.cytosolIndexs, :);
            this.freeTRNAs          = r.counts(r.matureIndexs(this.rna.matureTRNAIndexs), this.compartment.cytosolIndexs, :);
            this.freeTMRNA          = r.counts(r.matureIndexs(this.tmRNAGlobalIndex), this.compartment.cytosolIndexs, :);
            this.aminoacylatedTRNAs = r.counts(r.aminoacylatedIndexs(this.rna.matureTRNAIndexs), this.compartment.cytosolIndexs, :);
            this.aminoacylatedTMRNA = r.counts(r.aminoacylatedIndexs(this.tmRNAGlobalIndex), this.compartment.cytosolIndexs, :);
            this.boundTMRNA         = r.counts(r.boundIndexs(this.tmRNAGlobalIndex), this.compartment.cytosolIndexs, :);
        end

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

            r = this.rna;
            m = this.monomer;

            m.counts(m.nascentIndexs, this.compartment.cytosolIndexs, :)  = this.monomers;
            r.counts(r.matureIndexs(this.rna.matureTRNAIndexs), this.compartment.cytosolIndexs, :) = this.freeTRNAs;
            r.counts(r.matureIndexs(this.tmRNAGlobalIndex), this.compartment.cytosolIndexs, :) = this.freeTMRNA;
            r.counts(r.aminoacylatedIndexs(this.rna.matureTRNAIndexs), this.compartment.cytosolIndexs, :) = this.aminoacylatedTRNAs;
            r.counts(r.aminoacylatedIndexs(this.tmRNAGlobalIndex), this.compartment.cytosolIndexs, :) = this.aminoacylatedTMRNA;
            r.counts(r.boundIndexs(this.tmRNAGlobalIndex), this.compartment.cytosolIndexs, :) = this.boundTMRNA;
        end
    end

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

            this.monomers           = zeros(length(this.polypeptide.monomerLengths), 1, numTimePoints);
            this.mRNAs              = zeros(length(this.gene.mRNAIndexs), 1, numTimePoints);
            this.freeTRNAs          = zeros(length(this.rna.matureTRNAIndexs), 1, numTimePoints);
            this.freeTMRNA          = zeros(1, 1, numTimePoints);
            this.aminoacylatedTRNAs = zeros(length(this.rna.matureTRNAIndexs), 1, numTimePoints);
            this.aminoacylatedTMRNA = zeros(1, 1, numTimePoints);
            this.boundTMRNA         = zeros(1, 1, numTimePoints);
        end
    end

    %model
    methods
        %Calculate
        %- contribution to FBA objective
        %- minimum expression consistent with cell cycle length
        %
        %Contributions, by process, to protein synthesis:
        %                            AAs --(translation)----> nascent monomer + water
        %          nascent monomer + H2O --(processing)-----> processed protein + AAs
        %processed monomer + metabolites --(folding)--------> folded monomer
        %   folded monomer + metabolites --(modification)---> mature monomer + metabolites
        %              modified monomers --(complexation)---> nascent complexs
        % nascent complexs + metabolites --(folding)--------> mature complexes
        %          monomer/complex + H2O --(decay)----------> AAs + modified AAs
        %        phosphorylated AA + H2O --(catabolism)-----> AA + Pi
        %            lipoyl lysine + H2O --(catabolism)-----> lysine + lipoate
        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

            %breaking tRNA-peptidyl bond
            bmProd(this.substrateIndexs_water) = ...
                bmProd(this.substrateIndexs_water) + ...
                this.polypeptide.monomerLengths' * states.monomerProductions;
            byProd(this.substrateIndexs_hydrogen) = ...
                byProd(this.substrateIndexs_hydrogen) + ...
                this.polypeptide.monomerLengths' * states.monomerProductions;

            %peptide bond formation
            byProd(this.substrateIndexs_water) = ...
                byProd(this.substrateIndexs_water) + ...
                max(0, this.polypeptide.monomerLengths - 1)' * states.monomerProductions;

            %energy
            translationEnergy = ...
                2 * this.polypeptide.monomerLengths' * states.monomerProductions + ...  %elongation
                3 * sum(states.monomerProductions);                                     %initiation, termination
            bmProd(this.substrateIndexs_gtp)       = bmProd(this.substrateIndexs_gtp)       + translationEnergy;
            bmProd(this.substrateIndexs_water)     = bmProd(this.substrateIndexs_water)     + translationEnergy;
            byProd(this.substrateIndexs_gdp)       = byProd(this.substrateIndexs_gdp)       + translationEnergy;
            byProd(this.substrateIndexs_phosphate) = byProd(this.substrateIndexs_phosphate) + translationEnergy;
            byProd(this.substrateIndexs_hydrogen)  = byProd(this.substrateIndexs_hydrogen)  + translationEnergy;

            %% enzymes

            %Ribosome
            fractionInitiating = this.ribosomeElongationRate *  ...
                states.monomerProductions0' * (1 ./ this.polypeptide.monomerLengths) / ...
                sum(states.monomerProductions0); %fraction of ribosomes that are initiating
            aaPolymerized0 = (states.monomerProductions0' * this.polypeptide.monomerLengths);
            minEnzExp(this.enzymeIndexs_ribosome30S) = 1.65 * aaPolymerized0 / this.ribosomeElongationRate / (1 - fractionInitiating);
            minEnzExp(this.enzymeIndexs_ribosome50S) = 1.65 * aaPolymerized0 / this.ribosomeElongationRate / (1 - fractionInitiating);

            %initiation, elongation, termination factors
            minEnzExp(this.enzymeIndexs_initiationFactors)     = 3 * minEnzExp(this.enzymeIndexs_ribosome30S) *      fractionInitiating;
            minEnzExp(this.enzymeIndexs_elongationFactors)     = 3 * minEnzExp(this.enzymeIndexs_ribosome30S) * (1 - fractionInitiating);
            minEnzExp(this.enzymeIndexs_terminationFactor)     = 3;
            minEnzExp(this.enzymeIndexs_recyclingFactor)       = 3;
            minEnzExp(this.enzymeIndexs_elongationGFactor)     = 3 + minEnzExp(this.enzymeIndexs_elongationGFactor);

            %abortive translation enzymes
            minEnzExp(this.enzymeIndexs_tmRNA)                 = 3;
            minEnzExp(this.enzymeIndexs_tmRNABindingProtein)   = 3;
            minEnzExp(this.enzymeIndexs_peptidylTRNAHydrolase) = 3;
        end

        %Initialize ribosome, bound mRNAs, translation factor states
        %Assumptions:
        %- Metabolic cost of the transcription of these amino acids is negligible
        %- Probability that so many ribosomes bind a given mRNA that they can't
        %  pack along the mRNA without steric interactions is negligible. Thus we
        %  don't try to handle this case separately.
        function initializeState(this)
            %references to states
            rib = this.ribosome;
            pol = this.polypeptide;

            %Check ribosomes and polypeptides are synchronized
            if ...
                    ~isequal(rib.boundMRNAs, pol.boundMRNAs) || ...
                    ~isequal(rib.mRNAPositions, pol.nascentMonomerLengths) || ...
                    ~isequal(rib.tmRNAPositions, pol.proteolysisTagLengths)
                throw(MException('Translation:error', 'Ribosomes and polypeptides must be synchronized'));
            end

            %enzymes
            this.enzymes = this.enzymes + this.boundEnzymes;
            this.boundEnzymes(:) = 0;
            this.enzymes(this.enzymeIndexs_ribosome30S) = this.enzymes(this.enzymeIndexs_ribosome30S) + this.enzymes(this.enzymeIndexs_ribosome70S);
            this.enzymes(this.enzymeIndexs_ribosome50S) = this.enzymes(this.enzymeIndexs_ribosome50S) + this.enzymes(this.enzymeIndexs_ribosome70S);
            this.enzymes(this.enzymeIndexs_ribosome70S) = 0;
            this.enzymes(this.enzymeIndexs_ribosome30S) = this.enzymes(this.enzymeIndexs_ribosome30S) + this.enzymes(this.enzymeIndexs_ribosome30SIF3);
            this.enzymes(this.enzymeIndexs_initiationFactor3) = this.enzymes(this.enzymeIndexs_initiationFactor3) + this.enzymes(this.enzymeIndexs_ribosome30SIF3);
            this.enzymes(this.enzymeIndexs_ribosome30SIF3) = 0;

            freeTranslationFactors  = this.enzymes(this.enzymeIndexs_translationFactors);
            boundTranslationFactors = this.boundEnzymes(this.enzymeIndexs_translationFactors);
            ribosome30S             = this.enzymes(this.enzymeIndexs_ribosome30S);
            ribosome30SIF3          = this.enzymes(this.enzymeIndexs_ribosome30SIF3);
            ribosome50S             = this.enzymes(this.enzymeIndexs_ribosome50S);
            ribosome70S             = this.enzymes(this.enzymeIndexs_ribosome70S);

            boundRibosome70S = min([
                ribosome30S;
                ribosome50S;
                freeTranslationFactors(this.enzymeIndexs_elongationFactors);
                sum(this.mRNAs)]);
            ribosome30S = ribosome30S - boundRibosome70S;
            ribosome50S = ribosome50S - boundRibosome70S;
            freeTranslationFactors(this.enzymeIndexs_elongationFactors)  = freeTranslationFactors(this.enzymeIndexs_elongationFactors)  - boundRibosome70S;
            boundTranslationFactors(this.enzymeIndexs_elongationFactors) = boundTranslationFactors(this.enzymeIndexs_elongationFactors) + boundRibosome70S;

            %allocate
            rib.states                = repmat(rib.notExistValue, [2*boundRibosome70S, 1, 1]);
            pol.boundMRNAs            = zeros(2 * boundRibosome70S, 1, 1);
            pol.nascentMonomerLengths = zeros(2 * boundRibosome70S, 1, 1);
            pol.proteolysisTagLengths = zeros(2 * boundRibosome70S, 1, 1);

            %bind random mRNA
            for i = 1:boundRibosome70S
                pol.boundMRNAs(i) = this.randStream.randsample(numel(this.mRNAs), 1, true, this.mRNAs .* pol.monomerLengths);
            end

            %actively bound
            %- partition sequence among RNA polymerases bound to transcription unit
            %- randomly select position of RNA polymerases within partition
            for i = 1:length(this.mRNAs)
                boundRibosomes = find(pol.boundMRNAs == i);
                partitions = round(...
                    (0:length(boundRibosomes)) * ...
                    (pol.monomerLengths(i) + this.ribosomeElongationRate - 1) / length(boundRibosomes) -...
                    1/2*this.ribosomeElongationRate);

                for j = 1:length(boundRibosomes)
                    minState = partitions(j)   + 1/2 * this.ribosomeElongationRate;
                    maxState = partitions(j+1) - 1/2 * this.ribosomeElongationRate;
                    rib.states(boundRibosomes(j)) = rib.activeValue;
                    if maxState < minState
                        pol.nascentMonomerLengths(boundRibosomes(j)) = round(max(1, min(pol.monomerLengths(i), mean(partitions(j:j+1)))));
                    else
                        pol.nascentMonomerLengths(boundRibosomes(j)) = round(minState + this.randStream.rand * (maxState - minState));
                    end
                end
            end

            %synchronize ribosome and polypeptide states
            rib.boundMRNAs = pol.boundMRNAs;
            rib.mRNAPositions = pol.nascentMonomerLengths;
            rib.tmRNAPositions = pol.proteolysisTagLengths;

            %store enzymes
            this.enzymes(this.enzymeIndexs_translationFactors)      = freeTranslationFactors;
            this.boundEnzymes(this.enzymeIndexs_translationFactors) = boundTranslationFactors;
            this.enzymes(this.enzymeIndexs_ribosome30S)             = ribosome30S;
            this.enzymes(this.enzymeIndexs_ribosome30SIF3)          = ribosome30SIF3;
            this.enzymes(this.enzymeIndexs_ribosome50S)             = ribosome50S;
            this.enzymes(this.enzymeIndexs_ribosome70S)             = ribosome70S;
            this.boundEnzymes(this.enzymeIndexs_ribosome70S)        = boundRibosome70S;

            %decrement counts of protein
            this.copyToState();

            g = this.gene;
            rna = this.rna;
            mon = this.monomer;
            cpx = this.complex;
            pcComp = sum(cpx.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3);

            initMonCnts = mon.counts;
            initCpxCnts = cpx.counts;
            matureProteinWt = max(0, sum(mon.dryWeight + cpx.dryWeight) - pol.dryWeight);
            while sum(mon.dryWeight + cpx.dryWeight) > matureProteinWt
                idx = this.randStream.randsample(numel(mon.counts) + numel(cpx.counts), 1, false, [mon.counts(:); cpx.counts(:)]);
                if idx <= numel(mon.counts)
                    mon.counts(idx) = mon.counts(idx) - 1;
                else
                    idx = idx - numel(mon.counts);
                    [i, ~] = ind2sub(size(cpx.counts), idx);
                    i = mod(i - 1, numel(cpx.matureIndexs)) + 1;
                    if any(pcComp(:, i))
                        rna.counts(rna.matureIndexs(setdiff(1:end, rna.matureMRNAIndexs))) = ...
                            + rna.counts(rna.matureIndexs(setdiff(1:end, rna.matureMRNAIndexs))) ...
                            + pcComp(:, i);
                    end
                    cpx.counts(idx) = cpx.counts(idx) - 1;
                end
            end
            mon.counts = ...
                + mon.counts ...
                + mon.updateExternalState(-(initMonCnts - mon.counts), true);
            cpx.counts = ...
                + cpx.counts ...
                + cpx.updateExternalState(-(initCpxCnts - cpx.counts), true);

            this.copyFromState();
        end

        %resource requirements
        function result = calcResourceRequirements_Current(this)
            initiationFactor3 = this.enzymes(this.enzymeIndexs_initiationFactor3);
            ribosome30S       = this.enzymes(this.enzymeIndexs_ribosome30S);
            ribosome30SIF3    = this.enzymes(this.enzymeIndexs_ribosome30SIF3);
            ribosome50S       = this.enzymes(this.enzymeIndexs_ribosome50S);
            nRibosomes = ...
                this.enzymes(this.enzymeIndexs_ribosome70S) + ...
                this.boundEnzymes(this.enzymeIndexs_ribosome70S) + ...
                min(ribosome30SIF3 + min(ribosome30S, initiationFactor3), ribosome50S);

            nInitiations = nRibosomes;
            nElongations = min([
                sum(this.aminoacylatedTRNAs) + this.aminoacylatedTMRNA
                nRibosomes * this.ribosomeElongationRate]);
            nTerminations = nRibosomes;

            %energy for initiation, elongation, termination
            result = zeros(size(this.substrates));
            result(this.substrateIndexs_gtp)   = nInitiations + 2 * nElongations + 3 * nTerminations;
            result(this.substrateIndexs_water) = nInitiations + 2 * nElongations + 3 * nTerminations;

            %water to cleave peptidyl-tRNA bound
            result(this.substrateIndexs_water) = - nInitiations - nElongations;

            %water from peptide bond formation
            result(this.substrateIndexs_water) = + nElongations;
        end

        %simulation
        function evolveState(this)
            %% references to states
            rib = this.ribosome;
            pol = this.polypeptide;

            %% Check ribosomes and polypeptides are synchronized
            if ...
                    ~isequal(rib.boundMRNAs, pol.boundMRNAs) || ...
                    ~isequal(rib.mRNAPositions, pol.nascentMonomerLengths) || ...
                    ~isequal(rib.tmRNAPositions, pol.proteolysisTagLengths)
                throw(MException('Translation:error', 'Ribosomes and polypeptides must be synchronized'));
            end

            %% allocate memory and initialize variables

            % Numbers of enzymes
            freeTranslationFactors     = this.enzymes(this.enzymeIndexs_translationFactors);
            boundTranslationFactors    = this.boundEnzymes(this.enzymeIndexs_translationFactors);
            ribosome30S                = this.enzymes(this.enzymeIndexs_ribosome30S);
            ribosome30SIF3             = this.enzymes(this.enzymeIndexs_ribosome30SIF3);
            ribosome50S                = this.enzymes(this.enzymeIndexs_ribosome50S);
            ribosome70S                = this.enzymes(this.enzymeIndexs_ribosome70S);
            boundRibosome70S           = this.boundEnzymes(this.enzymeIndexs_ribosome70S);
            freeTMRNA                  = this.freeTMRNA; %#ok<*PROP>
            freeTMRNABindingProtein    = this.enzymes(this.enzymeIndexs_tmRNABindingProtein);
            freePeptidylTRNAHydrolase  = this.enzymes(this.enzymeIndexs_peptidylTRNAHydrolase);
            boundTMRNA                 = this.boundTMRNA;
            aminoacylatedTmRNA         = this.aminoacylatedTMRNA;

            % form ribosomes
            new_ribosome30SIF3 = min(ribosome30S, freeTranslationFactors(this.enzymeIndexs_initiationFactor3));
            ribosome30SIF3 = ribosome30SIF3 + new_ribosome30SIF3;
            ribosome30S    = ribosome30S - new_ribosome30SIF3;
            freeTranslationFactors(this.enzymeIndexs_initiationFactor3) = freeTranslationFactors(this.enzymeIndexs_initiationFactor3) - new_ribosome30SIF3;

            %allocate space to store states of new ribosomes
            if boundRibosome70S + ribosome30SIF3 > numel(rib.states)
                rib.states = [
                    rib.states;
                    repmat(rib.notExistValue, [boundRibosome70S + ribosome30SIF3 1])];
                pol.boundMRNAs = [
                    pol.boundMRNAs;
                    zeros([boundRibosome70S + ribosome30SIF3 1])];
                pol.nascentMonomerLengths = [
                    pol.nascentMonomerLengths;
                    zeros([boundRibosome70S + ribosome30SIF3 1])];
                pol.proteolysisTagLengths = [
                    pol.proteolysisTagLengths;
                    zeros([boundRibosome70S + ribosome30SIF3 1])];
            end

            %ribosomes
            elongationRate = this.ribosomeElongationRate;               %ribosome elongation rate (aa/s)
            actRibs = find(...                                          %indices of active ribosomes
                rib.states == rib.activeValue | ...
                rib.states == rib.stalledValue);
            notExistRibs = find(rib.states == rib.notExistValue);       %indices of free ribosomes 30SIF3
            elngRibs = false(size(rib.states, 1), 1);                   %indicate of ribosome actively translating
            elngSeqs = zeros(size(rib.states, 1), elongationRate);      %tRNA sequences of mRNAs segments that should be translated
            termRibs = false(size(rib.states, 1), 1);
            termRibs(rib.states == rib.activeValue, 1) = ...
                pol.monomerLengths(pol.boundMRNAs(rib.states == rib.activeValue, 1), 1) == ...
                pol.nascentMonomerLengths(rib.states == rib.activeValue, 1);
            termRibs = termRibs | pol.proteolysisTagLengths == pol.proteolysisTagLength;

            %mRNAs
            bndProbs = this.mRNAs;                                      %unnormalized probability of binding mRNAs

            %tRNAs
            nTRNAs = length(this.freeTRNAs);                            %number of tRNA species
            tRNAs  = this.freeTRNAs + this.aminoacylatedTRNAs;          %total numbers of tRNAs (free and aminoacylated)

            %energy
            energyAllocation = max(0, this.substrates(this.substrateIndexs_gtp));
            availableEnergy  = energyAllocation;                        %energy (GTP) allocated to translation

            %water
            availableWater = this.substrates(this.substrateIndexs_water);

            %% recycle elongation factors
            freeTranslationFactors(this.enzymeIndexs_elongationFactors) = ...
                freeTranslationFactors(this.enzymeIndexs_elongationFactors) + ...
                boundTranslationFactors(this.enzymeIndexs_elongationFactors);
            boundTranslationFactors(this.enzymeIndexs_elongationFactors) = 0;

            %% select bound mRNAs to elongate
            %elongate if all elongation factors available
            nElongatingRibosomes = min([
                numel(actRibs)
                boundRibosome70S
                floor(min(availableEnergy, availableWater) / 2)]);

            randomOrder = this.randStream.randperm(boundRibosome70S);
            for idx = 1:nElongatingRibosomes
                i = actRibs(randomOrder(idx));

                %Use elongation factors to elongate
                if ~termRibs(i)
                    if all(freeTranslationFactors(this.enzymeIndexs_elongationFactors))
                        freeTranslationFactors(this.enzymeIndexs_elongationFactors)  = freeTranslationFactors(this.enzymeIndexs_elongationFactors)  - 1;
                        boundTranslationFactors(this.enzymeIndexs_elongationFactors) = boundTranslationFactors(this.enzymeIndexs_elongationFactors) + 1;
                    else
                        continue;
                    end
                end

                %Extract portion of sequence to be translated
                switch rib.states(i)
                    case rib.activeValue
                        idxs = pol.nascentMonomerLengths(i) + 1:min([
                            pol.monomerLengths(pol.boundMRNAs(i));
                            pol.nascentMonomerLengths(i) + elongationRate]);
                        elngSeqs(i, :) = padarray(pol.monomerTRNASequences{pol.boundMRNAs(i)}(idxs)', [0 elongationRate-length(idxs)], 0, 'post');
                        elngRibs(i) = true;
                    case rib.stalledValue
                        idxs = pol.proteolysisTagLengths(i) + 1:min([
                            pol.proteolysisTagLength;
                            pol.proteolysisTagLengths(i) + elongationRate]);
                        elngSeqs(i, :) = padarray(pol.proteolysisTagTRNASequence(idxs)', [0 elongationRate-length(idxs)], 0, 'post');
                        elngRibs(i) = true;
                end
            end
            nElongatingRibosomes = sum(elngRibs);

            %% make/update free ribosomes, initiate translation
            %Ribosome 70S assembly proceeds as:
            %30S-IF3 + IF2 + IF1 --> 30S-IF1-IF2-IF3-mRNA --> 30S-IF1-IF2-IF3-mRNA-tRNA
            %--> 30S-50S-IF1-IF2-mRNA-tRNA --> 30S-50S-mRNA-tRNA
            %All of these steps occur quite rapidly, so in our model if enough factors,
            %subunits, RNAs exist, then a 70S ribosome can be formed and
            %translation initiated.
            %(tRNA is bound in elongation, not initiation.)

            %Initiate if
            %- 30S-IF3 and 50S available
            %- initiation factors IF1 IF2 available
            %- energy available
            %- ribosome binding factor available
            %- mRNAs available
            %- (already verified by knowledge base computeTRNASequences
            %   method that first codon is formylmethionine)
            nInitiatingRibosomes = min([
                ribosome30SIF3;
                ribosome50S;
                freeTranslationFactors(this.enzymeIndexs_initiationFactor1);
                freeTranslationFactors(this.enzymeIndexs_initiationFactor2);
                max(0, min(availableEnergy, availableWater) - 2 * nElongatingRibosomes);
                sum(this.mRNAs)]);
            if nInitiatingRibosomes > 0
                ribosome30SIF3   = ribosome30SIF3   - nInitiatingRibosomes;
                ribosome50S      = ribosome50S      - nInitiatingRibosomes;
                boundRibosome70S = boundRibosome70S + nInitiatingRibosomes;
                freeTranslationFactors(this.enzymeIndexs_initiationFactor3) = ...
                    freeTranslationFactors(this.enzymeIndexs_initiationFactor3) + ...
                    nInitiatingRibosomes;
                availableEnergy = availableEnergy - nInitiatingRibosomes;
                availableWater = availableWater - nInitiatingRibosomes;
                for idx = 1:nInitiatingRibosomes
                    i = notExistRibs(idx);

                    %randomly choose mRNA to bind
                    pol.boundMRNAs(i) = this.randStream.randsample(numel(bndProbs), 1, true, bndProbs);

                    %update state of ribosome
                    rib.states(i) = rib.activeValue;

                    %update binding probabilities to prevent additional ribosomes
                    %from initiating translation with same mRNA at same time
                    bndProbs(pol.boundMRNAs(i)) = max(0, bndProbs(pol.boundMRNAs(i)) - 1);
                end
            end

            %% Translate active sequences
            %compute progress of translation
            nInit = sum((pol.nascentMonomerLengths(elngRibs) + pol.proteolysisTagLengths(elngRibs)) == 0);
            availableEnergyWater = floor(min(availableEnergy, availableWater - nInit));
            [translationProgress, this.aminoacylatedTRNAs, translationCost, availableEnergyWater2] = ...
                edu.stanford.covert.cell.sim.util.polymerize(...
                elngSeqs(elngRibs, :), ...
                this.aminoacylatedTRNAs, 1:nTRNAs, 0, availableEnergyWater, 2, this.randStream);
            if any(translationProgress)
                %update water, gtp
                availableWater = ...
                    + availableWater ...
                    - (availableEnergyWater - availableEnergyWater2) ...
                    - sum((pol.nascentMonomerLengths(elngRibs) + pol.proteolysisTagLengths(elngRibs)) == 0 & translationProgress > 0);
                availableEnergy = availableEnergy - (availableEnergyWater - availableEnergyWater2);

                %advance positions of ribosomes
                pol.nascentMonomerLengths(elngRibs) = ...
                    pol.nascentMonomerLengths(elngRibs) + ...
                    translationProgress .* (rib.states(elngRibs) == rib.activeValue);
                if any(rib.states(elngRibs) == rib.stalledValue)
                    pol.proteolysisTagLengths(elngRibs) = ...
                        pol.proteolysisTagLengths(elngRibs) + ...
                        translationProgress .* (rib.states(elngRibs) == rib.stalledValue);
                end
            end

            %% Terminate, remaining in elongating state, or convert to stalled state
            %If these conditions are true
            %(1) ribosome has translated entire sequence,
            %(2) termination factor is available (pfrA)--dissociates
            %    peptide from complex needs 1GTP,
            %(3) ribosome recycling factor (frr) and EF-G (fusA) available
            %    (use 1GTP) to split off the 50S subunit,
            %(4) IF-3 (infC) available to split off mRNA and tRNA and
            %(5) energy is available (total 2GTPs)
            %
            %Then terminate translation:
            %- use energy to terminate
            %- release ribosomes which have completed translation
            %- increment expression of the corresponding protein
            %- IF3 remains bound to 30S

            randomOrder = this.randStream.randperm(length(elngRibs));
            tmpProg = zeros(size(elngRibs));
            tmpProg(elngRibs) = translationProgress;
            for idx = 1:length(elngRibs)
                i = randomOrder(idx);
                if elngRibs(i) && (tmpProg(i) || pol.nascentMonomerLengths(i) >= pol.monomerLengths(pol.boundMRNAs(i)))
                    if (pol.nascentMonomerLengths(i) >= pol.monomerLengths(pol.boundMRNAs(i)) || ...
                            pol.proteolysisTagLengths(i) >= pol.proteolysisTagLength) && ...
                            freeTranslationFactors(this.enzymeIndexs_terminationFactor) && ...
                            freeTranslationFactors(this.enzymeIndexs_recyclingFactor) && ...
                            (freeTranslationFactors(this.enzymeIndexs_elongationGFactor) || tmpProg(i)) && ...
                            freeTranslationFactors(this.enzymeIndexs_initiationFactor3) && ...
                            availableEnergy >= 2 && ...
                            availableWater >= 2

                        %use energy, (tRNAs have already been released in elongation)
                        availableEnergy = availableEnergy - 2;
                        availableWater = availableWater - 2;

                        switch rib.states(i)
                            case rib.activeValue
                                %update protein count
                                this.monomers(pol.boundMRNAs(i)) = this.monomers(pol.boundMRNAs(i)) + 1;
                            case rib.stalledValue
                                %update stalled polypeptide count
                                pol.abortedPolypeptides = [
                                    pol.abortedPolypeptides;
                                    pol.boundMRNAs(i) pol.nascentMonomerLengths(i) pol.proteolysisTagLength];

                                %release tmRNA
                                boundTMRNA = boundTMRNA - 1;
                                freeTMRNA  = freeTMRNA  + 1;
                        end

                        %release ribosome
                        rib.states(i) = rib.notExistValue;

                        %release mRNA
                        pol.nascentMonomerLengths(i) = 0;
                        pol.proteolysisTagLengths(i) = 0;
                        pol.boundMRNAs(i)            = 0;

                        %release ribosome, split into 30S and 50S
                        boundRibosome70S = boundRibosome70S - 1;
                        ribosome30S      = ribosome30S      + 1;
                        ribosome50S      = ribosome50S      + 1;
                    end
                elseif rib.states(i) == rib.activeValue && ...
                        pol.nascentMonomerLengths(i) > 0 && ...
                        aminoacylatedTmRNA > 0 && ...
                        freeTMRNABindingProtein > 0  && ...
                        freePeptidylTRNAHydrolase > 0 && ...
                        this.randStream.rand < this.tmRNABindingProbability

                    %set ribosome to stalled state
                    rib.states(i) = rib.stalledValue;

                    %increment bound tmRNAs
                    boundTMRNA = boundTMRNA + 1;

                    %decrement free aminoacylated tmRNAs
                    aminoacylatedTmRNA = aminoacylatedTmRNA - 1;

                    %increment proteolysis tag length by 1 for alanine
                    %bound to tRNA-like portion of tmRNA
                    pol.proteolysisTagLengths(i) = 1;
                end
            end

            %% synchronize ribosome and polypeptide states
            rib.boundMRNAs = pol.boundMRNAs;
            rib.mRNAPositions = pol.nascentMonomerLengths;
            rib.tmRNAPositions = pol.proteolysisTagLengths;

            %% store enzymes
            this.enzymes(this.enzymeIndexs_translationFactors)         = freeTranslationFactors;
            this.boundEnzymes(this.enzymeIndexs_translationFactors)    = boundTranslationFactors;
            this.enzymes(this.enzymeIndexs_ribosome30S)                = ribosome30S;
            this.enzymes(this.enzymeIndexs_ribosome30SIF3)             = ribosome30SIF3;
            this.enzymes(this.enzymeIndexs_ribosome50S)                = ribosome50S;
            this.enzymes(this.enzymeIndexs_ribosome70S)                = ribosome70S;
            this.boundEnzymes(this.enzymeIndexs_ribosome70S)           = boundRibosome70S;
            this.freeTMRNA                                             = freeTMRNA;
            this.enzymes(this.enzymeIndexs_tmRNABindingProtein)        = freeTMRNABindingProtein;
            this.enzymes(this.enzymeIndexs_peptidylTRNAHydrolase)      = freePeptidylTRNAHydrolase;
            this.boundTMRNA                                            = boundTMRNA;
            this.aminoacylatedTMRNA                                    = aminoacylatedTmRNA;

            %% acount for used substrates
            %aminoacylated tRNAs
            this.freeTRNAs = tRNAs - this.aminoacylatedTRNAs;

            %energy (energy-availableEnergy), water
            this.substrates(this.substrateIndexs_gtp)       = this.substrates(this.substrateIndexs_gtp)       - (energyAllocation - availableEnergy);
            this.substrates(this.substrateIndexs_water)     = availableWater;
            this.substrates(this.substrateIndexs_gdp)       = this.substrates(this.substrateIndexs_gdp)       + (energyAllocation - availableEnergy);
            this.substrates(this.substrateIndexs_phosphate) = this.substrates(this.substrateIndexs_phosphate) + (energyAllocation - availableEnergy);
            this.substrates(this.substrateIndexs_hydrogen)  = this.substrates(this.substrateIndexs_hydrogen)  + (energyAllocation - availableEnergy) + sum(translationCost);
        end
    end

    %get methods of dependent local state
    methods
        function value = getDryWeight(this)
            if size(this.monomers, 3) == 1
                value = this.getDryWeight@edu.stanford.covert.cell.sim.Process() + (...
                    this.monomer.molecularWeights(this.monomer.nascentIndexs)' * this.monomers           + ...
                    this.rna.molecularWeights(this.rna.matureIndexs(this.rna.matureTRNAIndexs))'          * this.freeTRNAs          + ...
                    this.rna.molecularWeights(this.rna.aminoacylatedIndexs(this.rna.matureTRNAIndexs))'   * this.aminoacylatedTRNAs + ...
                    this.rna.molecularWeights(this.rna.aminoacylatedIndexs(this.tmRNAGlobalIndex))'       * this.aminoacylatedTMRNA) ...
                    / edu.stanford.covert.util.ConstantUtil.nAvogadro;
            else
                value = this.getDryWeight@edu.stanford.covert.cell.sim.Process() + (...
                    permute(this.monomer.molecularWeights(this.monomer.nascentIndexs)' * permute(this.monomers,           [1 3 2]),[1 3 2]) + ...
                    permute(this.rna.molecularWeights(this.rna.matureIndexs(this.rna.matureTRNAIndexs))'        * permute(this.freeTRNAs,          [1 3 2]),[1 3 2]) + ...
                    permute(this.rna.molecularWeights(this.rna.aminoacylatedIndexs(this.rna.matureTRNAIndexs))' * permute(this.aminoacylatedTRNAs, [1 3 2]),[1 3 2]) + ...
                    permute(this.rna.molecularWeights(this.rna.aminoacylatedIndexs(this.tmRNAGlobalIndex))'     * permute(this.aminoacylatedTMRNA, [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:
% Load whole cell simulation from disk/database.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 1/12/2011
function [simulation, states, options, parameters, fittedConstants, metadata] = loadSimulation(directoryOrWID, timeInit, timeFin, downsampleStepSec, downsampleType)
%initialize
setWarnings();
setPath();
setPreferences();

%import
import edu.stanford.covert.cell.sim.Simulation;
import edu.stanford.covert.cell.sim.util.DatabaseLogger;
import edu.stanford.covert.cell.sim.util.DiskLogger;
import edu.stanford.covert.cell.sim.util.CachedSimulationObjectUtil;
import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
import edu.stanford.covert.db.MySQLDatabase;

%load metadata
if ischar(directoryOrWID)
    directoryOrWID = SimulationDiskUtil.getSimulation(directoryOrWID);
    metadata = DiskLogger.loadMetadata(directoryOrWID);
else
    if ~exist('database', 'var')
        database = MySQLDatabase(getConfig());
    end
    metadata = DatabaseLogger.loadMetadata(database, directoryOrWID);
end

%construct simulation object
simulation = CachedSimulationObjectUtil.load(metadata.revision);

%load options, parameters, fitted constants based on data stored in disk/database
if ischar(directoryOrWID)
    if nargin < 2
        timeInit = [];
    end
    if nargin < 3
        timeFin = [];
    end
    if nargin < 4
        downsampleStepSec = [];
    end
    if nargin < 5
        downsampleType = [];
    end
    [states, metadata, options, parameters, fittedConstants, randStreamStates] = ...
        DiskLogger.load(directoryOrWID, '-independent', timeInit, timeFin, downsampleStepSec, downsampleType);
else
    if ~exist('database', 'var')
        database = MySQLDatabase(getConfig());
    end
    if nargin > 1
        warning('WholeCell:warning', 'Additional options ignored');
    end
    [states, metadata, options, parameters, fittedConstants, randStreamStates] = ...
        DatabaseLogger.load(simulation, database, directoryOrWID);
end

%apply options, parameters, constants, time courses to simulation
simulation.applyOptions(options);
simulation.applyParameters(parameters);
simulation.applyFittedConstants(fittedConstants);
simulation.applyRandStreamStates(randStreamStates);
simulation.allocateMemoryForState(numel(states.Time.values));

for i = 1:numel(simulation.states)
    state = simulation.states{i};
    stateID = state.wholeCellModelID(7:end);
    for j = 1:numel(state.stateNames)
        name = state.stateNames{j};
        try
            state.(name) = states.(stateID).(name);
        catch %#ok<CTCH>
            warning('WholeCell:warning', 'Data not provided for %s %s state', state.name, name);
        end
    end
end

for i = 1:numel(simulation.processes)
    process = simulation.processes{i};
    process.copyFromState();
end

simulation.applyPerturbationsToConstants();

%cleanup
if exist('database', 'var')
    database.close();
end


text:
function bsv_test03 ( )

%*****************************************************************************80
%
%% BSV_TEST03 varies the left boundary value ALPHA in the Burgers equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BSV_TEST03:\n' );
  fprintf ( 1, '  Solution of steady viscous Burgers equation.\n' );
  fprintf ( 1, '  Vary the left boundary condition ALPHA around the value +1.\n' );

  a = -1.0;
  b = +1.0;
  alpha_test = [ 0.96, 0.98, 0.99, 0.995, 1.0, 1.005, 1.010, 1.02, 1.04 ];
  beta = -1.0;
  nu = 0.1;
  test_num = length ( alpha_test );
  n = 161;
  output = 0;

  u = zeros(n,test_num);

  fprintf ( 1, '\n' );

  for test = 1 : test_num

    alpha = alpha_test(test);

    fprintf ( 1, '  Using ALPHA = %g\n', alpha );

    u(:,test) = bsv ( a, b, alpha, beta, nu, n, output );

  end

  x = linspace ( a, b, n );

  figure ( 3 )
  clf
  hold on
  plot ( x, u, 'b-', 'LineWidth', 3 );
  plot ( [a,b], [0.0,0.0], 'r-', 'LineWidth', 2 )
  grid on
  title ( 'ALPHA = 0.96, 0.98, 0.99, 0.995, 1, 1.005, 1.01, 1.02, 1.04' )
  xlabel ( '<--- X --->' )
  ylabel ( '<---U(X) --->' )
  axis ( [ a, b, -1.5, alpha ] )
  hold off
  filename = 'bsv_test03.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot to file "%s".\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Complex Step Diff Checking
clc
format long g
F = @(x) exp(x)/((cos(x))^3 + (sin(x))^3)
G = symJac(F);
H = symHess(F);
x0 = pi/2;
true1d = G(x0)
true2d = H(x0)
I = sqrt(2)/2*(1i + 1);

k = 6;
cs1d = imag(F(x0 + 1i*10^-k)/(1*10^-k))
cs2d = imag((F(x0 + I*10^-k) + F(x0 - I*10^-k))/(10^-k)^2)

%% Multivariable 2nd der
clc
F = @(x) exp(x(1))/((cos(x(2)))^3 + (sin(x(2)))^3)
G = symJac(F);
H = symHess(F);
x0 = [2;2];
true1d = G(x0)
true2d = H(x0)

cstepHess(F,x0)

%% Multivariable 2nd der 4x4
clc
f = @(x) x(1)*x(4)*sum(x(1:3)) + x(3);
g = @(x) [ x(1)*x(4) + x(4)*sum(x(1:3))
                x(1)*x(4)
                x(1)*x(4) + 1
                x(1)*sum(x(1:3)) ];
c = @(x) [ prod(x); sum(x.^2) ];
j = @(x) [ prod(x)./x 2*x ];
H = @(x,sigma,lambda) sigma*[ 2*x(4)             0      0   0;
                      x(4)               0      0   0;
                      x(4)               0      0   0;
                      2*x(1)+x(2)+x(3)  x(1)  x(1)  0 ] +....
          lambda(1)*[    0          0         0         0;
                              x(3)*x(4)     0         0         0;
                              x(2)*x(4) x(1)*x(4)     0         0;
                              x(2)*x(3) x(1)*x(3) x(1)*x(2)     0  ] +...
          lambda(2)*diag([2 2 2 2]);

x0 = randn(4,1);
sig = randn(1,1); lambda = randn(2,1);
% x0 = [0.1;0.2;0.3;0.4];

%Gradient
cstepJac(f,x0) - g(x0)'
%Jacobian
cstepJac(c,x0) - j(x0)'
%Hessian Of Objective
cstepHess(f,x0,true) - H(x0,1,[0;0])
%Hessian Of Lagrangian
cstepHessLag(f,c,x0,sig,lambda,true) - H(x0,sig,lambda)

%% Speed Test
clc
La = 0.6604; Lh = 0.1778; m_b = 0.575; m_f = 0.575; Lw = 0.4699; m_w = 1.87; g = 9.81; Kf = 0.1188; u = [1;1];
F = @(x) (0.4e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * La ^ 2 * Lh ^ 3 * m_b * m_f ^ 2 + 0.4e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * La ^ 2 * Lh ^ 3 * m_b ^ 2 * m_f + 0.6e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * La ^ 2 * Lh * Lw ^ 2 * m_b * m_f * m_w + 0.4e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * La ^ 4 * Lh * m_b * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * La ^ 2 * Lh * Lw ^ 2 * m_f ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * Lh * Lw ^ 4 * m_b * m_w ^ 2 + 0.4e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * La ^ 4 * Lh * m_b ^ 2 * m_f + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * La ^ 2 * Lh * Lw ^ 2 * m_b ^ 2 * m_w + cos(x(2)) ^ 3 * cos(x(1)) * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.2e1 * cos(x(2)) * sin(x(2)) * Lh ^ 3 * Lw ^ 2 * m_b * m_f * m_w * cos(x(1)) ^ 2 + cos(x(2)) * sin(x(2)) * Lh ^ 3 * Lw ^ 2 * m_f ^ 2 * m_w * cos(x(1)) ^ 2 + cos(x(2)) * sin(x(2)) * Lh ^ 3 * Lw ^ 2 * m_b ^ 2 * m_w * cos(x(1)) ^ 2 + 0.4e1 * cos(x(2)) ^ 3 * cos(x(1)) * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * Lh * Lw ^ 4 * m_f * m_w ^ 2 - cos(x(2)) ^ 3 * cos(x(1)) * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w - 0.4e1 * cos(x(2)) ^ 3 * cos(x(1)) * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b ^ 2 * m_f) / Lh / cos(x(1)) ^ 2 / (0.4e1 * La ^ 4 * m_b ^ 2 * m_f + 0.4e1 * La ^ 4 * m_b * m_f ^ 2 + Lw ^ 4 * m_b * m_w ^ 2 + Lw ^ 4 * m_f * m_w ^ 2 + 0.6e1 * La ^ 2 * Lw ^ 2 * m_b * m_f * m_w + La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + 0.2e1 * Lh ^ 2 * Lw ^ 2 * m_b * m_f * m_w + Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w) * x(4) ^ 2 + ((0.8e1 * cos(x(2)) ^ 2 * cos(x(1)) * La ^ 2 * Lh ^ 3 * m_b ^ 2 * m_f + 0.8e1 * cos(x(2)) ^ 2 * cos(x(1)) * La ^ 2 * Lh ^ 3 * m_b * m_f ^ 2 + 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) * Lh ^ 3 * Lw ^ 2 * m_b ^ 2 * m_w + 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) * Lh ^ 3 * Lw ^ 2 * m_f ^ 2 * m_w + 0.8e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 3 * La ^ 4 * Lh * m_b ^ 2 * m_f + 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 3 * Lh * Lw ^ 4 * m_f * m_w ^ 2 + 0.8e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 3 * La ^ 4 * Lh * m_b * m_f ^ 2 + 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 3 * Lh * Lw ^ 4 * m_b * m_w ^ 2 - 0.2e1 * cos(x(1)) * La ^ 2 * Lh * Lw ^ 2 * m_f ^ 2 * m_w - 0.2e1 * cos(x(1)) * La ^ 2 * Lh * Lw ^ 2 * m_b ^ 2 * m_w - 0.4e1 * cos(x(1)) * Lh ^ 3 * Lw ^ 2 * m_b * m_f * m_w - 0.8e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b * m_f ^ 2 + 0.12e2 * cos(x(2)) ^ 2 * cos(x(1)) ^ 3 * La ^ 2 * Lh * Lw ^ 2 * m_b * m_f * m_w + 0.8e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b ^ 2 * m_f + 0.4e1 * cos(x(2)) ^ 2 * cos(x(1)) * Lh ^ 3 * Lw ^ 2 * m_b * m_f * m_w + 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 3 * La ^ 2 * Lh * Lw ^ 2 * m_b ^ 2 * m_w + 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 3 * La ^ 2 * Lh * Lw ^ 2 * m_f ^ 2 * m_w - 0.12e2 * cos(x(1)) * La ^ 2 * Lh * Lw ^ 2 * m_b * m_f * m_w - 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w - 0.2e1 * cos(x(1)) * Lh ^ 3 * Lw ^ 2 * m_b ^ 2 * m_w - 0.2e1 * cos(x(1)) * Lh ^ 3 * Lw ^ 2 * m_f ^ 2 * m_w - 0.8e1 * cos(x(1)) * La ^ 2 * Lh ^ 3 * m_b ^ 2 * m_f - 0.8e1 * cos(x(1)) * La ^ 2 * Lh ^ 3 * m_b * m_f ^ 2 - 0.8e1 * cos(x(1)) * La ^ 4 * Lh * m_b ^ 2 * m_f - 0.8e1 * cos(x(1)) * La ^ 4 * Lh * m_b * m_f ^ 2 - 0.2e1 * cos(x(1)) * Lh * Lw ^ 4 * m_b * m_w ^ 2 - 0.2e1 * cos(x(1)) * Lh * Lw ^ 4 * m_f * m_w ^ 2) / Lh / cos(x(1)) ^ 2 / (0.4e1 * La ^ 4 * m_b ^ 2 * m_f + 0.4e1 * La ^ 4 * m_b * m_f ^ 2 + Lw ^ 4 * m_b * m_w ^ 2 + Lw ^ 4 * m_f * m_w ^ 2 + 0.6e1 * La ^ 2 * Lw ^ 2 * m_b * m_f * m_w + La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + 0.2e1 * Lh ^ 2 * Lw ^ 2 * m_b * m_f * m_w + Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w) * x(6) + (0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) * sin(x(1)) * Lh ^ 3 * Lw ^ 2 * m_b ^ 2 * m_w + 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) * sin(x(1)) * Lh ^ 3 * Lw ^ 2 * m_f ^ 2 * m_w + 0.8e1 * cos(x(2)) ^ 2 * cos(x(1)) * sin(x(1)) * La ^ 2 * Lh ^ 3 * m_b ^ 2 * m_f + 0.4e1 * cos(x(2)) ^ 2 * cos(x(1)) * sin(x(1)) * Lh ^ 3 * Lw ^ 2 * m_b * m_f * m_w + 0.8e1 * cos(x(2)) ^ 2 * cos(x(1)) * sin(x(1)) * La ^ 2 * Lh ^ 3 * m_b * m_f ^ 2) / Lh / cos(x(1)) ^ 2 / (0.4e1 * La ^ 4 * m_b ^ 2 * m_f + 0.4e1 * La ^ 4 * m_b * m_f ^ 2 + Lw ^ 4 * m_b * m_w ^ 2 + Lw ^ 4 * m_f * m_w ^ 2 + 0.6e1 * La ^ 2 * Lw ^ 2 * m_b * m_f * m_w + La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + 0.2e1 * Lh ^ 2 * Lw ^ 2 * m_b * m_f * m_w + Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w) * x(5)) * x(4) + (0.2e1 * cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * Lh ^ 3 * Lw ^ 2 * m_b * m_f * m_w - cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * La ^ 2 * Lh * Lw ^ 2 * m_b ^ 2 * m_w + cos(x(2)) ^ 3 * cos(x(1)) ^ 3 * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + cos(x(2)) * cos(x(1)) * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w - 0.2e1 * cos(x(2)) * cos(x(1)) ^ 3 * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + 0.2e1 * cos(x(2)) * cos(x(1)) ^ 3 * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w - cos(x(2)) ^ 3 * cos(x(1)) ^ 3 * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w - cos(x(2)) * cos(x(1)) * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w - cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * La ^ 2 * Lh * Lw ^ 2 * m_f ^ 2 * m_w - 0.4e1 * cos(x(2)) * sin(x(2)) * Lh ^ 3 * Lw ^ 2 * m_b * m_f * m_w * cos(x(1)) ^ 2 - 0.8e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * La ^ 2 * Lh ^ 3 * m_b * m_f ^ 2 - 0.8e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * La ^ 2 * Lh ^ 3 * m_b ^ 2 * m_f - 0.2e1 * cos(x(2)) * sin(x(2)) * Lh ^ 3 * Lw ^ 2 * m_f ^ 2 * m_w * cos(x(1)) ^ 2 - 0.2e1 * cos(x(2)) * sin(x(2)) * Lh ^ 3 * Lw ^ 2 * m_b ^ 2 * m_w * cos(x(1)) ^ 2 - cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * Lh * Lw ^ 4 * m_f * m_w ^ 2 + 0.4e1 * cos(x(2)) ^ 3 * cos(x(1)) ^ 3 * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b ^ 2 * m_f - 0.4e1 * cos(x(2)) ^ 3 * cos(x(1)) ^ 3 * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b * m_f ^ 2 - 0.4e1 * cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * La ^ 4 * Lh * m_b ^ 2 * m_f - 0.4e1 * cos(x(2)) * cos(x(1)) * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b * m_f ^ 2 + 0.4e1 * cos(x(2)) * cos(x(1)) * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b ^ 2 * m_f - 0.8e1 * cos(x(2)) * cos(x(1)) ^ 3 * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b ^ 2 * m_f + cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * Lh ^ 3 * Lw ^ 2 * m_b ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * Lh ^ 3 * Lw ^ 2 * m_f ^ 2 * m_w + 0.8e1 * cos(x(2)) * cos(x(1)) ^ 3 * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b * m_f ^ 2 - 0.4e1 * cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * La ^ 4 * Lh * m_b * m_f ^ 2 + 0.4e1 * cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * La ^ 2 * Lh ^ 3 * m_b ^ 2 * m_f + 0.4e1 * cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * La ^ 2 * Lh ^ 3 * m_b * m_f ^ 2 - cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * Lh * Lw ^ 4 * m_b * m_w ^ 2 - 0.6e1 * cos(x(2)) * cos(x(1)) ^ 4 * sin(x(2)) * La ^ 2 * Lh * Lw ^ 2 * m_b * m_f * m_w) / Lh / cos(x(1)) ^ 2 / (0.4e1 * La ^ 4 * m_b ^ 2 * m_f + 0.4e1 * La ^ 4 * m_b * m_f ^ 2 + Lw ^ 4 * m_b * m_w ^ 2 + Lw ^ 4 * m_f * m_w ^ 2 + 0.6e1 * La ^ 2 * Lw ^ 2 * m_b * m_f * m_w + La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + 0.2e1 * Lh ^ 2 * Lw ^ 2 * m_b * m_f * m_w + Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w) * x(6) ^ 2 + (-0.2e1 * cos(x(2)) * cos(x(1)) ^ 3 * La * Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w - 0.2e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * Lh ^ 3 * Lw ^ 2 * m_f ^ 2 * m_w - 0.2e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * Lh ^ 3 * Lw ^ 2 * m_b ^ 2 * m_w - 0.8e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * La ^ 2 * Lh ^ 3 * m_b ^ 2 * m_f + 0.2e1 * cos(x(2)) * cos(x(1)) ^ 3 * La * Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w - 0.8e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * La ^ 2 * Lh ^ 3 * m_b * m_f ^ 2 - 0.4e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * Lh ^ 3 * Lw ^ 2 * m_b * m_f * m_w - 0.8e1 * cos(x(2)) * cos(x(1)) ^ 3 * La ^ 3 * Lh ^ 2 * m_b ^ 2 * m_f + 0.2e1 * cos(x(2)) * cos(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w - 0.2e1 * cos(x(2)) * cos(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.8e1 * cos(x(2)) * cos(x(1)) ^ 3 * La ^ 3 * Lh ^ 2 * m_b * m_f ^ 2 - 0.8e1 * cos(x(2)) * cos(x(1)) * La ^ 3 * Lh ^ 2 * m_b * m_f ^ 2 + 0.8e1 * cos(x(2)) * cos(x(1)) * La ^ 3 * Lh ^ 2 * m_b ^ 2 * m_f) / Lh / cos(x(1)) ^ 2 / (0.4e1 * La ^ 4 * m_b ^ 2 * m_f + 0.4e1 * La ^ 4 * m_b * m_f ^ 2 + Lw ^ 4 * m_b * m_w ^ 2 + Lw ^ 4 * m_f * m_w ^ 2 + 0.6e1 * La ^ 2 * Lw ^ 2 * m_b * m_f * m_w + La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + 0.2e1 * Lh ^ 2 * Lw ^ 2 * m_b * m_f * m_w + Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w) * x(5) * x(6) + (-0.4e1 * cos(x(2)) * cos(x(1)) * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b * m_f ^ 2 - cos(x(2)) * cos(x(1)) * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + cos(x(2)) * cos(x(1)) * sin(x(1)) * La * Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + 0.4e1 * cos(x(2)) * cos(x(1)) * sin(x(1)) * La ^ 3 * Lh ^ 2 * m_b ^ 2 * m_f) / Lh / cos(x(1)) ^ 2 / (0.4e1 * La ^ 4 * m_b ^ 2 * m_f + 0.4e1 * La ^ 4 * m_b * m_f ^ 2 + Lw ^ 4 * m_b * m_w ^ 2 + Lw ^ 4 * m_f * m_w ^ 2 + 0.6e1 * La ^ 2 * Lw ^ 2 * m_b * m_f * m_w + La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + 0.2e1 * Lh ^ 2 * Lw ^ 2 * m_b * m_f * m_w + Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w) * x(5) ^ 2 + (-0.3e1 * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_b ^ 2 * m_f + sin(x(2)) * sin(x(1)) * g * Lh * Lw ^ 3 * m_b * m_w ^ 2 + sin(x(2)) * sin(x(1)) * g * Lh * Lw ^ 3 * m_f * m_w ^ 2 - cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 3 * Lw * m_b ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 3 * Lw * m_f ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f - cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w - cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w - cos(x(2)) * cos(x(1)) ^ 2 * g * La * Lw ^ 3 * m_b * m_w ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * g * La * Lw ^ 3 * m_f * m_w ^ 2 + cos(x(2)) * cos(x(1)) * g * Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w - cos(x(2)) * cos(x(1)) * g * Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_b ^ 3 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_f ^ 3 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b ^ 3 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_f ^ 3 + cos(x(2)) * cos(x(1)) ^ 3 * g * La ^ 3 * Lw * m_b ^ 2 * m_w - cos(x(2)) * cos(x(1)) ^ 3 * g * La ^ 3 * Lw * m_f ^ 2 * m_w - 0.4e1 * cos(x(2)) * cos(x(1)) ^ 3 * g * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f - 0.3e1 * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_b * m_f ^ 2 - 0.3e1 * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b ^ 2 * m_f - 0.3e1 * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b * m_f ^ 2 + sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_b ^ 2 * m_w + sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_f ^ 2 * m_w + 0.4e1 * cos(x(2)) * cos(x(1)) ^ 3 * g * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 3 * g * La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w - cos(x(2)) * cos(x(1)) ^ 3 * g * La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 3 * g * La * Lw ^ 3 * m_b * m_w ^ 2 - cos(x(2)) * cos(x(1)) ^ 3 * g * La * Lw ^ 3 * m_f * m_w ^ 2 + 0.2e1 * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lw ^ 2 * u(1) * m_f * m_w - 0.2e1 * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lw ^ 2 * u(2) * m_b * m_w - 0.2e1 * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lw ^ 2 * u(2) * m_f * m_w - sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_b ^ 2 - sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_f ^ 2 - sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_b ^ 2 - sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_f ^ 2 - sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_b ^ 2 - sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_f ^ 2 - sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_b ^ 2 - sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_b ^ 2 - cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_b ^ 2 - cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_f ^ 2 + 0.2e1 * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lw ^ 2 * u(1) * m_b * m_w + 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(1) * m_b * m_f - 0.2e1 * cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(2) * m_b * m_f - cos(x(1)) * Kf * La ^ 2 * Lw ^ 2 * u(1) * m_b * m_w - cos(x(1)) * Kf * La ^ 2 * Lw ^ 2 * u(2) * m_b * m_w + 0.2e1 * cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_b * m_f - 0.2e1 * cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_b * m_f + cos(x(1)) * Kf * La ^ 2 * Lw ^ 2 * u(1) * m_f * m_w + cos(x(1)) * Kf * La ^ 2 * Lw ^ 2 * u(2) * m_f * m_w + cos(x(1)) * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_b ^ 2 * cos(x(2)) ^ 2 + cos(x(1)) * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_b ^ 2 * cos(x(2)) ^ 2 + cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_b ^ 2 * cos(x(1)) ^ 2 + cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_f ^ 2 * cos(x(1)) ^ 2 - cos(x(1)) * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_f ^ 2 * cos(x(2)) ^ 2 - cos(x(1)) * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_f ^ 2 * cos(x(2)) ^ 2 - cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_b ^ 2 * cos(x(1)) ^ 2 - cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_f ^ 2 * cos(x(1)) ^ 2 - cos(x(1)) * g * La * Lw ^ 3 * m_b * m_w ^ 2 * cos(x(2)) ^ 2 - 0.4e1 * cos(x(2)) * cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 - cos(x(1)) * g * La * Lh ^ 2 * Lw * m_f ^ 2 * m_w + 0.4e1 * cos(x(2)) * cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f + cos(x(1)) * g * La * Lh ^ 2 * Lw * m_b ^ 2 * m_w + cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f * cos(x(2)) ^ 2 + cos(x(1)) * g * La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w * cos(x(2)) ^ 2 + cos(x(1)) * g * La ^ 3 * Lw * m_f ^ 2 * m_w * cos(x(2)) ^ 2 + cos(x(1)) * g * La * Lw ^ 3 * m_f * m_w ^ 2 * cos(x(2)) ^ 2 - cos(x(1)) * g * La ^ 3 * Lw * m_b ^ 2 * m_w * cos(x(2)) ^ 2 - cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 * cos(x(2)) ^ 2 - cos(x(1)) * g * La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w * cos(x(2)) ^ 2 + 0.2e1 * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_b * m_f - 0.2e1 * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_b * m_f + Kf * Lh ^ 2 * Lw ^ 2 * u(1) * m_b * m_w + Kf * Lh ^ 2 * Lw ^ 2 * u(1) * m_f * m_w - Kf * Lh ^ 2 * Lw ^ 2 * u(2) * m_b * m_w - Kf * Lh ^ 2 * Lw ^ 2 * u(2) * m_f * m_w + cos(x(2)) * cos(x(1)) ^ 3 * g * Lw ^ 4 * m_b * m_w ^ 2 - cos(x(2)) * cos(x(1)) ^ 3 * g * Lw ^ 4 * m_f * m_w ^ 2 - sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_b ^ 3 - sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_f ^ 3 - sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b ^ 3 - sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_f ^ 3 + cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 4 * m_b ^ 2 * m_f - cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 4 * m_b * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 2 * Lh ^ 2 * m_b ^ 3 - cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 2 * Lh ^ 2 * m_f ^ 3 - cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(2) * m_b ^ 2 - cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(2) * m_f ^ 2 - 0.2e1 * cos(x(2)) ^ 2 * Kf * Lh ^ 4 * u(1) * m_b * m_f + 0.2e1 * cos(x(2)) ^ 2 * Kf * Lh ^ 4 * u(2) * m_b * m_f + 0.2e1 * cos(x(1)) ^ 2 * Kf * La ^ 4 * u(1) * m_b * m_f - 0.2e1 * cos(x(1)) ^ 2 * Kf * La ^ 4 * u(2) * m_b * m_f - 0.2e1 * cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(1) * m_b * m_f + 0.2e1 * cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(2) * m_b * m_f + cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(1) * m_b ^ 2 + cos(x(2)) ^ 2 * cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(1) * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 4 * u(1) * m_b ^ 2 - cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 4 * u(1) * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 4 * u(2) * m_b ^ 2 - cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 4 * u(2) * m_f ^ 2 - cos(x(1)) * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_b ^ 2 - cos(x(1)) * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_b ^ 2 - cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_b ^ 2 - cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_f ^ 2 + cos(x(1)) * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_f ^ 2 + cos(x(1)) * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_f ^ 2 + cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_b ^ 2 + cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_f ^ 2 - cos(x(1)) * Kf * La ^ 4 * u(2) * m_f ^ 2 * cos(x(2)) ^ 2 + cos(x(1)) * Kf * La ^ 4 * u(1) * m_b ^ 2 * cos(x(2)) ^ 2 + cos(x(1)) * Kf * La ^ 4 * u(2) * m_b ^ 2 * cos(x(2)) ^ 2 - cos(x(1)) * Kf * La ^ 4 * u(1) * m_f ^ 2 * cos(x(2)) ^ 2 - cos(x(1)) * g * La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w - cos(x(1)) * g * La ^ 3 * Lw * m_f ^ 2 * m_w - cos(x(1)) * g * La * Lw ^ 3 * m_f * m_w ^ 2 + cos(x(1)) * g * La ^ 3 * Lw * m_b ^ 2 * m_w + cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + cos(x(1)) * g * La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + cos(x(1)) * g * La * Lw ^ 3 * m_b * m_w ^ 2 + cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_b ^ 3 * cos(x(2)) ^ 2 + cos(x(1)) * g * La ^ 4 * m_b ^ 2 * m_f * cos(x(2)) ^ 2 - cos(x(1)) * g * La ^ 4 * m_b * m_f ^ 2 * cos(x(2)) ^ 2 - cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_f ^ 3 * cos(x(2)) ^ 2 - cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f + sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_b ^ 2 * m_w + sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_f ^ 2 * m_w - sin(x(2)) * sin(x(1)) * g * La * Lh * Lw ^ 2 * m_b ^ 2 * m_w - sin(x(2)) * sin(x(1)) * g * La * Lh * Lw ^ 2 * m_f ^ 2 * m_w + 0.2e1 * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_b * m_f * m_w - cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_b ^ 2 * m_w - cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_f ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 3 * g * La * Lh ^ 2 * Lw * m_b ^ 2 * m_w - cos(x(2)) * cos(x(1)) ^ 3 * g * La * Lh ^ 2 * Lw * m_f ^ 2 * m_w + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_b ^ 3 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_f ^ 3 - cos(x(2)) * cos(x(1)) ^ 2 * g * La * Lh ^ 2 * Lw * m_b ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 2 * g * La * Lh ^ 2 * Lw * m_f ^ 2 * m_w - cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_b ^ 2 * m_f - cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_b * m_f ^ 2 + 0.3e1 * cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b ^ 2 * m_f + 0.3e1 * cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b * m_f ^ 2 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_f ^ 3 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b ^ 3 - 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_f ^ 2 - 0.2e1 * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_b * m_f + 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_f ^ 2 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_f ^ 2 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_f ^ 2 - 0.2e1 * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_b * m_f + 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_b ^ 2 - 0.2e1 * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_b * m_f + 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_f ^ 2 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_b ^ 2 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_b ^ 2 + 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_b ^ 2 - 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_b ^ 2 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_b ^ 2 - 0.2e1 * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_b * m_f + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_f ^ 2 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_b ^ 2 + cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lw ^ 2 * u(1) * m_b * m_w - cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lw ^ 2 * u(1) * m_f * m_w + cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lw ^ 2 * u(2) * m_b * m_w - cos(x(2)) * cos(x(1)) ^ 2 * Kf * La ^ 2 * Lw ^ 2 * u(2) * m_f * m_w - 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_f ^ 2 - 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_b ^ 2 + cos(x(1)) * Kf * La ^ 2 * Lw ^ 2 * u(1) * m_b * m_w * cos(x(2)) ^ 2 + cos(x(1)) * Kf * La ^ 2 * Lw ^ 2 * u(2) * m_b * m_w * cos(x(2)) ^ 2 - 0.2e1 * cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(1) * m_b * m_f * cos(x(1)) ^ 2 + 0.2e1 * cos(x(2)) ^ 2 * Kf * La ^ 2 * Lh ^ 2 * u(2) * m_b * m_f * cos(x(1)) ^ 2 - cos(x(1)) * Kf * La ^ 2 * Lw ^ 2 * u(1) * m_f * m_w * cos(x(2)) ^ 2 - cos(x(1)) * Kf * La ^ 2 * Lw ^ 2 * u(2) * m_f * m_w * cos(x(2)) ^ 2 + cos(x(1)) * g * La * Lh ^ 2 * Lw * m_f ^ 2 * m_w * cos(x(2)) ^ 2 - cos(x(1)) * g * La * Lh ^ 2 * Lw * m_b ^ 2 * m_w * cos(x(2)) ^ 2 + Kf * Lh ^ 4 * u(1) * m_b ^ 2 + Kf * Lh ^ 4 * u(1) * m_f ^ 2 - Kf * Lh ^ 4 * u(2) * m_b ^ 2 - Kf * Lh ^ 4 * u(2) * m_f ^ 2 + cos(x(1)) * g * La ^ 4 * m_f ^ 3 - cos(x(1)) * g * La ^ 4 * m_b ^ 3 + 0.2e1 * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_b * m_f * m_w - 0.2e1 * sin(x(2)) * sin(x(1)) * g * La * Lh * Lw ^ 2 * m_b * m_f * m_w - cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_b ^ 2 * m_w - cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_f ^ 2 * m_w - 0.2e1 * cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_b * m_f * m_w - cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_b ^ 2 * m_f - cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La ^ 3 * Lh * m_b * m_f ^ 2 + 0.3e1 * cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b ^ 2 * m_f + 0.3e1 * cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b * m_f ^ 2 - cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_b ^ 2 * m_w - cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_f ^ 2 * m_w - 0.4e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b ^ 2 * m_f - 0.4e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La * Lh ^ 3 * m_b * m_f ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_b ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_f ^ 2 * m_w - sin(x(2)) * sin(x(1)) * Kf * La * Lh * Lw ^ 2 * u(1) * m_b * m_w - sin(x(2)) * sin(x(1)) * Kf * La * Lh * Lw ^ 2 * u(1) * m_f * m_w - sin(x(2)) * sin(x(1)) * Kf * La * Lh * Lw ^ 2 * u(2) * m_b * m_w - sin(x(2)) * sin(x(1)) * Kf * La * Lh * Lw ^ 2 * u(2) * m_f * m_w - 0.2e1 * cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_b * m_f - 0.2e1 * cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_b * m_f + 0.2e1 * cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_b * m_f + 0.2e1 * cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_b * m_f + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_b ^ 2 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_f ^ 2 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_b ^ 2 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_f ^ 2 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_b ^ 2 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_f ^ 2 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_b ^ 2 + cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_f ^ 2 - Kf * La ^ 2 * Lh ^ 2 * u(2) * m_b ^ 2 - Kf * La ^ 2 * Lh ^ 2 * u(2) * m_f ^ 2 + 0.2e1 * Kf * Lh ^ 4 * u(1) * m_b * m_f - 0.2e1 * Kf * Lh ^ 4 * u(2) * m_b * m_f + Kf * La ^ 2 * Lh ^ 2 * u(1) * m_b ^ 2 + Kf * La ^ 2 * Lh ^ 2 * u(1) * m_f ^ 2 - cos(x(1)) * g * La ^ 4 * m_f ^ 3 * cos(x(2)) ^ 2 + cos(x(1)) * g * La ^ 4 * m_b ^ 3 * cos(x(2)) ^ 2 - cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_b ^ 3 - cos(x(1)) * g * La ^ 4 * m_b ^ 2 * m_f + cos(x(1)) * g * La ^ 4 * m_b * m_f ^ 2 + cos(x(1)) * g * La ^ 2 * Lh ^ 2 * m_f ^ 3 + cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 4 * m_b ^ 3 - cos(x(2)) * cos(x(1)) ^ 2 * g * La ^ 4 * m_f ^ 3 - cos(x(1)) * Kf * La ^ 4 * u(1) * m_b ^ 2 - cos(x(1)) * Kf * La ^ 4 * u(2) * m_b ^ 2 + cos(x(1)) * Kf * La ^ 4 * u(1) * m_f ^ 2 + cos(x(1)) * Kf * La ^ 4 * u(2) * m_f ^ 2 - cos(x(2)) ^ 2 * Kf * Lh ^ 4 * u(1) * m_b ^ 2 - cos(x(2)) ^ 2 * Kf * Lh ^ 4 * u(1) * m_f ^ 2 + cos(x(2)) ^ 2 * Kf * Lh ^ 4 * u(2) * m_b ^ 2 + cos(x(2)) ^ 2 * Kf * Lh ^ 4 * u(2) * m_f ^ 2 + cos(x(1)) ^ 2 * Kf * La ^ 4 * u(1) * m_b ^ 2 + cos(x(1)) ^ 2 * Kf * La ^ 4 * u(1) * m_f ^ 2 - cos(x(1)) ^ 2 * Kf * La ^ 4 * u(2) * m_b ^ 2 - cos(x(1)) ^ 2 * Kf * La ^ 4 * u(2) * m_f ^ 2 - cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(1) * m_b ^ 2 - cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(1) * m_f ^ 2 + cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(2) * m_b ^ 2 + cos(x(1)) ^ 2 * Kf * Lh ^ 4 * u(2) * m_f ^ 2 + cos(x(1)) ^ 2 * Kf * Lw ^ 4 * u(1) * m_w ^ 2 - cos(x(1)) ^ 2 * Kf * Lw ^ 4 * u(2) * m_w ^ 2 + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_b ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_f ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La * Lh * Lw ^ 2 * m_b ^ 2 * m_w + cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La * Lh * Lw ^ 2 * m_f ^ 2 * m_w + 0.2e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_b * m_f * m_w + 0.2e1 * cos(x(2)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_b * m_f * m_w - cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_b ^ 2 * m_w - cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_f ^ 2 * m_w - 0.2e1 * cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * Lh ^ 3 * Lw * m_b * m_f * m_w + 0.2e1 * cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(1) * m_b * m_f + 0.2e1 * cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh ^ 3 * u(2) * m_b * m_f + 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh * Lw ^ 2 * u(1) * m_b * m_w - 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh * Lw ^ 2 * u(1) * m_f * m_w - 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh * Lw ^ 2 * u(2) * m_b * m_w + 0.2e1 * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La * Lh * Lw ^ 2 * u(2) * m_f * m_w - 0.2e1 * cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(1) * m_b * m_f - 0.2e1 * cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * Kf * La ^ 3 * Lh * u(2) * m_b * m_f - 0.2e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_b * m_f * m_w - 0.2e1 * cos(x(2)) * cos(x(1)) ^ 2 * sin(x(2)) * sin(x(1)) * g * La * Lh * Lw ^ 2 * m_b * m_f * m_w + 0.2e1 * cos(x(2)) * cos(x(1)) * sin(x(2)) * sin(x(1)) * g * La ^ 2 * Lh * Lw * m_b * m_f * m_w) / Lh / cos(x(1)) ^ 2 / (0.4e1 * La ^ 4 * m_b ^ 2 * m_f + 0.4e1 * La ^ 4 * m_b * m_f ^ 2 + Lw ^ 4 * m_b * m_w ^ 2 + Lw ^ 4 * m_f * m_w ^ 2 + 0.6e1 * La ^ 2 * Lw ^ 2 * m_b * m_f * m_w + La ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + La ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b ^ 2 * m_f + 0.4e1 * La ^ 2 * Lh ^ 2 * m_b * m_f ^ 2 + 0.2e1 * Lh ^ 2 * Lw ^ 2 * m_b * m_f * m_w + Lh ^ 2 * Lw ^ 2 * m_b ^ 2 * m_w + Lh ^ 2 * Lw ^ 2 * m_f ^ 2 * m_w);
x0 = [1;1;1;1;1;1];

tic
mj = mklJac(F,x0);
toc
tic
aj = autoJac(F,x0);
toc
tic
cj = cstepJac(F,x0');
toc

aj-cj

%% Vector Fun Test
clc
%Function
i = (1:50)';
fun = @(x) x(1)*exp(-x(2)*i) + x(3);
x0 = ones(3,1);

tic
mj = mklJac(fun,x0);
toc
tic
aj = autoJac(fun,x0);
toc
tic
cj = cstepJac(fun,x0);
toc

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function grad = mapDataTerm_gradImage(SR, model, LR, W, Wt)

    if isempty(model.photometricParams)
        % Get residual error for super-resolved estimate.
        r = getResidual(SR, LR, W, model.photometricParams);
        % Weight residual error by associated confidence map
        confidence = model.confidence;
        if isempty(confidence)
            % Default: equal confidence for each pixel
            confidence = ones(size(LR));
        end
        if ~isvector(confidence)
            confidence = imageToVector(confidence);
        end
        % Calculate the gradient w.r.t. image pixels of HR image
        if strcmp(model.errorModel, 'l2NormErrorModel')
            grad = - 2 * (Wt * (confidence .* r));
        elseif strcmp(model.errorModel, 'l1NormErrorModel')
            grad = - Wt * (confidence .* sign(r));
        else
            % User-defined error model
            [~, g] = model.errorModel(r);
            grad = - Wt * (confidence .* g);
        end
    else
        % Calculate the residual using photometric parameters.
        numFrames = size(model.photometricParams.mult, 3);
        numLRPixel = length(LR)/numFrames;
        if isvector(model.photometricParams.mult(:,:,1))
            % Global, affine photometric model
            bm = zeros(size(LR));
            ba = zeros(size(LR));
            for k = 1:numFrames
                bm( ((k-1)*numLRPixel + 1):(k*numLRPixel) ) = repmat(model.photometricParams.mult(k), numLRPixel, 1);
                ba( ((k-1)*numLRPixel + 1):(k*numLRPixel) ) = repmat(model.photometricParams.add(k), numLRPixel, 1);
            end
        else
            % Local (pixel-wise) photometric model
            bm = zeros(size(LR));
            ba = zeros(size(LR));
            for k = 1:numFrames
                bm( ((k-1)*numLRPixel + 1):(k*numLRPixel) ) = imageToVector(model.photometricParams.mult(:,:,k));
                ba( ((k-1)*numLRPixel + 1):(k*numLRPixel) ) = imageToVector(model.photometricParams.add(:,:,k));
            end
        end
        r = LR - bm .* (W*SR) - ba;
        % Weight residual error by associated confidence map
        confidence = model.confidence;
        if ~isvector(confidence)
            confidence = imageToVector(confidence);
        end
        if isempty(confidence)
            % Default: equal confidence for each pixel
            confidence = ones(size(LR));
        end

        % Calculate the gradient w.r.t. image pixels of HR image taking
        % photometric parameters into account.
%         grad = 0;
%         for k = 1:numFrames
%             m = numLRPixel*(k-1) + 1;
%             %Wkt = Wt(:, m:(m + numLRPixel-1));
%             Wk = W(m:(m + numLRPixel-1), :);
%             ck = confidence(m:(m + numLRPixel-1));
%             rk = r(m:(m + numLRPixel-1));
%
%             if ~isvector(model.photometricParams.mult(:,:,k))
%                 % Local (pixel-wise) photometric model
%                 bm = imageToVector( model.photometricParams.mult(:,:,k) );
%             else
%                 % Global, affine photometric model
%                 bm = model.photometricParams.mult(k);
%             end
%         end

            if strcmp(model.errorModel, 'l2NormErrorModel')
                %grad = grad - 2 * Wkt * (bm .* ck .* rk);
                grad = (-2 * (bm .* confidence .* r)' * W)';
                %grad = grad + g';
            elseif strcmp(model.errorModel, 'l1NormErrorModel')
                %grad = grad - Wkt * (bm .* ck .* sign(rk));
                grad = - Wt * (bm .* confidence .* sign(r));
                %grad = grad + g';
            else
                % User-defined error model
                [~, g] = model.errorModel(rk);
                %grad = grad - Wkt * (bm .* ck .* g);
                grad = - Wt * (bm .* confidence .* g);
                %grad = grad + g';
            end
        end
    end

author: thomas-koehler

repo: SupER


text:
function x = is_1d_regression(x)
x = isnumeric(x) && size(x,2) == 1;

repo: REST


text:
function test_bug472

% MEM 2gb
% WALLTIME 00:13:01
% DEPENDENCY buffer.mexa64 buffer.mexmaci buffer.mexw64 buffer.mexglx buffer.mexmaci64 buffer.mexmac buffer.mexw32

% start without a buffer
ft_destroy_buffer

% number of attempts
cnt = 1;

% use default url
url = 'buffer://localhost:1972';

% dummy header
hdr.Fs = 256;
hdr.nChans = 100;

stopwatch = tic;
% run the test for 10 minutes
% after a few attempts, MATLAB crashes
while (toc(stopwatch)<600)
    disp(['counter: ' num2str(cnt)]); cnt = cnt + 1;
    ft_create_buffer(1972);
    ft_write_data(url,rand(hdr.nChans,10000), 'header', hdr, 'dataformat', 'fcdc_buffer', 'append', false);
    ft_read_header(url)
    ft_destroy_buffer
    pause(0.1); % probably meaningless
end

disp('the test completed without detected problems');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function X = tdma(A,B,C,D)
%TriDiagonal Matrix Algorithm (TDMA) or Thomas Algorithm
% A_i*X_(i-1) + B_i*X_i + C_i*X_(i+1) = D_i (where A_1 = 0, C_n = 0)
% A,B,C,D are input vectors. X is the solution, also a vector.

% Copyright 2013 The MathWorks, Inc.

Cp = C;
Dp = D;
n = length(A);
X = zeros(n,1);
% Performs Gaussian elimination
Cp(1) = C(1)/B(1);
Dp(1) = D(1)/B(1);
for i = 2:n
Cp(i) = C(i)/(B(i)-Cp(i-1)*A(i));
Dp(i) = (D(i)-Dp(i-1)*A(i))/(B(i)-Cp(i-1)*A(i));
end
% Backward substitution, since X(n) is known first.
X(n) = Dp(n);
for i = n-1:-1:1
X(i) = Dp(i)-Cp(i)*X(i+1);
end

author: Sable

repo: mcbench-benchmarks


text:
classdef TestBlender
    %TestBlender

    methods (Static)
        function test_1
            btypes = {'FeatherBlender', 'MultiBandBlender'};
            for i=1:numel(btypes)
                obj = cv.Blender(btypes{i});
                typename = obj.typeid();
                validateattributes(typename, {'char'}, {'row', 'nonempty'});

                %{
                %TODO
                obj.prepare(corners, sizes);
                obj.feed(img, mask, tl);
                [dst, dst_mask] = obj.blend();
                %}
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [rgb,points3d]=read_3d_pts_general(depthInpaint,K,depthInpaintsize,imageName,crop)
% Convert depth to point cloud. Each pixel in the depth map corresponds a
% point in 3D. Provided by Shuran Song
% depthInpaint: input depth
% K: camera intrinsic
% depthInpantsize: size(depthInpaint)
% imageName: put it emtpy(as far as I've used this function)
% crop: let it be 1

    %K is [fx 0 cx; 0 fy cy; 0 0 1];
    %K = frames.K;
    if ~isempty(K)
        cx = K(1,3); cy = K(2,3);
        fx = K(1,1); fy = K(2,2);
    else
        fx = 5.19e+02;
        fy = 5.19e+02;
        cx = 320;
        cy = 240;
    end
    invalid = depthInpaint==0;
    if ~isempty(imageName)
        rgb = im2double(imageName);
    else
        rgb = double(cat(3,zeros(depthInpaintsize(1),depthInpaintsize(2)),...
            ones(depthInpaintsize(1),depthInpaintsize(2)),...
            zeros(depthInpaintsize(1),depthInpaintsize(2))));
    end
    rgb = reshape(rgb, [], 3);
    %3D points
    [x,y] = meshgrid((1:depthInpaintsize(2))+ crop(2)-1, (1:depthInpaintsize(1))+ crop(1)-1);
    x3 = (x-cx).*depthInpaint*1/fx;
    y3 = (y-cy).*depthInpaint*1/fy;
    z3 = depthInpaint;
    points3d = [x3(:) -y3(:) z3(:)];

    points3d(invalid(:),:) = NaN;
end

author: zhirongw


text:
function d = fd05 ( p )

%*****************************************************************************80
%
%% FD05 is a signed distance function for the cylinder with a hole.
%
%  Modified:
%
%    15 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P(N,3), one or more points.
%
%    Output, real D(N), the signed distance of each point to the boundary of the region.
%
  r = sqrt ( p(:,1).^2 + p(:,2).^2 );
  z = p(:,3);

  d1 = r - 1.0;
  d2 = z - 1.0;
  d3 = - z - 1.0;
  d4 = sqrt ( d1.^2 + d2.^2 );
  d5 = sqrt ( d1.^2 + d3.^2 );

  d = dintersect ( dintersect ( d1, d2 ), d3 );
  ix = ( 0.0 < d1 ) & ( 0.0 < d2 );
  d(ix) = d4(ix);
  ix = ( 0.0 < d1 ) & ( 0.0 < d3 );
  d(ix) = d5(ix);

  d = ddiff ( d, dsphere ( p, 0.0, 0.0, 0.0, 0.5 ) );

  return
end


text:
function b = int(a,x,L,U)
% function B = int(A,X,L,U)
%
% DESCRIPTION
%   Element-by-element integration of a polynomial with respect
%   to a single variable.
%
% INPUTS
%   A: polynomial
%   X: Scalar polynomial variable [Optional with default X = A.varname{1}]
%   L: Lower limit of definite integral
%   U: Upper limit of definite integral
%
% OUTPUTS
%   B: polynomial
%
% SYNTAX
%   B = int(A,X)
%     Indefinite integral of the polynomial, A, with respect to X.
%     X should be a polynomial variable or string.   Integration is done
%     element-by-element if A is a matrix.
%   B = int(A,X,L,U)
%     Definite integral of A with respect to X from lower limit L to
%     upper limit U.
%   B = int(A,X,[L U]);
%     Equivalent to B = diff(A,X,L,U)
%
% EXAMPLE
%   pvar x y z;
%   a = 2*x^3 - 2*x*z^2 + 5*y*z;
%   b = int(a,x)
%   diff(b,x)-a
%   c = int(a,[0 1])
%
% See also: diff, jacobian

% 12/6/2010 PJS  Initial Coding

% Process Inputs/ Error Checking
if nargin==1
    x = a.varname{1};
    L = []; U = [];
elseif nargin==2
    if isa(x,'double')
        % B = int(A,[L,U])
        L = x(1);
        U = x(2);
        x = a.varname{1};
    else
        % B = int(A,X)
        L = []; U = [];
    end
elseif nargin==3
    if isa(x,'double')
        % B = int(A,L,U)
        U = L;
        L = x;
        x = a.varname{1};
    else
        % B = int(A,X,[L U])
        U = L(2);
        L = L(1);
    end
elseif nargin~=4
    error(['Invalid syntax for the "int" command. ' ...
        'Type "help int" for more information.'])
end

if ispvar(x) && length(x)==1
    x = x.varname{1};
elseif ~ischar(x)
    error('X must be a single polynomial variable or a string');
end

% Get polynomial info about A
acoef = a.coefficient;
adeg = a.degmat;
avar = a.varname;
sza = a.matdim;
Naterms = size(acoef,1);

% Find variable we are differentiating with respect to.
varnumb = find( strcmp(avar,x) );

% Perform indefinite integral
szb = sza;
if isempty(varnumb)
    % int a(y) dx = a(y)*x
    bcoef = acoef;
    bdeg = [adeg ones(Naterms,1)];
    bvar = [avar; x];
else
    % int a(y)*x^n dx = a(y)*x^(n+1) / (n+1)
    bdeg = adeg;
    nplus1 = bdeg(:,varnumb)+1;
    bdeg(:,varnumb) = nplus1;

    bcoef = acoef;
    bcoef = lrscale( bcoef, 1./nplus1 , []);

    bvar = avar;
end

chkval = 0; % skip validity check
b = polynomial(bcoef,bdeg,bvar,szb,chkval);

% Evaluate definite integral
if ~isempty(L)
    bU = subs(b,{x},U);
    bL = subs(b,{x},L);
    b = bU-bL;
end

% Combine any common terms
b = combine(b);

author: yu-jiang

repo: radpbook

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


text:
function cpu_timing ( n )

%*****************************************************************************80
%
%% CPU_TIMING measures execution time with CPUTIME.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CPU_TIMING:\n' );
  fprintf ( 1, '  MATLAB has a built in command called CPUTIME,\n' );
  fprintf ( 1, '  which measures the CPU time elapsed so far.\n' );
  fprintf ( 1, '  The difference of two readings measures the work done lately.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TASK #1:\n' )
  fprintf ( 1, '  Multiply two random dense matices of order N\n' );
  fprintf ( 1, '  using MATLAB''s A*B operation.\n' );

  A = rand ( n, n );
  B = rand ( n, n );

  t1 = cputime ( );
  C = A * B;
  t2 = cputime ( );

  f = n * n * n;
  w = f / 1000000;
  s = t2 - t1;
  r = w / s;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  MegaFLOPs = %d\n', w );
  fprintf ( 1, '  Elapsed CPU time  = %f\n', s );
  fprintf ( 1, '  Rate  = %f\n', r );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TASK #2:\n' )
  fprintf ( 1, '  Multiply two random dense matices of order N\n' );
  fprintf ( 1, '  using explicit loops.\n' );

  A = rand ( n, n );
  B = rand ( n, n );

  t1 = cputime ( );
  for i = 1 : n
    for k = 1 : n
      C(i,k) = 0;
      for j = 1 : n
        C(i,k) = C(i,k) + A(i,j) * B(j,k);
      end
    end
  end
  t2 = cputime ( );

  f = n * n * n;
  w = f / 1000000;
  s = t2 - t1;
  r = w / s;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  MegaFLOPs = %d\n', w );
  fprintf ( 1, '  Elapsed CPU time  = %f\n', s );
  fprintf ( 1, '  Rate  = %f\n', r );

  return
end