In [1]:
import numpy as np
import theano
import theano.tensor as T
%pylab inline
import time
Populating the interactive namespace from numpy and matplotlib
In [2]:
rng = np.random
N = 4
dt= theano.shared(0.01, name="dt")
In [3]:
rn=lambda a:rng.randn(N)*0.0001+a
x = theano.shared(rn(0.1), name="x")
y = theano.shared(rn(0.1), name="y")
z = theano.shared(rn(0.), name="z")

Rossler方程式

In [4]:
#rossler                                                                                                                                                                                
a=T.dscalar("a")
b=T.dscalar("b")
c=T.dscalar("c")

rossler_param=[a,b,c]

xn=x+dt*(-y-z)
yn=y+dt*(x+a*y)
zn=z+dt*(b+z*(x-c))
In [5]:
rossler = theano.function(
                inputs=rossler_param,
                outputs=[xn],
                updates={x:xn,y:yn,z:zn},
                on_unused_input='warn'
                )
-c:5: UserWarning: The parameter 'updates' of theano.function() expects an OrderedDict, got <type 'dict'>. Using a standard dictionary here results in non-deterministic behavior. You should use an OrderedDict if you are using Python 2.7 (theano.compat.python2x.OrderedDict for older python), or use a list of (shared, update) pairs. Do not just convert your dictionary to this type before the call as the conversion will still be non-deterministic.
In [6]:
tn=20000

for文で実行した場合

In [7]:
tp=np.array([0]*N)
t0=time.clock()
for i in xrange(tn):
    tp=np.vstack((tp,rossler(0.1,0.1,14.)))
print time.clock()-t0
1.61605205877

scan(loop構文)での実行

In [8]:
x=T.dvector("X")
y=T.dvector("Y")
z=T.dvector("Z")

n=T.iscalar("n")

def rossler(x,y,z,a,b,c):
    xn=x+dt*(-y-z)
    yn=y+dt*(x+a*y)
    zn=z+dt*(b+z*(x-c))
    return [xn,yn,zn]

ros_result, ros_updates = theano.scan(fn=lambda x,y,z,aa,bb,cc:rossler(x,y,z,a,b,c),
                              sequences=None,
                              outputs_info=[x,y,z],
                              non_sequences=[a,b,c],
                              n_steps=n)

rosslers= theano.function(inputs=[x,y,z,a,b,c,n], outputs=ros_result, updates=ros_updates)
C:\Python27\lib\site-packages\theano\scan_module\scan_perform_ext.py:85: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility
  from scan_perform.scan_perform import *
In [9]:
an=rn(0.1)
t0=time.clock()
res=rosslers(an,an,an,0.1,0.1,14.,tn)
print time.clock()-t0
0.144549848153

for文と比べて圧倒的に速い

In [10]:
#CPU
an=rn(0.1)
t0=time.clock()
res=rosslers(an,an,an,0.1,0.1,14.,tn)
print time.clock()-t0
0.152729158485
In [153]:
#len(res)
an
#len(res[0])
Out[153]:
array([ 0.09995907,  0.10008567,  0.10000477,  0.09998015])

最後の値だけ残せば良い書き方

In [141]:
values, updates = theano.scan(lambda a,b,c:{x:xn,y:yn,z:zn},
                              non_sequences=[a,b,c],
                              n_steps=n)
rosslers = theano.function([a,b,c,n], values, updates=updates, allow_input_downcast=True)
In [142]:
tp=np.array([0]*N)
t0=time.clock()
res=rosslers(0.1,0.1,14.,tn)
print time.clock()-t0
0.000347134076947

結果のplot

In [260]:
cn=tp.T
for cc in cn:
    plt.plot(cc[1:])
In [264]:
for cc in cn:
    plt.plot(cc[:10000])

微妙の初期値の違いが増幅されている。カオス的振る舞い

In [17]:
rossler_param=T.dvector("param")
dt= theano.shared(0.01, name="dt")

x=T.dvector("X")
y=T.dvector("Y")
z=T.dvector("Z")

n=T.iscalar("n")

def dx(xx):
    return -xx[1]-xx[2]

def rossler(x,y,z,p):
    a=p[0]
    b=p[1]
    c=p[2]
    xx=[x,y,z]
#    xn=x+dt*(-y-z)
    xn=x+dt*dx(xx)
    yn=y+dt*(x+a*y)
    zn=z+dt*(b+z*(x-c))
    return [xn,yn,zn]

ros_result, ros_updates = theano.scan(rossler,
                              sequences=None,
                              outputs_info=[x,y,z],
                              non_sequences=[rossler_param],
                              n_steps=n)

rosslers= theano.function(inputs=[x,y,z,rossler_param,n],
                          outputs=ros_result, 
                          updates=ros_updates)
C:\Python27\lib\site-packages\theano\scan_module\scan_perform_ext.py:85: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility
  from scan_perform.scan_perform import *
In [24]:
an=rn(0.1)
t0=time.clock()
res=rosslers(an,an,an,[0.1,0.1,14.],tn)
print time.clock()-t0
0.149190887714
In [25]:
cn=res.T
for cc in cn:
    plt.plot(cc[1:])
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-25-a1334f50f276> in <module>()
----> 1 cn=res.T
      2 for cc in cn:
      3     plt.plot(cc[1:])

AttributeError: 'list' object has no attribute 'T'
In [32]:
res[0][0][0]
Out[32]:
0.097979097139283808
In [33]:
res[0][0]
Out[33]:
array([ 0.0979791 ,  0.09791782,  0.09794244,  0.09798657])
In [ ]: