diff --git a/.github/workflows/python-test.yml b/.github/workflows/python-test.yml
index 87ae6652d..e892f30bd 100644
--- a/.github/workflows/python-test.yml
+++ b/.github/workflows/python-test.yml
@@ -23,6 +23,7 @@ jobs:
           pip install z3-solver
           pip install exact
           pip install pysdd
+          pip install pychoco
           sudo snap install minizinc --classic
           pip install minizinc
       - name: Test with pytest
diff --git a/cpmpy/exceptions.py b/cpmpy/exceptions.py
index ced4d1be0..7ed0e2143 100644
--- a/cpmpy/exceptions.py
+++ b/cpmpy/exceptions.py
@@ -17,6 +17,12 @@ class MinizincNameException(CPMpyException):
 class MinizincBoundsException(CPMpyException):
     pass
 
+class ChocoTypeException(CPMpyException):
+    pass
+
+class ChocoBoundsException(CPMpyException):
+    pass
+
 class NotSupportedError(CPMpyException):
     pass
 
diff --git a/cpmpy/solvers/__init__.py b/cpmpy/solvers/__init__.py
index 376416057..e66a71d6c 100644
--- a/cpmpy/solvers/__init__.py
+++ b/cpmpy/solvers/__init__.py
@@ -18,6 +18,7 @@
         pysdd
         z3
         exact
+        choco
         utils
 
     ===============
@@ -33,6 +34,7 @@
         CPM_pysdd
         CPM_z3
         CPM_exact
+        CPM_choco
 
     =================
     List of functions
@@ -51,3 +53,5 @@
 from .pysdd import CPM_pysdd
 from .z3 import CPM_z3
 from .exact import CPM_exact
