From 7dc42cb1cc61f9d8c3164fd68eb8fcdd207a72f5 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:19:11 +0100
Subject: [PATCH] added snakemake context paragraph

---
 .gitignore            |  1 +
 Snakemake-Tutorial.md | 92 ++++++++++++++++++++++++++++++++++++++++++-
 simulation.py         | 22 +++++++++++
 3 files changed, 114 insertions(+), 1 deletion(-)
 create mode 100644 simulation.py

diff --git a/.gitignore b/.gitignore
index af25493..48713c3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
 venv
 Snakefile
+/.snakemake/
diff --git a/Snakemake-Tutorial.md b/Snakemake-Tutorial.md
index 398249e..84fbcd9 100644
--- a/Snakemake-Tutorial.md
+++ b/Snakemake-Tutorial.md
@@ -13,8 +13,24 @@ numerical integration scheme and the number of time steps.
 The following function solves the ODE numerically:
 
 ```python
+import numpy as np
+
 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
+
+    for i in range(Nsteps-1):
+        vpred = v[i] + dt / 2 * dv
+        dv = (1 - eta * vpred) / (1 + eta * dt / 2)
+        v[i+1] = vpred + dt / 2 * dv
+    return np.vstack((np.arange(Nsteps) * dt, v))
 ```
 
 It returns a single array where the first line is time and the second line is velocity.
@@ -77,3 +93,77 @@ happens. Let's see how to retrieve the wildcard values and output paths from the
 script file.
 
 ### Snakemake context
+
+A running script launched with Snakemake has access to the workflow context with
+the `snakemake.script.snakemake` object, which we can import in our
+`simulation.py` script:
+
+```python
+from snakemake.script import snakemake
+```
+
+You can print the following properties of the context object:
+
+```python
+print("input", snakemake.input)
+print("output", snakemake.output)
+print("wildcards", dict(snakemake.wildcards))
+print("parameters", dict(snakemake.params))
+```
+
+We can test this with:
+
+```
+snakemake -j1 free_fall,eta=1,dt=1e-1,total_time=10.txt
+
+< ... >
+
+input
+output free_fall,eta=1,dt=1e-1,total_time=10.txt
+wildcards {'eta': '1', 'dt': '1e-1', 'total_time': '10'}
+parameters {}
+```
+
+Note that the `snakemake.wildcards` object can be interpreted as a dictionary,
+with each value associated to a wildcard being a string (note that the values
+here are not numbers). This means numerical values need to be converted before
+being passed on to our function.
+
+Now that we know how to retrieve context values, we can use them to run our
+simulation and save the result in a text file:
+
+```python
+wildcards = snakemake.wildcards
+eta = float(wildcards.eta)
+dt = float(wildcards.dt)
+total_time = float(wildcards.total_time)
+
+results = solve_free_fall(eta, dt, total_time)
+np.savetxt(snakemake.output[0], results)
+```
+
+Note that `snakemake.output` is always a list. Since we only have one output, we
+take `snakemake.output[0]`.
+
+## New rule: plotting
+
+We now want to implement the second part of the workflow: plotting the results
+of a single simulation.
+
+The rule we will implement will have an input file, which is the same as the
+output of the previous rule. Snakemake makes this easy with a `rules` object
+which is used to access context from other rules (such as output files). One
+must be careful: **wildcards** need to be conserved between rules. The number of
+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.
+
+```python
+rule simulation_plot:
+    input:
+        data=rules.simulation.output[0],
+    output:
+        plot_file='time_series,eta={eta},dt={dt},total_time={total_time}.pdf'
+    script:
+        "simulation_plot.py"
+```
diff --git a/simulation.py b/simulation.py
new file mode 100644
index 0000000..9aedf22
--- /dev/null
+++ b/simulation.py
@@ -0,0 +1,22 @@
+import numpy as np
+from snakemake.script import snakemake
+
+def solve_free_fall(eta, dt, total_time):
+    Nsteps = int(total_time / dt)
+    v = np.zeros(Nsteps)
+    dv = 1
+
+    for i in range(Nsteps-1):
+        vpred = v[i] + dt / 2 * dv
+        dv = (1 - eta * vpred) / (1 + eta * dt / 2)
+        v[i+1] = vpred + dt / 2 * dv
+    return np.vstack((np.arange(Nsteps) * dt, v))
+
+
+wildcards = snakemake.wildcards
+eta = float(wildcards.eta)
+dt = float(wildcards.dt)
+total_time = float(wildcards.total_time)
+
+results = solve_free_fall(eta, dt, total_time)
+np.savetxt(snakemake.output[0], results)