;;; Unified Fusion-Based Push-Pull Arrays
;;; An implementation of lazy arrays with fusion.

(module fusion-arrays

(;;; Core API
 make-array array-from-list array-size array-ref array->list
 compute! force!
 
 ;;; Array operations (all lazy by default)
 array+ array- array* array/ array-pow array-min array-max
 array-scale array-add-scalar array-sub-scalar array-div-scalar
 array-negate array-abs array-sqrt array-exp array-log
 array-sin array-cos array-tan array-floor array-ceiling
 array> array< array= array-and array-or array-not
 
 ;;; Core functional operations
 array-map array-reduce array-fold array-slice array-concat
 array-filter array-zip array-take array-drop array-diff
 array-cumsum array-every? array-any? array-count array-select
 array-where array-unique array-group-by

 array-map* array-fold* array-reduce* array-filter* array-zip*

 ;;; Statistical operations
 array-sum array-mean array-var array-std array-min-val array-max-val
 array-moving-average array-moving-std array-moving-min array-moving-max
 
 ;;; Signal processing
 array-convolve array-gaussian-filter array-correlate
 array-fft array-ifft array-autocorrelation array-integrate
 
 ;;; Utilities
 array-pipeline array-linspace array-random-normal array-meshgrid
 array-zeros array-ones array-constant array-range
 
 ;;; Performance and debugging
 ;fusion-stats clear-fusion-stats! profile-computation
 array-complexity fusion-plan pp-fusion)

(import scheme (chicken base) (chicken foreign)
        (chicken memory) (chicken bitwise) (chicken locative)
        datatype matchable srfi-1 srfi-4 srfi-69
        (prefix random-mtzig random:))

        

;;;=============================================================================
;;; Core fusion array types
;;;=============================================================================

;;; All arrays are either materialized data or lazy fusion expressions
(define-datatype farray farray?
  ;; Materialized arrays (numeric data)
  (materialized-array
   (data (lambda (x) (or (vector? x) (f32vector? x) (f64vector? x)
                         (s32vector? x) (u32vector? x) (s64vector? x)))) ; typed vector (f64vector, f32vector, etc.)
   (size exact-integer?)    ; number of elements
   (element-type symbol?)   ; 'f64, 'f32, 's64, 's32, 'u32, 'generic
   (metadata list?))        ; optional metadata (sampling rate, units, etc.)
  
  ;; Lazy fusion expressions (no data until materialized)
  (fusion-expr
   (operation symbol?)      ; 'add, 'mul, 'sin, 'map, 'slice, 'concat, etc.
   (operands list?)         ; list of operand farrays
   (params list?)           ; operation parameters (scalars, functions, etc.)
   (result-type symbol?)    ; inferred result type
   (result-size exact-integer?) ; computed result size
   (metadata list?))        ; propagated metadata
  
  ;; Lazy reduction expressions (reduce array to scalar)
  (fusion-reduction
   (operation symbol?)      ; 'sum, 'mean, 'var, 'std, 'min, 'max, 'fold, etc.
   (operand farray?)        ; single operand array
   (params list?)           ; reduction parameters (function, initial value, etc.)
   (result-type symbol?)    ; scalar result type
   (metadata list?)))       ; propagated metadata

         
  
;;;=============================================================================
;;; Type inference
;;;=============================================================================

(define *type-hierarchy* '(s32 u32 s64 f32 f64 generic))
(define *type-ranks* '((s32 . 0) (u32 . 1) (s64 . 2) (f32 . 3) (f64 . 4) (generic . 5)))

(define (type-rank type)
  "Get numeric rank of type for promotion"
  (cdr (assoc type *type-ranks*)))

(define (promote-types type1 type2)
  "Find the promoted type for two operand types"
  (let ((rank1 (type-rank type1))
        (rank2 (type-rank type2)))
    (list-ref *type-hierarchy* (max rank1 rank2))))

(define (infer-result-type op operand-types)
  "Infer result type for operation on given operand types"
  (case op
    ;; Arithmetic operations preserve highest type
    ((add sub mul pow min max)
     (fold promote-types (car operand-types) (cdr operand-types)))
    
    ;; Division always promotes to floating point
    ((div)
     (let ((promoted (fold promote-types (car operand-types) (cdr operand-types))))
       (case promoted
         ((s32 u32 s64) 'f64)
         ((f32) 'f32)
         (else promoted))))
    
    ;; Transcendental functions always produce floating point
    ((exp log sin cos tan)
     (let ((input-type (car operand-types)))
       (case input-type
         ((s32 u32 s64) 'f64)
         ((f32) 'f32)
         (else 'f64))))

    ;; sqrt can return floating point or complex
    ((sqrt) 'generic)
     
    
    ;; Comparison operations produce boolean (represented as f64)
    ((gt lt eq and or not) 'f64)
    
    ;; Statistical operations
    ((sum mean) (car operand-types))
    ((var std) (case (car operand-types)
                 ((s32 u32 s64) 'f64)
                 (else (car operand-types))))
    
    ;; Map preserves or promotes based on function result
    ((map) (car operand-types))  ; Conservative estimate
    
    ;; Default: preserve first operand type
    (else (car operand-types))))

;;;=============================================================================
;;; Core fusion engine
;;;=============================================================================

;;; Evaluates expressions index-by-index without creating intermediate arrays

(define (evaluate-fusion-at-index expr index)
  "Evaluate fusion expression at specific index"
   
  (match expr
    ;; Base case: materialized array
    (($ farray 'materialized-array data size element-type _)
     (typed-vector-ref data index element-type))
    
    ;; Arithmetic operations
    (($ farray 'fusion-expr 'add operands _ _ _ _)
     (apply + (map (lambda (op) (evaluate-fusion-at-index op index)) operands)))
    
    (($ farray 'fusion-expr 'sub (operand1 operand2) _ _ _ _)
     (- (evaluate-fusion-at-index operand1 index)
        (evaluate-fusion-at-index operand2 index)))
    
    (($ farray 'fusion-expr 'mul operands _ _ _ _)
     (apply * (map (lambda (op) (evaluate-fusion-at-index op index)) operands)))
    
    (($ farray 'fusion-expr 'div (operand1 operand2) _ _ _ _)
     (/ (evaluate-fusion-at-index operand1 index)
        (evaluate-fusion-at-index operand2 index)))
    
    (($ farray 'fusion-expr 'pow (operand1 operand2) _ _ _ _)
     (expt (evaluate-fusion-at-index operand1 index)
           (evaluate-fusion-at-index operand2 index)))
    
    ;; Scalar operations
    (($ farray 'fusion-expr 'scale (operand) (scalar) _ _ _)
     (* (evaluate-fusion-at-index operand index) scalar))
    
    (($ farray 'fusion-expr 'add-scalar (operand) (scalar) _ _ _)
     (+ (evaluate-fusion-at-index operand index) scalar))
    
    (($ farray 'fusion-expr 'sub-scalar (operand) (scalar) _ _ _)
     (- (evaluate-fusion-at-index operand index) scalar))
    
    (($ farray 'fusion-expr 'div-scalar (operand) (scalar) _ _ _)
     (/ (evaluate-fusion-at-index operand index) scalar))
    
    ;; Unary operations
    (($ farray 'fusion-expr 'negate (operand) _ _ _ _)
     (- (evaluate-fusion-at-index operand index)))
    
    (($ farray 'fusion-expr 'abs (operand) _ _ _ _)
     (abs (evaluate-fusion-at-index operand index)))
    
    (($ farray 'fusion-expr 'sqrt (operand) _ _ _ _)
     (sqrt (evaluate-fusion-at-index operand index)))
    
    (($ farray 'fusion-expr 'exp (operand) _ _ _ _)
     (exp (evaluate-fusion-at-index operand index)))
    
    (($ farray 'fusion-expr 'log (operand) _ _ _ _)
     (log (evaluate-fusion-at-index operand index)))
    
    (($ farray 'fusion-expr 'sin (operand) _ _ _ _)
     (sin (evaluate-fusion-at-index operand index)))
    
    (($ farray 'fusion-expr 'cos (operand) _ _ _ _)
     (cos (evaluate-fusion-at-index operand index)))
    
    (($ farray 'fusion-expr 'tan (operand) _ _ _ _)
     (tan (evaluate-fusion-at-index operand index)))
    
    (($ farray 'fusion-expr 'floor (operand) _ _ _ _)
     (floor (evaluate-fusion-at-index operand index)))
    
    (($ farray 'fusion-expr 'ceiling (operand) _ _ _ _)
     (ceiling (evaluate-fusion-at-index operand index)))
    
    ;; Comparison operations
    (($ farray 'fusion-expr 'gt (operand1 operand2) _ _ _ _)
     (if (> (evaluate-fusion-at-index operand1 index)
            (evaluate-fusion-at-index operand2 index))
         1.0 0.0))
    
    (($ farray 'fusion-expr 'lt (operand1 operand2) _ _ _ _)
     (if (< (evaluate-fusion-at-index operand1 index)
            (evaluate-fusion-at-index operand2 index))
         1.0 0.0))
    
    (($ farray 'fusion-expr 'eq (operand1 operand2) _ _ _ _)
     (if (= (evaluate-fusion-at-index operand1 index)
            (evaluate-fusion-at-index operand2 index))
         1.0 0.0))

    ;; Logical operations (0.0 = false, anything else = true)
    (($ farray 'fusion-expr 'and (operand1 operand2) _ _ _ _)
     (if (and (not (= (evaluate-fusion-at-index operand1 index) 0.0))
              (not (= (evaluate-fusion-at-index operand2 index) 0.0)))
         1.0 0.0))

    (($ farray 'fusion-expr 'or (operand1 operand2) _ _ _ _)
     (if (or (not (= (evaluate-fusion-at-index operand1 index) 0.0))
             (not (= (evaluate-fusion-at-index operand2 index) 0.0)))
         1.0 0.0))

    (($ farray 'fusion-expr 'not (operand) _ _ _ _)
     (if (= (evaluate-fusion-at-index operand index) 0.0) 1.0 0.0))

    ;; Element-wise min/max (distinct from reduction min/max)
    (($ farray 'fusion-expr 'min (operand1 operand2) _ _ _ _)
     (min (evaluate-fusion-at-index operand1 index)
          (evaluate-fusion-at-index operand2 index)))

    (($ farray 'fusion-expr 'max (operand1 operand2) _ _ _ _)
     (max (evaluate-fusion-at-index operand1 index)
          (evaluate-fusion-at-index operand2 index)))
    
    ;; Map operation - apply function to operand at index
    (($ farray 'fusion-expr 'map (operand) (func) _ _ _)
     (func (evaluate-fusion-at-index operand index)))
    
    ;; Enhanced slice with step support
    (($ farray 'fusion-expr 'slice (source) (start end step) _ _ _)
     (let ((source-index (+ start (* index step))))
       (if (and (>= source-index start) (< source-index end))
           (evaluate-fusion-at-index source source-index)
           0.0)))  ; or handle out-of-bounds as needed
    
    ;; Concatenation - map index to appropriate source array
    (($ farray 'fusion-expr 'concat operands sizes _ _ _)
     (concat-evaluate-at-index operands sizes index))

    ;; Filter operation
    ;; TODO: needs special handling for sparse results
    (($ farray 'fusion-expr 'filter (operand) (predicate) _ _ _)
     (let ((source-val (evaluate-fusion-at-index operand index)))
       (if (predicate source-val)
           source-val
           0.0)))  ; filtered out values become 0 (simplified)
    
    ;; Zip operation - combine multiple arrays element-wise with function
    (($ farray 'fusion-expr 'zip operands (func) _ _ _)
     (apply func (map (lambda (op) (evaluate-fusion-at-index op index)) operands)))
    
    ;; Difference operation
    (($ farray 'fusion-expr 'diff (operand) (n) _ _ _)
     (if (< index n)
         0.0  ; pad with zeros
         (- (evaluate-fusion-at-index operand index)
            (evaluate-fusion-at-index operand (- index n)))))
    
    ;; Moving window operations - these require special handling
    (($ farray 'fusion-expr 'moving-average (source) (window-size) _ _ _)
     (moving-average-at-index source index window-size))
    
    (($ farray 'fusion-expr 'moving-std (source) (window-size) _ _ _)
     (moving-std-at-index source index window-size))
    
    ;; Convolution - evaluated using direct convolution at each point
    (($ farray 'fusion-expr 'convolve (signal kernel) (mode) _ _ _)
     (convolve-at-index signal kernel index mode))
    
    (else
     (error "Unknown fusion expression" expr))))

;;; Helper function for concatenation evaluation
(define (concat-evaluate-at-index operands sizes index)
  "Evaluate concatenation at specific index by mapping to correct source array"
  (let loop ((ops operands) (szs sizes) (offset 0))
    (cond
      ((null? ops) 0.0)  ; out of bounds
      ((< index (+ offset (car szs)))
       (evaluate-fusion-at-index (car ops) (- index offset)))
      (else
       (loop (cdr ops) (cdr szs) (+ offset (car szs)))))))

;;; Type-specific fusion materialization
(define (materialize-fusion expr target-type)
  "Materialize fusion expression with type-specific optimizations"
  (let* ((size (array-size expr))
         (result-data (allocate-typed-vector target-type size)))
    
    ;; Generate optimized loop based on target type
    (case target-type
      ((f64)
       (do ((i 0 (+ i 1)))
           ((= i size))
         (f64vector-set! result-data i 
                        (exact->inexact (evaluate-fusion-at-index expr i)))))
      
      ((f32)
       (do ((i 0 (+ i 1)))
           ((= i size))
         (f32vector-set! result-data i 
                        (exact->inexact (evaluate-fusion-at-index expr i)))))
      
      ((s64)
       (do ((i 0 (+ i 1)))
           ((= i size))
         (s64vector-set! result-data i 
                        (inexact->exact (evaluate-fusion-at-index expr i)))))
      
      ((s32)
       (do ((i 0 (+ i 1)))
           ((= i size))
         (s32vector-set! result-data i 
                        (inexact->exact (evaluate-fusion-at-index expr i)))))
      
      ((u32)
       (do ((i 0 (+ i 1)))
           ((= i size))
         (u32vector-set! result-data i 
                        (inexact->exact (evaluate-fusion-at-index expr i)))))
      
      (else
       (do ((i 0 (+ i 1)))
           ((= i size))
         (vector-set! result-data i (evaluate-fusion-at-index expr i)))))
    
    (materialized-array result-data size target-type 
                       (fusion-expr-metadata expr))))

;;; Reduction materialization
(define (materialize-reduction reduction target-type)
  "Materialize reduction to scalar value"
  (match reduction
    (($ farray 'fusion-reduction op operand params result-type _)
     (case op
       ((sum)
        (let ((size (array-size operand)))
          (do ((i 0 (+ i 1))
               (acc 0.0 (+ acc (evaluate-fusion-at-index operand i))))
              ((= i size) acc))))
       
       ((mean)
        (let ((sum-val (materialize-reduction 
                        (fusion-reduction 'sum operand '() result-type '()) 
                        target-type))
              (size (array-size operand)))
          (/ sum-val size)))
       
       ((min)
        (let ((size (array-size operand)))
          (if (= size 0)
              +inf.0
              (do ((i 1 (+ i 1))
                   (acc (evaluate-fusion-at-index operand 0)
                        (min acc (evaluate-fusion-at-index operand i))))
                  ((= i size) acc)))))
       
       ((max)
        (let ((size (array-size operand)))
          (if (= size 0)
              -inf.0
              (do ((i 1 (+ i 1))
                   (acc (evaluate-fusion-at-index operand 0)
                        (max acc (evaluate-fusion-at-index operand i))))
                  ((= i size) acc)))))
       
       ((var)
        (let* ((mean-val (materialize-reduction 
                          (fusion-reduction 'mean operand '() result-type '()) 
                          target-type))
               (size (array-size operand))
               (ddof (if (null? params) 1 (car params))))
          (when (<= size ddof)
            (error "Not enough data points for variance"))
          (let ((sum-sq-diff
                 (do ((i 0 (+ i 1))
                      (acc 0.0 (let ((diff (- (evaluate-fusion-at-index operand i) mean-val)))
                                 (+ acc (* diff diff)))))
                     ((= i size) acc))))
            (/ sum-sq-diff (- size ddof)))))
       
       ((std)
        (sqrt (materialize-reduction 
               (fusion-reduction 'var operand params result-type '()) 
               target-type)))
       
       ((fold)
        (let ((func (car params))
              (init (cadr params))
              (size (array-size operand)))
          (do ((i 0 (+ i 1))
               (acc init (func (evaluate-fusion-at-index operand i) acc)))
              ((= i size) acc))))
       
       (else
        (error "Unknown reduction operation" op))))))

;;;=============================================================================
;;; Array construction and basic operations
;;;=============================================================================

(define (make-array type size #!optional (fill-value 0))
  "Create a new materialized array"
  (let ((data (allocate-typed-vector type size fill-value)))
    (materialized-array data size type '())))

(define (array-from-list lst #!optional (target-type #f))
  "Create array from list with automatic type inference"
  (let* ((size (length lst))
         (inferred-type (or target-type (infer-list-type lst)))
         (data (allocate-typed-vector inferred-type size)))
    
    ;; Fill the array
    (do ((i 0 (+ i 1))
         (lst lst (cdr lst)))
        ((null? lst))
      (typed-vector-set! data i (car lst) inferred-type))
    
    (materialized-array data size inferred-type '())))

(define (array-size arr)
  "Get size of array (works for materialized, fusion expressions, and reductions)"
  (match arr
    (($ farray 'materialized-array _ size _ _) size)
    (($ farray 'fusion-expr _ _ _ _ size _) size)
    (($ farray 'fusion-reduction _ _ _ _ _) 1)
    (else (error 'array-size "no matching pattern for array" arr)))
  )  ; reductions produce scalars

(define (array-ref arr index)
  "Get element at index (forces evaluation for fusion expressions)"
  (match arr
    (($ farray 'materialized-array data _ element-type _)
     (typed-vector-ref data index element-type))
    (($ farray 'fusion-expr _ _ _ _ _ _)
     (evaluate-fusion-at-index arr index))
    (($ farray 'fusion-reduction _ _ _ _ _)
     (compute! arr))  ; reductions are scalars
    (else (error 'array-ref "no matching pattern for array" arr)))
  )
(define (array->list arr)
  "Convert array to list (forces full materialization)"
  (match arr
    (($ farray 'fusion-reduction _ _ _ _ _)
     (list (compute! arr)))  ; reductions produce single values
    (else
     (let ((size (array-size arr)))
       (do ((i (- size 1) (- i 1))
            (acc '() (cons (array-ref arr i) acc)))
           ((< i 0) acc))))))

;;; Core materialization functions
(define (compute! expr #!optional (target-type #f))
  "Force computation of fusion expression"
  (match expr
    ;; Already materialized
    (($ farray 'materialized-array _ _ _ _) 
     (if target-type
         (convert-array-type expr target-type)
         expr))
    
    ;; Fusion expression - materialize with true fusion
    (($ farray 'fusion-expr _ _ _ result-type _ _)
     (materialize-fusion expr (or target-type result-type)))
    
    ;; Reduction expression - compute scalar result
    (($ farray 'fusion-reduction _ _ _ result-type _)
     (materialize-reduction expr (or target-type result-type)))
    (else (error 'compute! "no matching pattern for array" expr))))

(define (force! expr)
  "Alias for compute!"
  (compute! expr))

;;;=============================================================================
;;; Lazy array operations. All create fusion expressions.
;;;=============================================================================

;;; Binary arithmetic operations
(define (array+ arr1 arr2)
  "Lazy element-wise addition"
  (create-binary-fusion 'add arr1 arr2))

(define (array- arr1 arr2)
  "Lazy element-wise subtraction"
  (create-binary-fusion 'sub arr1 arr2))

(define (array* arr1 arr2)
  "Lazy element-wise multiplication"
  (create-binary-fusion 'mul arr1 arr2))

(define (array/ arr1 arr2)
  "Lazy element-wise division"
  (create-binary-fusion 'div arr1 arr2))

(define (array-pow arr1 arr2)
  "Lazy element-wise exponentiation"
  (create-binary-fusion 'pow arr1 arr2))

(define (array-min arr1 arr2)
  "Lazy element-wise minimum"
  (create-binary-fusion 'min arr1 arr2))

(define (array-max arr1 arr2)
  "Lazy element-wise maximum"
  (create-binary-fusion 'max arr1 arr2))

;;; Scalar operations
(define (array-scale arr scalar)
  "Lazy scalar multiplication"
  (create-scalar-fusion 'scale arr scalar))

(define (array-add-scalar arr scalar)
  "Lazy scalar addition"
  (create-scalar-fusion 'add-scalar arr scalar))

(define (array-sub-scalar arr scalar)
  "Lazy scalar subtraction"
  (create-scalar-fusion 'sub-scalar arr scalar))

(define (array-div-scalar arr scalar)
  "Lazy scalar division"
  (create-scalar-fusion 'div-scalar arr scalar))

;;; Unary operations
(define (array-negate arr)
  "Lazy element-wise negation"
  (create-unary-fusion 'negate arr))

(define (array-abs arr)
  "Lazy element-wise absolute value"
  (create-unary-fusion 'abs arr))

(define (array-sqrt arr)
  "Lazy element-wise square root"
  (create-unary-fusion 'sqrt arr))

(define (array-exp arr)
  "Lazy element-wise exponential"
  (create-unary-fusion 'exp arr))

(define (array-log arr)
  "Lazy element-wise natural logarithm"
  (create-unary-fusion 'log arr))

(define (array-sin arr)
  "Lazy element-wise sine"
  (create-unary-fusion 'sin arr))

(define (array-cos arr)
  "Lazy element-wise cosine"
  (create-unary-fusion 'cos arr))

(define (array-tan arr)
  "Lazy element-wise tangent"
  (create-unary-fusion 'tan arr))

(define (array-floor arr)
  "Lazy element-wise floor"
  (create-unary-fusion 'floor arr))

(define (array-ceiling arr)
  "Lazy element-wise ceiling"
  (create-unary-fusion 'ceiling arr))

;;; Comparison operations
(define (array> arr1 arr2)
  "Lazy element-wise greater than"
  (create-binary-fusion 'gt arr1 arr2))

(define (array< arr1 arr2)
  "Lazy element-wise less than"
  (create-binary-fusion 'lt arr1 arr2))

(define (array= arr1 arr2)
  "Lazy element-wise equality"
  (create-binary-fusion 'eq arr1 arr2))

(define (array-and arr1 arr2)
  "Lazy element-wise logical AND (0.0 = false, non-zero = true)"
  (create-binary-fusion 'and arr1 arr2))

(define (array-or arr1 arr2)
  "Lazy element-wise logical OR (0.0 = false, non-zero = true)"
  (create-binary-fusion 'or arr1 arr2))

(define (array-not arr)
  "Lazy element-wise logical NOT (0.0 → 1.0, non-zero → 0.0)"
  (create-unary-fusion 'not arr))

;;;=============================================================================
;;; Core functional operations
;;;=============================================================================

;;; Map operation
(define (array-map func arr)
  "Lazy map operation - applies function to each element"
  (let* ((size (array-size arr))
         (arr-type (get-array-type arr))
         ;; Try to infer result type by applying function to a sample
         (sample-result (func (array-ref arr 0)))
         (result-type (classify-number sample-result)))
    (fusion-expr 'map (list arr) (list func)
                 (or result-type 'generic) size '())))

;;; Reduce operations
(define (array-reduce op arr #!optional (params '()))
  "Create lazy reduction operation"
  (let* ((arr-type (get-array-type arr))
         (result-type (case op
                        ((sum mean) arr-type)
                        ((var std) (case arr-type
                                     ((s32 s64 u32) 'f64)
                                     (else arr-type)))
                        ((min max) arr-type)
                        ((fold) 'generic)
                        (else 'generic))))
    (fusion-reduction op arr params result-type '())))

(define (array-fold func init arr)
  "Lazy fold operation with function and initial value"
  (array-reduce 'fold arr (list func init)))

;;; Slice operation
(define (array-slice arr start end #!optional (step 1))
  "Create enhanced lazy slice with step support"
  (let* ((arr-size (array-size arr))
         (clamped-start (max 0 (min start arr-size)))
         (clamped-end (max clamped-start (min end arr-size)))
         (result-size (max 0 (inexact->exact (ceiling (/ (- clamped-end clamped-start) step)))))
         (result-type (get-array-type arr)))
    (when (<= step 0)
      (error "Step must be positive" step))
    (fusion-expr 'slice (list arr) (list clamped-start clamped-end step)
                 result-type result-size '())))

;;; Lazy concatenation
(define (array-concat . arrays)
  "Create lazy concatenation of multiple arrays"
  (when (null? arrays)
    (error "Cannot concatenate empty list of arrays"))
  (let* ((sizes (map array-size arrays))
         (total-size (apply + sizes))
         (result-type (fold promote-types (get-array-type (car arrays))
                           (map get-array-type (cdr arrays)))))
    (fusion-expr 'concat arrays sizes result-type total-size '())))

;;; Filter operation (creates sparse representation)
(define (array-filter predicate arr)
  "Lazy filter operation - creates sparse array representation"
  (let* ((size (array-size arr))
         (arr-type (get-array-type arr)))
    ;; Note: This creates a sparse representation where filtered elements are 0
    ;; For true filtering, we'd need to compute the actual result size
    (fusion-expr 'filter (list arr) (list predicate) arr-type size '())))

;;; Zip operation
(define (array-zip func . arrays)
  "Lazy zip operation - combines multiple arrays with function"
  (when (null? arrays)
    (error "Cannot zip empty list of arrays"))
  (let* ((sizes (map array-size arrays))
         (min-size (apply min sizes))
         (types (map get-array-type arrays))
         (result-type (fold promote-types (car types) (cdr types))))
    ;; Check that all arrays have the same size
    (unless (apply = sizes)
      (print "Warning: Arrays have different sizes, using minimum size"))
    (fusion-expr 'zip arrays (list func) result-type min-size '())))

;;;=============================================================================
;;; Pipeline-compatible (data-first) variants of data-last operations
;;;=============================================================================

(define (array-map* arr func)
  "Data-first variant of array-map for use in array-pipeline.
   (array-map* arr func) == (array-map func arr)"
  (array-map func arr))

(define (array-fold* arr func init)
  "Data-first variant of array-fold for use in array-pipeline.
   (array-fold* arr func init) == (array-fold func init arr)"
  (array-fold func init arr))

(define (array-reduce* arr op params)
  "Data-first, non-variadic variant of array-reduce for use in array-pipeline.
   (array-reduce* arr op params) == (array-reduce op arr params)"
  (array-reduce op arr params))

(define (array-filter* arr predicate)
  "Data-first variant of array-filter for use in array-pipeline.
   (array-filter* arr predicate) == (array-filter predicate arr)"
  (array-filter predicate arr))

(define (array-zip* arr func . other-arrays)
  "Data-first variant of array-zip for use in array-pipeline.
   (array-zip* arr func arr2 arr3 ...) == (array-zip func arr arr2 arr3 ...)
   In a pipeline, arr is supplied by the macro; func and other-arrays are given explicitly."
  (apply array-zip func arr other-arrays))

;;;=============================================================================
;;; Statistical operations
;;;=============================================================================

(define (array-sum arr)
  "Compute sum using lazy reduction"
  (compute! (array-reduce 'sum arr)))

(define (array-mean arr)
  "Compute mean using lazy reduction"
  (compute! (array-reduce 'mean arr)))

(define (array-var arr #!optional (ddof 1))
  "Compute variance using lazy reduction"
  (compute! (array-reduce 'var arr (list ddof))))

(define (array-std arr #!optional (ddof 1))
  "Compute standard deviation using lazy reduction"
  (compute! (array-reduce 'std arr (list ddof))))

(define (array-min-val arr)
  "Find minimum value using lazy reduction"
  (compute! (array-reduce 'min arr)))

(define (array-max-val arr)
  "Find maximum value using lazy reduction"
  (compute! (array-reduce 'max arr)))

;;; Moving window operations (lazy)
(define (array-moving-average arr window-size)
  "Lazy moving average"
  (let* ((size (array-size arr))
         (result-size size)
         (result-type (get-array-type arr)))
    (fusion-expr 'moving-average (list arr) (list window-size)
                 result-type result-size '())))

(define (array-moving-std arr window-size)
  "Lazy moving standard deviation"
  (let* ((size (array-size arr))
         (result-size size)
         (result-type (promote-types (get-array-type arr) 'f64)))
    (fusion-expr 'moving-std (list arr) (list window-size)
                 result-type result-size '())))

(define (array-moving-min arr window-size)
  "Lazy moving minimum"
  (let* ((size (array-size arr))
         (result-size (max 0 (+ (- size window-size) 1)))
         (result-type (get-array-type arr)))
    (fusion-expr 'moving-min (list arr) (list window-size)
                 result-type result-size '())))

(define (array-moving-max arr window-size)
  "Lazy moving maximum"
  (let* ((size (array-size arr))
         (result-size (max 0 (+ (- size window-size) 1)))
         (result-type (get-array-type arr)))
    (fusion-expr 'moving-max (list arr) (list window-size)
                 result-type result-size '())))

;;;=============================================================================
;;; Additional functional operations
;;;=============================================================================

(define (array-take arr n)
  "Take first n elements"
  (array-slice arr 0 (min n (array-size arr))))

(define (array-drop arr n)
  "Drop first n elements"
  (array-slice arr (min n (array-size arr)) (array-size arr)))

(define (array-diff arr #!optional (n 1))
  "Discrete difference using fusion"
  (let* ((size (array-size arr))
         (result-size (max 0 (- size n)))
         (result-type (get-array-type arr)))
    (fusion-expr 'diff (list arr) (list n) result-type result-size '())))

(define (array-every? predicate arr)
  "Test if predicate is true for all elements"
  (let ((filtered (array-filter predicate arr)))
    (= (array-size filtered) (array-size arr))))

(define (array-any? predicate arr)
  "Test if predicate is true for any element"
  (let ((filtered (array-filter predicate arr)))
    (> (array-size filtered) 0)))

(define (array-count predicate arr)
  "Count elements satisfying predicate"
  (compute! (array-reduce 'sum (array-map (lambda (x) (if (predicate x) 1 0)) arr))))

(define (array-select arr indices)
  "Select elements at given indices (lazy)"
  (array-map (lambda (idx) (array-ref arr (inexact->exact idx))) 
             (if (farray? indices) indices (array-from-list indices))))

(define (array-where predicate arr)
  "Return indices where predicate is true"
  (let ((size (array-size arr)))
    (let loop ((i 0) (indices '()))
      (if (= i size)
          (array-from-list (reverse indices))
          (loop (+ i 1)
                (if (array-ref arr i)
                    (cons i indices)
                    indices))))
    ))

(define (array-unique arr)
  "Return unique elements (forces materialization for hash table)"
  (let ((seen (make-hash-table))
        (size (array-size arr)))
    (let loop ((i 0) (unique-list '()))
      (if (= i size)
          (array-from-list (reverse unique-list))      
          (let* ((val (array-ref arr i))
                 (seen? (hash-table-exists? seen val)))
            (if (not seen?) (hash-table-set! seen val #t))
            (loop (+ i 1) (if (not seen)
                              (cons val unique-list)
                              unique-list))))
      ))
  )

(define (array-group-by key-func arr)
  "Group array elements by key function (forces materialization)"
  (let ((groups (make-hash-table))
        (size (array-size arr)))
    (do ((i 0 (+ i 1)))
        ((= i size))
      (let* ((val (array-ref arr i))
             (key (key-func val)))
        (hash-table-update!/default groups key
                                   (lambda (group) (cons val group))
                                   '())))
    ;; Convert to association list of arrays
    (map (lambda (key-group)
           (cons (car key-group) 
                 (array-from-list (reverse (cdr key-group)))))
         (hash-table->alist groups))))

(define (array-cumsum arr)
  "Cumulative sum using fold with accumulation"
  (let* ((size (array-size arr))
         (result-data (allocate-typed-vector (get-array-type arr) size)))
    (let loop ((i 0) (running-sum 0.0))
      (if (= i size)
          (materialized-array result-data size (get-array-type arr) '())
          (begin
            (typed-vector-set! result-data i running-sum (get-array-type arr))
            (loop (+ i 1) (+ running-sum (array-ref arr i))))
        ))
    ))

;;;=============================================================================
;;; Signal processing
;;;=============================================================================


  #>
  #include "fft_serial.h"
  <#

;; Foreign function bindings (C implementation from fft_serial.c)
(define cffti
  (foreign-lambda void "cffti"
                  int              ; n         array length (power of 2)
                  (c-pointer double))) ; w[]   twiddle-factor table (length n)

(define cfft2
  (foreign-lambda void "cfft2"
                  int              ; n
                  (c-pointer double) ; x[]    input  (interleaved complex, 2n doubles)
                  (c-pointer double) ; y[]    output (interleaved complex, 2n doubles)
                  (c-pointer double) ; w[]    twiddle-factor table
                  double))           ; sgn    +1.0 forward, -1.0 inverse

(define (power-of-2? n)
  "True iff n is a positive power of two."
  (and (> n 0) (= 0 (bitwise-and n (- n 1)))))

(define (fill-complex-buffer! buf arr n)
  "Fill interleaved-complex f64vector BUF from real fusion array ARR.
   BUF must have length 2*n.  Imaginary parts are initialised to 0."
  (do ((i 0 (+ i 1)))
      ((= i n))
    (f64vector-set! buf (* i 2) (exact->inexact (array-ref arr i)))))

(define (extract-complex-result buf n scale)
  "Split interleaved-complex f64vector BUF into two materialized f64 arrays,
   multiplying every value by SCALE.  Returns (values real-arr imag-arr)."
  (let ((re (make-f64vector n))
        (im (make-f64vector n)))
    (do ((i 0 (+ i 1)))
        ((= i n))
      (f64vector-set! re i (* scale (f64vector-ref buf (* i 2))))
      (f64vector-set! im i (* scale (f64vector-ref buf (+ (* i 2) 1)))))
    (values (materialized-array re n 'f64 '())
            (materialized-array im n 'f64 '()))))

(define (array-fft arr)
  "Forward FFT of a real-valued fusion array.
   Returns (values real-part imag-part) as materialized f64 arrays.
   The input is treated as purely real; its size must be a power of 2.

   Complexity: O(n log n).  Forces materialization of ARR."
  (let* ((n (array-size arr)))
    (unless (power-of-2? n)
      (error "array-fft: input size must be a power of 2" n))
    (let ((w (make-f64vector n   0.0))   ; twiddle factors  - length n
          (x (make-f64vector (* 2 n) 0.0))   ; interleaved input  - 2*n
          (y (make-f64vector (* 2 n) 0.0)))  ; interleaved output - 2*n
      (fill-complex-buffer! x arr n)
      (cffti n (make-locative w))
      ;; cfft2 always writes its result into the Y buffer
      (cfft2 n (make-locative x) (make-locative y) (make-locative w) +1.0)
      (extract-complex-result y n 1.0))))

(define (array-ifft real-arr imag-arr)
  "Inverse FFT from separate real and imaginary part arrays.
   Returns (values real-part imag-part) scaled by 1/n.
   Both input arrays must have equal size, which must be a power of 2.

   For the IFFT of a real signal's FFT the imaginary output will be ~0
   (up to floating-point rounding).  Forces materialization of both inputs."
  (let* ((n (array-size real-arr)))
    (unless (power-of-2? n)
      (error "array-ifft: input size must be a power of 2" n))
    (unless (= n (array-size imag-arr))
      (error "array-ifft: real and imaginary arrays must have equal size"
             n (array-size imag-arr)))
    (let ((w (make-f64vector n      0.0))
          (x (make-f64vector (* 2 n) 0.0))
          (y (make-f64vector (* 2 n) 0.0)))
      ;; Interleave real and imaginary inputs into x
      (do ((i 0 (+ i 1)))
          ((= i n))
        (f64vector-set! x (* i 2)       (exact->inexact (array-ref real-arr i)))
        (f64vector-set! x (+ (* i 2) 1) (exact->inexact (array-ref imag-arr i))))
      (cffti n (make-locative w))
      ;; sgn = -1.0 selects the inverse transform; C code does NOT divide by n
      (cfft2 n (make-locative x) (make-locative y) (make-locative w) -1.0)
      ;; Normalise: divide every output element by n
      (extract-complex-result y n (/ 1.0 n)))))

(define (array-convolve signal kernel #!optional (mode 'valid))
  "Lazy convolution operation"
  (let* ((sig-size (array-size signal))
         (ker-size (array-size kernel))
         (sig-type (get-array-type signal))
         (ker-type (get-array-type kernel))
         (result-type (promote-types sig-type ker-type))
         (result-size (case mode
                        ((valid) (max 0 (+ (- sig-size ker-size) 1)))
                        ((same) sig-size)
                        ((full) (+ sig-size ker-size -1)))))
    (fusion-expr 'convolve (list signal kernel) (list mode)
                 result-type result-size '())))

(define (array-gaussian-filter arr sigma)
  "Apply Gaussian filter using lazy convolution"
  (let* ((kernel-size (min 31 (round (inexact->exact (+ (* 6 sigma) 0.5)))))
         (kernel-half (quotient kernel-size 2))
         (kernel-data (make-f64vector kernel-size))
         (kernel-sum
          (let loop ((i 0) (sum 0.0))
            (if (= i kernel-size) sum
                (let* ((x (- i kernel-half))
                       (val (exp (- (/ (* x x) (* 2 sigma sigma))))))
                  (f64vector-set! kernel-data i val)
                  (loop (+ i 1)
                        (+ sum val))))))
         )
                
    
    ;; Create normalized Gaussian kernel
    (do ((i 0 (+ i 1)))
        ((= i kernel-size))
      (f64vector-set! kernel-data i (/ (f64vector-ref kernel-data i) kernel-sum)))
    
    (let ((kernel (materialized-array kernel-data kernel-size 'f64 '())))
      (array-convolve arr kernel 'same))))

(define (array-correlate arr1 arr2 #!optional (mode 'full))
  "Cross-correlation of two arrays"
  (case mode
    ((full)
     ;; Full correlation: convolve arr1 with reversed arr2
     (let ((rev-arr2 (array-slice arr2 (- (array-size arr2) 1) -1 -1)))
       (array-convolve arr1 rev-arr2 'full)))
    ((valid)
     (array-convolve arr1 arr2 'valid))
    ((same)
     (array-convolve arr1 arr2 'same))
    (else
     (error "Unknown correlation mode" mode))))

(define (array-autocorrelation arr)
  "Autocorrelation of array"
  (array-correlate arr arr 'full))

(define (array-integrate arr dx)
  "Numerical integration using trapezoidal rule"
  (let* ((size (array-size arr))
         (integral (array-scale (array-cumsum arr) dx)))
    ;; Apply trapezoidal correction
    (array+ (array-scale (array-slice integral 0 (- size 1)) 0.5)
            (array-scale (array-slice integral 1 size) 0.5))))

;;;=============================================================================
;;; Utility functions
;;;=============================================================================

  
(define (array-linspace start end n)
  "Create linearly spaced array"
  (when (<= n 1)
    (error "Need at least 2 points for linspace"))
  (let* ((step (/ (- end start) (- n 1)))
         (data (make-f64vector n)))
    (do ((i 0 (+ i 1)))
        ((= i n))
      (f64vector-set! data i (+ start (* i step))))
    (materialized-array data n 'f64 '())))

(define (array-meshgrid x-arr y-arr)
  "Create coordinate grids for 2D evaluation.
   Returns (values XX YY), each a flat f64 array of length nx*ny
   in row-major order.  Analogous to NumPy meshgrid with indexing='xy'."
  (let* ((nx (array-size x-arr))
         (ny (array-size y-arr))
         (total (* nx ny))
         (xx (make-f64vector total))
         (yy (make-f64vector total)))
    (do ((i 0 (+ i 1)))           ; row  -> y dimension
        ((= i ny))
      (do ((j 0 (+ j 1)))         ; col  -> x dimension
          ((= j nx))
        (let ((idx (+ (* i nx) j)))
          (f64vector-set! xx idx (exact->inexact (array-ref x-arr j)))
          (f64vector-set! yy idx (exact->inexact (array-ref y-arr i))))))
    (values
     (materialized-array xx total 'f64 '())
     (materialized-array yy total 'f64 '()))))

(define random-state (make-parameter (random:init)))
(define (set-random-seed! x) (random-state (random:init x)))
(define (random-uniform)
  (random:randu! (random-state)))

(define (array-random-normal n #!optional (mean 0.0) (std 1.0) (seed #f))
  "Generate normally distributed random array"
  (when seed (set-random-seed! seed))
  (let ((data (make-f64vector n)))
    (do ((i 0 (+ i 2)))
        ((>= i n))
      ;; Box-Muller transform
      (let* ((u1 (random-uniform))
             (u2 (random-uniform))
             (z0 (* std (sqrt (* -2.0 (log u1))) (cos (* 2.0 3.14159 u2))))
             (z1 (* std (sqrt (* -2.0 (log u1))) (sin (* 2.0 3.14159 u2)))))
        (f64vector-set! data i (+ mean z0))
        (when (< (+ i 1) n)
          (f64vector-set! data (+ i 1) (+ mean z1)))))
    (materialized-array data n 'f64 '())))

(define (array-zeros n #!optional (type 'f64))
  "Create array of zeros"
  (make-array type n 0))

(define (array-ones n #!optional (type 'f64))
  "Create array of ones"
  (make-array type n 1))

(define (array-constant n value #!optional (type 'f64))
  "Create array filled with constant value"
  (make-array type n value))

(define (array-range n #!optional (start 0) (step 1) (type 'f64))
  "Create range array"
  (let ((data (allocate-typed-vector type n)))
    (do ((i 0 (+ i 1)))
        ((= i n))
      (typed-vector-set! data i (+ start (* i step)) type))
    (materialized-array data n type '())))

;;; Pipeline macro for readable chaining

(define-syntax array-pipeline
  (syntax-rules (lambda)
    ((_ arr) arr)
    ;; Handle the case where op is a lambda expression
    ((_ arr ((lambda args body) . more) rest ...)
     (array-pipeline ((lambda args body) arr . more) rest ...))
    ((_ arr (lambda args body) rest ...)
     (array-pipeline ((lambda args body) arr) rest ...))
    ((_ arr (op . args) rest ...)
     (array-pipeline (op arr . args) rest ...))
    ))

;;;=============================================================================
;;; Helpers routines
;;;=============================================================================


(define (fusion-expr-metadata x)
  (match x
    (($ farray 'fusion-expr operations operands _ _ _ metadata)
     metadata)
    (else (error 'fusion-expr-metdata "argument is not a fusion expression" x))))

(define (create-binary-fusion op arr1 arr2)
  "Create binary fusion expression"
  (let* ((size1 (array-size arr1))
         (size2 (array-size arr2))
         (type1 (get-array-type arr1))
         (type2 (get-array-type arr2)))
    (unless (= size1 size2)
      (error "Array sizes must match" size1 size2))
    (let ((result-type (infer-result-type op (list type1 type2))))
      (fusion-expr op (list arr1 arr2) '() result-type size1 '()))))

(define (create-unary-fusion op arr)
  "Create unary fusion expression"
  (let* ((size (array-size arr))
         (type (get-array-type arr))
         (result-type (infer-result-type op (list type))))
    (fusion-expr op (list arr) '() result-type size '())))

(define (create-scalar-fusion op arr scalar)
  "Create scalar fusion expression"
  (let* ((size (array-size arr))
         (type (get-array-type arr))
         (result-type (infer-result-type op (list type))))
    (fusion-expr op (list arr) (list scalar) result-type size '())))

(define (get-array-type arr)
  "Get the element type of an array"
  (match arr
    (($ farray 'materialized-array _ _ element-type _) element-type)
    (($ farray 'fusion-expr _ _ _ result-type _ _) result-type)
    (($ farray 'fusion-reduction _ _ _ result-type _) result-type)
    (else (error 'get-array-type "no matching pattern for array" arr))
    ))

(define (allocate-typed-vector type size #!optional (fill 0))
  "Allocate vector of specific type"
  (case type
    ((f64) (make-f64vector size (exact->inexact fill)))
    ((f32) (make-f32vector size (exact->inexact fill)))
    ((s64) (make-s64vector size (inexact->exact fill)))
    ((s32) (make-s32vector size (inexact->exact fill)))
    ((u32) (make-u32vector size (inexact->exact fill)))
    (else (make-vector size fill))))

(define (typed-vector-ref vec index type)
  "Reference element from typed vector"
  (case type
    ((f64) (f64vector-ref vec index))
    ((f32) (f32vector-ref vec index))
    ((s64) (s64vector-ref vec index))
    ((s32) (s32vector-ref vec index))
    ((u32) (u32vector-ref vec index))
    (else (vector-ref vec index))))

(define (typed-vector-set! vec index value type)
  "Set element in typed vector"
  (case type
    ((f64) (f64vector-set! vec index (exact->inexact value)))
    ((f32) (f32vector-set! vec index (exact->inexact value)))
    ((s64) (s64vector-set! vec index (inexact->exact value)))
    ((s32) (s32vector-set! vec index (inexact->exact value)))
    ((u32) (u32vector-set! vec index (inexact->exact value)))
    (else (vector-set! vec index value))))

(define (infer-list-type lst)
  "Infer appropriate type for list of values"
  (let ((types (filter identity (map classify-number lst))))
    (if (null? types)
        'generic
        (fold promote-types (car types) (cdr types)))))

(define (classify-number val)
  "Classify a number by type"
  (cond
   ((not (number? val)) #f)
   ((complex? val) 'generic)
   ((real? val) 'f64)
   ((and (number? val) (exact? val) (integer? val))
    (cond
     ((and (>= val -2147483648) (<= val 2147483647)) 's32)
     ((and (>= val 0) (<= val 4294967295)) 'u32)
     (else 's64)))
   (else 'generic)))

;;; Special fusion operations that require window access
(define (moving-average-at-index source index window-size)
  "Compute moving average at specific index"
  (let ((start (max 0 (- index (quotient window-size 2))))
        (end (min (array-size source) (+ index (quotient window-size 2) 1))))
    (let loop ((i start) (sum 0.0) (count 0))
      (if (= i end)
          (/ sum count)
          (loop (+ i 1)
                (+ sum (evaluate-fusion-at-index source i))
                (+ count 1))))
    ))

(define (moving-std-at-index source index window-size)
  "Compute moving standard deviation at specific index"
  (let ((start (max 0 (- index (quotient window-size 2))))
        (end (min (array-size source) (+ index (quotient window-size 2) 1))))
    (let loop ((i start) (mean 0.0) (m2 0.0) (count 0))
      (if (= i end)
          (if (> count 1)
              (sqrt (/ m2 (- count 1)))
              0.0)
          (let* ((x (evaluate-fusion-at-index source i))
                 (delta (- x mean)))
            (loop (+ i 1)
                  (+ mean (/ delta count))
                  (+ m2 (* delta (- x mean)))
                  (+ count 1))
            ))
      ))
  )
      

(define (convolve-at-index signal kernel index mode)
  "Compute convolution at specific index"
  (let* ((sig-size (array-size signal))
         (ker-size (array-size kernel)))
    (let loop ((j 0) (sum 0.0))
      (if (= j ker-size) sum
          (let ((sig-idx (case mode
                           ((valid) (+ index j))
                           ((same) (+ index j (- (quotient ker-size 2))))
                           ((full) (+ index j)))))
            (loop (+ j 1)
                  (if (and (>= sig-idx 0) (< sig-idx sig-size))
                      (+ sum (* (evaluate-fusion-at-index signal sig-idx)
                                (evaluate-fusion-at-index kernel j)))
                      sum))
            )))
    ))

;;; Performance monitoring and debugging
(define *fusion-stats* (make-hash-table))

(define (fusion-stats)
  "Get fusion performance statistics"
  (hash-table->alist *fusion-stats*))

(define (clear-fusion-stats!)
  "Clear fusion statistics"
  (hash-table-clear! *fusion-stats*))

(define (array-complexity expr)
  "Calculate computational complexity of expression"
  (match expr
    (($ farray 'materialized-array _ _ _ _) 1)
    (($ farray 'fusion-expr _ operands _ _ _ _)
     (+ 1 (apply + (map array-complexity operands))))
    (($ farray 'fusion-reduction _ operand _ _ _)
     (+ 1 (array-complexity operand)))
    (else (error 'array-complexity "no matching pattern for array" expr))
    ))


(define (fusion-plan expr)
  "Show the fusion plan for an expression"
  (match expr
    (($ farray 'materialized-array _ size type _)
     `(materialized ,type ,size))
    (($ farray 'fusion-expr op operands params _ size _)
     `(fusion ,op ,(map fusion-plan operands) ,params ,size))
    (($ farray 'fusion-reduction op operand params result-type _)
     `(reduction ,op ,(fusion-plan operand) ,params ,result-type))
    (else (error 'fusion-plan "no matching pattern for array" expr))
    ))

(define (pp-fusion expr)
  "Pretty print fusion expression tree"
  (define (print-tree tree depth)
    (let ((indent (make-string (* depth 2) #\space)))
      (match tree
        (('materialized type size)
         (print indent "DATA[" type "," size "]"))
        (('fusion op operands params size)
         (print indent "FUSE[" op "," size "]")
         (for-each (lambda (operand) (print-tree operand (+ depth 1))) operands))
        (('reduction op operand params result-type)
         (print indent "REDUCE[" op "→" result-type "]")
         (print-tree operand (+ depth 1))))))
  (print-tree (fusion-plan expr) 0))

;;; Convert array type
(define (convert-array-type arr target-type)
  "Convert array to target type if needed"
  (if (eq? (get-array-type arr) target-type)
      arr
      (let* ((size (array-size arr))
             (new-data (allocate-typed-vector target-type size)))
        (do ((i 0 (+ i 1)))
            ((= i size))
          (typed-vector-set! new-data i (array-ref arr i) target-type))
        (materialized-array new-data size target-type '()))))

) ; end module