+from .choco import CPM_choco
+
diff --git a/cpmpy/solvers/choco.py b/cpmpy/solvers/choco.py
new file mode 100644
index 000000000..f161b1360
--- /dev/null
+++ b/cpmpy/solvers/choco.py
@@ -0,0 +1,545 @@
+#!/usr/bin/env python
+## choco.py
+##
+"""
+    Interface to Choco solver's Python API
+
+    ...
+
+    Documentation of the solver's own Python API:
+    https://pypi.org/project/pychoco/
+    https://pychoco.readthedocs.io/en/latest/
+
+    ===============
+    List of classes
+    ===============
+
+    .. autosummary::
+        :nosignatures:
+
+        CPM_choco
+"""
+import time
+
+import numpy as np
+
+from ..transformations.normalize import toplevel_list
+from .solver_interface import SolverInterface, SolverStatus, ExitStatus
+from ..expressions.core import Expression, Comparison, Operator, BoolVal
+from ..expressions.globalconstraints import DirectConstraint
+from ..expressions.variables import _NumVarImpl, _IntVarImpl, _BoolVarImpl, NegBoolView, intvar
+from ..expressions.globalconstraints import GlobalConstraint
+from ..expressions.utils import is_num, is_int, is_boolexpr, is_any_list, get_bounds
+from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.get_variables import get_variables
+from ..transformations.flatten_model import flatten_constraint, flatten_objective
+from ..transformations.comparison import only_numexpr_equality
+from ..transformations.linearize import canonical_comparison
+from ..transformations.reification import only_bv_reifies, reify_rewrite
+from ..exceptions import ChocoBoundsException, ChocoTypeException, NotSupportedError
+
+
+class CPM_choco(SolverInterface):
+    """
+    Interface to the Choco solver python API
+
+    Requires that the 'pychoco' python package is installed:
+    $ pip install pychoco
+
+    See detailed installation instructions at:
+    https://pypi.org/project/pychoco/
+    https://pychoco.readthedocs.io/en/latest/
+
+    Creates the following attributes (see parent constructor for more):
+    chc_model: the pychoco.Model() created by _model()
+    chc_solver: the choco Model().get_solver() instance used in solve()
+
+    """
+
+    @staticmethod
+    def supported():
+        # try to import the package
+        try:
+            import pychoco as chc
+            return True
+        except ImportError:
+            return False
+
+    def __init__(self, cpm_model=None, subsolver=None):
+        """
+        Constructor of the native solver object
+
+        Requires a CPMpy model as input, and will create the corresponding
+        choco model and solver object (chc_model and chc_solver)
+
+        chc_model and chc_solver can both be modified externally before
+        calling solve(), a prime way to use more advanced solver features
+
+        Arguments:
+        - cpm_model: Model(), a CPMpy Model() (optional)
+        - subsolver: None
+        """
+        if not self.supported():
+            raise Exception("Install the python 'pychoco' package to use this solver interface")
+
+        import pychoco as chc
+
+        assert (subsolver is None), "Choco does not support any subsolver"
+
+        # initialise the native solver objects
+        self.chc_model = chc.Model()
+
+        # for the objective
+        self.obj = None
+        self.minimize_obj = None
+        self.helper_var = None
+        # for solving with assumption variables, TO-CHECK
+
+        # initialise everything else and post the constraints/objective
+        super().__init__(name="choco", cpm_model=cpm_model)
+
+    def solve(self, time_limit=None, **kwargs):
+        """
+            Call the Choco solver
+
+            Arguments:
+            - time_limit:  maximum solve time in seconds (float, optional)
+            - kwargs:      any keyword argument, sets parameters of solver object
+
+        """
+
+        # call the solver, with parameters
+        self.chc_solver = self.chc_model.get_solver()
+
+        start = time.time()
+
+        if time_limit is not None:
+            self.chc_solver.limit_time(str(time_limit) + "s")
+
+        if self.has_objective():
+            sol = self.chc_solver.find_optimal_solution(maximize= not self.minimize_obj,
+                                                        objective=self.solver_var(self.obj),
+                                                        **kwargs)
+        else:
+            sol = self.chc_solver.find_solution()
+        end = time.time()
+
+        # new status, get runtime
+        self.cpm_status = SolverStatus(self.name)
+        self.cpm_status.runtime = end - start
+
+        # translate exit status
+        if sol is not None:
+            if time_limit is None or self.cpm_status.runtime < time_limit: # solved to optimality
+                self.cpm_status.exitstatus = ExitStatus.OPTIMAL
+            else: # solved, but optimality not proven
+                self.cpm_status.exitstatus = ExitStatus.FEASIBLE
+        elif time_limit is None or self.cpm_status.runtime < time_limit: # proven unsat
+            self.cpm_status.exitstatus = ExitStatus.UNSATISFIABLE
+        else:
+            self.cpm_status.exitstatus = ExitStatus.UNKNOWN  # can happen when timeout is reached...
+
+
+        # True/False depending on self.chc_status
+        has_sol = sol is not None
+
+        # translate solution values (of user specified variables only)
+        self.objective_value_ = None
+        if has_sol:
+            # fill in variable values
+            for cpm_var in self.user_vars:
+                value = sol.get_int_val(self.solver_var(cpm_var))
+                if isinstance(cpm_var, _BoolVarImpl):
+                    cpm_var._value = bool(value)
+                else:
+                    cpm_var._value = value
+
+            # translate objective
+            if self.has_objective():
+                self.objective_value_ = sol.get_int_val(self.solver_var(self.obj))
+
+        return has_sol
+
+    def solveAll(self, display=None, time_limit=None, solution_limit=None, call_from_model=False, **kwargs):
+        """
+            Compute all (optimal) solutions, map them to CPMpy and optionally display the solutions.
+
+            Arguments:
+                - display: either a list of CPMpy expressions, OR a callback function, called with the variables after value-mapping
+                        default/None: nothing displayed
+                - solution_limit: stop after this many solutions (default: None)
+                - time_limit:  maximum solve time in seconds (float, default: None)
+
+            Returns: number of solutions found
+        """
+
+        if time_limit is not None:
+            self.chc_solver.limit_time(str(time_limit) + "s")
+
+        self.chc_solver = self.chc_model.get_solver()
+        start = time.time()
+        if self.has_objective():
+            sols = self.chc_solver.find_all_optimal_solutions(maximize=not self.minimize_obj,
+                                                              solution_limit=solution_limit,
+                                                              objective=self.solver_var(self.obj),
+                                                              **kwargs)
+        else:
+            sols = self.chc_solver.find_all_solutions(solution_limit=solution_limit)
+        end = time.time()
+
+        # new status, get runtime
+        self.cpm_status = SolverStatus(self.name)
+        self.cpm_status.runtime = end - start
+
+        # display if needed
+        if display is not None:
+            for sol in sols:
+                # map the solution to user vars
+                for cpm_var in self.user_vars:
+                    value = sol.get_int_val(self.solver_var(cpm_var))
+                    if isinstance(cpm_var, _BoolVarImpl):
+                        cpm_var._value = bool(value)
+                    else:
+                        cpm_var._value = value
+                # print the desired display
+                if isinstance(display, Expression):
+                    print(display.value())
+                elif isinstance(display, list):
+                    print([v.value() for v in display])
+                else:
+                    display()  # callback
+
+        return len(sols)
+
+    def solver_var(self, cpm_var):
+        """
+            Creates solver variable for cpmpy variable
+            or returns from cache if previously created
+        """
+        if is_num(cpm_var):  # shortcut, eases posting constraints
+            if not is_int(cpm_var):
+                raise ValueError(f"Choco only accepts integer constants, got {cpm_var} of type {type(cpm_var)}")
+            if cpm_var < -2147483646 or cpm_var > 2147483646:
+                raise ChocoBoundsException(
+                    "Choco does not accept integer literals with bounds outside of range (-2147483646..2147483646)")
+            return int(cpm_var)
+
+        # special case, negative-bool-view
+        # work directly on var inside the view
+        if isinstance(cpm_var, NegBoolView):
+            return self.chc_model.bool_not_view(self.solver_var(cpm_var._bv))
+
+        # create if it does not exist
+        if cpm_var not in self._varmap:
+            if isinstance(cpm_var, _BoolVarImpl):
+                revar = self.chc_model.boolvar(name=str(cpm_var.name))
+            elif isinstance(cpm_var, _IntVarImpl):
+                if cpm_var.lb < -2147483646 or cpm_var.ub > 2147483646:
+                    raise ChocoBoundsException(
+                        "Choco does not accept variables with bounds outside of range (-2147483646..2147483646)")
+                revar = self.chc_model.intvar(cpm_var.lb, cpm_var.ub, name=str(cpm_var.name))
+            else:
+                raise NotImplementedError("Not a known var {}".format(cpm_var))
+            self._varmap[cpm_var] = revar
+
+        return self._varmap[cpm_var]
+
+    def objective(self, expr, minimize):
+        """
+            Post the given expression to the solver as objective to minimize/maximize
+
+            - expr: Expression, the CPMpy expression that represents the objective function
+            - minimize: Bool, whether it is a minimization problem (True) or maximization problem (False)
+
+            'objective()' can be called multiple times, only the last one is stored
+
+            (technical side note: constraints created during conversion of the objective
+            are premanently posted to the solver. Choco accepts variables to maximize or minimize
+            so it is needed to post constraints and create auxiliary variables)
+        """
+
+        # make objective function non-nested
+        obj_var = intvar(*get_bounds(expr))
+        self += obj_var == expr
+
+        self.obj = obj_var
+        self.minimize_obj = minimize  # Choco has as default to maximize
+
+    def has_objective(self):
+        return self.obj is not None
+
+
+    def _to_var(self, val):
+        from pychoco.variables.intvar import IntVar
+        if is_int(val):
+            # Choco accepts only int32, not int64
+            if val < -2147483646 or val > 2147483646:
+                raise ChocoBoundsException(
+                    "Choco does not accept integer literals with bounds outside of range (-2147483646..2147483646)")
+            return self.chc_model.intvar(int(val), int(val))  # convert to "variable"
+        elif isinstance(val, _NumVarImpl):
+            return self.solver_var(val)  # use variable
+        else:
+            raise ValueError(f"Cannot convert {val} of type {type(val)} to Choco variable, expected int or NumVarImpl")
+
+        # elif isinstance(val, IntVar):
+        #     return val
+        # return None
+
+    def _to_vars(self, vals):
+        if is_any_list(vals):
+            return [self._to_vars(v) for v in vals]
+        return self._to_var(vals)
+
+    def transform(self, cpm_expr):
+        """
+            Transform arbitrary CPMpy expressions to constraints the solver supports
+
+            Implemented through chaining multiple solver-independent **transformation functions** from
+            the `cpmpy/transformations/` directory.
+
+            See the 'Adding a new solver' docs on readthedocs for more information.
+
+        :param cpm_expr: CPMpy expression, or list thereof
+        :type cpm_expr: Expression or list of Expression
+
+        :return: list of Expression
+        """
+
+        cpm_cons = toplevel_list(cpm_expr)
+        supported = {"min", "max", "abs", "count", "element", "alldifferent", "alldifferent_except0", "allequal",
+                     "table", "InDomain", "cumulative", "circuit", "gcc", "inverse", "nvalue"}
+        supported_reified = supported # choco supports reification of any constraint
+        cpm_cons = decompose_in_tree(cpm_cons, supported, supported_reified)
+        cpm_cons = flatten_constraint(cpm_cons)  # flat normal form
+        cpm_cons = canonical_comparison(cpm_cons)
+        cpm_cons = reify_rewrite(cpm_cons, supported = supported_reified | {"sum", "wsum"})  # constraints that support reification
+        cpm_cons = only_numexpr_equality(cpm_cons, supported=frozenset(["sum", "wsum", "sub"]))  # support >, <, !=
+
+        return cpm_cons
+
+    def __add__(self, cpm_expr):
+        """
+            Eagerly add a constraint to the underlying solver.
+
+            Any CPMpy expression given is immediately transformed (through `transform()`)
+            and then posted to the solver in this function.
+
+            This can raise 'NotImplementedError' for any constraint not supported after transformation
+
+            The variables used in expressions given to add are stored as 'user variables'. Those are the only ones
+            the user knows and cares about (and will be populated with a value after solve). All other variables
+            are auxiliary variables created by transformations.
+
+        :param cpm_expr: CPMpy expression, or list thereof
+        :type cpm_expr: Expression or list of Expression
+
+        :return: self
+        """
+        # add new user vars to the set
+        get_variables(cpm_expr, collect=self.user_vars)
+
+        # transform and post the constraints
+        for con in self.transform(cpm_expr):
+            c = self._get_constraint(con)
+            if c is not None: # Reification constraints are not posted
+                c.post()
+
+        return self
+
+    def _get_constraint(self, cpm_expr):
+        """
+        Get a solver's constraint by a supported CPMpy constraint
+
+        :param cpm_expr: CPMpy expression
+        :type cpm_expr: Expression
+
+        """
+
+        # Operators: base (bool), lhs=numexpr, lhs|rhs=boolexpr (reified ->)
+        if isinstance(cpm_expr, Operator):
+            # 'and'/n, 'or'/n, '->'/2
+            if cpm_expr.name == 'and':
+                return self.chc_model.and_(self.solver_vars(cpm_expr.args))
+            elif cpm_expr.name == 'or':
+                return self.chc_model.or_(self.solver_vars(cpm_expr.args))
+
+            # elif cpm_expr.name == "->": # prepared for if pychoco releases if_then addition
+            #     cond, subexpr = cpm_expr.args
+            #     if isinstance(cond, _BoolVarImpl) and isinstance(subexpr, _BoolVarImpl):
+            #         return self.chc_model.or_(self.solver_vars([~cond, subexpr]))
+            #     elif isinstance(cond, _BoolVarImpl):
+            #         return self._get_constraint(subexpr).implied_by(self.solver_var(cond))
+            #     elif isinstance(subexpr, _BoolVarImpl):
+            #         return self._get_constraint(cond).implies(self.solver_var(subexpr))
+            #     else:
+            #         ValueError(f"Unexpected implication: {cpm_expr}")
+
+            elif cpm_expr.name == '->':
+                cond, subexpr = cpm_expr.args
+                if isinstance(cond, _BoolVarImpl) and isinstance(subexpr, _BoolVarImpl): # bv -> bv
+                    chc_cond, chc_subexpr = self.solver_vars([cond, subexpr])
+                elif isinstance(cond, _BoolVarImpl): # bv -> expr
+                    chc_cond = self.solver_var(cond)
+                    chc_subexpr = self._get_constraint(subexpr).reify()
+                elif isinstance(subexpr, _BoolVarImpl): # expr -> bv
+                    chc_cond = self._get_constraint(cond).reify()
+                    chc_subexpr = self.solver_var(subexpr)
+                else:
+                    raise ValueError(f"Unexpected implication {cpm_expr}")
+
+                return self.chc_model.or_([self.chc_model.bool_not_view(chc_cond), chc_subexpr])
+
+            else:
+                raise NotImplementedError("Not a known supported Choco Operator '{}' {}".format(
+                    cpm_expr.name, cpm_expr))
+
+        # Comparisons: both numeric and boolean ones
+        # numexpr `comp` bvar|const
+        elif isinstance(cpm_expr, Comparison):
+            lhs, rhs = cpm_expr.args
+            op = cpm_expr.name if cpm_expr.name != "==" else "="
+
+            if is_boolexpr(lhs) and is_boolexpr(rhs): #boolean equality -- Reification
+                # # prepared for if pychoco releases reify_with addition
+                # if isinstance(lhs, _BoolVarImpl) and isinstance(lhs, _BoolVarImpl):
+                #     return self.chc_model.all_equal(self.solver_vars([lhs, rhs]))
+                # elif isinstance(lhs, _BoolVarImpl):
+                #     return self._get_constraint(rhs).reify_with(self.solver_var(lhs))
+                # elif isinstance(rhs, _BoolVarImpl):
+                #     return self._get_constraint(lhs).reify_with(self.solver_var(rhs))
+                # else:
+                #     raise ValueError(f"Unexpected reification {cpm_expr}")
+
+                if isinstance(lhs, _BoolVarImpl) and isinstance(lhs, _BoolVarImpl):
+                    chc_var, bv = self.solver_vars([lhs, rhs])
+                elif isinstance(lhs, _BoolVarImpl):
+                    bv = self._get_constraint(rhs).reify()
+                    chc_var = self.solver_var(lhs)
+                elif isinstance(rhs, _BoolVarImpl):
+                    bv = self._get_constraint(lhs).reify()
+                    chc_var = self.solver_var(rhs)
+                else:
+                    raise ValueError(f"Unexpected reification {cpm_expr}")
+                return self.chc_model.all_equal([chc_var, bv])
+
+            elif isinstance(lhs, _NumVarImpl):
+                return self.chc_model.arithm(self.solver_var(lhs), op, self.solver_var(rhs))
+            elif isinstance(lhs, Operator) and lhs.name in {'sum','wsum','sub'}:
+                if lhs.name == 'sum':
+                    return self.chc_model.sum(self.solver_vars(lhs.args), op, self.solver_var(rhs))
+                elif lhs.name == "sub":
+                    a, b = self.solver_vars(lhs.args)
+                    return self.chc_model.arithm(a, "-", b, op, self.solver_var(rhs))
+                elif lhs.name == 'wsum':
+                    wgt, x = lhs.args
+                    w = np.array(wgt).tolist()
+                    x = self.solver_vars(lhs.args[1])
+                    return self.chc_model.scalar(x, w, op, self.solver_var(rhs))
+
+            elif cpm_expr.name == '==':
+
+                chc_rhs = self._to_var(rhs) # result is always var
+                all_vars = {"min", "max", "abs", "div", "mod", "element", "nvalue"}
+                if lhs.name in all_vars:
+
+                    chc_args = self._to_vars(lhs.args)
+
+                    if lhs.name == 'min': # min(vars) = var
+                        return self.chc_model.min(chc_rhs, chc_args)
+                    elif lhs.name == 'max': # max(vars) = var
+                        return self.chc_model.max(chc_rhs, chc_args)
+                    elif lhs.name == 'abs': # abs(var) = var
+                        assert len(chc_args) == 1, f"Expected one argument of abs constraint, but got {chc_args}"
+                        return self.chc_model.absolute(chc_rhs, chc_args[0])
+                    elif lhs.name == "div": # var / var = var
+                        dividend, divisor = chc_args
+                        return self.chc_model.div(dividend, divisor, chc_rhs)
+                    elif lhs.name == 'mod': # var % var = var
+                        dividend, divisor = chc_args
+                        return self.chc_model.mod(dividend, divisor, chc_rhs)
+                    elif lhs.name == "element": # varsvar[var] = var
+                        # TODO: actually, Choco also supports ints[var] = var, but no mix of var and int in array
+                        arr, idx = chc_args
+                        return self.chc_model.element(chc_rhs, arr, idx)
+                    elif lhs.name == "nvalue": # nvalue(vars) = var
+                        # TODO: should look into leaving nvalue <= arg so can post atmost_nvalues here
+                        return self.chc_model.n_values(chc_args, chc_rhs)
+
+                elif lhs.name == 'count': # count(vars, var/int) = var
+                    arr, val = lhs.args
+                    return self.chc_model.count(self.solver_var(val), self._to_vars(arr), chc_rhs)
+                elif lhs.name == 'mul': # var * var/int = var/int
+                    a,b = self.solver_vars(lhs.args)
+                    if isinstance(a, int):
+                        a,b = b,a # int arg should always be second
+                    return self.chc_model.times(a,b, self.solver_var(rhs))
+                elif lhs.name == 'pow': # var ^ int = var
+                    chc_rhs = self._to_var(rhs)
+                    return self.chc_model.pow(*self.solver_vars(lhs.args),chc_rhs)
+
+                raise NotImplementedError(
+                    "Not a known supported Choco left-hand-side '{}' {}".format(lhs.name, cpm_expr))
+
+        # base (Boolean) global constraints
+        elif isinstance(cpm_expr, GlobalConstraint):
+
+            # many globals require all variables as arguments
+            if cpm_expr.name in {"alldifferent", "alldifferent_except0", "allequal", "circuit", "inverse"}:
+                chc_args = self._to_vars(cpm_expr.args)
+                if cpm_expr.name == 'alldifferent':
+                    return self.chc_model.all_different(chc_args)
+                elif cpm_expr.name == 'alldifferent_except0':
+                    return self.chc_model.all_different_except_0(chc_args)
+                elif cpm_expr.name == 'allequal':
+                    return self.chc_model.all_equal(chc_args)
+                elif cpm_expr.name == "circuit":
+                    return self.chc_model.circuit(chc_args)
+                elif cpm_expr.name == "inverse":
+                    return self.chc_model.inverse_channeling(*chc_args)
+
+            # but not all
+            elif cpm_expr.name == 'table':
+                assert (len(cpm_expr.args) == 2)  # args = [array, table]
+                array, table = self.solver_vars(cpm_expr.args)
+                return self.chc_model.table(array, table)
+            elif cpm_expr.name == 'InDomain':
+                assert len(cpm_expr.args) == 2  # args = [array, list of vals]
+                expr, table = self.solver_vars(cpm_expr.args)
+                return self.chc_model.member(expr, table)
+            elif cpm_expr.name == "cumulative":
+                start, dur, end, demand, cap = cpm_expr.args
+                # start, end, demand and cap should be var
+                start, end, demand, cap = self._to_vars([start, end, demand, cap])
+                # duration can be var or int
+                dur = self.solver_vars(dur)
+                # Create task variables. Choco can create them only one by one
+                tasks = [self.chc_model.task(s, d, e) for s, d, e in zip(start, dur, end)]
+                return self.chc_model.cumulative(tasks, demand, cap)
+            elif cpm_expr.name == "gcc":
+                vars, vals, occ = cpm_expr.args
+                return self.chc_model.global_cardinality(*self.solver_vars([vars, vals]), self._to_vars(occ))
+            else:
+                raise NotImplementedError(f"Unknown global constraint {cpm_expr}, should be decomposed! If you reach this, please report on github.")
+
+        # unlikely base case: Boolean variable
+        elif isinstance(cpm_expr, _BoolVarImpl):
+            return self.chc_model.and_([self.solver_var(cpm_expr)])
+
+        # unlikely base case: True or False
+        elif isinstance(cpm_expr, BoolVal):
+            # Choco does not allow to post True or False. Post "certainly True or False" constraints instead
+            if cpm_expr.args[0] is True:
+                return None
+            else:
+                if self.helper_var is None:
+                    self.helper_var = self.chc_model.intvar(0, 0)
+                return self.chc_model.arithm(self.helper_var, "<", 0)
+
+        # a direct constraint, pass to solver
+        elif isinstance(cpm_expr, DirectConstraint):
+            c = cpm_expr.callSolver(self, self.chc_model)
+            return c
+
+        # else
+        raise NotImplementedError(cpm_expr)  # if you reach this... please report on github
\ No newline at end of file
diff --git a/cpmpy/solvers/utils.py b/cpmpy/solvers/utils.py
index 2464a0cd6..f9768d782 100644
--- a/cpmpy/solvers/utils.py
+++ b/cpmpy/solvers/utils.py
@@ -28,6 +28,7 @@
 from .z3 import CPM_z3
 from .pysdd import CPM_pysdd
 from .exact import CPM_exact
