diff --git a/prism/benchmark.org b/prism/benchmark.org
index 37958ba0c6cdf61bc1cedfc8db0512230d5fecd4..7a49140b1984e3ad62b842a1ebd59310b1c2fa71 100644
--- a/prism/benchmark.org
+++ b/prism/benchmark.org
@@ -87,7 +87,6 @@
    models = []
    for t in models_table:
        models.append({'name': t[0], 'path': t[1], 'type': t[2], 'constants': t[3]})
-   models
    #+END_SRC
 
 ** Tooling
@@ -98,6 +97,7 @@
    import os
    import subprocess
    import json
+   import re
    #+END_SRC
 
    The following environment definitions probably need to be adapted for
@@ -190,7 +190,7 @@
        path_with_consts = const_base_path(path, const_assignments)
        coalgebra_file = path_with_consts + ".coalgebra"
        out = subprocess.run([config['ma_cmd'], 'refine', '--stats-json', coalgebra_file],
-                            capture_output=True)
+                            stdout=subprocess.DEVNULL, stderr=subprocess.PIPE)
        print(out.stderr.decode('utf-8'))
        return json.loads(out.stderr.decode('utf-8'))
 
@@ -198,27 +198,87 @@
        export_model(model, const_assignments)
        convert_model(model, const_assignments, initial_partition)
        return partition_model(model, const_assignments)
+   #+END_SRC
+
+** Valmari's implementation
+
+   We also compare the running time of our implementation against a direct
+   implementation of markov chain lumping that also supports MDPs. This is
+   expected to be a lot faster since it's written in hand optimized C++ and
+   specialized for its use case instead of being generic like ours.
 
-   def run_benchmark(model, const_assignments, initial_partition):
-       stats = run_all(model, const_assignments, initial_partition)
-       return [model, const_assignments, initial_partition,
-               stats['states'],
-               stats['edges'],
-               stats['explicit-states'],
-               stats['initial-partition-size'],
-               stats['final-partition-size'],
-               stats['overall-duration'],
-               stats['parse-duration'],
-               stats['refine-duration']]
+   First, we need to define where the executable for this implementation is:
+
+   #+BEGIN_SRC python :session :eval never-export :results output silent
+   config['valmari_cmd'] = '../../valmari_cpp/mdpmin'
    #+END_SRC
 
-   #+NAME: bench
-   #+BEGIN_SRC python :session :eval no-export :var table="" :colnames no :hlines yes
-   [['Model', 'Consts', 'Partition', "States", "Edges", "Sort 0", "I", "Q", 't(s)', 't_p(s)', 't_r(s)'], None] + \
-       list(map(lambda t: run_benchmark(t[0], t[1], t[2]), table[2:]))
+   And then some python to generate input files and run it.
+
+   #+BEGIN_SRC python :session :eval never-export :results output silent
+   def convert_model_valmari(name, const_assignments, initial_partition):
+       model = find_model(name)
+       path, const_assignments = replace_consts_in_path(model['path'], const_assignments)
+       path_with_consts = const_base_path(path, const_assignments)
+       tra_file = "%s.tra" % path_with_consts
+       sta_file = "%s.sta" % path_with_consts
+       output_file = path_with_consts + ".valmari"
+       if os.path.isfile(output_file):
+           return
+       if initial_partition:
+           subprocess.run([config['prism_converter_cmd'],
+                           '--model-type', model['type'],
+                           '--output-format', 'valmari',
+                           '--states-file', sta_file,
+                           "--partition-on-variables", initial_partition,
+                           tra_file, output_file])
+       else:
+           subprocess.run([config['prism_converter_cmd'],
+                           '--model-type', model['type'],
+                           '--output-format', 'valmari',
+                           '--states-file', sta_file,
+                           tra_file, output_file])
+
+   def run_valmari(name, const_assignments):
+       model = find_model(name)
+       path, const_assignments = replace_consts_in_path(model['path'], const_assignments)
+       path_with_consts = const_base_path(path, const_assignments)
+       valmari_file = open(path_with_consts + ".valmari")
+       out = subprocess.run(['perf', 'stat', config['valmari_cmd']],
+                            stdin=valmari_file, stdout=subprocess.DEVNULL,
+                            stderr=subprocess.PIPE, env={'LC_ALL': 'C'})
+       m = re.search(r'(\d+.\d+)\W+seconds time elapsed', out.stderr.decode('utf-8'))
+       return m[1]
+
+   def run_all_valmari(name, const_assignments, initial_partition):
+       convert_model_valmari(name, const_assignments, initial_partition)
+       return run_valmari(name, const_assignments)
    #+END_SRC
 
 * Benchmarks
+
+  #+BEGIN_SRC python :session :eval never-export :results output silent :exports none
+  def run_benchmark(model, const_assignments, initial_partition):
+      stats = run_all(model, const_assignments, initial_partition)
+      valmari = run_all_valmari(model, const_assignments, initial_partition)
+      return [model, const_assignments, initial_partition,
+              stats['states'],
+              stats['edges'],
+              stats['explicit-states'],
+              stats['initial-partition-size'],
+              stats['final-partition-size'],
+              stats['overall-duration'],
+              stats['parse-duration'],
+              stats['refine-duration'],
+              valmari]
+  #+END_SRC
+
+  #+NAME: bench
+  #+BEGIN_SRC python :session :eval no-export :var table="" :colnames no :hlines yes
+  [['Model', 'Consts', 'Partition', "States", "Edges", "Sort 0", "I", "Q", 't(s)', 't_p(s)', 't_r(s)', 't_valmari(s)'], None] + \
+      list(map(lambda t: run_benchmark(t[0], t[1], t[2]), table[2:]))
+  #+END_SRC
+
 ** Bounded Retransmission Protocol
 
    #+BEGIN_QUOTE