;;;; -*- mode:scheme -*- ;;;;

;;; Neural Net Training.  Follows terminology in Winston Chap. 22.
;;; In particular:
;;; delta_j = o_j * (1 - o_j) * beta_j
;;; y_j = o_j
;;; y* = desired

;;; Sample training data Note that this is a bit different from the
;;; usual representation, so we need a simple conversion function.
;;; See neural-convert-sample below.

(define *xor-samples*  '(((1 0) (1))
			 ((0 1) (1))
			 ((0 0) (0))
			 ((1 1) (0))))

;;; Controls drawing
(define *neural-draw-interval-while-training* 500)

;;; Weights used during training
(define *neural-weights* '())
(define *neural-initial-weight-magnitude* 0.01)
(define *neural-deltas* '())

;;; High-Level Interface for two-class classification

;;; Rate Constant
(define *neural-rate* 1.0)
(define *neural-momentum* 0.1)
(define *neural-max-epochs* 1000)
(define *neural-min-rms* 0.01)

;; The actual names of the classes - will be converted to 0 and 1
(define *neural-class0* #f)			; set by neural-train
(define *neural-class1* #f)

(define (neural-train data . hidden)
  ;; a minimal network for binary classification, this gets called
  ;; with the structure of the hidden units.  If no arguments are
  ;; provided, then there are no hidden units.  If it gets called
  ;; with, for example, (neural-train *data* 3 2) then two layers of
  ;; hidden units will be created, the first with 3 units and the next
  ;; with 2 units.
  (let ((classes (map car (class-counts data))))
    (cond ((= (length classes) 2)
	   (cond ((or (equal? classes '(0 1))
		      (equal? classes '(0.0 1.0))
		      (equal? classes '(1 0))
		      (equal? classes '(1.0 0.0)))
		  (set! *neural-class0* 0)
		  (set! *neural-class1* 1))
		 (else
		  (set! *neural-class0* (first classes))
		  (set! *neural-class1* (second classes)))))
	  (else 
	   (error "Only two classes allowed, but we have" classes))))
  (backprop *neural-max-epochs*
	    (neural-convert-training-data data *neural-class0* *neural-class1*)
	    (initialize-weights
	     (append
	      (list (length (data-point-features (first data))))
	      hidden
	      '(1)))
	    *neural-min-rms*)
  'done)

(define (neural-continue-train data)
  (backprop *neural-max-epochs*
	    (neural-convert-training-data data *neural-class0* *neural-class1*)
	    *neural-weights*		; use the current weights
	    *neural-min-rms*))

(define (neural-classify data-point)
  ;; This makes a prediction of the class, it needs to convert to
  ;; neural format and map the 0-1 prediction into the output classes.
  (let* ((converted 
	  (neural-convert-sample (normalize-data-point data-point)
				 *neural-class1*))
	 (out (make-neural-prediction (first converted)))
	 (prediction (if (>= out 0.5) *neural-class1* *neural-class0*)))
    (display* "The prediction is " prediction "(" out ")"
	      ".  Correct is " (data-point-class data-point)
	      "."
	      )
    prediction))

(define (make-neural-prediction features)
  ;; return firt output.
  (first (compute-final-outputs 
	  (forward-propagate-output features *neural-weights*))))

;;; Convert usual training data to neural net representation

(define (neural-convert-training-data data class-0 class-1)
  (random-reorder 
   (map (lambda (x) (neural-convert-sample x class-1))
	(normalize-data data))))

(define (neural-convert-sample x class-1)
  (list (data-point-features x)
	(list (if (eq? (data-point-class x) class-1)
		  0.9 0.1))))

;;; Train

;;;  Purpose:	Train the net
;;;  Arguments:	Maximum number of epochs, samples, RMS to stop

(define (backprop epoch-limit samples initial-weights min-rms)
  (newline)
  (set! *neural-weights* initial-weights) ; will be changed during training
  (set! *neural-deltas* #f)		; keeps previous deltas
  (let ((sample-count (length samples))
	(drawing-samples
	 (map (lambda (s) (make-data-point (second s) (first s))) samples)))
    (do ((step 0 (+ 1 step))
	 (limit (* epoch-limit sample-count)))
	((cond ((= step limit) #t)
	       ((zero? (modulo step (* 10 sample-count)))
		(let ((rms-error (compute-rms-error samples)))
		  (print-average-error (/ step sample-count) rms-error)
		  (if (< rms-error min-rms) #t #f)))
	       (#t #f)))
      (if (and *draw* *neural-draw-interval-while-training*)
	  (if (zero? (remainder step (* sample-count *neural-draw-interval-while-training*)))
	      (if (= step 0)
		  (draw-classifier neural-get-point drawing-samples)
		  (draw-classifier neural-get-point drawing-samples *window*))))
      (let* ((sample (list-ref samples (modulo step sample-count)))
	     (sample-inputs (first sample))
	     (desired-outputs (second sample)))
	(single-step step sample-inputs desired-outputs)))
    (if *draw* (draw-classifier neural-get-point drawing-samples *window*))
    'done))

;;;  Purpose:	Perform a single training step
;;;  Arguments:	Step number
;;;  		Sample inputs
;;;  		Desired outputs for those inputs
;;;  Remarks:	Weights changed by side effect
(define (single-step step input desired-outputs)
  (let* ((outputs (forward-propagate-output input *neural-weights*))
	 (final-outputs (compute-final-outputs outputs))
	 (output-betas (output-layer-betas final-outputs
					   desired-outputs))
	 (betas (reverse (backward-propagate-beta (reverse outputs)
						  output-betas
						  (reverse *neural-weights*))))
	 (partials (invert-partials (compute-partials betas (cons input outputs)))))
    (if *neural-momentum*		; remember the deltas
	(set! *neural-deltas* partials))
    (set! *neural-weights* 
	  (add-by-layers *neural-weights* 
			 (if (and *neural-momentum* *neural-deltas*)
			     ;; momentum
			     (add-by-layers
			      (multiply-by-layers *neural-rate* partials)
			      (multiply-by-layers *neural-momentum* *neural-deltas*))
			     ;; no momentum
			     (multiply-by-layers *neural-rate* partials))
			 ))))

;;; Forward Propagation

;;;  Purpose:	Propagate values from inputs to outputs
;;;  Arguments:	Inputs to net
;;;  Returns:	Outputs of net
(define (forward-propagate-output inputs-to-k remaining-weights)
  (if (null? remaining-weights)
      '()
    (let* ((weights-to-k (first remaining-weights))
	   (output-k (forward-propagate-output-one-layer
		      (cons -1.0 inputs-to-k)
		      weights-to-k)))
      (cons output-k
	    (forward-propagate-output
	     output-k
	     (rest remaining-weights))))))

;;;  Arguments:	Outputs of leftward layer, weights to this layer
;;;  Returns:	Outputs of this layer
(define (forward-propagate-output-one-layer outputs-from-j weights-to-k)
  (map (lambda (weights-j-to-k)
	      (sigmoid (vector-dot-product outputs-from-j weights-j-to-k)))
	  weights-to-k))

;;;  Purpose:	Fetch final outputs from layer-by-layer outputs
;;;  Arguments:	Layer-by-layer outputs
;;;  Returns:	Final outputs
(define (compute-final-outputs outputs)
  (first (last-pair outputs)))

;;; Backward Propagation

;;;  Arguments:	Obvious from names
;;;  Returns:	Betas associated with output layer
(define (output-layer-betas actual-outputs desired-outputs)
  (vector-difference desired-outputs actual-outputs))

;;;  Arguments:	Weights to rightward layer
;;;  		Outputs of rightward layer
;;;   		Betas of rightward layer
;;;  Returns:	Betas of this layer
(define (backward-propagate-beta-one-layer weights-to-k outputs-k betas-k)
  (map (lambda (weights)
	      (reduce + 0.0 (map * 
				  weights
				  (map derivative-of-sigmoid outputs-k)
				  betas-k)))
	  (transform-weights weights-to-k)))

;;;  Purpose:	Propagate betas backward from outputs to inputs
;;;  Arguments:	Outputs of each layer, from outputs to first layer
;;;  		Betas for each layer
;;;  		Weights of each layer, from outputs to inputs
;;;  Returns:	All betas
(define (backward-propagate-beta reversed-outputs betas-k reversed-weights)
  (if (null? reversed-weights)
      '()
    (let* ((outputs-k (first reversed-outputs))
	   (weights-k (first reversed-weights))
	   (betas-j (backward-propagate-beta-one-layer weights-k
						       outputs-k
						       betas-k)))
      (cons betas-k
	    (backward-propagate-beta (rest reversed-outputs)
				     (rest betas-j)
				     (rest reversed-weights))))))
   
;;;  Purpose:	Convert left-to-right weight description to right-to-left
;;;  Arguments:	Left-to-right weight description
;;;  Returns:	Right-to-left weight description
;;;  Remarks:	Transforms a single layer
(define (transform-weights weights)
  (let ((result '()) (upper-bound (length (first weights))))
    (do ((n 0 (+ 1 n)))
	((= n upper-bound) (reverse result))
      (set! result (cons (map (lambda (x) (list-ref x n)) weights) result)))))


;;;  Purpose:	Convert right-to-left partials description to left-to-right
;;;  Arguments:	Right-to-left partials description
;;;  Returns:	Left-to-right partials description
;;;  Remarks:	Inverts all layers
(define (invert-partials layers)
  (map transform-weights layers))

;;; Compute Partials

;;;  Arguments:	Betas for this layer, outputs from this layer
;;;  Returns:	partials for all weights in this layer
(define (compute-partials-one-layer betas outputs)
  (let ((outputs-i (first outputs))
	(outputs-j (second outputs))
	(betas-j (first betas)))
    (map (lambda (oi)
	   (map (lambda (oj bj)
		  ;; This is deltaj
		  (* oi (derivative-of-sigmoid oj) bj))
		outputs-j
		betas-j))
	 (cons -1 outputs-i))))

;;;  Purpose:	Compute partials for all weights
;;;  Arguments:	All betas, all outputs
;;;  Returns:	All partials
(define (compute-partials betas outputs)
  (if (null? betas)
      '()
    (cons 
     (compute-partials-one-layer betas outputs)
     (compute-partials (rest betas) (rest outputs)))))

;;; Error Computation

;;;  Purpose:	Computs average error for all samples, 
;;;  		where error for one sample is rms error over all outputs
;;;  Arguments:	All samples
;;;  Returns:	Average rms error
(define (compute-rms-error samples)
  (/ (reduce + 0.0
	     (map
	      (lambda (sample) 
		  (vector-rms-differences
		   (compute-final-outputs
		     (forward-propagate-output (first sample) *neural-weights*))
		   (second sample)))
		  samples))
     (length samples)))

;;; Manipulate Weights

;;;  Purpose:	Takes two weight descriptions and adds corresponding elements
;;;  Arguments:	Typically, the current weights and changes
;;;  Returns:	Typically, new weights
(define (add-by-layers weights partials)
  (if (null? weights)
      '()
    (if (number? weights)
	(+ weights partials)
      (cons (add-by-layers (first weights) (first partials))
	    (add-by-layers (rest weights) (rest partials))))))

;;;  Purpose:	Mulitplies each weight by a multiplier
;;;  Arguments:	Typically, a rate and an expression containing partials
;;;  Returns:	Typically, changes to be made
(define (multiply-by-layers multiplier weights)
  (if (null? weights)
      '()
    (if (number? weights)
	(* multiplier weights)
      (cons (multiply-by-layers multiplier (first weights))
	    (multiply-by-layers multiplier (rest weights))))))

;;; Inform User

;;;  Purpose:	Generate formated progress information
;;;  Arguments:	Current epoch number, set of samples
(define (print-average-error epoch rms-error)
  (display* "Epochs: " epoch
	    " Average rms error: " rms-error))

;;; Initial Weight Computation

;;;  Purpose:	Produce initialized weights
;;;  Arguments:	A list of layer sizes, from inputs to outputs
;;;  Returns:	Initialized weights
(define (initialize-weights dimensions)
  (define (make-layer from-count to-count)
    (define (another-random-number counter)
      (if (zero? counter)
	  '()
	  (cons (- (random (* *neural-initial-weight-magnitude* 2)) 
		   *neural-initial-weight-magnitude*)
		(another-random-number (- counter 1)))))
    (if (zero? to-count)
	'()
	(cons (another-random-number (+ 1 from-count))
	      (make-layer from-count (- to-count 1)))))
  (if (null? (rest dimensions))
      '()
      (cons
       (make-layer (first dimensions) (second dimensions))
       (initialize-weights (rest dimensions)))))

;;; Vector Operations

(define (vector-dot-product v1 v2)
  (reduce + 0.0 (map * v1 v2)))

(define (vector-difference v1 v2)
  (map - v1 v2))

(define (vector-squared-differences v1 v2)
  (let ((v (vector-difference v1 v2)))
    (vector-dot-product v v)))

(define (vector-rms-differences v1 v2)
  (let ((v (vector-difference v1 v2)))
    (sqrt (/ (vector-dot-product v v)
	     (length v)))))

;;; Sigmoid Computations

;;;  Arguments:	Input of sigma function
;;;  Returns:	Output of sigmoid function
(define (sigmoid input)
  (if (> input 50.0)			; just clamp it to 1
      1.0
      (if (< input -50.0)		; just clamp it to 0
	  0.0
	  (/ 1.0 (+ 1 (exp (- input)))))))

;;;  Arguments:	Output of sigmoid function
;;;  Returns:	Derivitive of output of sigmoid function wrt output variable
(define (derivative-of-sigmoid output)
  (* output (- 1 output)))