What is causing type instability in this struct?

General Usage
1

I have the following struct and I’m a little confused why using @code_warntype on it’s construct says there’s unstable types. The output from that command is a little verbose but I’ll copy one of the red parts in case that is useful. Any help would be appreciated, thanks!

mutable struct System{T, M, MS, N, L, K}
    coords::Matrix{T}
    coords_initial::Matrix{T}
    coords_unwrapped::Matrix{T}
    masses::SVector{N,M}
    masses_sqrt::SVector{N,MS}
    box_size::L
    n_particles::Integer
    kB::K
end

function System(r0, masses, box_size, N_atoms, kB)
    masses_sqrt = sqrt.(masses)
    
    return System{eltype(r0), eltype(masses), eltype(masses_sqrt), N_atoms, typeof(box_size), typeof(kB)}(
        deepcopy(r0), deepcopy(r0), deepcopy(r0), masses, masses_sqrt, box_size, N_atoms, kB)
end

│ %22 = Core.apply_type(Main.System, %17, %18, %19, N_atoms, %20, %21)::Type{System{Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothi
ng}}, Quantity{Float64, 𝐌 𝐍^-1, Unitful.FreeU
nits{(g, mol^-1), 𝐌 𝐍^-1, nothing}}, Quantity
{Float64, 𝐌^1/2 𝐍^-1/2, Unitful.FreeUnits{(g^
1/2, mol^-1/2), 𝐌^1/2 𝐍^-1/2, nothing}}, _A,
Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), �
, nothing}}, Quantity{Float64, 𝐋^2 𝐌 𝐍^-1 �
^^-1 𝐓^-2, Unitful.FreeUnits{(kcal, K^-1, mol^
-1), 𝐋^2 𝐌 𝐍^-1 𝚯^-1 𝐓^-2, nothing}}}} whe
re _A

2

Integer is an abstract type. It comprises all possible integers, such as Int32 or Int64. If you want a concrete type, use one of the types just mentioned or use Int which defaults to your system integer width, usually Int64 on most computers.

3

That didn’t seem to fix it.

In the code_warntype output it goes through all of the parameters and n_particles was not flagged. My feeling is that the issue is with the SVectors (maybe mass_sqrt) but I can’t really interpret the output too well.

4

That’s hard to tell without seeing code. Can you provide a minimum example? (Have to leave the pc now so won’t be able to help further.)

5

No rush,

this is similar to how I initialize and it also gives type instability in the print out

L =  20.9872u"Å"
N_atoms = 256
kB =  1.9872e-3u"kcal * mol^-1 * K^-1"
masses = @SVector ones(N_atoms)
r0 = zeros(N_atoms,3) * unit(L)
sys = @code_warntype System(r0, masses, L, N_atoms, kB)
1 Like
6

Not sure it helps but deepcopy is usually type-unstable

1 Like
7

It will be type unstable as you are creating a struct using an integer in the type, but this information is only available at runtime, causing the function to be unstable. To be type stable, the compiler needs all the information necessary to deduce the concrete types of all objects at compile time, which only has access to the types of the variables, not the contents.

You could fix this by using Val types for the integers. E.g.

function System(r0, masses, box_size, N_atoms::Val{N}, kB) where {N}
    masses_sqrt = sqrt.(masses)
    
    return System{eltype(r0), eltype(masses), eltype(masses_sqrt), N, typeof(box_size), typeof(kB)}(
        deepcopy(r0), deepcopy(r0), deepcopy(r0), masses, masses_sqrt, box_size, N, kB)
end

And modify the call like this:

@code_warntype System(r0, masses, L, Val(N_atoms), kB)
8

@code_warntype System(r0, r0, r0, masses, deepcopy(masses), L, N_atoms, kB)
returns perfectly blue output, seems like everything can be inferred?

1 Like
9

Yeah so passing 3 deep copies to the constructor and changing Integer to Int64 works, but why does calling deepcopy in the constructor differ from doing it outside the constructor.

Even though that works, I do think its kind of ugly and bad style to pass 3 versions of the same array into the constructor and from a user POV they should not need to know the struct needs 3 copies of the same array. That should be handled in the constructor.

1 Like
10

Yeah, so the problem is that you are using a value, N_atoms, as a type variable, even though it is not in the type domain.

But I also don’t get why you need N_atoms as a type parameter. Isn’t that already encoded in the type parameter N, which is part of SVector{N,M}? It seems redundant and possibly error prone to duplicate this. Furthermore, N_atoms is also encoded in the field n_particle, which means it’s encoded three different places in your struct, and since it is mutable, it is possible for the field value to come out of sync with the type parameter(s). (Edit: Actually, the number is encoded in six different places, since it’s also in the size of the fields coords, coords_initial and coords_unwrapped.)

But there are other things.

  • StaticArrays are intended for small arrays, and 256 is not small, you might as well just use regular arrays.
  • You don’t need deepcopy, just regular copy
  • Encoding both masses and masses_sqrt is redundant, and prone to becoming inconsistent, if one array is modified and not the other.
  • Coordinates can be nicely encoded by a Vector{SVector{3, Float64}}.
  • I wonder what all the different matrices are for anyway.

But, that’s exactly what you did, you handled it in the constructor. What do you mean here?

11

Maybe you could make a Particle type, and then defined System somewhat like this:

struct System{P<:Particle, L, K}
    particles::Vector{P}
    box_size::L
    kB::K
end

And then n_particles(sys::System) = length(sys.particles), etc.

12

Yeah the struct isn’t the prettiest, its mostly just data storage. The coordinates are without explaining necessary to track separately (and copy should work), but initialized from the same array. The mass and mass_sqrt are both stored to avoid re-calculating. I end up using both arrays a LOT and its cheaper to just store the sqrt version of the masses here.

I’ll remove the SVectors like you suggested. They are probably unnecessary. Guess I overcomplicated things lol.

Thanks all.

13

What about making a Particle type that contains positions and mass? That seems tidier to me.

For example, if you want to add or remove a particle from the system, you could do

add!(sys::System, p::Particle) = push!(sys.particles, p)

etc.

Adding a particle to your current System is quite complex.

All sorts of operations are likely to become easier to handle.

2 Likes
14

Yeah that’s probably something I should do, I was focused on getting the code working first I guess.

I could implement something like to get the masses…
masses(sys::System) = mass.(sys.particles)