Deep Belief Nets.go

DBN.go

package main

import (
	"fmt"
	"math"
	"math/rand"
)


func uniform(min float64, max float64) float64 {
	return rand.Float64() * (max - min) + min
}

func binomial(n int, p float64) int {
	if p < 0 || p > 1 { return 0 }

	c := 0
	var r float64
	
	for i := 0; i < n; i++ {
		r = rand.Float64()		
		if r < p { c++ }
	}

	return c
}

func sigmoid(x float64) float64 {
	return 1.0 / (1.0 + math.Exp(-x))
}


type DBN struct {
	N int
	n_ins int
	hidden_layer_sizes []int
	n_outs int
	n_layers int
	sigmoid_layers []HiddenLayer
	rbm_layers []RBM
	log_layer LogisticRegression
}

type HiddenLayer struct {
	N int
	n_in int
	n_out int
	W [][]float64
	b []float64
}

type RBM struct {
	N int
	n_visible int
	n_hidden int
	W [][]float64
	hbias []float64
	vbias []float64
}

type LogisticRegression struct {
	N int
	n_in int
	n_out int
	W [][]float64
	b []float64
}


// DBN
func DBN__construct(this *DBN, N int, n_ins int, hidden_layer_sizes []int, n_outs int, n_layers int) {
	var input_size int
	
	this.N = N
	this.n_ins = n_ins
	this.hidden_layer_sizes = hidden_layer_sizes
	this.n_outs = n_outs
	this.n_layers = n_layers

	this.sigmoid_layers = make([]HiddenLayer, n_layers)
	this.rbm_layers = make([]RBM, n_layers)
	
	// construct multi-layer
	for i := 0; i < n_layers; i++ {
		if i == 0 {
			input_size = n_ins
		} else {
			input_size = hidden_layer_sizes[i-1]
		}

		// construct sigmoid_layer
		HiddenLayer__construct(&(this.sigmoid_layers[i]), N, input_size, hidden_layer_sizes[i], nil, nil)

		// construct rbm_layer
		RBM__construct(&(this.rbm_layers[i]), N, input_size, hidden_layer_sizes[i], this.sigmoid_layers[i].W, this.sigmoid_layers[i].b, nil)
	}

	// layer for output using LogisticRegression
	LogisticRegression__construct(&(this.log_layer), N, hidden_layer_sizes[n_layers-1], n_outs)
}

func DBN_pretrain(this *DBN, train_X [][]int, lr float64, k int, epochs int){
	var (
		layer_input []int
		prev_layer_input_size int
		prev_layer_input []int
	)


	for i := 0; i < this.n_layers; i++ {	// layer-wise
		for epoch := 0; epoch < epochs; epoch++ {	 // training epochs
			for n := 0; n < this.N; n++ {	 // input x1...xN

				// layer input
				for l := 0; l <= i; l++ {
					if l == 0 {
						layer_input = make([]int, this.n_ins)
						for j := 0; j < this.n_ins; j++ { layer_input[j] = train_X[n][j] }
					} else {
						if l == 1 {
							prev_layer_input_size = this.n_ins
						} else {
							prev_layer_input_size = this.hidden_layer_sizes[l-2]
						}

						prev_layer_input = make([]int, prev_layer_input_size)
						for j := 0; j < prev_layer_input_size; j++ { prev_layer_input[j] = layer_input[j] }

						layer_input = make([]int, this.hidden_layer_sizes[l-1])

						HiddenLayer_sample_h_given_v(&(this.sigmoid_layers[l-1]), prev_layer_input, layer_input)
					}
				}

				RBM_contrastive_divergence(&(this.rbm_layers[i]), layer_input, lr, k)
			}
		}
	}
}

