;;;; -*- Scheme -*-
;;
;;  Copyright (c) 1997 Massachusetts Institute of Technology
;;  
;;  This material was developed by the Amorphous Computing project at the
;;  Massachusetts Institute of Technology Artificial Intelligence
;;  Laboratory.  Permission to copy this software, to redistribute it, and
;;  to use it for any purpose is granted, subject to the following
;;  restrictions and understandings.
;;  
;;  1. Any copy made of this software must include this copyright notice
;;  in full.
;;  
;;  2. Users of this software agree to make their best efforts (a) to
;;  return to the MIT Amorphous Computing project any improvements or
;;  extensions that they make, so that these may be included in future
;;  releases; and (b) to inform MIT of noteworthy uses of this software.
;;  
;;  3. All materials developed as a consequence of the use of this
;;  software shall duly acknowledge such use, in accordance with the usual
;;  standards of acknowledging credit in academic research.
;;  
;;  4. MIT has made no warrantee or representation that the operation of
;;  this software will be error-free, and MIT is under no obligation to
;;  provide any services, by way of maintenance, update, or otherwise.
;;  
;;  5. In conjunction with products arising from the use of this material,
;;  there shall be no use of the name of the Massachusetts Institute of
;;  Technology nor of any adaptation thereof in any advertising,
;;  promotional, or sales literature without prior written consent from
;;  MIT in each case.


;; Interface and Utilities for HLSIM
;;-----------------------------------
;; Interface Functions for making a simulation object
;; Simulation object remembers the communication radius (simulation.get sim 'radius)
;; Different types of layouts (random, smoothed, hex, rectangular)
;; Functions for saving and reusing the same layout (for experiments, debugging)
;; Functions for Display management, color-me routines, and drawing list of processors
;; Utilities such as nbrhood size -> radius conversions, find-nearest-processor,
;; compute-nbrhood-stats, pi, etc
;; Many functions copied from dder's utils file
;;---------------------------------------------------------------

(declare (usual-integrations))

;; Generic Simulation Creation
;;-----------------------------

(define (make-random-sim n r)
  (let ((x (make-initialized-vector n (lambda (i) i (random 1.0))))
	(y (make-initialized-vector n (lambda (i) i (random 1.0))))
	(z (make-vector n 0.0)))
    (return-sim-with-layout x y z n r)))

(define (return-sim-with-layout x y z n r)
  (let ((sim (%make-simulation)))
    (do ((i 0 (+ i 1))) ((= i n)) (simulation.add-processor sim))
    (set-simulation.x! sim x)
    (set-simulation.y! sim y)
    (set-simulation.z! sim z)
    (set-simulation.neighbours! sim (point-neighbours x y z r))
    (simulation.define sim 'radius r)
    sim))

;; Save Current layout 
;;--------------------
; For reusing the same layout with different programs,
; Allows experiments using the same processor locations, 
; Also allows changing the radius for the same layout (nbrhoods are recalculated)
; note: fasdump/fasload does not work on sim

(define (save-layout sim)
  (list (vector-copy (simulation.x sim))
	(vector-copy (simulation.y sim))
	(vector-copy (simulation.z sim))))

; create a simulation using a saved layout
(define (make-new-sim-with-layout layout r)
  (let ((x (vector-copy (first layout)))
	(y (vector-copy (second layout)))
	(z (vector-copy (third layout)))
	(n (vector-length (first layout))))
    (return-sim-with-layout x y z n r)))

;; Simulation Display
;;--------------------

(define (simulation.display! sim flag)
  (set-simulation.display!
   sim
   (cond ((not flag)  #F)
	 ((graphics-device? flag) flag)
	 (else (make-display)))))

; can use to create a base window
(define (make-display)
  (let ((g  (make-graphics-device)))
    (graphics-operation g 'set-foreground-color "black")
    (graphics-operation g 'set-background-color "white")
    (graphics-set-coordinate-limits g -.01 -.01 1.01 1.01)
    (graphics-clear g)
    g))

; simple color-me function      
(define (%color-me sim i color)
  ; (define d 0.007)
  ; instead use a size proportional to the number of processors

  (define (draw color)
    (let ((g (simulation.display sim))
	  (x (vector-ref (simulation.x sim) i))
	  (y (vector-ref (simulation.y sim) i))
	  (d (* .5 (/ 1 (sqrt (simulation.n-processors sim)))) ))
      (graphics-operation g 'set-foreground-color color)
      (graphics-operation g 'fill-circle x y d)))

  (cond ((string? color)
	 (draw color))
	((number? color)
	 (draw (spectrum-color color)))
	(else 
	 (write-line (list "error color is not a string or number" color))
	 )))

(define (spectrum-color i)
  ; red to black
  (let ((spectrum-vec 
	 '#("#F02020" "#E02020" "#C02020" "#A02020" "#902020" "#202020")))
    (if (< i (vector-length spectrum-vec))
	(if (>= i 0)
	    (vector-ref spectrum-vec i)
	    (vector-ref spectrum-vec 0))
	(vector-ref spectrum-vec (- (vector-length spectrum-vec) 1))
	)
    ))

; some helpful coloring functions

(define (draw-processors sim proc-list color)
  (do ((i 0 (+ i 1))) ((= i (simulation.n-processors sim)))
    (if (member i proc-list)
	(processor.eval (vector-ref (simulation.processors sim) i)
			`(color-me ,color))))
  'done)

(define (color-neighbours sim i color)
  (draw-processors sim 
	(vector->list (vector-ref (simulation.neighbours sim) i)) color))

(define (draw-comms-circle sim i color)
  (graphics-operation (simulation.display sim) 'set-foreground-color color)
  (graphics-operation (simulation.display sim) 'draw-circle
		      (vector-ref (simulation.x sim) i)
		      (vector-ref (simulation.y sim) i)
		      (simulation.get sim 'radius)))

;; Different Layouts
;;-------------------

; RECTANGULAR
(define (make-rect-sim n r)  (make-jiggly-sim n r 0))

; PSEUDO-RANDOM (JIGGLY)
; Make jiggly simulation.  The j argument is the jigglyness, 0 <= j < 1
; Default j=0.5 NB If n is not a square number, n processors will be created,
; but there will be gaps in the grid near the top right. (dder)

(define (make-jiggly-sim n r #!optional j)
  (let ((j (if (default-object? j) 0.5 j))
	(d (ceiling (sqrt n))) ; number of divisions on axes of unit square
	(x (make-vector n))
	(y (make-vector n))
	(z (make-vector n 0.0)))
    (let loop ((p-i 0) (x-i 0) (y-i 0))
      (cond ( (= p-i n) 'done )
	    ( (= x-i d) (loop p-i 0 (1+ y-i)) )
	    ( else (vector-set! x p-i (/ (+ x-i (* j (random 1.0))) d))
		   (vector-set! y p-i (/ (+ y-i (* j (random 1.0))) d))
		   (loop (1+ p-i) (1+ x-i) y-i) )))
    ; remaining code is identical to make-sim/1
    (return-sim-with-layout x y z n r)))

; HEXAGONAL
; #columns= #rows, therefore the result is not a unit square surface
; columns = xdivisions = sqrt(total), therefore to avoid gaps in upper right
; corner choose n to be a perfect square.

(define (hex-total->columns total)
  (ceiling->exact (sqrt total)))

(define (hex-columns->u cols)
  (/ 1 (+ cols .5)))

(define (hex-columns->rowheight cols)
  (* (sqrt 3) (hex-columns->u cols) .5))

(define (make-hex-sim n r)
  (let* ((cols (hex-total->columns n))
	 (xunit (hex-columns->u cols))
	 (yunit (hex-columns->rowheight cols))
	 (offset (/ xunit 2))
	 ; return x,y,z vectors
	 (x (make-vector n))
	 (y (make-vector n))
	 (z (make-vector n 0.0))
	 )
    ;debug
    (write-line (list 'unit-distance xunit))

    (let loop ((p-i 0) (x-i 0) (y-i 0))
      (cond ; completed all processors
	    ((= p-i n) 'done )
	    ; end of a row, start new row
	    ((= x-i cols) 
	     (loop p-i 0 (1+ y-i)))
	    ; otherwise assign x,y 
	    (else (if (even? y-i)	; even & odd rows are offset
		      (vector-set! x p-i (exact->inexact (* x-i xunit)))
		      (vector-set! x p-i (+ offset (exact->inexact (* x-i xunit)))))
		  (vector-set! y p-i (exact->inexact (* y-i yunit)))
		  (loop (1+ p-i) (1+ x-i) y-i)
		  )))

    ; remaining code is identical to make-sim/1
    (return-sim-with-layout x y z n r)))

; Different version of hex
; specify the number of processors you want along the x axis (number of columns)
; creates enough rows to fill the unit square
(define (hex-nc->approx-total nc) (* (/ 2 (sqrt 3)) (square nc)))

(define (make-hex-square-sim nc r)
    (let* ((xunit  (exact->inexact (/ 1 (- nc .5))))
	   (offset (/ xunit 2.0))
	   (yunit  (* (/ (sqrt 3) 2) xunit))
	   (nr     (inexact->exact (ceiling (/ 1 yunit))))
	   (n      (* nc nr))
	   ; coordinate vectors
	   (x (make-vector n))
	   (y (make-vector n))
	   (z (make-vector n 0.0))
	   )

      (write-line (list "total processors" n))
      (do ((row 0 (+ row 1))) ((= row nr))
	(let ((ycoord (* row yunit)) (p (* row nc)))
	  ; place row
	  (do ((i 0 (+ i 1))) ((= i nc))
	    (vector-set! y (+ p i) ycoord)
	    (vector-set! x (+ p i) (+ (if (even? row) offset 0.0)
				      (exact->inexact (* i xunit))))
	    )
	))

      ; remaining code is identical to make-sim/1
      (return-sim-with-layout x y z n r)))


; SMOOTHED
; Creates a smoothed 2D layout First compute a reasonable distance.
; Then place points one by one, and check to make sure it doesn't
; conflict with other points.  For reasonable dist, placing N points
; on a square grid would require 1/sqrt(N) spacing. Some fraction of
; that (default 50%) is reasonable (dder)
;
; very slow O(n^2), rewrote to use spatial grids (rad)

(define (square x) (* x x))

(define (make-smooth-sim n r #!optional f)
  (let ((mindist (square (* (if (default-object? f) 0.5 f)
			    (/ 1 (sqrt n)))))
	(vx (make-vector n 0.0))
	(vy (make-vector n 0.0))
	(vz (make-vector n 0.0)))

    ; Spatial Grid (2D)
    (let* ((grid-divs (inexact->exact (floor (sqrt n))))
	   (grid-buckets (make-initialized-vector grid-divs (lambda (i) 
			    (make-initialized-vector grid-divs (lambda (j) '()))))))

      (define (get-bucket i j) (vector-ref (vector-ref grid-buckets i) j))
      (define (set-bucket i j val) (vector-set! (vector-ref grid-buckets i) j val))

      (define (grid-row x y) (inexact->exact (floor (* y grid-divs))))
      (define (grid-col x y) (inexact->exact (floor (* x grid-divs))))

      (define (insert-into-grid x y)
	(set-bucket (grid-row x y) (grid-col x y)
		    (cons (cons x y) (get-bucket (grid-row x y) (grid-col x y)))))

      (define (compute-dist x y pt)
	(+ (square (- x (car pt))) (square (- y (cdr pt)))))

      ; Check point against all other pts in the bucket
      ; return #f if point is too close to some point
      (define (check-grid-bucket x y r c)
	; check bucket range and return #t if no such bucket
	(if (or (< r 0) (>= r grid-divs) (< c 0) (>= c grid-divs))
	    #t
	    (let checklp ((ptlist (get-bucket r c)))
	      (if (null? ptlist) #t
		  (if (< (compute-dist x y (first ptlist)) mindist)
		      #f
		      (checklp (cdr ptlist)) ))
	      )))

      ; Point must check out against all nbring grid buckets
      (define (check-point x y)
	(let ((r (grid-row x y)) (c (grid-col x y)))
	  (and (check-grid-bucket x y (- r 1) (- c 1)) 
	       (check-grid-bucket x y (- r 1) c)
	       (check-grid-bucket x y (- r 1) (+ c 1))
	       (check-grid-bucket x y r (- c 1)) 
	       (check-grid-bucket x y r c)
	       (check-grid-bucket x y r (+ c 1))
	       (check-grid-bucket x y (+ r 1) (- c 1)) 
	       (check-grid-bucket x y (+ r 1) c)
	       (check-grid-bucket x y (+ r 1) (+ c 1))
	     )))

      ; Assign points one at a time, and check with previous points
      (write-line `(Creating smooth layout: mindist ,(sqrt mindist)))

      (let lp ((i 0) (c 0))
	(if (>= i n) 
	    'done
	    (let ((x0 (random 1.0)) (y0 (random 1.0)))
	      (vector-set! vx i x0) 
	      (vector-set! vy i y0)
	      (if (check-point x0 y0)
		  (begin 
		    (if (> c 0) (write-line (list "repeated node" i c "times")))
		    (insert-into-grid x0 y0) (lp (+ i 1) 0))
		  ; try to place the point 50 times before giving up
		  ; c keeps track of the number of trys
		  (if (< c 50)
		      (lp i (+ c 1))
		      (begin (insert-into-grid x0 y0) (lp (+ i 1) 0)) ))
	      )))
      ) ; end spatial grid let

    ; remaining code is identical to make-sim/1
    (return-sim-with-layout vx vy vz n r)))

#|
(define (make-smooth-sim n r #!optional f)
  (let ((mindist (square (* (if (default-object? f) 0.5 f)
			    (/ 1 (sqrt n)))))
	(x (make-vector n 0.0))
	(y (make-vector n 0.0))
	(z (make-vector n 0.0)))

    (define (compute-dist i j)
      (+ (square (- (vector-ref x i) (vector-ref x j)))
	 (square (- (vector-ref y i) (vector-ref y j)))
	 (square (- (vector-ref z i) (vector-ref z j)))))

    ; Check point against all other points
    (define (check-point i)
      (let checklp ((j 0))
	(if (>= j i) 
	    #t
	    (if (< (compute-dist i j) mindist)
		#f
		(checklp (+ j 1))))))

    ; Assign points one at a time, and check with previous points
    (write-line `(Creating smooth layout: mindist ,(sqrt mindist)))

    (let lp ((i 0) (c 0))
      (if (>= i n) 'done
	  (begin (vector-set! x i (random 1.0))
		 (vector-set! y i (random 1.0))
		 (if (check-point i)
		     (lp (+ i 1) 0)
		     ; try to place the point 50 times before giving up
		     ; c keeps track of the number of trys
		     (if (< c 50)
			 (begin (write-line (list "repeated node" i c "times"))
				(lp i (+ c 1)))
			 (lp (+ i 1) 0))))
	  ))

    ; remaining code is identical to make-sim/1
    (return-sim-with-layout x y z n r)))
|#

;; General Utility Functions
;;---------------------------

(define pi (* 4 (atan 1 1)))
; nbr = * pi r^2 total
(define (radius->nbrhood r total)     (* pi r r total))
(define (nbrhood->radius nbr total)   (sqrt (/ nbr (* pi total))))
(define (nbrhood->total nbrs r)	      (ceiling->exact (/ nbrs (* pi r r))))

; simulation.find-processor-nearest is a useful library routine
; which searches an existing simulation for the processor closest
; to the given coordinates. (dder)

(define (simulation.find-processor-nearest sim x y #!optional z)
  (define (sq x) (* x x))
  (let ((z (if (default-object? z) 0 z))
	(x-vector (simulation.x sim))
	(y-vector (simulation.y sim))
	(z-vector (simulation.z sim))
	(r2-min 1e99)   ; the smallest square of distance found (start > 3)
	(min-i #F))     ; index of the processor at r2-min (else #F)
    (do ((i 0 (+ i 1))) ((= i (simulation.n-processors sim)))
      (let ((r2 (+ (sq (- x (vector-ref x-vector i)))
		   (sq (- y (vector-ref y-vector i)))
		   (sq (- z (vector-ref z-vector i))))))
	(if (< r2 r2-min)
	    (begin
	      (set! min-i i)
	      (set! r2-min r2)))))
    min-i))

; simulation.find-processors-in-region is a useful library routine which 
; searches an existing simulation for all the processors which satisfy a 
; predicate (specified as a function of the processor coordinates).
; NB An alternative approach to specifying arbitrary regions is to use
; a drawing package to produce a PNM image for the HLSIM sensor facility.
; (dder)

(define (simulation.find-processors-in-region sim in-region?)
  (let ((x-vector (simulation.x sim))
	(y-vector (simulation.y sim))
	(z-vector (simulation.z sim))
	(n (simulation.n-processors sim)))
    (let loop ((i 0) (result '()))
	 (cond ((= i n) result)
	       ((in-region? (vector-ref x-vector i)
			    (vector-ref y-vector i)
			    (vector-ref z-vector i)) 
		(loop (1+ i) (cons i result)))
	       (else (loop (1+ i) result))))))

; The make-circle-predicate procedure returns a predicate for use in
; simulation.find-processors-in-region.  The predicate describes a
; circular region with origin at (x, y) and radius r (excluding points
; with radius equal to r).  For example, to find processors in a circle
; centered at (0.5, 0.5) with radius 0.1, use
; (simulation.find-processors-in-region sim (make-circle-predicate .5 .5 .1))

(define (make-circle-predicate x0 y0 r)
  (let ((sq (lambda (a) (* a a)))
	(r2 (* r r)))
       (lambda (x y z)
	       (< (+ (sq (- x x0)) (sq (- y y0))) r2))))

; Compute Statistics
; average and std dev in neighborhood size

(define (compute-nbrhood-stats sim)
  (let ((lst (vector->list (simulation.neighbours sim))))
    (compute-avg-stddev (map vector-length lst))))

(define (compute-avg-stddev lst)
  (let ((n (length lst)) (avg 0) (stddev 0))
    (set! avg (/ (fold-right + 0 lst) n))
    (let ((tmp (map (lambda (x) (square (- x avg))) lst)))
      (set! stddev (sqrt (/ (fold-right +  0 tmp) n))))
    (cons (exact->inexact avg) (exact->inexact stddev))))

; Turn off broadcast notices
; to undo reload comms*
(define (notify-broadcast sim my-id message) #F)
(define (notify-delivery sim my-id message) #F)

; To turn on compiler help
; (set! compiler:generate-type-checks? #T)
; (set! compiler:generate-range-checks? #T)