Using TensorVariable into interpolated function with numpy

Hello everyone.
I hope I can show the problem I have using mainly (I think) pytensor and numpy.

I want to make inference over some physical parameter x which is a variable of a complicated physical model.
I have tried to make the inference using the physical model as it is but the graphs to be constructed are so complex that my desktop computer runs out of memory.

Therefore I have thought that if I generate an interpolating polynomial for the physical model and use it as an alternative I could solve the problem. But here a new problem arose: when I want to generate samples of the interpolated model, there is an incompatibility between the tensor variable used by pymc and the numerical input expected by the interpolating polynomial generated with numpy (the same problem occurs using scipy interpolators).

I suspect that the solution can be given by translating the tensor variable to a float variable, but I have not found the correct way to do it.

Below I copy a minimal outline of the code I have implemented. I hope someone knows how to deal with this kind of issues.

Thank you very much.

import pymc as pm
import pytensor
import pytensor.tensor as pt

import numpy as np

from multiprocessing import Pool

1. Create the interpolated function evaluate a complicated physical model in a 1D-grid using p.map


 >>> def f_model(x):
    	      # a complicated function of x and another (fixed for now) parameters
    	     return f
>>> x_list = [ ]
>>> xRange = np.linspace(xmin, xmax, nx) # for some values of xmin, xmax and nx

>>> for j in range(nx):
    	    x_list.append(xRange[j])
  
>>> p = Pool(6)
>>> f_sampled = p.map(f_model,x_list)

# generate the interpolation

>>> x_new = np.linspace(xmin, xmax, Nx) # same values for xmin and xmax; but now Nx much bigger than nx

>>> f_interp = np.interp(x_new, xRange, f_sampled)

here I plot together f_interp and f_sampled and the interpolation seems to be satisfactory. Also, I print some evaluation for the interpolation and I get

>>> f_out = np.interp(x_value,xRange,f_sampled) # for some numerical x_value
>>> print(str(f_out),type(f_out))

f_value <class 'numpy.float64'>
**2. Generate the inference for x given observed data**

## Define the deterministic model and the prior for x


>>> def data_model(x):

        S1 = np.interp(x,xRange,f_sampled)
        Dim = 1
	    a = pt.zeros(Dim) 
	    a = S1
	    return a
>>> print('some evaluation: ', str(data_model(x_value))+'  '+str(type(data_model(x_value))))

some evaluation:  numerical_value  <class 'numpy.float64'>

## Then I define the Model:

>>> with pm.Model() as model_bayes:
    	x_in = pm.Uniform('x', xmin, xmax)

    	def sampled_model(x_rnd = x_in):
    	    return data_model(x_rnd)
        
>>> print('10 draws for x: ', str(pm.draw(x_in,10)))

>>> print('evaluation of the physical model:', str(sampled_model(pm.draw(x_in,1))))
>>> print('kind of variable: ', type(sampled_model(pm.draw(x_in,1))))

10 draws for x:  [x0 x1 x2 x3 x4 x5 x6 x7 x8 x9]

evaluation of the physical model: a_numerical_value
kind of variable:  <class 'numpy.float64'>

 And here comes the error

>>> with model_bayes:
    	function_pm = pm.Deterministic('modelo', sampled_model())


TypeError                                 Traceback (most recent call last)
<ipython-input-41-cd6021848196> in <cell line: 10>()
      9 
     10 with model_bayes:
---> 11     function_pm = pm.Deterministic('modelo',sampled_model())

1 frames
<ipython-input-38-3c75629263ab> in data_model(x)
     12     
     13 
---> 14     x_in = float(x)
     15 
     16     S1 = np.interp(x_in,xRange,f_sampled)

**TypeError: float() argument must be a string or a real number, not 'TensorVariable'**Preformatted text

pytensor tensors are not numbers, they are symbolic objects used to construct a computational graph. You cannot use them with numpy or scipy, or any other package for that matter. See here for an introduction to what pytensor is and how PyMC uses it.

You have two options:

  1. Write your interpolation function using pytensor functions. We don’t have a pt.interp available, but looking at the numpy source code, we should have all the tools to do it using pytensor functions.
  2. Wrap the numpy/scipy code in an pytensor Op, following e.g. this example.

More discussion on this here.

Hello Jesse!

Thank you very much for your recommendations. I have followed the post on using a “Using a black box …” and it seems to work fine for my case. I was able to obtain a Deterministic function by defining a ‘polynomial Operator’ in a similar way to the ‘LogLike’ example in the post. And also I was able to sample the model using pm.sample, but when I tried to use pm.smc.sample_smc, I encountered the following error:

UnusedInputError: Pytensor.function was asked to create a function computing outputs given certain inputs, but the provided input variable at index 0 is not part of the computational graph needed to compute the outputs: joined_inputs. To convert this error into a warning, you can pass the parameter on_unused_input='warn' to Pytensor.function. To disable it completely, use on_unused_input='ignore'.

I will continue reading about this issue. Thanks again for your suggestions.

Best wishes,

Mariano