From 8788d45fe06eaa22516ee3b66953c5c18ee25d4c Mon Sep 17 00:00:00 2001 From: Dilawar Singh Date: Mon, 15 Nov 2021 14:49:29 +0530 Subject: [PATCH 1/2] Added first version. --- examples/S15_python/change_env.py | 21 +-- examples/S15_python/integrate_with_moose.py | 2 +- examples/S15_python/template.py | 27 ++-- examples/S7_surfaces/stick2.py | 14 +- examples/S99_more/maze/maze.txt | 134 ++++++++++++++++++++ examples/S99_more/maze/maze_solve.py | 53 ++++++++ source/python/smoldyn/smoldyn.py | 33 ++++- source/python/smoldyn/types.py | 13 +- 8 files changed, 261 insertions(+), 36 deletions(-) create mode 100644 examples/S99_more/maze/maze.txt create mode 100644 examples/S99_more/maze/maze_solve.py diff --git a/examples/S15_python/change_env.py b/examples/S15_python/change_env.py index 1b347bf4..0df5b044 100644 --- a/examples/S15_python/change_env.py +++ b/examples/S15_python/change_env.py @@ -8,6 +8,7 @@ r1_ = None bouton_ = None + def generate_spike(t, args): x = random.random() if x < 0.1: @@ -19,7 +20,7 @@ def update_kf(val): global r1_, bouton_ r1_.kf = 1000 ** val c = min(1.0, val) - bouton_.setStyle('both', color=[c, c, c]) + bouton_.setStyle("both", color=[c, c, c]) def build_model_smoldyn(): @@ -31,7 +32,9 @@ def build_model_smoldyn(): """ global r1_, bouton_ s = smoldyn.Simulation(low=[-500, -500], high=[1500, 1500]) - sv = s.addSpecies("SV", difc=dict(all=2400, front=10), color="blue", display_size=10) + sv = s.addSpecies( + "SV", difc=dict(all=2400, front=10), color="blue", display_size=10 + ) sv.addToSolution(100, lowpos=(0, 0), highpos=(1000, 1000)) # add a reaction with generates the sv at a fixed rate. Its better to @@ -47,17 +50,17 @@ def build_model_smoldyn(): decay = s.addReaction("decay", subs=[trans], prds=[], rate=math.log(2) / 20e-3) # BOUTON - path = s.addPath2D((1000, 0), (1000, 1000), (0, 1000), (0, 0)) + path = s.addPath2D([(1000, 0), (1000, 1000), (0, 1000), (0, 0)]) bouton_ = s.addSurface("bouton", panels=path.panels) - bouton_.setStyle('both', color="blue") - bouton_.setAction('both', [sv], "reflect") + bouton_.setStyle("both", color="blue") + bouton_.setAction("both", [sv], "reflect") # this is the bottom surface of bouton. This is sticky for synaptic # vesciles rect1 = smoldyn.Rectangle(corner=(0, 0), dimensions=[1000], axis="+y") bottom = s.addSurface("boutonBottom", panels=[rect1]) - bottom.setStyle('both', color="red") - bottom.setAction('back', trans, "reflect") # but it reflect neurotranmitter + bottom.setStyle("both", color="red") + bottom.setAction("back", trans, "reflect") # but it reflect neurotranmitter # SV stick to bottom of the bouton and also detach back with a smaller # rate. @@ -72,14 +75,14 @@ def build_model_smoldyn(): s.connect(generate_spike, update_kf, step=20) s.addGraphics("opengl", iter=10, text_display="time") - print('[INFO] Starting simulation ...') + print("[INFO] Starting simulation ...") s.run(stop=20, dt=0.001) print("Done") + def main(): build_model_smoldyn() if __name__ == "__main__": main() - diff --git a/examples/S15_python/integrate_with_moose.py b/examples/S15_python/integrate_with_moose.py index 52a91e5b..bc8862e2 100644 --- a/examples/S15_python/integrate_with_moose.py +++ b/examples/S15_python/integrate_with_moose.py @@ -86,7 +86,7 @@ def update_kf(val): s.addReaction("decay", subs=[trans], prds=[], rate=math.log(2)/20e-3) # BOUTON -path = s.addPath2D((1000, 0), (550, 1000), (450, 1000), (0, 0)) +path = s.addPath2D([(1000, 0), (550, 1000), (450, 1000), (0, 0)]) bouton_ = s.addSurface("bouton", panels=path.panels) bouton_.setStyle('both', color="blue") bouton_.setAction('both', [sv], "reflect") diff --git a/examples/S15_python/template.py b/examples/S15_python/template.py index 6736903c..d649275a 100644 --- a/examples/S15_python/template.py +++ b/examples/S15_python/template.py @@ -13,12 +13,12 @@ __author__ = "Dilawar Singh" __email__ = "dilawars@ncbs.res.in" -import smoldyn +import smoldyn # Model parameters -K_FWD = 0.001 # substrate-enzyme association reaction rate -K_BACK = 1 # complex dissociation reaction rate -K_PROD = 1 # complex reaction rate to product +K_FWD = 0.001 # substrate-enzyme association reaction rate +K_BACK = 1 # complex dissociation reaction rate +K_PROD = 1 # complex reaction rate to product # Simulation starts with declaring a Simulation object with the system boundaries. s = smoldyn.Simulation(low=[-1, -1], high=[1, 1]) @@ -36,30 +36,35 @@ # as well. Here, `all` molecules reflect at both surface faces. sph1 = smoldyn.Sphere(center=(0, 0), radius=1, slices=50) membrane = s.addSurface("membrane", panels=[sph1]) -membrane.setAction('both', [S, E, P, ES], "reflect") -membrane.setStyle('both', color="black", thickness=1) +membrane.setAction("both", [S, E, P, ES], "reflect") +membrane.setStyle("both", color="black", thickness=1) # Define a compartment, which is region inside the 'membrane' surface. inside = s.addCompartment(name="inside", surface=membrane, point=[0, 0]) # Chemical reactions. Here, E + S <-> ES -> P r1 = s.addBidirectionalReaction( - "r1", subs=[(E,"front"), (S,"bsoln")], prds=[(ES,"front")], kf=K_FWD, kb=K_BACK) + "r1", subs=[(E, "front"), (S, "bsoln")], prds=[(ES, "front")], kf=K_FWD, kb=K_BACK +) r1.reverse.productPlacement("pgemmax", 0.2) r2 = s.addReaction( - "r2", subs=[(ES, "front")], prds=[(E, "front"), (P, "bsoln")], rate=K_PROD) + "r2", subs=[(ES, "front")], prds=[(E, "front"), (P, "bsoln")], rate=K_PROD +) # Place molecules for initial condition inside.addMolecules(S, 500) membrane.addMolecules((E, "front"), 100) # Output and other run-time commands -s.setOutputFile('templateout.txt', True) +s.setOutputFile("templateout.txt", True) s.addCommand("molcountheader templateout.txt", "B") s.addCommand("molcount templateout.txt", "N", step=10) s.setGraphics( - "opengl_good", bg_color="white", frame_thickness=1, - text_display=["time", S, (E, "front"), (ES, "front"), P] ) + "opengl_good", + bg_color="white", + frame_thickness=1, + text_display=["time", S, (E, "front"), (ES, "front"), P], +) s = s.run(stop=10, dt=0.01) diff --git a/examples/S7_surfaces/stick2.py b/examples/S7_surfaces/stick2.py index 2b02003d..a02188fd 100644 --- a/examples/S7_surfaces/stick2.py +++ b/examples/S7_surfaces/stick2.py @@ -11,7 +11,9 @@ s = smoldyn.Simulation(low=[0, 0], high=[100, 100]) red = s.addSpecies("red", color="red", difc=dict(all=3, front=0), display_size=5) -yellow = s.addSpecies("yellow", color="black", difc=dict(soln=3, back=1), display_size=5) +yellow = s.addSpecies( + "yellow", color="black", difc=dict(soln=3, back=1), display_size=5 +) blue = s.addSpecies("blue", color="blue", difc=3, display_size=5) red.addToSolution(100) @@ -19,18 +21,18 @@ blue.addToSolution(50, pos=(20, 20)) # Construct a closed path in 2D. -p = s.addPath2D((0, 0), (100, 0), (100, 100), (0, 100), closed=True) +p = s.addPath2D([(0, 0), (100, 0), (100, 100), (0, 100)], closed=True) walls = s.addSurface("walls", panels=p.panels) -walls.setAction('both', [red, yellow, blue], "reflect") -walls.setStyle('both', color="black") +walls.setAction("both", [red, yellow, blue], "reflect") +walls.setStyle("both", color="black") sph = smoldyn.Sphere(center=(50, 50), radius=20, slices=20) surf = s.addSurface("stick", panels=[sph]) surf.setRate(red, "fsoln", "front", rate=1, revrate=0.1) surf.setRate(yellow, "bsoln", "back", rate=1, revrate=0.1) surf.setRate(blue, "fsoln", "bsoln", rate=1, revrate=0.1) -surf.setStyle('front', color=(1, 0.7, 0)) -surf.setStyle('back', color=(0.6, 0, 0.6)) +surf.setStyle("front", color=(1, 0.7, 0)) +surf.setStyle("back", color=(0.6, 0, 0.6)) surf.addMolecules((red, "front"), 100) s.setGraphics("opengl") diff --git a/examples/S99_more/maze/maze.txt b/examples/S99_more/maze/maze.txt new file mode 100644 index 00000000..598f4524 --- /dev/null +++ b/examples/S99_more/maze/maze.txt @@ -0,0 +1,134 @@ +-0.5 -0.5 -0.5 9.5 +-0.5 9.5 9 9.5 +9.5 -0.5 9.5 9.5 +0 -0.5 9.5 -0.5 +0.5 -0.5 0.5 0.5 +0.5 2.5 -0.5 2.5 +-0.5 2.5 0.5 2.5 +0.5 4.5 -0.5 4.5 +-0.5 4.5 0.5 4.5 +0.5 5.5 0.5 6.5 +0.5 7.5 -0.5 7.5 +-0.5 7.5 0.5 7.5 +0.5 7.5 0.5 8.5 +1.5 -0.5 1.5 0.5 +0.5 0.5 0.5 -0.5 +1.5 1.5 0.5 1.5 +0.5 1.5 1.5 1.5 +1.5 5.5 0.5 5.5 +1.5 6.5 0.5 6.5 +1.5 5.5 1.5 6.5 +0.5 5.5 1.5 5.5 +0.5 6.5 0.5 5.5 +0.5 6.5 1.5 6.5 +0.5 8.5 0.5 7.5 +1.5 7.5 1.5 8.5 +1.5 0.5 1.5 -0.5 +2.5 2.5 1.5 2.5 +1.5 2.5 1.5 1.5 +2.5 3.5 1.5 3.5 +2.5 2.5 2.5 3.5 +1.5 2.5 2.5 2.5 +1.5 3.5 2.5 3.5 +2.5 3.5 2.5 4.5 +1.5 6.5 1.5 5.5 +1.5 8.5 1.5 7.5 +3.5 0.5 2.5 0.5 +3.5 1.5 2.5 1.5 +3.5 0.5 3.5 1.5 +2.5 0.5 3.5 0.5 +2.5 1.5 3.5 1.5 +3.5 2.5 2.5 2.5 +2.5 3.5 2.5 2.5 +2.5 2.5 3.5 2.5 +3.5 2.5 3.5 3.5 +3.5 3.5 3.5 4.5 +2.5 4.5 2.5 3.5 +3.5 4.5 3.5 5.5 +3.5 5.5 2.5 5.5 +3.5 6.5 2.5 6.5 +2.5 5.5 3.5 5.5 +3.5 7.5 2.5 7.5 +2.5 6.5 3.5 6.5 +2.5 7.5 3.5 7.5 +3.5 8.5 3.5 9.5 +3.5 1.5 3.5 0.5 +4.5 1.5 3.5 1.5 +3.5 1.5 4.5 1.5 +4.5 3.5 3.5 3.5 +3.5 3.5 3.5 2.5 +4.5 3.5 4.5 4.5 +3.5 4.5 3.5 3.5 +3.5 3.5 4.5 3.5 +3.5 5.5 3.5 4.5 +4.5 5.5 4.5 6.5 +4.5 7.5 3.5 7.5 +3.5 7.5 4.5 7.5 +3.5 9.5 3.5 8.5 +5.5 0.5 5.5 1.5 +5.5 1.5 4.5 1.5 +4.5 1.5 5.5 1.5 +5.5 3.5 4.5 3.5 +5.5 2.5 5.5 3.5 +4.5 4.5 4.5 3.5 +4.5 3.5 5.5 3.5 +5.5 4.5 4.5 4.5 +5.5 3.5 5.5 4.5 +4.5 4.5 5.5 4.5 +4.5 6.5 4.5 5.5 +5.5 6.5 5.5 7.5 +5.5 8.5 4.5 8.5 +4.5 8.5 5.5 8.5 +6.5 0.5 5.5 0.5 +6.5 -0.5 6.5 0.5 +5.5 0.5 5.5 -0.5 +6.5 1.5 5.5 1.5 +6.5 0.5 6.5 1.5 +5.5 1.5 5.5 0.5 +5.5 0.5 6.5 0.5 +5.5 1.5 6.5 1.5 +6.5 2.5 5.5 2.5 +5.5 3.5 5.5 2.5 +5.5 2.5 6.5 2.5 +6.5 3.5 5.5 3.5 +5.5 4.5 5.5 3.5 +5.5 3.5 6.5 3.5 +6.5 5.5 5.5 5.5 +5.5 5.5 6.5 5.5 +6.5 7.5 5.5 7.5 +5.5 7.5 5.5 6.5 +5.5 7.5 6.5 7.5 +7.5 -0.5 7.5 0.5 +6.5 0.5 6.5 -0.5 +6.5 1.5 6.5 0.5 +7.5 0.5 7.5 1.5 +7.5 2.5 6.5 2.5 +7.5 3.5 6.5 3.5 +6.5 2.5 7.5 2.5 +6.5 3.5 7.5 3.5 +7.5 6.5 7.5 7.5 +7.5 7.5 7.5 8.5 +8.5 -0.5 8.5 0.5 +7.5 0.5 7.5 -0.5 +8.5 0.5 7.5 0.5 +7.5 1.5 7.5 0.5 +7.5 0.5 8.5 0.5 +8.5 2.5 8.5 3.5 +8.5 4.5 7.5 4.5 +8.5 3.5 8.5 4.5 +7.5 4.5 8.5 4.5 +7.5 7.5 7.5 6.5 +8.5 8.5 7.5 8.5 +7.5 8.5 7.5 7.5 +7.5 8.5 8.5 8.5 +8.5 8.5 8.5 9.5 +8.5 0.5 8.5 -0.5 +9.5 1.5 8.5 1.5 +8.5 1.5 9.5 1.5 +9.5 2.5 8.5 2.5 +8.5 3.5 8.5 2.5 +8.5 2.5 9.5 2.5 +8.5 4.5 8.5 3.5 +9.5 5.5 8.5 5.5 +8.5 5.5 9.5 5.5 +8.5 9.5 8.5 8.5 diff --git a/examples/S99_more/maze/maze_solve.py b/examples/S99_more/maze/maze_solve.py new file mode 100644 index 00000000..39182127 --- /dev/null +++ b/examples/S99_more/maze/maze_solve.py @@ -0,0 +1,53 @@ +# +# Solve a maze using brute force +# + +__author__ = "Dilawar Singh" +__email__ = "dilawar.s.rajput@gmail.com" + +import math +import smoldyn + +from pathlib import Path +import typing as T + +def add_maze(s, mazefile) -> T.List[smoldyn.Panel]: + """Add maze""" + panels = [] + offset, scale = 1.5, 10 + with open(mazefile, "r") as f: + for i, ll in enumerate(f.read().strip().split("\n")): + l = [scale * (offset + float(x)) for x in ll.split()] + panels += s.addPath2D([(l[0], l[1]), (l[2], l[3])]).panels + + # close the entry notch + r = s.addRectangle(corner=[10, 10], dimensions=[30], axis="+y", name="r") + panels.append(r) + return panels + + +s = smoldyn.Simulation(low=[0, 0], high=[120, 120]) + +A = s.addSpecies("A", difc=2, color="blue") +A.addToSolution(1000, pos=[15, 15]) + +B = s.addSpecies("B", difc=0, color='black', display_size=3) +B.addToSolution(10, lowpos=[105,112], highpos=[110,112]) + +panels = add_maze(s, Path(__file__).parent / "maze.txt") +maze = s.addSurface("maze", panels=panels) +maze.setAction("front", [A], "reflect") +maze.setAction("back", [A], "reflect") +maze.setStyle("both", thickness=1, color="red") + +# When any A meets B, maze is solved. So we set a reaction between A and B +s.addReaction("done", [A, B], [], rate=10) +s.addCommand("ifless B 10 text ", "E") + +# Stop when A + B -> ∅ happens i.e. maze is solved. +s.addCommand("ifless B 10 pause", "E") + +# TODO: Show the solution. + +s.addGraphics("opengl", iter=50) +s.run(3000, dt=0.01) diff --git a/source/python/smoldyn/smoldyn.py b/source/python/smoldyn/smoldyn.py index c436679a..fae93f49 100644 --- a/source/python/smoldyn/smoldyn.py +++ b/source/python/smoldyn/smoldyn.py @@ -35,8 +35,7 @@ from typing import Union, Tuple, List, Dict, Optional, Sequence - -from smoldyn.types import Color, BoundType, ColorType, DiffConst +from smoldyn.types import Color, BoundType, ColorType, DiffConst, PointType from smoldyn import _smoldyn # Path of model file. @@ -845,8 +844,32 @@ def setAction( assert k == _smoldyn.ErrorCode.ok +def _axis(p1: PointType, p2: PointType) -> str: + (x1, y1), (x2, y2) = p1, p2 + theta = math.atan2(y2 - y1, x2 - x1) + __logger__.debug(f"theta={theta} {p1} and {p2}") + if theta in [0.0, math.pi, math.pi / 2.0, -math.pi / 2.0]: + if math.isclose(theta, 0.0): + axis = "+y" + elif math.isclose(theta, math.pi): + axis = "-y" + elif math.isclose(theta, math.pi / 2.0): + axis = "-x" + elif math.isclose(theta, -math.pi / 2.0): + axis = "+x" + else: + raise RuntimeError( + "Should not be here Python3 numerical computation is broken!" + ) + + class Path2D(object): - def __init__(self, *points, simulation: _smoldyn.Simulation, closed: bool = False): + def __init__( + self, + points: List[PointType], + simulation: _smoldyn.Simulation, + closed: bool = False, + ): """Construct a 2D path from given points. A Path2D consists of `Rectangle` and `Triangle`. @@ -2268,8 +2291,8 @@ def connect(self, func, target, step: int, args: List[float] = []): """ return super().connect(func, target, step, args) - def addPath2D(self, *points, closed: bool = False): - return Path2D(*points, simulation=super(), closed=closed) + def addPath2D(self, points: List[PointType], closed: bool = False): + return Path2D(points, simulation=super(), closed=closed) def addPort( self, diff --git a/source/python/smoldyn/types.py b/source/python/smoldyn/types.py index 99834db9..8ff1b323 100644 --- a/source/python/smoldyn/types.py +++ b/source/python/smoldyn/types.py @@ -9,9 +9,14 @@ # as [1,1,0,1] specifying RGBA value. ColorType = Union[str, Tuple[float, float, float], Tuple[float, float, float, float]] +Point2 = Tuple[float, float] +Point3 = Tuple[float, float, float] + +PointType = Union[Point2, Point3] + """Diffusion coefficient can be a single numerical value (applies to all states of molecules), or it could be a dictionary of molecule state and correponding -diffusion coefficient. +diffusion coefficient. """ DiffConst = Union[float, Dict[str, float]] @@ -19,9 +24,10 @@ given, it applies to all dimentions.""" BoundType = Union[str, List[str]] + class Color: - """Color class. - """ + """Color class.""" + def __init__(self, color): assert not isinstance(color, dict) self.name = color if not isinstance(color, Color) else color.name @@ -35,4 +41,3 @@ def _toRGBA(self): def __str__(self): return str(self.name) - From a4cca184e3190c66692a561772ff317228ad571f Mon Sep 17 00:00:00 2001 From: Dilawar Singh Date: Mon, 15 Nov 2021 17:08:44 +0530 Subject: [PATCH 2/2] temp: control is NOT returned back to python Wrote to Steve. --- examples/S99_more/maze/maze_solve.py | 17 ++++++++++------- source/Smoldyn/smolsim.cpp | 10 +++++++++- source/python/CMakeLists.txt | 1 - 3 files changed, 19 insertions(+), 9 deletions(-) diff --git a/examples/S99_more/maze/maze_solve.py b/examples/S99_more/maze/maze_solve.py index 39182127..556b96b1 100644 --- a/examples/S99_more/maze/maze_solve.py +++ b/examples/S99_more/maze/maze_solve.py @@ -11,6 +11,7 @@ from pathlib import Path import typing as T + def add_maze(s, mazefile) -> T.List[smoldyn.Panel]: """Add maze""" panels = [] @@ -31,8 +32,8 @@ def add_maze(s, mazefile) -> T.List[smoldyn.Panel]: A = s.addSpecies("A", difc=2, color="blue") A.addToSolution(1000, pos=[15, 15]) -B = s.addSpecies("B", difc=0, color='black', display_size=3) -B.addToSolution(10, lowpos=[105,112], highpos=[110,112]) +B = s.addSpecies("B", difc=0, color="black", display_size=3) +B.addToSolution(10, lowpos=[105, 112], highpos=[110, 112]) panels = add_maze(s, Path(__file__).parent / "maze.txt") maze = s.addSurface("maze", panels=panels) @@ -42,12 +43,14 @@ def add_maze(s, mazefile) -> T.List[smoldyn.Panel]: # When any A meets B, maze is solved. So we set a reaction between A and B s.addReaction("done", [A, B], [], rate=10) -s.addCommand("ifless B 10 text ", "E") # Stop when A + B -> ∅ happens i.e. maze is solved. -s.addCommand("ifless B 10 pause", "E") - -# TODO: Show the solution. +s.addOutputData('posdata') +s.addCommand("molpos B posdata", "E") +s.addCommand("ifless B 10 stop", "E") s.addGraphics("opengl", iter=50) -s.run(3000, dt=0.01) +s.run(3000, dt=0.01, quit_at_end=True) + +print(s.getOutputData("posdata")) +print('all done') diff --git a/source/Smoldyn/smolsim.cpp b/source/Smoldyn/smolsim.cpp index 0c59a857..b8139ddc 100644 --- a/source/Smoldyn/smolsim.cpp +++ b/source/Smoldyn/smolsim.cpp @@ -2533,8 +2533,16 @@ int simulatetimestep(simptr sim) { sim->time+=sim->dt; // --- end of time step --- simsetvariable(sim,"time",sim->time); er=simdocommands(sim); - if(er) return er; +#ifdef COMPILE_AS_PY_MODULE + if(er == 7) { + fprintf(stdout, "command stop"); + sim->tmax = sim->time; + sim->tbreak = sim->time; + er = 0; + } +#endif + if(er) return er; if(sim->time>=sim->tmax) return 1; if(sim->time>=sim->tbreak) return 10; diff --git a/source/python/CMakeLists.txt b/source/python/CMakeLists.txt index dab7daaa..a2ad6d6e 100644 --- a/source/python/CMakeLists.txt +++ b/source/python/CMakeLists.txt @@ -11,7 +11,6 @@ target_compile_options(_pysmoldyn PRIVATE -DENABLE_PYTHON_CALLBACK) # especially where scanf is used for user input. target_compile_options(_pysmoldyn PRIVATE -DCOMPILE_AS_PY_MODULE) - target_include_directories(_pysmoldyn PRIVATE ${Python3_INCLUDE_DIRS} ${PYBIND11_SOURCE_DIR}/include) # This is suggested by pybind11: reduces binary size and disables some warnings.