
# ------------------------------------------------------------------------------------------
# # Julia 很快
#
# 基准数值测试程序（benchmarks）通常被用来比较编程语言的性能。
# 这些 benchmarks 可能引起很多的讨论，首先是作为基准测试的到底是什么，再有就是究竟是什么导致不同语言间的差异。
# 这些简单的问题有时比你乍想之下复杂得多。
#
# 这个 notebook 目的在于让你实现一个简单的基准测试。
# 通过这个 notebook 可以看到作者的配有四核 Macbook Pro 表现如何，或者你可以自己跑跑看。
#
# （材料来自于 MIT 的 Steven Johnson 的一场演讲：[18S096/Boxes-and-registers.ipynb at master ·
# mitmath/18S096](https://github.com/mitmath/18S096/blob/master/lectures/lecture1/Boxes-and-
# registers.ipynb)）
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 目录
#
# - [`sum` 函数的定义](#sum:-一个简单易懂的函数)
# - [不同语言中的各种基准测试](#不同语言中的各种基准测试)
#     1. [C (手写)](#1.-C-语言)
#     2. [C (手写 + `-ffast-math`)](#2.-使用--ffast-math-的-C-语言)
#     3. [python (内置)](#3.-Python-的内置-sum)
#     4. [python (numpy)](#4.-Python:-numpy)
#     5. [python (手写)](#5.-Python（手写）)
#     6. [Julia (内置)](#6.-Julia（内置）)
#     7. [Julia (手写)](#7.-Julia（手写）)
#     8. [Julia (手写 + SIMD)](#8.-Julia（手写-+--@simd）)
# - [性能测试小结](#性能测试小结)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## `sum`: 一个简单易懂的函数
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 考虑**求和**函数 `sum(a)`，它是用来计算
# $$
# \mathrm{sum}(a) = \sum_{i=1}^n a_i
# $$
# 其中 $n$ 是 `a` 的长度。
# ------------------------------------------------------------------------------------------

a = rand(10^7)  # 在[0,1)上均匀分布的随机数组成的1维向量

sum(a)

# ------------------------------------------------------------------------------------------
# 预期结果应该在 $0.5 \times 10^7$ 左右，因为所有元素的平均值应该是 0.5 左右
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 不同语言中的各种基准测试
# ------------------------------------------------------------------------------------------

@time sum(a)

@time sum(a)

@time sum(a)

# ------------------------------------------------------------------------------------------
# `@time` 宏会产生受干扰的结果，所以不是做基准测试的最佳选择！
#
# 好在 Julia 有个 `BenchmarkTools.jl` 包来简单并准确地进行基准测试：
# ------------------------------------------------------------------------------------------

using Pkg
Pkg.add("BenchmarkTools")

using BenchmarkTools

# ------------------------------------------------------------------------------------------
# ###  1. C 语言
#
# C 语言对人类不友好但是对机器友好，经常作为一个黄金标准。
# 在 C 语言的两倍以内通常来说就比较理想了。
# 然而，即使仅考虑 C 语言，不同的 C 语言程序员也可能采用或不采用不同的优化方式。
#
# 原作者不会 C 语言，所以他看不懂下面这个代码块中的代码，
# 但是他知道 Julia 中可以编译、运行 C 语言代码。注意 `"""` 所包围的是个多行字符串。
# ------------------------------------------------------------------------------------------

using Libdl
C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()  # 创建一个临时文件


# 通过将 C_code 导入 gcc 来编译成共享链接库
# (只有装了 gcc 才行):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end

# 定义一个调用这个 C 语言函数的 Julia 函数：
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_sum(a)

c_sum(a) ≈ sum(a)  # 输入\approx后按<TAB>键得到≈符号

c_sum(a) - sum(a)  

≈  # `isapprox` 函数的别名

?isapprox

# ------------------------------------------------------------------------------------------
# 现在我们可以在 Julia 中直接对 C 语言代码进行基准测试了：
# ------------------------------------------------------------------------------------------

c_bench = @benchmark c_sum($a)

println("C: 最快 $(minimum(c_bench.times) / 1e6) 毫秒")

d = Dict()  # 字典，也就是关联数组（键值对）
d["C"] = minimum(c_bench.times) / 1e6  # 毫秒
d

using Plots
gr()

using Statistics  # 导入统计包以计算标准差
t = c_bench.times / 1e6  # 毫秒
m, σ = minimum(t), std(t)

histogram(t, bins=500,
    xlim=(m - 0.01, m + σ),
    xlabel="milliseconds", ylabel="count", label="")

# ------------------------------------------------------------------------------------------
# ### 2. 使用 `-ffast-math` 的 C 语言
#
# 如果我们允许 C 重新安排浮点操作，它就会使用 SIMD 指令(single instruction, multiple data)进行矢量化。
# ------------------------------------------------------------------------------------------

const Clib_fastmath = tempname()   # 创建临时文件

# 同上但是增加了 -ffast-math 标志
open(`gcc -fPIC -O3 -msse3 -xc -shared -ffast-math -o $(Clib_fastmath * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end

# 定义一个调用这个 C 语言函数的 Julia 函数：
c_sum_fastmath(X::Array{Float64}) = ccall(("c_sum", Clib_fastmath), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_fastmath_bench = @benchmark $c_sum_fastmath($a)

d["C -ffast-math"] = minimum(c_fastmath_bench.times) / 1e6  # 毫秒

# ------------------------------------------------------------------------------------------
# ### 3. Python 的内置 `sum`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `PyCall` 包给 Julia 提供了 Python 接口：
# ------------------------------------------------------------------------------------------

using Pkg; Pkg.add("PyCall")
using PyCall

# 获得 Python 的内置 "sum" 函数：
pysum = pybuiltin("sum")

pysum(a)

pysum(a) ≈ sum(a)

py_list_bench = @benchmark $pysum($a)

d["Python built-in"] = minimum(py_list_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 4. Python: `numpy`
#
# **条件允许的情况下利用硬件 "SIMD"**
#
# `numpy` 是一个 Python 的优化过的 C 语言库。
# 在 Julia 版 Python 中的安装方法为：
# ------------------------------------------------------------------------------------------

using Pkg; Pkg.add("Conda")
using Conda

Conda.add("numpy")

numpy_sum = pyimport("numpy")["sum"]

py_numpy_bench = @benchmark $numpy_sum($a)

numpy_sum(a)

numpy_sum(a) ≈ sum(a)

d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 5. Python（手写）
# ------------------------------------------------------------------------------------------

py"""
def py_sum(A):
    s = 0.0
    for a in A:
        s += a
    return s
"""

sum_py = py"py_sum"

py_hand = @benchmark $sum_py($a)

sum_py(a)

sum_py(a) ≈ sum(a)

d["Python hand-written"] = minimum(py_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 6. Julia（内置）
#
# **是直接用 Julia 写的，不是 C！**
# ------------------------------------------------------------------------------------------

@which sum(a)

j_bench = @benchmark sum($a)

d["Julia built-in"] = minimum(j_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 7. Julia（手写）
# ------------------------------------------------------------------------------------------

function mysum(A)   
    s = 0.0  # s = zero(eltype(a))
    for a in A
        s += a
    end
    s
end

j_bench_hand = @benchmark mysum($a)

d["Julia hand-written"] = minimum(j_bench_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ### 8. Julia（手写 +  `@simd`）
# ------------------------------------------------------------------------------------------

function mysum_simd(A)   
    s = 0.0  # s = zero(eltype(A))
    @simd for a in A
        s += a
    end
    s
end

j_bench_hand_simd = @benchmark mysum_simd($a)

mysum_simd(a)

d["Julia hand-written simd"] = minimum(j_bench_hand_simd.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# ## 性能测试小结
# ------------------------------------------------------------------------------------------

for (key, value) in sort(collect(d), by=last)
    println(rpad(key, 25, "."), lpad(round(value; digits=1), 6, "."))
end