+from .choco import CPM_choco
 
 def param_combinations(all_params, remaining_keys=None, cur_params=None):
     """
@@ -76,6 +77,7 @@ def base_solvers(cls):
                 ("pysat", CPM_pysat),
                 ("pysdd", CPM_pysdd),
                 ("exact", CPM_exact),
+                ("choco", CPM_choco),
                ]
 
     @classmethod
@@ -134,7 +136,7 @@ def lookup(cls, name=None):
 
 # using `builtin_solvers` is DEPRECATED, use `SolverLookup` object instead
 # Order matters! first is default, then tries second, etc...
-builtin_solvers = [CPM_ortools, CPM_gurobi, CPM_minizinc, CPM_pysat, CPM_exact]
+builtin_solvers = [CPM_ortools, CPM_gurobi, CPM_minizinc, CPM_pysat, CPM_exact, CPM_choco]
 def get_supported_solvers():
     """
         Returns a list of solvers supported on this machine.
diff --git a/cpmpy/transformations/comparison.py b/cpmpy/transformations/comparison.py
index 6f4f40b91..49e24f41d 100644
--- a/cpmpy/transformations/comparison.py
+++ b/cpmpy/transformations/comparison.py
@@ -1,8 +1,9 @@
 import copy
 
 from .flatten_model import get_or_make_var
