Following the theano Tutorial


In [9]:
%matplotlib inline

In [1]:
from theano import *
import theano.tensor as T

Baby Steps - Algebra

Adding two Scalars


In [2]:
import numpy
from theano import function

In [3]:
x = T.dscalar('x')
y = T.dscalar('y')
z = x+y
f = function([x,y],z)

In [7]:
print type(x), type(y), type(z), type(f) # good to know what these 
# new classes are in theano


<class 'theano.tensor.var.TensorVariable'> <class 'theano.tensor.var.TensorVariable'> <class 'theano.tensor.var.TensorVariable'> <class 'theano.compile.function_module.Function'>

In [8]:
f(2,3)


Out[8]:
array(5.0)

In [9]:
numpy.allclose(f(16.3,12.1),28.4)


Out[9]:
True

In [10]:
x.type


Out[10]:
TensorType(float64, scalar)

In [11]:
T.dscalar


Out[11]:
TensorType(float64, scalar)

In [12]:
x.type is T.dscalar


Out[12]:
True

In [ ]:

"Prefer constructors like matrix, vector and scalar to dmatrix, dvector and dscalar because the former will give you float32 variables when floatX=float32." - cf. Using the GPU Theano tutorial


In [13]:
xf = T.scalar('xf')
yf = T.scalar('yf')
zf = xf + yf
ff = function([xf,yf],zf)

In [14]:
from theano import pp

In [15]:
print(pp(z))


(x + y)

In [16]:
print(pp(zf))


(xf + yf)

In [17]:
x = T.dmatrix('x')
y = T.dmatrix('y')
z = x + y 
f = function([x,y],z)

In [18]:
f([[1,2],[3,4]], [[10,20],[30,40]])


Out[18]:
array([[ 11.,  22.],
       [ 33.,  44.]])

In [19]:
f(numpy.array([[1,2],[3,4]]),numpy.array([[10,20],[30,40]]))


Out[19]:
array([[ 11.,  22.],
       [ 33.,  44.]])

In [20]:
xf = T.matrix('xf')
xy = T.matrix('yf')
zf = xf + yf
ff = function([xf,yf],zf)

In [22]:
ff([[1,2],[3,4]], [[10,20],[30,40]])


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-22-38fbc8fb8a3b> in <module>()
----> 1 ff([[1,2],[3,4]], [[10,20],[30,40]])

/home/topolo/Public/anaconda2/lib/python2.7/site-packages/theano/compile/function_module.pyc in __call__(self, *args, **kwargs)
    784                         s.storage[0] = s.type.filter(
    785                             arg, strict=s.strict,
--> 786                             allow_downcast=s.allow_downcast)
    787 
    788                     except Exception as e:

/home/topolo/Public/anaconda2/lib/python2.7/site-packages/theano/tensor/type.pyc in filter(self, data, strict, allow_downcast)
    175             raise TypeError("Wrong number of dimensions: expected %s,"
    176                             " got %s with shape %s." % (self.ndim, data.ndim,
--> 177                                                         data.shape))
    178         if not data.flags.aligned:
    179             try:

TypeError: ('Bad input argument to theano function with name "<ipython-input-20-f711db7fb599>:4"  at index 1(0-based)', 'Wrong number of dimensions: expected 0, got 2 with shape (2, 2).')

In [25]:
a = theano.tensor.vector()
b = theano.tensor.vector()
out = a**2 + b**2 + 2 * a * b
f = theano.function([a,b],out)
print(f([1,2],[4,5]))


[ 25.  49.]

"At this point it would be wise to begin familiarizing yourself more systematically with Theano’s fundamental objects and operations by browsing this section of the library: Basic Tensor Functionality." cf. More Examples

Custom tensor types


In [28]:
dtensor5 = T.TensorType('float64', (False,)*5)

In [29]:
x = dtensor5()
z = dtensor5('z')

In [32]:
my_dmatrix = T.TensorType('float64', (False,)*2)
x = my_dmatrix()
my_dmatrix == T.dmatrix


Out[32]:
True

Converting from Python Objects


In [33]:
x = shared(numpy.random.randn(3,4))

In [34]:
x = T.dmatrix('x')
s = 1 / ( 1 + T.exp(-x))
logistic = theano.function([x],s)
logistic([[0,1],[-1,-2]])


Out[34]:
array([[ 0.5       ,  0.73105858],
       [ 0.26894142,  0.11920292]])
\begin{gathered} s(x) = \frac{1}{1+\exp{-x} } = \frac{1+\tanh{(x/2) } }{2} \end{gathered}

In [35]:
s2 = (1 + T.tanh(x/2))/2
logistic2 = theano.function([x],s2)
logistic2([[0,1],[-1,-2]])


Out[35]:
array([[ 0.5       ,  0.73105858],
       [ 0.26894142,  0.11920292]])

Computing More than one Thing at the Same Time (!!!)


In [36]:
a,b = T.dmatrices('a','b')
diff = a-b
abs_diff = abs(diff)
diff_squared = diff**2
f = theano.function([a,b],[diff,abs_diff,diff_squared])

In [37]:
f([[1,1],[1,1]], [[0,1],[2,3]])


Out[37]:
[array([[ 1.,  0.],
        [-1., -2.]]), array([[ 1.,  0.],
        [ 1.,  2.]]), array([[ 1.,  0.],
        [ 1.,  4.]])]

Setting a Default Value for an Argument


In [38]:
from theano import In
from theano import function
x,y = T.dscalars('x','y')
z = x + y
f= function([x,In(y,value=1)],z)

In [39]:
f(33)


Out[39]:
array(34.0)

In [40]:
f(33,2)


Out[40]:
array(35.0)

"Inputs with default values must follow inputs without default values (like Python’s functions). There can be multiple inputs with default values. These parameters can be set positionally or by name, as in standard Python:"


In [42]:
x,y,w = T.dscalars('x', 'y', 'w')
z = (x+y)*w
f = function([x,In(y,value=1),In(w,value=2,name='w_by_name')],z)
f(33)


Out[42]:
array(68.0)

In [43]:
f(33,2)


Out[43]:
array(70.0)

In [44]:
f(33,0,1)


Out[44]:
array(33.0)

In [45]:
f(33,w_by_name=1)


Out[45]:
array(34.0)

In [46]:
f(33,w_by_name=1,y=0)


Out[46]:
array(33.0)

Using Shared Variables


In [47]:
from theano import shared
state = shared(0)
inc = T.iscalar('inc')
accumulator = function([inc],state,updates=[(state,state+inc)])

In [48]:
print(state.get_value())


0

In [49]:
accumulator(1)


Out[49]:
array(0)

In [50]:
print(state.get_value())


1

In [51]:
accumulator(300)


Out[51]:
array(1)

In [52]:
print(state.get_value())


301

"It is possible to reset the state. Just use the .set_value() method:"


In [53]:
state.set_value(-1)

In [54]:
accumulator(3)


Out[54]:
array(-1)

In [55]:
print(state.get_value())


2

In [56]:
decrementor = function([inc],state, updates=[(state,state-inc)])

In [57]:
decrementor(2)


Out[57]:
array(2)

In [58]:
print(state.get_value())


0

"Also, Theano has more control over where and how shared variables are allocated, which is one of the important elements of getting good performance on the GPU."


In [60]:
fn_of_state = state * 2 + inc
# The type of foo must match the shared variable we are replacing
# with the "givens"
foo = T.scalar(dtype=state.dtype)
skip_shared = function([inc,foo], fn_of_state, givens=[(state,foo)])
skip_shared(1,3)


Out[60]:
array(7)

In [61]:
print(state.get_value())


0

Copying functions


In [64]:
inc = T.iscalar('inc')
accumulator = theano.function([inc],state, updates=[(state,state+inc)])
accumulator(10)


Out[64]:
array(0)

In [65]:
print(state.get_value())


10

"We can use copy() to create a similar accumulator but with its own internal state using the swap parameter, which is a dictionary of shared variables to exchange:"


In [66]:
new_state = theano.shared(0)
new_accumulator = accumulator.copy(swap={state:new_state})
new_accumulator(100)


Out[66]:
[array(0)]

In [67]:
print(new_state.get_value())


100

In [68]:
print(state.get_value())


10

In [69]:
null_accumulator = accumulator.copy(delete_updates=True)
null_accumulator(9000)


---------------------------------------------------------------------------
UnusedInputError                          Traceback (most recent call last)
<ipython-input-69-9c62592c8b05> in <module>()
----> 1 null_accumulator = accumulator.copy(delete_updates=True)
      2 null_accumulator(9000)

/home/topolo/Public/anaconda2/lib/python2.7/site-packages/theano/compile/function_module.pyc in copy(self, share_memory, swap, delete_updates, name, profile)
    719                                 # can contain inplace. DebugMode check
    720                                 # that.
--> 721                                 accept_inplace=True,
    722                                 ).create(input_storage,
    723                                          storage_map=new_storage_map)

/home/topolo/Public/anaconda2/lib/python2.7/site-packages/theano/compile/function_module.pyc in __init__(self, inputs, outputs, mode, accept_inplace, function_builder, profile, on_unused_input, fgraph, output_keys)
   1413 
   1414         # Check if some input variables are unused
-> 1415         self._check_unused_inputs(inputs, outputs, on_unused_input)
   1416 
   1417         # Make a list of (SymbolicInput|SymblicInputKits, indices,

