(require 'array)
(require 'subarray)
(require 'array-for-each)
(require 'eps-graph)
(require 'printf)
(require 'Peano-fill)
(require 'Hilbert-fill)

(define (euclidean-distance p1 p2)
  (define offs (map - p1 p2))
  (sqrt (apply + (map * offs offs))))

(define (profile-closeness proc rank samples start end)
  (define len (- end start))
  (let ((pra (make-array (vector) len))
	(results (make-array (A:floR64b) samples 2)))
    (do ((idx start (+ 1 idx))
	 (jdx 0 (+ 1 jdx)))
	((>= idx end))
      (array-set! pra (proc idx rank) jdx))
    (do ((off samples (+ -1 off)))
	((< off 1) results)
      (let ((acc 0))
	(array-for-each (lambda (p1 p2)
			  (define d12 (euclidean-distance p1 p2))
			  (set! acc (+ acc d12)))
			(subarray pra (list 0 (- len off 1)))
			(subarray pra (list off (+ -1 len))))
	(array-set! results off (+ -1 off) 0)
	(array-set! results (/ acc (- len off) off) (+ -1 off) 1)))))

(define (profile-log-closeness proc rank samples start end)
  (define len (- end start))
  (let ((pra (make-array (vector) len))
	(results (make-array (A:floR64b) samples 2)))
    (do ((idx start (+ 1 idx))
	 (jdx 0 (+ 1 jdx)))
	((>= idx end))
      (array-set! pra (proc idx rank) jdx))
    (do ((off samples (+ -1 off)))
	((< off 1) results)
      (let ((acc 0))
	(array-for-each (lambda (p1 p2)
			  (define d12 (euclidean-distance p1 p2))
			  (set! acc (+ acc d12)))
			(subarray pra (list 0 (- len off 1)))
			(subarray pra (list off (+ -1 len))))
	(array-set! results (real-log10 off) (+ -1 off) 0)
	(array-set! results (real-log10 (/ acc (- len off))) (+ -1 off) 1)))))

;;; Defaults from graph parameters
(define graph:dimensions '(550 350))
(define graph:font "Times")
(define graph:font-size 13)

(define (profile name title proc distance start end)
  (let ((xrange (list 1 distance)))
    (apply create-postscript-graph
	   (string-append name ".eps") graph:dimensions
	   (set-linewidth 0)
	   (set-font graph:font graph:font-size)
	   (set-margin-templates "++1.0" "++1.0")
	   (whole-page)
	   (setup-plot xrange '(0 1))
	   (title-top (sprintf #f "%s; %d points"
			       (or title name) (- end start))
		      "Rank 2, 3, 4, 6, and 9")
	   (outline-rect plotrect)
	   (rule-vertical leftedge "Euclidean-distance / scalar-distance" 10)
	   (rule-horizontal bottomedge "scalar distance" 5)
	   (grid-horizontals)
	   (grid-verticals)
	   (clip-to-rect plotrect)
	   ;;(set-color 'heath)
	   (map (lambda (dim)
		  (plot-column (profile-closeness proc dim distance start end)
			       0 1 'line))
		'(2 3 4 6 9)))
    ;;(eps:viewer name)
    (eps:convert name)
    ))

(define (profile-log name title proc distance start end)
  (let ((xrange (list 0 (real-log10 distance))))
    (apply create-postscript-graph
	   (string-append name ".eps") graph:dimensions
	   (set-linewidth 0)
	   (set-font graph:font graph:font-size)
	   (set-margin-templates "++1.0" "++1.0")
	   (whole-page)
	   (setup-plot xrange '(0 1.25))
	   (title-top (sprintf #f "%s; %d points"
			       (or title name) (- end start))
		      "Rank 2, 3, 4, 6, and 9")
	   (outline-rect plotrect)
	   (rule-vertical leftedge "log_10(Euclidean-distance)" 10)
	   (rule-horizontal bottomedge "log_10(scalar-distance)" 5)
	   (grid-horizontals)
	   (grid-verticals)
	   (clip-to-rect plotrect)
	   ;;(set-color 'heath)
	   (map (lambda (dim)
		  (plot-column (profile-log-closeness proc dim distance start end)
			       0 1 'line))
		'(2 3 4 6 9)))
    ;;(eps:viewer name)
    (eps:convert name)
    ))

(define (eps:viewer path)
  (system (string-append "gv '" path ".eps'")))

(define (eps:convert path)
  (and (system (string-append "convert '" path ".eps' '" path ".png'"))
       (delete-file (string-append path ".eps"))))

(define (bitwise-delaminate count k)
  (define nbs (* count (+ 1 (quotient (integer-length k) count))))
  (do ((kdx (- nbs count) (- kdx count))
       (lst (vector->list (make-vector count 0))
	    (map (lambda (k bool) (+ (if bool 1 0) (arithmetic-shift k 1)))
		 lst
		 (integer->list (arithmetic-shift k (- kdx)) count))))
      ((negative? kdx) lst)))

(profile "Z-curve-r" #f
	 (lambda (idx rank) (bitwise-delaminate rank idx))
	 15 0 (expt 2 12))
(profile "Peano-r" "Peano Space-Filling Curve"
	 natural->peano-coordinates
	 15 0 (expt 3 6))
(profile "Hilbert-r" "Hilbert Space-Filling Curve"
	 integer->hilbert-coordinates
	 15 0 (expt 2 12))

(profile-log "Z-curve-l" #f
	     (lambda (idx rank) (bitwise-delaminate rank idx))
	     100 0 (expt 2 12))
(profile-log "Peano-l" "Peano Space-Filling Curve"
	     natural->peano-coordinates
	     100 0 (expt 3 6))
(profile-log "Hilbert-l" "Hilbert Space-Filling Curve"
	     integer->hilbert-coordinates
	     100 0 (expt 2 12))