-from ..expressions.core import Comparison
-from ..expressions.variables import _NumVarImpl
+from ..expressions.core import Comparison, Operator
+from ..expressions.utils import is_boolexpr
+from ..expressions.variables import _NumVarImpl, _BoolVarImpl
 
 """
   Transformations regarding Comparison constraints.
@@ -28,17 +29,54 @@ def only_numexpr_equality(constraints, supported=frozenset()):
     # shallow copy (could support inplace too this way...)
     newcons = copy.copy(constraints)
 
-    for i,con in enumerate(newcons):
-        if isinstance(con, Comparison) and con.name != '==':
-            # LHS <op> IV    with <op> one of !=,<,<=,>,>=
-            lhs = con.args[0]
-            if not isinstance(lhs, _NumVarImpl) and not lhs.name in supported:
-                # LHS is unsupported for LHS <op> IV, rewrite to `(LHS == A) & (A <op> IV)`
-                (lhsvar, lhscons) = get_or_make_var(lhs)
-                # replace comparison by A <op> IV
-                newcons[i] = Comparison(con.name, lhsvar, con.args[1])
-                # add lhscon(s), which will be [(LHS == A)]
-                assert(len(lhscons) == 1), "only_numexpr_eq: lhs surprisingly non-flat"
-                newcons.insert(i, lhscons[0])
+    for i,cpm_expr in enumerate(newcons):
+
+        if isinstance(cpm_expr, Operator) and cpm_expr.name == "->":
+            cond, subexpr = cpm_expr.args
+            if not isinstance(cond, _BoolVarImpl): # expr -> bv
+                res = only_numexpr_equality([cond], supported)
+                if len(res) > 1:
+                    newcons[i] = res[1].implies(subexpr)
+                    newcons.insert(i, res[0])
+
+            elif not isinstance(subexpr, _BoolVarImpl):  # expr -> bv
+                res = only_numexpr_equality([subexpr], supported)
+                if len(res) > 1:
+                    newcons[i] = cond.implies(res[1])
+                    newcons.insert(i, res[0])
+            else: #bv -> bv
+                pass
+
+
+        if isinstance(cpm_expr, Comparison):
+            lhs, rhs = cpm_expr.args
+
+            if cpm_expr.name == "==" and is_boolexpr(lhs) and is_boolexpr(rhs): # reification, check recursively
+
+                if not isinstance(lhs, _BoolVarImpl):  # expr == bv
+                    res = only_numexpr_equality([lhs], supported)
+                    if len(res) > 1:
+                        newcons[i] = res[1] == rhs
+                        newcons.insert(i, res[0])
+
+                elif not isinstance(rhs, _BoolVarImpl):  # bv == expr
+                    res = only_numexpr_equality([rhs], supported)
+                    if len(res) > 1:
+                        newcons[i] = lhs == res[1]
+                        newcons.insert(i, res[0])
+                else:  # bv == bv
+                    pass
+
+            elif cpm_expr.name != "==":
+                # LHS <op> IV    with <op> one of !=,<,<=,>,>=
+                lhs = cpm_expr.args[0]
+                if not isinstance(lhs, _NumVarImpl) and lhs.name not in supported:
+                    # LHS is unsupported for LHS <op> IV, rewrite to `(LHS == A) & (A <op> IV)`
+                    (lhsvar, lhscons) = get_or_make_var(lhs)
+                    # replace comparison by A <op> IV
+                    newcons[i] = Comparison(cpm_expr.name, lhsvar, cpm_expr.args[1])
+                    # add lhscon(s), which will be [(LHS == A)]
+                    assert(len(lhscons) == 1), "only_numexpr_eq: lhs surprisingly non-flat"
+                    newcons.insert(i, lhscons[0])
 
     return newcons
diff --git a/cpmpy/transformations/linearize.py b/cpmpy/transformations/linearize.py
index d18e7d048..e0fea4c50 100644
--- a/cpmpy/transformations/linearize.py
+++ b/cpmpy/transformations/linearize.py
@@ -269,8 +269,10 @@ def canonical_comparison(lst_of_expr):
             elif isinstance(lhs, Comparison):
                 lhs = canonical_comparison(lhs)[0]
                 newlist.append(lhs.implies(rhs))
+            else:
+                newlist.append(cpm_expr)
 
-        if isinstance(cpm_expr, Comparison):
+        elif isinstance(cpm_expr, Comparison):
             lhs, rhs = cpm_expr.args
             if isinstance(lhs, Comparison) and cpm_expr.name == "==":  # reification of comparison
                 lhs = canonical_comparison(lhs)[0]
diff --git a/tests/test_constraints.py b/tests/test_constraints.py
index 246e39ed2..40b7e1590 100644
--- a/tests/test_constraints.py
+++ b/tests/test_constraints.py
@@ -17,6 +17,7 @@
                   "pysat": {"circuit", "element","min","max","count", "nvalue", "allequal","alldifferent","cumulative"},
                   "pysdd": {"circuit", "element","min","max","count", "nvalue", "allequal","alldifferent","cumulative",'xor'},
                   "exact": {},
+                  "choco": {}
                   }
 
 # Exclude certain operators for solvers.
diff --git a/tests/test_direct.py b/tests/test_direct.py
index 63da961ed..71b94d78b 100644
--- a/tests/test_direct.py
+++ b/tests/test_direct.py
@@ -4,7 +4,7 @@
 import pytest
 
 from cpmpy import *
-from cpmpy.solvers import CPM_gurobi, CPM_pysat, CPM_minizinc, CPM_pysdd, CPM_z3, CPM_exact
+from cpmpy.solvers import CPM_gurobi, CPM_pysat, CPM_minizinc, CPM_pysdd, CPM_z3, CPM_exact, CPM_choco
 
 
 class TestDirectORTools(unittest.TestCase):
@@ -127,4 +127,17 @@ def test_direct_poly(self):
 
         self.assertEqual(y.value(), poly_val)
 
+@pytest.mark.skipif(not CPM_choco.supported(),
+                    reason="pychoco not installed")
+class TestDirectChoco(unittest.TestCase):
+
+    def test_direct_global(self):
+        iv = intvar(1,9, shape=3)
+
+        model = SolverLookup.get("choco")
+
+        model += DirectConstraint("increasing", iv)
+        model += iv[1] < iv[0]
+
+        self.assertFalse(model.solve())
 
diff --git a/tests/test_solveAll.py b/tests/test_solveAll.py
index 88b31a3e5..a8d0dd82b 100644
--- a/tests/test_solveAll.py
+++ b/tests/test_solveAll.py
@@ -29,8 +29,7 @@ def test_solveall_no_obj(self):
                 count = solver.solveAll(solution_limit=1000, display=add_sol)
             self.assertEqual(3, count)
             self.assertSetEqual(sols,
-                                {"[True, True]","[True, False]","[False, True]"})
-
+                                {"[True, True]", "[True, False]", "[False, True]"})
 
     def test_solveall_with_obj(self):
 
diff --git a/tests/test_solvers.py b/tests/test_solvers.py
index e8f06849b..0297afc6c 100644
--- a/tests/test_solvers.py
+++ b/tests/test_solvers.py
@@ -8,6 +8,8 @@
 from cpmpy.solvers.minizinc import CPM_minizinc
 from cpmpy.solvers.gurobi import CPM_gurobi
 from cpmpy.solvers.exact import CPM_exact
+from cpmpy.solvers.choco import CPM_choco
+
 
 from cpmpy.exceptions import MinizincNameException
 
@@ -526,6 +528,118 @@ def test_false(self):
             if cls.supported():
                 self.assertFalse(m.solve(solver=name))
 
+    @pytest.mark.skipif(not CPM_choco.supported(),
+                        reason="pychoco not installed")
+    def test_choco(self):
+        bv = cp.boolvar(shape=3)
+        iv = cp.intvar(0, 9, shape=3)
+        # circular 'bigger then', UNSAT
+        m = cp.Model([
+            bv[0].implies(iv[0] > iv[1]),
+            bv[1].implies(iv[1] > iv[2]),
+            bv[2].implies(iv[2] > iv[0])
+        ])
+        m += sum(bv) == len(bv)
+        s = cp.SolverLookup.get("choco", m)
+
+        self.assertFalse(s.solve())
+
+        m = cp.Model(~(iv[0] != iv[1]))
+        s = cp.SolverLookup.get("choco", m)
+        self.assertTrue(s.solve())
+
+        m = cp.Model((iv[0] == 0) & ((iv[0] != iv[1]) == 0))
+        s = cp.SolverLookup.get("choco", m)
+        self.assertTrue(s.solve())
+
+        m = cp.Model([~bv, ~((iv[0] + abs(iv[1])) == sum(iv))])
+        s = cp.SolverLookup.get("choco", m)
+        self.assertTrue(s.solve())
+
+    @pytest.mark.skipif(not CPM_choco.supported(),
+                        reason="pychoco not installed")
+    def test_choco_element(self):
+
+        # test 1-D
+        iv = cp.intvar(-8, 8, 3)
+        idx = cp.intvar(-8, 8)
+        # test directly the constraint
+        constraints = [cp.Element(iv, idx) == 8]
+        model = cp.Model(constraints)
+        s = cp.SolverLookup.get("choco", model)
+        self.assertTrue(s.solve())
+        self.assertTrue(iv.value()[idx.value()] == 8)
+        self.assertTrue(cp.Element(iv, idx).value() == 8)
+        # test through __get_item__
+        constraints = [iv[idx] == 8]
+        model = cp.Model(constraints)
+        s = cp.SolverLookup.get("choco", model)
+        self.assertTrue(s.solve())
+        self.assertTrue(iv.value()[idx.value()] == 8)
+        self.assertTrue(cp.Element(iv, idx).value() == 8)
+        # test 2-D
+        iv = cp.intvar(-8, 8, shape=(3, 3))
+        idx = cp.intvar(0, 3)
+        idx2 = cp.intvar(0, 3)
+        constraints = [iv[idx, idx2] == 8]
+        model = cp.Model(constraints)
+        s = cp.SolverLookup.get("choco", model)
+        self.assertTrue(s.solve())
+        self.assertTrue(iv.value()[idx.value(), idx2.value()] == 8)
+
+    @pytest.mark.skipif(not CPM_choco.supported(),
+                        reason="pychoco not installed")
+    def test_choco_gcc_alldiff(self):
+
+        iv = cp.intvar(-8, 8, shape=5)
+        occ = cp.intvar(0, len(iv), shape=3)
+        val = [1, 4, 5]
+        model = cp.Model([cp.GlobalCardinalityCount(iv, val, occ)])
+        solver = cp.SolverLookup.get("choco", model)
+        self.assertTrue(solver.solve())
+        self.assertTrue(cp.GlobalCardinalityCount(iv, val, occ).value())
+        self.assertTrue(all(cp.Count(iv, val[i]).value() == occ[i].value() for i in range(len(val))))
+        occ = [2, 3, 0]
+        model = cp.Model([cp.GlobalCardinalityCount(iv, val, occ), cp.AllDifferent(val)])
+        solver = cp.SolverLookup.get("choco", model)
+        self.assertTrue(solver.solve())
+        self.assertTrue(cp.GlobalCardinalityCount(iv, val, occ).value())
+        self.assertTrue(all(cp.Count(iv, val[i]).value() == occ[i] for i in range(len(val))))
+        self.assertTrue(cp.GlobalCardinalityCount([iv[0],iv[2],iv[1],iv[4],iv[3]], val, occ).value())
+
+    @pytest.mark.skipif(not CPM_choco.supported(),
+                        reason="pychoco not installed")
+    def test_choco_inverse(self):
+        from cpmpy.solvers.ortools import CPM_ortools
+
+        fwd = cp.intvar(0, 9, shape=10)
+        rev = cp.intvar(0, 9, shape=10)
+
+        # Fixed value for `fwd`
+        fixed_fwd = [9, 4, 7, 2, 1, 3, 8, 6, 0, 5]
+        # Inverse of the above
+        expected_inverse = [8, 4, 3, 5, 1, 9, 7, 2, 6, 0]
+
+        model = cp.Model(cp.Inverse(fwd, rev), fwd == fixed_fwd)
+
+        solver = cp.SolverLookup.get("choco", model)
+        self.assertTrue(solver.solve())
+        self.assertEqual(list(rev.value()), expected_inverse)
+
+    @pytest.mark.skipif(not CPM_choco.supported(),
+                        reason="pychoco not installed")
+    def test_choco_objective(self):
+        iv = cp.intvar(0,10, shape=2)
+        m = cp.Model(iv >= 1, iv <= 5, maximize=sum(iv))
+        s = cp.SolverLookup.get("choco", m)
+        self.assertTrue( s.solve() )
+        self.assertEqual( s.objective_value(), 10)
+
+        m = cp.Model(iv >= 1, iv <= 5, minimize=sum(iv))
+        s = cp.SolverLookup.get("choco", m)
+        self.assertTrue( s.solve() )
+        self.assertEqual(s.objective_value(), 2)
+
     @pytest.mark.skipif(not CPM_gurobi.supported(),
                         reason="Gurobi not installed")
     def test_gurobi_element(self):