Skip to content

Instantly share code, notes, and snippets.

@tanakahx tanakahx/bm.lisp
Last active Sep 27, 2015

Embed
What would you like to do?
Simple Boltzmann Machine Simulation
(in-package :cl-user)
(eval-when (:compile-toplevel :load-toplevel :execute)
(ql:quickload '(:cl-num-utils :lla :array-operations) :silent t))
(defpackage :bm
(:use :cl)
(:import-from :clnu
:e*
:transpose)
(:import-from :clnu.mx
:mx)
(:import-from :lla
:mm)
(:import-from :aops
:flatten
:reshape
:reshape-col
:nrow
:ncol
:sub)
(:export :bm))
(in-package :bm)
(defparameter *node-count* 3)
(defparameter *ws* (clnu.mx:mx t
( 0 -2 -2)
(-2 0 -2)
(-2 -2 0)))
(defparameter *bs* (clnu.mx:mx t
(-5)
(-5)
(-5)))
(defparameter *temperature* 1)
(defparameter *update-count* 100)
(defun sigmoid (s temperature)
(/ 1 (+ 1 (exp (/ (- s) temperature)))))
(defun prob-sigmoid (s temperature)
"シグモイド関数を使った確率的2値モデル。"
(if (< (random 1d0) (sigmoid s temperature)) 1 0))
(defun append1 (lst x)
(append lst (list x)))
(defun power-set (n)
"n bitの状態変数の集合をリストとして生成する。"
(labels ((rec (n acc)
(if (> n 0)
(rec (1- n)
(append (mapcar #'(lambda (x) (append1 x 0)) acc)
(mapcar #'(lambda (x) (append1 x 1)) acc)))
acc)))
(rec n (list nil))))
(defun power-set-array (n)
"n bitの状態変数の集合をベクタとして生成する。"
(mapcar #'(lambda (x) (aops:reshape-col (coerce x 'vector))) (power-set n)))
(defun energy (ws bs xs)
"エネルギー関数。状態変数の存在確率の理論値を求めるために使用。"
(+ (reduce #'+ (aops:flatten (clnu:e* -1/2 ws (lla:mm xs (clnu:transpose xs)))))
(aref (lla:mm (clnu:transpose bs) xs) 0 0)))
(defun softmax (ws bs xs temperature)
"重みws, 閾値bs, 温度temperatureにおける状態変数xsの存在確率を求める。"
(loop for s in (power-set-array (aops:nrow xs))
sum (exp (- (/ (energy ws bs s) temperature))) into total
finally (return (/ (exp (- (/ (energy ws bs xs) temperature))) total))))
(defun weighted-sum (ws xs bs row)
(- (lla:mm (aops:sub ws row) (aops:flatten xs)) (aref bs row 0)))
(defun update-nodes! (ws bs xs temperature)
"状態変数xsを重みws, 閾値bsにより1要素ずつ順番に更新する。xsは変更されることに注意。"
(loop for i below (aops:nrow xs)
for s = (weighted-sum ws xs bs i)
for x = (prob-sigmoid s temperature)
do (setf (aref xs i 0) x)
finally (return xs)))
(defun make-nodes (n)
"0/1にランダムに初期化した状態変数を作成する。"
(let ((xs (make-array (list n 1))))
(dotimes (i n)
(setf (aref xs i 0) (random 2)))
xs))
(defun make-hist (n)
(make-array (expt 2 n)))
(defun hist-count (hist bin)
"ヒストグラムのbinをインクリメントする。"
(incf (aref hist bin)))
(defun vector-to-number (xs &key (radix 10))
"vectorをradix進数の数値として解釈する。状態変数をヒストグラムのbinにエンコードするために使用する。"
(parse-integer (format nil "~{~a~}"
(let ((*print-base* radix))
(map 'list #'princ-to-string
(reverse (aops:flatten xs)))))
:radix radix))
(defun bm (ws bs &key (update-count *update-count*) (temperature *temperature*))
"ボルツマンマシンをシミュレートして状態の存在確率の理論値と実測値を表示する。
エネルギー関数から求めた重みws, 閾値bsを与える。"
(assert (= (aops:ncol ws) (aops:nrow bs)))
(let* ((node-count (aops:nrow ws))
(xs (make-nodes node-count))
(hist (make-hist node-count)))
(loop repeat update-count
for bin = (vector-to-number (update-nodes! ws bs xs temperature) :radix 2)
do (hist-count hist bin))
(let ((expected (mapcar #'(lambda (xs) (softmax ws bs xs temperature)) (power-set-array 3)))
(actual (map 'list #'(lambda (x) (float (/ x update-count))) hist)))
(format t "理論値: ~{~1,3f~^ ~}~%" expected)
(format t "実測値: ~{~1,3f~^ ~}~%" actual))))
#+ (or) (bm:bm (clnu.mx:mx t (0 -2 -2) (-2 0 -2) (-2 -2 0))
(clnu.mx:mx t (-5) (-5) (-5))
:update-count 1000
:temperature 1)
#+ (or) (bm:bm (clnu.mx:mx t (0 -2 -2) (-2 0 -2) (-2 -2 0))
(clnu.mx:mx t (-1) (-1) (-1))
:update-count 1000
:temperature 1)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.