from pylab import *
from scipy.ndimage import filters
def F(a): return array(a,'f')
figsize(12,6)

temp = F(mean(imread("page.png"),2))
roi = (slice(500,1000),slice(500,1000))
target = F(1-(temp[roi]>0.5))
gray()
subplot(121); imshow(target)
image = maximum(filters.gaussian_filter(target+0.2*randn(*target.shape),3.0),(rand(*image.shape)>0.95))
image -= amin(image); image /= amax(image)
subplot(122); imshow(image)
h,w = image.shape

def sigmoid(x): return 1/(1+exp(-x))

N = 8
l = arange(2*r+1)
r = 3
fs = array([0.01*ones((2*r+1,2*r+1)) for i in range(8)])
fs[0,:,r] = .1
fs[1,r,:] = .1
fs[2,r-1:r+2,r-1:r+2] = .1
fs[3,:2,:2] = .1; fs[3,:2,-2:] = .1; fs[3,-2:,:2] = .1; fs[3,-2:,-2:] = .1
fs[4,l,l] = .1
fs[5,l,l[::-1]] = .1
fs[6,r-1:r+2,r] = .1
fs[7,r,r-1:r+2] = .1
for i in range(N):
    fs -= mean(fs)
    fs[i] /= sum(fs**2)**.5
thetas = zeros(N)
ws = sin(1.7*arange(N))
theta2 = 0.0
figsize(12,6)
for i in range(N):
    subplot(2,4,i+1); imshow(fs[i])

hs = array([sigmoid(filters.convolve(image,f)+thetas[i]) for i,f in enumerate(fs)])
for i in range(N):
    subplot(2,4,i+1); imshow(hs[i])
hs.shape

pred = sigmoid(dot(ws,hs.reshape(N,h*w)).reshape(h,w)+theta2)
imshow(pred)

def train(niter,eta=0.1):
    global N,r,fs,thetas,ws,theta2
    for iter in range(niter):
        hs = array([sigmoid(filters.convolve(image,f)+thetas[i]) for i,f in enumerate(fs)])
        pred = sigmoid(dot(ws,hs.reshape(N,h*w)).reshape(h,w)+theta2)
        err = sum((pred-target)**2)/(w*h)
        if iter%10==0: sys.stderr.write("\n")
        sys.stderr.write("%d:%.3f "%(iter,err))
        delta = (pred-target)*pred*(1-pred)
        delta /= h*w
        dw = array([sum(delta*hs[i]) for i,f in enumerate(fs)])
        for k,f in enumerate(fs):
            delta1 = delta*ws[k]*hs[k]*(1-hs[k])
            df = array([[sum(delta1*roll(roll(image,i,0),j,1)) for j in range(-r,r+1)] for i in range(-r,r+1)])
            fs[k] -= eta*df
            thetas[k] -= eta*sum(delta1)
        ws -= eta*dw
        theta2 -= eta*sum(delta)
    

def forward(image):
    hs = array([sigmoid(filters.convolve(image,f)+thetas[i]) for i,f in enumerate(fs)])
    pred = sigmoid(dot(ws,hs.reshape(N,h*w)).reshape(h,w)+theta2)
    return pred

fs.shape

train(500,eta=5.0)
pred = forward(image)
imshow(pred)

figsize(12,12)
subplot(131); imshow(target)
subplot(132); imshow(image)
subplot(133); imshow(pred)

hs = array([sigmoid(filters.convolve(image,f)+thetas[i]) for i,f in enumerate(fs)])
for i in range(N):
    subplot(2,4,i+1); imshow(hs[i])
hs.shape

for i in range(N):
    subplot(2,4,i+1); imshow(fs[i])

print thetas
print ws
print theta2