From 96f0734dffc0c62d195353f1584c816ecf364406 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Lucas=20Fr=C3=A9rot?= <lucas.frerot@sorbonne-universite.fr>
Date: Wed, 14 Jan 2026 14:31:58 +0100
Subject: [PATCH] added plot step

---
 Snakemake-Tutorial.md | 34 ++++++++++++++++++++++++++++++++++
 simulation.py         |  8 ++++++++
 simulation_plot.py    | 11 +++++++++++
 3 files changed, 53 insertions(+)
 create mode 100644 simulation_plot.py

diff --git a/Snakemake-Tutorial.md b/Snakemake-Tutorial.md
index 84fbcd9..d11d64f 100644
--- a/Snakemake-Tutorial.md
+++ b/Snakemake-Tutorial.md
@@ -138,6 +138,8 @@ eta = float(wildcards.eta)
 dt = float(wildcards.dt)
 total_time = float(wildcards.total_time)
 
+assert eta >= 0 and dt > 0 and total_time > 0
+
 results = solve_free_fall(eta, dt, total_time)
 np.savetxt(snakemake.output[0], results)
 ```
@@ -158,12 +160,44 @@ wildcards in the output of a rule can never be smaller than the number of
 wildcards in the input (although it can be larger). This means we have to reuse
 the wildcards in the filename for the plot.
 
+Let us also say we want to change the color of the plot. Since this is not
+really a simulation parameter, we will add a `color` parameter to the rule
+itself. The advantage is that to change the parameter value one can quickly look
+at the rule, there is no need to go read the source code.
+
 ```python
 rule simulation_plot:
     input:
         data=rules.simulation.output[0],
     output:
         plot_file='time_series,eta={eta},dt={dt},total_time={total_time}.pdf'
+    params:
+        plot_color='red'
     script:
         "simulation_plot.py"
 ```
+
+Note that we have used *named* attributes for input, output, etc. This means we
+can use the names `data`, `plot_file` and `plot_color` in the script, like so:
+
+```python
+import matplotlib.pyplot as plt
+import numpy as np
+
+from snakemake.script import snakemake
+
+data = np.loadtxt(snakemake.input.data)
+
+plt.plot(data[0], data[1], color=snakemake.params.plot_color)
+plt.xlabel('$t$')
+plt.ylabel('$v(t)$')
+plt.savefig(snakemake.output.plot_file)
+```
+
+You can generate the plot with:
+
+```
+snakemake -j1 time_series,eta=1,dt=1e-1,total_time=10.pdf
+```
+
+Notice that only the plot step gets executed if the simulation data is up-to-date.
diff --git a/simulation.py b/simulation.py
index 9aedf22..cca2de4 100644
--- a/simulation.py
+++ b/simulation.py
@@ -2,6 +2,12 @@ import numpy as np
 from snakemake.script import snakemake
 
 def solve_free_fall(eta, dt, total_time):
+    """
+    Solve free fall equation with centered FD scheme.
+
+    Returns time and velocity values as Numpy arrays.
+    """
+
     Nsteps = int(total_time / dt)
     v = np.zeros(Nsteps)
     dv = 1
@@ -18,5 +24,7 @@ eta = float(wildcards.eta)
 dt = float(wildcards.dt)
 total_time = float(wildcards.total_time)
 
+assert eta >= 0 and dt > 0 and total_time > 0
+
 results = solve_free_fall(eta, dt, total_time)
 np.savetxt(snakemake.output[0], results)
diff --git a/simulation_plot.py b/simulation_plot.py
new file mode 100644
index 0000000..14cf251
--- /dev/null
+++ b/simulation_plot.py
@@ -0,0 +1,11 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+from snakemake.script import snakemake
+
+data = np.loadtxt(snakemake.input.data)
+
+plt.plot(data[0], data[1], color=snakemake.params.plot_color)
+plt.xlabel('$t$')
+plt.ylabel('$v(t)$')
+plt.savefig(snakemake.output.plot_file)