diff --git a/jastadd-mquat-benchmark/results/.gitignore b/jastadd-mquat-benchmark/results/.gitignore
index b60117a22d9f06d568c670cc8852bfeeab2b0746..e9ad181c876940b67d218dc9a4d6a52253ca0444 100644
--- a/jastadd-mquat-benchmark/results/.gitignore
+++ b/jastadd-mquat-benchmark/results/.gitignore
@@ -1,6 +1,7 @@
 *
 !.gitignore
 !jastadd-mquat-plots.ipynb
+!plot.py
 !vis_scenarios.ipynb
 !to-html.sh
 !fr
diff --git a/jastadd-mquat-benchmark/results/plot.py b/jastadd-mquat-benchmark/results/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..949c39afafaab3e576bb7112525b6dd365c361a6
--- /dev/null
+++ b/jastadd-mquat-benchmark/results/plot.py
@@ -0,0 +1,109 @@
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from matplotlib import colors as mcolors
+
+# constants
+bar_width = 0.2
+solver_names = ['ilp-direct', 'simple']
+patterns = ['❌', '✔']
+
+colors = [c[4:] for c in (sorted(mcolors.TABLEAU_COLORS.keys())) if 'dark' + c[4:] in mcolors.CSS4_COLORS]
+colors.reverse()
+
+
+def to_label(name):
+    if '-' in name:
+        tokens = name.split('-')
+        return '{} ({})'.format(tokens[0].upper(), tokens[1].title())
+    else:
+        return name.title()
+
+
+def load(name, show_head=False):
+    data = pd.read_csv(name)
+    data['tool'] = data.name.astype(str)
+    data['Gen'] = (data.Gen.astype(int) > 0) * data.Gen.astype(int)
+    data['total'] = data.Gen.astype(int) + data.Solved.astype(int)
+    if show_head:
+        data.head()
+    return data
+
+
+def load_scenario_file(name, show_head=False):
+    for filename in os.listdir("scenarios"):
+        if filename.endswith(name):
+            return load("scenarios/" + filename, show_head)
+
+
+def plot_range(names, labels, ax):
+    # constants
+
+    # colors = plt.get_cmap("rainbow")
+    #     colors = list(mcolors.to_rgba_array(color))
+    N = len(solver_names)
+    width = 0.05
+    ind = np.arange(N)
+
+    barpos = []
+    for i in range(len(names)):
+        barpos += list((ind + i * (N + 0.5)) * width)
+
+    labelpos = []
+    for i in range(len(names)):
+        labelpos += [(0.5 + i * (N + 0.5)) * width]
+
+    for (name, i) in zip(names, range(len(names))):
+
+        # load data
+        data = load_scenario_file(name)
+        means = data.total
+
+        rect = ax.bar((ind + i * (N + 0.5)) * width, means, width=width, color=colors)
+
+        for r, isValid, color in zip(rect, data.Valid, colors):
+            ax.text(r.get_x() + r.get_width() / 2, max(5, r.get_y() + r.get_height()),
+                    patterns[isValid],
+                    fontname='symbola', fontsize=16, ha='center', va='bottom', color='black')
+            ax.text(r.get_x() + r.get_width() / 2, max(5, r.get_y() + r.get_height()),
+                    '%d' % int(r.get_height()),
+                    fontsize=10, ha='center', va='bottom', color='black')
+
+
+
+    ax.set_xticks(labelpos)
+    ax.set_xticklabels(labels)
+    # ax.legend(rect, [to_label(n) for n in solver_names])
+    ax.set_yticks([1, 1000, 1000000])
+
+
+def create_grid():
+    variants = (2, 4)
+    requests = (1, 2, 3)
+    depth = (1, 2, 3)
+    resources = [15, 30]
+
+
+    fig, axs = plt.subplots(len(variants), len(requests), figsize=(12, 12), subplot_kw=dict(yscale="log"))
+
+    for v, i in zip(variants, range(len(variants))):
+        for q, j in zip(requests, range(len(requests))):
+            plot_range(['size_v%d_q%d_d%d_r15.csv' % (v, q, d) for d in depth], ['depth = %d' % d for d in depth],
+                       axs[i, j])
+
+    for i in range(len(variants)):
+        axs[i][-1].set_ylabel('Implementation\nVariants = ' + str(variants[i]), rotation=0)
+        axs[i][-1].yaxis.set_label_coords(1.2, 0.5)
+
+    for i in range(len(requests)):
+        axs[0][i].set_xlabel('Requests = ' + str(requests[i]), rotation=0)
+        axs[0][i].xaxis.set_label_coords(0.5, 1.05)
+
+    fig.patch.set_facecolor('w')
+
+    plt.show()
+
+
+create_grid()