func DBN_finetune(this *DBN, train_X [][]int, train_Y [][]int, lr float64, epochs int) {
	var (
		layer_input []int
		prev_layer_input []int
	)

	for epoch := 0; epoch < epochs; epoch++ {
		for n := 0; n < this.N; n++ {	 // input x1...xN

			// layer input
			for i := 0; i < this.n_layers; i++ {
				if i == 0 {
					prev_layer_input = make([]int, this.n_ins)
					for j := 0; j < this.n_ins; j++ { prev_layer_input[j] = train_X[n][j] }
				} else {
					prev_layer_input = make([]int, this.hidden_layer_sizes[i-1])
					for j:= 0; j < this.hidden_layer_sizes[i-1]; j++ { prev_layer_input[j] = layer_input[j] }
				}

				layer_input = make([]int, this.hidden_layer_sizes[i])
				HiddenLayer_sample_h_given_v(&(this.sigmoid_layers[i]), prev_layer_input, layer_input)
			}

			LogisticRegression_train(&(this.log_layer), layer_input, train_Y[n], lr)
		}
		// lr *= 0.95
	}
}

func DBN_predict(this *DBN, x []int, y []float64) {
	var (
		layer_input []float64
	)	
	prev_layer_input := make([]float64, this.n_ins)
	for j := 0; j < this.n_ins; j++ { prev_layer_input[j] = float64(x[j]) }

	
	// layer activation
	for i := 0; i < this.n_layers; i++ {
		layer_input = make([]float64, this.sigmoid_layers[i].n_out)

		for k := 0; k < this.sigmoid_layers[i].n_out; k++ {
			linear_outuput := 0.0

			for j := 0; j < this.sigmoid_layers[i].n_in; j++ {
				linear_outuput += this.sigmoid_layers[i].W[k][j] * prev_layer_input[j]
			}
			linear_outuput += this.sigmoid_layers[i].b[k]
			layer_input[k] = sigmoid(linear_outuput)
		}

		if i < this.n_layers-1 {
			prev_layer_input = make([]float64, this.sigmoid_layers[i].n_out)

			for j := 0; j < this.sigmoid_layers[i].n_out; j++ {
				prev_layer_input[j] = layer_input[j]
			}
		}
	}

	for i := 0; i < this.log_layer.n_out; i++ {
		y[i] = 0
		for j := 0; j < this.log_layer.n_in; j++ {
			y[i] += this.log_layer.W[i][j] * layer_input[j]
		}
		y[i] += this.log_layer.b[i]
	}

	LogisticRegression_softmax(&(this.log_layer), y)
}


// HiddenLayer
func HiddenLayer__construct(this *HiddenLayer, N int, n_in int, n_out int, W [][]float64, b []float64) {
	a := 1.0 / float64(n_in)

	this.N = N
	this.n_in = n_in
	this.n_out = n_out

	if W == nil {
		this.W = make([][]float64, n_out)
		for i := 0; i < n_out; i++ { this.W[i] = make([]float64, n_in) }
		
		for i := 0; i < n_out; i++ {
			for j := 0; j < n_in; j++ {
				this.W[i][j] = uniform(-a, a)
			}
		}
	} else {
		this.W = W
	}

	if b == nil {
		this.b = make([]float64, n_out)
	} else {
		this.b = b
	}
}

func HiddenLayer_output(this *HiddenLayer, input []int, w []float64, b float64) float64 {
	linear_output := 0.0

	for j := 0; j < this.n_in; j++ {
		linear_output += w[j] * float64(input[j])
	}
	linear_output += b

	return sigmoid(linear_output)
}

func HiddenLayer_sample_h_given_v(this *HiddenLayer, input []int, sample []int) {
	for i := 0; i < this.n_out; i++ {
		sample[i] = binomial(1, HiddenLayer_output(this, input, this.W[i], this.b[i]))
	}
}


// RBM
func RBM__construct(this *RBM, N int, n_visible int, n_hidden int, W [][]float64, hbias []float64, vbias []float64) {
	a := 1.0 / float64(n_visible)

	this.N = N
	this.n_visible = n_visible
	this.n_hidden = n_hidden

	if W == nil {
		this.W = make([][]float64, n_hidden)
		for i := 0; i < n_hidden; i++ { this.W[i] = make([]float64, n_visible) }

		for i := 0; i < n_hidden; i++ {
			for j := 0; j < n_visible; j++ {
				this.W[i][j] = uniform(-a, a)
			}
		}
	} else {
		this.W = W
	}

	if hbias == nil {
		this.hbias = make([]float64, n_hidden)
	} else {
		this.hbias = hbias
	}

	if vbias == nil {
		this.vbias = make([]float64, n_visible)
	} else {
		this.vbias = vbias
	}
}