/home/topolo/Public/anaconda2/lib/python2.7/site-packages/theano/compile/function_module.pyc in _check_unused_inputs(self, inputs, outputs, on_unused_input)
   1551                 elif on_unused_input == 'raise':
   1552                     raise UnusedInputError(msg % (inputs.index(i),
-> 1553                                                   i.variable, err_msg))
   1554                 else:
   1555                     raise ValueError("Invalid value for keyword "

UnusedInputError: theano.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: inc.
To make this error into a warning, you can pass the parameter on_unused_input='warn' to theano.function. To disable it completely, use on_unused_input='ignore'.

In [70]:
print(state.get_value())


10

Using Random Numbers


In [71]:
from theano.tensor.shared_randomstreams import RandomStreams
from theano import function
srng = RandomStreams(seed=234)
rv_u = srng.uniform((2,2)) # represents a random stream of 2x2 matrices
rv_n = srng.normal((2,2))
f = function([], rv_u)
g = function([], rv_n, no_default_updates=True) # Not updating rv_n.rng
nearly_zeros = function([],rv_u+rv_u - 2 * rv_u)

"The RandomStream only work on the CPU, MRG31k3p work on the CPU and GPU. CURAND only work on the GPU." cf. http://deeplearning.net/software/theano/tutorial/examples.html#other-implementations


In [74]:
from theano.sandbox.rng_mrg import MRG_RandomStreams

In [76]:
from theano.sandbox.cuda import CURAND_RandomStreams

In [77]:
f_val0 = f()

In [78]:
f_val1 = f()

"When we add the extra argument no_default_updates=True to function (as in g), then the random number generator state is not affected by calling the returned function. So, for example, calling g multiple times will return the same numbers."


In [79]:
g_val0 = g() # different numbers from f_val0 and f_val1

In [80]:
g_val1 = g()

"An important remark is that a random variable is drawn at most once during any single function execution. So the nearly_zeros function is guaranteed to return approximately 0 (except for rounding error) even though the rv_u random variable appears three times in the output expression."

Seeding Streams


In [82]:
rng_val = rv_u.rng.get_value(borrow=True) # Get the ring for rv_u
rng_val.seed(89234) # seeds the generator
rv_u.rng.set_value(rng_val, borrow=True) # Assign back seeded rng

In [83]:
srng.seed(902340) # seeds rv_u and rv_n with different seeds each

Sharing Streams Between Functions


In [84]:
state_after_v0 = rv_u.rng.get_value().get_state()

In [85]:
nearly_zeros()  # this affects rv_u's generator


Out[85]:
array([[ 0.,  0.],
       [ 0.,  0.]])

In [86]:
v1 = f()
rng = rv_u.rng.get_value(borrow=True)
rng.set_state(state_after_v0)
rv_u.rng.set_value(rng,borrow=True)

In [87]:
v2 =f() # v2 != v1

In [88]:
v3=f() # v3 == v1

In [90]:
v2.view()


Out[90]:
array([[ 0.33919835,  0.85344878],
       [ 0.14881562,  0.79659413]])

In [91]:
v1.view()


Out[91]:
array([[ 0.5025809 ,  0.99544429],
       [ 0.75073355,  0.17926032]])

In [92]:
v3.view()


Out[92]:
array([[ 0.5025809 ,  0.99544429],
       [ 0.75073355,  0.17926032]])

Copying Random State Between Theano Graphs


In [93]:
from __future__ import print_function
from theano.sandbox.rng_mrg import MRG_RandomStreams
from theano.tensor.shared_randomstreams import RandomStreams

In [94]:
class Graph():
    def __init__(self, seed=123):
        self.rng = RandomStreams(seed)
        self.y = self.rng.uniform(size=(1,))

In [96]:
g1 = Graph(seed=123)
f1 = theano.function([], g1.y)
g2 = Graph(seed=987)
f2 = theano.function([], g2.y)

# By default, the two functions are out of sync.
f1()


Out[96]:
array([ 0.72803009])

In [97]:
f2()


Out[97]:
array([ 0.55056769])

In [98]:
def copy_random_state(g1,g2):
    if isinstance(g1.rng, MRG_RandomStreams):
        g2.rng.rstate = g1.rng.rstate
    for (su1, su2) in zip(g1.rng.state_updates, g2.rng.state_updates):
        su2[0].set_value(su1[0].get_value())

In [99]:
# We now copy the state of the theano random number generators.
copy_random_state(g1, g2)
f1()


Out[99]:
array([ 0.59044123])

In [100]:
f2()


Out[100]:
array([ 0.59044123])

Other Random Distributions

are found here at other distributions implemented

A Real Example: Logistic Regression


In [101]:
import numpy
import theano
import theano.tensor as T
rng = numpy.random

N = 400      # training sample size
feats = 784  # number of input variables

# generate a dataset: D = (input_values, target_class)
D = (rng.randn(N, feats), rng.randint(size=N, low=0, high=2))
training_steps = 10000

# Declare Theano symbolic variables
x = T.dmatrix("x")
y = T.dvector("y")

# initialize the weight vector w randomly
#
# this and the following bias variable b
# are shared so they keep their values
# between training iterations (updates)
w = theano.shared(rng.randn(feats), name="w")

# initialize the bias term
b = theano.shared(0., name="b")

print("Initial model:")
print(w.get_value())
print(b.get_value())

# Construct Theano expression graph
p_1 = 1 / (1 + T.exp(-T.dot(x, w) - b))  # Probability that target =  1
prediction = p_1 > 0.5                   # The prediction thresholded
xent = -y * T.log(p_1) - (1-y * T.log(1-p_1)) # Cross-entropy loss function
cost = xent.mean() + 0.01 * ( w** 2).sum() # The cost to minimize
gw, gb = T.grad(cost, [w,b])             # Compute the gradient of the cost
                                    # w.r.t weight vector w and 
                                     # bias term b
                                    # (we shall return to this in a 
                                    # following section of this tutorial)

# Compile
train = theano.function(
            inputs=[x,y],
            outputs=[prediction, xent],
            updates=((w,w-0.1 *gw), (b,b-0.1 * gb)))
predict = theano.function(inputs=[x], outputs=prediction)

# Train
for i in range(training_steps):
    pred, err = train(D[0], D[1])
    
print("Final model:")
print(w.get_value())
print(b.get_value())
print("target values for D:")
print(D[1])
print("prediction on D:")
print(predict(D[0]))


Initial model:
[  1.02090375e+00   8.73947142e-01   1.95614041e+00   7.07377428e-01
  -6.53816172e-01   7.42545996e-01  -9.80754517e-01   5.61128882e-01
   1.78775715e+00   2.55831279e-01  -8.67444321e-02   7.94001501e-01
  -6.70527935e-01  -1.96668538e+00   2.54200412e+00   8.89567452e-01
   9.74970709e-01   8.48005589e-01  -6.09636316e-01   8.50633037e-01
   1.68271257e+00   4.93346994e-01  -5.51691867e-02   2.37368331e-01
   1.62746146e-01   1.55918000e+00  -8.34594814e-01   1.38944861e+00
  -1.23368104e+00  -3.93048716e-01   1.05566395e+00  -1.31731246e+00
   7.08793200e-01  -2.69903059e-01   1.59489513e+00  -5.41472996e-01
  -1.51657778e+00  -1.74491648e+00  -4.87264078e-01  -1.40829581e+00
   9.16277590e-01   1.13288831e-01  -8.15545784e-01   3.63391354e-01
   4.59170798e-01  -4.92534498e-01   1.59924371e-01  -5.79592495e-01
   8.28797390e-01   1.01823000e+00   9.47734065e-01   1.07973174e+00
   4.30089731e-02  -8.06497075e-02   8.95636221e-01  -1.18478568e-02
  -5.81466361e-01   4.24213296e-01  -9.76244057e-01   9.77803579e-01
  -1.68833204e-01   7.43143151e-01  -1.45011683e+00  -4.39116887e-01
   2.02524479e-01  -1.48128730e-01  -1.57771132e+00   1.50352107e-01
  -3.28331728e-01   6.71897989e-02  -1.19704026e+00   1.52688204e+00
  -1.31571397e+00   8.90164259e-01  -2.23393343e-01   4.99356497e-01
  -1.29504971e+00   7.02042318e-01  -3.24683905e-01  -6.27943493e-01
   3.14126673e-01   1.12213995e-01  -9.21144105e-01  -2.43648932e-01
  -2.33448211e-01  -1.00637079e+00   6.33923303e-01   6.56206738e-01
   1.45152861e+00   1.33058182e+00   5.10506128e-01   1.16950441e+00
  -7.71856464e-01   4.80736974e-01   2.31723934e-01   1.49188536e+00
  -7.33637390e-01  -3.60145054e-01  -6.62364149e-01   1.15997382e+00
   4.85467482e-02  -2.49511990e-01  -2.16758659e+00  -8.67892909e-01
   9.63635308e-01   5.77244652e-03   1.95557853e+00   6.59938526e-01
  -5.92379578e-01  -1.42131361e+00  -6.06090501e-01  -1.44791105e-01
  -4.07448243e-01   1.08839737e+00  -6.09163649e-01  -3.83439487e-01
  -5.52671444e-02  -7.45281579e-01   2.36207307e-01   6.73715161e-01
  -5.06680836e-01  -3.89571760e-01   5.43189401e-01   5.18757744e-01
  -8.18926229e-01  -1.31366774e-01   1.24943697e+00  -2.04555766e+00
   2.30376888e-01  -1.82374680e+00  -1.49620886e+00   7.66641416e-01
   4.93167826e-01   4.79110747e-01  -1.86664191e+00  -1.09007317e+00
   1.54655862e+00  -4.69143883e-01   2.97751272e-01  -1.89391756e+00
   3.17911849e-01   1.95615980e+00   9.86954201e-02   1.05970291e+00
  -2.14845759e+00   1.16862930e-01  -5.98409143e-01   5.44959433e-01
   1.69291338e+00  -2.03597258e-01   5.04499321e-02  -1.08120920e+00
   4.22101164e-01   1.82562928e+00   3.55452430e-01  -6.40125897e-01
   3.88448160e-01   7.04888678e-01  -1.51880862e+00   8.31411212e-01
  -1.27244052e+00   1.18952583e+00   1.08579966e-01  -9.32388052e-01
  -4.74739101e-01   1.52918600e+00  -5.25556385e-01  -4.14615352e-01
   3.18785478e-01  -9.88229458e-01  -8.35711414e-01  -5.70864096e-01
  -4.84650587e-01   5.56635563e-01   1.28846167e+00   5.66758002e-02
   1.55206322e+00   1.38991673e+00   7.79811803e-01   1.64131920e+00
   8.81968840e-01   9.24259617e-01  -1.11325013e-01   5.30073122e-01
   7.39416638e-01  -3.30580630e-01   5.22214454e-01   1.38711119e-01
   2.70007039e+00   1.13349621e+00   2.03398657e+00   4.02464985e-02
   5.30575756e-01  -7.19422036e-01   1.61156869e-01   7.02390558e-01
  -1.94458803e-01  -7.07393385e-01  -6.16683694e-01  -7.03852782e-01
  -2.13252703e-01  -1.26214924e+00   7.86406827e-01   9.48912220e-01
   9.61491556e-01   1.13303970e-01  -7.11180757e-01   1.33364223e+00
   7.51945657e-01   1.36873393e+00  -7.33127552e-01   1.54788909e-01
  -8.81978191e-01   2.39732742e-01   3.16681446e-01   1.25647817e+00
   9.86339419e-01  -8.00005460e-01  -1.58500220e+00   2.25529851e+00
  -2.57658643e-01  -8.35908448e-01   5.37293722e-01   1.46937651e+00
   4.22615189e-01  -9.25028307e-01  -1.08224378e+00  -4.17584822e-01
  -1.65624617e+00  -1.20110360e+00  -1.71466593e+00   1.06068512e+00
  -7.10531288e-01  -2.45182989e-01  -8.73764848e-01  -3.14402273e-01
   2.84017903e-01   1.04617480e+00  -2.00541051e-01  -6.76963341e-01
   6.42284249e-01  -7.27256802e-01   5.00257949e-01   1.77173328e-01
  -3.35200393e-01   4.54306408e-01   9.71647231e-01  -1.58096887e+00
  -1.62698537e+00  -1.20757309e-01  -4.84239650e-01  -1.01253549e+00
  -4.55594449e-01  -1.24059703e+00   1.45024757e+00  -6.36831201e-01
   6.66481267e-01  -1.20226156e+00   6.26097700e-01   4.36649752e-01
  -1.90033855e-01  -2.16185075e-01  -1.45658116e-01  -5.42178438e-01
  -9.51822596e-01   6.03622035e-01   3.47322478e+00   3.80057626e+00
   4.30018315e-01  -1.18529132e+00   1.05549365e+00   1.08000352e+00
   1.14564527e+00   9.32282009e-01   3.22643679e-01  -1.29549056e-01
  -5.48090410e-01   2.07474158e-01   1.23735892e+00  -1.40366607e-01
  -8.79224625e-02   1.05414726e+00   1.10139014e+00  -1.10684451e+00
  -8.03057691e-01  -1.36056405e+00  -6.74985268e-01   5.99045433e-01
  -9.44934312e-01   1.08325628e+00   6.61105015e-01   2.00477600e+00
  -6.21931145e-01  -4.42610971e-01   1.55512565e+00  -2.34940342e-01
  -1.80660809e+00  -1.04444304e+00  -1.31868758e-01  -1.06363300e+00
   1.99452744e-01  -5.16354069e-01  -7.17794152e-01  -4.93125496e-01
   7.30863966e-01  -3.97392088e-02  -8.07285651e-01   3.62122093e-02
   2.87212919e-01  -3.42469655e-01   1.69802421e-01  -1.20004403e+00
   1.80603626e-01  -2.77886862e-01   1.28721740e+00   2.20948332e-01
   1.07498766e+00   7.61151978e-01  -1.14772422e+00  -4.10463090e-01
   1.41356006e+00  -1.39942092e+00  -9.41572246e-01  -1.95133217e+00
   5.81492929e-01  -1.04094464e+00   6.84228468e-02  -9.13995323e-01
   1.27109853e+00   4.38327007e-01  -1.32285865e+00   4.19381907e-01
  -7.58841541e-01  -6.31574470e-01   6.11066011e-01  -1.31629659e+00
  -1.12835248e+00   4.44155282e-01   2.03144065e-01  -7.41700262e-01
  -1.69316743e+00  -5.31673350e-01   1.31910625e+00   2.39838001e-01
  -1.92196504e+00  -2.35910409e+00  -3.14737549e-01   1.87122084e+00
  -5.72948937e-02  -9.73029924e-01  -1.85086017e-01   1.10281882e+00
  -6.35725636e-01   1.42455637e+00   1.34764904e+00  -1.07256518e+00
  -1.27960447e-01  -3.85691781e-01   2.00689202e-01   5.30360204e-01
  -5.84661295e-02   5.59501939e-01   1.71760346e-01   3.21997310e-01
  -6.26850186e-01  -3.45337736e-01  -7.58126939e-01  -4.70202100e-02
   7.53177730e-01  -8.95705512e-01  -3.60578077e-01  -3.22447867e-01
   8.53417880e-02   5.78246985e-02   1.23998693e+00  -2.06135804e+00
  -1.55010859e+00  -5.80190662e-01   6.14965392e-01  -1.35307374e-01
   4.13915303e-01   2.13702341e+00  -1.59648172e+00  -5.82890700e-01
  -2.24164635e+00   1.33868557e+00  -6.80063461e-01   6.06424504e-01
  -3.45359685e-01   1.79120094e+00  -2.08000677e-01  -1.34963171e-01
   1.22894951e-01  -3.12841995e-02   2.14765435e+00   1.08531374e+00
   6.60823630e-01   1.13268657e+00  -5.77218119e-01  -6.23931014e-01
  -1.45827522e-01  -4.25202420e-02  -1.11147074e-01  -2.31996228e+00
   3.50586416e-01  -9.45301488e-01  -2.49751518e-01  -6.04875064e-01
  -1.29692191e+00  -4.98995297e-01  -1.63612825e-01  -8.66242100e-01
   1.66359117e+00   8.74237357e-01  -7.47449334e-01  -3.04265302e-01
   7.32225773e-02   7.00115207e-01  -7.41867476e-01  -1.41263745e-01
  -3.96428795e-01  -9.69013991e-01   3.29063868e-01  -1.27983409e-01
  -5.83896926e-01  -1.67696037e+00  -7.63331188e-01   1.18874060e-01
   8.48481636e-01  -1.63397512e+00  -8.19186762e-01  -3.32855184e-01
   1.16852423e+00   1.39392843e+00  -5.38042416e-01  -1.07356083e+00
  -1.14886365e+00   5.21086479e-01  -2.12885678e+00  -4.12849981e-01
   7.91528403e-01  -8.23170497e-01   1.87441979e+00  -1.36993704e+00
   4.76834305e-01  -3.26370969e-01  -3.67854749e-01  -7.26938399e-01
   8.09801448e-01   2.51198920e+00   3.56915229e-01  -2.54983451e-01
  -1.44171443e+00   4.60024234e-01   2.63685057e-01  -5.94825247e-01
  -9.84242448e-01  -4.97994804e-02   2.35485879e+00  -1.93900450e-01
  -6.53016302e-01   8.87233726e-01  -2.80150239e+00   9.05762006e-01
   4.08678677e-01  -6.81497772e-01   1.17321980e+00   8.19168014e-01
  -9.40753945e-01   1.60284566e+00  -1.52814257e+00  -1.27721158e+00
   6.81123620e-01   1.71472764e-01   7.65446208e-02  -9.53235010e-01
  -8.71715727e-01   3.15672308e-01   4.06340239e-01   9.83262660e-01
  -5.88616524e-01   1.21725679e+00  -7.86854635e-01  -1.28339733e+00
   8.22405802e-01  -5.28495086e-01  -8.00970266e-01  -1.75614387e+00
  -2.05405368e+00  -1.27566391e-01   2.36984594e-01  -9.01924717e-01
  -4.30518484e-01  -2.91784938e-01   6.54276508e-01  -9.51377711e-02
  -3.98568543e-01   6.60975269e-01  -1.41035066e+00   1.04368355e+00
  -4.11681786e-01  -1.43639928e+00   1.19240559e+00   4.95506687e-01
   1.57903626e+00  -2.79574310e-01   1.01947458e+00   2.78850004e-01
   5.75906939e-01  -3.66377257e-01   2.73652690e-01  -1.49404684e+00
  -7.01688625e-01  -1.10984134e+00   1.39178673e+00   6.08530075e-02
  -1.52852451e+00   7.87616683e-01  -6.44778055e-02  -3.45015409e-01
  -1.35845365e+00   2.12514267e-01  -7.25537962e-01  -4.56041157e-01
   1.89184553e-01  -1.51740959e+00   3.88779144e-01   1.41803863e+00
  -5.28035622e-01  -1.59373290e+00   1.36580111e+00  -3.97282334e-01
   5.31464020e-01   1.27122602e+00   1.96755527e-01   2.71952742e+00
  -5.46571756e-01   9.19082325e-01   8.29534026e-01   5.74567062e-01
  -1.66778996e+00   9.55768627e-01   9.18737099e-01   7.21911661e-01
  -4.55460821e-01  -6.22525028e-01  -1.53614233e+00  -1.11358462e+00
   3.67407611e-01   4.56837263e-01   5.05896647e-01  -3.82801562e-01
   1.89874303e+00   1.03889449e+00   4.23131063e-01  -8.11990652e-01
  -3.03945375e-02  -1.99952646e+00   4.46676030e-01   1.95721820e+00
   1.08873300e+00  -8.37814242e-03  -2.07361839e-01   3.60056055e-01
  -1.67468643e+00  -7.06049428e-01   9.65039385e-01   7.12357184e-01
   2.42874285e-01  -3.36268160e-01   5.06770876e-01   7.23624267e-01
   1.45186075e+00  -1.61955497e+00  -1.03411671e-01   5.59007686e-01
   5.20425979e-02   1.12046422e+00  -5.01872995e-01   1.06420808e+00
  -1.22003112e+00   7.89863655e-01  -5.85431950e-01  -2.24825569e+00
   4.22184348e-01  -2.06462345e-01  -2.26357030e+00   6.99624901e-01
   1.03726395e+00   4.53993244e-03   4.51026021e-01   1.29149321e-01
   4.79950159e-01   9.42474221e-01  -7.17021320e-01   8.08699977e-01
   1.29705448e+00  -7.29463525e-01   7.61622706e-01   1.33106913e+00
   1.26861623e+00   7.29526959e-01  -5.95726986e-01   5.03512395e-01
  -1.29418546e+00  -5.90249629e-01  -3.73911118e-01   1.89357574e+00
  -1.28354545e-01  -1.16177445e+00  -2.30600786e+00  -2.04517667e-01
  -3.76288692e-01   1.01317587e-01  -6.96236526e-01   1.90057656e-01
  -1.21489214e-01  -5.57939235e-01   1.87873908e+00   5.07359236e-01
   3.06841944e-01   6.65237276e-01   1.33444863e+00   1.44620678e+00
  -1.28550235e+00   1.13688401e+00  -5.34448063e-01  -1.32370673e-01
  -1.36045000e-01  -6.60853337e-01   2.81913662e-01  -1.03592075e+00
  -1.41860719e+00   1.81524496e+00   1.41390111e+00   4.23515124e-01
   1.02554970e+00   9.63881472e-01   7.12236894e-01  -3.71863624e-02
   1.76545480e+00   3.26268427e-01  -1.29671739e+00  -7.86767462e-01
  -1.23347607e+00   9.52597569e-01  -3.51206854e-01  -5.63237921e-01
  -5.69407473e-02   1.52866118e+00  -1.00496869e+00   1.63372644e+00
  -6.50946328e-01  -7.16273718e-01  -5.90126661e-01   1.30985556e+00
  -3.79768720e-01  -9.99915425e-01  -1.71145067e+00  -1.57329194e+00
  -7.30696372e-02   7.62151506e-01   1.48689646e+00   4.31745960e-01
   2.63882808e-01   1.72206428e-01   2.00851731e-01  -1.18314347e+00
  -1.50402279e+00   1.95153249e+00  -2.78003872e-01  -1.93854878e-03
   1.16582803e+00  -3.20238544e-01   9.36585278e-01   1.48863154e+00
  -1.39264019e+00   1.55666358e+00  -4.20662970e-01  -5.50893690e-02
  -7.10504077e-01  -8.73934840e-01  -6.26132410e-01  -1.29781690e-01
  -2.22777306e+00   3.73831654e-01  -3.82432449e-01   5.37121996e-02
  -9.32183728e-01   3.53740442e+00  -1.05598008e+00  -7.47679114e-03
  -1.06551955e+00  -3.46521140e-01   3.88344981e-01  -3.64222377e-01
  -1.37048480e+00   1.01306470e-01   1.64514228e+00  -3.61445671e-01
   2.31614366e-01  -1.25968487e+00  -6.62537224e-01  -5.61715767e-01
  -1.13596993e+00  -1.36075508e+00   7.31594663e-01   2.67213081e-01
   7.13141895e-02  -2.46130956e-01  -1.39377306e+00   1.29318009e-01
  -1.92097772e-01  -1.15100304e+00   1.02016481e-01  -1.93068490e+00
   9.69566668e-01   9.75694115e-01  -7.57001670e-01  -7.60531329e-01
  -3.03552280e-01  -3.67388243e-01   7.60265551e-02   4.42746285e-01
   1.92679575e-01   1.04068792e+00  -1.11139888e+00  -1.72269547e+00
  -2.30986045e-01  -2.83927273e-01  -7.21803150e-01  -1.37631072e+00
  -1.10935589e+00   8.09954622e-01   5.83783993e-01  -7.31810290e-01
   3.06739636e-01   1.27915913e+00   2.73785857e-01   9.65830330e-01
  -6.15660861e-01   1.64785898e-01  -4.61886474e-02   3.25744723e-01
   3.67606054e-01  -6.68566804e-01  -3.20837511e-01  -2.09469223e-01
  -6.60728681e-01  -3.49847761e-01   7.84861730e-01   1.26622383e+00
  -1.80302074e+00   1.53355423e+00  -1.22745466e+00  -7.02346088e-01
   3.17239216e-01  -1.25357347e+00   2.39430740e-01   7.85747146e-01
   3.81228699e-01   5.72990733e-01  -4.98818013e-01   1.42759290e+00
   1.69600421e-01  -3.16755901e-01  -4.16914206e-01   1.22881192e+00
  -5.58348844e-02  -6.26676967e-01   7.38764150e-02   3.32789760e-01
   1.67321192e+00   1.74966096e-01   2.65341416e-01  -3.29541011e-01
  -1.86929223e+00  -1.14107067e+00  -7.53407493e-01   6.99519962e-01
   4.29125580e-01   2.19232831e-01   7.42525121e-01  -5.42444700e-01]
0.0
Final model:
[  3.77556006e+00   1.67658791e+00  -3.29136306e+00  -1.30945858e-01
  -1.18248243e+00   1.78664688e+00  -1.25736948e+00  -2.02258094e-01
  -9.98941721e-01  -7.63579832e-01   1.57014674e+00  -3.94171743e-01
   1.15549836e+00   3.86136052e-01   1.49428543e-01  -3.36443795e-01
  -1.06508201e+00   3.64432700e-01  -5.03306943e-01  -2.56198989e+00
   5.18392678e-02   1.94736399e+00  -5.96506785e-01   5.43487120e-01
  -4.42454966e+00   1.09235863e+00   1.58964361e+00  -6.50332550e-01
   1.11307437e+00   4.90929707e-01  -1.11052327e+00   1.28429702e+00
  -2.09541131e-01   3.77961492e+00   7.58902758e-01  -9.27180732e-01
  -8.22309817e-01   3.13218677e+00   8.31525921e-01  -1.47034985e+00
  -4.92172262e-01   1.01060797e+00  -4.05332960e-01  -1.41643633e+00
   2.69651131e+00  -3.55972129e-01  -2.48351450e+00   2.49434037e-01
  -2.67932735e-01   9.79001494e-01  -9.92249171e-01  -5.96254783e-01
  -7.99429663e-01  -7.35963162e-02  -2.40140608e-01  -2.12372700e+00
  -1.44504208e+00   3.75125978e+00  -1.26714397e+00  -3.18828201e-01
  -4.49929404e+00  -1.28865508e+00  -1.28861373e+00  -1.32248557e+00
  -1.57317700e+00  -2.12905251e+00  -2.64333717e-01  -2.01285336e+00
  -1.61006020e-01   2.27637572e+00   1.27586244e+00  -3.91120851e-01
  -1.36317244e+00  -3.82788078e-01  -1.92953169e-01  -2.88078903e-01
  -3.88464164e-01  -1.81768162e+00   1.69723438e+00   4.78324656e-01
  -1.57736636e-01  -2.51055129e+00  -7.41662774e-01  -8.76968592e-01
   4.73273895e-02  -3.97794058e+00   3.18338634e-01   7.68762163e-01
  -1.67187179e+00   6.14911653e-01   1.10559059e+00   2.47633722e+00
  -4.73221220e-01   1.35373891e+00  -6.50581076e-02   9.55459296e-02
   2.03508078e-01  -2.53910629e+00   6.10795696e-01   1.38991127e+00
   9.15067374e-01   3.63101705e+00   8.26695293e-01  -6.68788818e-01
  -1.74436254e+00  -3.33700192e+00  -1.53928722e+00  -2.96420698e-01
   1.57903269e+00   1.79884911e+00  -2.50262001e-01  -1.08973660e+00
  -1.85498194e+00  -1.19668614e+00  -4.26989972e-01  -7.91529732e-01
   1.31224993e+00   3.63791777e+00  -2.94546634e-01  -2.45556566e+00
  -1.04590869e+00   6.08167811e-01   2.77936798e+00  -9.12923966e-01
  -7.83906049e-01   3.56631289e-01  -3.70281016e-01   4.54090725e-01
   1.60936498e+00   2.15636982e+00  -1.89672027e+00   6.21229996e-02
  -1.50861299e+00   1.02692098e+00   1.04198787e+00  -7.97214484e-01
   1.14241240e+00   4.56888918e-01   7.37187678e-01   2.60242725e+00
   3.95684108e+00  -2.15633608e+00   4.70875935e-02   1.83190139e+00
   1.36581873e+00  -3.86096041e-01   3.01203822e+00   2.76356831e+00
  -1.98179848e+00   6.40089101e-01   2.20144846e+00  -1.06899827e+00
   3.09186005e+00  -7.11234569e-01  -3.78276117e+00   1.83545092e-01
  -1.21826498e+00  -5.10159929e-01   3.25947346e+00  -1.58883034e+00
  -2.49492079e+00   1.15062631e+00   2.46314738e+00   1.83675178e+00
   1.90239257e+00  -2.91436561e+00  -3.27309202e+00  -1.15415531e+00
   2.15120857e+00   2.61265590e+00   1.04777203e+00  -6.99138307e-01
   1.87069022e+00   4.36184386e-01   3.33629754e-01  -1.50581420e+00
   2.94741903e+00  -2.34902225e+00   8.38337601e-01   5.00249453e+00
  -2.38272140e-01  -3.01384569e-01   1.89878373e-01  -2.19144302e+00
   3.58792779e+00  -3.52910136e+00   1.26603845e+00   7.75531918e-01
   2.00485484e+00  -3.53271420e+00   8.10974186e-01   1.84989351e+00
   2.51069243e+00  -4.26123595e-01  -8.50810390e-01  -4.10358099e-01
  -9.42958105e-01  -4.61438083e-01  -4.31821036e+00  -6.02840390e-01
  -2.51884238e+00  -7.50987410e-01  -1.77051049e-01  -6.47304330e-02
  -1.97060390e+00   6.03047738e-01   1.55041662e+00  -2.68797620e+00
  -5.15198874e-01   3.97884570e-01   2.98066051e+00   4.96007616e+00
   1.80033955e+00  -1.29148046e-01  -1.65518882e+00  -3.11845247e+00
  -1.11648300e+00   1.73962809e+00  -1.99331693e+00  -4.71585789e-01
  -1.02846099e+00  -1.25604064e+00  -2.03443821e+00  -1.34359639e+00
  -1.47179647e+00   2.02529621e-01   3.97979332e-03   9.49734864e-01
  -9.27013373e-01  -2.79633499e+00   1.74473437e+00  -1.71349731e+00
   2.74689851e-01   2.39719164e+00   8.98615247e-01   9.54759734e-01
  -7.86490056e-01   2.61971730e+00  -5.77213868e-01   1.10154531e+00
  -2.31061307e+00   2.45251463e+00  -2.58008735e+00  -6.03860022e-01
  -3.52720049e-01   2.74259886e+00  -6.69666767e-01   5.29239930e-01
  -6.22485624e-02   2.60974732e-01  -8.06554943e-01  -4.04441241e-02
   3.12796271e-01   3.59310410e-01   2.50710253e+00  -1.13593713e+00
   2.29305634e+00  -6.16724464e-01  -1.81585729e+00  -2.30747339e-01
  -8.54934300e-01  -2.53499110e-01   2.49450311e+00   1.20811077e+00
  -2.54033571e+00  -1.36962728e+00  -3.09588677e+00   2.45733519e+00
   1.10845535e+00   1.30030015e+00   6.36532470e-01   1.14534142e+00
  -1.35369898e+00  -8.31991520e-01  -2.85235698e+00  -3.19790257e+00
  -1.08782735e+00  -5.41336738e-01  -2.87775857e+00  -1.04789255e+00
   1.62996278e+00   5.17471118e-01  -4.63886043e-01  -1.50424824e+00
  -4.18059294e-01   1.64772733e+00  -2.47286106e+00   3.92516794e+00
   3.70956959e+00   3.31061339e-01   3.21938139e-01  -1.98907732e+00
  -1.62564497e-01   9.25237287e-01  -2.39858458e+00  -1.47865761e+00
  -1.07760350e+00  -2.12954216e-01   3.56462897e+00  -1.15009932e+00
  -9.37380338e-01   3.83415308e+00  -9.34323105e-02   7.63654725e-02
   1.11119989e+00   1.07718664e+00  -1.03859306e+00  -1.57055699e+00
  -4.46823604e-01  -1.06687172e+00   8.21890243e-01  -4.42575478e-03
  -4.48360060e+00   7.43566499e-01  -6.17821929e-01   6.03409669e-01
  -3.66864713e-01   4.48455869e-03   2.39363432e+00   2.33992814e-01
   4.21215749e+00   4.55885777e+00   3.92714524e-01  -8.59827680e-01
   6.31566745e-01  -3.34226322e+00  -1.03161005e-01  -2.74638014e+00
   1.23743057e+00  -7.47138590e-01   1.12158829e+00   2.10546219e+00
   4.10013374e+00  -9.54526634e-01   8.29485521e-01  -3.39853916e+00
  -1.70827359e+00  -4.33272812e-01  -8.00704300e-01   3.32479096e+00
  -2.64906932e+00   5.94967233e-01   1.81417239e+00  -2.83409352e+00
  -2.76620738e-01   1.19838737e+00  -7.31473471e-01  -2.80155572e-01
   2.15023089e+00  -1.17006423e+00   1.14970105e+00   2.51563880e-01
  -4.38913007e-01   2.04797482e+00  -7.94636949e-02   3.43728270e-01
  -6.02989618e-02   1.22478848e+00  -1.16759381e+00   3.23341798e-01
   1.53320439e+00  -4.26366857e-01   1.04518067e+00  -6.36737194e-01
   3.25811916e-01  -3.30528851e+00   4.11734855e-01   2.50890638e+00
  -2.21787212e+00  -7.17371728e-01   3.42200962e-01  -1.17867176e+00
   1.20111575e+00  -1.01109980e+00   1.20831464e+00   1.24711079e-03
  -1.71053384e+00   1.73891705e+00   4.44005955e-01  -1.77678764e+00
   3.66508508e+00  -7.04048127e-01  -2.68029967e+00   1.75221936e-01
   6.00537286e-01  -2.73095435e+00   2.73826244e+00  -9.13690757e-01
   1.26268560e+00  -1.97674163e-01   5.21461602e-02   7.08639718e-01
   1.46285279e+00   3.32089390e+00  -1.41666032e+00   5.40736688e-01
  -2.05106750e+00  -2.34978746e+00   1.55974139e+00  -1.91425957e+00
  -2.10754689e+00   7.22699816e-01   1.29999777e+00  -1.65573513e-01
   3.70725978e+00   3.55281708e+00  -2.01923747e+00  -2.89179581e+00
   2.59579559e+00  -1.75259386e+00  -1.08206430e-01   7.45312239e-03
   5.85764039e-02   6.14869059e-01  -1.29978952e-02  -2.92607897e+00
   8.14497304e-01   1.32776160e+00  -2.09535483e+00  -6.23154792e-01
  -1.02933599e+00  -6.79472063e-01   1.36643972e-01   1.30739153e+00
  -7.70938238e-02   1.00573421e+00   9.51384392e-01  -1.53817711e+00
  -2.35552187e+00  -7.56852770e-01  -1.26638116e+00  -2.36581588e+00
   1.87100888e-01  -1.80800625e+00   2.28044384e+00   1.43620211e-01
   1.41963125e+00   8.60386903e-01   2.23154240e+00   2.68496183e+00
  -1.64176977e+00   1.45565317e+00   3.93317919e-01   4.63312083e+00
  -2.22928479e-01   1.43214324e+00  -6.36949102e-01  -7.63117543e-01
   5.71798558e-01  -9.99673642e-01  -1.85985574e+00  -3.71322109e+00
  -6.49851612e-01  -1.56570988e+00   2.12323186e+00  -1.20302700e-02
  -1.78012832e+00   7.68680703e-01  -1.32861308e+00   2.48974714e+00
  -3.93162518e-01   1.43291241e+00  -1.89201627e-01  -1.03040512e+00
   2.59071790e+00  -2.12592921e-01   1.61041461e+00   1.17219257e+00
  -1.46006051e+00  -6.91233414e-01  -6.22940430e-01  -8.66083482e-01
   7.55935471e-01   1.37146465e+00  -2.16174006e+00  -2.81248199e+00
   6.54337758e-02   9.78157136e-01  -3.96411245e-01   1.96117462e-01
   3.04485846e+00   1.73779365e+00   6.50505894e-01  -2.92031943e+00
  -1.02245870e+00   2.53395754e+00  -1.21661383e+00  -1.81445843e+00
  -2.41390310e-01   2.29891887e+00   7.72474664e-02   1.66204376e-01
  -1.32964716e+00   1.09600867e+00   1.63308666e+00   9.43247657e-01
   1.28885703e+00   2.82512813e+00  -5.84870327e-01  -1.91773348e+00
  -1.49308557e+00  -8.16129941e-01  -7.46493838e-01  -1.53153033e+00
   1.81696140e-01  -5.07009044e-02   1.99327335e+00   2.58423667e+00
  -5.77897925e-01   3.56759945e+00   3.98155760e-01  -5.07832445e-01
   1.81953549e-01   1.85387295e+00   1.38552765e+00   2.46977981e-01
   2.36358716e-01  -7.73858776e-01   1.07066315e+00   2.01714544e+00
  -6.71629691e-01  -3.88064770e+00  -1.44425362e-01  -1.81274366e+00
  -1.46691251e+00  -3.43959920e+00   2.80986444e+00  -1.17730142e+00
  -2.15336320e+00   1.55093198e+00   1.15522980e+00   1.31141415e+00
  -1.48348157e+00  -1.44346424e+00  -2.68688220e+00   4.84821223e-01
   1.53581802e+00   2.26427221e+00   9.72914780e-02  -1.42469335e+00
  -1.82068986e+00  -1.86450827e+00   3.02468900e-01  -4.79596616e-01
   9.47111398e-01   1.14024438e+00  -8.47782688e-01  -1.14826786e+00
  -2.70538807e-01   2.64888030e-01   9.43489510e-01   4.65450048e-01
  -1.77615395e+00   1.50161604e-01   1.03315021e+00  -2.39435007e+00
   1.40142814e+00   8.15688739e-01  -1.08414515e+00   7.60111618e-01
   1.41953427e+00  -1.25001136e+00   1.30342412e-02  -1.60291024e+00
   2.15567458e+00  -1.58323187e+00   3.25596645e-01   7.01505666e-01
   8.59101035e-01   1.80383869e+00   2.92749652e+00   6.09916419e-01
  -2.71513471e+00   5.60586660e-01   2.25575384e-01   6.28246042e-01
   9.32571286e-01   2.93902667e+00   1.72953108e+00  -2.95012867e-01
   2.24469224e+00  -2.71902436e-01  -1.10800541e+00   7.32789542e-01
   7.02507751e-01   7.39097377e-01   4.14989958e+00   2.67674147e+00
   4.14115952e-01  -3.84743988e-01  -1.27175024e+00  -2.19024992e+00
  -2.18140078e+00  -1.05260452e+00  -3.49783833e+00   5.57055340e-01
  -1.67155251e+00  -3.34448800e-01   6.28179616e-01  -2.74330916e-02
  -2.06374038e+00   6.35278593e-01   1.18046664e+00   5.60642386e-01
   1.35393263e+00  -2.61514826e+00   1.77800056e+00   1.09163009e+00
   2.05349909e+00   1.89596602e+00  -9.02341407e-01  -5.93588299e-01
   3.01763579e+00  -4.37601720e+00  -2.27288981e+00  -5.79656815e+00
   7.64501306e-02   6.17753732e-01  -4.78503966e-01   2.11041028e+00
   9.12261551e-01  -6.46346074e-01   1.28761655e+00  -1.81290786e+00
  -3.46692744e-01  -1.12915892e+00   1.14566399e+00  -8.74047326e-02
   4.75776147e-01   2.89152270e+00  -1.12603767e+00   1.81324941e+00
  -2.37493794e-01   1.22815705e+00   5.48781676e-01  -8.47839017e-01
  -6.79308765e-01   1.84570285e+00  -4.84912772e-01   1.49939217e+00
  -1.42955765e+00   8.99106720e-01   2.09687546e+00   2.32900062e-01
  -2.16848707e-01  -4.04704275e+00   2.35142329e+00  -1.18962640e+00
   8.94147748e-01   2.93820180e+00   2.07633335e+00   8.07215486e-01
  -4.65727597e-01  -1.03032045e+00   2.74614894e+00  -8.16323188e-01
   3.51113518e+00   2.05167186e+00  -1.55314194e+00  -9.40403259e-01
  -4.09650743e-02  -2.36659383e+00   5.53580697e-01   1.09338078e+00
  -3.87102251e-01  -1.39174144e+00  -1.80277537e+00  -2.10924254e+00
   1.54615036e+00   8.56128405e-01  -1.33996909e+00  -1.71725473e+00
  -1.57777016e+00   2.65809549e+00   8.54861156e-01   2.19105117e+00
  -1.94545749e+00  -1.31990906e+00  -2.48547338e+00  -1.21667731e+00
   1.12321564e-01   7.03858613e-01  -1.31534567e+00  -3.01909282e+00
   1.76111359e+00   1.03464584e+00   1.64381967e+00   2.16596567e+00
  -1.93396820e+00   3.86992118e+00   2.65106268e-01  -1.41634383e+00
  -1.71242017e+00   6.57565704e-01  -2.79514824e+00   1.10829177e-01
   2.06668705e+00   4.07656873e-01   1.97741136e+00  -9.54814686e-01
   1.37971688e+00   1.06623143e+00  -2.89450811e+00  -2.35645806e+00
   3.92589659e+00   2.83499175e+00   9.80082282e-01   1.68067343e+00
  -1.06050525e-01  -1.74825733e+00  -1.51233969e+00   2.79399091e-01
   1.61218267e+00   6.79050628e-01  -7.60576055e-01   1.07710598e+00
   8.68093576e-01   7.96236864e-01  -3.37133751e+00   1.90457592e+00
  -3.59889053e-01  -1.55153960e+00   1.28538465e+00   2.24184560e+00
   1.50305971e+00  -1.88109843e+00   6.17140621e-02  -9.47180228e-01
  -2.18103376e+00  -1.24519494e+00   3.45694238e+00   2.55172108e+00
   1.40859574e-02  -9.53010206e-01  -1.92951266e+00  -2.29692512e+00
   3.45593292e+00   1.08317253e+00  -3.38025123e+00   1.59851381e+00
   1.74317752e+00  -2.45672994e+00  -1.78299865e+00  -2.62810516e+00
   1.46304208e+00  -8.25624629e-01  -3.11702318e+00   2.20752930e+00
  -2.64307436e+00  -2.50232719e+00   1.34700372e+00  -3.25327855e+00
  -1.38651861e+00   6.84753858e-01  -4.39893291e-01   2.41137859e+00
   4.02819568e+00   2.59001772e+00  -4.89821103e-01   3.59276693e+00
  -1.69033166e+00  -9.00221564e-01   7.14153357e-01   3.66780469e-01
   3.59698167e+00  -1.12640871e+00   3.03631173e+00   1.48208320e+00
   1.27545162e+00   3.33220799e+00   2.06530977e-01   3.00507382e+00
   1.55880239e+00   1.20344288e+00   2.17279620e+00  -1.49944096e+00
  -1.25077507e+00  -2.21842681e+00   1.11603549e+00   4.58609893e-01
  -4.98923219e-01   1.70758638e-01   1.72125168e+00   5.63973756e-03]
497.500000002
target values for D:
[0 1 1 1 1 1 1 0 1 0 1 0 0 1 0 1 0 1 0 1 1 0 1 1 0 1 1 1 0 0 1 1 1 0 1 1 0
 1 0 1 1 0 0 0 1 0 1 1 1 0 1 1 0 1 0 0 1 1 1 1 1 1 0 0 0 1 0 0 1 1 1 1 0 0
 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 1 0 1 0 1 0 0 1 1 0 1 1 1 0 1 1 1 1 0 1
 0 1 1 1 0 0 1 0 1 1 0 0 1 0 1 1 1 1 0 1 0 1 1 1 0 0 0 1 1 0 0 0 0 0 1 0 0
 1 1 1 1 0 1 0 1 0 0 1 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 1 1 1 0 1 0 0 1 1
 0 1 0 1 0 1 1 1 0 1 1 0 0 1 0 0 0 0 0 1 0 1 1 0 0 0 0 1 1 0 1 1 1 1 0 0 0
 0 0 1 1 0 1 0 0 1 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 0 1 0 1 1 0 0 0 0 0
 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 1 0 0 0 1 0 1 0 1 0 0 0 1 1 1
 0 0 1 0 1 0 0 1 0 1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 1 0 1 1 0 0 1 1 1 1 1 0 1
 0 0 0 1 1 1 1 0 1 0 0 1 0 1 0 0 1 0 1 0 1 1 1 0 1 1 1 0 1 0 0 1 0 0 0 1 0
 1 0 1 1 0 1 1 1 1 1 1 0 0 1 1 1 1 0 1 0 1 0 1 0 0 1 0 0 1 0]
prediction on D:
[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]

Derivatives in Theano

Computing Gradients


In [1]:
import numpy
import theano
import theano.tensor as T
from theano import pp

For this

$ \begin{gathered} \frac{d (x^2) }{ dx} = 2 \cdot x \end{gathered} $


In [2]:
x = T.dscalar('x')
y = x ** 2
gy = T.grad(y,x)
pp(gy) # print out the gradient prior to optimization


Out[2]:
'((fill((x ** TensorConstant{2}), TensorConstant{1.0}) * TensorConstant{2}) * (x ** (TensorConstant{2} - TensorConstant{1})))'

fill((x ** 2), 1.0) means to make a matrix of the same shape as x ** 2 and fill it with 1.0


In [3]:
f = theano.function([x],gy)
f(4)


Out[3]:
array(8.0)

In [4]:
numpy.allclose(f(94.2), 188.4)


Out[4]:
True

A plot of the gradient of the logistic function, with x on the x-axis and $ds(x)/dx$ on the y-axis


In [6]:
x = T.dmatrix('x')
s = T.sum(1 / (1 + T.exp(-x)))
gs = T.grad(s, x)
dlogistic = theano.function([x], gs)
dlogistic([[0, 1], [-1, -2]])


Out[6]:
array([[ 0.25      ,  0.19661193],
       [ 0.19661193,  0.10499359]])

Computing the Jacobian


In [7]:
x = T.dvector('x')
y = x ** 2
J, updates = theano.scan(lambda i, y, x: T.grad(y[i], x), sequences=T.arange(y.shape[0]), non_sequences=[y,x] )
f = theano.function([x], J, updates=updates)
f([4, 4])


Out[7]:
array([[ 8.,  0.],
       [ 0.,  8.]])

Computing the Hessian


In [8]:
x = T.dvector('x')
y = x ** 2
cost = y.sum()
gy = T.grad(cost, x)
H, updates = theano.scan(lambda i, gy, x : T.grad(gy[i], x), sequences=T.arange(gy.shape[0]), non_sequences=[gy, x] )
f = theano.function([x], H, updates=updates)
f([4,4])


Out[8]:
array([[ 2.,  0.],
       [ 0.,  2.]])

R-operator


In [9]:
W = T.dmatrix('W')
V = T.dmatrix('V')
x = T.dvector('x')
y = T.dot(x,W)
JV = T.Rop(y, W, V)
f = theano.function([W,V,x],JV)
f([[1,1], [1,1]],[[2,2],[2,2]],[0,1])


Out[9]:
array([ 2.,  2.])

L-operator


In [10]:
W = T.dmatrix('W')
v = T.dvector('v')
x = T.dvector('x')
y = T.dot(x,W)
VJ = T.Lop(y,W,v)
f = theano.function([v,x],VJ)
f([2,2],[0,1])


Out[10]:
array([[ 0.,  0.],
       [ 2.,  2.]])

Hessian times a Vector


In [11]:
x = T.dvector('x')
v = T.dvector('v')
y = T.sum(x ** 2)
gy = T.grad(y, x)
vH = T.grad(T.sum( gy * v), x)
f= theano.function([x,v], vH)
f([4,4], [2,2])


Out[11]:
array([ 4.,  4.])

or, making use of the R-operator:


In [12]:
x = T.dvector('x')
v = T.dvector('v')
y = T.sum( x ** 2)
gy = T.grad(y,x)
Hv = T.Rop(gy, x, v)
f = theano.function([x,v], Hv)
f([4,4],[2,2])


Out[12]:
array([ 4.,  4.])

Conditions


In [13]:
from theano import tensor as T
from theano.ifelse import ifelse
import theano, time, numpy

a,b = T.scalars('a','b')
x,y = T.matrices('x','y')

In [14]:
z_switch = T.switch(T.lt(a,b), T.mean(x), T.mean(y))
z_lazy   = ifelse(T.lt(a, b), T.mean(x), T.mean(y))

In [3]:
rng = numpy.random

In [4]:
N=400
feats=784

In [5]:
D = (rng.randn(N, feats).astype(theano.config.floatX), rng.randint(size=N,low=0, high=2).astype(theano.config.floatX))

In [11]:
print(D[0].shape); D[0]


(400, 784)
Out[11]:
array([[-0.51540321,  0.45720622, -1.07207155, ..., -0.14092861,
        -1.29015946,  0.79136688],
       [ 0.41795   ,  0.58462614, -0.13168216, ...,  0.65973812,
        -0.46558085, -0.15008399],
       [ 0.09026403, -1.53019106, -0.43091467, ...,  0.42072779,
         0.85275686, -2.31214309],
       ..., 
       [ 1.0708015 , -1.14941835,  2.08980322, ..., -0.02039757,
        -0.42613024, -1.25479817],
       [-1.50268137, -0.1785637 , -1.19366574, ...,  0.53540504,
        -0.3872523 ,  1.20521164],
       [ 0.76399541,  0.65020931, -1.35631859, ...,  0.5690015 ,
         0.69514877,  0.35484901]], dtype=float32)

In [10]:
print(D[1].shape); D[1]


(400,)
Out[10]:
array([ 0.,  0.,  1.,  0.,  0.,  0.,  1.,  1.,  0.,  1.,  0.,  0.,  1.,
        0.,  1.,  1.,  0.,  0.,  1.,  1.,  1.,  0.,  0.,  0.,  0.,  1.,
        0.,  0.,  0.,  1.,  1.,  0.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,
        1.,  1.,  1.,  0.,  0.,  0.,  1.,  1.,  0.,  1.,  1.,  1.,  1.,
        1.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,  1.,
        0.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  1.,  1.,  1.,  0.,  0.,
        1.,  1.,  1.,  0.,  1.,  1.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,
        1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  1.,  0.,  0.,  1.,  1.,  1.,  1.,  0.,  1.,  0.,  1.,  0.,
        0.,  1.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  1.,  0.,  0.,  1.,
        1.,  0.,  0.,  1.,  1.,  1.,  0.,  1.,  1.,  1.,  0.,  0.,  0.,
        0.,  0.,  1.,  1.,  0.,  1.,  0.,  1.,  1.,  0.,  1.,  0.,  1.,
        1.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  1.,  0.,  0.,  1.,
        0.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,  1.,  0.,  1.,  1.,  0.,
        1.,  1.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  1.,  1.,
        1.,  0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,  1.,  0.,
        1.,  0.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  1.,  0.,  1.,
        1.,  1.,  1.,  0.,  0.,  1.,  1.,  0.,  0.,  0.,  1.,  1.,  1.,
        1.,  1.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  1.,
        0.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,
        1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,
        1.,  1.,  1.,  1.,  0.,  0.,  1.,  0.,  1.,  1.,  0.,  0.,  0.,
        0.,  1.,  0.,  1.,  0.,  0.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  1.,  0.,  1.,  1.,  1.,  1.,  0.,  1.,  1.,  0.,
        0.,  1.,  0.,  0.,  1.,  1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.,
        0.,  0.,  1.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,
        1.,  0.,  0.,  0.,  1.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,
        0.,  1.,  1.,  0.,  0.,  0.,  1.,  1.,  1.,  0.,  1.,  0.,  1.,
        0.,  1.,  1.,  0.,  0.,  1.,  1.,  1.,  1.,  0.,  1.,  1.,  0.,
        1.,  0.,  0.,  0.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,  1.,  0.,
        1.,  1.,  0.,  1.,  1.,  1.,  1.,  0.,  1.,  0.], dtype=float32)

In [8]:
training_steps = 10000

In [12]:
# Declare Theano symbolic variables
x = theano.shared(D[0], name="x")
y = theano.shared(D[1], name="y")
w = theano.shared(rng.randn(feats).astype(theano.config.floatX),name="w")
b = theano.shared(numpy.asarray(0., dtype=theano.config.floatX),name="b")

In [13]:
# Setting the tag.test_value attribute gives the variable its test value, i.e.
# provide Theano with a default test-value
x.tag.test_value = D[0]
y.tag.text_value = D[1]

In [17]:
# Construct Theano expression graph
p_1 = 1 / ( 1+ T.exp( - T.dot( x,w) - b)) # Probabilty of having a 1
prediction = p_1 > 0.5  # the prediction that is done: 0 or 1
xent = -y * T.log(p_1) - (1-y) * T.log(1-p_1) # Cross-entropy
cost = T.cast( xent.mean(), theano.config.floatX) +  0.01 * (w ** 2).sum()  # the cost to optimize
gw, gb = T.grad(cost, [w,b])

In [18]:
# compile exprression to functions
train = theano.function( 
                        inputs=[],
                        outputs=[prediction, xent],
                        updates=[(w, w - 0.01 * gw), (b, b - 0.01 * gb)],
                        name="train")
predict = theano.function(inputs=[], outputs=prediction, name="predict")

In [20]:
if any([n.op.__class__.__name__ in ['Gemv', 'CGemv', 'Gemm', 'CGemm'] for n in train.maker.fgraph.toposort()]):
    print('Used the cpu')

In [22]:
# elif
if any([n.op.__class__.__name__ in ['GpuGemm', 'GpuGemv'] for n in train.maker.fgraph.toposort()]):
    print('Used the gpu')


Used the gpu

In [23]:
for i in range(training_steps):
    pred, err = train()

In [29]:
print("target values for D")
print(D[1])


target values for D
[ 0.  0.  1.  0.  0.  0.  1.  1.  0.  1.  0.  0.  1.  0.  1.  1.  0.  0.
  1.  1.  1.  0.  0.  0.  0.  1.  0.  0.  0.  1.  1.  0.  1.  0.  1.  1.
  1.  0.  1.  1.  1.  1.  0.  0.  0.  1.  1.  0.  1.  1.  1.  1.  1.  1.
  0.  1.  1.  1.  0.  1.  1.  1.  0.  1.  1.  0.  1.  1.  1.  1.  1.  0.
  0.  1.  1.  1.  0.  0.  1.  1.  1.  0.  1.  1.  0.  1.  1.  1.  1.  1.
  0.  1.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.
  1.  1.  1.  1.  0.  1.  0.  1.  0.  0.  1.  1.  0.  1.  0.  1.  0.  1.
  1.  0.  0.  1.  1.  0.  0.  1.  1.  1.  0.  1.  1.  1.  0.  0.  0.  0.
  0.  1.  1.  0.  1.  0.  1.  1.  0.  1.  0.  1.  1.  0.  1.  1.  1.  1.
  1.  0.  0.  1.  0.  0.  1.  0.  1.  0.  1.  1.  1.  0.  1.  1.  0.  1.
  1.  0.  1.  1.  0.  0.  0.  1.  1.  0.  0.  0.  0.  1.  1.  1.  0.  0.
  0.  0.  1.  1.  1.  1.  1.  0.  1.  0.  1.  0.  0.  1.  0.  1.  0.  1.
  0.  1.  1.  0.  1.  1.  1.  1.  0.  0.  1.  1.  0.  0.  0.  1.  1.  1.
  1.  1.  0.  0.  1.  0.  0.  0.  0.  1.  0.  0.  1.  0.  1.  0.  1.  1.
  1.  0.  1.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  1.
  1.  0.  0.  1.  1.  1.  1.  0.  0.  1.  0.  1.  1.  0.  0.  0.  0.  1.
  0.  1.  0.  0.  1.  0.  1.  0.  0.  0.  0.  0.  0.  0.  1.  0.  1.  1.
  1.  1.  0.  1.  1.  0.  0.  1.  0.  0.  1.  1.  0.  0.  1.  0.  0.  1.
  0.  0.  0.  1.  1.  0.  1.  0.  0.  0.  0.  0.  0.  1.  1.  0.  0.  0.
  1.  0.  1.  1.  0.  0.  0.  0.  0.  0.  1.  1.  0.  0.  0.  1.  1.  1.
  0.  1.  0.  1.  0.  1.  1.  0.  0.  1.  1.  1.  1.  0.  1.  1.  0.  1.
  0.  0.  0.  1.  0.  1.  0.  0.  0.  0.  1.  0.  1.  1.  0.  1.  1.  1.
  1.  0.  1.  0.]

In [32]:
print("prediction on D")
print(predict().astype(theano.config.floatX))


prediction on D
[ 0.  0.  1.  0.  0.  0.  1.  1.  0.  1.  0.  0.  1.  0.  1.  1.  0.  0.
  1.  1.  1.  0.  0.  0.  0.  1.  0.  0.  0.  1.  1.  0.  1.  0.  1.  1.
  1.  0.  1.  1.  1.  1.  0.  0.  0.  1.  1.  0.  1.  1.  1.  1.  1.  1.
  0.  1.  1.  1.  0.  1.  1.  1.  0.  1.  1.  0.  1.  1.  1.  1.  1.  0.
  0.  1.  1.  1.  0.  0.  1.  1.  1.  0.  1.  1.  0.  1.  1.  1.  1.  1.
  0.  1.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.
  1.  1.  1.  1.  0.  1.  0.  1.  0.  0.  1.  1.  0.  1.  0.  1.  0.  1.
  1.  0.  0.  1.  1.  0.  0.  1.  1.  1.  0.  1.  1.  1.  0.  0.  0.  0.
  0.  1.  1.  0.  1.  0.  1.  1.  0.  1.  0.  1.  1.  0.  1.  1.  1.  1.
  1.  0.  0.  1.  0.  0.  1.  0.  1.  0.  1.  1.  1.  0.  1.  1.  0.  1.
  1.  0.  1.  1.  0.  0.  0.  1.  1.  0.  0.  0.  0.  1.  1.  1.  0.  0.
  0.  0.  1.  1.  1.  1.  1.  0.  1.  0.  1.  0.  0.  1.  0.  1.  0.  1.
  0.  1.  1.  0.  1.  1.  1.  1.  0.  0.  1.  1.  0.  0.  0.  1.  1.  1.
  1.  1.  0.  0.  1.  0.  0.  0.  0.  1.  0.  0.  1.  0.  1.  0.  1.  1.
  1.  0.  1.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  1.
  1.  0.  0.  1.  1.  1.  1.  0.  0.  1.  0.  1.  1.  0.  0.  0.  0.  1.
  0.  1.  0.  0.  1.  0.  1.  0.  0.  0.  0.  0.  0.  0.  1.  0.  1.  1.
  1.  1.  0.  1.  1.  0.  0.  1.  0.  0.  1.  1.  0.  0.  1.  0.  0.  1.
  0.  0.  0.  1.  1.  0.  1.  0.  0.  0.  0.  0.  0.  1.  1.  0.  0.  0.
  1.  0.  1.  1.  0.  0.  0.  0.  0.  0.  1.  1.  0.  0.  0.  1.  1.  1.
  0.  1.  0.  1.  0.  1.  1.  0.  0.  1.  1.  1.  1.  0.  1.  1.  0.  1.
  0.  0.  0.  1.  0.  1.  0.  0.  0.  0.  1.  0.  1.  1.  0.  1.  1.  1.
  1.  0.  1.  0.]

scan - Looping in Theano


In [39]:
print(y.get_value().shape)
y.get_value()


(400,)
Out[39]:
array([ 0.,  0.,  1.,  0.,  0.,  0.,  1.,  1.,  0.,  1.,  0.,  0.,  1.,
        0.,  1.,  1.,  0.,  0.,  1.,  1.,  1.,  0.,  0.,  0.,  0.,  1.,
        0.,  0.,  0.,  1.,  1.,  0.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,
        1.,  1.,  1.,  0.,  0.,  0.,  1.,  1.,  0.,  1.,  1.,  1.,  1.,
        1.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,  1.,
        0.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  1.,  1.,  1.,  0.,  0.,
        1.,  1.,  1.,  0.,  1.,  1.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,
        1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  1.,  0.,  0.,  1.,  1.,  1.,  1.,  0.,  1.,  0.,  1.,  0.,
        0.,  1.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  1.,  0.,  0.,  1.,
        1.,  0.,  0.,  1.,  1.,  1.,  0.,  1.,  1.,  1.,  0.,  0.,  0.,
        0.,  0.,  1.,  1.,  0.,  1.,  0.,  1.,  1.,  0.,  1.,  0.,  1.,
        1.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  1.,  0.,  0.,  1.,
        0.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,  1.,  0.,  1.,  1.,  0.,
        1.,  1.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  1.,  1.,
        1.,  0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,  1.,  0.,
        1.,  0.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  1.,  0.,  1.,
        1.,  1.,  1.,  0.,  0.,  1.,  1.,  0.,  0.,  0.,  1.,  1.,  1.,
        1.,  1.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  1.,
        0.,  1.,  0.,  1.,  1.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,
        1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,
        1.,  1.,  1.,  1.,  0.,  0.,  1.,  0.,  1.,  1.,  0.,  0.,  0.,
        0.,  1.,  0.,  1.,  0.,  0.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  1.,  0.,  1.,  1.,  1.,  1.,  0.,  1.,  1.,  0.,
        0.,  1.,  0.,  0.,  1.,  1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.,
        0.,  0.,  1.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,
        1.,  0.,  0.,  0.,  1.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,
        0.,  1.,  1.,  0.,  0.,  0.,  1.,  1.,  1.,  0.,  1.,  0.,  1.,
        0.,  1.,  1.,  0.,  0.,  1.,  1.,  1.,  1.,  0.,  1.,  1.,  0.,
        1.,  0.,  0.,  0.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,  1.,  0.,
        1.,  1.,  0.,  1.,  1.,  1.,  1.,  0.,  1.,  0.], dtype=float32)

Simple loop with accumulation: computing $A^k$


In [9]:
import numpy as np

In [60]:
k = theano.shared(np.int32(1),"k")
A = theano.shared(np.array(range(10)),"A")

In [61]:
result, updates = theano.scan(fn=lambda prior_result, A: prior_result * A, 
                                 outputs_info=T.ones_like(A),
                                 non_sequences=A,
                                 n_steps=k)

In [62]:
result, updates = theano.scan(fn=lambda prior_result, A: prior_result * A, 
                                 outputs_info=T.ones_like(A),
                                 non_sequences=A,
                                 n_steps=2)

In [63]:
final_result = result[-1]

In [64]:
power=theano.function(inputs=[A,k], outputs=final_result, updates=updates)


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-64-3d27e2cadcf2> in <module>()
----> 1 power=theano.function(inputs=[A,k], outputs=final_result, updates=updates)

/home/topolo/PropD/Theano/theano/compile/function.pyc in function(inputs, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input)
    324                    on_unused_input=on_unused_input,
    325                    profile=profile,
--> 326                    output_keys=output_keys)
    327     # We need to add the flag check_aliased inputs if we have any mutable or
    328     # borrowed used defined inputs

