You need to wrap your python code in a pytensor Op. See here for an example that focuses on likelihood functions. Your case will be the same in terms of constructing an Op, but you will call the function on symbolic inputs rather than put it in a CustomDist logp. Here is another example where I wrap scipy.optimize.root in an Op and use it in a model, which is close to how you will use your external functions I guess.