diff --git a/examples/chess-recognition.py b/examples/chess-recognition.py
index 1c21f9249..39fbb4a78 100644
--- a/examples/chess-recognition.py
+++ b/examples/chess-recognition.py
@@ -119,17 +119,11 @@ def softmax(vector):
     (Count(board[7], p) == 0),
     (Count(board[7], P) == 0),
     # The number of pieces of each piece type without promotion
-    ((Count(board, p) == 8) | (Count(board, P) == 8)).implies(
-        (Count(board, b) <= 2) &
-        (Count(board, r) <= 2) &
-        (Count(board, n) <= 2) &
-        (Count(board, q) <= 1)
+    (
+        (Count(board, b) + Count(board, r) + Count(board, n) + Count(board, q) <= 2+2+1+1 + 8-Count(board, p))
     ),
-    ((Count(board, P) == 8) | (Count(board, p) == 8)).implies(
-        (Count(board, B) <= 2) &
-        (Count(board, R) <= 2) &
-        (Count(board, N) <= 2) &
-        (Count(board, Q) <= 1)
+    (
+        (Count(board, B) + Count(board, R) + Count(board, N) + Count(board, Q) <= 2+2+1+1 + 8-Count(board, P))
     ),
     # Bishops can't have moved if the pawns are still in starting position
     ((board[1, 1] == p) & (board[1, 3] == p) & (
diff --git a/examples/sudoku_chaos_killer.py b/examples/sudoku_chaos_killer.py
new file mode 100644
index 000000000..36079c742
--- /dev/null
+++ b/examples/sudoku_chaos_killer.py
@@ -0,0 +1,119 @@
+import cpmpy as cp
+import numpy as np
+import time
+
+# This CPMpy example solves a sudoku by KNT, which can be found on https://logic-masters.de/Raetselportal/Raetsel/zeigen.php?id=0009KE
+
+SIZE = 9
+
+# sudoku cells
+cell_values = cp.intvar(1,SIZE, shape=(SIZE,SIZE))
+# decision variables for the regions
+cell_regions = cp.intvar(1,SIZE, shape=(SIZE,SIZE))
+# regions cardinals
+region_cardinals = cp.intvar(1,SIZE, shape=(SIZE,SIZE))
+
+
+def killer_cage(idxs, values, regions, total):
+    # the sum of the cells in the cage must equal the total
+    # all cells in the cage must be in the same region
+    constraints = []
+    constraints.append(cp.sum([values[r, c] for r,c in idxs]) == total)
+    constraints.append(cp.AllEqual([regions[r, c] for r,c in idxs]))
+    return constraints
+
+def get_neighbours(r, c):
+    # a cell must be orthogonally adjacent to a cell in the same region
+
+    # check if on top row
+    if r == 0:
+        if c == 0:
+            return [(r, c+1), (r+1, c)]
+        elif c == SIZE-1:
+            return [(r, c-1), (r+1, c)]
+        else:
+            return [(r, c-1), (r, c+1), (r+1, c)]
+    # check if on bottom row
+    elif r == SIZE-1:
+        if c == 0:
+            return [(r, c+1), (r-1, c)]
+        elif c == SIZE-1:
+            return [(r, c-1), (r-1, c)]
+        else:
+            return [(r, c-1), (r, c+1), (r-1, c)]
+    # check if on left column
+    elif c == 0:
+        return [(r-1, c), (r+1, c), (r, c+1)]
+    # check if on right column
+    elif c == SIZE-1:
+        return [(r-1, c), (r+1, c), (r, c-1)]
+    # check if in the middle
+    else:
+        return [(r-1, c), (r+1, c), (r, c-1), (r, c+1)]
+
+    
+
+m = cp.Model(
+    killer_cage(np.array([[0,0],[0,1],[1,1]]), cell_values, cell_regions, 18),
+    killer_cage(np.array([[1,0],[2,0],[2,1]]), cell_values, cell_regions, 8),
+    killer_cage(np.array([[3,0],[3,1],[3,2],[2,2],[1,2]]), cell_values, cell_regions, 27),
+    killer_cage(np.array([[4,0],[4,1]]), cell_values, cell_regions, 17),
+    killer_cage(np.array([[5,0],[5,1]]), cell_values, cell_regions, 8),
+    killer_cage(np.array([[1,3],[1,4]]), cell_values, cell_regions, 6),
+    killer_cage(np.array([[0,7],[0,8]]), cell_values, cell_regions, 4),
+    killer_cage(np.array([[2,3],[3,3],[3,4],[4,3]]), cell_values, cell_regions, 12),
+    killer_cage(np.array([[1,5],[2,5],[3,5],[2,4]]), cell_values, cell_regions, 28),
+    killer_cage(np.array([[4,4],[4,5],[5,5],[4,6]]), cell_values, cell_regions, 16),
+    killer_cage(np.array([[2,8],[3,8]]), cell_values, cell_regions, 15),
+    killer_cage(np.array([[3,6],[3,7],[4,7],[4,8]]), cell_values, cell_regions, 17),
+    killer_cage(np.array([[6,2],[6,3],[7,3]]), cell_values, cell_regions, 21),
+    killer_cage(np.array([[7,2],[8,2]]), cell_values, cell_regions, 5),
+    killer_cage(np.array([[8,3],[8,4]]), cell_values, cell_regions, 15),
+    killer_cage(np.array([[6,4],[7,4]]), cell_values, cell_regions, 8),
+    killer_cage(np.array([[6,5],[7,5],[6,6]]), cell_values, cell_regions, 19),
+    killer_cage(np.array([[8,5],[8,6]]), cell_values, cell_regions, 11),
+    killer_cage(np.array([[7,6],[7,7],[8,7]]), cell_values, cell_regions, 10),
+)
+
+for i in range(cell_values.shape[0]):
+    m += cp.AllDifferent(cell_values[i,:])
+    m += cp.AllDifferent(cell_values[:,i])
+
+
+total_cells = SIZE**2
+for i in range(total_cells-1):
+    r1 = i // SIZE
+    c1 = i % SIZE
+    
+    neighbours = get_neighbours(r1, c1)
+    
+    # at least one neighbour must be in the same region with a smaller cardinal number or the cell itself must have cardinal number 1; this forces connectedness of regions
+    m += cp.any([cp.all([cell_regions[r1,c1] == cell_regions[r2,c2], region_cardinals[r1, c1] > region_cardinals[r2, c2]]) for r2,c2 in neighbours]) | (region_cardinals[r1,c1] == 1)
+        
+        
+
+for r in range(1, SIZE+1):
+    # Enforce size for each region
+    m += cp.Count(cell_regions, r) == SIZE  # each region must be of size SIZE
+    
+    
+# Create unique IDs for each (value, region) pair to enforce all different
+m += cp.AllDifferent([(cell_values[i,j]-1)*SIZE+cell_regions[i,j]-1 for i in range(SIZE) for j in range(SIZE)])
+# Only 1 cell with cardinal number 1 per region
+for r in range(1, SIZE+1):
+    m += cp.Count([(region_cardinals[i,j])*(cell_regions[i,j] == r) for i in range(SIZE) for j in range(SIZE)], 1) == 1
+
+    
+# Symmetry breaking for the regions
+# fix top-left and bottom-right region to reduce symmetry
+m += cell_regions[0,0] == 1
+m += cell_regions[SIZE-1, SIZE-1] == SIZE
+
+sol = m.solve()
+
+print("The solution is:")
+print(cell_values.value())
+print("The regions are:")
+print(cell_regions.value())
+print("The region cardinals are:")
+print(region_cardinals.value())
\ No newline at end of file
diff --git a/examples/sudoku_schrodingers_rat.py b/examples/sudoku_schrodingers_rat.py
new file mode 100644
index 000000000..bfb7a5761
--- /dev/null
+++ b/examples/sudoku_schrodingers_rat.py
@@ -0,0 +1,347 @@
+import time
+import cpmpy as cp
+import numpy as np
+
+from cpmpy.tools.explain import optimal_mus, mus
+from cpmpy.tools.explain.utils import make_assump_model
+from cpmpy.transformations.get_variables import get_variables
+
+solver = "ortools"
+
+# This cpmpy example solves a sudoku by Zanno, which can be found on https://logic-masters.de/Raetselportal/Raetsel/zeigen.php?id=000N39
+
+# sudoku cells
+cells = cp.intvar(1,9, shape=(9,9), name="values")
+# path indices 0 if not on path else the var indicates at what point the rat passes this cell
+path = cp.intvar(0, 81, shape=(9,9), name="path")
+# list of cells (before, after) on path with indices and value, used for emitter constraint. (Due to the structure of the puzzle, we know there will never be more than 20 cells between two emitters on the path. This in itself is hard to know, but for efficiency reasons we reduce this amount using 'expert knowledge'.)
+sequence = cp.intvar(-1, 9, shape=(9,9,6,20), name="sequence")  # for each cell: before value, before row, before column, after value, after row, after column | -1 if not applicable
+
+# inducers for diagonal walls (prevent diagonal self-crossings)
+inducers = cp.intvar(0, 80, shape=(8,8), name="inducers")
+
+# givens
+# vertical walls
+V_WALLS = np.array([[0,0,0,0,0,0,1,0],
+                    [0,0,0,1,0,1,0,1],
+                    [0,0,0,0,1,1,1,0],
+                    [1,0,1,0,1,0,0,1],
+                    [1,0,0,0,1,0,0,0],
+                    [0,0,0,1,0,0,0,0],
+                    [0,1,0,0,0,1,0,0],
+                    [0,1,0,0,0,0,0,1],
+                    [0,0,0,0,0,0,1,0]])
+# horizontal walls
+H_WALLS = np.array([[0,1,0,0,0,0,0,0,0],
+                    [0,0,1,1,0,0,0,0,0],
+                    [0,0,1,0,0,0,0,1,0],
+                    [0,0,0,0,0,1,0,0,0],
+                    [0,0,0,0,0,0,0,1,1],
+                    [0,0,0,1,0,0,1,1,0],
+                    [0,0,0,0,0,0,0,0,0],
+                    [0,0,0,1,1,0,0,0,0],
+                    [0,0,0,0,0,0,0,0,0]])
+# corners walls (to block diagonal movement and future proofing for puzzles that only include walls just made out of a corner)
+C_WALLS = np.array([[1,1,0,1,0,1,1,1],
+                    [1,1,1,1,1,1,1,1],
+                    [1,1,1,1,1,1,1,1],
+                    [1,1,1,0,1,1,0,1],
+                    [1,0,1,1,1,0,1,1],
+                    [0,1,1,1,0,1,1,1],
+                    [0,1,0,0,0,1,0,1],
+                    [0,1,1,1,1,0,1,1]])
+
+# emitters 1 = renban, 2 = nabner, 3 = modular, 4 = entropy, 5 = region_sum, 6 = ten_sum
+EMITTERS = {(0, 5): 6, (2, 1): 5, (2, 4): 3, (2, 8): 1, (4, 2): 4, (4, 6): 6, (5, 0): 3, (5, 6): 2, (6, 5): 2, (6, 6): 4, (7, 5): 5, (8, 2): 1}
+
+def gate(idx1, idx2):
+    # the path can not pass from idx2 to idx1, and the value in idx1 must be greater than the value in idx2
+    r1, c1 = idx1
+    r2, c2 = idx2
+    return cp.all([cells[r1,c1] > cells[r2,c2], path[r1,c1] - 1 != path[r2,c2]])
+
+def get_reachable_neighbours(row, column):
+    # from a cell get the indices of its reachable neighbours
+    reachable_neighbours = []
+    if row != 0:
+        if H_WALLS[row-1, column] == 0:
+            reachable_neighbours.append([row - 1, column])
+        if column != 0:
+            if C_WALLS[row-1, column-1] == 0:
+                reachable_neighbours.append([row - 1, column - 1])
+        if column != 8:
+            if C_WALLS[row-1, column] == 0:
+                reachable_neighbours.append([row - 1, column + 1])
+    if row != 8:
+        if H_WALLS[row, column] == 0:
+            reachable_neighbours.append([row + 1, column])
+        if column != 0:
+            if C_WALLS[row, column-1] == 0:
+                reachable_neighbours.append([row + 1, column - 1])
+        if column != 8:
+            if C_WALLS[row, column] == 0:
+                reachable_neighbours.append([row + 1, column + 1])
+    if column != 0:
+        if V_WALLS[row, column-1] == 0:
+            reachable_neighbours.append([row, column - 1])
+    if column != 8:
+        if V_WALLS[row, column] == 0:
+            reachable_neighbours.append([row, column + 1])
+    return reachable_neighbours
+
+def path_valid(path):
+    constraints = []
+    for r in range(path.shape[0]):
+        for c in range(path.shape[0]):
+            neighbours = get_reachable_neighbours(r, c)
+            non_emitter_neighbours = [n for n in neighbours if tuple(n) not in EMITTERS]
+            emitter_neighbours = [n for n in neighbours if tuple(n) in EMITTERS]
+            if (r,c) == (2,8):
+                # The path starts on emitter. It doesn't have any previous cells so the first 3 vectors in the sequence must be fully -1. As for the the last 3, they are taken from the neighbour. Vice versa also applies.
+                constraints.append(cp.all([cp.all(sequence[r,c,:3].flatten() == -1)]))
+                constraints.append(cp.any([cp.all([path[nr,nc] == 2, sequence[r,c,3,0] == cells[nr,nc], sequence[r,c,4,0] == nr, sequence[r,c,5,0] == nc, cp.all(sequence[r,c,3,1:] == sequence[nr,nc,3,:19]), cp.all(sequence[r,c,4,1:] == sequence[nr,nc,4,:19]), cp.all(sequence[r,c,5,1:] == sequence[nr,nc,5,:19]), sequence[nr,nc,0,0] == cells[r,c], sequence[nr,nc,1,0] == r, sequence[nr,nc,2,0] == c, cp.all(sequence[nr,nc,0,1:] == sequence[r,c,0,:19]), cp.all(sequence[nr,nc,1,1:] == sequence[r,c,1,:19]), cp.all(sequence[nr,nc,2,1:] == sequence[r,c,2,:19])]) for nr, nc in neighbours]))
+            elif (r,c) == (6,6):
+                # nothing comes after this emitter on the path
+                constraints.append(cp.all([cp.all(sequence[r,c,3:].flatten() == -1)]))
+                constraints.append(path[r,c] == cp.max(path))
+            else:
+                # for any pathcell, the next pathcell must always be reachable
+                constraints.append((path[r,c] != 0).implies(cp.any([path[neighbour[0], neighbour[1]] == path[r,c] + 1 for neighbour in neighbours])))
+                if (r,c) not in EMITTERS:
+                    # carry over sequence values from non-emitter to non-emitter
+                    constraints.append(cp.all([((path[r,c] != 0) & (path[r,c] + 1 == path[nr,nc])).implies(cp.all([sequence[r,c,3,0] == cells[nr,nc], sequence[r,c,4,0] == nr, sequence[r,c,5,0] == nc, cp.all(sequence[r,c,3,1:] == sequence[nr,nc,3,:19]), cp.all(sequence[r,c,4,1:] == sequence[nr,nc,4,:19]), cp.all(sequence[r,c,5,1:] == sequence[nr,nc,5,:19]), sequence[nr,nc,0,0] == cells[r,c], sequence[nr,nc,1,0] == r, sequence[nr,nc,2,0] == c, cp.all(sequence[nr,nc,0,1:] == sequence[r,c,0,:19]), cp.all(sequence[nr,nc,1,1:] == sequence[r,c,1,:19]), cp.all(sequence[nr,nc,2,1:] == sequence[r,c,2,:19])])) for nr,nc in non_emitter_neighbours]))
+                    # carry over sequence values from non-emitter to emitter
+                    constraints.append(cp.all([((path[r,c] != 0) & (path[r,c] + 1 == path[nr,nc])).implies(cp.all([sequence[r,c,3,0] == cells[nr,nc], sequence[r,c,4,0] == nr, sequence[r,c,5,0] == nc, cp.all(sequence[r,c,3,1:] == -1), cp.all(sequence[r,c,4,1:] == -1), cp.all(sequence[r,c,5,1:] == -1), sequence[nr,nc,0,0] == cells[r,c], sequence[nr,nc,1,0] == r, sequence[nr,nc,2,0] == c, cp.all(sequence[nr,nc,0,1:] == sequence[r,c,0,:19]), cp.all(sequence[nr,nc,1,1:] == sequence[r,c,1,:19]), cp.all(sequence[nr,nc,2,1:] == sequence[r,c,2,:19])])) for nr,nc in emitter_neighbours]))
+                else:
+                    # carry over sequence values from emitter to non-emitter
+                    constraints.append(cp.all([((path[r,c] != 0) & (path[r,c] + 1 == path[nr,nc])).implies(cp.all([sequence[r,c,3,0] == cells[nr,nc], sequence[r,c,4,0] == nr, sequence[r,c,5,0] == nc, cp.all(sequence[r,c,3,1:] == sequence[nr,nc,3,:19]), cp.all(sequence[r,c,4,1:] == sequence[nr,nc,4,:19]), cp.all(sequence[r,c,5,1:] == sequence[nr,nc,5,:19]), sequence[nr,nc,0,0] == cells[r,c], sequence[nr,nc,1,0] == r, sequence[nr,nc,2,0] == c, cp.all(sequence[nr,nc,0,1:] == -1), cp.all(sequence[nr,nc,1,1:] == -1), cp.all(sequence[nr,nc,2,1:] == -1)])) for nr,nc in non_emitter_neighbours]))
+                    # carry over sequence values from emitter to emitter
+                    constraints.append(cp.all([((path[r,c] != 0) & (path[r,c] + 1 == path[nr,nc])).implies(cp.all([sequence[r,c,3,0] == cells[nr,nc], sequence[r,c,4,0] == nr, sequence[r,c,5,0] == nc, cp.all(sequence[r,c,3,1:] == -1), cp.all(sequence[r,c,4,1:] == -1), cp.all(sequence[r,c,5,1:] == -1), sequence[nr,nc,0,0] == cells[r,c], sequence[nr,nc,1,0] == r, sequence[nr,nc,2,0] == c, cp.all(sequence[nr,nc,0,1:] == -1), cp.all(sequence[nr,nc,1,1:] == -1), cp.all(sequence[nr,nc,2,1:] == -1)])) for nr,nc in emitter_neighbours]))
+
+            # for any non-pathcell, its sequence must be fully -1
+            constraints.append((path[r,c] == 0).implies(cp.all(sequence[r,c] == -1)))
+
+            # if the path moves diagonally, it induces a C_WALL
+            for nr, nc in neighbours:
+                if abs(nr - r) == 1 and abs(nc - c) == 1:
+                    wall_r = min(r, nr)
+                    wall_c = min(c, nc)
+                    constraints.append((cp.all([path[r,c] != 0, path[r,c] + 1 == path[nr,nc]])).implies(inducers[wall_r, wall_c] == path[r,c])) # this forces no diagonal crossings
+    return constraints
+
+def same_region(r1, c1, r2, c2):
+    # check if two cells are in the same 3x3 region
+    return cp.all([(r1 // 3 == r2 // 3), (c1 // 3 == c2 // 3)])
+
+def renban(array, rs, cs):
+    # the line must form a set of consecutive non repeating digits
+    cons = []
+    renban_emitters = [k for k, v in EMITTERS.items() if v == 1]
+
+    for i in range(len(array) - 1):
+        cons.append((cp.all([array[i] != -1, array[i+1] == -1])).implies(cp.all([cp.any([rs[i] != er, cs[i] != ec]) for (er, ec) in renban_emitters])))  # if the next cell is -1, the current cell can't be a renban emitter
+    cons.append(cp.all([cp.AllDifferentExceptN(array, -1), cp.max(array) - cp.min([a + 10*(a == -1) for a in array]) + 1 == len(array) - cp.Count(array, -1)]))
+    return cons
+
+def nabner(array, rs, cs):
+    # no two digits are consecutive and no digit repeats
+    cons = []
+    nabner_emitters = [k for k, v in EMITTERS.items() if v == 2]
+
+    for i in range(len(array) - 1):
+        cons.append((cp.all([array[i] != -1, array[i+1] == -1])).implies(cp.all([cp.any([rs[i] != er, cs[i] != ec]) for (er, ec) in nabner_emitters])))  # if the next cell is -1, the current cell can't be a nabner emitter
+
+    for i in range(len(array) - 1):
+        cons.append((array[i] != -1).implies(cp.all([cp.abs(array[i] - array[j]) > 1 for j in range(i+1,len(array))])))
+    cons.append(cp.AllDifferentExceptN(array, -1))
+    return cons
+
+def modular(array, rs, cs):
+    # every set of 3 consecutive digits must have a different value mod 3
+    cons = []
+    modular_emitters = [k for k, v in EMITTERS.items() if v == 3]
+
+    for i in range(len(array) - 1):
+        cons.append((cp.all([array[i] != -1, array[i+1] == -1])).implies(cp.all([cp.any([rs[i] != er, cs[i] != ec]) for (er, ec) in modular_emitters])))  # if the next cell is -1, the current cell can't be a modular emitter
+    arr = cp.intvar(-1, 2, shape=len(array))
+    cons.append(cp.all([cp.all(arr == (array % 3)), cp.all([cp.AllDifferentExceptN(arr[i:i+3], -1) for i in range(len(arr) - 2)])]))
+    return cons
+
+def entropy(array, rs, cs):
+    # every set of 3 consecutive digit must have 1 low digit (1-3) and 1 middle digit (4-6) and 1 high digit (7-9)
+    cons = []
+    entropy_emitters = [k for k, v in EMITTERS.items() if v == 4]
+
+    for i in range(len(array) - 1):
+        cons.append((cp.all([array[i] != -1, array[i+1] == -1])).implies(cp.all([cp.any([rs[i] != er, cs[i] != ec]) for (er, ec) in entropy_emitters])))  # if the next cell is -1, the current cell can't be an entropy emitter
+    cons.append(cp.all([cp.AllDifferentExceptN([(a+2) // 3 for a in array[i:i+3]], 0) for i in range(len(array) - 2)]))
+    return cons
+
+def region_sum(array, rs, cs, order):
+    # box borders divide the line into segments of the same sum
+    running_sums = cp.intvar(-1, 45, shape=len(array), name=f"running_sums_{rs[0]}_{cs[0]}_{order}")
+    region_sum = cp.intvar(1, 45, name=f"region_sum_{rs[0]}_{cs[0]}_{order}")
+    cons = []
+    region_sum_emitters = [k for k, v in EMITTERS.items() if v == 5]
+
+    cons.append(region_sum == running_sums[0])
+    cons.append(running_sums[-1] == array[-1])
+    for i in range(len(array)-1):
+        cons.append((cp.all([array[i+1] != -1, ~same_region(rs[i], cs[i], rs[i+1], cs[i+1])])).implies(cp.all([running_sums[i] == array[i], region_sum == running_sums[i+1]])))
+        cons.append((cp.all([(array[i+1] == -1), (array[i] != -1)])).implies(running_sums[i] == array[i]))
+        cons.append((cp.all([array[i+1] != -1, same_region(rs[i], cs[i], rs[i+1], cs[i+1])])).implies(cp.all([running_sums[i] == running_sums[i+1] + array[i]])))
+        cons.append((array[i] == -1).implies(running_sums[i] == array[i]))
+        # the last cell can not be another region_sum emitter
+        cons.append((cp.all([array[i] != -1, array[i+1] == -1])).implies(cp.all([cp.any([rs[i] != er, cs[i] != ec]) for (er, ec) in region_sum_emitters])))  # if the next cell is -1, the current cell can't be a region_sum emitter
+    cons.append(cp.Count(running_sums, region_sum) >= 2)  # at least two segments
+    return cons
+
+def ten_sum(array, rs, cs, order):
+    # the line can be divided into one or more non-overlapping segments that each sum to 10
+    running_sums = cp.intvar(-1, 45, shape=len(array), name=f"ten_running_sums_{rs[0]}_{cs[0]}_{order}")
+    splits = cp.boolvar(shape=len(array)-1, name=f"splits_{rs[0]}_{cs[0]}_{order}")
+    region_sum = 10
+    cons = []
+    ten_sum_emitters = [k for k, v in EMITTERS.items() if v == 6]
+
+    cons.append(region_sum == running_sums[0])
+    cons.append(running_sums[-1] == array[-1])
+    for i in range(len(array)-1):
+        cons.append((cp.all([array[i+1] != -1, splits[i]])).implies(cp.all([running_sums[i] == array[i], region_sum == running_sums[i+1]])))
+        cons.append((cp.all([(array[i+1] == -1), (array[i] != -1)])).implies(running_sums[i] == array[i]))
+        cons.append((cp.all([array[i+1] != -1, ~splits[i]])).implies(cp.all([running_sums[i] == running_sums[i+1] + array[i]])))
+        cons.append((cp.all([array[i] != -1, array[i+1] == -1])).implies(cp.all([cp.any([rs[i] != er, cs[i] != ec]) for (er, ec) in ten_sum_emitters])))  # if the next cell is -1, the current cell can't be a ten_sum emitter
+    return cons
+
+def activate_lines(sequence):
+    # Iterate over all emitters and activate their respective constraints based on if they are on the path
+    constraints = []
+
+    for er,ec in EMITTERS:
+        line = EMITTERS[(er,ec)]
+        before = np.concatenate(([cells[er,ec]], sequence[er,ec,0,:]))  # before emitter
+        after = np.concatenate(([cells[er,ec]], sequence[er,ec,3,:]))   # after emitter
+        before_rs = np.concatenate(([er], sequence[er,ec,1,:]))  # before emitter rows
+        before_cs = np.concatenate(([ec], sequence[er,ec,2,:]))  # before emitter columns
+        after_rs = np.concatenate(([er], sequence[er,ec,4,:]))   # after emitter rows
+        after_cs = np.concatenate(([ec], sequence[er,ec,5,:]))   # after emitter columns
+        
+        cons_before = []
+        cons_after = []
+
+        if line == 1:
+            cons_before = renban(before, before_rs, before_cs)
+            cons_after = renban(after, after_rs, after_cs)
+        elif line == 2:
+            cons_before = nabner(before, before_rs, before_cs)
+            cons_after = nabner(after, after_rs, after_cs)
+        elif line == 3:
+            cons_before = modular(before, before_rs, before_cs)
+            cons_after = modular(after, after_rs, after_cs)
+        elif line == 4:
+            cons_before = entropy(before, before_rs, before_cs)
+            cons_after = entropy(after, after_rs, after_cs)
+        elif line == 5:
+            cons_before = region_sum(before, before_rs, before_cs, "before")
+            cons_after = region_sum(after, after_rs, after_cs, "after")
+        elif line == 6:
+            cons_before = ten_sum(before, before_rs, before_cs, "before")
+            cons_after = ten_sum(after, after_rs, after_cs, "after")
+        for c in cons_before:
+            constraints.append((path[er,ec] > 1).implies(c))
+        for c in cons_after:
+            constraints.append((cp.all([0 < path[er,ec],path[er,ec] < cp.max(path)])).implies(c))
+    return constraints
+
+m = cp.Model(
+
+    # path givens
+    path[2,8] == 1,
+    cp.max(path) == path[6,6],
+
+    # no duplicate path indices
+    cp.AllDifferentExcept0(path),
+
+    # gates
+    gate((0,4), (1,4)),
+    gate((2,3), (2,2)),
+    gate((6,6), (7,6)),
+    
+    # emitter constraints
+    activate_lines(sequence),
+    # general path constraints
+    path_valid(path)
+)
+
+def regroup_to_blocks(grid):
+    # Create an empty list to store the blocks
+    blocks = [[] for _ in range(9)]
+
+    for row_index in range(9):
+        for col_index in range(9):
+            # Determine which block the current element belongs to
+            block_index = (row_index // 3) * 3 + (col_index // 3)
+            # Add the element to the appropriate block
+            blocks[block_index].append(grid[row_index][col_index])
+
+    return blocks
+
+
+blocks = regroup_to_blocks(cells)
+
+for i in range(cells.shape[0]):
+    m += cp.AllDifferent(cells[i,:])
+    m += cp.AllDifferent(cells[:,i])
+    m += cp.AllDifferent(blocks[i])
+
+def print_grid(grid):
+    for r in range(grid.shape[0]*2+1):
+        row = ""
+        if r == 0 or r == grid.shape[0]*2:
+            for c in range(grid.shape[1]*2+1):
+                if c % 2 == 0:
+                    row += "+"
+                else:
+                    row += "---"
+        elif r % 2 == 0:
+            for c in range(grid.shape[1]*2+1):
+                if c == 0 or c == grid.shape[1]*2:
+                    row += "+"
+                elif c % 2 == 0:
+                    if C_WALLS[r//2 - 1, c//2 - 1] == 1:
+                        row += "+"
+                    else:
+                        row += " "
+                else:
+                    if H_WALLS[r//2 - 1, c//2] == 1:
+                        row += "---"
+                    else:
+                        row += "   "
+        else:
+            for c in range(grid.shape[1]*2+1):
+                if c == 0 or c == grid.shape[1]*2:
+                    row += "|"
+                elif c % 2 == 0:
+                    if V_WALLS[r//2, c//2 - 1] == 1:
+                        row += "|"
+                    else:
+                        row += " "
+                else:
+                    if (r//2, c//2) not in EMITTERS:
+                        row += "   "
+                    else:
+                        row += " " + str(EMITTERS[(r//2, c//2)]) + " "
+        print(row)
+    print("")
+
+print("The puzzle is:")
+print_grid(cells)
+
+print("Number of constraints:", len(m.constraints))
+
+sol = m.solve(solver=solver)
+
+print("The solution is:")
+print(cells.value())
+print("The path is (0 if not on the path):")
+print(path.value())
\ No newline at end of file
diff --git a/examples/sudoku_zeroville.py b/examples/sudoku_zeroville.py
new file mode 100644
index 000000000..34c05379a
--- /dev/null
+++ b/examples/sudoku_zeroville.py
@@ -0,0 +1,124 @@
+import cpmpy as cp
+import numpy as np
+
+# This cpmpy example solves a sudoku by marty_sears, which can be found on https://logic-masters.de/Raetselportal/Raetsel/zeigen.php?id=000ONO
+
+# w, h of block
+w, h = 3, 2
+
+# sudoku cells
+cells = cp.intvar(0,6, shape=(9,9))
+
+# orientations of blocks
+orientations = cp.boolvar(shape=6)
+
+# starts and ends
+starts = cp.intvar(0, 7, shape=(6, 2))
+ends = cp.intvar(2, 9, shape=(6, 2))
+# block shapes
+blocks = cp.intvar(2, 3, shape=(6, 2))
+
+# intvars for cell to block mapping
+cell_blocks = cp.intvar(0, 6, shape=(9, 9))
+
+def same_difference(array):
+    # adjacent cells on the line must all have the same difference
+    diff = cp.intvar(0,6, shape=1)
+    return cp.all([abs(a-b) == diff for a,b in zip(array, array[1:])])
+
+def zipper(array):
+    # equidistant cells from the middle, sum to the value in the value in the middle
+    assert len(array) % 2 == 1
+    mid = len(array) // 2
+    return cp.all([array[i] + array[len(array)-1-i] == array[mid] for i in range(mid)])
+
+def parity(array):
+    # adjacent cells on the line must have different parity
+    return cp.all([a % 2 != b % 2 for a,b in zip(array, array[1:])])
+
+def whisper(array):
+    # adjacent cells on the line must have a difference of at least five (5)
+    return cp.all([abs(a-b) >= 5 for a,b in zip(array, array[1:])])
+
+def palindrome(array):
+    # the line must read the same forwards and backwards
+    return cp.all([array[i] == array[len(array)-1-i] for i in range(len(array)//2)])
+
+def slow_thermo(array):
+    # digits on a slow thermo must increase or stay equal
+    return cp.all([a <= b for a,b in zip(array, array[1:])])
+
+def region_sum(args):
+    # box borders divide the line into segments of the same sum
+    # use cell_blocks vars, track sum of segment until adjacent digits belong to different blocks
+    # the input in this case is actually the indices which can be used both for the cells vars and the cell_blocks vars
+
+    running_sums = cp.intvar(0, 21, shape=len(args))
+    region_sum = cp.intvar(0, 21)
+    cons = []
+    cons.append(region_sum == running_sums[0])
+    cons.append(running_sums[-1] == cells[args[-1,0], args[-1,1]])
+    for i in range(len(args)-1):
+        cons.append((cell_blocks[args[i, 0], args[i, 1]] != cell_blocks[args[i+1, 0], args[i+1, 1]]).implies(cp.all([running_sums[i] == cells[args[i, 0], args[i, 1]], region_sum == running_sums[i+1]])))
+        cons.append((cell_blocks[args[i, 0], args[i, 1]] == cell_blocks[args[i+1, 0], args[i+1, 1]]).implies(cp.all([running_sums[i] == running_sums[i+1] + cells[args[i,0], args[i,1]]])))
+
+    return cp.all(cons)
+
+
+m = cp.Model(
+    # add line constraints
+    parity(cells[0,:3]),
+    same_difference(np.concatenate((cells[0,3:5], [cells[1,5]], cells[0,6:8]))),
+    zipper(np.concatenate((cells[3,5:], cells[2::-1,8]))),
+    parity(np.array([cells[1,6], cells[2,6], cells[2,7]])),
+    whisper(np.array([cells[1,3], cells[2,2], cells[3,3], cells[2,3]])),
+    palindrome(np.array([cells[5,3], cells[6,3], cells[6,4]])),
+    palindrome(np.array([cells[5,4], cells[4,4], cells[4,5]])),
+    slow_thermo(np.concatenate((cells[8,2::-1], cells[7:4:-1,0], [cells[5,1]], [cells[4,2]], [cells[3,2]]))),
+    slow_thermo(np.array([cells[6,2], cells[7,3], cells[8,4]])),
+    zipper(np.array([cells[5,6], cells[6,7], cells[7,6]])),
+    same_difference(np.array([cells[4,8], cells[5,7], cells[6,8]])),
+    same_difference(cells[8,6:]),
+    region_sum(np.array([(1,0), (2,0), (3,0), (3,1), (2,1), (1,1)]))
+)
+
+for i in range(cells.shape[0]):
+    m += cp.AllDifferentExcept0(cells[i,:])
+    m += cp.AllDifferentExcept0(cells[:,i])
+    
+for i in range(6):
+    # block orientation
+    m += (orientations[i] == 0).implies(blocks[i,0] == w)
+    m += (orientations[i] == 1).implies(blocks[i,0] == h)
+    m += (orientations[i] == 0).implies(blocks[i,1] == h)
+    m += (orientations[i] == 1).implies(blocks[i,1] == w)
+    # block starts and ends
+    m += starts[i,0] + blocks[i,0] == ends[i,0]
+    m += starts[i,1] + blocks[i,1] == ends[i,1]
+
+# block constraints
+# Create mapping of value, block to unique ID, then all different. This also ensures there is just 6 cells for each block
+m += cp.AllDifferentExceptN([cell_blocks[i,j]*6 + cells[i,j]-1 for i in range(9) for j in range(9)], -1)
+
+# Use start and end to restrict cell_blocks values
+for i in range(9):
+    for j in range(9):
+        for b in range(6):
+            in_block = cp.boolvar()
+            in_block = cp.all([
+                cp.all([i >= starts[b,0], i < ends[b,0]]),
+                cp.all([j >= starts[b,1], j < ends[b,1]])]
+            )
+            
+            # this in also enforces the no_overlap2d constraint in a non-global way, but it is useful to actually have the mapping
+            m += (in_block).implies(cp.all([cell_blocks[i,j] == b+1, cells[i,j] != 0]))
+            m += (~in_block).implies(cell_blocks[i,j] != b+1)
+        m += (cell_blocks[i,j] == 0).implies(cells[i,j] == 0)
+
+
+# print(m)
+sol = m.solve()
+print("The solution is:")
+print(cells.value())
+print("The blocks are mapped like this:")
+print(cell_blocks.value())
\ No newline at end of file