/home/topolo/PropD/Theano/theano/compile/pfunc.pyc in pfunc(params, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input, output_keys)
    447                                          rebuild_strict=rebuild_strict,
    448                                          copy_inputs_over=True,
--> 449                                          no_default_updates=no_default_updates)
    450     # extracting the arguments
    451     input_variables, cloned_extended_outputs, other_stuff = output_vars

/home/topolo/PropD/Theano/theano/compile/pfunc.pyc in rebuild_collect_shared(outputs, inputs, replace, updates, rebuild_strict, copy_inputs_over, no_default_updates)
    174             raise TypeError(('Cannot use a shared variable (%s) as explicit '
    175                              'input. Consider substituting a non-shared'
--> 176                              ' variable via the `givens` parameter') % v)
    177 
    178     # Fill update_d and update_expr with provided updates

TypeError: Cannot use a shared variable (A) as explicit input. Consider substituting a non-shared variable via the `givens` parameter

In [65]:
power=theano.function(inputs=[],outputs=final_result, updates=updates)

In [66]:
power()


Out[66]:
array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])

In [67]:
power()


Out[67]:
array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])

In [68]:
A.get_value()


Out[68]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

In [69]:
result, updates = theano.scan(fn=lambda prior_result, A: prior_result * A, 
                                 outputs_info=T.ones_like(A),
                                 non_sequences=A,
                                 n_steps=k)

