Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Flux diagram for chosen species #80

Merged
merged 2 commits into from
Jun 13, 2012
Merged

Flux diagram for chosen species #80

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
70e0335
Flux diagram creation extended focused on a centralized species
Jun 13, 2012
7e65555
Merge remote-tracking branch 'origin/master' into fluxDiagram_Chosen_…
Jun 13, 2012
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 48 additions & 24 deletions generateFluxDiagram.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,15 +47,15 @@

################################################################################

def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, outputDirectory, speciesDirectory=None, settings = None):
def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, outputDirectory, centralSpecies = None, speciesDirectory=None, settings = None):
"""
For a given `reactionModel` and simulation results stored as arrays of
`times`, species `concentrations`, and `reactionRates`, generate a series
of flux diagrams as frames of an animation, then stitch them together into
a movie. The individual frames and the final movie are saved on disk at
`outputDirectory.`
"""

global maximumNodeCount, maximumEdgeCount, timeStep, concentrationTolerance, speciesRateTolerance
# Allow user defined settings for flux diagram generation if given
if settings:
maximumNodeCount = settings['maximumNodeCount']
Expand All @@ -70,6 +70,13 @@ def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, out
reactionList = reactionModel.core.reactions[:]
numReactions = len(reactionList)

#search index of central species:
if centralSpecies is not None:
for i, species in enumerate(speciesList):
if species.label == centralSpecies:
centralSpeciesIndex = i
break

# Compute the rates between each pair of species (big matrix warning!)
speciesRates = numpy.zeros((len(times),numSpecies,numSpecies), numpy.float64)
for index, reaction in enumerate(reactionList):
Expand All @@ -92,22 +99,40 @@ def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, out

# Determine the nodes and edges to keep
nodes = []; edges = []
for i in range(numSpecies*numSpecies):
productIndex, reactantIndex = divmod(speciesIndex[-i-1], numSpecies)
if reactantIndex > productIndex:
# Both reactant -> product and product -> reactant are in this list,
# so only keep one of them
continue
if maxSpeciesRates[reactantIndex, productIndex] == 0:
break
if reactantIndex not in nodes and len(nodes) < maximumNodeCount: nodes.append(reactantIndex)
if productIndex not in nodes and len(nodes) < maximumNodeCount: nodes.append(productIndex)
if len(nodes) > maximumNodeCount:
break
edges.append([reactantIndex, productIndex])
if len(edges) >= maximumEdgeCount:
break

if centralSpecies is None:
for i in range(numSpecies*numSpecies):
productIndex, reactantIndex = divmod(speciesIndex[-i-1], numSpecies)
if reactantIndex > productIndex:
# Both reactant -> product and product -> reactant are in this list,
# so only keep one of them
continue
if maxSpeciesRates[reactantIndex, productIndex] == 0:
break
if reactantIndex not in nodes and len(nodes) < maximumNodeCount: nodes.append(reactantIndex)
if productIndex not in nodes and len(nodes) < maximumNodeCount: nodes.append(productIndex)
if len(nodes) > maximumNodeCount:
break
edges.append([reactantIndex, productIndex])
if len(edges) >= maximumEdgeCount:
break
else:
nodes.append(centralSpeciesIndex)
for index, reaction in enumerate(reactionList):
for reactant, product in reaction.pairs:
reactantIndex = speciesList.index(reactant)
productIndex = speciesList.index(product)
if maxSpeciesRates[reactantIndex, productIndex] == 0:
break
if len(nodes) > maximumNodeCount or len(edges) >= maximumEdgeCount:
break
if reactantIndex == centralSpeciesIndex:
if productIndex not in nodes:
nodes.append(productIndex)
edges.append([reactantIndex, productIndex])
if productIndex == centralSpeciesIndex:
if reactantIndex not in nodes:
nodes.append(reactantIndex)
edges.append([reactantIndex, productIndex])
# Create the master graph
# First we're going to generate the coordinates for all of the nodes; for
# this we use the thickest pen widths for all nodes and edges
Expand Down Expand Up @@ -224,7 +249,7 @@ def simulate(reactionModel, reactionSystem, settings = None):
Generate and return a set of core and edge species and reaction fluxes
by simulating the given `reactionSystem` using the given `reactionModel`.
"""

global maximumNodeCount, maximumEdgeCount, timeStep, concentrationTolerance, speciesRateTolerance
# Allow user defined settings for flux diagram generation if given
if settings:
maximumNodeCount = settings['maximumNodeCount']
Expand Down Expand Up @@ -487,7 +512,7 @@ def loadRMGPyJob(inputFile, chemkinFile=None, speciesDict=None):

################################################################################

def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = False, settings = None, chemkinOutput = ''):
def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = False, settings = None, chemkinOutput = '', centralSpecies = None):
"""
Generates the flux diagram based on a condition 'inputFile', chemkin.inp chemkinFile,
a speciesDict txt file, plus an optional chemkinOutput file.
Expand All @@ -510,7 +535,7 @@ def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = Fals
time, coreSpeciesConcentrations, coreReactionRates, edgeReactionRates = loadChemkinOutput(chemkinOutput, rmg.reactionModel)

print 'Generating flux diagram for chemkin output...'
generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '1'), speciesPath, settings)
generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '1'), centralSpecies, speciesPath, settings)

else:
# Generate a flux diagram video for each reaction system
Expand All @@ -532,9 +557,8 @@ def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = Fals
time, coreSpeciesConcentrations, coreReactionRates, edgeReactionRates = simulate(rmg.reactionModel, reactionSystem, settings)

print 'Generating flux diagram for reaction system {0:d}...'.format(index+1)
generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates,
os.path.join(savePath, '{0:d}'.format(index+1)), speciesPath, settings)

generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '{0:d}'.format(index+1)),
centralSpecies, speciesPath, settings)
################################################################################

if __name__ == '__main__':
Expand Down