func RBM_contrastive_divergence(this *RBM, input []int, lr float64, k int) {
	ph_mean := make([]float64, this.n_hidden)
	ph_sample := make([]int, this.n_hidden)
	nv_means := make([]float64, this.n_visible)
	nv_samples := make([]int, this.n_visible)
	nh_means := make([]float64, this.n_hidden)
	nh_samples := make([]int, this.n_hidden)

	/* CD-k */
	RBM_sample_h_given_v(this, input, ph_mean, ph_sample)

	for step := 0; step < k; step++ {
		if step == 0 {
			RBM_gibbs_hvh(this, ph_sample, nv_means, nv_samples, nh_means, nh_samples)
		} else {
			RBM_gibbs_hvh(this, nh_samples, nv_means, nv_samples, nh_means, nh_samples)
		}
	}

	for i := 0; i < this.n_hidden; i++ {
		for j := 0; j < this.n_visible; j++ {
			this.W[i][j] += lr * (ph_mean[i] * float64(input[j]) - nh_means[i] * float64(nv_samples[j])) / float64(this.N)
		}
		this.hbias[i] += lr * (float64(ph_sample[i]) - nh_means[i]) / float64(this.N)
	}

	for i := 0; i < this.n_visible; i++ {
		this.vbias[i] += lr * float64(input[i] - nv_samples[i]) / float64(this.N)
	}
}

func RBM_sample_h_given_v(this *RBM, v0_sample []int, mean []float64, sample []int) {
	for i := 0; i < this.n_hidden; i++ {
		mean[i] = RBM_propup(this, v0_sample, this.W[i], this.hbias[i])
		sample[i] = binomial(1, mean[i])
	}
}

func RBM_sample_v_given_h(this *RBM, h0_sample []int, mean []float64, sample []int) {
	for i := 0; i < this.n_visible; i++ {
		mean[i] = RBM_propdown(this, h0_sample, i, this.vbias[i])
		sample[i] = binomial(1, mean[i])
	}
}

func RBM_propup(this *RBM, v []int, w []float64, b float64) float64 {
	pre_sigmoid_activation := 0.0
	
	for j := 0; j < this.n_visible; j++ {
		pre_sigmoid_activation += w[j] * float64(v[j])
	}
	pre_sigmoid_activation += b
	
	return sigmoid(pre_sigmoid_activation)
}

func RBM_propdown(this *RBM,	h []int, i int, b float64) float64 {
	pre_sigmoid_activation := 0.0
	
	for j := 0; j < this.n_hidden; j++ {
		pre_sigmoid_activation += this.W[j][i] * float64(h[j])
	}
	pre_sigmoid_activation += b

	return sigmoid(pre_sigmoid_activation)
}

func RBM_gibbs_hvh(this *RBM, h0_sample []int, nv_means []float64, nv_samples []int, nh_means []float64, nh_samples []int) {
	RBM_sample_v_given_h(this, h0_sample, nv_means, nv_samples)
	RBM_sample_h_given_v(this, nv_samples, nh_means, nh_samples)
}

func RBM_reconstruct(this *RBM, v []int, reconstructed_v []float64) {
	h := make([]float64, this.n_hidden)
	var pre_sigmoid_activation float64

	for i := 0; i < this.n_hidden; i++ {
		h[i] = RBM_propup(this, v, this.W[i], this.hbias[i])
	}

	for i := 0; i < this.n_visible; i++ {
		pre_sigmoid_activation = 0.0
		for j := 0; j < this.n_hidden; j++ {
			pre_sigmoid_activation += this.W[j][i] * h[j]
		}
		pre_sigmoid_activation += this.vbias[i]

		reconstructed_v[i] = sigmoid(pre_sigmoid_activation)
	}
}

