调试函数时，发现NaN判断不了

综合讨论区 Julia入门
1

我发现下面的函数测试可以，换了个场景运行不了。于是进行测试：

function getWakeVoltage(zvec::Vector{Float64}, ρvec::Vector{Float64}, funz2wakez::Union{Interpolations.Extrapolation,Function}, Δz, Nb, elementcharge=1.6021766208e-19)
	volvec = zeros(size(zvec, 1))
	println("1.0=", volvec[1:10])
	for i in 1:length(zvec)
		s = 0.0
		@fastmath for j in 1:length(zvec)
			s += funz2wakez((i-j)*Δz) * ρvec[j] * Δz
		end
		volvec[i] = s

		if isnan(s)
			println(i)
			break
		else
			println("s=", s, ", isnan?: ", isnan(s), ", isequal?: ", isequal(s, NaN))
		end
	end
	println("1.1=", volvec[1:10])
	-Nb*elementcharge*volvec
end

这次报错，发现volvec最后都变成NaN。于是尝试检查s的值。结果：

s=NaN, isnan?: false, isequal?: false

也就是s是NaN，但是没有判断出来。不知道为啥。
如果有大佬能够一眼看出发现我鲁棒性差的原因，就更好了。

2

直接测试 s = NaN 无法复现。

versioninfo() 以及系统环境说一下。
另外最好给一个能运行、可复现 bug 的示例。
尽量多化简一下，不然不太好调试。

3

我的版本信息

julia> versioninfo()
Julia Version 1.6.7
Commit 3b76b25b64 (2022-07-19 15:11 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: AMD Ryzen 5 2600 Six-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, znver1)
Environment:
  JULIA_NUM_THREADS = 8
  JULIA_PKG_SERVER = https://mirrors.bfsu.edu.cn/julia/static

然后是我这里在REPL中赋值给s，是能够用isnan判断的

julia> s=NaN
NaN

julia> isnan(s)
true

我提供下我测试的代码：

computername = ENV["COMPUTERNAME"]
if computername == "家里电脑"
    push!(LOAD_PATH, "E:/Documents/JuliaProgram/MyModules/")
elseif computername == "宿舍电脑"
    push!(LOAD_PATH, "F:/Documents/JuliaProgramme/MyModules")
elseif computername=="DESKTOP-FDIN0I7"
    push!(LOAD_PATH, "D:/Documents/Julia/MyModules")
elseif computername=="TINA"
    push!(LOAD_PATH, "D:/JuliaProgramme/MyModules")
end

using DoubleRFs, CSV, DataFrames, Interpolations, BenchmarkTools, Plots

function getWakeFunction(wakepath::String)
    df=CSV.read(wakepath, DataFrame)
	LinearInterpolation(df[!,"z"],1e12*df[!,"wakez"])
end

function main()
    funwakez = getWakeFunction("./wakezlin.csv")
    σz = 0.00252
    Nb = 5e-9/1.6e-19
    v1 = 116.66328e6
    v2 = 0.0
    r = v2/v1
    ϕs = 2.8155729
    ϕ2s = 0.0
    h1 = 216820
    h = 3
    ηp = 1.429249e-5
    σδ = 3.836791e-4
    Δz = 0.01σz
    centerenergy = 45.5e9
    circum=9.9635e4
    Σ = 0.5
    χ = 0.01
    denerr = 1e-12
    
    gr()

    # 全程不变的量
    zvec = collect(-Σ:Δz:Σ)
    potenVecCavity = getPotenCavity(zvec, ϕs, ϕ2s, v1, r, h1, h, circum, centerenergy)
    zR, criPotential = DoubleRFs.getRBound(zvec, potenVecCavity, ϕs, ϕ2s, v1, r, h1, h, circum, centerenergy, σz) # 头部粒子存在与否的边界，及该点临界potential
    
    zL = DoubleRFs.getLBound(zvec, potenVecCavity, criPotential, σz)
    #    由于右侧(头部)不会受到尾场影响，因此只需要求一次右侧z和临界的Hamiltonian即可。
    #    左侧（尾部）则为了保险起见每次更新。
    ρvec=DoubleRFs.getDistFromPotential(zvec, potenVecCavity, zL, zR, ηp, σδ) # 计算初始密度分布

    # gr()
    # plt1 = plot(zvec, potenVecCavity, label="Cavity Potential")
    # plot!(plt1, zvec, getPotenWake(zvec, ρvec, funwakez, Nb, Δz, centerenergy, circum), label="")
    # xlims!(plt1, -0.04, 0.04)
    # ylims!(plt1, -1e-11, 1e-10)
    # display(plt1)

    plt=plot(zvec, ρvec, label="")
    title!(plt, "χ="*string(χ))
    xlims!(plt, -0.015, 0.015)

    # 第一次迭代
    ρvecpre = ρvec
    println("第1次调用")
    ρvecnow = DoubleRFs.getρVecNow(ρvec, zvec, potenVecCavity, funwakez, zR, criPotential, Nb, centerenergy, circum, ηp, σδ, σz, Δz)
    
    ρvec = (1-χ).*ρvecpre .+ χ.*ρvecnow   # 这是利用默认权重χ， 直接迭代的分布， 下一行则是相应的初始判据
    cri = DoubleRFs.convergecriterion(ρvec, ρvecpre, zvec, χ, zL, zR)

    # 开始循环迭代
    while χ >= 0.001
        ρvecpre = ρvec
        println("第2次调用")
        ρvecnow = DoubleRFs.getρVecNow(ρvec, zvec, potenVecCavity, funwakez, zR, criPotential, Nb, centerenergy, circum, ηp, σδ, σz, Δz)

        # 找到合适的χ，使分布收敛
        suc = false
        while suc == false
            ρvec = (1-χ).*ρvecpre .+ χ.*ρvecnow
            tmpcri = DoubleRFs.convergecriterion(ρvec, ρvecpre,  zvec, χ, zL, zR)
            if tmpcri < cri     # 找到的χ能否让判据收敛
                suc = true
                cri = tmpcri
                println("χ=", χ, ", criterion=", cri)
            elseif χ <= 0.001   # 如果χ过小，则不用再寻找合适的χ了
                break
            else
                χ = χ*0.99
            end
        end
    end
    # zvec, ρvec
    plot!(plt, zvec, ρvec, label="")
    display(plt)
    nothing
