~dieggsy/numchi

34960576e81b6f68f511e6c1a239cdd4628764f2 — dieggsy 1 year, 10 months ago 8e00449
Add tests, fix erroring trig functions
2 files changed, 265 insertions(+), 157 deletions(-)

M math.scm
M tests/run.scm
M math.scm => math.scm +41 -24
@@ 29,6 29,8 @@
 add
 conj)

;; TODO: probably return real output for real arrays, like NumPy

@(==== "Trigonometric functions")

(define (::sin a #!key dtype)


@@ 81,10 83,11 @@
       (array-map asin a)
       (or dtype
           (cond ((gena? a) gen)
                 ((or (complex-array? a)
                      (flonum-array? a))
                 ((complex-array? a)
                  (array-dtype a))
                 (else f64))))
                 ((f32a? a) c64)
                 ((f64a? a) c128)
                 (else c128))))
      (asin a)))

(define arcsin


@@ 100,10 103,11 @@
       (array-map acos a)
       (or dtype
           (cond ((gena? a) gen)
                 ((or (complex-array? a)
                      (flonum-array? a))
                 ((complex-array? a)
                  (array-dtype a))
                 (else f64))))
                 ((f32a? a) c64)
                 ((f64a? a) c128)
                 (else c128))))
      (acos a)))

(define arccos


@@ 136,6 140,7 @@
       (@to "(or array number)"))
  ::atan)

;; TODO: add note about complex numbers and arc* procedures
@("Standard trigonometric functions."

  "Similar: [[https://numpy.org/doc/stable/reference/routines.math.html#trigonometric-functions|Trigonometric functions]]")


@@ 190,7 195,7 @@
         ((or (c64a? a) (c64a? b)) c64)
         ((or (f32a? a) (f32a? b)) f32)
         (else f64))))
