how to fit a method belonging to an instance with pymc3?

Question

I failed to fit a method belonging to an instance of a class, as a Deterministic function, with PyMc3. Can you show me how to do that ?

For simplicity, my case is summarised below with a simple example. In reality, my constraint is that everything is made through a GUI and actions like ‘find_MAP’ should be inside methods linked to pyqt buttons.

I want to fit the function ‘FunctionIWantToFit’ over the data points. Problem, the following code:

import numpy as np
import pymc3 as pm3
from scipy.interpolate import interp1d
import theano.tensor as tt
import theano.compile

class cprofile:
    def __init__(self):
        self.observed_x = np.array([0.3,1.4,3.1,5,6.8,9,13.4,17.1])
        self.observations = np.array([6.25,2.75,1.25,1.25,1.5,1.75,1.5,1])
        self.x = np.arange(0,18,0.5)

    @theano.compile.ops.as_op(itypes=[tt.dscalar,tt.dscalar,tt.dscalar],
                              otypes=[tt.dvector])
    def FunctionIWantToFit(self,t,y,z):
        # can be complicated but simple in this example
        # among other things, this FunctionIWantToFit depends on a bunch of 
        # variables and methods that belong to this instance of the class cprofile,
        # so it cannot simply be put outside the class ! (like in the following example)
        val=t+y*self.x+z*self.x**2
        interp_values = interp1d(self.x,val)
        return interp_values(self.observed_x)

    def doMAP(self):
        model = pm3.Model()
        with model:
            t = pm3.Uniform("t",0,5)
            y = pm3.Uniform("y",0,5)
            z = pm3.Uniform("z",0,5)
            MyModel = pm3.Deterministic('MyModel',self.FunctionIWantToFit(t,y,z))
            obs = pm3.Normal('obs',mu=MyModel,sd=0.1,observed=self.observations)
            start = pm3.find_MAP()
            print('start: ',start)

test=cprofile()
test.doMAP()

gives the following error:

Traceback (most recent call last):

  File "<ipython-input-15-3dfb7aa09f84>", line 1, in <module>
    runfile('/Users/steph/work/profiles/GUI/pymc3/so.py', wdir='/Users/steph/work/profiles/GUI/pymc3')

  File "/Users/steph/anaconda/lib/python3.5/site-packages/spyder/utils/site/sitecustomize.py", line 866, in runfile
    execfile(filename, namespace)

  File "/Users/steph/anaconda/lib/python3.5/site-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/Users/steph/work/profiles/GUI/pymc3/so.py", line 44, in <module>
    test.doMAP()

  File "/Users/steph/work/profiles/GUI/pymc3/so.py", line 38, in doMAP
    MyModel = pm3.Deterministic('MyModel',self.FunctionIWantToFit(x,y,z))

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gof/op.py", line 668, in __call__
    required = thunk()

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gof/op.py", line 912, in rval
    r = p(n, [x[0] for x in i], o)

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/compile/ops.py", line 522, in perform
    outs = self.__fn(*inputs)

TypeError: FunctionIWantToFit() missing 1 required positional argument: 'z'

What’s wrong ?

remark 1: I systematically get an error message concerning the last parameter of ‘FunctionIWantToFit’. here it’s ‘z’ but if I remove z from the signature, the error message concerns ‘y’ (identical except from the name of the variable). if I add a 4th variable ‘w’ in the signature, the error message concerns ‘w’ (identical except from the name of the variable).

rk2: it looks like I missed something very basic in ‘theano’ or ‘pymc3’, because when I put ‘FunctionIWantToFit’ outside the class, it works. See the following example.

class cprofile:
    def __init__(self):
        self.observations = np.array([6.25,2.75,1.25,1.25,1.5,1.75,1.5,1])

    def doMAP(self):
        model = pm3.Model()
        with model:
            t = pm3.Uniform("t",0,5)
            y = pm3.Uniform("y",0,5)
            z = pm3.Uniform("z",0,5)
            MyModel = pm3.Deterministic('MyModel',FunctionIWantToFit(t,y,z))
            obs = pm3.Normal('obs',mu=MyModel,sd=0.1,observed=self.observations)
            start = pm3.find_MAP()
            print('start: ',start)

@theano.compile.ops.as_op(itypes=[tt.dscalar,tt.dscalar,tt.dscalar],
                              otypes=[tt.dvector])
def FunctionIWantToFit(t,y,z):
        observed_x = np.array([0.3,1.4,3.1,5,6.8,9,13.4,17.1])
        x = np.arange(0,18,0.5)
        val=t+y*x+z*x**2
        interp_values = interp1d(x,val)
        return interp_values(observed_x)

test=cprofile()
test.doMAP()

gives:

Warning: gradient not available.(E.g. vars contains discrete variables). MAP estimates may not be accurate for the default parameters. Defaulting to non-gradient minimization fmin_powell.
WARNING:pymc3:Warning: gradient not available.(E.g. vars contains discrete variables). MAP estimates may not be accurate for the default parameters. Defaulting to non-gradient minimization fmin_powell.
Optimization terminated successfully.
         Current function value: 1070.673818
         Iterations: 4
         Function evaluations: 179
start:  {'t_interval_': array(-0.27924150484602733), 'y_interval_': array(-9.940000425802811), 'z_interval_': array(-12.524909223913992)}

Except that I don’t know how to do that without big modifications in several modules, since the real ‘FunctionIWantToFit’ depends on a bunch of variables and methods that belong to this instance of the class profile.

In fact I 'm not even sure I know how to do that since ‘FunctionIWantToFit’ should then have objects in arguments (that I currently use via self) and I'm not sure how to do that with the theano decorator.

So I would prefer to avoid this solution... unless necessary. then I need explanations on how to implement it...

---

added on april 9, 2017:

Even without the interpolation question, it doesn't work because I must have missed something obvious with theano and/or pymc3. Please can you explain the problem ? I just want to compare model and data. First, it's such a shame being stuck to pymc2. ; second, I'm sure I'm not the only one with such a basic problem.

For example, let's consider variations around this very basic code:

import numpy as np
import theano
import pymc3
theano.config.compute_test_value = 'ignore'
theano.config.on_unused_input = 'ignore'

class testclass:
    x = np.arange(0,18,0.5)
    observed_x = np.array([0.3,1.4,3.1,5,6.8,9,13.4,17.1])
    observations = np.array([6.25,2.75,1.25,1.25,1.5,1.75,1.5,1])

    def testfunc(self,t,y,z):
        t2 = theano.tensor.dscalar('t2')
        y2 = theano.tensor.dscalar('y2')
        z2 = theano.tensor.dscalar('z2')
        val = t2 + y2 * self.observed_x + z2 * self.observed_x**2
        f = theano.function([t2,y2,z2],val)
        return f

test=testclass()
model = pymc3.Model()
with model:
    t = pymc3.Uniform("t",0,5)
    y = pymc3.Uniform("y",0,5)
    z = pymc3.Uniform("z",0,5)

with model:
   MyModel = pymc3.Deterministic('MyModel',test.testfunc(t,y,z))

with model:
   obs = pymc3.Normal('obs',mu=MyModel,sd=0.1,observed=test.observations)

this code fails at the last line with the error message: TypeError: unsupported operand type(s) for -: 'TensorConstant' and 'Function'

if I change 'testfunc' into:

def testfunc(self,t,y,z):
    t2 = theano.tensor.dscalar('t2')
    y2 = theano.tensor.dscalar('y2')
    z2 = theano.tensor.dscalar('z2')
    val = t2 + y2 * self.observed_x + z2 * self.observed_x**2
    f = theano.function([t2,y2,z2],val)
    fval = f(t,y,z,self.observed_x)
    return fval

the code fails at the 'MyModel =' line with error TypeError: ('Bad input argument to theano function with name "/Users/steph/work/profiles/GUI/pymc3/theanotest170409.py:32" at index 0(0-based)', 'Expected an array-like object, but found a Variable: maybe you are trying to call a function on a (possibly shared) variable instead of a numeric array?')

if I go back to the original 'testfunc' but change the last 'with model' lines with:

with model:
   fval = test.testfunc(t,y,z)
   obs = pymc3.Normal('obs',mu=fval,sd=0.1,observed=test.observations)

the error is the same as the first one.

I presented here only 3 tries but I would like to underline that I tried many many combinations, simpler and simpler until these ones, during hours. I have the feeling pymc3 shows a huge change of spirit, compared to pymc2, that I didn't get and is poorly documented...

Cromer

Answer 1

I finally converged toward the successful code below:

import numpy as np
import theano
from scipy.interpolate import interp1d
import pymc3 as pm3
theano.config.compute_test_value = 'ignore'
theano.config.on_unused_input = 'ignore'

class cprofile:
    observations = np.array([6.25,2.75,1.25,1.25,1.5,1.75,1.5,1])
    x = np.arange(0,18,0.5)
    observed_x = np.array([0.3,1.4,3.1,5,6.8,9,13.4,17.1])    

    def doMAP(self):
        model = pm3.Model()
        with model:
            t = pm3.Uniform("t",0,5)
            y = pm3.Uniform("y",0,5)
            z = pm3.Uniform("z",0,5)
            obs=pm3.Normal('obs',
              mu=FunctionIWantToFit(self)(t,y,z),
              sd=0.1,observed=self.observations)
            start = pm3.find_MAP()
            print('start: ',start)

class FunctionIWantToFit(theano.gof.Op):
    itypes=[theano.tensor.dscalar,
            theano.tensor.dscalar,
            theano.tensor.dscalar]
    otypes=[theano.tensor.dvector]

    def __init__(self, cp):
        self.cp = cp # note cp is an instance of the 'cprofile' class

    def perform(self,node, inputs, outputs):
        t, y, z = inputs[0], inputs[1], inputs[2]

        xxx = self.cp.x
        temp = t+y*xxx+z*xxx**2
        interpolated_concentration = interp1d(xxx,temp)   
        outputs[0][0] = interpolated_concentration(self.cp.observed_x)

testcp=cprofile()
testcp.doMAP()

thanks to the answer by Dario because I was too slow to understand the first answer by myself. I get it retrospectively but I strongly think the pymc3 doc is painfully unclear. It should contain very simple and illustrative examples.

However I didn’t succed in doing anything that work following the comment by Chris. Could anyone explain and/or give an example ?

One more thing: I don’t know whether my example above is efficient or could be simplified. In particular it gives me the impression the instance ‘testcp’ is copied twice in memory. More comments/answers are welcome to go further.

Cromer

Answer 2

Ok, let's do this by parts. First I'll explain the error messages that you got, and then I'll tell you how I would proceed.

On the first question, the direct reason why you're getting a complaint on the missing parameters is because your function, defined inside the class, takes as input (self, t, y, z), while you're declaring it in the op decorator as having only three inputs (t, y, z). You would have to declare the inputs as being four in your decorator to account for the class instance itself.

On "added on april 9, 2017:", the first code will not work because the output of test.testfunc(t,y,z) is a theano function itself. pymc3.Deterministic is expecting it to output theano variables (or python variables). Instead, make test.testfun output val = t2 + y2 * self.observed_x + z2 * self.observed_x**2 directly.

Then, on "if I change 'testfunc' into:", you get that error because of the way pymc3 is trying to work with theano functions. Long story short, the problem is that when pymc3 is making use of this function, it will send it theano variables, while fval is expecting numerical variables (numpy arrays or other). As in the previous paragraph, you just need to output val directly: no need to compile any theano function for this.

As for how I would proceed, I would try to declare the class instance as input to the theano decorator. Unfortunately, I can't find any documentation on how to do this and it might actually be impossible (see this old post, for example).

Then I would try to pass everything the function needs as inputs and define it outside of the class. This could be quite cumbersome and if it needs methods as input, then you run into additional problems.

Another way of doing this is to create a child class of theano.gof.Op whose init method takes your class (or rather an instance of it) as input and then define your perform() method. This would look something like this:

class myOp(theano.gof.Op):
    """ These are the inputs/outputs you used in your as_op
    decorator.
    """
    itypes=[tt.dscalar,tt.dscalar,tt.dscalar]
    otypes=[tt.dvector]
    def __init__(self, myclass):
        """ myclass would be the class you had from before, which
        you called cprofile in your first block of code."""
        self.myclass = myclass
    def perform(self,node, inputs, outputs):
        """ Here you define your operations, but instead of
        calling everyting from that class with self.methods(), you
        just do self.myclass.methods().

        Here, 'inputs' is a list with the three inputs you declared
        so you need to unpack them. 'outputs' is something similar, so
        the function doesn't actually return anything, but saves all
        to outputs. 'node' is magic juice that keeps the world
        spinning around; you need not do anything with it, but always
        include it.
        """
        t, y, z = inputs[0][0], inputs[0][1], inputs[0][2]
        outputs[0][0] = t+y*self.myclass.x+z*self.myclass.x**2
myop = myOp(myclass)

Once you have done this, you can use myop as your Op for the rest of your code. Note that some parts are missing. You can check my example for more details.

As for the example, you do not need to define the grad() method. Because of this, you can do all operations in perform() in pure python, if that helps.

Alternatively, and I say this with a smirk on my face, if you have access to the definition of the class you're using, you can also make it inherit from theano.gof.Op, create the perform() method (as in my other example, where you left a message) and try to use it like that. It could create conflicts with whatever else you're doing with that class and it's probably quite hard to get right, but might be fun to try.

Monticule

Answer 3

theano.compile.ops.as_op is just a short-hand for defining simple Theano Ops. If you want to code more involved ones, it is better to define it in a separate class. Objects of this class could of course take a reference to an instance of your cprofile if that really is necessary.

http://deeplearning.net/software/theano/extending/extending_theano.html

Gravure

thank you for your answer. I was away and will try this now. – Silicosis

sorry but it remains unclear to me. I tried many combinations of a home-made theano op with pymc3 but nothing works and the docs are painfully blurred... Could you provide an example of code ? – Silicosis

If you write your own interpolation function using Theano, you should not need an op at all. 1D interpolation should be straightforward. – Troublesome