#!/usr/bin/env python3

import sys


# ----- grammar; this is new, needs lots of testing

from pyparsing import StringEnd,Literal,Word,ZeroOrMore,OneOrMore,Optional,Forward,alphas,nums

def label(s):
  return lambda x:[[s]+list(x)]

plus = Literal('+').suppress()
minus = Literal('-').suppress()
equals = Literal('=').suppress()
at = Literal('@').suppress()
superscript = (Literal('^')|Literal('**')).suppress()
arrow = Literal('->').suppress()
lparen = Literal('(').suppress()
rparen = Literal(')').suppress()

digits = Word(nums)
pointfloat = Optional(digits)+Literal('.')+digits|digits+Literal('.')
exponentfloat = (pointfloat|digits)+(Literal('e')|Literal('E'))+Optional(Literal('+')|Literal('-'))+digits
number = (exponentfloat|pointfloat|digits).setParseAction(lambda x:[['number',float(''.join(x))]])

name = Word(alphas,alphas+nums+'_').setParseAction(label('name'))

expr = Forward()
atom = lparen+expr+rparen|name|number
factor = Forward()
pmfactor = ((Literal('+')|Literal('-'))+factor).setParseAction(label('pmfactor'))
power = (atom+Optional(superscript+factor)).setParseAction(label('power'))
factor << (pmfactor|power)
term = (factor+ZeroOrMore((Literal('*')|Literal('/'))+factor)).setParseAction(label('term'))
expr << (term+ZeroOrMore((Literal('+')|Literal('-'))+term)).setParseAction(label('expr'))

compartment = (Literal('compartment').suppress()+name+equals+expr).setParseAction(label('compartment'))
parameter = (Literal('parameter').suppress()+name+equals+expr).setParseAction(label('parameter'))

inputs = (name+ZeroOrMore(plus+name)).setParseAction(label('inputs'))
outputs = (name+ZeroOrMore(plus+name)).setParseAction(label('outputs'))
reaction = (inputs+arrow+outputs+at+expr).setParseAction(label('reaction'))

command = (compartment|parameter|reaction)+StringEnd()


# ----- convert parsed strings to differential equations

P = {} # P[p] is value of parameter p (number, or None if undefined)
C = {} # C[c] is initial value in compartment c (number, or None if undefined)
D = {} # D[c] is derivative of compartment c (expression)

def evaluate(x):
  # convert reactions to D for inputs and outputs
  # plug in parameters
  # leave compartment names as symbols
  # could also fold constants here for speed later, but not bothering

  if x[0] == 'number':
    assert len(x) == 2
    return x

  if x[0] == 'name':
    assert len(x) == 2
    if x[1] in P:
      if P[x[1]] == None: raise Exception('parameter %s undefined'%x[1])
      return P[x[1]]
    if x[1] not in C: raise Exception('name %s undefined'%x[1])
    return x

  if x[0] == 'pmfactor':
    assert len(x) == 3
    return [x[0],x[1],evaluate(x[2])]

  if x[0] == 'power':
    if len(x) == 2:
      return evaluate(x[1])
    assert len(x) == 3
    return ['arith',evaluate(x[1]),'^',evaluate(x[2])]

  if x[0] in ['term','expr']:
    if len(x) == 2: return evaluate(x[1])
    operations = ['*','/'] if x[0] == 'term' else ['+','-']
    assert len(x) >= 2
    result = ['arith',evaluate(x[1])]
    for j in range(2,len(x),2):
      assert x[j] in operations
      assert len(x) >= j+2
      result += [x[j],evaluate(x[j+1])]
    return result

  if x[0] == 'compartment':
    assert len(x) == 3
    c = x[1]
    assert len(c) == 2
    assert c[0] == 'name'
    c = c[1]
    if c in P: raise Exception('parameter %s already defined'%c)
    if c in C: raise Exception('compartment %s already defined'%c)
    C[c] = evaluate(x[2])
    D[c] = ['number',float(0)]
    return

  if x[0] == 'parameter':
    assert len(x) == 3
    p = x[1]
    assert len(p) == 2
    assert p[0] == 'name'
    p = p[1]
    if p in P: raise Exception('parameter %s already defined'%p)
    if p in C: raise Exception('compartment %s already defined'%p)
    P[p] = evaluate(x[2])
    return

  if x[0] in ['inputs','outputs']:
    assert len(x) >= 2
    return [x[0]]+[evaluate(xj) for xj in x[1:]]

  if x[0] == 'reaction':
    assert len(x) == 4
    assert x[1][0] == 'inputs'
    assert x[2][0] == 'outputs'
    assert x[3][0] == 'expr'
    inputs = evaluate(x[1])[1:]
    outputs = evaluate(x[2])[1:]
    speed = evaluate(x[3])

    speed = ['arith',speed]
    for c in inputs: speed += ['*',c]
    if len(inputs) >= 2:
      population = ['arith']
      speed += ['/',['arith',['population'],'^',['number',len(inputs)-1]]]

    coeff = {c:0 for c in C}

    for c in inputs:
      assert len(c) == 2
      assert c[0] == 'name'
      c = c[1]
      if c not in C: raise Exception('compartment %s undefined'%c)
      coeff[c] -= 1

    for c in outputs:
      assert len(c) == 2
      assert c[0] == 'name'
      c = c[1]
      if c not in C: raise Exception('compartment %s undefined'%c)
      coeff[c] += 1

    for c in C:
      if coeff[c] == 1:
        D[c] = ['arith',D[c],'+',speed]
      elif coeff[c] == -1:
        D[c] = ['arith',D[c],'-',speed]
      elif coeff[c] != 0:
        D[c] = ['arith',D[c],'+',['arith',speed,'*',['number',float(coeff[c])]]]

    return

  raise Exception('undefined operation %s'%x)


