Clojure の mean を高速化できなかった

# まとめ
Clojureの速度が気になったのでJuliaを使ってみたが、Clojureが忘れられず、`partition`や`reductions`をJuliaで書いていた。そこに手間をかけるのであれば、Clojureの高速化にひと手間かけても良いのではないかと考えた。

しかし、シーケンスの平均値を取る操作を長いベクタに対して行った際、その実行時間を満足いくまで縮めることができなかった。更に縮めるためには恐らくデータ構造を変える必要があり、それはClojureの良さを殺してしまうように思われるので、やはりJuliaを使おう思う。

5000万個のランダム値に対する平均を取る操作の比較が下記。Clojureが4.7秒に対し、Juliaは0.37秒。素のPythonが26秒で、Numpyは0.39秒。Common Lispは大体1.75秒。

```clojure
;; Clojure
user=> (time (mean (for [_ (range 50000000)] (rand))))
"Elapsed time: 4762.215565 msecs"
```

```julia
# Julia
julia> @time mean([rand() for _ in 1:50000000])
  0.377605 seconds (52.42 k allocations: 383.969 MiB, 22.62% gc time, 9.38% compilation time)
```

```python
# Python
from time import time
from random import random
from statistics import mean
>>> t = time(); mean([random() for _ in range(0, 50000000)]); time() - t
25.926154136657715

# Numpy
import numpy as np
>>> t = time(); np.random.rand(50000000).mean(); time() - t
0.3986964225769043
```

```common-lisp
; Common Lisp
; 5000万個だと落ちるので、500万個
* (time (alexandria:mean     
    (let ((l nil))
      (dotimes (i 5000000 l)
        (push (random 1.0) l)))))
Evaluation took:
  0.175 seconds of real time
  0.172530 seconds of total run time (0.153032 user, 0.019498 system)
  [ Run times consist of 0.037 seconds GC time, and 0.136 seconds non-GC time. ]
  98.86% CPU
  469,180,265 processor cycles
  79,995,616 bytes consed
```

# Clojure高速化メモ

## REPL
大きな出力をさせようとすると、REPLがネックになることが多々ある。心当たりがある場合は、`take`等で出力範囲を区切る。

## プロファイル
まずは、コードのどこが遅いかを見付ける必要がある。選択肢は下記。

### `time`
時間を図りたい部分を`time`で括るだけ。結果の表示がうるさい/遅い場合は、`do`/`nil`も併用する。遅延評価部分は実行されないので、必要に応じて`doall`も併用する。
```clojure
(time (for [_ (range 10000)] (range 10000)))
; (out) "Elapsed time: 0.185629 msecs"
; REPLに結果を表示させようとするため、実際はしばらく待たないと結果が表示されず、REPLも使えない。

(time (do (for [_ (range 10000)] (range 10000)) nil))
; (out) "Elapsed time: 0.180818 msecs"
; (do ... nil) で括ると出力が無くなるため、実行後すぐにREPLを使える。

(time (do (doall (for [_ (range 10000)] (range 10000))) nil))
; (out) "Elapsed time: 2.578813 msecs"
; doallを使うと、遅延評価部分が現実化されるため、実行時間が長くなる。
; REPLに表示させる場合は、常にこの現実化+REPL用の処理時間がかかる。
```
下記のような関数にまとめてしまっても良いかもしれない。
```clojure
(defn dotime [& something]
  (time (apply doall something))
  nil)

(dotime (for [_ (range 10000)] (range 10000)))
; (out) "Elapsed time: 3.307578 msecs"
```