In [70]:
power=theano.function(inputs=[],outputs=final_result, updates=updates)

In [71]:
power()


Out[71]:
array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])

In [72]:
k.get_value()


Out[72]:
array(1, dtype=int32)

In [76]:
k.set_value( np.int32(3) )

In [77]:
power()


Out[77]:
array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])

From this link, Loop, which for some reason a Google search overlooks most of the time, I then found this cf. good ipython notebook with explanation and more examples, a scan tutorial written by Pierre Luc Carrier

Example 1 : As simple as it gets


In [19]:
import numpy as np

In [20]:
from theano import sandbox

In [21]:
X1 = T.vector('vector1')
X2 = T.vector('vector2')

The parameter fn receives a function or lambda expression that expresses computation to do at every iteration.

Since we wish to iterate over both X1 and X2 simultaneously, provide them as sequences. This means that every iteration will operate on 2 inputs; an element from X1 and the corresponding element from X2.


In [22]:
output, updates = theano.scan(fn=lambda a, b : a * b, sequences=[X1,X2])

output contains outputs of fn from every timestep concatenated into a tensor. In our case, the output of a single timestep is a scalar so output is a vector where output[i] is the output of the ith iteration.

updates details if and how the execution of scan updates any shared variable in the graph. It should be provided as an argument when compiling the Theano function.


