sbml simulation with groovy

sbml.groovy

/* settings */
def filename = './mapk.xml'
def duration = 4000.0
def dt = 0.1
/************/

def doc = new XmlParser().parse(filename)
def model = doc.model[0]
def species = [:]
def parameters = [:]
def reactions = [:]
def speciesIdList = null

model.listOfSpecies[0].species.each { species[it.@id] = new Double(it.@initialAmount) }
speciesIdList = species.keySet()

model.listOfReactions[0].reaction.each { re ->
  re.kineticLaw[0].listOfParameters?.getAt(0).parameter.each { p ->
    parameters[p.@id] = new Double(p.@value)
  }
  def reactants = re.listOfReactants[0].speciesReference*.@species
  def products = re.listOfProducts[0].speciesReference*.@species
  def math = { node ->
    def name = node.name().getLocalPart()
    switch (name) {
      case "ci":
        def key = node.value()[0]
        if (parameters[key] != null)
          return [new Tuple("constant", parameters[key])]
        else
          return [new Tuple("species", key)]
      case "cn":
        return [new Tuple("constant", new Double(node.value()[0]))]
      case ["plus", "minus", "times", "divide", "power"]:
        return [new Tuple("operator", name)]
      case "apply":
        def op = call(node.children()[0])
        def list = call(node.children()[1])
        (2..<node.children().size()).each { list += owner.call(node.children()[it]) + op }
        return list
      default:
        assert false
    }
  }(re.kineticLaw[0].math[0].children()[0])
  reactions[re.@id] = [reactants: reactants, products: products, math: math]
}

print "# time"
speciesIdList.each { print " ${it}" }
println ''

def opmap = [
  plus:   { r, l -> r + l },
  minus:  { r, l -> r - l },
  times:  { r, l -> r * l },
  divide: { r, l -> r / l },
  power:  { r, l -> r ** l }
]
for (t = 0.0; t <= duration; t += dt) {
  reactions.each { id, info ->
    def stack = []
    info.math.each { tuple ->
      def type = tuple.first()
      switch (type) {
        case "operator":
          def left = stack.pop()
          def right = stack.pop()
          stack += opmap[tuple.last()](right, left)
          break
        case "species":
          stack += species[tuple.last()]
          break
        case "constant":
          stack += tuple.last()
          break
        default:
          assert false
          break
      }
    }

    def diff = stack.pop() * dt
    info.reactants.each { species_next[it] -= diff }
    info.products.each { species_next[it] += diff }
  }
  species = species_next

  print t
  speciesIdList.each { print " ${species[it]}" }
  println ''
}