### tufte
[参考記事](https://qiita.com/totakke/items/e0464da6a269cfc5de6b)では[criterium](https://github.com/hugoduncan/criterium) (☆1.1k) が利用されているが、2021年11月から1年半ほど動きがないので、[tufte](https://github.com/ptaoussanis/tufte)を試す。こちらは2022年12月更新の、☆486。READMEを見る限り、こちらの方がcriteriumよりも細かく情報を取れそう。
ターミナルからREPLで`clj`している場合は良いが、[Conjure](https://github.com/Olical/conjure)を利用している場合、Conjure出力画面ではなく裏で動いているターミナルに結果が表示されるので注意。
```clojure
;; deps.edn
  :deps {
    com.taoensso/tufte {:mvn/version "2.4.5"}
  }
```
計測した部分を`(p ...)`で括って、最初の引数にキーワードを渡してあげるだけ。最後に全体を`(profile {} ...)`で括る。
```clojure
(ns user
  (:require [taoensso.tufte :as tufte :refer [defnp p profile]]))

(tufte/add-basic-println-handler!
  {:format-pstats-opts {:columns [:n-calls :min :max :mean :mad :clock :total]}})

(profile
  {}
  (let [a (p :a (range 10000))]
    (p :body (range 10000))))
;; (省略)
;; pId           nCalls        Min        Max       Mean   MAD      Clock  Total
;; 
;; :a                 1     1.60μs     1.60μs     1.60μs   ±0%     1.60μs     1%
;; :body              1   387.00ns   387.00ns   387.00ns   ±0%   387.00ns     0%
;; 
;; Accounted                                                       1.99μs     1%
;; Clock                                                         180.80μs   100%
```
なお、`profile`の後ろの`{}`には、プロファイルのオプションを入れることができる。そしてtufteについても、関数を作っておく。
```clojure
(defn doprof [& something]
  (profile {} (apply doall something) nil))
```

## 高速化

### `mean`
Clojureは標準で`mean`が実装されていないので、自分で簡単なものを書いていたが、これがボトルネックだった。なので、逐次的に値を求めるような書き方に変更をした。
下記の通りいくつかのパターンを試した (実際はそれぞれ何回かサンプルを取った) ところ、ライブラリのものが少し他より速かった。
```clojure
;; これが遅い
(defn mean [coll]
  (/ (apply + coll) (count coll)))
; (time (mean (for [_ (range 1000000)] (rand))))
; (out) "Elapsed time: 106.225616 msecs"

;; reduceで一回だけcollを眺めれば良いように変更
(defn mean [coll]
  (let [r (reduce (fn [mem val] {:S (+ (mem :S) val) :n (inc (mem :n))})
                  {:S 0 :n 0}
                  coll)]
    (/ (r :S) (r :n))))
; (out) "Elapsed time: 70.677517 msecs"

;; recurバージョン
(defn mean 
  ([coll] (mean coll 0 0))
  ([coll s n]
    (if (empty? coll)
      (/ s n)
      (recur (rest coll) (+ s (first coll)) (inc n)))))
; (out) "Elapsed time: 79.161107 msecs"

;; recurベクタバージョン
(defn mean 
  ([coll] (mean (vec coll) 0 0))
  ([coll s n]
    (if (empty? coll)
      (/ s n)
      (recur (pop coll) (+ s (peek coll)) (inc n)))))
; (out) "Elapsed time: 207.029417 msecs"

;; ライブラリ
;; net.mikera/core.matrix {:mvn/version "0.63.0"}
(ns user
  (:require [clojure.core.matrix.stats :refer [mean]]))
; (out) "Elapsed time: 70.736877 msecs"
```
わざわざ書くほどでもないが、メモ。$M$=平均、$S$=合計。
$$
\begin{align*}
M_n &= \frac{a_1 + a_2 + \cdots + a_n}{n} = \frac{S_n}{n} \\\\
S_n &= S_{n-1}+a_n \\\\
M_n &= \frac{S_{n-1}+a_n}{n}
\end{align*}
$$

上記の速度はJuliaに比べると、十分に遅いように思われる。
```julia
# Julia
using Statistics
@time mean([rand() for _ in 1:1000000])
# 0.033709 seconds (52.42 k allocations: 10.131 MiB, 83.37% compilation time)
```