1# -*- coding: utf-8 -*-
2# Dilawar Singh <dilawar.s.rajput@gmail.com>, 2020-04-29
3# Python implementation of box.txt
4
5import smoldyn as sm
6import numpy as np
7import os.path
8
9def test_data_output():
10    s = sm.Simulation(low=[0, 0, 0], high=[100, 100, 100], boundary_type="ppp")
11
12    s.addBox(size=10)
13
14    # declaration of species A, B, and C with attributes.
15    a = s.addSpecies("A", state="soln", difc=1, color="red", mol_list="Alist")
16    b = s.addSpecies("B", state="soln", color="green", difc=1, mol_list="Blist")
17    c = s.addSpecies("C", state="soln", difc=1.0, color="blue", mol_list="Clist")
18    s.addMolecules(a, 1000)
19    s.addMolecules(b, 1000)
20
21    s.addBidirectionalReaction("r1", subs=[c], prds=(a, b), kf=0.1, kb=100)
22
23    # FIXME: Prints only upto 10 (2 iterations rather than 100)
24    s.setOutputFile("box.dat")
25    # s.setGraphics("opengl")
26    c = s.addCommand("molcount box.dat", "i", on=0, off=100, step=10)
27    s.run(100, dt=0.1, overwrite=True)
28    print("Simulation over")
29
30    # Now read the box.dat and verify the results.
31    # On OSX, it may not work. Homebrew python is a shell script which chdir to
32    # script dir before running it. I.e., box.dat would not be found easily.
33    if not os.path.exists('box.dat'):
34        quit(0)
35    data = np.loadtxt("box.dat")
36    assert data.shape == (11, 4), data.shape
37    expected = [50.036, 590.727, 590.727, 409.272]
38    print(data.mean(axis=0), expected)
39    assert np.isclose(
40        data.mean(axis=0), expected, atol=1e-1, rtol=1e-1
41    ).all(), data.mean(axis=0)
42    quit(0)
43
44
45def main():
46    test_data_output()
47
48
49if __name__ == "__main__":
50    main()
51