In [23]:
f = theano.function(inputs=[X1,X2], outputs=output, updates=updates)

If updates is omitted, the state of any shared variables modified by Scan will not be updated properly. Random number sampling, for instance, relies on shared variables. If updates is not provided, the state of the random number generator won't be updated properly and the same numbers might be sampled repeatedly. Always provide updates when compiling Theano functions.


In [10]:
X1_value = np.arange(0,5).astype(theano.config.floatX)  # [0,1,2,3,4]
X2_value = np.arange(1,6).astype(theano.config.floatX)  # [1,2,3,4,5]

In [11]:
f(X1_value,X2_value)


Out[11]:
array([  0.,   2.,   6.,  12.,  20.], dtype=float32)

In [12]:
f.maker.fgraph.toposort()


Out[12]:
[Shape_i{0}(vector1),
 GpuFromHost(vector1),
 Shape_i{0}(vector2),
 GpuFromHost(vector2),
 Elemwise{minimum,no_inplace}(Shape_i{0}.0, Shape_i{0}.0),
 Elemwise{lt,no_inplace}(Elemwise{minimum,no_inplace}.0, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}(Elemwise{lt,no_inplace}.0, Elemwise{minimum,no_inplace}.0, Shape_i{0}.0, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}[(0, 1)](Elemwise{lt,no_inplace}.0, Elemwise{minimum,no_inplace}.0, Shape_i{0}.0, TensorConstant{0}),
 Elemwise{le,no_inplace}(Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, TensorConstant{0}),
 Elemwise{le,no_inplace}(Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}[(0, 1)].0, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)](Elemwise{le,no_inplace}.0, TensorConstant{0}, Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, Shape_i{0}.0),
 Elemwise{switch,no_inplace}(Elemwise{le,no_inplace}.0, TensorConstant{0}, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)](Elemwise{le,no_inplace}.0, TensorConstant{0}, Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}[(0, 1)].0, Shape_i{0}.0),
 Elemwise{switch,no_inplace}(Elemwise{le,no_inplace}.0, TensorConstant{0}, TensorConstant{0}),
 ScalarFromTensor(Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)].0),
 ScalarFromTensor(Elemwise{switch,no_inplace}.0),
 ScalarFromTensor(Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)].0),
 ScalarFromTensor(Elemwise{switch,no_inplace}.0),
 GpuSubtensor{int64:int64:int8}(GpuFromHost.0, ScalarFromTensor.0, ScalarFromTensor.0, Constant{1}),
 GpuSubtensor{int64:int64:int8}(GpuFromHost.0, ScalarFromTensor.0, ScalarFromTensor.0, Constant{1}),
 GpuElemwise{Mul}[(0, 0)](GpuSubtensor{int64:int64:int8}.0, GpuSubtensor{int64:int64:int8}.0),
 HostFromGpu(GpuElemwise{Mul}[(0, 0)].0)]

