nn2.lisp

;;; -*- coding:utf-8; mode:lisp -*-

;; This is a secondary work for d.hatena.ne.jp/masatoi/20160106/1452057909.
;; See for the article regarding the copyright.
;;                                                  2015-01-09 @nfunato

;; オリジナルに対する本質的な変更は無く、
;; 実質的な変更点は、以下の手法で焦点でない細部を見えないようにしたこと
;;   - 数値ベクトルへのarefアクセスは、calc- という関数に分離した
;;   - 数値ベクトル以外のarefアクセスはstructure accessor経由にした
;;   - そこそこ高階の繰り返しオペレータを使うようにした
;;
;; 高階オペレータの一つとして、make-random-weight の書き換えに metatilities libraryの
;; maparray! を使った(といっても十分に一般的で簡単な定義だが)ので、
;; 下位関数の linearize-arrayと共に コピーなりインポートなりが必要
;; (https://github.com/gwkkwg/metatilities-base/blob/master/dev/l0-arrays.lisp)

(in-package :cl-user)
(defpackage wiz-nn (:use :cl))
(in-package :wiz-nn)

;;; Structures

(defstruct layer
  in-dim
  out-dim
  w-mat
  u-vec
  z-vec
  delta-vec
  activation-func
  activation-func-diff)

(defstruct nn
  n-of-layers
  layer-vec
  learning-rate)

(defun nn-nth-layer  (nn n) (aref (nn-layer-vec nn) n))
(defun nn-last-layer (nn)   (nn-nth-layer nn (1- (nn-n-of-layers nn))))

;;; Constructors

(defun rand-by-delta (delta)
  ;; return a random number between -delta and +delta
  (- (random (* delta 2)) delta))

(defun make-dfarray (dim)
  (make-array dim :element-type 'double-float :initial-element 0d0))

(defun make-random-weight (in-dim out-dim)
  ;; using maparray! in metatilities
  (maparray! (make-array (list out-dim in-dim) :element-type 'double-float)
             (lambda (i) (declare (ignore i)) (rand-by-delta 0.1d0))))

(defun make-random-layer (in-dim out-dim activation-func activation-func-diff)
  (make-layer :in-dim  in-dim
              :out-dim out-dim
              :w-mat (make-random-weight in-dim out-dim)
              :u-vec (make-dfarray out-dim)
              :z-vec (make-dfarray out-dim)
              :delta-vec (make-dfarray out-dim)
              :activation-func activation-func
              :activation-func-diff activation-func-diff))

(defun make-random-nn (dim-list act-fn-pair-list &optional (rate 0.01d0))
  (assert (= (1- (length dim-list)) (length act-fn-pair-list)))
  (let ((layers (loop for (in-d out-d)   on dim-list
                      for (afn afn-diff) in act-fn-pair-list
                      collect (make-random-layer in-d out-d afn afn-diff))))
    (make-nn :n-of-layers  (1- (length dim-list))
             :layer-vec    (coerce layers 'vector)
             :learning-rate rate)))

;;; Activation functions

;; RLF; Rectified Linear Function
(defun RLF (u)
  (if (> u 0d0) u 0d0))

(defun RLF-diff (u)
  (if (>= u 0d0) 1d0 0d0))

;; Identical function
;; Differntial of identity function
(defun one (x)
  (declare (ignore x))
  1d0)

;; Logistic function
(defun logistic (u)
  (/ 1d0 (+ 1d0 (exp (- u)))))

(defun logistic-diff (u)
  (let ((f (logistic u)))
    (* f (- 1d0 f))))

;; Hyperbolic tangent
(defun tanh-diff (u)
  (let ((tanh-u (tanh u)))
    (- 1d0 (* tanh-u tanh-u))))

;;; Feed-forward

(defun calc-u-vec (in-vec layer)
  (loop for j from 0 below (layer-out-dim layer) do
    (setf (aref (layer-u-vec layer) j)
          (loop for i from 0 below (layer-in-dim layer)
                summing (* (aref in-vec i) (aref (layer-w-mat layer) j i)))))
  (layer-u-vec layer))

(defun calc-z-vec (layer)
  ;; the dimensions of z-vec and u-vec are equal by definition
  (map-into (layer-z-vec layer)
            (layer-activation-func layer) (layer-u-vec layer))
  (layer-z-vec layer))

(defun calc-u-vec/z-vec (in-vec layer)
  (calc-u-vec in-vec layer)
  (calc-z-vec layer))

(defun forward (in-vec nn)
  (loop for vec   =      in-vec then (layer-z-vec layer)
        for layer across (nn-layer-vec nn)
        do (calc-u-vec/z-vec vec layer)))

;;; Back-propagation

(defun calc-delta-vec (curr next)
  (loop for j from 0 below (layer-in-dim next) do
    (setf (aref (layer-delta-vec curr) j)
          (* (funcall (layer-activation-func-diff curr)
                      (aref (layer-u-vec curr) j))
             (loop for k from 0 below (layer-out-dim next)
                   summing
                   (* (aref (layer-delta-vec next) k)
                      (aref (layer-w-mat     next) k j)))))))

(defun backward (train-vec nn)
  (let ((last-layer (nn-last-layer nn)))
    ;; calculate last layer's delta
    (map-into (layer-delta-vec last-layer)
              #'- (layer-z-vec last-layer) train-vec)
    ;; calculate other deltas
    (loop for l from (- (nn-n-of-layers nn) 2) downto 0
          for curr = (nn-nth-layer nn l)
          for next = (nn-nth-layer nn (1+ l))
          do (calc-delta-vec curr next))))

(defun predict (in-vec nn)
  (forward in-vec nn)
  (layer-z-vec (nn-last-layer nn)))

(defun calc-w-mat (layer prev-layer-z-vec learning-rate)
  (loop for i from 0 below (layer-in-dim layer) do
    (loop for j from 0 below (layer-out-dim layer) do
      (setf (aref (layer-w-mat layer) j i)
            (- (aref (layer-w-mat layer) j i)
               (* learning-rate
                  (aref prev-layer-z-vec        i)
                  (aref (layer-delta-vec layer) j)))))))

(defun update (in-vec train-vec nn)
  (forward in-vec nn)
  (backward train-vec nn)
  ;; update first layer
  (let ((rate  (nn-learning-rate nn))
        (first (nn-nth-layer nn 0)))
    (calc-w-mat first in-vec rate)
    ;; update other layer
    (loop for l      from 1 below (nn-n-of-layers nn)
          for prev = (nn-nth-layer nn (1- l))
          for curr = (nn-nth-layer nn l)
          do (calc-w-mat curr (layer-z-vec prev) rate))))

nn3.lisp

;;; -*- coding:utf-8; mode:lisp -*-

;;; This is a secondary work for d.hatena.ne.jp/masatoi/20160106/1452057909.
;;; See for the article regarding the copyright.
;;;                                                 2015-01-11 by @nfunato

(in-package :cl-user)
(defpackage wiz-nn (:use :cl))
(in-package :wiz-nn)

;;; General Utilities

(defun rand-by-delta (delta)            ; return a random number in [-d,+d)
  (- (random (* delta 2)) delta))

(defun mat-dims (mat) (array-dimensions mat))
(defun vec-dim  (vec) (car (array-dimensions vec)))

(defun array-map1! (dest-arr fn arr)    ; dimension independent
  (loop for i from 0 below (array-total-size arr) do
    (setf (row-major-aref dest-arr i) (funcall fn (row-major-aref arr i))))
  dest-arr)

(defun apply-mat! (result-mat fn mat row-in-vec column-in-vec)
  #+:test
  (assert (and (equal (mat-dims mat) (mat-dims result-mat))
               (equal (mat-dims mat)
                      (list (vec-dim column-in-vec) (vec-dim row-in-vec)))))
  (loop for j from 0 for jv across column-in-vec do
    (loop for i from 0 for iv across row-in-vec do
      (let ((mv (aref mat j i)))
        (setf (aref result-mat j i) (funcall fn mv iv jv)))))
  result-mat)

(defun mult-mat-vec! (ovec imat ivec &optional transpose-p)
  #+:test
  (assert (equal (mat-dims imat)
                 (if transpose-p
                     (list (vec-dim ivec) (vec-dim ovec))
                     (list (vec-dim ovec) (vec-dim ivec)))))
  (flet ((inner-product (j)
           (loop for vv across ivec
                 for i  from 0
                 ;; if transpose-p, slicing is done for vertical direction
                 for mv = (if transpose-p (aref imat i j) (aref imat j i))
                 sum (* mv vv))))
    (loop for j from 0 below (length ovec) do
      (setf (aref ovec j) (inner-product j)))
    ovec))

;;; Data Structures

;; dump/undump し易いように、dimensionsは 2重持ちせず、weighting-matrix だけに留める
(defclass layer ()
  ((weighting-matrix :reader w-mat :initarg :w-mat)
   (u-vec :reader u-vec :initarg :u-vec)
   (z-vec :reader z-vec :initarg :z-vec)
   (delta-vec :reader d-vec :initarg :d-vec)
   (activation-fn :reader act-fn :initarg :act-fn)
   (activation-differential-fn
    :reader act-diff-fn :initarg :act-diff-fn)))

(defmethod in-dim  ((w-mat array)) (cadr (array-dimensions w-mat)))
(defmethod out-dim ((w-mat array)) (car  (array-dimensions w-mat)))
(defmethod in-dim  ((l layer))     (in-dim  (w-mat l)))
(defmethod out-dim ((l layer))     (out-dim (w-mat l)))

(defun make-df-vector (dim)
  (make-array dim :element-type 'double-float :initial-element 0d0))

(defun make-random-df-matrix (in-dim out-dim)
  (flet ((fn (x) (declare (ignore x)) (rand-by-delta 0.1d0)))
    (let ((mat (make-array (list out-dim in-dim) :element-type 'double-float)))
      (array-map1! mat #'fn mat))))

(defun make-layer (w-mat act-fn act-diff-fn)
  (let ((out-dim (out-dim w-mat)))
    (make-instance 'layer
                   :w-mat w-mat
                   :act-fn act-fn
                   :act-diff-fn act-diff-fn
                   :u-vec (make-df-vector out-dim)
                   :z-vec (make-df-vector out-dim)
                   :d-vec (make-df-vector out-dim))))

(defun make-random-weight-layer (in-dim out-dim act-fn act-diff-fn)
  (make-layer (make-random-df-matrix in-dim out-dim) act-fn act-diff-fn))

;; ネットワーク全体はレイヤーの(ベクトルでなく)リストということにする。
;; (後で出てくるが、loop 文を on で使いたいので)

(defclass neural-net ()
  ((layer-list :reader layer-list :initarg :layer-list)
   (learning-rate :reader learning-rate :initarg :learning-rate)))

(defmethod last-layer ((nn neural-net)) (car (last (layer-list nn))))

(defun make-neural-net (layers learning-rate)
  (make-instance 'neural-net :layer-list layers :learning-rate learning-rate))

;;; Top-level Constructor

(defun make-random-neural-net (dim-list fn-pair-list &optional (rate 0.01d0))
  (assert (= (1- (length dim-list)) (length fn-pair-list)))
  (let ((layers (loop for (in-d out-d) on dim-list
                      for (f diff-f)   in fn-pair-list
                      collect (make-random-weight-layer in-d out-d f diff-f))))
    (make-neural-net layers rate)))

;;; Activation Functions

;; RLF; Rectified Linear Function
(defun RLF (u)
  (if (> u 0d0) u 0d0))

(defun RLF-diff (u)
  (if (>= u 0d0) 1d0 0d0))

;; Identical function
;; Differntial of identity function
(defun one (x)
  (declare (ignore x))
  1d0)

;; Logistic function
(defun logistic (u)
  (/ 1d0 (+ 1d0 (exp (- u)))))

(defun logistic-diff (u)
  (let ((f (logistic u)))
    (* f (- 1d0 f))))

;; Hyperbolic tangent
(defun tanh-diff (u)
  (let ((tanh-u (tanh u)))
    (- 1d0 (* tanh-u tanh-u))))

;;; Feed-forward

(defmethod calc-u-vec/z-vec ((curr layer) (prev-z-vec vector))
  (let ((u-vec (u-vec curr))
        (z-vec (z-vec curr)))
    (mult-mat-vec! u-vec (w-mat curr) prev-z-vec)
    (map-into z-vec (act-fn curr) u-vec)))

(defmethod calc-u-vec/z-vec ((curr layer) (prev layer))
  (calc-u-vec/z-vec curr (z-vec prev)))

(defmethod forward ((nn neural-net) in-vec)
  (loop for (prev curr) on (cons in-vec (layer-list nn))
        while curr do (calc-u-vec/z-vec curr prev)))

;;; Back-propagation

;; when curr is not last layer
(defmethod calc-delta-vec ((curr layer) (next layer))
  (flet ((differentiate (d u) (* d (funcall (act-diff-fn curr) u))))
    (let ((d-vec (d-vec curr)))
      ;; currently we accept looping twice, i.e. mult-mat-vec! and map-into
      (mult-mat-vec! d-vec (w-mat next) (d-vec next) t)
      (map-into d-vec #'differentiate d-vec (u-vec curr)))))

;; only when curr is last layer
(defmethod calc-delta-vec ((last layer) (train-vec vector))
  (map-into (d-vec last) #'- (z-vec last) train-vec))

(defmethod backward ((nn neural-net) train-vec)
  (loop for (next curr) on (cons train-vec (reverse (layer-list nn)))
        while curr do (calc-delta-vec curr next)))

;;; Learning (update) and Predicting

(defmethod calc-w-mat ((curr layer) (prev-z-vec vector) learning-rate)
  (flet ((fn (mat-v pz d) (- mat-v (* learning-rate pz d))))
    (let ((w-mat (w-mat curr)))
      (apply-mat! w-mat #'fn w-mat prev-z-vec (d-vec curr)))))

(defmethod calc-w-mat ((curr layer) (prev layer) learning-rate)
  (calc-w-mat curr (z-vec prev) learning-rate))

(defmethod update ((nn neural-net) in-vec train-vec)
  (forward nn in-vec)
  (backward nn train-vec)
  (loop for (prev curr) on (cons in-vec (layer-list nn))
        while curr do (calc-w-mat curr prev (learning-rate nn))))

(defun predict (nn in-vec)
  (forward nn in-vec)
  (z-vec (last-layer nn)))

;;; Sample Usage

(defun make-input-data ()
  (loop repeat 100
        ;; vector の 第2要素(1d0)は バイアス
        collect (vector (rand-by-delta pi) 1d0)))

(defun make-train-data (input-data)
  (mapcar (lambda (x) (vector (sin (aref x 0)))) input-data))

(defvar *input-data* (make-input-data))
(defvar *train-data* (make-train-data *input-data*))

(defun make-sample-nn2 ()
  (make-random-neural-net
    '(2 50 1)
    (list (list #'RLF #'RLF-diff)
          (list #'identity #'one))
    0.05d0))

(defun make-sample-nn3 ()
  (make-random-neural-net
    '(2 20 10 1)
    (list (list #'logistic #'logistic-diff)
          (list #'RLF #'RLF-diff)
          (list #'identity #'one))
    0.05d0))

(defparameter *nn2* (make-sample-nn2))
(defparameter *nn3* (make-sample-nn3))

(defun learn (nn &optional (cycle 1000))
  (dotimes (i cycle)
    (mapc (lambda (in out) (update nn in out))
          *input-data*
          *train-data*))
  nil)

(defvar *saved-data* nil)

(defmethod dump-nn-w-mats (nn)
  (mapcar #'w-mat (layer-list nn)))

(defmethod undump-nn-w-mats (nn w-mats)
  (flet ((reinitialize-layer (layer w-mat)
           (let ((out-dim (out-dim w-mat)))
             (with-slots (weighting-matrix u-vec z-vec delta-vec) layer
               (setf weighting-matrix w-mat
                     u-vec (make-df-vector out-dim)
                     z-vec (make-df-vector out-dim)
                     delta-vec (make-df-vector out-dim))))))
    (loop for l  in (layer-list nn)
          for wm in w-mats
          do (reinitialize-layer l wm))))

(defun plot-data (nn)
  (loop for x from -3.14 to 3.14 by 0.01
        collect (aref (predict nn (vector (* x 1d0) 1d0))
                      0)))

;;; Debugging stuff

;; see http://power.mech.eng.osaka-cu.ac.jp/~takiyama/gnuplot-howto.html
;; in short,
;;   > set datafile separator ","
;;   > plot 'nn_tmp.csv'
(defun fileout-plot-data (nn &optional (path #p"./nn_tmp.csv"))
  (with-open-file (st path :direction :output :if-exists :supersede)
    (loop for x from -3.14 to 3.14 by 0.01
          for y in (plot-data nn)
          do (format st "~a,~a~%" x y))))

(defvar *saved-w-mats*
  nil
  )

(defvar *saved-input*
  nil
  )

(defun undump-data (nn &optional (w-mats *saved-w-mats*) (input *saved-input*))
  (undump-nn-w-mats nn w-mats)
  (setf *input-data* input)
  (setf *train-data* (make-train-data input))
  nil
  )