# ----- read input

while True:
  line = sys.stdin.readline()
  if not line: break
  if '#' in line: line = line[:line.index('#')]
  line = line.strip()
  if line == '': continue

  evaluate(command.parseString(line)[0])

for c in C:
  print(c,D[c])


# ----- solve differential equations (approximately)

def fullyevaluate(x):
  # return value of x as float (not as ['number',float])
  # uses current[c] defined below

  if x[0] == 'number': return x[1]
  if x[0] == 'name': return current[x[1]]

  if x[0] == 'population': return sum(current[c] for c in C)

  if x[0] == 'pmfactor':
    assert x[1] in ['+','-']
    if x[1] == '-': return -fullyevaluate(x[2])
    return fullyevaluate(x[2])

  if x[0] == 'arith':
    result = fullyevaluate(x[1])
    for j in range(2,len(x),2):
      assert x[j] in ['+','-','*','/','^']
      if x[j] == '+':
        result = result+fullyevaluate(x[j+1])
      if x[j] == '-':
        result = result-fullyevaluate(x[j+1])
      if x[j] == '*':
        result = result*fullyevaluate(x[j+1])
      if x[j] == '/':
        result = result/fullyevaluate(x[j+1])
      if x[j] == '^':
        result = result**fullyevaluate(x[j+1])
    return result

  raise Exception('undefined operation %s'%x)

current = {c:fullyevaluate(C[c]) for c in C}

history = {c:[] for c in C}
historyday = []
maxpopulation = 0

day = 0
stopday = 365.25
while day < stopday:

  historyday += [day]
  for c in C: history[c] += [current[c]]

  population = sum(current[c] for c in C)
  maxpopulation = max(population,maxpopulation)

  change = {c:0 for c in current}

  for c in C:
    d = fullyevaluate(D[c])
    change[c] += d

  daydelta = 1/72.0

  assert daydelta > 0
  day += daydelta
  for c in current:
    current[c] += change[c]*daydelta


# ----- plotting

outputfile = 'epi.pdf'
  
import datetime
xmin = datetime.datetime(2020,1,1)
xmax = datetime.datetime(2020,1,1)+datetime.timedelta(days=stopday)
days = [xmin+datetime.timedelta(days=x) for x in historyday]

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from matplotlib.backends.backend_pdf import PdfPages

fig,ax1 = plt.subplots()
# ax2 = ax1.twinx() # in case plots are desired on other scales
# ax3 = ax2.twinx()
fig.set_size_inches(8,4)

ax1.set_xlim([xmin,xmax])
# ax2.set_xlim([xmin,xmax])
# ax3.set_xlim([xmin,xmax])

ax1.set_ylim([0,maxpopulation])
ax1.tick_params(axis='y',labelcolor='black',labelleft=False,labelright=True)
for c in C:
  ax1.plot(days,[y for y in history[c]],color='black')

ax1.xaxis.set_major_locator(mdates.MonthLocator())
ax1.xaxis.set_major_formatter(mdates.DateFormatter('%b'))
plt.setp(ax1.get_xticklabels(),rotation=90)

fig.tight_layout()
with PdfPages(outputfile) as pdf:
  pdf.savefig()
  plt.close()