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

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

(module (math-utils vector)

(;export
  vector-sum
  absolute-magnitude
  cosine-similarity
  vector-apply
  vector-reduce
  vector-compare vector-compare2
  dot-product
  cross-product
  softmax softmax*)

(import scheme)
(import (chicken base))
(import (chicken fixnum))
(import (chicken flonum))
(import (chicken type))
(import math-utils)

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

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

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

(: vector-sum (number-vector -> number))

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

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

(: vector-reduce ((number number -> number) number number-vector #!optional fixnum fixnum -> number))
(: vector-apply (nary-numeric #!rest number-vector -> number-vector))

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

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

(: softmax (number-vector -> number-vector))
(: softmax* (number-vector #!optional real -> number-vector))

;;

(define-inline (sqrt-exact x) (inexact->exact (sqrt x)))

(define (sumsqr acc elm) (+ (sqr elm) acc))

;;

(define (vector-reduce f id v #!optional (st 0) (ed (vector-length v)))
  (let loop ((i st) (acc id))
    (if (fx= i ed) acc
        (loop (fx+ i 1) (f acc (vector-ref v i)))) ) )

(define (vector-sum v) (vector-reduce + 0 v))

(define (absolute-magnitude v) (sqrt-exact (vector-reduce sumsqr 0 v)))

;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) (sumsqr maga ai) (sumsqr magb bi)) ) ) ) ) )

; assumes all same length
(define-inline (*vector-apply/1 f a)
  (let* ((len (vector-length a))
         (r (make-vector len)) )
    (do ((i 0 (fx+ i 1)))
         ((fx= i len) r)
      (set! (vector-ref r i) (f (vector-ref a i))) ) ) )

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

(define (vector-compare2 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 (fx= i alen) 0
              (let ((d (- (vector-ref a i) (vector-ref b i))))
                (if (not (= 0 d)) d
                    (loop (fx+ i 1)) ) ) ) ) ) ) )

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

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

;0 -> error
;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)))
    (cond ((or (fx= 0 alen) (fx= 0 blen))
            0)
          ((fx= alen blen)
            (case alen
              ((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) ) ) )
          (else
            (error 'cross-product "undefined for varying length" alen blen)) ) ) )

;1 ~= summ(vmax(V)) but 1 ?= summ(V), normalize
(define (softmax v)
  (let* ((ve (*vector-apply/1 exp v))
         (sum (vector-sum ve)) )
    (*vector-apply/1 (cut / <> sum) ve) ) )

#|
We can also use T values less than 1 to INCREASE the distance between the
values entering softmax, and in the limit as T->0, we max the maximum value be
infinitely far from the rest as it enters the exponent function. In code, this
must be handled as a special case that finds the maximum output instead of
selecting an output probabilistically.

So the legal values of T range from 0 (pick the maximum) through infinity
(where all values would have equal probabilities because they are all 0 as
they enter softmax).

temp = 1/alpha
|#

(define-constant EPSILON (inexact->exact flonum-epsilon))
(define-inline (near? m n) (or (= m n) (<= (abs (- m n)) EPSILON)))

(define (softmax* v #!optional (t 1))
  (assert (not (near? 0 t)))
  (softmax (if (= t 1) v (*vector-apply/1 (cut / <> t) v))) )

) ;(math-utils vector)