window_regres2.py

"""
x[t] = fp(u)[t]
     = mean( v[t-T+1:t+1] ) = ( v[t-T+1] + v[t-T+2] + ... + v[t-1] + v[t] )/T
v[t] = tanh(u[t])

E = sum( 0.5 * (x - y)**2 )
dE/dx[t] = x[t] - y[t]
dE/dv[t] = ( dE/dx[t+T-1] + dE/dx[t+T-2] + ... + dE/dx[t+1] + dE/dx[t] )/T
dE/du[t] = dE/dv[t] * (1 - x[t]**2)

result: http://www.flickr.com/photos/arataka/tags/up20091213143236/show/
"""

import pylab
import numpy
tanh = numpy.tanh

class WindowRegres(object):

    def __init__(self, steps, wlen, _learn_rate=5, title=''):
        self.t          = numpy.arange(steps+wlen)
        self.x          = numpy.zeros((steps,), dtype=numpy.float)
        self.x          = numpy.zeros((steps,), dtype=numpy.float)
        self.v          = numpy.zeros((steps+wlen,), dtype=numpy.float)
        self.u          = numpy.zeros((steps+wlen,), dtype=numpy.float)
        self.dEdv       = numpy.zeros((steps+wlen,), dtype=numpy.float)
        self.dlist      = range(wlen+1)
        self.wlen       = wlen
        self.learn_rate = float(_learn_rate)/steps
        self.title = title

    def fp(self):
        u = self.u
        v = self.v = tanh(u)
        x = self.x
        x.fill(0)
        steps = self.x.shape[0]
        wlen = self.wlen
        for d in self.dlist:
            x += v[d:steps+d]
        x /= (wlen+1)

    def bp(self):
        y = self.y
        v = self.v
        x = self.x
        dEdv = self.dEdv
        steps = self.x.shape[0]
        wlen = self.wlen
        self.E = numpy.sum( 0.5*(x-y)**2 )/steps
        dEdx = x-y
        dEdv.fill(0)
        for d in self.dlist:
            dEdv[d:steps+d] += dEdx
        #dEdv /= (wlen+1)
        dEdu = dEdv * (1-v**2)

        self.u -= self.learn_rate * dEdu

    def learn(self, epochs, draw_step=None):
        self.rec = numpy.zeros((epochs,1), dtype=numpy.float)
        for ep in xrange(epochs):
            self.fp()
            self.bp()
            self.rec[ep,0] = self.E
            if draw_step != None and ep % draw_step == 0:
                self.plot()
        self.plot()

    def plot(self):
        wlen = self.wlen
        pylab.clf()
        
        pylab.subplot(311)
        pylab.semilogx(self.rec[:,0])
        pylab.ylim(0,None)
        pylab.title(self.title)
        pylab.ylabel('Error')
        
        pylab.subplot(312)
        pylab.plot(self.t[wlen:], self.y, linewidth=0.5, label='y')
        pylab.plot(self.t[wlen:], self.x, linewidth=2.0, label='x')
        pylab.plot(self.t       , self.v, linewidth=0.5, label='v')
        pylab.ylim(-1.05,1.05)
        pylab.legend()
        pylab.ylabel('y, x, v')
        
        pylab.subplot(313)
        pylab.plot(self.t, self.u, label='u')
        pylab.legend()
        pylab.ylabel('u')
        
        pylab.draw()
        pylab.show()


tch_len = 500
for wlen in range(0,60,10):
    wr = WindowRegres(tch_len,wlen,10, title='T=%d'%wlen)
    t = 2 * numpy.pi * numpy.arange(tch_len) / 50
    wr.y = numpy.sin(t) * 0.8 + numpy.random.randn(tch_len) * 0.5
    wr.learn(10000)
    pylab.savefig('T%03d.png'%wlen)