In [16]:
output, updates = theano.scan(fn=lambda a, b : sandbox.cuda.basic_ops.gpu_from_host( a * b ), sequences=[X1,X2])
f = theano.function(inputs=[X1,X2], outputs=output, updates=updates)

In [17]:
f(X1_value,X2_value)


Out[17]:
array([  0.,   2.,   6.,  12.,  20.], dtype=float32)

In [18]:
f.maker.fgraph.toposort()


Out[18]:
[Shape_i{0}(vector2),
 GpuFromHost(vector2),
 Shape_i{0}(vector1),
 GpuFromHost(vector1),
 Elemwise{minimum,no_inplace}(Shape_i{0}.0, Shape_i{0}.0),
 Elemwise{lt,no_inplace}(Elemwise{minimum,no_inplace}.0, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}(Elemwise{lt,no_inplace}.0, Elemwise{minimum,no_inplace}.0, Shape_i{0}.0, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}(Elemwise{lt,no_inplace}.0, Elemwise{minimum,no_inplace}.0, Shape_i{0}.0, TensorConstant{0}),
 Elemwise{le,no_inplace}(Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, TensorConstant{0}),
 Elemwise{le,no_inplace}(Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)](Elemwise{le,no_inplace}.0, TensorConstant{0}, Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, Shape_i{0}.0),
 Elemwise{switch,no_inplace}(Elemwise{le,no_inplace}.0, TensorConstant{0}, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)](Elemwise{le,no_inplace}.0, TensorConstant{0}, Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, Shape_i{0}.0),
 Elemwise{switch,no_inplace}(Elemwise{le,no_inplace}.0, TensorConstant{0}, TensorConstant{0}),
 ScalarFromTensor(Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)].0),
 ScalarFromTensor(Elemwise{switch,no_inplace}.0),
 ScalarFromTensor(Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)].0),
 ScalarFromTensor(Elemwise{switch,no_inplace}.0),
 GpuSubtensor{int64:int64:int8}(GpuFromHost.0, ScalarFromTensor.0, ScalarFromTensor.0, Constant{1}),
 GpuSubtensor{int64:int64:int8}(GpuFromHost.0, ScalarFromTensor.0, ScalarFromTensor.0, Constant{1}),
 GpuElemwise{Mul}[(0, 0)](GpuSubtensor{int64:int64:int8}.0, GpuSubtensor{int64:int64:int8}.0),
 for{gpu,scan_fn}(Elemwise{minimum,no_inplace}.0, GpuElemwise{Mul}[(0, 0)].0, Elemwise{minimum,no_inplace}.0),
 HostFromGpu(for{gpu,scan_fn}.0)]

An interesting thing is that we never explicitl told Scan how many iterations to run. It was automatically inferred. When given sequences, Scan will run as many iterations as length of the shortest sequence.


In [24]:
f(X1_value, X2_value[:4])


Out[24]:
array([  0.,   2.,   6.,  12.], dtype=float32)

In [26]:
def vec_addition(a,b):
    return sandbox.cuda.basic_ops.gpu_from_host( a + b )

In [27]:
output, updates = theano.scan(fn=vec_addition, sequences=[X1,X2])
f = theano.function(inputs=[X1,X2], outputs=output, updates=updates)

In [28]:
f(X1_value,X2_value)


Out[28]:
array([ 1.,  3.,  5.,  7.,  9.], dtype=float32)

In [29]:
f.maker.fgraph.toposort()


Out[29]:
[Shape_i{0}(vector2),
 GpuFromHost(vector2),
 Shape_i{0}(vector1),
 GpuFromHost(vector1),
 Elemwise{minimum,no_inplace}(Shape_i{0}.0, Shape_i{0}.0),
 Elemwise{lt,no_inplace}(Elemwise{minimum,no_inplace}.0, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}(Elemwise{lt,no_inplace}.0, Elemwise{minimum,no_inplace}.0, Shape_i{0}.0, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}(Elemwise{lt,no_inplace}.0, Elemwise{minimum,no_inplace}.0, Shape_i{0}.0, TensorConstant{0}),
 Elemwise{le,no_inplace}(Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, TensorConstant{0}),
 Elemwise{le,no_inplace}(Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)](Elemwise{le,no_inplace}.0, TensorConstant{0}, Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, Shape_i{0}.0),
 Elemwise{switch,no_inplace}(Elemwise{le,no_inplace}.0, TensorConstant{0}, TensorConstant{0}),
 Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)](Elemwise{le,no_inplace}.0, TensorConstant{0}, Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}.0, Shape_i{0}.0),
 Elemwise{switch,no_inplace}(Elemwise{le,no_inplace}.0, TensorConstant{0}, TensorConstant{0}),
 ScalarFromTensor(Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)].0),
 ScalarFromTensor(Elemwise{switch,no_inplace}.0),
 ScalarFromTensor(Elemwise{Composite{Switch(i0, i1, minimum(i2, i3))}}[(0, 2)].0),
 ScalarFromTensor(Elemwise{switch,no_inplace}.0),
 GpuSubtensor{int64:int64:int8}(GpuFromHost.0, ScalarFromTensor.0, ScalarFromTensor.0, Constant{1}),
 GpuSubtensor{int64:int64:int8}(GpuFromHost.0, ScalarFromTensor.0, ScalarFromTensor.0, Constant{1}),
 GpuElemwise{Add}[(0, 0)](GpuSubtensor{int64:int64:int8}.0, GpuSubtensor{int64:int64:int8}.0),
 for{gpu,scan_fn}(Elemwise{minimum,no_inplace}.0, GpuElemwise{Add}[(0, 0)].0, Elemwise{minimum,no_inplace}.0),
 HostFromGpu(for{gpu,scan_fn}.0)]

In [31]:
X1.get_value


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-31-fba24b6c4b82> in <module>()
----> 1 X1.get_value

AttributeError: 'TensorVariable' object has no attribute 'get_value'

So we can do the following with scan:
$\forall \, X_1,X_2 \in \mathbb{R}^d$,
$$ + : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d \\ \verb|output|[i] = X_1[i] + X_2[i], \qquad \, \forall \, i = 0,1, \dots d-1 $$

$$ \odot : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d \\ \verb|output|[i] = X_1[i] * X_2[i], \qquad \, \forall \, i = 0,1, \dots d-1 $$

In [33]:
alpha_expn = T.scalar('alpha')  # \alpha
Xs = T.vector('Xs')  # Xs \equiv X[i]

In [38]:
lambda x, k : x**k


Out[38]:
<function __main__.<lambda>>

In [43]:
output, updates = theano.scan(fn=lambda x, k : x**k ,sequences=[Xs])


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-43-e543b36f2512> in <module>()
----> 1 output, updates = theano.scan(fn=lambda x, k : x**k ,sequences=[Xs])

/home/topolo/PropD/Theano/theano/scan_module/scan.pyc in scan(fn, sequences, outputs_info, non_sequences, n_steps, truncate_gradient, go_backwards, mode, name, profile, allow_gc, strict, return_list)
    771     # and outputs that needs to be separated
    772 
--> 773     condition, outputs, updates = scan_utils.get_updates_and_outputs(fn(*args))
    774     if condition is not None:
    775         as_while = True

TypeError: <lambda>() takes exactly 2 arguments (1 given)

In [45]:
output, updates = theano.scan(fn=lambda x, k : x**k ,sequences=[Xs,alpha_expn])


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-45-0a6d243197ac> in <module>()
----> 1 output, updates = theano.scan(fn=lambda x, k : x**k ,sequences=[Xs,alpha_expn])

/home/topolo/PropD/Theano/theano/scan_module/scan.pyc in scan(fn, sequences, outputs_info, non_sequences, n_steps, truncate_gradient, go_backwards, mode, name, profile, allow_gc, strict, return_list)
    501                 # If not we need to use copies, that will be replaced at
    502                 # each frame by the corresponding slice
--> 503                 actual_slice = seq['input'][k - mintap]
    504                 _seq_val = tensor.as_tensor_variable(seq['input'])
    505                 _seq_val_slice = _seq_val[k - mintap]

/home/topolo/PropD/Theano/theano/tensor/var.pyc in __getitem__(self, args)
    577                     self, *theano.tensor.subtensor.Subtensor.collapse(
    578                         args,
--> 579                         lambda entry: isinstance(entry, Variable)))
    580 
    581     def take(self, indices, axis=None, mode='raise'):

/home/topolo/PropD/Theano/theano/gof/op.pyc in __call__(self, *inputs, **kwargs)
    602         """
    603         return_list = kwargs.pop('return_list', False)
--> 604         node = self.make_node(*inputs, **kwargs)
    605 
    606         if config.compute_test_value != 'off':

/home/topolo/PropD/Theano/theano/tensor/subtensor.pyc in make_node(self, x, *inputs)
    479                 len(idx_list), x.type.ndim))
    480             exception.subtensor_invalid = True
--> 481             raise exception
    482 
    483         input_types = Subtensor.collapse(idx_list,

ValueError: The index list is longer (size 1) than the number of dimensions of the tensor(namely 0). You are asking for a dimension of the tensor that does not exist! You might need to use dimshuffle to add extra dimension to your tensor.

What about reduce?


In [46]:
def vec_addition(a,b):
    return sandbox.cuda.basic_ops.gpu_from_host( a + b )

In [55]:
output, updates = theano.reduce(fn=vec_addition, sequences=[X1,X2],outputs_info=[None,])

In [56]:
f = theano.function(inputs=[X1,X2], outputs=output, updates=updates)

In [57]:
f(X1_value, X2_value)


Out[57]:
array(9.0, dtype=float32)

In [58]:
def vec_elem_mult(a,b):
    return sandbox.cuda.basic_ops.gpu_from_host( a*b)

In [59]:
output, updates = theano.reduce(fn=vec_elem_mult, sequences=[X1,X2],outputs_info=[None,])

In [60]:
f = theano.function(inputs=[X1,X2], outputs=output, updates=updates)

In [61]:
f(X1_value, X2_value)


Out[61]:
array(20.0, dtype=float32)

What about matrices, tensors?


In [62]:
X1 = T.matrix('matrix1')
X2 = T.matrix('matrix2')

In [63]:
output, updates = theano.scan(fn=lambda a, b : sandbox.cuda.basic_ops.gpu_from_host( a * b ), sequences=[X1,X2])

In [64]:
f = theano.function(inputs=[X1,X2], outputs=output, updates=updates)

In [65]:
X1_value = np.arange(0,6).reshape(2,3).astype(theano.config.floatX); print(X1_value)
X2_value = np.arange(1,7).reshape(2,3).astype(theano.config.floatX); print(X2_value)


[[ 0.  1.  2.]
 [ 3.  4.  5.]]
[[ 1.  2.  3.]
 [ 4.  5.  6.]]

In [66]:
f(X1_value,X2_value)


Out[66]:
array([[  0.,   2.,   6.],
       [ 12.,  20.,  30.]], dtype=float32)

In [67]:
X1 = T.tensor3('tensor1')
X2 = T.tensor3('tensor2')

In [68]:
output, updates = theano.scan(fn=vec_addition, sequences=[X1,X2])
f = theano.function(inputs=[X1,X2], outputs=output, updates=updates)

In [69]:
X1_value = np.arange(0,24).reshape(2,3,4).astype(theano.config.floatX); print(X1_value)
X2_value = np.arange(1,25).reshape(2,3,4).astype(theano.config.floatX); print(X2_value)


[[[  0.   1.   2.   3.]
  [  4.   5.   6.   7.]
  [  8.   9.  10.  11.]]

 [[ 12.  13.  14.  15.]
  [ 16.  17.  18.  19.]
  [ 20.  21.  22.  23.]]]
[[[  1.   2.   3.   4.]
  [  5.   6.   7.   8.]
  [  9.  10.  11.  12.]]

 [[ 13.  14.  15.  16.]
  [ 17.  18.  19.  20.]
  [ 21.  22.  23.  24.]]]

In [70]:
f(X1_value,X2_value)


Out[70]:
array([[[  1.,   3.,   5.,   7.],
        [  9.,  11.,  13.,  15.],
        [ 17.,  19.,  21.,  23.]],

       [[ 25.,  27.,  29.,  31.],
        [ 33.,  35.,  37.,  39.],
        [ 41.,  43.,  45.,  47.]]], dtype=float32)

In [76]:
X_t = T.vector('X_t')
output,updates=theano.scan(fn=lambda x:x**3,sequences=[X_t,])

In [78]:
f=theano.function(inputs=[X_t,],outputs=output,updates=updates)

In [79]:
f.maker.fgraph.toposort()


Out[79]:
[Shape_i{0}(X_t),
 GpuFromHost(X_t),
 ScalarFromTensor(Shape_i{0}.0),
 Elemwise{Composite{Switch(LE(i0, i1), i1, i2)}}[(0, 0)](Shape_i{0}.0, TensorConstant{0}, TensorConstant{0}),
 ScalarFromTensor(Elemwise{Composite{Switch(LE(i0, i1), i1, i2)}}[(0, 0)].0),
 GpuSubtensor{int64:int64:int8}(GpuFromHost.0, ScalarFromTensor.0, ScalarFromTensor.0, Constant{1}),
 GpuElemwise{Pow}[(0, 0)](GpuSubtensor{int64:int64:int8}.0, CudaNdarrayConstant{[ 3.]}),
 HostFromGpu(GpuElemwise{Pow}[(0, 0)].0)]

In [80]:
f=theano.function(inputs=[X_t,],outputs=sandbox.cuda.basic_ops.gpu_from_host(output),updates=updates)

In [81]:
f.maker.fgraph.toposort()


Out[81]:
[Shape_i{0}(X_t),
 GpuFromHost(X_t),
 ScalarFromTensor(Shape_i{0}.0),
 Elemwise{Composite{Switch(LE(i0, i1), i1, i2)}}[(0, 0)](Shape_i{0}.0, TensorConstant{0}, TensorConstant{0}),
 ScalarFromTensor(Elemwise{Composite{Switch(LE(i0, i1), i1, i2)}}[(0, 0)].0),
 GpuSubtensor{int64:int64:int8}(GpuFromHost.0, ScalarFromTensor.0, ScalarFromTensor.0, Constant{1}),
 GpuElemwise{Pow}[(0, 0)](GpuSubtensor{int64:int64:int8}.0, CudaNdarrayConstant{[ 3.]})]

In [84]:
X1_value = np.arange(0,5).astype(theano.config.floatX)  # [0,1,2,3,4] 
print(X1_value)


[ 0.  1.  2.  3.  4.]

In [83]:
f(X1_value)


Out[83]:
CudaNdarray([  0.   1.   8.  27.  64.])

So this is the mathematical equivalent to
$ t=0,1,\dots T-1$, $t\in \mathbb{Z}^+$,
$X\in \mathbb{R}^T$
$$ \forall \, t = 0, 1, \dots T-1, \\ f:\mathbb{R} \to \mathbb{R} \\ X(t) \mapsto f(X(t)) = (X(t))^3 \qquad \, (\text{for example}) \\ $$


In [86]:
output,updates=theano.reduce(fn=lambda x:x**4,sequences=[X_t,],outputs_info=[None,])
f=theano.function(inputs=[X_t,],outputs=output,updates=updates)

In [87]:
f(X1_value)


Out[87]:
array(256.0, dtype=float32)

Example 2: Non-sequences

We need some variables to be available "as is" at every iteration of the loop. We don't want scan to iterate over them and give only part of them at every iteration.


In [71]:
X = T.matrix('X')
W = T.matrix('W')
b = T.vector('b')

For the sake of variety (and so lambda is the same as defining a Python function), define computation to be done at every iteration of the loop using step(), instead of lambda expression.

To have $W$ and $b$ be available at every iteration, use the argument non_sequences. Contrary to sequences, non-sequences are not iterated upon by Scan.

This means step() function will need to operate on 3 symbolic inputs; 1 for our sequences $X$, 1 for each non-sequences $W$ and $b$.

The inputs that correspond to the non-sequences are always last and in same order at the non-sequences provided to Scan. This means correspondence between inputs of the step() function and arguments to scan() is the following:

  • $v$ : individual element of the sequence $X$
  • $W,b$ : non-sequences $W,b$, respectively

In [72]:
def step(v,W,b):
    return T.dot(v,W) + b

output,updates=theano.scan(fn=step,
                            sequences=[X],
                              non_sequences=[W,b])
print(updates)


OrderedUpdates()

In [73]:
f=theano.function(inputs=[X,W,b],
                 outputs=sandbox.cuda.basic_ops.gpu_from_host( output),
                 updates=updates)

In [74]:
X_value = np.arange(-3,3).reshape(3,2).astype(theano.config.floatX)
W_value = np.eye(2).astype(theano.config.floatX)
b_value=np.arange(2).astype(theano.config.floatX)

In [88]:
print(X_value); print(W_value); print(b_value)


[[-3. -2.]
 [-1.  0.]
 [ 1.  2.]]
[[ 1.  0.]
 [ 0.  1.]]
[ 0.  1.]

In [75]:
f(X_value,W_value,b_value)


Out[75]:
CudaNdarray([[-3. -1.]
 [-1.  1.]
 [ 1.  3.]])

Notice how scan is on the first dimension (or, counting from 0, the 0th dimension), always. So 1 way to think about it is discretized time $t\in \mathbb{R} \xrightarrow{ \text{ discretize } } t\in\mathbb{Z}^+$

$$ X:\mathbb{Z}^+ \to \mathbb{R}^d \text{ or } \text{Mat}_{\mathbb{R}}(N_1,N_2) \text{ or } \tau^r_s(V), \text{ space of tensors of type } (r,s), \tau_s^r(V) = \lbrace (r+s)-\text{linear maps}, \underbrace{V\times \dots \times V}_{r} \times \underbrace{V^* \times \dots \times V^*}_{s} \to \mathbb{F} \rbrace $$

$$ \forall \, t = 0 , 1, \dots T-1, \\ \text{ e.g. } \theta \in \text{Mat}_{\mathbb{R}}(d,N_2) \\ \qquad b \in \mathbb{R}^{N_2} $$$$ f_{\theta,b}:\mathbb{R}^d \to \mathbb{R}^{N_2} \\ X(t) \mapsto f(X(t)) = X(t) \cdot \theta + b $$

and so nonsequences help to return the functional

$$ f:\text{Mat}_{\mathbb{R}}(d,N_2) \times \mathbb{R}^{N_2} \to \text{Hom}(\mathbb{R}^d, \mathbb{R}^{N_2} ) \\ f:(\theta,b) \mapsto f_{\theta,b} $$

Notice the right action of $\theta$ onto $X(t)$, as necessitated by how the size dimensions are defined. This should be duly noted, as scan only iterates across the first dimension. So to write the equivalent step, with the usual left action,


In [132]:
X = T.tensor3('X')
W = T.matrix('W')
b = T.vector('b')

In [133]:
#def step_left(v,W,b):
def step_left(v,W):
    return T.dot(W,v)  # + b  

#output, updates=theano.scan(fn=step_left,sequences=[X],non_sequences=[W,b])  
#f=theano.function(inputs=[X,W,b],outputs=output,updates=updates)

output, updates=theano.scan(fn=step_left,sequences=[X],non_sequences=[W])  
f=theano.function(inputs=[X,W],outputs=output,updates=updates)

In [134]:
X_value = np.arange(-3,3).reshape(3,2,1).astype(theano.config.floatX)
W_value = np.arange(1,5).reshape(2,2).astype(theano.config.floatX)  
#b_value=np.arange(2).reshape(2,1).astype(theano.config.floatX)

In [135]:
test_result_left =f(X_value,W_value) #,b_value)

In [136]:
test_result_left


Out[136]:
array([[[ -7.],
        [-17.]],

       [[ -1.],
        [ -3.]],

       [[  5.],
        [ 11.]]], dtype=float32)

In [127]:
test_result_left[0].shape


Out[127]:
(2, 1)

In [140]:
def step_left(v,W,b):
    return T.dot(W,v)   + b  

output, updates=theano.scan(fn=step_left,sequences=[X],non_sequences=[W,b])  
f=theano.function(inputs=[X,W,b],outputs=output,updates=updates)

In [143]:
b_value=np.arange(2).reshape(2).astype(theano.config.floatX)
test_result_left =f(X_value,W_value, b_value)

In [144]:
test_result_left


Out[144]:
array([[[ -7.,  -6.],
        [-17., -16.]],

       [[ -1.,   0.],
        [ -3.,  -2.]],

       [[  5.,   6.],
        [ 11.,  12.]]], dtype=float32)

In [151]:
def step_left(v,W,b):
    return ( T.dot(W,v)   + b )

In [158]:
output, updates=theano.scan(fn=step_left,sequences=[X],outputs_info=[None],non_sequences=[W,b])

In [159]:
test_result_left =f(X_value,W_value, b_value)

In [160]:
test_result_left


Out[160]:
array([[[ -7.,  -6.],
        [-17., -16.]],

       [[ -1.,   0.],
        [ -3.,  -2.]],

       [[  5.,   6.],
        [ 11.,  12.]]], dtype=float32)

In [108]:
np.dot(W_value,X_value)


Out[108]:
array([[[ -7.],
        [ -1.],
        [  5.]],

       [[-17.],
        [ -3.],
        [ 11.]]], dtype=float32)

In [109]:
W_value


Out[109]:
array([[ 1.,  2.],
       [ 3.,  4.]], dtype=float32)

In [110]:
np.dot(W_value,X_value[0])


Out[110]:
array([[ -7.],
       [-17.]], dtype=float32)

In [137]:
X_value


Out[137]:
array([[[-3.],
        [-2.]],

       [[-1.],
        [ 0.]],

       [[ 1.],
        [ 2.]]], dtype=float32)

In [139]:
b_value


Out[139]:
array([[ 0.],
       [ 1.]], dtype=float32)

In [156]:
T.zeros_like


Out[156]:
<function theano.tensor.basic.zeros_like>

Example 3: Reusing outputs from the previous iterations


In [161]:
def step(m_row, cumulative_sum):
    return m_row + cumulative_sum

The trick part is informing Scan that our step function expects as input the output of a previous iteration. A new parameter, outputs_info, achieves this. This parameter is used to tell Scan how we intend to use each of the ouputs that are computer at each iteration.

This parameter can be omitted (like we have done so far) when the step function doesn't depend on any output of a previous iteration.

outputs_info takes a sequence with 1 element for every output of the step() function:

  • For a non-recurrent output, element should be None.
  • For simple recurrent output, (iteration $t$ depends on value of iteration at $t-1$, say), the element must be a tensor. Scan will interpret it as being an initial state for a recurrent output, and give it as input to the first iteration, pretending it's the output value from a previous iteration.

The step() expects 1 additional input for each simple recurrent output: these inputs correspond to outputs from previous iteration and are always after inputs that correspond to sequences, but before those that correspond to non-sequences.


In [162]:
M=T.matrix('X')
s=T.vector('s')  # initial value for the cumulative sum

In [163]:
output, updates = theano.scan(fn=step,
                                sequences=[M],
                                 outputs_info=[s])

In [165]:
f=theano.function(inputs=[M, s], 
                     outputs=output,
                     updates=updates)

In [166]:
M_value=np.arange(9).reshape(3,3).astype(theano.config.floatX)
s_value=np.zeros((3,),dtype=theano.config.floatX)

In [167]:
f(M_value,s_value)


Out[167]:
array([[  0.,   1.,   2.],
       [  3.,   5.,   7.],
       [  9.,  12.,  15.]], dtype=float32)

In [168]:
M_value


Out[168]:
array([[ 0.,  1.,  2.],
       [ 3.,  4.,  5.],
       [ 6.,  7.,  8.]], dtype=float32)

In [170]:
print(M_value[0])
print(M_value[1])


[ 0.  1.  2.]
[ 3.  4.  5.]

For input $\begin{aligned} & X:\mathbb{Z}^+ \to R-\text{Module} \\ & t\mapsto X(t) \end{aligned} $, e.g. $R$-Module such as $\mathbb{R}^d, V, \text{Mat}_{\mathbb{R}}(N_1,N_2), \tau_s^r(V)$, and a function $f$ that acts at each iteration, i.e. $\forall \, t=0,1,\dots T$, $$ \begin{aligned} & f: R-\text{Module} \times R-\text{Module} \to R-\text{Module} \\ & (X_1,X_0) \mapsto X_1 + X_0 \end{aligned} $$, then we want to express $$ f(X(t),X(t-1)) = X(t)+X(t-1) \qquad \, \forall \, t=0,1\dots T-1, $$

In the end, we should get $$ X\in (R-\text{Module})^T \xrightarrow{ \text{ scan } } Y \in (R-\text{Module})^T $$ $Y \equiv $ output


In [171]:
# Further example
X=T.vector('X')
x_0=T.scalar('x_0') # initial value for the cumulative sum  
output, updates = theano.scan(fn=step, sequences=[X],outputs_info=[x_0])

In [173]:
f=theano.function(inputs=[X,x_0],outputs=output,updates=updates)

In [181]:
X_value =np.arange(9).astype(theano.config.floatX); print(X_value)
x_0_val = np.cast[theano.config.floatX](0.)


[ 0.  1.  2.  3.  4.  5.  6.  7.  8.]

In [182]:
f(X_value,x_0_val)


Out[182]:
array([  0.,   1.,   3.,   6.,  10.,  15.,  21.,  28.,  36.], dtype=float32)

This is classically what (parallel) scan should do.

Indeed,


In [183]:
output, updates = theano.scan(fn=lambda x_1,x_0 : x_1*x_0, sequences=[X],outputs_info=[x_0])
f=theano.function(inputs=[X,x_0],outputs=output,updates=updates)

In [185]:
X_value =np.arange(1,11).astype(theano.config.floatX); print(X_value)
x_0_val = np.cast[theano.config.floatX](1.)


[  1.   2.   3.   4.   5.   6.   7.   8.   9.  10.]

In [186]:
f(X_value,x_0_val)


Out[186]:
array([  1.00000000e+00,   2.00000000e+00,   6.00000000e+00,
         2.40000000e+01,   1.20000000e+02,   7.20000000e+02,
         5.04000000e+03,   4.03200000e+04,   3.62880000e+05,
         3.62880000e+06], dtype=float32)

In summary, the dictionary between the mathematics and the Python theano code for scan seems to be the following:

If $k=1$, $\forall \, t = 0,1,\dots T-1$, \begin{equation} \begin{gathered} \begin{aligned} & F:R-\text{Module}\times R-\text{Module} \to R-\text{Module} \\ & F(X(t),X(t-1)) \mapsto X(t) \end{aligned} \Longleftrightarrow \text{Python function (object) or Python lambda expression } \mapsto \verb|scan(fn= )| \\ (X(0),X(1),\dots X(T-1)) \in (R-\text{Module})^T \Longleftrightarrow \verb|scan(sequences= )| \\ X(-1) \in R-\text{Module} \Longleftrightarrow \verb|scan(outputs_info=[ ] )| \end{gathered} \end{equation}

Example 4: Reusing outputs from multiple past iterations

"...Since every example so far had only 1 output at every iteration of the loop, we will also compute, at each time step, the ratio between the new term of the Fibonacci sequence and the previous term." I think what Pierre means is that we're going to try to implement for the output, ratio, as a non-recurrent term, just for the sake of this pedagogical example.


In [187]:
def step(f_minus2,f_minus1):
    new_f = f_minus2+f_minus1
    ratio=new_f/f_minus1
    return new_f,ratio

Defining the value of outputs_info:

Recall that, for non-recurrent outputs, the value is None, and, for simple recurrent outputs, the value is a single initial state. For general recurrent outputs, where iteration $t$ may depend on multiple past values, the value is a dictionary.

That dictionary has 2 values:

  • taps : list declaring which previous values of that output every iteration will need, e.g. [-2,-1]
  • initial : tensor of initial values. If every initial value has $n$ dimensions, initial will be a single tensor of $n+1$ dimensions with as many initial values as the oldest requested tap. In the case of Fibonacci sequence, individual initial values are scalars, so initial will be a vector.

In [189]:
f_init = T.vector()
outputs_info = [dict(initial=f_init, taps=[-2,-1]),None]

In [190]:
output, updates=theano.scan(fn=step,outputs_info=outputs_info,n_steps=10)
next_fibonacci_terms=output[0]
ratios_between_terms=output[1]

In [193]:
f=theano.function(inputs=[f_init],
                    outputs=[next_fibonacci_terms,ratios_between_terms],updates=updates)

In [196]:
out=f([1,1]) 
print(len(out))
print(out[0])
print(out[1])


2
[   2.    3.    5.    8.   13.   21.   34.   55.   89.  144.]
[ 2.          1.5         1.66666663  1.60000002  1.625       1.61538458
  1.61904764  1.61764705  1.61818182  1.6179775 ]

EY : 20170324 note. Notice the n_steps parameter that's utilized now in scan. I'll try to explain it.

$\forall \, t = 0, 1, \dots T-1$, $T\Longleftrightarrow $ n_steps$=T$,

Consider \begin{equation} \begin{aligned} & F:(R-\text{Module})^k \to R-\text{Module} \\ & F(X(t-k),X(t-(k-1)),\dots X(t-1)) = X(t) \end{aligned} \Longleftrightarrow \verb|fn| \in \text{Python function (object) } \end{equation}

If $k=1$, we'll need to be given $X(0) \in R-\text{Module}$. Perhaps consider $\forall \, t=-k,-(k-1), \dots -1,0,1\dots T-1$ (``in full'').

For $k>1$, we'll need to be given (or declare) $\lbrace X(-k),X(-(k-1)),\dots X(-1)\rbrace$.

So for $k=1$, $X(-1) \in R-\text{Module}$ needed $\Longleftrightarrow $ e.g. T.scalar() if $R$-Module $=\mathbb{R}$.
for $k>1$, $(X(-k),X(-(k-1)),\dots X(-1)) \in (R-\text{Module})^k \Longleftrightarrow $ e.g. T.vector(), into ''initial'' of a Python dict, if $R$-Module $=\mathbb{R}$.

Also, for $k>1$,

$$ (-k,-(k-1), \dots -1) \Longleftrightarrow \verb|taps| = [-k,-(k-1),\dots -1] (\text{a Python}\verb| list|) $$

scan, essentially, does this: \begin{equation} \begin{gathered} (X(-k),X(-(k-1)),\dots X(-1)) \mapsto (X(0),X(1)\dots X(T-1)) \\ F(X(t-k), X(t-(k-1))\dots X(t-1)) = X(t), \qquad \, \forall \, t=0,1,\dots T-1 \end{gathered} \end{equation} given $F:(R-\text{Module})^k \to R-\text{Module}$, with $T=$ n_steps.


In [197]:
def F(Xtm2,Xtm1):
    """
    cf. https://www.math.cmu.edu/~af1p/Teaching/Combinatorics/Slides/Generating-Functions.pdf
    How many ways to spend n dollars, where everyday, you can only spend it on 1 dollar for a bun 
    OR 2 dollars for an ice cream 
    OR 2 dollars for a pastry?
    """
    new_f = Xtm2 * 2 + Xtm1 
    return new_f  

X_init = T.vector()
outputs_info=[dict(initial=X_init, taps=[-2,-1])]

In [198]:
output,updates = theano.scan(fn=F, outputs_info=outputs_info,n_steps=20)

In [199]:
f=theano.function(inputs=[X_init],outputs=output,updates=updates)

In [200]:
f([1,1])


Out[200]:
array([  3.00000000e+00,   5.00000000e+00,   1.10000000e+01,
         2.10000000e+01,   4.30000000e+01,   8.50000000e+01,
         1.71000000e+02,   3.41000000e+02,   6.83000000e+02,
         1.36500000e+03,   2.73100000e+03,   5.46100000e+03,
         1.09230000e+04,   2.18450000e+04,   4.36910000e+04,
         8.73810000e+04,   1.74763000e+05,   3.49525000e+05,
         6.99051000e+05,   1.39810100e+06], dtype=float32)

Exercises

Exercise 1 - Computing a polynomial


In [201]:
coefficients = T.vector("coefficients")
x=T.scalar("x")

In [202]:
def step(coeff,power,free_var):
    return coeff * free_var ** power

In [203]:
# Generate the components of the polynomial
max_coefficients_supported = 10000
full_range=T.arange(max_coefficients_supported)

In [204]:
components, updates = theano.scan(fn=step, 
                                 outputs_info=None,
                                sequences=[coefficients, full_range],
                                 non_sequences=x)

In [205]:
polynomial = components.sum()

In [206]:
calculate_polynomial = theano.function(inputs=[coefficients,x],outputs=polynomial,updates=updates)

In [207]:
test_coeff=np.asarray([1,0,2],dtype=theano.config.floatX)

In [208]:
calculate_polynomial(test_coeff,3)


Out[208]:
array(19.0)

Solution

cf. scan_ex1_solution.py


In [209]:
coefficients = T.vector("coefficients")
x=T.scalar("x")
max_coefficients_supported = 10000

In [214]:
def step(coeff,power,prior_value,free_var):
    return prior_value + (coeff * (free_var ** power))

In [216]:
# Generate the components of the polynomial
full_range = T.arange(max_coefficients_supported,dtype=theano.config.floatX)
outputs_info = np.zeros((),dtype=theano.config.floatX)

In [217]:
print(outputs_info)
components, updates = theano.scan(fn=step,
                                 sequences=[coefficients,full_range],
                                 outputs_info=outputs_info,
                                 non_sequences=x)


0.0

In [218]:
polynomial=components[-1]

In [219]:
calculate_polynomial=theano.function(inputs=[coefficients,x], outputs=polynomial, updates=updates)

In [220]:
test_coeff=np.asarray([1,0,2],dtype=theano.config.floatX)

In [222]:
calculate_polynomial(test_coeff,3)


Out[222]:
array(19.0, dtype=float32)

Exercise 2 - Sampling without replacement

  • takes as input a vector of probabilities and a scalar

In [223]:
probabilities = T.vector()
nb_samples = T.iscalar()  

rng = T.shared_randomstreams.RandomStreams(1234)

In [231]:
def sample_from_pvect(pvect):
    """ Provided utility function: given a symbolic vector of 
    probabilities (which MUST sum to 1), sample one element 
    and return its index
    """
    onehot_sample = rng.multinomial(n=1,pvals=pvect)
    sample = onehot_sample.argmax()
    return sample  # sample \in \mathbb{Z}^+, i.e. sample = 0,1,...K-1, with K= total number of possible outcomes

def set_p_to_zero(pvect, i):
    """ Provided utility function: given a symbolic vector of 
    probabilities and an index 'i', set the probability of the 
    i-th element to 0 and renormalize the probabilities so they
    sum to 1.
    """  
    new_pvect = T.set_subtensor(pvect[i], 0.)
    new_pvect = new_pvect / new_pvect.sum()  
    return new_pvect

EY:20170325 notes: raw_random - Low-level random numbers sample $n$ times from a multinomial distribution, defined by probabilities pvals. It's this formula :

$$ \binom{N}{k} p^k(1-p)^{N-k} $$

for a binomial distribution, probability of picking out the outcome corresponding to probability $p$, $k$ times out of $N$.

My solution, without looking at the author's beforehand

In [242]:
def sample_step(pvect):
    """ sample_step - sample without replacement, at 1 given iteration 
    """
    sample = sample_from_pvect(pvect) # \in \mathbb{Z}^+, i.e. sample=0,1,...K-1, with K=total number of possible outcomes
    new_pvect = set_p_to_zero(pvect,sample)  # we had to remove the drawn sample out of the equation
    return new_pvect

In [243]:
# this line is not needed, since we want the inputted "hard", numerical values to be the initial values
#pvect_0 = T.vector() # the initial set of probabilities for all $K$ possibilities (outcomes)

In [244]:
output,update=theano.scan(fn=sample_step,outputs_info=[probabilities],n_steps=nb_samples)

In [245]:
# Compiling the function
f = theano.function(inputs=[probabilities, nb_samples], outputs=output,updates=update)

In [246]:
# Testing the function
test_probs = np.asarray([0.6,0.3,0.1], dtype=theano.config.floatX)

In [247]:
for i in range(10):
    print(f(test_probs, 2))


[[ 0.    0.75  0.25]
 [ 0.    0.    1.  ]]
[[ 0.85714281  0.          0.14285713]
 [ 0.          0.          1.        ]]
[[ 0.    0.75  0.25]
 [ 0.    1.    0.  ]]
[[ 0.    0.75  0.25]
 [ 0.    0.    1.  ]]
[[ 0.66666669  0.33333334  0.        ]
 [ 0.          1.          0.        ]]
[[ 0.    0.75  0.25]
 [ 0.    1.    0.  ]]
[[ 0.    0.75  0.25]
 [ 0.    1.    0.  ]]
[[ 0.    0.75  0.25]
 [ 0.    1.    0.  ]]
[[ 0.85714281  0.          0.14285713]
 [ 0.          0.          1.        ]]
[[ 0.    0.75  0.25]
 [ 0.    0.    1.  ]]

Solution from author (Pierre Le Duc?)

cf. scan_ex2_solution.py


In [236]:
def step(p):
    sample=sample_from_pvect(p)
    new_p=set_p_to_zero(p,sample)
    return new_p,sample

In [237]:
output,updates = theano.scan(fn=step,outputs_info=[probabilities,None],n_steps=nb_samples)

In [238]:
modified_probabilities,samples=output

In [239]:
f=theano.function(inputs=[probabilities,nb_samples],outputs=[samples],updates=updates)

In [240]:
# Testing the function
test_probs=np.asarray([0.6,0.3,0.1],dtype=theano.config.floatX)
for i in range(10):
    print(f(test_probs,2))


[array([0, 2])]
[array([1, 0])]
[array([0, 2])]
[array([0, 2])]
[array([1, 0])]
[array([2, 0])]
[array([0, 1])]
[array([0, 1])]
[array([2, 0])]
[array([0, 1])]

Using theano's crossentropy


In [10]:
import numpy as np
import matplotlib.pyplot as plt

In [11]:
import sympy
from sympy import Symbol
from sympy.plotting import plot

In [12]:
plot(sympy.log( Symbol('x')) )


Out[12]:
<sympy.plotting.plot.Plot at 0x7f5538232b50>

In [95]:
m_test = 2  # total number of examples, m 
y_test = np.random.randint(2, size=m_test)
y_predicted_test = np.ones(m_test)*0.8
print(y_test)
print(y_predicted_test)


[1 0]
[ 0.8  0.8]

In [96]:
-(y_test*np.log(y_predicted_test)+(1.-y_test)*np.log(1.-y_predicted_test)).mean()


Out[96]:
0.91629073187415511

In [97]:
y_test = theano.shared(y_test.astype(theano.config.floatX))
y_predicted_test = theano.shared(y_predicted_test.astype(theano.config.floatX))

In [98]:
J_binary = T.nnet.binary_crossentropy(y_predicted_test, y_test).mean()
J_categorical = T.nnet.categorical_crossentropy(y_predicted_test, y_test).mean()

In [99]:
print( theano.function( inputs=[], outputs=J_categorical)() )


0.223143551314

In [100]:
print( theano.function( inputs=[], outputs=J_binary)() )


0.916290731874

making an example for categorical crossentropy


In [67]:
y3=np.zeros((m_test,3))

In [68]:
y3_predicted_cls = np.random.randint(3,size=m_test)
print(y3_predicted_cls)
for i in range(m_test):
    y3[i][y3_predicted_cls[i]] = 1.
print(y3)


[1 2]
[[ 0.  1.  0.]
 [ 0.  0.  1.]]

In [80]:
Kclses = 3
y3_predicted = np.array([np.random.dirichlet(np.ones(Kclses),size=1).flatten() for i in range(m_test)])
print(y3_predicted)


[[ 0.2914426   0.48511168  0.22344572]
 [ 0.69105175  0.22599138  0.08295687]]

In [76]:
y3_predicted[:,0].sum()


Out[76]:
0.73285247918113439

In [77]:
y3[


Out[77]:
array([[ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

In [85]:
-(y3[0][0]*np.log(y3_predicted[0][0])+(1-y3[0][0])*np.log(y3_predicted[0][np.arange(Kclses)!=0].sum()))


Out[85]:
0.34452420410160056

In [83]:
y3_predicted[0][np.arange(Kclses)!=0].sum()


Out[83]:
0.70855740198393524

In [88]:
categorical_results=[]
for i in range(m_test):
    example_row=[]
    for cls in range(Kclses):
        entropy= \
            (y3[i][cls]*np.log(y3_predicted[i][cls])+(1-y3[i][cls])*np.log(y3_predicted[i][np.arange(Kclses)!=cls].sum()))
        example_row.append(-entropy)
    categorical_results.append(example_row)
categorical_results = np.array( categorical_results )

In [89]:
categorical_results


Out[89]:
array([[ 0.3445242 ,  0.72337615,  0.25288874],
       [ 1.17458149,  0.25617226,  2.48943439]])

In [90]:
categorical_results.sum(1)


Out[90]:
array([ 1.32078909,  3.92018814])

In [91]:
categorical_results.sum(1).mean()


Out[91]:
2.6204886194141581

In [92]:
y3_test=theano.shared(y3.astype(theano.config.floatX))
y3_predicted_test=theano.shared(y3_predicted.astype(theano.config.floatX))

In [93]:
J_categorical=T.nnet.categorical_crossentropy(y3_predicted_test,y3_test).mean()

In [94]:
print( theano.function(inputs=[],outputs=J_categorical)() )


1.60640527128

In [ ]: