In [1]:
%matplotlib inline
from matplotlib import pyplot as plt
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
import pandas as pd
import tensorflow as tf
import tfgraphviz as tfg

from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
import numpy as np
import matplotlib.pyplot as plt
In [2]:
class ReluNet():
    def __init__(self, hidden_num = 1, batch_size = 100, training_epochs = 1000, learning_rate= 0.01):
        # Parameter initiation
        self.hidden_num = hidden_num
        self.batch_size = batch_size
        self.training_epochs = training_epochs
        self.learning_rate= learning_rate
        # clear previous graph
        tf.reset_default_graph()
        # Network Defination
        self.X = tf.placeholder("float")
        self.Y = tf.placeholder("float")
        
        self.hidden_num = hidden_num
        self.create_relu_net()

    def create_relu_net(self):
        # hidden layer        
        self.W = []; self.b = []; 
        self.h = []; self.h.append(self.X)       
        subnet_arch = list(np.hstack([np.array(1), np.array(self.hidden_num)]))
        for i in range(len(subnet_arch)-1):
            if i==0:
                self.W.append(tf.get_variable("Relu_W" + str(i+1), initializer= \
                            np.ones((subnet_arch[i], subnet_arch[i+1]), dtype = np.float32), trainable=False))
                self.b.append(tf.Variable(-np.linspace(0, 1, self.hidden_num[0], dtype = np.float32)))
                self.h.append(tf.nn.relu(tf.matmul(self.h[i], self.W[i]) + self.b[i]))
            else:
                self.W.append(tf.get_variable("Relu_W" + str(i+1), initializer= \
                            tf.random_normal([subnet_arch[i], subnet_arch[i+1]], mean=0.0, stddev=0.1)))
                self.b.append(tf.Variable(tf.random_normal([subnet_arch[i+1]], mean=0.0, stddev=0.001), name = "Relu_b" + str(i+1)))
                self.h.append(tf.matmul(self.h[i], self.W[i]) + self.b[i])

        self.W.append(tf.get_variable("Relu_W_Output", initializer= \
                                 tf.random_normal([subnet_arch[-1], 1], mean=0.0, stddev=0.1)))
        self.b.append(tf.Variable(tf.random_normal([1], mean=0.0, stddev=0.001), name = "Relu_b_Output"))
        self.h.append(tf.matmul(self.h[-1], self.W[-1]) + self.b[-1])

        self.y_model = self.h[-1]
        self.cost = tf.reduce_mean(tf.square(self.Y - self.y_model)) 

        # optimization
        self.train_op = tf.train.AdamOptimizer(self.learning_rate).minimize(self.cost)
        # initilize
        self.sess = tf.Session()
        self.init_op = tf.group(
                tf.local_variables_initializer(),
                tf.global_variables_initializer()
                )
        self.sess.run(self.init_op)

    def train(self, data):
        train_size = len(data)
        train_y = MinMaxScaler().fit_transform(data)
        train_x = np.linspace(1,train_size,train_size).reshape([-1,1])/train_size
        
        self.err_train = []
        training_data = np.hstack([train_x,train_y])
        for epoch in range(self.training_epochs):
            mini_batches = [training_data[k:k+self.batch_size] for k in range(0, train_size, self.batch_size)]
            for i in range(len(mini_batches)):
                self.sess.run(self.train_op, feed_dict={self.X: mini_batches[i][:,:1], self.Y: mini_batches[i][:,1:]})
            if (epoch + 1) % 100 == 0:
                err_train = self.sess.run([self.cost], feed_dict={self.X: train_x, self.Y: train_y})
                self.err_train.append(err_train)
                print("Epoch ", epoch + 1, " Loss:", err_train)
        pred_train, err_train = self.sess.run([self.y_model, self.cost], feed_dict={self.X: train_x, self.Y: train_y})
        plt.figure(1, figsize=(8, 5))
        plt.plot(self.err_train)
        return train_y, pred_train 
        
    def predict(self, x):
        for i in range(self.hidden_layer):
            print("The " + str(i+1) +" Layer:",np.round(self.sess.run(self.hs[i], feed_dict={self.X: x}),3))
        pred = self.sess.run([self.y_model], feed_dict={self.X: x})
        print("The Final Output:",pred)
        return pred

Test Example: Simple S-Curve

In [3]:
def doppler(DataNum, noise_level):
    x = np.linspace(0, 1, DataNum)
    y0 = 1/(1 + np.exp(-25*(x-0.5)))
    y = y0 + noise_level*np.random.normal(0, 1, DataNum) 
    return x, y0, y 

DataNum = 2000
x, y0, y = doppler(DataNum, 0.1)
plt.figure(1, figsize=(8, 5))
plt.plot(x, y, 'o', MarkerSize=2)
plt.plot(x, y0, color="red", LineWidth=3)
plt.show()
In [4]:
X = x.reshape([-1,1])
Y = y.reshape([-1,1])
nnpar = [20, 2]
model = ReluNet(hidden_num = nnpar, batch_size = 100, training_epochs = 1000, learning_rate= 0.001)
train_y, pred_train = model.train(Y)

Epoch  100  Loss: [0.006833778]
Epoch  200  Loss: [0.0044063623]
Epoch  300  Loss: [0.0040119058]
Epoch  400  Loss: [0.0039584157]
Epoch  500  Loss: [0.0039392863]
Epoch  600  Loss: [0.003923848]
Epoch  700  Loss: [0.003910889]
Epoch  800  Loss: [0.0039049562]
Epoch  900  Loss: [0.003901668]
Epoch  1000  Loss: [0.0038983335]
In [5]:
plt.figure(1, figsize=(8, 5))
plt.plot(train_y, '.', MarkerSize=3)
plt.plot(pred_train, color='red')
plt.title('ReLU Structure: ' + str(nnpar), fontsize=15, fontweight="bold")
plt.show()
In [6]:
W1 = model.sess.run(model.W[0]); b1 = model.sess.run(model.b[0])
W2 = model.sess.run(model.W[1]); b2 = model.sess.run(model.b[1])
W3 = model.sess.run(model.W[2]); b3 = model.sess.run(model.b[2])

fig, axs = plt.subplots(1, 2, figsize=(12, 5))
for j in range(W1.shape[1]):
    if (-b1[j] > 0) & (-b1[j] < 1):
        axs[0].plot(model.sess.run(model.h[1][:,j], feed_dict={model.X: X})) 
axs[0].set_title('ReLU Bases in Layer 1', fontweight="bold")
for k in range(W2.shape[1]):
    axs[1].plot(W3[k]*model.sess.run(model.h[2][:,k], feed_dict={model.X: X}))
axs[1].set_title('Weighted ReLU Ridges in Layer 2', fontweight="bold")
plt.show()

Test Example: Doppler Curve

In [7]:
def doppler(DataNum, noise_level):
    x = np.linspace(0, 1, DataNum)
    # y0 = np.sin((2.1*np.pi)/(x+0.25))
    y0 = np.sqrt(x*(1-x))*np.sin((2.1*np.pi)/(x+0.05))
    y = y0 + noise_level*np.random.normal(0, 1, DataNum) 
    return x, y0, y 

DataNum = 2000
x, y0, y = doppler(DataNum, 0.1)

plt.figure(1, figsize=(8, 5))
plt.plot(x, y, 'o', MarkerSize=2)
plt.plot(x, y0, color="red", LineWidth=3)
plt.show()
In [8]:
X = x.reshape([-1,1])
Y = y.reshape([-1,1])
nnpar = [40, 4]
model = ReluNet(hidden_num = nnpar, batch_size = 500, training_epochs = 5000, learning_rate= 0.001)
train_y, pred_train = model.train(Y)

Epoch  100  Loss: [0.04001194]
Epoch  200  Loss: [0.032852326]
Epoch  300  Loss: [0.028393814]
Epoch  400  Loss: [0.026127491]
Epoch  500  Loss: [0.024917228]
Epoch  600  Loss: [0.023910755]
Epoch  700  Loss: [0.022973744]
Epoch  800  Loss: [0.022091903]
Epoch  900  Loss: [0.021295909]
Epoch  1000  Loss: [0.020579666]
Epoch  1100  Loss: [0.019936418]
Epoch  1200  Loss: [0.019334577]
Epoch  1300  Loss: [0.01879752]
Epoch  1400  Loss: [0.018310769]
Epoch  1500  Loss: [0.017870773]
Epoch  1600  Loss: [0.017404065]
Epoch  1700  Loss: [0.016990392]
Epoch  1800  Loss: [0.016659135]
Epoch  1900  Loss: [0.016368749]
Epoch  2000  Loss: [0.01612094]
Epoch  2100  Loss: [0.015902413]
Epoch  2200  Loss: [0.015719015]
Epoch  2300  Loss: [0.015564656]
Epoch  2400  Loss: [0.015437052]
Epoch  2500  Loss: [0.016679242]
Epoch  2600  Loss: [0.015185165]
Epoch  2700  Loss: [0.015088987]
Epoch  2800  Loss: [0.014997373]
Epoch  2900  Loss: [0.014921343]
Epoch  3000  Loss: [0.01485786]
Epoch  3100  Loss: [0.014749292]
Epoch  3200  Loss: [0.014666612]
Epoch  3300  Loss: [0.0146089075]
Epoch  3400  Loss: [0.014556049]
Epoch  3500  Loss: [0.014508543]
Epoch  3600  Loss: [0.014549146]
Epoch  3700  Loss: [0.014400864]
Epoch  3800  Loss: [0.014363024]
Epoch  3900  Loss: [0.014327561]
Epoch  4000  Loss: [0.014296406]
Epoch  4100  Loss: [0.014252818]
Epoch  4200  Loss: [0.014226784]
Epoch  4300  Loss: [0.01420406]
Epoch  4400  Loss: [0.01418352]
Epoch  4500  Loss: [0.014149501]
Epoch  4600  Loss: [0.014130714]
Epoch  4700  Loss: [0.014116798]
Epoch  4800  Loss: [0.01422654]
Epoch  4900  Loss: [0.014070169]
Epoch  5000  Loss: [0.014056598]
In [9]:
plt.figure(1, figsize=(8, 5))
plt.plot(train_y, '.', MarkerSize=3)
plt.plot(pred_train, color='red')
plt.title('ReLU Structure: ' + str(nnpar), fontsize=15, fontweight="bold")
plt.show()
In [10]:
W1 = model.sess.run(model.W[0]); b1 = model.sess.run(model.b[0])
W2 = model.sess.run(model.W[1]); b2 = model.sess.run(model.b[1])
W3 = model.sess.run(model.W[2]); b3 = model.sess.run(model.b[2])

fig, axs = plt.subplots(1, 2, figsize=(12, 5))
for j in range(W1.shape[1]):
    if (-b1[j] > 0) & (-b1[j] < 1):
        axs[0].plot(model.sess.run(model.h[1][:,j], feed_dict={model.X: X})) 
axs[0].set_title('ReLU Bases in Layer 1', fontweight="bold")
for k in range(W2.shape[1]):
    axs[1].plot(W3[k]*model.sess.run(model.h[2][:,k], feed_dict={model.X: X}))
axs[1].set_title('Weighted ReLU Ridges in Layer 2', fontweight="bold")
plt.show()