;;;; math-utils.vector.scm  -*- Scheme -*-
;;;; Kon Lovett, Jul '23
;;;; Kon Lovett, Feb '24

;;;; Issues
;;
;; - Functor (VECTOR FIELD) - f32vector & fL | ...

(module (math-utils vector)

(;export
  absolute-magnitude
  cosine-similarity
  vector-apply
  vector-mul vector-div
  vector-sum vector-dif
  vector-min vector-max vector-lcm vector-gcd
  vector-compare
  dot-product
  cross-product)

(import scheme)
(import (chicken base))
(import (chicken fixnum))
(import (chicken type))

;(define-type real (or integer ratnum float))

(define-type boolean-vector (vector-of boolean))
(define-type number-vector (vector-of number))

(define-type nary (* * #!rest * --> *))
#; ;UNUSED
(define-type nary-numeric (number number #!rest number --> *))

(: absolute-magnitude (number-vector --> number))

(: cosine-similarity (number-vector number-vector --> number))

(: vector-apply (nary #!rest vector --> vector))

(: vector-mul (number-vector number-vector #!rest number-vector --> number-vector))
(: vector-div (number-vector number-vector #!rest number-vector --> number-vector))
(: vector-sum (number-vector number-vector #!rest number-vector --> number-vector))
(: vector-dif (number-vector number-vector #!rest number-vector --> number-vector))
(: vector-min (number-vector number-vector #!rest number-vector --> number-vector))
(: vector-max (number-vector number-vector #!rest number-vector --> number-vector))
(: vector-lcm (number-vector number-vector #!rest number-vector --> number-vector))
(: vector-gcd (number-vector number-vector #!rest number-vector --> number-vector))

(: vector-compare (number-vector number-vector #!rest number-vector --> number))

(: dot-product (number-vector number-vector --> number))
(: cross-product (number-vector number-vector --> (or number number-vector)))

;;

(define-inline (sqrt-exact x) (inexact->exact (sqrt x)))
(define-inline (sqr a) (* a a))

(define (absolute-magnitude v)
  (let ((len (vector-length v)))
    (let loop ((i 0) (acc 0))
      (if (fx= i len) (sqrt-exact acc)
          (loop (fx+ i 1) (+ (sqr (vector-ref v i)) acc))) ) ) )

;https://en.wikipedia.org/wiki/Cosine_similarity
(define (cosine-similarity a b)
  ;unroll
  (let ((len (vector-length a)))
    (let loop ((i 0) (dot 0) (maga 0) (magb 0))
      (if (fx= i len) (/ dot (sqrt-exact (* maga magb)))
          (let ((ai (vector-ref a i)) (bi (vector-ref b i)))
            (loop (fx+ i 1) (+ (* ai bi) dot) (+ (sqr ai) maga) (+ (sqr bi) magb)) ) ) ) ) )

; assumes all same length
(define (vector-apply f . vs)
  (let* ((a (car vs))
         (n (vector-length a))
         (r (make-vector n)) )
    (do ((i 0 (fx+ i 1)))
        ((fx= i n) r)
      (set! (vector-ref r i)
            (apply f (map (cut vector-ref <> i) vs))) ) ) )

(define (vector-mul . vs) (apply vector-apply * vs))
(define (vector-div . vs) (apply vector-apply / vs))
(define (vector-sum . vs) (apply vector-apply + vs))
(define (vector-dif . vs) (apply vector-apply - vs))
(define (vector-min . vs) (apply vector-apply min vs))
(define (vector-max . vs) (apply vector-apply max vs))
(define (vector-lcm . vs) (apply vector-apply lcm vs))
(define (vector-gcd . vs) (apply vector-apply gcd vs))

(define vector-compare
  (case-lambda
    ((a b)
      (let* ((alen (vector-length a))
             (blen (vector-length b))
             (dlen (fx- alen blen)) )
        (if (not (fx= 0 dlen))
            dlen
            (let loop ((i 0))
              (if (= i alen)
                  0
                  (let ((d (- (vector-ref a i) (vector-ref b i))))
                    (if (= 0 d)
                        (loop (fx+ i 1))
                        d ) ) ) ) ) ) )
    (vs
      (let loop ((vs vs))
        (if (or (null? vs) (null? (cdr vs)))
            0
            (let ((d (vector-compare (car vs) (cadr vs))))
              (if (= 0 d)
                  (loop (cdr vs))
                  d ) ) ) ) ) ) )

#;
(define (vector-compare . vs)
  (let ((vs vs))
    (if (or (null? vs) (null? (cdr vs)))
        0
        (let ((d (vector-compare (car vs) (cadr vs))))
          (if (= 0 d)
              (loop (cdr vs))
              d ) ) ) ) )

(define (dot-product a b)
  (let ((len (vector-length a)))
    (let loop ((i 0) (acc 0))
      (if (fx= i len) acc
          (loop (fx+ i 1) (+ (* (vector-ref a i) (vector-ref b i)) acc))) ) ) )

;0 1 2 -> X, 3 -> #(X Y Z), 4 -> #(X Y Z W), 8 -> #(X Y Z W I J K L), * -> Error
(define (cross-product a b)
  (let ((alen (vector-length a)) (blen (vector-length b)))
    (if (not (fx= alen blen))
        (error 'cross-product "undefined for varying length" alen blen)
        (case alen
          ((0)
           0 )
          ((1)
            (* (vector-ref a 0) (vector-ref b 0)) )
          ((2)
            (- (* (vector-ref a 0) (vector-ref b 1))
               (* (vector-ref a 1) (vector-ref b 0))) )
          ((3)
            (vector
              ;0
              (- (* (vector-ref a 1) (vector-ref b 2))
                 (* (vector-ref a 2) (vector-ref b 1)))
              ;1
              (- (* (vector-ref a 2) (vector-ref b 0))
                 (* (vector-ref a 0) (vector-ref b 2)))
              ;2
              (- (* (vector-ref a 0) (vector-ref b 1))
                 (* (vector-ref a 1) (vector-ref b 0)))) )
          ((4)
            ;{ {a[0] * b[0]} - {a[1] * b[1]} - {a[2] * b[2]} - {a[3] * b[3]} }
            ;{ {a[0] * b[1]} + {a[1] * b[0]} + {a[2] * b[3]} - {a[3] * b[2]} }
            ;{ {a[0] * b[2]} - {a[1] * b[3]} + {a[2] * b[0]} + {a[3] * b[1]} }
            ;{ {a[0] * b[3]} + {a[1] * b[2]} - {a[2] * b[1]} + {a[3] * b[0]} }
            (vector
              ;0
              (- (* (vector-ref a 0) (vector-ref b 0))
                 (* (vector-ref a 1) (vector-ref b 1))
                 (* (vector-ref a 2) (vector-ref b 2))
                 (* (vector-ref a 3) (vector-ref b 3)))
              ;1
              (- (+ (+ (* (vector-ref a 0) (vector-ref b 1))
                       (* (vector-ref a 1) (vector-ref b 0)))
                    (* (vector-ref a 2) (vector-ref b 3)))
                 (* (vector-ref a 3) (vector-ref b 2)))
              ;2
              (+ (+ (- (* (vector-ref a 0) (vector-ref b 2))
                       (* (vector-ref a 1) (vector-ref b 3)))
                    (* (vector-ref a 2) (vector-ref b 0)))
                 (* (vector-ref a 3) (vector-ref b 1)))
              ;3
              (+ (- (+ (* (vector-ref a 0) (vector-ref b 3))
                       (* (vector-ref a 1) (vector-ref b 2)))
                    (* (vector-ref a 2) (vector-ref b 1)))
                 (* (vector-ref a 3) (vector-ref b 0)))) )
          ((8)
            ;{ {a[0] * b[0]} - {a[1] * b[1]} - {a[2] * b[2]} - {a[3] * b[3]} - {b[4] * a[4]} - {b[5] * a[5]} - {b[6] * a[6]} - {b[7] * a[7]} }
            ;{ {a[0] * b[1]} + {a[1] * b[0]} + {a[2] * b[3]} - {a[3] * b[2]} - {b[4] * a[5]} + {b[5] * a[4]} + {b[6] * a[7]} - {b[7] * a[6]} }
            ;{ {a[0] * b[2]} - {a[1] * b[3]} + {a[2] * b[0]} + {a[3] * b[1]} - {b[4] * a[6]} - {b[5] * a[7]} + {b[6] * a[4]} + {b[7] * a[5]} }
            ;{ {a[0] * b[3]} + {a[1] * b[2]} - {a[2] * b[1]} + {a[3] * b[0]} - {b[4] * a[7]} + {b[5] * a[6]} - {b[6] * a[5]} + {b[7] * a[4]} }
            ;{ {b[4] * a[0]} - {b[5] * a[1]} - {b[6] * a[2]} - {b[7] * a[3]} + {a[4] * b[0]} + {a[5] * b[1]} + {a[6] * b[2]} + {a[7] * b[3]} }
            ;{ {b[4] * a[1]} + {b[5] * a[0]} + {b[6] * a[3]} - {b[7] * a[2]} - {a[4] * b[1]} + {a[5] * b[0]} - {a[6] * b[3]} + {a[7] * b[2]} }
            ;{ {b[4] * a[2]} - {b[5] * a[3]} + {b[6] * a[0]} + {b[7] * a[1]} - {a[4] * b[2]} + {a[5] * b[3]} + {a[6] * b[0]} - {a[7] * b[1]} }
            ;{ {b[4] * a[3]} + {b[5] * a[2]} - {b[6] * a[1]} + {b[7] * a[0]} - {a[4] * b[3]} - {a[5] * b[2]} + {a[6] * b[1]} + {a[7] * b[0]} }

            (vector
              ;0
              (- (* (vector-ref a 0) (vector-ref b 0))
                 (* (vector-ref a 1) (vector-ref b 1))
                 (* (vector-ref a 2) (vector-ref b 2))
                 (* (vector-ref a 3) (vector-ref b 3))
                 (* (vector-ref b 4) (vector-ref a 4))
                 (* (vector-ref b 5) (vector-ref a 5))
                 (* (vector-ref b 6) (vector-ref a 6))
                 (* (vector-ref b 7) (vector-ref a 7)))
              ;1
              (- (+ (+ (- (- (+ (+ (* (vector-ref a 0) (vector-ref b 1))
                                   (* (vector-ref a 1) (vector-ref b 0)))
                                (* (vector-ref a 2) (vector-ref b 3)))
                             (* (vector-ref a 3) (vector-ref b 2)))
                          (* (vector-ref b 4) (vector-ref a 5)))
                       (* (vector-ref b 5) (vector-ref a 4)))
                    (* (vector-ref b 6) (vector-ref a 7)))
                 (* (vector-ref b 7) (vector-ref a 6)))
              ;2
              (+ (+ (- (- (+ (+ (- (* (vector-ref a 0) (vector-ref b 2))
                                   (* (vector-ref a 1) (vector-ref b 3)))
                                (* (vector-ref a 2) (vector-ref b 0)))
                             (* (vector-ref a 3) (vector-ref b 1)))
                          (* (vector-ref b 4) (vector-ref a 6)))
                       (* (vector-ref b 5) (vector-ref a 7)))
                    (* (vector-ref b 6) (vector-ref a 4)))
                 (* (vector-ref b 7) (vector-ref a 5)))
              ;3
              (+ (- (+ (- (+ (- (+ (* (vector-ref a 0) (vector-ref b 3))
                                   (* (vector-ref a 1) (vector-ref b 2)))
                                (* (vector-ref a 2) (vector-ref b 1)))
                             (* (vector-ref a 3) (vector-ref b 0)))
                          (* (vector-ref b 4) (vector-ref a 7)))
                       (* (vector-ref b 5) (vector-ref a 6)))
                    (* (vector-ref b 6) (vector-ref a 5)))
                 (* (vector-ref b 7) (vector-ref a 4)))
              ;4
              (+ (+ (+ (+ (- (- (- (* (vector-ref b 4) (vector-ref a 0))
                                   (* (vector-ref b 5) (vector-ref a 1)))
                                (* (vector-ref b 6) (vector-ref a 2)))
                             (* (vector-ref b 7) (vector-ref a 3)))
                          (* (vector-ref a 4) (vector-ref b 0)))
                       (* (vector-ref a 5) (vector-ref b 1)))
                    (* (vector-ref a 6) (vector-ref b 2)))
                 (* (vector-ref a 7) (vector-ref b 3)))
              ;5
              (+ (- (+ (- (- (+ (+ (* (vector-ref b 4) (vector-ref a 1))
                                   (* (vector-ref b 5) (vector-ref a 0)))
                                (* (vector-ref b 6) (vector-ref a 3)))
                             (* (vector-ref b 7) (vector-ref a 2)))
                          (* (vector-ref a 4) (vector-ref b 1)))
                       (* (vector-ref a 5) (vector-ref b 0)))
                    (* (vector-ref a 6) (vector-ref b 3)))
                 (* (vector-ref a 7) (vector-ref b 2)))
              ;6
              (- (+ (+ (- (+ (+ (- (* (vector-ref b 4) (vector-ref a 2))
                                   (* (vector-ref b 5) (vector-ref a 3)))
                                (* (vector-ref b 6) (vector-ref a 0)))
                             (* (vector-ref b 7) (vector-ref a 1)))
                          (* (vector-ref a 4) (vector-ref b 2)))
                       (* (vector-ref a 5) (vector-ref b 3)))
                    (* (vector-ref a 6) (vector-ref b 0)))
                 (* (vector-ref a 7) (vector-ref b 1)))
              ;7
              (+ (+ (- (- (+ (- (+ (* (vector-ref b 4) (vector-ref a 3))
                                   (* (vector-ref b 5) (vector-ref a 2)))
                                (* (vector-ref b 6) (vector-ref a 1)))
                             (* (vector-ref b 7) (vector-ref a 0)))
                          (* (vector-ref a 4) (vector-ref b 3)))
                       (* (vector-ref a 5) (vector-ref b 2)))
                    (* (vector-ref a 6) (vector-ref b 1)))
                 (* (vector-ref a 7) (vector-ref b 0)))) )
          (else
           (error 'cross-product "undefined for this length" alen blen) ) ) ) ) )

) ;(math-utils vector)