Translating python loop to Theano

Question

I want to transform a P vector of n elements into a Z vector of n+1 elements using the following formula:

I was able to do that in python with:

import numpy as np

def p_to_Z(p):
    k = len(p)
    Z = np.zeros(k+1)
    Z[0] = p[0]

    Z[1:k] = [p[i] / (1 - sum(p[range(i)])) for i in range(1,k)]
    Z[-1] = 1
    return Z

Example:

p_ex = np.array([0.3,0.4,0.1,0.2])
p_to_Z(p_ex)
output:
array([ 0.3, 0.57142857, 0.33333333,  1. , 1. ])

and for a matrix with:

M = np.array([[0.2, 0.05, 0.1],
                     [0.5, 2.0, 0.2],
                     [0.1, 0.2, 1.0]])

np.apply_along_axis(p_to_Z, 1, M)

array([[ 0.2 , 0.0625 , 0.13333333, 1.],
       [ 0.5 , 4. ,-0.13333333, 1.],
       [ 0.1 , 0.22222222,  1.42857143, 1.]])

I have been trying to translate this into theano but I have been getting my ass kicked by the scan function. Could some one help me define these functions in theano? I would be forever grateful.

EDIT: Failed attempts

import theano.tensor as tt
import numpy as np
import theano
from theano.compile.ops import as_op

p = tt.vector('p')
results, updates = theano.scan(lambda i: p[i] / (1 - tt.sum(p[:i-1])), sequences=tt.arange(1,N))
fn = theano.function([p,N],results)

p_ex = np.array([0.3,0.4,0.1,0.2])
Z = fn(p_ex,tt.get_vector_length(p_ex))

and with decorator:

@as_op(itypes=[tt.dvector],otypes=[tt.dvector])
def p_to_Z_theno(p):
    N = tt.get_vector_length(p)
    Z= theano.scan(lambda i: p[i] / (1 - tt.sum(p[:i - 1])), sequences=tt.arange(1, N))
    fn = theano.function([p, N], results)
    Z[0]=p[0]
    Z[-1]=Z
    return(Z)

The first one obviously do not generate the expected results, while the second gives:

Traceback (most recent call last):   File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2882, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)   File "<ipython-input-21-211efc53b449>", line 1, in <module>
    p_to_Z_theno(p_ex)   File "/usr/local/lib/python3.6/dist-packages/theano/gof/op.py", line 615, in __call__
    node = self.make_node(*inputs, **kwargs)   File "/usr/local/lib/python3.6/dist-packages/theano/gof/op.py", line 983, in make_node
    if not all(inp.type == it for inp, it in zip(inputs, self.itypes)):   File "/usr/local/lib/python3.6/dist-packages/theano/gof/op.py", line 983, in <genexpr>
    if not all(inp.type == it for inp, it in zip(inputs, self.itypes)): AttributeError: 'numpy.ndarray' object has no attribute 'type'

When you say "translte to Theano", do you mean use vectorized numpy operations? — Mad Physicist, Aug 10 '18 at 18:16
Also, your formula only gives k elements, not k+1. Did you mean `for i = 2, ..., k, and Z_k+1 = 1`? — Mad Physicist, Aug 10 '18 at 18:18
Sorry, I edited the post. k = n+1. So, what I am doing actually is to implement a bayesian model in pymc3. Pymc3 works with theano, and I have defined a prior on P, but the likelihood function must be on Z. Thus, with this function, I want to be able to input a theano variable and return a theano variable, as well as being able to extend this to a matrix. — user2597079, Aug 10 '18 at 18:27
"I will not even put what I have tried, since I have failed miserably.". Please post your attempt, and a detailed explanation of the failure. — Mad Physicist, Aug 10 '18 at 18:40
I'll wait for you to post the attempt to see what you did. I am not super familiar with Theano, but I came up with something very simple that seems to work based on the docs. — Mad Physicist, Aug 10 '18 at 18:50
I don't think it's even worth it. I am both newbie in python and theano, and what I did wasn't even remotely worth of improvement. I think it would be nice to see what you have come up with. — user2597079, Aug 10 '18 at 18:55
When I say not worth of improvement, I mean the fucntion didn't even run. — user2597079, Aug 10 '18 at 18:55
That's fine. We're all here to learn. Post what you have and show the error. It won't hurt and will probably benefit you. I'll wait. — Mad Physicist, Aug 10 '18 at 18:56
Take a look at my answer BTW. You almost certainly don't need scan. Rather, you need something like cumsum, since Theano is at least as good at doing element-wise operations as numpy. — Mad Physicist, Aug 10 '18 at 18:58
There, I have put my attempts. Yeah, with your answer I am now seeing the problem in a different way. I will give a try with cumsum — user2597079, Aug 10 '18 at 19:13
Sorry, would you mind formatting the error as code please? It will make it much easier to read. — Mad Physicist, Aug 10 '18 at 19:22
I don't feel like delving into Theano enough to fix what you did, but I hope that the next person looking at this knows what went wrong at a glance. In the meantime, I've updated my answer based purely on the docs (which are excellent as I've found). I've never used Theano before. — Mad Physicist, Aug 10 '18 at 19:42
Unless you have a reason that compels you to use Theano specifically, I would strongly recommend using some other library. Theano is no longer under development. — PMende, Aug 10 '18 at 19:47

Mad Physicist · Accepted Answer · 2018-08-10T19:44:59.947

Let's start by thinking about this in terms of vectorized operations. Your initial equation can be written very simply if we use numpy's cumsum function:

def p2Z_numpy(p):
    Z = np.r_[p, 1]
    Z[1:-1] /= np.cumsum(p[:-1])
    return Z

For example:

>>> p2Z_numpy([0.3, 0.4, 0.1, 0.2])
array([0.3       , 1.33333333, 0.14285714, 0.25      , 1.        ])

The reason that is an important starting point is that Theano allows for all the same sort of operations, but on a much larger scale. So let's do the same division, concatenation and cumulative sum in Theano:

import theano.tensor as tt
import theano.tensor.extra_ops as te
from theano import function

p = tt.vector()
denom = 1 - te.cumsum(p[:-1])
frac = p[1:] / denom
Z = tt.concatentat([p[0:1], frac, [1]])

p2Z = function([p], Z)

Testing it out:

>>> p2Z([0.3, 0.4, 0.1, 0.2])
array([0.3       , 0.57142857, 0.33333333, 1.        , 1.        ])

While I've broken the function into steps, you could actually turn it into a massively ugly one-liner:

p = tt.vector()
Z = tt.concatenate([p[0:1], p[1:] / (1 - te.cumsum(p[:-1])), [1]])
p2Z = function([p], Z)

Or even

p2Z = function([p], tt.concatenate([p[0:1], p[1:] / (1 - te.cumsum(p[:-1])), [1]]))

Thats nice, thank you. However, I need a function that would work with theano variables. — user2597079, Aug 10 '18 at 18:30

Translating python loop to Theano

1 Answers1