#!/usr/bin/env python
# From http://arctrix.com/nas/python/bpnn.py

import time, math, random
random.seed(0)

def sigmoid(x):
    return math.tanh(x)

def dsigmoid(y):
    return 1.0 - y**2

def rand(a, b):
    return (b-a)*random.random() + a

def makeMatrix(I, J, fill=0.0):
    m = []
    for i in range(I):
        m.append([fill]*J)
    return m

class NN:
    def __init__(self, ni, nh, no):
        # number of input, hidden, and output nodes
        self.ni = ni + 1 # +1 for bias node
        self.nh = nh
        self.no = no

        # activations for nodes
        self.ai = [1.0]*self.ni
        self.ah = [1.0]*self.nh
        self.ao = [1.0]*self.no

        # create weights
        self.wi = makeMatrix(self.ni, self.nh)
        self.wo = makeMatrix(self.nh, self.no)
        # set them to random vaules
        for i in range(self.ni):
            for j in range(self.nh):
                self.wi[i][j] = rand(-0.2, 0.2)
        for j in range(self.nh):
            for k in range(self.no):
                self.wo[j][k] = rand(-2.0, 2.0)

        # last change in weights for momentum
        self.ci = makeMatrix(self.ni, self.nh)
        self.co = makeMatrix(self.nh, self.no)

    def update(self, inputs):
        if len(inputs) != self.ni-1:
            raise ValueError('wrong number of inputs')

        # input activations
        for i in range(self.ni-1):
            #self.ai[i] = sigmoid(inputs[i])
            self.ai[i] = inputs[i]

        # hidden activations
        for j in range(self.nh):
            sum = 0.0
            for i in range(self.ni):
                sum = sum + self.ai[i] * self.wi[i][j]
            self.ah[j] = sigmoid(sum)

        # output activations
        for k in range(self.no):
            sum = 0.0
            for j in range(self.nh):
                sum = sum + self.ah[j] * self.wo[j][k]
            self.ao[k] = sigmoid(sum)

        return self.ao[:]


    def backPropagate(self, targets, N, M):
        if len(targets) != self.no:
            raise ValueError('wrong number of target values')

        # calculate error terms for output
        output_deltas = [0.0] * self.no
        for k in range(self.no):
            error = targets[k]-self.ao[k]
            output_deltas[k] = dsigmoid(self.ao[k]) * error

        # calculate error terms for hidden
        hidden_deltas = [0.0] * self.nh
        for j in range(self.nh):
            error = 0.0
            for k in range(self.no):
                error = error + output_deltas[k]*self.wo[j][k]
            hidden_deltas[j] = dsigmoid(self.ah[j]) * error

        # update output weights
        for j in range(self.nh):
            for k in range(self.no):
                change = output_deltas[k]*self.ah[j]
                self.wo[j][k] = self.wo[j][k] + N*change + M*self.co[j][k]
                self.co[j][k] = change
                #print N*change, M*self.co[j][k]

        # update input weights
        for i in range(self.ni):
            for j in range(self.nh):
                change = hidden_deltas[j]*self.ai[i]
                self.wi[i][j] = self.wi[i][j] + N*change + M*self.ci[i][j]
                self.ci[i][j] = change

        # calculate error
        error = 0.0
        for k in range(len(targets)):
            error = error + 0.5*(targets[k]-self.ao[k])**2
        return error


    def test(self, patterns, threshold=0.5):
        for p in patterns:
            res = self.update(p[0])
            res = list(map(lambda x:-1 if x < -threshold else 1 if x > threshold else 0,res))
            print(p[0], '->', "+" if res == p[1] else "-", "Got: ", res, "Expected: ", p[1])

    def weights(self):
        print('Input weights:')
        for i in range(self.ni):
            print(self.wi[i])
        print()
        print('Output weights:')
        for j in range(self.nh):
            print(self.wo[j])

    def train(self, patterns, iterations=1000, N=0.5, M=0.1):
        # N: learning rate
        # M: momentum factor
        for i in range(iterations):
            error = 0.0
            for p in patterns:
                inputs = p[0]
                targets = p[1]
                self.update(inputs)
                error = error + self.backPropagate(targets, N, M)
            if i % 100 == 0:
                print('error %-.5f' % error)

def fizzbuzz():
    CONSTANTS=[1] # bias node
    def int2bits(n, digits):
        res = [0]*digits
        i = 0
        while n > 0:
            res[i] = int(n) & 1
            i = i + 1
            n = n // 2
        return res

    def int2quadbits(n, quadbits):
        res = [0,0,0,0,] * quadbits
        i = 0
        while n > 0:
            res[i:i+4]=int2bits(n%15,4)
            i += 4
            n=n//15
        return CONSTANTS + res
    def fizzbuzz(j):
        F=0
        if j % 15 == 0: return [1,1,F]
        if j % 3 == 0: return [1,F,F]
        if j % 5 == 0: return [F,1,F]
        return [F,F,1]

    assert int2bits(11, 4) == [1, 1, 0, 1]
    assert int2quadbits(11, 1) == CONSTANTS+[1, 1, 0, 1]
    assert int2quadbits(37, 1) == CONSTANTS+[1, 1, 1, 0,  0,1,0,0]

    pat100 = [[int2quadbits(j,2), fizzbuzz(j)] for j in range(1,101)]
    pat200 = [[int2quadbits(j,2), fizzbuzz(j)] for j in range(101,201)]

    inputs  = len(pat200[0][0])
    outputs = len(pat200[0][1])

    bp = NN(inputs, inputs+outputs, outputs)
    bp.train(pat200, N=0.05,M=0.2)
    bp.test(pat100)


def benchmark():
    # XOR
    patterns = [
        [[-1,-1], [-1]],
        [[-1,1], [1]],
        [[1,-1], [1]],
        [[1,1], [-1]],
    ]

    bp = NN(2, 3, 1)
    bp.train(patterns, 10000)
    bp.test(patterns)

if __name__ == '__main__':
    start = time.clock()
    #benchmark()
    fizzbuzz()
    end = time.clock()
    print
    print (end - start)