NB. Standard Deviation Benchmark: Native J vs Futhark
NB. Compares performance on a compute-intensive operation

load '/path/to/futhark.ijs'

NB. ===============================================================
NB. Standard Deviation Implementations
NB. ===============================================================

NB. Native J standard deviation (population)
NB. stddev = sqrt(mean((x - mean(x))^2))
j_stddev =: 3 : 0
  mean =. (+/ y) % # y
  sq_diffs =. *: y - mean
  %: (+/ sq_diffs) % # y
)

NB. Futhark standard deviation
NB. More compute-intensive: mean calculation, element-wise ops, reduction, sqrt
fut_stddev =: 'entry f (xs: []f64) : f64 = let n = f64.i64 (length xs) in let mean = reduce (+) 0 xs / n in let sq_diffs = map (\x -> (x - mean) * (x - mean)) xs in f64.sqrt (reduce (+) 0 sq_diffs / n)' fut

NB. ===============================================================
NB. Run Benchmarks
NB. ===============================================================

echo ''
echo '=================================================================='
echo '     STANDARD DEVIATION BENCHMARK: Native J vs Futhark'
echo '=================================================================='
echo ''
echo '  stddev(x) = sqrt( mean( (x - mean(x))^2 ) )'
echo ''
echo '  Operations per element:'
echo '    - Compute mean (reduction over all elements)'
echo '    - Subtract mean from each element'
echo '    - Square each difference'
echo '    - Sum all squared differences (reduction)'
echo '    - Divide by n and take sqrt'
echo ''
echo '  This is MORE COMPUTE-INTENSIVE than dot product because:'
echo '    - Multiple passes over data'
echo '    - More arithmetic operations per element'
echo '    - Better opportunity for parallel speedup'
echo ''
echo '=================================================================='
echo ''
echo 'Backend: ', FUTHARK_BACKEND, '  Threads: ', ": FUTHARK_THREADS
echo ''

NB. Test sizes - from small to very large
sizes =: 1000 10000 100000 500000 1000000 2000000 5000000 10000000 20000000 50000000 100000000 

echo 'Size          J (ms)      Futhark (ms)  Speedup    Match'
echo '----------    --------    ------------  ---------  -----'

run_bench =: 3 : 0
  n =. y

  NB. Generate random data with known statistical properties
  data =. n ?@$ 0

  NB. Number of runs for averaging (fewer for large arrays)
  runs =. 3 >. <. 10000000 % n

  NB. J timing
  j_total =. 0
  for_r. i. runs do.
    t0 =. 6!:1 ''
    j_result =. j_stddev data
    t1 =. 6!:1 ''
    j_total =. j_total + 1000 * t1 - t0
  end.
  j_ms =. j_total % runs

  NB. Futhark timing
  fut_total =. 0
  for_r. i. runs do.
    t0 =. 6!:1 ''
    fut_result =. fut_stddev data
    t1 =. 6!:1 ''
    fut_total =. fut_total + 1000 * t1 - t0
  end.
  fut_ms =. fut_total % runs

  NB. Check results match (within floating point tolerance)
  rel_err =. | 1 - fut_result % j_result
  match =. rel_err < 1e_9

  NB. Calculate speedup (>1 means Futhark is faster)
  speedup =. j_ms % fut_ms >. 0.001

  NB. Format output
  size_s =. 10 {. ": n
  j_s =. 10 {. 8j2 ": j_ms
  fut_s =. 12 {. 8j2 ": fut_ms
  spd_s =. 9 {. (5j2 ": speedup), 'x'
  match_s =. > match { 'FAIL';'OK  '

  echo size_s, '    ', j_s, '  ', fut_s, '  ', spd_s, ' ', match_s

  n, j_ms, fut_ms, speedup, match
)

results =: run_bench"0 sizes

echo ''
echo '=================================================================='
echo 'RESULTS ANALYSIS'
echo '=================================================================='
echo ''

NB. Find crossover point
speedups =. 3 {"1 results
crossover =. 1 i.~ speedups >: 1

analyze =: 3 : 0
  'sizes speedups crossover' =. y
  if. crossover < # sizes do.
    echo 'Crossover point: ~', (": crossover { sizes), ' elements'
    echo '  Futhark becomes faster than J at this size'
  else.
    echo 'Crossover point: Not reached'
    echo '  J was faster for all tested sizes'
  end.
  echo ''

  NB. Show max speedup
  max_speedup =. >./ speedups
  max_idx =. speedups i. max_speedup
  echo 'Maximum speedup: ', (5j2 ": max_speedup), 'x at ', (": max_idx { sizes), ' elements'
  echo ''

  NB. Summary statistics
  echo 'Average speedup: ', (5j2 ": (+/ speedups) % # speedups), 'x'
  echo ''
)

analyze sizes;speedups;crossover

echo 'Key observations:'
echo '  - Speedup > 1.0x means Futhark is faster'
echo '  - Standard deviation is compute-bound (multiple ops per element)'
echo '  - Futhark can parallelize both reductions and element-wise ops'
echo '  - J''s implementation is sequential but highly optimized'
echo ''
echo '=================================================================='

exit 0