;;;; -*- mode:scheme -*- ;;;;
;;;; Created: 7 Nov 95

;;;; UTILITIES

(load-option 'format)

(define (first x) (car x))

(define (second x) (cadr x))

(define (third x) (caddr x))

(define (rest x) (cdr x))

(define (reduce function initial-value items)
  "
  Purpose:	   Use function to combine list items with an initial value
  Example:         (reduce + 1 '(2 3 4))
  Value:           10
  "
  (if (null? items)
      initial-value
      (reduce function
              (function initial-value (first items)) ; New initial value
              (rest items))))			     ; New list of items

(define (nth n list-of-elements) (list-ref list-of-elements n))

;;; Test Data

(define *samples*  '(((1 0) (0.0))
		   ((0 1) (0.0))
		   ((0 0) (1.0))
		   ((1 1) (1.0))))

(define *initial-weights* 
      '(((0.255775571 -0.698110938 0.2516861)
         (0.3975036 0.906125069 0.360087872))
        ((0.129490018 -0.5424665 -0.582934))))

;;;fix bug from lack of (random)
(define (fix-by factor weight-thingy)
  (if (number? weight-thingy)
      (* factor weight-thingy)
      (if (null? weight-thingy)
	  '()
	  (cons (fix-by factor (car weight-thingy))
		(fix-by factor (cdr weight-thingy))))))
	  
(define *large-initial-weights*
  (fix-by 100 *initial-weights*))

(define *tiny-initial-weights*
  (fix-by 0.01 *initial-weights*))


;;; Rate Constant

(define *rate* 2.0)

;;; 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)))))

;;; Sigma Computations

(define (sigma input)
  "
  Arguments:	Input of sigma function
  Returns:	Output of sigma function
  "
  (if (> input 50.0)
      1.0
      (if (< input -50.0)
	  0.0
	  (/ 1.0 (+ 1 (exp (- input)))))))

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

;;; Forward Propagation

(define (forward-propagate-output inputs-to-k remaining-weights)
  "
  Purpose:	Propagate values from inputs to outputs
  Arguments:	Inputs to net
  Returns:	Outputs of net
  "
  (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))))))

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

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

;;; Backward Propagation

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

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

(define (backward-propagate-beta reversed-outputs betas-k reversed-weights)
  "
  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
  "
  (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))))))
   
(define (transform-weights 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
  "
  (let ((result '()) (upper-bound (length (first weights))))
    (do ((n 0 (+ 1 n)))
	((= n upper-bound) (reverse result))
      (set! result (cons (map (lambda (x) (nth n x)) weights) result)))))

(define (invert-partials layers)
  "
  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
  "
  (map transform-weights layers))

;;; Compute Partials

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

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

;;; Train

(define (train epoch-limit)
  "
  Purpose:	Train the net
  Arguments:	Maximum number of epochs
  "

  (format #t "~%")
  (let ((sample-count (length *samples*)))
    (do ((step 0 (+ 1 step))
	 (limit (* epoch-limit sample-count)))
	((cond ((= step limit) #t)
	       ((zero? (modulo step (* 25 sample-count)))
		(let ((rms-error (compute-rms-error *samples*)))
		  (print-average-error (/ step sample-count) rms-error)
		  (if (< rms-error 0.1) #t #f)))
	       (#t #f)))
      (let* ((sample (nth (modulo step sample-count) *samples*))
	     (sample-inputs (first sample))
	     (desired-outputs (second sample)))
	(single-step step sample-inputs desired-outputs)))))

(define (single-step step input desired-outputs)
  "
  Purpose:	Perform a single training step
  Arguments:	Step number
  		Sample inputs
  		Desired outputs for those inputs
  Remarks:	Weights changed by side effect
  "
  (let* ((outputs (forward-propagate-output input *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 *weights*))))
	 (partials (invert-partials (compute-partials betas (cons input outputs)))))
    (set! *weights* (add-by-layers *weights* 
				   (multiply-by-layers *rate* partials)))))

;;; Error Computation

(define (compute-rms-error samples)
  "
  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
  "
  (/ (reduce + 0.0
	     (map
	      (lambda (sample) 
		  (vector-rms-differences
		   (compute-final-outputs
		     (forward-propagate-output (first sample) *weights*))
		   (second sample)))
		  samples))
     (length samples)))

;;; Manipulate Weights

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

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

;;; Inform User

(define (print-average-error epoch rms-error)
  "
  Purpose:	Generate formated progress information
  Arguments:	Current epoch number, set of samples
  "
  (format #t "Epochs: ~a  Average rms error: ~a~%"
	  epoch
	  rms-error))

;;; Initial Weight Computation

(define *initial-weight-magnitude* 1.0)

(define (initialize-weights dimensions)
  "
  Purpose:	Produce initialized weights
  Arguments:	A list of layer sizes, from inputs to outputs
  Returns:	Initialized weights
  "
  (if (null? (rest dimensions))
      '()
    (cons
     (make-layer (first dimensions) (second dimensions))
     (initialize-weights (rest dimensions)))))

(define (make-layer from-count to-count)
  "
  Purpose:	Helper
  "
  (if (zero? to-count)
      '()
    (cons (another-random-number (+ 1 from-count))
	  (make-layer from-count (- to-count 1)))))

(define (another-random-number counter)
  "
  Purpose:	Helper
  Remark:       Random returns a random number between 0 and the argument;
                All such numbers are equally likely---that is, the 
                distribution is uniform
  "
  (if (zero? counter)
      '()
    (cons (- (random (* *initial-weight-magnitude* 2)) 
	     *initial-weight-magnitude*)
	  (another-random-number (- counter 1)))))