end

main()

而这里的DoubleRFs库，被我上传到了Github，大佬可以下载下来像我一样push!调用。
这次遇到的问题函数是./src/Voltage.jl文件中的getWakeVoltage的第一个方法。
我测试通过的代码和记录也被我放在./test/中的TestVoltage.jl和TestVoltage.md中。

4

呃呃，大佬，我找到了我鲁棒性差的问题，是因为我的ρvec变[NaN NaN ... NaN]了。

但是isnan()无法判断的问题就交给你了。

5

你似乎没有提供这个文件。

6

抱歉。

链接: 百度网盘 请输入提取码 提取码: 1234 复制这段内容后打开百度网盘手机App，操作更方便哦

其实这个的问题是我读取这个wakelin.csv文件时，忘了和之前有所不同。所以保留了1e12这个系数，最终导致ρvec都变成NaN了。

那么问题可以简化一下：为什么

function getWakeVoltage(zvec::Vector{Float64}, ρvec::Vector{Float64}, funz2wakez::Union{Interpolations.Extrapolation,Function}, Δz, Nb, elementcharge=1.6021766208e-19)
	volvec = zeros(size(zvec, 1))
	println("1.0=", volvec[1:10])
	for i in 1:length(zvec)
		s = 0.0
		@fastmath for j in 1:length(zvec)
			s += funz2wakez((i-j)*Δz) * ρvec[j] * Δz
		end
		volvec[i] = s

		if isnan(s)
			println(i)
			break
		else
			println("s=", s, ", isnan?: ", isnan(s), ", isequal?: ", isequal(s, NaN))
		end
	end
	println("1.1=", volvec[1:10])
	-Nb*elementcharge*volvec
end

中ρvec为[NaN NaN ... NaN]时，我isnan(s)没判断出来？

7

盲猜一下，可能和这个宏有关。

julia> @fastmath isnan(NaN) |> println
false

julia> isnan(NaN) |> println
true

QQ图片20220911085356
QQ图片20220911085356597×548 31.6 KB

都给你优化了，你就说快不快吧 :rofl:

@fastmath 是以放弃精度为代价，去追求极致的速度。
如果你是做科学计算的，建议在任何语言中都不要用这个选项。

I mean, the whole point of fast-math is trading off speed with correctness. If fast-math was to give always the correct results, it wouldn’t be fast-math, it would be the standard way of doing math.

其中有一个

Allow optimizations for floating-point arithmetic that assume that arguments and results are not NaNs or ±Infs.

看起来很像我们这里的情况。

https://simonbyrne.github.io/notes/fastmath/

目前还剩一个问题，是这个优化是怎么传播到下方的 block 中的。
isnan 并不在优化的范围内。

2 个赞
8

:joy:确实快

好这就取消@fastmath。

找Julia bug的问题，就交给你们专业人士了。