(define-generic2 (%atan2 (<array> a) (<number> x))
(define-generic (%atan2 (<array> a) (<number> x))
  (maybe-copy
   (array-map
    (lambda (n)


@@ 201,12 206,24 @@
         ((cplxnum? x) c128)
         ((f32a? a) f32)
         (else f64))))
(define-generic (%atan2 (<number> x) (<array> a))
  (maybe-copy
   (array-map
    (lambda (n)
      (atan x n))
    a)
   (cond ((gena? a) gen)
         ((complex-array? a) (array-body a))
         ((cplxnum? x) c128)
         ((f32a? a) f32)
         (else f64))))
(define-generic (%atan2 (<number> x) (<number> y))
  (atan x y))

(define arctan2
  @(fn (@sig (arctan2 a b #!key type))
       (@to "(or array number)")) atan2)
       (@to "(or array number)"))
  atan2)

@("Computes (angle (make-rectangular a b))."
  "Similar: [[https://numpy.org/doc/stable/reference/generated/numpy.atan2.html|numpy.atan2]]")


@@ 316,17 333,17 @@
  (log (+ x (sqrt (+ (expt x 2) 1)))))

(define (asinh a #!key dtype)
  @((@sig (asin a #!key dtype))
    (@to "(or array number)"))
  @((@to "(or array number)"))
  (if (array? a)
      (maybe-copy
       (array-map %asinh a)
       (or dtype
           (cond ((gena? a) gen)
                 ((or (complex-array? a)
                      (flonum-array? a))
                 ((complex-array? a)
                  (array-dtype a))
                 (else f64))))
                 ((f32a? a) c64)
                 ((f64a? a) c128)
                 (else c128))))
      (%asinh a)))

(define arcsinh


@@ 338,17 355,17 @@
  (log (+ x (sqrt (- (expt x 2) 1)))))

(define (acosh a #!key dtype)
  @((@sig (acos a #!key dtype))
    (@to "(or array number)"))
  @((@to "(or array number)"))
  (if (array? a)
      (maybe-copy
       (array-map %acosh a)
       (or dtype
           (cond ((gena? a) gen)
                 ((or (complex-array? a)
                      (flonum-array? a))
                 ((complex-array? a)
                  (array-dtype a))
                 (else f64))))
                 ((f32a? a) c64)
                 ((f64a? a) c128)
                 (else c128))))
      (%acosh a)))

(define arccosh


@@ 360,23 377,23 @@
  (/ (log (/ (+ 1 x) (- 1 x))) 2))

(define (atanh a #!key dtype)
  @((@sig (atan a [b] #!key dtype))
    (@to "(or array number)"))
  @((@to "(or array number)"))
  (if (array? a)
      (maybe-copy
       (array-map %atanh a)
       (or dtype
           (cond ((gena? a) gen)
                 ((or (complex-array? a)
                      (flonum-array? a))
                 ((complex-array? a)
                  (array-dtype a))
                 (else f64))))
                 ((f32a? a) c64)
                 ((f64a? a) c128)
                 (else c128))))
      (%atanh a)))

(define arctanh
  @(fn (@sig (arctanh a #!key dtype))
       (@to "(or array number)"))
  %atanh)
  atanh)

@("Standard hyperbolic functions."


M tests/run.scm => tests/run.scm +224 -133
@@ 1,9 1,13 @@
;; TODO: dtype checking

;; TODO: for every procedure, probably should make sure to check behavior for:
;; flonum array, integer array, complex array at the very least. Should
;; probably test negative numbers too....
(module main ()
  (import r7rs
          (chicken base)
          pyffi
          numchi
          (prefix numchi nc:)
          test
          srfi-1
          chicken.format


@@ 14,10 18,10 @@
          srfi-13
          srfi-133)

  (current-test-verbosity #f)
  (current-test-verbosity #t)
  (define random-test-number 50)

  (array-pretty-print #f)
  (nc:array-pretty-print #f)

  (py-start)
  (py-import "numpy")


@@ 37,23 41,23 @@
    (syntax-rules ()
      [(_ name a1 a2)
       (test-assert name
         (equal? (array->list* a1)
         (equal? (nc:array->list* a1)
                 (ndarray->list a2)))]
      [(_ a1 a2)
       (test-assert
           (equal? (array->list* a1)
           (equal? (nc:array->list* a1)
                   (ndarray->list a2)))]))

  (define (arrays-fl= a1 a2)
    (and
     (equal? (array-shape a1)
     (equal? (nc:array-shape a1)
             (ndarray-shape a2))
     (every
     identity
     (map (lambda (a b)
            (< (- a b) (current-test-epsilon)))
          (flatten (array->list* a1))
          (flatten (ndarray->list a2))))))
      identity
      (map (lambda (a b)
             (< (- a b) (current-test-epsilon)))
           (flatten (nc:array->list* a1))
           (flatten (ndarray->list a2))))))

  (define-syntax test-array/fl
    (syntax-rules ()


@@ 75,13 79,13 @@
           (base-shape (list-tabulate
                        dims
                        (lambda i (fx+ (pseudo-random-integer max-dim-elems) 2))))
           (shape (domain base-shape))
           (shape (nc:domain base-shape))
           (base-lst (list-tabulate (apply * base-shape)
                                    (lambda i (fx+ (if nonzero 1 0)
                                                 (pseudo-random-integer
                                                  (fx-
                                                   (if small 10 500)
                                                   (if nonzero 1 0)))))))
                                    (lambda _ (fx+ (if nonzero 1 0)
                                                   (pseudo-random-integer
                                                    (fx-
                                                     (if small 10 500)
                                                     (if nonzero 1 0)))))))
           (python-base-shape (string-append "("
                                             (string-join (map ->string base-shape) ",")
                                             ")"))


@@ 89,11 93,28 @@
                                           (string-join (map ->string base-lst) ",")
                                           "]")))
      (values
       (list->array base-lst shape s64-storage-class)
       (nc:list->array base-lst shape nc:s64)
       (py-eval (format "numpy.array(~a).reshape(~a)"
                        python-base-lst
                        python-base-shape)))))

  (define (random-like a #!key small nonzero)
    (let* ((shape (nc:array-shape a))
           (arr (nc:fromfunction
                 (lambda _ (fx+ (if nonzero 1 0)
                                (pseudo-random-integer
                                 (fx-
                                  (if small 10 500)
                                  (if nonzero 1 0)))))
                 shape
                 #:dtype (nc:array-dtype a))))
      (values
       arr
       (py-eval
        (format "numpy.array(~a).reshape(~a)"
                (string-append "[" (string-join (map ->string (nc:array->list arr)) ",") "]")
                (string-append "(" (string-join (map ->string (vector->list shape)) ",") ")"))))))

  (test-begin "numchi")

  (test-group "Basic test"


@@ 103,129 124,199 @@

    ;; Test array->list*, since this'll be the basis for array comparison across
    ;; the two libraries
    (test-array "basic" (list->array '(1 2 3 4) (domain #(2 2)))
    (test-array "basic" (nc:list->array '(1 2 3 4) (nc:domain #(2 2)))
                (py-eval "numpy.array([1,2,3,4]).reshape((2,2))"))

    (test-array "basic" (list->array '(1 2 3 4 5 6 7 8 9 10 11 12) (domain #(3 2 2)))
    (test-array "basic" (nc:list->array '(1 2 3 4 5 6 7 8 9 10 11 12) (nc:domain #(3 2 2)))
                (py-eval "numpy.array([1,2,3,4,5,6,7,8,9,10,11,12]).reshape((3,2,2))")))


  (test-group "Array procedures"
    (define-pyfun ("numpy.transpose" ndarray-transpose) arr axes)
    (define-pyfun ("numpy.diagonal" ndarray-diagonal) arr offset axis1 axis2)
    (define-pyfun ("numpy.trace" ndarray-trace) arr offset axis1 axis2)
    (define-pyfun ("numpy.clip" ndarray-clip) arr min max)
    (define-pyfun ("numpy.ravel" ndarray-flatten) arr)
    (define-pyfun ("numpy.reshape" ndarray-reshape) arr shape)
  (test-group "array.scm"
    (test-group "Array procedures"
      (define-pyfun ("numpy.transpose" ndarray-transpose) arr axes)
      (define-pyfun ("numpy.diagonal" ndarray-diagonal) arr offset axis1 axis2)
      (define-pyfun ("numpy.trace" ndarray-trace) arr offset axis1 axis2)
      (define-pyfun ("numpy.clip" ndarray-clip) arr min max)
      (define-pyfun ("numpy.ravel" ndarray-flatten) arr)
      (define-pyfun ("numpy.reshape" ndarray-reshape) arr shape)

      (do ((i 0 (add1 i)))
          ((fx= i random-test-number))
        (let*-values ([(arr py-arr) (random-array)]
                      [(ra1 ra2) (let ((ra1 (pseudo-random-integer (nc:array-ndim arr))))
                                   (let loop ((ra2 (pseudo-random-integer (nc:array-ndim arr))))
                                     (if (fx= ra2 ra1)
                                         (loop (pseudo-random-integer (nc:array-ndim arr)))
                                         (values ra1 ra2))))]
                      [(random-offset)
                       (fx- (pseudo-random-integer
                             (sub1 (fx+ (vector-ref (nc:array-shape arr) ra1)
                                        (vector-ref (nc:array-shape arr) ra2))))
                            (sub1 (vector-ref (nc:array-shape arr) ra1)))]
                      [(ndim) (nc:array-ndim arr)]
                      [(random-permutation)
                       (random-permutation!
                        (vector-unfold
                         (lambda (i x)
                           (values x (fx- x 1)))
                         ndim
                         (fx- ndim 1)))])
          (test-array "transpose"
                      (nc:array-transpose arr
                                          axes:
                                          random-permutation)
                      (ndarray-transpose py-arr
                                         random-permutation))
          (test-array "flatten"
                      (nc:array-flatten arr)
                      (ndarray-flatten py-arr))
          (let* ((bounds (nc:interval-upper-bounds->vector (nc:array-domain arr))))
            (vector-set! bounds 1 (* (vector-ref bounds 0)
                                     (vector-ref bounds 1)))
            (test-array "reshape"
                        ;; Map identity to the array so that we're testing
                        ;; numchi's general reshape and not
                        ;; specialized-array-reshape
                        (nc:array-reshape (nc:array-map identity arr)
                                          (subvector bounds 1))
                        (ndarray-reshape py-arr
                                         (subvector bounds 1))))
          (test-array "diagonal"
                      (nc:array-diagonal arr
                                         offset: random-offset
                                         axis1: ra1
                                         axis2: ra2)
                      (ndarray-diagonal py-arr
                                        random-offset
                                        ra1
                                        ra2))
          (let ((trace (nc:array-trace arr
                                       offset: random-offset
                                       axis1: ra1
                                       axis2: ra2))
                (py-trace (ndarray-trace py-arr
                                         random-offset
                                         ra1
                                         ra2)))
            (if (number? trace)
                (test "trace" trace (numpy-item py-trace))
                (test-array "trace" trace py-trace )))))

      (let*-values ([(arr py-arr) (random-array)])
        (test-array "array-clip"
                    (nc:array-clip arr min: 10 max: 40)
                    (ndarray-clip py-arr 10 40))))

    (do ((i 0 (add1 i)))
        ((fx= i random-test-number))
    (test-group "axiswise procedures"
      (define-pyfun ("numpy.amax" ndarray-max) arr axis)
      (define-pyfun ("numpy.amin" ndarray-min) arr axis)
      (define-pyfun ("numpy.sum" ndarray-sum) arr axis)
      (define-pyfun ("numpy.mean" ndarray-mean) arr axis)
      (define-pyfun ("numpy.var" ndarray-var) arr axis)
      (define-pyfun ("numpy.std" ndarray-std) arr axis)
      (define-pyfun ("numpy.prod" ndarray-prod) arr axis)
      (define-pyfun ("numpy.all" ndarray-all) arr axis)
      (define-pyfun ("numpy.any" ndarray-any) arr axis)
      (do ((i 0 (add1 i)))
          ((fx= i random-test-number))
        (let*-values ([(arr py-arr) (random-array)]
                      [(random-axis) (pseudo-random-integer (nc:array-ndim arr))])
          (test-array "max"
                      (nc:array-max arr axis: random-axis)
                      (ndarray-max py-arr random-axis))
          (test-array "min"
                      (nc:array-min arr axis: random-axis)
                      (ndarray-min py-arr random-axis))
          (test-array "sum"
                      (nc:array-sum arr axis: random-axis)
                      (ndarray-sum py-arr random-axis))
          (test-array/fl "mean"
                         (nc:array-map inexact (nc:array-mean arr axis: random-axis))
                         (ndarray-mean py-arr random-axis))
          (test-array/fl "var"
                         (nc:array-map inexact (nc:array-var arr axis: random-axis))
                         (ndarray-var py-arr random-axis))
          (test-array/fl "std"
                         (nc:array-map inexact (nc:array-std arr axis: random-axis))
                         (ndarray-std py-arr random-axis)))
        (let*-values ([(arr py-arr) (random-array small: #t)]
                      [(random-axis) (pseudo-random-integer (nc:array-ndim arr))])
          (test-array "prod"
                      (nc:array-prod arr axis: random-axis)
                      (ndarray-prod py-arr random-axis))))
      ;; This test might be kinda useless, but it could at least catch silly
      ;; mistakes.
      (let*-values ([(arr py-arr) (random-array nonzero: #t)]
                    [(random-axis) (pseudo-random-integer (nc:array-ndim arr))])
        (test-array "all"
                    (nc:array-all arr axis: random-axis)
                    (ndarray-all py-arr random-axis))
        (test-array "any"
                    (nc:array-any* arr axis: random-axis)
                    (ndarray-any py-arr random-axis)))))

  ;; (test-group "array-creation.scm")

  (test-group "math.scm"
    ;; Again, these are just to ensure we catch silly mistakes, no crazy
    ;; dimension manipulation here
    (test-group "Trigonometric functions"
      (define-pyfun ("numpy.sin" ndarray-sin) arr)
      (define-pyfun ("numpy.cos" ndarray-cos) arr)
      (define-pyfun ("numpy.tan" ndarray-tan) arr)
      (define-pyfun ("numpy.arcsin" ndarray-arcsin) arr)
      (define-pyfun ("numpy.arccos" ndarray-arccos) arr)
      (define-pyfun ("numpy.arctan" ndarray-arctan) arr)
      (define-pyfun ("numpy.arctan2" ndarray-arctan2) arr brr)
      (define-pyfun ("numpy.degrees" ndarray-degrees) arr)
      (define-pyfun ("numpy.radians" ndarray-radians) arr)
      (let*-values ([(arr py-arr) (random-array)]
                    [(ra1 ra2) (let ((ra1 (pseudo-random-integer (array-ndim arr))))
                                 (let loop ((ra2 (pseudo-random-integer (array-ndim arr))))
                                   (if (fx= ra2 ra1)
                                       (loop (pseudo-random-integer (array-ndim arr)))
                                       (values ra1 ra2))))]
                    [(random-offset)
                     (fx- (pseudo-random-integer
                           (sub1 (fx+ (vector-ref (array-shape arr) ra1)
                                      (vector-ref (array-shape arr) ra2))))
                          (sub1 (vector-ref (array-shape arr) ra1)))]
                    [(ndim) (array-ndim arr)]
                    [(random-permutation)
                      (random-permutation!
                       (vector-unfold
                        (lambda (i x)
                          (values x (fx- x 1)))
                        ndim
                        (fx- ndim 1)))])
        (test-array "transpose"
                    (array-transpose arr
                                     axes:
                                     random-permutation)
                    (ndarray-transpose py-arr
                                      random-permutation))
        (test-array "flatten"
                    (array-flatten arr)
                    (ndarray-flatten py-arr))
        (let* ((bounds (interval-upper-bounds->vector (array-domain arr))))
          (vector-set! bounds 1 (* (vector-ref bounds 0)
                                   (vector-ref bounds 1)))
          (test-array "reshape"
                      ;; Map identity to the array so that we're testing
                      ;; numchi's general reshape and not
                      ;; specialized-array-reshape
                      (array-reshape (array-map identity arr)
                                     (subvector bounds 1))
                      (ndarray-reshape py-arr
                                       (subvector bounds 1))))
        (test-array "diagonal"
                    (array-diagonal arr
                                    offset: random-offset
                                    axis1: ra1
                                    axis2: ra2)
                    (ndarray-diagonal py-arr
                                      random-offset
                                      ra1
                                      ra2))
        (let ((trace (array-trace arr
                                  offset: random-offset
                                  axis1: ra1
                                  axis2: ra2))
              (py-trace (ndarray-trace py-arr
                                       random-offset
                                       ra1
                                       ra2)))
          (if (number? trace)
              (test "trace" trace (numpy-item py-trace))
              (test-array "trace" trace py-trace )))))

    (let*-values ([(arr py-arr) (random-array)])
      (test-array "array-clip" (array-clip arr min: 10 max: 40)
                  (ndarray-clip py-arr 10 40))))

  (test-group "axiswise procedures"
    (define-pyfun ("numpy.amax" ndarray-max) arr axis)
    (define-pyfun ("numpy.amin" ndarray-min) arr axis)
    (define-pyfun ("numpy.sum" ndarray-sum) arr axis)
    (define-pyfun ("numpy.mean" ndarray-mean) arr axis)
    (define-pyfun ("numpy.var" ndarray-var) arr axis)
    (define-pyfun ("numpy.std" ndarray-std) arr axis)
    (define-pyfun ("numpy.prod" ndarray-prod) arr axis)
    (define-pyfun ("numpy.all" ndarray-all) arr axis)
    (define-pyfun ("numpy.any" ndarray-any) arr axis)
    (do ((i 0 (add1 i)))
        ((fx= i random-test-number))
                    [(brr py-brr) (random-like arr)]
                    [(x) (pseudo-random-integer 500)]
                    [(y) (pseudo-random-integer 500)])
        (test-array/fl "sin" (nc:sin arr) (ndarray-sin py-arr))
        (test "sin (number)" (nc:sin x) (numpy-item (ndarray-sin x)))
        (test-array/fl "cos" (nc:cos arr) (ndarray-cos py-arr))
        (test "cos (number)" (nc:cos x) (numpy-item (ndarray-cos x)))
        (test-array/fl "tan" (nc:tan arr) (ndarray-tan py-arr))
        (test "tan (number)" (nc:tan x) (numpy-item (ndarray-tan x)))
        ;; Pyffi can't handle complex numbers, I think, so let's hold off on
        ;; these for now, just make sure they don't error on arrays for now
        (test-assert "arcsin" (nc:arcsin arr))
        (test-assert "arccos" (nc:arccos arr))
        ;; (test-array/fl "arcsin" (nc:arcsin arr) (ndarray-arcsin py-arr))
        ;; (test "arcsin (number)" (nc:arcsin x) (numpy-item (ndarray-arcsin x)))
        ;; (test-array/fl "arccos" (nc:arccos arr) (ndarray-arccos py-arr))
        ;; (test "arccos (number)" (nc:arccos x) (numpy-item (ndarray-arccos x)))
        (test-array/fl "arctan" (nc:arctan arr) (ndarray-arctan py-arr))
        (test "arctan (number)" (nc:arctan x) (numpy-item (ndarray-arctan x)))
        (test-array/fl "arctan2" (nc:arctan arr brr) (ndarray-arctan2 py-arr py-brr))
        (test-array/fl "arctan2 (array number)" (nc:arctan arr x) (ndarray-arctan2 py-arr x))
        (test-array/fl "arctan2 (number array)" (nc:arctan x arr) (ndarray-arctan2 x py-arr))
        (test "arctan2 (number number)" (nc:arctan x y) (numpy-item (ndarray-arctan2 x y)))
        (test-array/fl "degrees" (nc:degrees arr) (ndarray-degrees py-arr))
        (test-array/fl "radians" (nc:radians arr) (ndarray-radians py-arr))))

    (test-group "Hyperbolic functions"
      (define-pyfun ("numpy.sinh" ndarray-sinh) arr)
      (define-pyfun ("numpy.cosh" ndarray-cosh) arr)
      (define-pyfun ("numpy.tanh" ndarray-tanh) arr)
      (define-pyfun ("numpy.arcsinh" ndarray-arcsinh) arr)
      (define-pyfun ("numpy.arccosh" ndarray-arccosh) arr)
      (define-pyfun ("numpy.arctanh" ndarray-arctanh) arr)
      (let*-values ([(arr py-arr) (random-array)]
                    [(random-axis) (pseudo-random-integer (array-ndim arr))])
        (test-array "max" (array-max arr axis: random-axis)
                    (ndarray-max py-arr random-axis))
        (test-array "min" (array-min arr axis: random-axis)
                    (ndarray-min py-arr random-axis))
        (test-array "sum" (array-sum arr axis: random-axis)
                    (ndarray-sum py-arr random-axis))
        (test-array/fl "mean"
                       (array-map inexact (array-mean arr axis: random-axis))
                       (ndarray-mean py-arr random-axis))
        (test-array/fl "var"
                       (array-map inexact (array-var arr axis: random-axis))
                       (ndarray-var py-arr random-axis))
        (test-array/fl "std"
                       (array-map inexact (array-std arr axis: random-axis))
                       (ndarray-std py-arr random-axis)))
      (let*-values ([(arr py-arr) (random-array small: #t)]
                    [(random-axis) (pseudo-random-integer (array-ndim arr))])
        (test-array "prod"
                    (array-prod arr axis: random-axis)
                    (ndarray-prod py-arr random-axis))))
    (let*-values ([(arr py-arr) (random-array nonzero: #t)]
                  [(random-axis) (pseudo-random-integer (array-ndim arr))])
      (test-array "all"
                  (array-all arr axis: random-axis)
                  (ndarray-all py-arr random-axis))
      (test-array "any"
                  (array-any* arr axis: random-axis)
                  (ndarray-any py-arr random-axis))))
                    [(x) (pseudo-random-integer 500)])
        (test-array/fl "sinh" (nc:sinh arr) (ndarray-sinh py-arr))
        (test "sinh (number)" (nc:sinh x) (numpy-item (ndarray-sinh x)))
        (test-array/fl "cosh" (nc:cosh arr) (ndarray-cosh py-arr))
        (test "cosh (number)" (nc:cosh x) (numpy-item (ndarray-cosh x)))
        (test-array/fl "tanh" (nc:tanh arr) (ndarray-tanh py-arr))
        (test "tanh (number)" (nc:tanh x) (numpy-item (ndarray-tanh x)))

        ;; Ensure these don't error, we'll have to handle complex stuff later.
        (test-assert "arcsinh" (nc:arcsinh arr))
        (test-assert "arccosh" (nc:arccosh arr))
        (test-assert "arctanh" (nc:arctanh arr)))))

  (test-end "numchi")
  (py-stop)