Single precision vs double precision in create_state_from_snapshot()

[trvsst]

Hi hoomd community,

I noticed that when using create_state_from_snapshot() the resulting state is saved in single precision.

Below is a snippet, where I get the charge from a snapshot I create and np.sum(ch1) gives 1e-15. The same charge after create-from_snapshot (np.sum(ch2) gives 2e-5. Just confirming that this is the expected behaviour. I checked and the same loss of precision happens with the mass. The only issue is that when the charges add up to 2e-5 hoomd issues a warning in that the charge is not zero. Any comments are welcome.

.. code-block:: python
ch1 = snap.particles.charge
ms1 = snap.particles.mass
print('charges', np.sum(ch1))
print('masses', np.sum(ms1))
sim = hoomd.Simulation(device=dv, seed=sd)
sim.create_state_from_snapshot(snap)
ch2 = sim.state.get_snapshot().particles.charge
ms2 = sim.state.get_snapshot().particles.mass
print('charges 2', np.sum(ch2))
print('masses 2', np.sum(ms2))

[joaander]

I am not able to reproduce the behavior you describe when I run your script, I get

NameError: name 'snap' is not defined

Provided that you compiled HOOMD-blue with HOOMD_LONGREAL_SIZE=64 (the default) and snap is created via hoomd.Snapshot or state.get_snapshot then I do not expect the behavior you describe. Snapshot creates a double precision snapshot:

hoomd-blue/hoomd/snapshot.py

Line 95 in ec0b84c

self._cpp_obj = _hoomd.SnapshotSystemData_double()

hoomd-blue/hoomd/SnapshotSystemData.cc

Lines 169 to 171 in ec0b84c

	pybind11::class_<SnapshotSystemData<double>, std::shared_ptr<SnapshotSystemData<double>>>(
	m,
	"SnapshotSystemData_double")

and get_snapshot returns a double precision snapshot:

hoomd-blue/hoomd/SystemDefinition.cc

Line 299 in ec0b84c

.def("takeSnapshot_double", &SystemDefinition::takeSnapshot<double>)

hoomd-blue/hoomd/state.py

Line 308 in ec0b84c

cpp_snapshot = self._cpp_sys_def.takeSnapshot_double()

You can confirm this by checking the numpy array type of the snapshot:

>>> import hoomd
>>> s = hoomd.Snapshot()
>>> s.particles.N = 1
>>> s.particles.charge.dtype
dtype('float64')

You can also confirm that double precision values survive a round trip through create_state_from_snapshot/set_snapshot/get_snapshot.

>>> import hoomd
>>> import numpy
>>> 
>>> s = hoomd.Snapshot()
>>> s.particles.N = 1
>>> s.particles.types = ['A']
>>> s.particles.charge[0] = 1 + numpy.finfo(float).eps
>>> print(s.particles.charge[0])
1.0000000000000002
>>> 
>>> simulation = hoomd.Simulation(device=hoomd.device.CPU())
>>> simulation.create_state_from_snapshot(s)
>>> print(simulation.state.get_snapshot().particles.charge[0])
1.0000000000000002
>>> 
>>> simulation.state.set_snapshot(s)
>>> print(simulation.state.get_snapshot().particles.charge[0])
1.0000000000000002
>>> 

Note that GSD frames store values in single precision: https://gsd.readthedocs.io/en/v3.4.2/schema-hoomd.html

[trvsst]

Hi Josh,

Thanks for the detailed answer. I was expecting this would be a much trivial issue, so I thought with a snippet it would be enough.

Here is a more detailed answer. I have generated a snapshot with another program (hoodlt), so what I have done is to pickle the object,
which is the file

snap.p

The script

test_sp.py

does the following:
unpickle snap.p, which gives the original Frame object in double precision.

use the snapshot in create_from_snapshot and you can see that the result is in single precision
I also save it as gsd and the result is also in single precision, but not the same as create_from_snapshot

The script is quite verbose so as to make it self-explanatory. It is puzzling to me how it happens that the original snapshot ends up in single precision. Is it possible that the original Frame object (as read from the pickle file) has some flag that triggers hoomd or gsd to save it in single precision? I really have no clue what is it that can be happening.

files.zip

[joaander]

A float64 can represent any float32 value exactly. Therefore, your individual float32 inputs from snap_p are preserved all the way through the chain:

$ python -i test_sp.py
...
>>> np.max(np.abs(np.array(snap_p.particles.charge, dtype=np.float64) - snap.particles.charge))
np.float64(0.0)

Your sums over many quantities do not show that the precision is not maintained. What is happening is that your float32 sum experiences round-off error. The float64 sum experiences less round-off error and is therefore closer to to the correct sum:

$ python -i test_sp.py
...
>>> np.sum(np.array(snap_p.particles.charge, dtype=np.float64)) - np.sum(snap_p.particles.charge)
np.float64(3.1948089599609375e-05)

As an exercise, try permuting the order of the float32 array and performing the sum in float32 - you should find that the sum depends on the order of the elements.

You should prepare your input masses and charges with double precision computations if you would like them to have a specific sum.

[trvsst]

Thanks a lot Josh,

I am sure you are right, but I am a little confused by your answer. When I print the output in my machine, it shows that the pickled snap

snap_p.particles.charge

is of type float64

What do you mean by "Therefore, your individual float32 inputs from snap_p", as far as I see, there are no float32 inputs.

It is snap that somehow becomes float32.

I added a couple of lines to the script explicitly showing this:

snap_v2.zip

[joaander]

I am sure you are right, but I am a little confused by your answer. When I print the output in my machine, it shows that the pickled snap

I am confused by your statement:

$ python -i test_sp.py
>>> snap_p
<gsd.hoomd.Frame object at 0x11d82c980>
>>> snap_p.particles.mass.dtype
dtype('float32')
>>> snap_p.particles.charge.dtype
dtype('float32')

snap_p is clearly a gsd.hoomd.Frame object which I already stated should be single precision. You also confirmed in words that you expect a Frame object in the pickle file. The code block above was executed with your original snap.zip.

In your v2 script, you print the numpy dtype in snap_p earlier in the script and it is float64. I did not expect this. GSD stores single precision values. You should never place values in a Frame that you want to stay in float64 as GSD converts them to the expected type:

https://github.com/glotzerlab/gsd/blob/9c6868705dd4e4f256b1fe84cf614f47049eb56a/gsd/hoomd.py#L223-L228

This code is in the validate method which is called on the Frame in create_state_from_snapshot.

In your v2 script, you write a gsd file and read it back in:

print('Save and then read a gsd file, result has the same precision as create from snapshot, but it is not the same')
f = gsd.hoomd.open(name='snap.gsd', mode='w')
f.append(snap_p)
f.close()
with gsd.hoomd.open(name='snap.gsd', mode='r') as fi:
    snap = fi[0]
print('snap numpy type is', snap.particles.charge.dtype)

Again, the float32 output is expected. GSD SAVES FILES IN SINGLE PRECISION. I gave this answer in the first reply. This is expected behavior. If you would like to store values in double, use a different file format.

[trvsst]

Thanks a lot Josh. My confusion was because snap_p.particles.charge came as float64 when it should be single precision. Now I understand.
Thanks a lot again and have a happy Thanksgiving.