Undesired and inconsistent returning of NaN when instantiating a custom structure several times

General Usage
1

Hi,
This is my first post to the boards about a behavior I haven’t seen yet in a forum (appreciate any direction to a post if it exists!). I was attempting to build a b-spline structure atop a polynomial structure, and to do this I have to define a family of polynomials recursively. What I’m finding is that when I call the code from REPL, eg regbn(7), I get seemingly random occurences of NaN in the polynomial coefficients. Here is my attempt at a much simplified MWE demonstrating the issue:

""" Auxilliary structure routines """
struct Poly1d{T}
	deg::Int64
	c::Vector{T}

	function Poly1d(c::Vector)
		deg = convert(Int64,length(c))-1;
		S = typeof(c[1]);

		return new{S}(deg,c)
	end
end
function Base.:+(p₁::Poly1d,p₂::Poly1d)
	deg = p₁.deg<=p₂.deg ? p₂.deg : p₁.deg;
	T = typeof(p₁.c[1]+p₂.c[1]);
	c = Vector{T}(undef,deg+1);

	@inbounds for i=1:(deg+1)
		c[i] *= 0;
		c[i] += (p₁.deg+1)<=i ? p₁.c[i] : 0;
		c[i] += (p₂.deg+1)<=i ? p₂.c[i] : 0;
	end

	return Poly1d(c)
end

""" Recursion relation """
function regbⁿ(P::Vector{Poly1d{Float64}})
	n = length(P);
	Q = fill(Poly1d([0.0]),n+1);

	@inbounds for k=1:n+1
		if k!=n+1
			Q[k] += Poly1d([0.0,1.0])+P[k];
		end
		if k!=1
			Q[k] += Poly1d([-1.0,1.0])+P[k-1];
		end
	end
	return Q
end
function regbⁿ(deg::Integer)
	P = Vector{Vector{Poly1d{Float64}}}(undef,deg+1);
	P[1] = [Poly1d([1.0])];
	if deg==0
		return P[1]
	elseif deg>0
		@inbounds for k=1:deg
			flagnan = true
			while flagnan
				P[k+1] = regbⁿ(P[k])
				# Catch the erroneous NaN's
				ct = 0
				for i=1:length(P[k+1])
					flag = isnan.(P[k+1][i].c).==true
					if sum(flag)!=0
						ct+=1;
					end
				end
				if ct==0
					flagnan=false
				else
					@warn "At P[k+1]=P[$k+1] there were $ct many polys with nan's before continued" 
				end
			end
		end
		return P[end]
	end	
end

You’ll notice I’ve wrapped the part producing random NaN’s in a while loop that only terminates when no NaN’s are produced and warns the user how many times NaN’s occurred. I’ve also attached a screenshot of my REPL showing how the output changes despite being given the exact same input. I believe I have avoided any accidental argument mutation. The coefficients that do appear seem to be correct, and sometimes the code does run without returning any NaN’s. On the other hand, at other times calls will produce many more NaN’s than what the screenshot shows.

mew_issue_v1.7.3
mew_issue_v1.7.3886×1023 54.1 KB

In the screenshot you’ll see in the first and third calls, the NaN’s appeared in different polynomials (k=4 vs k=7) and in the second call no NaN was produced.

The while loop-fix may be a workable solution, but I wish I understood what was happening better to guard against in the future.

Appreciate any insights!

2

The vector c is initialised as undefined so it may indeed sometimes be NaN. Should you be setting it to zero instead? Either each entry in the loop, or using zeros when the vector is initialised.

1 Like
3

It’d be better to write this like

function Poly1d(c::Vector{T}) where T
    new{T}(length(c)-1, c) # `length` already returns an `Int`
end

This gives you a Poly1d{T} with the same type parameter as the element type of the vector. The same goes for your Base.:+ - you can use the type parameters here again:

function Base.:+(p1::Poly1d{T1}, p2::Poly1d{T2}) where {T1, T2}
    deg = max(p1.deg, p2.deg)
    rettype = promote_type(T1, T2)
    c = zeros(rettype, deg+1)

    for idx in eachindex(c)
        c[idx] += (p1.deg + 1) <= idx ? p1.c[idx] : zero(rettype)
        c[idx] += (p2.deg + 1) <= idx ? p2.c[idx] : zero(rettype)
    end

    return Poly1d(c)
end

Notice that I replaced the Vector{T}(undef, deg+1) with zeros(rettype, deg+1) - this is to prevent getting NaN from uninitialized memory and instead get a zero initialized array. You could also replace the c[i] *= 0 with just c[i] = zero(rettype), since NaN*0 == NaN, following IEEE floating point NaN semantics:

julia> NaN * 0
NaN
2 Likes
4

I think this still reads out of bounds. (p₁.deg+1)<=i ? p₁.c[i] : 0 tests that i is the last element of p₁.c or later, and then chooses to read from it. Deleting @inbounds (or using @Sukera’s version):

julia> Poly1d([1,2,3]) + Poly1d([4,5,6,7,8])
ERROR: BoundsError: attempt to access 3-element Vector{Int64} at index [4]

That may give you NaN (or segfaults) in addition to isnan(NaN * 0) as mentioned.

Edit: If the intention is to pad with zeros, then one tidy way is:

Base.:+(p1::Poly1d, p2::Poly1d) = Poly1d(padadd(p1.c, p2.c))

function padadd(x::AbstractVector, y::AbstractVector)
    length(x) >= length(y) || return padadd(y, x)
    [x[i] + get(y, i, zero(eltype(y))) for i in eachindex(x)]
end
2 Likes
5

I did have a bad feeling about it, thanks! This should be better:

for idx in eachindex(c)
    c[idx] += checkbounds(Bool, p1.c, idx) ? p1.c[idx] : zero(rettype)
    c[idx] += checkbounds(Bool, p2.c, idx) ? p2.c[idx] : zero(rettype)
end
6

Thanks everyone! Much clearer now. I hadn’t realized undef initialization can also generate NaN’s. And also @Sukera thanks for the style tips. I’m still building intuition for where-keyword uses as well as (custom) type promotions. I hadn’t directly used zeros (when I probably should have) because I didn’t know the cleanest way of initializing a zero element of a custom type. The comments help with that too.

Cheers

1 Like