z3-python: How to correctly handle substitution of candidate functions in Z3 for LIA problems?

Question

I'm new to z3, and have been working on a trying to solve LIA problems using Z3 (Z3 version 4.13.0 - 64 bit ) + python (3.12), and I think I need to substitute candidate expressions (which can be either expressions or integers) into a set of constraints. What I'm really trying to do is generate python functions which could take numbers of arguments e.g. x and y, and test them with z3 to see if they satisfy the constraints of the problem and are valid. These generated python functions could return an Int e.g.

def func_1(*args): 
    return 1

or something slightly more complex e.g.

def max_function(*args):
    if len(args) != 2:
        raise ValueError("max_function expects exactly 2 arguments.")
    x, y = values
    return If(x <= y, y, x)

A sample problem (but not limited to just this one) I'm trying to solve is the max of 2 ints:

(set-logic LIA)
(declare-fun f (Int Int) Int)
(declare-fun x () Int)
(declare-fun y () Int)
(assert (= (f x y) (f y x)))
(assert (and (<= x (f x y)) (<= y (f x y))))

As far as I can see, I need to substitute the candidate functions into the constraints in the solver. From previous answers I've understood that this approach is possibly the only way to do the following:

(declare-fun f (Int Int) Int) needs to become the following for func_1 (define-fun f (Int Int) Int(1)) and needs to become for func_2 (define-fun f (Int Int) Int (ite (<= x y) y x)) as I can't find a way to redefine the function i.e. go from declare-fun to define-fun, from what I've gathered, I need to specifically use python functions (not as z3 "Functions"), and substitute the candidate function/expression into the constraints. Which works for simple problems e.g. if I hardcode it, but I cannot seem to get it working for more complex ones in a more dynamic way i.e. with different constraints or numbers of variables (I've searched stackoverflow and tried to also incorporate a previous question's solution here). I've tried many different ways of doing this so can only assume that it is a fundamental misunderstanding. Here is one iteration that fails to substitute correctly for the commutativity constraint in the sample problem as a reproducible example:

def make_expression(possible_expression):
    if is_expr(possible_expression):
        return possible_expression
    elif isinstance(possible_expression, int):
        return IntVal(possible_expression)
    else:
        raise Exception("Cannot convert: %s" % possible_expression.__repr__())

def substitute_constraints(constraints: Collection[z3.ExprRef], func: z3.FuncDeclRef,
                           candidate_expression: typing.Union[z3.FuncDeclRef, z3.QuantifierRef, typing.Callable, int, z3.IntNumRef]) -> typing.List[z3.ExprRef]:
    def reconstruct_expression(expr: z3.ExprRef) -> z3.ExprRef:
        if is_app(expr) and expr.decl() == func:
            new_args = [reconstruct_expression(arg) for arg in expr.children()]
            if callable(candidate_expression):
                return candidate_expression(*new_args)
            elif isinstance(candidate_expression, (FuncDeclRef, QuantifierRef)):
                return candidate_expression(*new_args)
            else:
                return make_expression(candidate_expression)
        elif is_app(expr):
            new_args = [reconstruct_expression(arg) for arg in expr.children()]
            return expr.decl()(*new_args)
        else:
            return expr

    if isinstance(candidate_expression, (int, z3.IntNumRef)):
        substitution_pairs = [(func, make_expression(candidate_expression))]
    elif callable(candidate_expression):
        substitution_pairs = [(func, lambda *args: candidate_expression(*[make_expression(arg) for arg in args]))]
    else:
        substitution_pairs = [(func, candidate_expression)]

    return [substitute(c, *substitution_pairs) for c in constraints]

x, y = Ints('x y')
expr = x + y

def functionFromExpr(args, expr):
    return lambda *fargs: substitute(expr, *zip(args, [make_expression(e) for e in fargs]))

f = functionFromExpr([x, y], expr)
# These work ok
print(f(x, y))          # Should print: x + y
print(f(x, IntVal(4)))  # Should print: x + 4
print(f(IntVal(3), IntVal(4)))  # Should print: 3 + 4

# This does not work
x, y = Ints('x y')
f = Function('f', IntSort(), IntSort(), IntSort())

commutativity_constraint = f(x, y) == f(y, x)
print("Original commutativity constraint:")
print(commutativity_constraint)
def func_2(a, b):
    return If(a >= b, a, b)

substituted_constraints = substitute_constraints([commutativity_constraint], f, func_2)
for constraint in substituted_constraints:
    print(constraint)

substituted_constraints = substitute_constraints([commutativity_constraint], f, 42)
for constraint in substituted_constraints:
    print(constraint)

I usually get some variant of: z3.z3types.Z3Exception: Z3 invalid substitution, expression pairs expected.

I've tried using z3.Substitute, z3.Substitute_funs, z3.substitute_vars etc as documented here

I would like to avoid trying to do this with string manipulation in SMT-LIB if possible (generating candidate programs and resetting the solver and parsing the problem from a string), as I understand the recommendation is to not mix SMT-LIB and python i.e. stick to one or the other, and to try and use python functions for this kind of problem.

Please let me know if any of the assertions or assumptions I have made are incorrect (I've gathered them from previous answers along the way, though I may have misunderstood), hopefully that makes sense, any help would be greatly appreciated.