;;;; "rmdsff.scm" Space-filling functions and their inverses.
;;; Copyright (C) 2013, 2014 Aubrey Jaffer
;
;Permission to copy this software, to modify it, to redistribute it,
;to distribute modified versions, 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.  I have made no warranty or representation that the operation of
;this software will be error-free, and I am under no obligation to
;provide any services, by way of maintenance, update, or otherwise.
;
;3.  In conjunction with products arising from the use of this
;material, there shall be no use of my name in any advertising,
;promotional, or sales literature without prior written consent in
;each case.

(require 'array)

;;@code{(require 'space-filling)}
;;@ftindex space-filling

;;@ The algorithms and cell properties are described in
;;@url{http://people.csail.mit.edu/jaffer/Geometry/RMDSFF.pdf}

;;; A cell is an object encapsulating the information about a
;;; Hamiltonian path on a rectangular grid of side^rank nodes.
;;; Here are the accessors for a cell:
(define cell-type caar)
(define cell-side cadar)
(define cell-rank caddar)
(define (cell-index cell crds) (apply array-ref (cadadr cell) crds))
(define (cell-coords cell t) (vector-ref (caadr cell) t))
(define (cell-entry cell t) (vector-ref (caddr cell) t))
(define (cell-exit cell t) (vector-ref (cadddr cell) t))
(define (cell-rotation cell t) (vector-ref (cadddr (cdr cell)) t))

;;@args type rank side precession
;;@args type rank side
;;@args type rank
;;
;;@1 must be the symbol @code{diagonal}, @code{adjacent}, or
;;@code{centered}.  @2 must be an integer larger than 1.  @3, if
;;present, must be an even integer larger than 1 if @1 is
;;@code{adjacent} or an odd integer larger than 2 otherwise; @3
;;defaults to the smallest value.  @4, if present, must be an integer
;;between 0 and @3^@2-1; it is relevant only when @1 is
;;@code{diagonal} or @code{centered}.
;;
;;@args Hamiltonian-path-vector precession
;;@args Hamiltonian-path-vector
;;
;;@1 must be a vector of @var{side}^@var{rank} lists of @var{rank} of
;;integers encoding the coordinate positions of a Hamiltonian path on
;;the @var{rank}-dimensional grid of points starting and ending on
;;corners of the grid.  The starting corner must be the origin
;;(all-zero coordinates). If the side-length is even, then the ending
;;corner must be non-zero in only one coordinate; otherwise, the
;;ending corner must be the furthest diagonally opposite corner from
;;the origin.
;;
;;@code{make-cell} returns a data object suitable for passing as the
;;first argument to @code{integer->coordinates} or
;;@code{coordinates->integer}.
(define (make-cell arg1 . args)
  (define (make-serpentine-path rank s)
    (let loop ((path '(())) (rnk (+ -1 rank)))
      (if (negative? rnk) path
	  (loop (let iloop ((seq '()) (sc (+ -1 s)))
		  (if (negative? sc)
		      seq
		      (iloop (append (map (lambda (coords) (cons sc coords))
					  (if (odd? sc) (reverse path) path))
				     seq)
			     (+ -1 sc))))
		(+ -1 rnk)))))
  (if (list? arg1) (set! arg1 (list->vector arg1)))
  (cond
   ((> (length args) 3)
    (slib:error 'make-cell 'extra 'arguments 'not 'handled args))
   ((vector? arg1)
    (let ((path arg1)
	  (precession (and (not (null? args)) (car args))))
      (define frst (vector-ref path 0))
      (define len (vector-length path))
      (define s-1 (apply max (apply append (vector->list path))))
      (let* ((len-1 (+ -1 len))
	     (last (vector-ref path len-1))
	     (d (length frst)))
	;; returns index of first non-zero in LST
	(define (first-non-zero lst)
	  (define pos (lambda (n lst)
			(cond ((zero? (car lst)) (pos (+ 1 n) (cdr lst)))
			      (else n))))
	  (pos 0 lst))
	;; returns the traversal direction of the sub-path.
	(define (U_e N t)
	  (if (= t len-1)
	      (map - (vector-ref path t) (vector-ref path (- t 1)))
	      (let ((dH_i+1 (map - (vector-ref path (+ t 1)) (vector-ref path t))))
		(define dotpr (apply + (map * N dH_i+1)))
		(define csum (apply + dH_i+1))
		(if (or (and (zero? dotpr) (= 1 csum)) (= dotpr csum -1))
		    dH_i+1
		    (map - (vector-ref path t) (vector-ref path (- t 1)))
		    ))))
	(define (flip-direction dir cdir)
	  (map (lambda (px c) (modulo (+ px c) 2)) dir cdir))
	(define (path-diag? path)
	  (define prev frst)
	  (for-each (lambda (lst)
		      (if (not (= d (length lst)))
			  (slib:error 'non-uniform 'ranks frst lst)))
		    (vector->list path))
	  (for-each (lambda (cs)
		      (if (not (= 1 (apply + (map abs (map - prev cs)))))
			  (slib:error 'bad 'step prev cs))
		      (set! prev cs))
		    (cdr (vector->list path)))
	  (cond ((not (zero? (apply + frst))) (slib:error 'strange 'start frst))
		((not (= d (length last))) (slib:error 'non-uniform 'lengths path))
		((apply = s-1 last) #t)
		((and (= s-1 (apply + last))) #f)
		(else (slib:error 'strange 'net-travel frst last))))
	(define diag? (path-diag? path))
	(define entries (make-vector len (vector-ref path 0)))
	(define exits (make-vector
		       len (if diag?
			       (map (lambda (c) (quotient c s-1)) last)
			       (vector-ref path 1))))
	(define rotations (make-vector len 0))
	(define ipath (apply make-array
			     (A:fixZ32b -1)
			     (vector->list (make-vector d (+ 1 s-1)))))
	(define ord 0)
	(for-each (lambda (coords)
		    (apply array-set! ipath ord coords)
		    (set! ord (+ 1 ord)))
		  (vector->list path))
	(let lp ((t 1)
		 (prev-X (if diag?
			     (map (lambda (c) 1) (vector-ref entries 0))
			     (vector-ref path 1))))
	  (cond ((> t len-1)
		 (if (not (equal? (vector-ref exits len-1)
				  (map (lambda (c) (quotient c s-1)) last)))
		     (slib:warn (list (if diag? 'diagonal 'adjacent)
				      (+ 1 s-1) d precession)
				'bad 'last 'exit
				(vector-ref exits len-1) 'should 'be
				(map (lambda (c) (quotient c s-1)) last)))
		 (let ((ord 0)
		       (h (first-non-zero last)))
		   (for-each
		    (lambda (coords)
		      (vector-set! rotations
				   ord
				   (modulo
				    (cond ((not diag?)
					   (- h (first-non-zero
						 (map -
						      (vector-ref exits ord)
						      (vector-ref entries ord)))))
					  (precession (+ precession ord))
					  (else 0))
				    d))
		      (set! ord (+ 1 ord)))
		    (vector->list path)))
		 (list (list (if diag? 'diagonal 'adjacent)
			     (+ 1 s-1)
			     d
			     precession)
		       (list path ipath)
		       entries
		       exits
		       rotations
		       ))
		(else
		 (let ((N (flip-direction
			   prev-X
			   (map - (vector-ref path t)
				(vector-ref path (- t 1))))))
		   (define X (if diag?
				 (map (lambda (tn) (- 1 tn)) N)
				 (flip-direction N (U_e N t))))
		   (vector-set! entries t N)
		   (vector-set! exits t X)
		   (lp (+ 1 t) X))))))))
   ((< (car args) 2)
    (slib:error 'make-cell 'rank 'too 'small (car args)))
   (else
    (case arg1
      ((center centered)
       (let ((cell (make-cell (make-serpentine-path
			       (car args)
			       (if (null? (cdr args)) 3 (cadr args)))
			      (if (= 3 (length args)) (caddr args) #f))))
	 (if (not (eq? 'diagonal (cell-type cell)))
	     (slib:error 'make-cell 'centered 'must 'be 'diagonal (car cell)))
	 (set-car! (car cell) 'centered)
	 cell))
      ((diagonal opposite)
       (make-cell (make-serpentine-path
		   (car args)
		   (if (null? (cdr args)) 3 (cadr args)))
		  (if (= 3 (length args)) (caddr args) #f)))
      ((adjacent)
       (make-cell (make-serpentine-path
		   (car args)
		   (if (null? (cdr args)) 2 (cadr args)))
		  (if (= 3 (length args)) (caddr args) #f)))
      (else
       (slib:error 'make-cell 'unknown 'cell 'type arg1))))))

;;@ Hilbert, Peano, and centered Peano cells are generated
;;respectively by:
;;@example
;;(make-cell 'adjacent @var{rank} 2)   ; Hilbert
;;(make-cell 'diagonal @var{rank} 3)   ; Peano
;;(make-cell 'centered @var{rank} 3)   ; centered Peano
;;@end example

;; Positive k rotates left
(define (rotate-list lst k)
  (define len (length lst))
  (cond ((<= len 1) lst)
	(else
	 (set! k (modulo k len))
	 (if (zero? k)
	     lst
	     (let ((ans (list (car lst))))
	       (do ((ls (cdr lst) (cdr ls))
		    (tail ans (cdr tail))
		    (k (+ -1 k) (+ -1 k)))
		   ((<= k 0)
		    (append ls ans))
		 (set-cdr! tail (list (car ls)))))))))

;;@ In the conversion procedures, if the cell is @code{diagonal} or
;;@code{adjacent}, then the coordinates and scalar must be nonnegative
;;integers.  If @code{centered}, then the integers can be negative.

;;@body
;;@0 converts the integer @2 to a list of coordinates according to @1.
(define (integer->coordinates cell u)
  (define umag (case (cell-type cell)
		 ((centered) (* (abs u) 2))
		 (else u)))
  (define d (cell-rank cell))
  (define s (cell-side cell))
  (let* ((s^d (expt s d))
	 (s^d^2 (expt s^d d)))
    (define (align V t sde)
      (map (lambda (Vj Nj) (if (zero? Nj) Vj (- sde Vj)))
	   (rotate-list V (cell-rotation cell t))
	   (cell-entry cell t)))
    (define (rec u m w)
      (if (positive? w)
	  (let ((t (quotient u m)))
	    (map +
		 (map (lambda (y) (* y w)) (cell-coords cell t))
		 (align (rec (modulo u m) (quotient m s^d) (quotient w s))
			t
			(- w 1))))
	  (cell-coords cell 0)))
    (do ((uscl 1 (* uscl s^d^2))
	 (cscl 1 (* cscl s^d)))
	((> uscl umag)
	 (case (cell-type cell)
	   ((centered)
	    (let ((cscl/2 (quotient cscl 2)))
	      (map (lambda (c) (- c cscl/2))
		   (rec (+ u (quotient uscl 2))
			(quotient uscl s^d)
			(quotient cscl s)))))
	   (else (rec u (quotient uscl s^d) (quotient cscl s))))))))

;;@body
;;@0 converts the list of coordinates @2 to an integer according to @1.
(define (coordinates->integer cell V)
  (define maxc (case (cell-type cell)
		 ((centered) (* 2 (apply max (map abs V))))
		 (else (apply max V))))
  (define d (cell-rank cell))
  (define s (cell-side cell))
  (let* ((s^d (expt s d))
	 (s^d^2 (expt s^d d)))
    (define (align^-1 V t sde)
      (rotate-list (map (lambda (Vj Nj) (if (zero? Nj) Vj (- sde Vj)))
			V
			(cell-entry cell t))
		   (- (cell-rotation cell t))))
    (define (rec u V w)
      (if (positive? w)
	  (let ((dig (cell-index cell (map (lambda (c) (quotient c w)) V))))
	    (rec (+ dig (* s^d u))
		 (align^-1 (map (lambda (cx) (modulo cx w)) V)
			   dig
			   (- w 1))
		 (quotient w s)))
	  u))
    (do ((uscl 1 (* uscl s^d^2))
	 (cscl 1 (* cscl s^d)))
	((> cscl maxc)
	 (case (cell-type cell)
	   ((centered)
	    (let ((cscl/2 (quotient cscl 2)))
	      (- (rec 0
		      (map (lambda (c) (+ c cscl/2)) V)
		      (quotient cscl s))
		 (quotient uscl 2))))
	   (else (rec 0 V (quotient cscl s))))))))

;;@var{coordinates->integer} and @var{integer->coordinates} are
;;inverse functions when passed the same @var{cell} argument.