Skip to main content
Global

Stub 5 Henry

  • Page ID
    37832
  • \( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)\(\newcommand{\AA}{\unicode[.8,0]{x212B}}\)

    !pip install python-copasi
    # -*- coding: utf-8 -*-
    # Copyright (C) 2017 by Pedro Mendes, Virginia Tech Intellectual 
    # Properties, Inc., University of Heidelberg, and University of 
    # of Connecticut School of Medicine. 
    # All rights reserved. 
    
    # Copyright (C) 2010 - 2016 by Pedro Mendes, Virginia Tech Intellectual 
    # Properties, Inc., University of Heidelberg, and The University 
    # of Manchester. 
    # All rights reserved. 
    
    
    
    
    #  
    #  This is an example on how to calculate and output the Jacobian matrix
    #  in COPASI
    #  
    from COPASI import *
    import sys
    
    
    MODEL_STRING = "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n <!-- Created by COPASI version 4.5.31 (Debug) on 2010-05-11 13:40 with libSBML version 4.1.0-b3. -->\n <sbml xmlns=\"http://www.sbml.org/sbml/level2/version4\" level=\"2\" version=\"4\">\n <model metaid=\"COPASI1\" id=\"Model_1\" name=\"New Model\">\n <listOfUnitDefinitions>\n <unitDefinition id=\"volume\" name=\"volume\">\n <listOfUnits>\n <unit kind=\"litre\" scale=\"-3\"/>\n </listOfUnits>\n </unitDefinition>\n <unitDefinition id=\"substance\" name=\"substance\">\n <listOfUnits>\n <unit kind=\"mole\" scale=\"-3\"/>\n </listOfUnits>\n </unitDefinition>\n <unitDefinition id=\"unit_0\">\n <listOfUnits>\n <unit kind=\"second\" exponent=\"-1\"/>\n </listOfUnits>\n </unitDefinition>\n </listOfUnitDefinitions>\n <listOfCompartments>\n <compartment id=\"compartment_1\" name=\"compartment\" size=\"1\"/>\n </listOfCompartments>\n <listOfSpecies>\n <species metaid=\"COPASI2\" id=\"species_1\" name=\"A\" compartment=\"compartment_1\" initialConcentration=\"1\"/>\n <species metaid=\"COPASI3\" id=\"species_2\" name=\"B\" compartment=\"compartment_1\" initialConcentration=\"0\"/>\n <species metaid=\"COPASI4\" id=\"species_3\" name=\"C\" compartment=\"compartment_1\" initialConcentration=\"0\"/>\n </listOfSpecies>\n <listOfReactions>\n <reaction metaid=\"COPASI5\" id=\"reaction_1\" name=\"reaction_1\" reversible=\"false\">\n <listOfReactants>\n <speciesReference species=\"species_1\"/>\n </listOfReactants>\n <listOfProducts>\n <speciesReference species=\"species_2\"/>\n </listOfProducts>\n <kineticLaw>\n <math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n <apply>\n <times/>\n <ci> compartment_1 </ci>\n <ci> k1 </ci>\n <ci> species_1 </ci>\n </apply>\n </math>\n <listOfParameters>\n <parameter id=\"k1\" name=\"k1\" value=\"0.2\" units=\"unit_0\"/>\n </listOfParameters>\n </kineticLaw>\n </reaction>\n <reaction metaid=\"COPASI6\" id=\"reaction_2\" name=\"reaction_2\" reversible=\"false\">\n <listOfReactants>\n <speciesReference species=\"species_2\"/>\n </listOfReactants>\n <listOfProducts>\n <speciesReference species=\"species_3\"/>\n </listOfProducts>\n <kineticLaw>\n <math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n <apply>\n <times/>\n <ci> compartment_1 </ci>\n <ci> k1 </ci>\n <ci> species_2 </ci>\n </apply>\n </math>\n <listOfParameters>\n <parameter id=\"k1\" name=\"k1\" value=\"0.1\" units=\"unit_0\"/>\n </listOfParameters>\n </kineticLaw>\n </reaction>\n </listOfReactions>\n </model>\n </sbml>\n"
    MODEL_STRING = "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n <!-- Created by COPASI version 4.5.31 (Debug) on 2010-05-11 13:40 with libSBML version 4.1.0-b3. -->\n <sbml xmlns=\"http://www.sbml.org/sbml/level2/version4\" level=\"2\" version=\"4\">\n <model metaid=\"COPASI1\" id=\"Model_1\" name=\"New Model\">\n <listOfUnitDefinitions>\n <unitDefinition id=\"volume\" name=\"volume\">\n <listOfUnits>\n <unit kind=\"litre\" scale=\"-3\"/>\n </listOfUnits>\n </unitDefinition>\n <unitDefinition id=\"substance\" name=\"substance\">\n <listOfUnits>\n <unit kind=\"mole\" scale=\"-3\"/>\n </listOfUnits>\n </unitDefinition>\n <unitDefinition id=\"unit_0\">\n <listOfUnits>\n <unit kind=\"second\" exponent=\"-1\"/>\n </listOfUnits>\n </unitDefinition>\n </listOfUnitDefinitions>\n <listOfCompartments>\n <compartment id=\"compartment_1\" name=\"compartment\" size=\"1\"/>\n </listOfCompartments>\n <listOfSpecies>\n <species metaid=\"COPASI2\" id=\"species_1\" name=\"A\" compartment=\"compartment_1\" initialConcentration=\"1\"/>\n <species metaid=\"COPASI3\" id=\"species_2\" name=\"B\" compartment=\"compartment_1\" initialConcentration=\"0\"/>\n <species metaid=\"COPASI4\" id=\"species_3\" name=\"C\" compartment=\"compartment_1\" initialConcentration=\"0\"/>\n </listOfSpecies>\n <listOfReactions>\n <reaction metaid=\"COPASI5\" id=\"reaction_1\" name=\"reaction_1\" reversible=\"false\">\n <listOfReactants>\n <speciesReference species=\"species_1\"/>\n </listOfReactants>\n <listOfProducts>\n <speciesReference species=\"species_2\"/>\n </listOfProducts>\n <kineticLaw>\n <math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n <apply>\n <times/>\n <ci> compartment_1 </ci>\n <ci> k1 </ci>\n <ci> species_1 </ci>\n </apply>\n </math>\n <listOfParameters>\n <parameter id=\"k1\" name=\"k1\" value=\"0.2\" units=\"unit_0\"/>\n </listOfParameters>\n </kineticLaw>\n </reaction>\n <reaction metaid=\"COPASI6\" id=\"reaction_2\" name=\"reaction_2\" reversible=\"false\">\n <listOfReactants>\n <speciesReference species=\"species_2\"/>\n </listOfReactants>\n <listOfProducts>\n <speciesReference species=\"species_3\"/>\n </listOfProducts>\n <kineticLaw>\n <math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n <apply>\n <times/>\n <ci> compartment_1 </ci>\n <ci> k1 </ci>\n <ci> species_2 </ci>\n </apply>\n </math>\n <listOfParameters>\n <parameter id=\"k1\" name=\"k1\" value=\"0.1\" units=\"unit_0\"/>\n </listOfParameters>\n </kineticLaw>\n </reaction>\n </listOfReactions>\n </model>\n </sbml>\n"
    
    
    def main():
        assert CRootContainer.getRoot() != None
        # create a new datamodel
        dataModel = CRootContainer.addDatamodel()
        assert dataModel != None
        assert CRootContainer.getDatamodelList().size() == 1
        # next we import a simple SBML model from a string
    
        # clear the message queue so that we only have error messages from the import in the queue
        CCopasiMessage.clearDeque()
        result=True
        try:
          result = dataModel.importSBMLFromString(MODEL_STRING)
        except:
            sys.stderr.write("Import of model failed miserably.\n")
            return
        # check if the import was successful
        mostSevere = CCopasiMessage.getHighestSeverity()
        # if it was a filtered error, we convert it to an unfiltered type
        # the filtered error messages have the same value as the unfiltered, but they
        # have the 7th bit set which basically adds 128 to the value
        mostSevere = mostSevere & 127
    
        # we assume that the import succeeded if the return value is True and
        # the most severe error message is not an error or an exception
        if result != True and  mostSevere < CCopasiMessage.ERROR:
            sys.stderr.write("Sorry. Model could not be imported.\n")
            return
    
        #
        # now we tell the model object to calculate the jacobian
        #
        model = dataModel.getModel()
        assert model != None
    
        if model != None:
            # running a task, e.g. a trajectory will automatically make sure that
            # the initial values are transferred to the current state before the calculation begins.
            # If we use low level calculation methods like the one to calculate the jacobian, we
            # have to make sure the the initial values are applied to the state
            model.applyInitialValues()
            # we need an array that stores the result
            # the size of the matrix does not really matter because
            # the calculateJacobian autoamtically resizes it to the correct
            # size
            jacobian=FloatMatrix()
            # the first parameter to the calculation function is a reference to
            # the matrix where the result is to be stored
            # the second parameter is the derivationFactor for the calculation
            # it basically represents a relative delta value for the calculation of the derivatives
            # the third parameter is a boolean indicating whether the jacobian should
            # be calculated from the reduced (true) or full (false) system
            model.getMathContainer().calculateJacobian(jacobian, 1e-12, False)
            # now we print the result
            # the jacobian stores the values in the order they are
            # given in the user order in the state template so it is not really straight
            # forward to find out which column/row corresponds to which species
            stateTemplate = model.getStateTemplate()
            # and we need the user order
            userOrder = stateTemplate.getUserOrder()
            # from those two, we can construct an new vector that contains
            # the names of the entities in the jacobian in the order in which they appear in
            # the jacobian
            nameVector=[]
            entity = None
            status=-1
    
            for i in range(0,userOrder.size()):
                entity = stateTemplate.getEntity(userOrder.get(i))
                assert entity != None
                # now we need to check if the entity is actually
                # determined by an ODE or a reaction
                status = entity.getStatus()
    
                if (status == CModelEntity.Status_ODE or (status == CModelEntity.Status_REACTIONS and entity.isUsed())):
                    nameVector.append(entity.getObjectName())
    
            assert len(nameVector) == jacobian.numRows()
            # now we print the matrix, for this we assume that no
            # entity name is longer then 5 character which is a save bet since
            # we know the model
            print ("Jacobian Matrix:")
            print ("")
            row = "%7s" % (" ")
    
            for i in range(0,len(nameVector)):
                row = row + " %7s" % (nameVector[i])
    
            print (row)
    
            for i in range(0,len(nameVector)):
                row = "%7s" % (nameVector[i])
    
                for j in range(0,len(nameVector)):
                    row = row + " %7.3f" % (jacobian.get(i,j))
    
                print (row)
    
            # we can also calculate the jacobian of the reduced system
            # in a similar way
            model.getMathContainer().calculateJacobian(jacobian, 1e-12, True)
            # this time generating the output is actually simpler because the rows
            # and columns are ordered in the same way as the independent variables of the state temple
            print ("")
            print ("")
            print ("Reduced Jacobian Matrix:")
            print ("")
            row = "%7s" % (" ")
            
            iMax = stateTemplate.getNumIndependent()
            
            for i in range(0,iMax):
              row = row + " %7s" % (stateTemplate.getIndependent(i).getObjectName())
    
            print (row)
    
            for i in range(0,iMax):
                row = "%7s" % (stateTemplate.getIndependent(i).getObjectName())
    
                for j in range(0,iMax):
                    row = row + " %7.3f" % (jacobian.get(i,j))
    
                print (row)
    
    if(__name__ == '__main__'):
       main()