Implements the up/down tuple structured objects used by scmutils/SICM, for Racket

structures.rkt

#lang racket

(require
  srfi/1
  math
  racket/struct
  ;; syntax stuff
  (for-syntax
   syntax/parse
   syntax/parse/lib/function-header))

;; for function name experimantation
(define-syntax (named-lambda stx)
  (syntax-parse stx
    [(_named-lambda (f:id . args:formals) b ...)
     (syntax/loc stx
       (procedure-rename (λ args b ...) 'f))]))

;; ================================================

;; All "structures" are implemented as applicable structs
;; ======================================================
;; viz. matrices, up/down tuples and power series...later
;; (however see note about native Racket vectors, below)

;; Apply... "backwards"
(define (apply-struct S val)
  (define ((reverse-apply val) func) (func val))
  (tup-val-map (reverse-apply val) S))
      
;; Tuple Structures
;; ================

; the up/down tuple structure

(struct tuple (type value)
  #:transparent
  #:methods gen:custom-write
  [(define write-proc
     (make-constructor-style-printer
      (λ (t) (tuple-type t))
      (λ (t) (tuple->list t))))]
  #:property prop:custom-print-quotable 'never
  #:property prop:procedure
  (λ (f x)
    (match-define (tuple up/down contents) f)
    (tuple up/down (apply-struct contents x))))

;; Template functions
;; ------------------
;; Allow us to parameterize the containers which hold the value of the tuple
;; For instance, scmutils used 'vector', but other types can also be used,
;; for instance the Math Library Array.
;; GJS himself states that their decision to use 'vector' might change in the
;; future so hopefully this won't cause any problems.
;;
;; 'vector' is applicable in scmutils, but only by being promoted by default
;; to a structure, i.e. ((vector sin cos tan) 1) => (up .84147 .54030 1.5575)
;; I do not want to implement 'vector' as a struct (OBVIOUSLY), 
;; so if the above case ever arises a simple transformation 
;; ((vector ...) arg) => ((up ...) arg)  should handle it.

#|
; template for using Array as container
; (define aa (structure:expt (up 1 (down 2 4 6) 3 5) 9)) takes >300"
(define tup-val-map array-map)
(define tup-val->list array->list)
(define tup-val<-list list->array)
(define tup-val->vector array->vector)
(define tup-val<-vector vector->array)
(define tup-val-size array-size)
(define tup-val-ref array-ref)
(define tup-val-index vector)
(define tup-val-build build-array)
(define tup-val-proc (λ(proc)(compose proc first vector->list)))
|#


; template for using vector as container
; (define aa (structure:expt (up 1 (down 2 4 6) 3 5) 9)) takes 30"
; (s:dimension aa) => 10077696  ; acceptable
(define tup-val-map vector-map)
(define tup-val->list vector->list)
(define tup-val<-list list->vector)
(define tup-val->vector identity)
(define tup-val<-vector identity)
(define tup-val-size vector-length)
(define tup-val-ref vector-ref)
(define tup-val-index identity)
(define tup-val-build build-vector)
(define tup-val-proc identity)

; The two instances

(define (up . args)
  (list->tuple 'up args))
       
(define (down . args)
  (list->tuple 'down args))


;; Functions  
;; ---------
;; (tuple->list t)
;; (list->tuple type li)
;; (structure? s)
;; (up? s)
;; (down s)
;; (s:structure up/down array)
;; (vector->up v)
;; (vector->down v)
;; (s:->vector s)
;; (up->vector t)
;; (down->vector t)
;; (s:same s)
;; (s:opposite s)
;; (s:length s)
;; (s:dimension s)
;; (s:ref s i)
;; (s:generate n type proc)
;; (s:for-all proc s)
;; (s:fringe s1)
;; (s:for-each proc s)
;; (structure:expt s n)
;; (s:multiply s1 s2)
;; (s:* s1 s2)
;; (s:compatible-for-contraction? v1 v2)
;; (s:compatible-elements? v1 v2)
;; (scalar*structure s v)
;; (structure*scalar v s)
;; (structure/scalar v s)
;; (s:square)
;; (s:dot-product)
;; (s:outer-product s1 s2)
;; (s:map)
;; (s:map/l)
;; (s:map/r)
;; (s:map/r/l)
;; (s:elementwise proc)
;; (structure:elementwise proc)


; lists <-> tuple conversions

(define (tuple->list t)
  (tup-val->list (tuple-value t)))

(define (list->tuple type li)
  (tuple type (tup-val<-list li)))

; structure operations

(define structure? tuple?)

(define (up? t)
  (and (tuple? t) (eq? (tuple-type t) 'up)))

(define (down? t)
  (and (tuple? t) (eq? (tuple-type t) 'down)))

(define (s:structure up/down contents)  ;; that's a math library array
  (case up/down
    ((up contravariant vector)
     (tuple 'up contents))
    ((down covariant covector)
     (tuple 'down contents))
    (else
     (error "Bad up/down spec -- s:structure" up/down contents))))  

(define (vector->up v)
  (tuple 'up  (tup-val<-vector v)))

(define (vector->down v)
  (tuple 'down (tup-val<-vector v)))

(define (s:->vector t)
  (if (or (up? t) (down? t))
      (tup-val->vector (tuple-value t))
      (error "Bad structure -- s:->vector" t)))

(define (up->vector s)
  (if (up? s) (s:->vector s)
      (error "not up tuple -- up->vector")))
      
(define (down->vector s)
  (if (down? s) (s:->vector s)
      (error "not down tuple -- down->vector")))

(define (s:same t)
  (cond ((up?  t) 'up)
	((down? t) 'down)
	(else (error "Bad structure -- s:same" t))))

(define (s:opposite t)
  (cond ((up? t) 'down)
	((down? t) 'up)
	(else (error "Bad structure -- s:opposite" t))))

(define (s:length t)
  (if (structure? t)
      (tup-val-size (tuple-value t))
      1))

(define (s:dimension s)
  (if (structure? s)
      (reduce +
	      0
	      (map s:dimension
		   (tup-val->list (tuple-value s))))
      1))

(define (s:ref t . ilist)
  (define (s-ref str args)
    (let* ([num (car args)]
           [tail (cdr args)]
           [elem (tup-val-ref (tuple-value str) num)])
    (if (empty? tail) elem (s-ref elem tail))))
  (cond
    [(not (structure? t)) (if (equal? ilist '(0)) t
                              (error "Bad structure -- s:ref" t ilist))]
    [(empty? ilist) t]
    [else (s-ref t ilist)]))

(define (s:generate n up/down proc)  ; 0-based
  (tuple up/down
         (tup-val-build (tup-val-index n) (tup-val-proc proc))))

(define (s:forall proc s)
  (let ((n (s:length s)))
    (let loop ([i 1] [ans (proc (s:ref s 0))])
      (cond ((eq? i n) ans)
	    ((not ans) ans)
	    (else
	     (loop (add1 i) (proc (s:ref s i))))))))  

;;; S:FRINGE recursively traverses a structure, making up a list of
;;; the terminal elements.

(define (s:fringe s)
  (define (walk s ans)
    (if (structure? s)
	(let ([n (s:length s)])
	  (let loop ([i 0] [ans ans])
	    (if (eq? i n)
		ans
		(loop (add1 i)
                      (walk (s:ref s i) ans)))))
        (cons s ans)))
  (walk s '()))


(define (s:foreach proc s)  ;; for side-effects and mutable
  (define (walk s)
    (if (structure? s)
	(let ((n (s:length s)))
	  (let loop ((i 0))
	    (if (eq? i n)
		'done
		(begin (walk (s:ref s i))
		       (loop (add1 i))))))
	(proc s)))
  (walk s))


;; arithmetic on structures

(define (structure:expt t n)
  (cond [(zero? n) (error "Undefined power of zero for structures")]
        [(eq? n 2) (s:square t)] ; special case
        [else
         (define (s:expt t n)
           (cond [(eq? n 1) t]
                 [(> n 1) (s:* t (s:expt t (sub1 n)))]  ;; gen *
                 [else (error "Cannot: " `(expt ,t ,n))]))
         (s:expt t n)]))

(define (s:multiply s1 s2)  ;; include special case where identical, and square
  (cond [(s:compatible-for-contraction? s1 s2)
	 (s:dot-product s1 s2)]
	[(or *allowing-incompatible-multiplication*
	     (and (or (and (down? s1) (down? s2))
		      (and (up? s1) (up? s2)))
		  (s:forall (lambda (c)
			      (s:compatible-for-contraction? s1 c))
			    s2)))
         (s:generate (s:length s2) (s:same s2)
		     (lambda (i)
                       ;; my temp soln for recursive non-compatibles
                       (let ([factor (s:ref s2 i)])
                         (if (structure? factor)
                             (s:multiply s1 factor)
                             (structure*scalar s1 factor)))))]
	[else
	 (error "Incompatible multiplication" s1 s2)]))

(define s:* s:multiply)

(define *allowing-incompatible-multiplication* #t)

(define (s:compatible-for-contraction? v1 v2)
  (or (and (down? v1) (up? v2)
	   (s:compatible-elements? v1 v2))
      (and (up? v1) (down? v2)
	   (s:compatible-elements? v1 v2))))


(define (s:compatible-elements? v1 v2)
  (let ((n (s:length v1)))
    (and (eq? n (s:length v2))
	 (let loop ([i 0])
	   (cond ((eq? i n) #t)
		 ((or (not (structure? (s:ref v1 i)))
		      (not (structure? (s:ref v2 i))))
		  (loop (add1 i)))
		 ((s:compatible-for-contraction? (s:ref v1 i)
						 (s:ref v2 i))
		  (loop (add1 i)))
		 (else #f))))))

(define (scalar*structure x t) 
  (s:map/r (curry * x) t))

(define (structure*scalar t x)
  (scalar*structure x t)) 

(define (structure/scalar t x)
  (s:map/r (λ(y) (/ y x) t)))



#|
;; ***AG** (the code below is from scmutils)
;; might work better, since recursive -- check this
;;; Is this redundant with (s:dot-product v v)?
(define (s:square v)
  (let ((vv (s:->vector v)))
    (let ((n (vector-length vv)))
      (if (fix:= n 0)
	  :zero
	  (let lp ((i 1) (sum (g:square (vector-ref vv 0))))
	    (if (fix:= i n)
		sum
		(lp (fix:+ i 1)
		    (g:+ sum (g:square (vector-ref vv i))))))))))
|#

;; ***AG*** the following two seem simplistic!!
;; check/analyse this...

;; make fully-recursive for future, allowing polymorphic *, operators, etc
(define (s:square s)
  (let ([flat-elems (s:fringe s)])
    (apply + (map * flat-elems flat-elems))))

;; make fully-recursive for future, allowing polymorphic *, operators, etc
(define (s:dot-product t1 t2);   ***AG*** poor checks.  Improve!
  (cond  [(eq? (s:same t1) (s:same t2))
          (error "Incompatible structures -- s:dot-product" t1 t2)]
         [(= (s:dimension t1) (s:dimension t2))
          (apply + (map * (s:fringe t1) (s:fringe t2)))]
         [else
          (error "Incompatible dimensions -- s:dot-product" t1 t2)]))  ;; generalise +,*


;;; Given two structures their outer product makes a structure
(define (s:outer-product t2 t1)
  (s:map/r (lambda (s1)
	     (s:map/r (lambda (s2) (* s1 s2)) t2) ;; this * must become generic
             t2)
           t1))


;;; --- MAPPERS ---
;;; The following mappers only make sense if, when there is more than
;;; one structure they are all isomorphic.

(define (s:map proc . structures)
  (s:map/l proc structures))

(define (s:map/l proc structures)
  (if (structure? (car structures))
      (s:generate (s:length (car structures))
		  (s:same (car structures))
		  (λ (i)
		    (apply proc
			   (map (λ (s) (s:ref s i)) structures))))
      (apply proc structures)))

(define (s:map/r proc . structures)
  (s:map/r/l proc structures))

(define (s:map/r/l proc structures)
  (s:map/l (lambda elements
	     (if (structure? (car elements))
		 (s:map/r/l proc elements)
		 (apply proc elements)))
	   structures))

(define ((s:elementwise proc) . structures)
  (s:map/l proc structures))

(define structure:elementwise s:elementwise)


;; ToDos
;; =====
;; language syntax for polymorphism

#|
;; Testing....
;; -----------

;; apply structure
((up sin cos tan) 1)

;; simple struture map - unary function
(let* ([st1 (up 2 3 4)]
       [st2 (s:map sqrt st1)])
    (list st1 st2))

;; recursive struture map - n-ary function
(let ([s1 (up 1 (down 2 4 8) 3)]
        [s2 (up 7 (down 4 5 6) 9)]
        [s3 (up 7 (down 1 8 9) 5)])
    (s:map/r + s1 s2 s3))

;; contraction of compatible structures as per manual
(s:* (up (up 2 3) (down 5 7 11)) (down (down 13 17) (up 19 23 29)))

;; demo of non-commutative outer product as per manual
(s:* (up 2 3) (up 5 7 11))
(s:* (up 5 7 11) (up 2 3))

;; tuple exponentation (special case for 2)
(structure:expt (up 1 2 3) 2)
(structure:expt (up 1 2 3) 3)
|#