;;; Usage Examples for Fusion Arrays

(import fusion-arrays)
(import scheme (chicken base) (chicken time) (chicken gc) srfi-1 srfi-4)

;;;=============================================================================
;;; Example with all operations
;;;=============================================================================

(define (example-fusion-system)
  "Array fusion system examples"
  
  ;; Basic functional operations
  (print "Map, Reduce, Slice, Concat Operations:")
  (let* ((a (array-from-list '(1.0 2.0 3.0 4.0 5.0)))
         (b (array-from-list '(2.0 4.0 6.0 8.0 10.0)))
         
         ;; Map operation - all lazy!
         (squared (array-map (lambda (x) (* x x)) a))
         (doubled (array-map (lambda (x) (* x 2)) b))
         
         ;; Complex expression with map
         (complex-expr (array+ squared doubled))
         
         ;; Reduction operations
         (sum-a (array-reduce 'sum a))
         (mean-b (array-reduce 'mean b))
         
         ;; Slice operations
         (slice-a (array-slice a 1 4))
         (every-other (array-slice b 0 5 2))
         
         ;; Concatenation
         (concat-ab (array-concat a b)))
    
    (print "  Arrays: a=" (array->list a) " b=" (array->list b))
    (print "  Squared a: " (array->list (compute! squared)))
    (print "  Doubled b: " (array->list (compute! doubled)))
    (print "  Complex (a^2 + 2b): " (array->list (compute! complex-expr)))
    (print "  Sum of a: " (compute! sum-a))
    (print "  Mean of b: " (compute! mean-b))
    (print "  Slice a[1:4]: " (array->list (compute! slice-a)))
    (print "  Every other b: " (array->list (compute! every-other)))
    (print "  Concatenated: " (array->list (compute! concat-ab)))
    (print "  Fusion plan for complex expression:")
    (pp-fusion complex-expr)
    (print))
  
  ;; Performance demonstration with large datasets
  (print "Large-Scale Performance:")
  (let* ((n 100000)
         (large-a (array-random-normal (* 2 n)))
         (large-b (array-random-normal n))
         (large-c (array-random-normal n))
         
         ;; Mega-complex operation using all features
         (mega-expr 
          (array-reduce 'mean
                        (array-map (lambda (x) (* x x))
                                   (array+ (array-slice large-a 0 n)
                                           (array* (array-scale large-b 2.0)
                                                   (array-sqrt (array-abs large-c))))))))
    
    (print "  Large arrays: " n " elements each")
    (print "  Operation: mean(map(x^2, slice(a) + (2b * sqrt(|c|))))")
    (print "  Expression complexity: " (array-complexity mega-expr))
    
    (let ((start (current-process-milliseconds))
          (mem-before (memory-statistics)))
      (let ((result (compute! mega-expr)))
        (let ((end (current-process-milliseconds))
              (mem-after (memory-statistics)))
          (print "  Result: " result)
          (print "  Fusion time: " (- end start) " ms")
          (print "  Memory before: " mem-before " after: " mem-after))))
    (print))
  
  ;; Functional programming patterns
  (print "Functional Patterns:")
  (let* ((data (array-range 100 1 1 'f64))  ; 1 to 100
         
         ;; Map-reduce pattern: sum of squares of even numbers
         (even-squares-sum
          (array-reduce 'sum
                       (array-map (lambda (x) 
                                   (if (= (modulo (inexact->exact x) 2) 0)
                                       (* x x)
                                       0))
                                 data)))
         
         ;; Complex transformation pipeline
         (transformed
          (array-pipeline data
                         ;; Filter range
                         (array-slice 10 90)
                         ;; Transform
                         (array-map* (lambda (x) (sin (* x 0.1))))
                         ;; Moving average
                         (array-moving-average 5)
                         ;; Final scaling
                         (array-scale 100.0)))
         

         ;; Multi-array operations
         (combined (array-zip + 
                             (array-slice data 0 50)
                             (array-slice data 50 100))))
    
    (print "  Data: 1 to 100")
    (print "  Even squares sum: " (compute! even-squares-sum))
    (print "  Transformed pipeline size: " (array-size transformed))
    (print "  Combined array size: " (array-size combined))
    (print))
  
  )

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

(define (bandpass-filter signal low-freq high-freq sampling-rate)
  "Bandpass filter using fusion (simplified version)"
  (let* ((nyquist (/ sampling-rate 2.0))
         (low-norm (/ low-freq nyquist))
         (high-norm (/ high-freq nyquist))
         
         ;; Create filter kernels (simplified Gaussian approximation)
         (low-sigma (/ 1.0 (* 2.0 3.14159 low-norm)))
         (high-sigma (/ 1.0 (* 2.0 3.14159 high-norm)))
         
         ;; Apply filters using fusion
         (lowpass (array-gaussian-filter signal low-sigma))
         (highpass (array- signal (array-gaussian-filter signal high-sigma))))
    
    ;; Combine filters (this creates a fused expression)
    (array- lowpass (array- lowpass highpass))))

(define (spectral-power-density signal window-size overlap)
  "Compute power spectral density using overlapping windows"
  (let* ((signal-size (array-size signal))
         (step-size (inexact->exact (* window-size (- 1.0 overlap))))
         (n-windows (inexact->exact (floor (/ (- signal-size window-size) step-size))))
         (power-spectra '()))
    
    ;; Process each window
    (do ((i 0 (+ i 1))
         (start 0 (+ start step-size)))
        ((or (= i n-windows) (>= (+ start window-size) signal-size)))
      
      (let* (;; Extract window (lazy slice)
             (window (array-slice signal start (+ start window-size)))
             
             ;; Apply window function (fused with slice!)
             (windowed (array* window (hann-window window-size)))
             
             ;; Compute magnitude squared (would need FFT for real implementation)
             (power (array* windowed windowed)))
        
        (set! power-spectra (cons (compute! power) power-spectra))))
    
    ;; Average across windows
    (if (null? power-spectra)
        (array-zeros window-size)
        (let ((sum-spectrum (apply array+ power-spectra)))
          (array-scale sum-spectrum (/ 1.0 (length power-spectra)))))))

(define (hann-window n)
  "Create Hann window function"
  (let ((data (make-f64vector n)))
    (do ((i 0 (+ i 1)))
        ((= i n))
      (f64vector-set! data i (* 0.5 (- 1.0 (cos (/ (* 2.0 3.14159 i) (- n 1)))))))
    (materialized-array data n 'f64 '())))

;;;=============================================================================
;;; Performance benchmarks
;;;=============================================================================

(define (benchmark-fusion-performance)
  "Fusion performance benchmarks"
  (print "=== Fusion Performance Benchmarks ===\n")
  
  ;; Simple vs Complex Expressions
  (print "Expression Complexity Impact:")
  (let ((n 1000000))
    (let* ((a (array-random-normal n))
           (b (array-random-normal n))
           (c (array-random-normal n))
           
           ;; Simple expression
           (simple (array+ a b))
           
           ;; Complex expression  
           (complex (array+ (array* a b)
                            (array/ (array-sqrt c) 
                                    (array+ a (array-constant n 1.0))))))
      
      (print "  Array size: " n " elements")
      (print "  Simple expression complexity: " (array-complexity simple))
      (print "  Complex expression complexity: " (array-complexity complex))
      
      ;; Time simple expression
      (let ((start (current-milliseconds)))
        (let ((result1 (compute! simple)))
          (let ((mid (current-milliseconds)))
            ;; Time complex expression  
            (let ((result2 (compute! complex)))
              (let ((end (current-milliseconds)))
                (print "  Simple time: " (- mid start) " ms")
                (print "  Complex time: " (- end mid) " ms")
                ))))
        (print))))
    
  
  ;; Memory Usage Comparison
  (print "Memory Usage Analysis:")
  (let ((n 500000))  ; Half million elements
    (let* ((a (array-random-normal n))
           (b (array-random-normal n))
           (c (array-random-normal n))
           
           ;; Create deeply nested expression
           (nested-expr (array/ (array+ (array* a b) 
                                       (array- c a))
                               (array+ (array-sqrt a) 
                                      (array-exp (array-scale b 0.1))))))
      
      (print "  Expression: ((a*b) + (c-a)) / (sqrt(a) + exp(b*0.1))")
      (print "  Array size: " n " elements (" (* n 8) " bytes each)")
      (print "  Full instantiation approach would create:")
      (print "    temp1 = a * b         (" (* n 8) " bytes)")
      (print "    temp2 = c - a         (" (* n 8) " bytes)")
      (print "    temp3 = temp1 + temp2 (" (* n 8) " bytes)")
      (print "    temp4 = sqrt(a)       (" (* n 8) " bytes)")
      (print "    temp5 = b * 0.1       (" (* n 8) " bytes)")
      (print "    temp6 = exp(temp5)    (" (* n 8) " bytes)")
      (print "    temp7 = temp4 + temp6 (" (* n 8) " bytes)")
      (print "    result = temp3 / temp7(" (* n 8) " bytes)")
      (print "    Total memory: " (* 8 n 8) " bytes (" (exact->inexact (/ (* 8 n 8) 1024 1024)) " MB)") 
      (print)
      (print "  Fusion plan:")
      (pp-fusion nested-expr)
      ;; Actual computation
      (let ((start-mem (vector-ref (memory-statistics) 1))
            (start-time (current-process-milliseconds)))
        (let ((result (compute! nested-expr)))
          (gc)
          (let ((end-mem (vector-ref (memory-statistics) 1))
                (end-time (current-process-milliseconds)))
            (print "  Actual computation:")
            (print "    Time: " (- end-time start-time) " ms")
            (print "    Memory delta: " (- end-mem start-mem) " bytes")
            (print "    Result mean: " (array-mean result)))))
      (print)))
  
  
  ;; Cache Performance
  (print "Cache Performance Benefits:")
  (let ((n 1000000))
    (print "  Large arrays that don't fit in cache: " n " elements")
    (let* ((a (array-range n 0 1 'f64))
           (b (array-range n 1 1 'f64))
           
           ;; Expression that would require multiple cache-missing passes
           (multi-pass-expr (array+ (array+ (array+ a b) 
                                           (array-scale a 2.0))
                                   (array-sqrt b))))
      
      (print "  Expression: ((a + b) + (a * 2)) + sqrt(b)")
      (print "  Fusion plan:")
      (pp-fusion multi-pass-expr)
      (let ((start (current-process-milliseconds)))
        (let ((result (compute! multi-pass-expr)))
          (let ((end (current-process-milliseconds)))
            (print "  Fusion time: " (- end start) " ms")
            (print "  Result size: " (array-size result)))))
      ))
  
)

(define (example-signal-processing raw-eeg channels sampling-rate)
  "Signal preprocessing pipeline with fusion"
  (let* ((n-channels (length channels))
         (n-samples (array-size (car raw-eeg)))
         
         ;; Common average reference (fused across channels)
         (common-average 
          (array-scale (apply array+ raw-eeg) (/ 1.0 n-channels)))
         
         ;; Preprocessed channels
         (preprocessed
          (map (lambda (channel)
                 (array-pipeline channel
                                ;; Remove common average
                                (array- common-average)
                                ;; High-pass filter
                                (bandpass-filter 1.0 100.0 sampling-rate)
                                ;; Artifact removal (simple threshold)
                                (array-min (array-constant n-samples 100.0))
                                (array-max (array-constant n-samples -100.0))))
               raw-eeg)))
    
    `((raw . ,raw-eeg)
      (common-average . ,(compute! common-average))
      (preprocessed . ,(map (lambda (ch) (compute! ch)) preprocessed)))))

;;;=============================================================================
;;; Profiling
;;;=============================================================================

(define (profile-expression expr description)
  "Profile expression with detailed breakdown"
  (print "Profiling: " description)
  (print "Expression plan:")
  (pp-fusion expr)
  (print "Complexity: " (array-complexity expr))
  (print "Size: " (array-size expr) " elements")
  
  (gc)
  (let ((start-mem (vector-ref (memory-statistics) 1))
        (start-time (current-process-milliseconds)))
    (let ((result (compute! expr)))
      (gc)
      (let ((end-mem (vector-ref (memory-statistics) 1))
            (end-time (current-process-milliseconds)))
        (print "Computation time: " (- end-time start-time) " ms")
        (print "Memory delta: " (- end-mem start-mem) " bytes")
        (print "Elements/second: " 
               (/ (array-size expr) 
                  (/ (- end-time start-time) 1000.0)))
        (print "Result statistics:")
        (print "  Mean: " (array-mean result))
        (print "  Std: " (array-std result))
        (print "  Range: [" (array-min-val result) ", " (array-max-val result) "]")
        result))))

(define (compare-fusion-vs-eager expr)
  "Compare fusion vs eager evaluation (simulation)"
  (let* ((complexity (array-complexity expr))
         (size (array-size expr))
         (estimated-intermediates (- complexity 1)))
    
    (print "Fusion vs Eager Comparison:")
    (print "Expression complexity: " complexity)
    (print "Array size: " size " elements")
    (print "Estimated intermediates: " estimated-intermediates)
    (print)
    
    ;; Simulate eager evaluation memory usage
    (print "Eager evaluation (simulated):")
    (print "  Intermediate arrays: " estimated-intermediates)
    (print "  Memory usage: " (* estimated-intermediates size 8) " bytes")
    (print "  Cache passes: " complexity)
    (print)
    
    ;; Actual fusion
    (print "Fusion evaluation (actual):")
    (let ((result (profile-expression expr "Fused computation")))
      (print "  Intermediate arrays: 0")
      (print "  Cache passes: 1")
      result)))

(define (example-advanced-usage)
  "Show advanced usage patterns and performance"
  
  ;; Signal processing example
  (print "Signal Processing Pipeline:")
  (let* ((signal (array+ (array-sin (array-scale (array-range 1000) 0.1))
                        (array-scale (array-random-normal 1000) 0.1)))
         (filtered (bandpass-filter signal 0.05 0.3 10.0)))
    (profile-expression filtered "Bandpass filter")
    (print))
  
  ;; Performance benchmark
  (print "Performance Benchmarks:")
  (benchmark-fusion-performance)
  
  )


(define (example-logical-operations)
  "Demonstrate the logical operations"
  
  ;; Basic logical operations
  (print "Basic Logical Operations:")
  (let* ((a (array-from-list '(0.0 1.0 0.0 1.0 0.0)))
         (b (array-from-list '(0.0 0.0 1.0 1.0 1.0)))
         
         ;; Logical operations
         (and-result (array-and a b))
         (or-result (array-or a b))
         (not-a (array-not a)))
    
    (print "  a = " (array->list a))
    (print "  b = " (array->list b))
    (print "  a AND b = " (array->list (compute! and-result)))
    (print "  a OR b = " (array->list (compute! or-result)))
    (print "  NOT a = " (array->list (compute! not-a)))
    (print))
  
  ;; Complex logical expressions
  (print "Complex Logical Expressions:")
  (let* ((data (array-from-list '(-2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0)))
         
         ;; Create conditions
         (positive (array> data (array-constant (array-size data) 0.0)))
         (less-than-3 (array< data (array-constant (array-size data) 3.0)))
         
         ;; Combine conditions: positive AND less than 3
         (in-range (array-and positive less-than-3))
         
         ;; Negate: outside range
         (out-of-range (array-not in-range)))
    
    (print "  Data: " (array->list data))
    (print "  Positive: " (array->list (compute! positive)))
    (print "  Less than 3: " (array->list (compute! less-than-3)))
    (print "  In range (0, 3): " (array->list (compute! in-range)))
    (print "  Out of range: " (array->list (compute! out-of-range)))
    (print))
  
  ;; Neuroscience application: spike detection with multiple criteria
  (print "Neuroscience Application - Multi-Criteria Spike Detection:")
  (let* ((voltage (array-from-list '(-70.0 -65.0 -30.0 -25.0 -70.0 -35.0 -20.0 -68.0)))
         (derivative (array-from-list '(0.5 5.0 10.0 -5.0 -10.0 5.0 15.0 -2.0)))
         
         ;; Spike criteria
         (threshold -40.0)
         (min-slope 3.0)
         
         ;; Conditions
         (above-threshold (array> voltage (array-constant (array-size voltage) threshold)))
         (steep-rise (array> derivative (array-constant (array-size derivative) min-slope)))
         
         ;; Spike detection: voltage above threshold AND steep rise
         (spike-detected (array-and above-threshold steep-rise)))
    
    (print "  Voltage: " (array->list voltage))
    (print "  Derivative: " (array->list derivative))
    (print "  Above threshold (-40mV): " (array->list (compute! above-threshold)))
    (print "  Steep rise (>3): " (array->list (compute! steep-rise)))
    (print "  Spikes detected: " (array->list (compute! spike-detected)))
    (print))
  
  ;; Combining logical operations in complex pipelines
  (print "Complex Pipeline with Logical Operations:")
  (let* ((signal (array-random-normal 100 0.0 1.0))
         
         ;; Complex filtering: keep values in range [-1, 1] AND not near zero
         (in-main-range (array-and 
                         (array> signal (array-constant (array-size signal) -1.0))
                         (array< signal (array-constant (array-size signal) 1.0))))
         (not-near-zero (array-or
                         (array> signal (array-constant (array-size signal) 0.1))
                         (array< signal (array-constant (array-size signal) -0.1))))
         (keep-mask (array-and in-main-range not-near-zero))
         
         ;; Apply mask (simplified - multiply to zero out unwanted values)
         (filtered (array* signal keep-mask)))
    
    (print "  Signal size: " (array-size signal))
    (print "  Conditions: in [-1,1] AND (> 0.1 OR < -0.1)")
    (print "  Kept elements: " (array-count (lambda (x) (not (= x 0.0))) 
                                            (compute! filtered)))
    (print "  First 10 filtered: " (take (array->list (compute! filtered)) 10))
    (print))
  
  (print "Fusion Verification:")
  (let* ((a (array-from-list '(1.0 2.0 3.0 4.0 5.0)))
         (b (array-from-list '(0.0 1.0 0.0 1.0 0.0)))
         
         ;; Complex fused expression with logical operations
         (complex-expr (array-and 
                        (array> a (array-constant (array-size a) 2.0))
                        (array-or b (array-not b)))))
    
    (print "  Complex logical expression:")
    (print "  (a > 2) AND (b OR (NOT b))")
    (print "  Fusion plan:")
    (pp-fusion complex-expr)
    (print "  Result: " (array->list (compute! complex-expr)))
    (print)))

(define (example-meshgrid)
  "Examples with the meshgrid operation"
  
  ;; Basic 2D meshgrid
  (print "Basic 2D Meshgrid:")
  (let* ((x-coords (array-linspace -2.0 2.0 5))
         (y-coords (array-linspace -1.0 1.0 3)))
    
    (print "  X coordinates: " (array->list x-coords))
    (print "  Y coordinates: " (array->list y-coords))
    
    (receive (x-grid y-grid) (array-meshgrid x-coords y-coords)
      (print "  X grid (flattened): " (array->list x-grid))
      (print "  Y grid (flattened): " (array->list y-grid))
      (print "  Grid dimensions: " (array-size x-coords) " * " (array-size y-coords))
      (print "  Total points: " (array-size x-grid))
      (print)))
  
  ;; Computing function over meshgrid
  (print "Function Evaluation over Grid:")
  (let* ((x-coords (array-linspace -1.0 1.0 11))
         (y-coords (array-linspace -1.0 1.0 11)))
    
    (receive (x-grid y-grid) (array-meshgrid x-coords y-coords)
      ;; Compute z = x^2 + y^2 (distance squared from origin)
      (let ((z-values (array+ (array* x-grid x-grid) 
                             (array* y-grid y-grid))))
        
        (print "  Grid: 11 * 11 points")
        (print "  Function: z = x^2 + y^2")
        (print "  Z values computed: " (array-size z-values))
        (print "  Min z: " (array-min-val (compute! z-values)))
        (print "  Max z: " (array-max-val (compute! z-values)))
        (print "  First 10 z values: " (take (array->list (compute! z-values)) 10))
        (print))))
  
  ;; Neuroscience application: receptive field mapping
  (print "Neuroscience Application - Receptive Field Grid:")
  (let* ((x-positions (array-linspace -10.0 10.0 21))  ; degrees of visual angle
         (y-positions (array-linspace -10.0 10.0 21)))
    
    (receive (x-grid y-grid) (array-meshgrid x-positions y-positions)
      ;; Compute Gaussian receptive field centered at origin
      (let* ((sigma 3.0)
             (response (array-exp 
                        (array-scale 
                         (array+ (array* x-grid x-grid) 
                                (array* y-grid y-grid))
                         (/ -1.0 (* 2.0 sigma sigma))))))
        
        (print "  Receptive field grid: 21 * 21 positions")
        (print "  Gaussian RF with sigma = " sigma "deg")
        (print "  Total stimulus positions: " (array-size response))
        (print "  Peak response: " (array-max-val (compute! response)))
        (print "  Mean response: " (array-mean (compute! response)))
        (print "  Fusion plan:")
        (pp-fusion response)
        (print))))
  
  ;; Complex computation with meshgrid
  (print "Wave Interference Computation:")
  (let* ((x-coords (array-linspace 0.0 10.0 100))
         (y-coords (array-linspace 0.0 10.0 100)))
    
    (receive (x-grid y-grid) (array-meshgrid x-coords y-coords)
      ;; Two wave sources creating interference pattern
      ;; z = sin(r1) + sin(r2) where r1, r2 are distances from sources
      (let* ((source1-x 3.0) (source1-y 3.0)
             (source2-x 7.0) (source2-y 7.0)
             
             ;; Distances from sources
             (dx1 (array-sub-scalar x-grid source1-x))
             (dy1 (array-sub-scalar y-grid source1-y))
             (r1 (array-sqrt (array+ (array* dx1 dx1) (array* dy1 dy1))))
             
             (dx2 (array-sub-scalar x-grid source2-x))
             (dy2 (array-sub-scalar y-grid source2-y))
             (r2 (array-sqrt (array+ (array* dx2 dx2) (array* dy2 dy2))))
             
             ;; Interference pattern
             (pattern (array+ (array-sin r1) (array-sin r2))))
        
        (print "  Grid: 100 * 100 = " (array-size pattern) " points")
        (print "  Computing interference of two wave sources")
        (print "  Expression complexity: " (array-complexity pattern))
        
        (let ((start (current-process-milliseconds)))
          (let ((result (compute! pattern)))
            (let ((end (current-process-milliseconds)))
              (print "  Computation time: " (- end start) " ms")
              (print "  Pattern range: [" (array-min-val result) 
                     ", " (array-max-val result) "]")
              )))
        (print))))
  
  )

;; Run examples
(example-fusion-system)
(example-advanced-usage)

(example-logical-operations)
(example-meshgrid)