// LogisticRegression
func LogisticRegression__construct(this *LogisticRegression, N int, n_in int, n_out int) {
	this.N = N
	this.n_in = n_in
	this.n_out = n_out

	this.W = make([][]float64, n_out)
	for i := 0; i < n_out; i++ { this.W[i] = make([]float64, n_in) }
	
	this.b = make([]float64, n_out)
}

func LogisticRegression_train(this *LogisticRegression, x []int, y []int, lr float64) {
	p_y_given_x := make([]float64, this.n_out)
	dy := make([]float64, this.n_out)
	
	for i := 0; i < this.n_out; i++ {
		p_y_given_x[i] = 0
		for j := 0; j < this.n_in; j++ {
			p_y_given_x[i] += this.W[i][j] * float64(x[j])
		}
		p_y_given_x[i] += this.b[i]
	}
	LogisticRegression_softmax(this, p_y_given_x)
	
	for i := 0; i < this.n_out; i++ {
		dy[i] = float64(y[i]) - p_y_given_x[i]
		
		for j := 0; j < this.n_in; j++ {
			this.W[i][j] += lr * dy[i] * float64(x[j]) / float64(this.N)
		}

		this.b[i] += lr * dy[i] / float64(this.N)
	}
	
}

func LogisticRegression_softmax(this *LogisticRegression, x []float64) {
	var (
		max float64
		sum float64
	)

	for i := 0; i < this.n_out; i++ { if max < x[i] {max = x[i]} }
	for i := 0; i < this.n_out; i++ {
		x[i] = math.Exp(x[i] - max)
		sum += x[i]
	}

	for i := 0; i < this.n_out; i++ { x[i] /= sum }
}

func LogisticRegression_predict(this *LogisticRegression, x []int, y []float64) {
	for i := 0; i < this.n_out; i++ {
		y[i] = 0
		for j := 0; j < this.n_in; j++ {
			y[i] += this.W[i][j] * float64(x[j])
		}
		y[i] += this.b[i]
	}

	LogisticRegression_softmax(this, y)
}



func test_dbn() {
	rand.Seed(0)

	pretrain_lr := 0.1
	pretraining_epochs := 1000
	k := 1
	fintune_lr := 0.1
	fintune_epochs := 500

	train_N := 6
	test_N := 4
	n_ins := 6
	n_outs := 2
	hidden_layer_sizes := []int {3, 3}
	n_layers := len(hidden_layer_sizes)


	// training data
	train_X := [][]int {
		{1, 1, 1, 0, 0, 0},
		{1, 0, 1, 0, 0, 0},
		{1, 1, 1, 0, 0, 0},
		{0, 0, 1, 1, 1, 0},
		{0, 0, 1, 1, 0, 0},
		{0, 0, 1, 1, 1, 0},
	}

	train_Y := [][]int {
		{1, 0},
		{1, 0},
		{1, 0},
		{0, 1},
		{0, 1},
		{0, 1},
	}

	// construct DBN
	var dbn DBN
	DBN__construct(&dbn, train_N, n_ins, hidden_layer_sizes, n_outs, n_layers)

	// pretrain
	DBN_pretrain(&dbn, train_X, pretrain_lr, k, pretraining_epochs)

	// finetune
	DBN_finetune(&dbn, train_X, train_Y, fintune_lr, fintune_epochs)

	// test data
	test_X := [][]int {
		{1, 1, 0, 0, 0, 0},
		{1, 1, 1, 1, 0, 0},
		{0, 0, 0, 1, 1, 0},
		{0, 0, 1, 1, 1, 0},
	}

	test_Y := make([][]float64, test_N)
	for i := 0; i < test_N; i++ { test_Y[i] = make([]float64, n_outs)}

	// test
	for i := 0; i < test_N; i++ {
		DBN_predict(&dbn, test_X[i], test_Y[i])
		for j := 0; j < n_outs; j++ {
			fmt.Printf("%.5f ", test_Y[i][j])
		}
		fmt.Printf("\n")
	}
}



func main() {
	test_dbn()
}