Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove inaccurate validation tests. #1428

Merged
merged 7 commits into from
Nov 18, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,8 @@ types of particle simulations.
Additional information and publications.
- [HOOMD-blue benchmark scripts](https://github.com/glotzerlab/hoomd-benchmarks):
Scripts to evaluate the performance of HOOMD-blue simulations.
- [HOOMD-blue validation tests](https://github.com/glotzerlab/hoomd-validation):
Scripts to validate that HOOMD-blue performs accurate simulations.

## Related tools

Expand Down
104 changes: 0 additions & 104 deletions hoomd/conftest.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,6 @@
import atexit
import os
import numpy
import math
import warnings
from hoomd.logging import LoggerCategories
from hoomd.snapshot import Snapshot
from hoomd import Simulation
Expand Down Expand Up @@ -548,108 +546,6 @@ def autotuned_kernel_parameter_check(instance, activate, all_optional=False):
assert instance.kernel_parameters == initial_kernel_parameters


class BlockAverage:
"""Block average method for estimating standard deviation of the mean.

Args:
data: List of values
"""

def __init__(self, data):
# round down to the nearest power of 2
N = 2**int(math.log(len(data)) / math.log(2))
if N != len(data):
warnings.warn(
"Ignoring some data. Data array should be a power of 2.")

block_sizes = []
block_mean = []
block_variance = []

# take means of blocks and the mean/variance of all blocks, growing
# blocks by factors of 2
block_size = 1
while block_size <= N // 8:
num_blocks = N // block_size
block_data = numpy.zeros(num_blocks)

for i in range(0, num_blocks):
start = i * block_size
end = start + block_size
block_data[i] = numpy.mean(data[start:end])

block_mean.append(numpy.mean(block_data))
block_variance.append(numpy.var(block_data) / (num_blocks - 1))

block_sizes.append(block_size)
block_size *= 2

self._block_mean = numpy.array(block_mean)
self._block_variance = numpy.array(block_variance)
self._block_sizes = numpy.array(block_sizes)
self.data = numpy.array(data)

# check for a plateau in the relative error before the last data point
block_relative_error = numpy.sqrt(self._block_variance) / numpy.fabs(
self._block_mean)
relative_error_derivative = (numpy.diff(block_relative_error)
/ numpy.diff(self._block_sizes))
if numpy.all(relative_error_derivative > 0):
warnings.warn("Block averaging failed to plateau, run longer")

def get_hierarchical_errors(self):
"""Get details on the hierarchical errors."""
return (self._block_sizes, self._block_mean, self._block_variance)

@property
def standard_deviation(self):
"""float: The error estimate on the mean."""
if numpy.all(self.data == self.data[0]):
return 0

return numpy.sqrt(numpy.max(self._block_variance))

@property
def mean(self):
"""float: The mean."""
return self._block_mean[-1]

@property
def relative_error(self):
"""float: The relative error."""
return self.standard_deviation / numpy.fabs(self.mean)

def assert_close(self,
reference_mean,
reference_deviation,
z=6,
max_relative_error=0.02):
"""Assert that the distribution is constent with a given reference.

Also assert that the relative error of the distribution is small.
Otherwise, test runs with massive fluctuations would likely lead to
passing tests.

Args:
reference_mean: Known good mean value
reference_deviation: Standard deviation of the known good value
z: Number of standard deviations
max_relative_error: Maximum relative error to allow
"""
sample_mean = self.mean
sample_deviation = self.standard_deviation

assert sample_deviation / sample_mean <= max_relative_error

# compare if 0 is within the confidence interval around the difference
# of the means
deviation_diff = ((sample_deviation**2
+ reference_deviation**2)**(1 / 2.))
mean_diff = math.fabs(sample_mean - reference_mean)
deviation_allowed = z * deviation_diff
assert mean_diff <= deviation_allowed


class ListWriter(hoomd.custom.Action):
"""Log a single quantity to a list.

Expand Down
1 change: 0 additions & 1 deletion hoomd/hpmc/pytest/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@ set(files __init__.py
test_small_box_2d.py
test_small_box_3d.py
conftest.py
test_nec.py
depletants_sphere.py
)

Expand Down
112 changes: 0 additions & 112 deletions hoomd/hpmc/pytest/test_nec.py

This file was deleted.

1 change: 0 additions & 1 deletion hoomd/md/pytest/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@ set(files __init__.py
test_dihedral.py
test_flags.py
test_kernel_parameters.py
test_lj_equation_of_state.py
test_potential.py
test_pppm_coulomb.py
test_manifolds.py
Expand Down
120 changes: 0 additions & 120 deletions hoomd/md/pytest/test_alj.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,133 +2,13 @@
# Part of HOOMD-blue, released under the BSD 3-Clause License.

import numpy as np
import numpy.testing as npt
import pytest

import hoomd
import hoomd.conftest
from hoomd import md


@pytest.mark.validate
def test_conservation(simulation_factory, lattice_snapshot_factory):
# For test, use a unit area hexagon.
particle_vertices = np.array([[6.20403239e-01, 0.00000000e+00, 0],
[3.10201620e-01, 5.37284966e-01, 0],
[-3.10201620e-01, 5.37284966e-01, 0],
[-6.20403239e-01, 7.59774841e-17, 0],
[-3.10201620e-01, -5.37284966e-01, 0],
[3.10201620e-01, -5.37284966e-01, 0]])
area = 1.0
circumcircle_radius = 0.6204032392788702
incircle_radius = 0.5372849659264116
num_vertices = len(particle_vertices)

circumcircle_diameter = 2 * circumcircle_radius

# Just initialize in a simple cubic lattice.
sim = simulation_factory(
lattice_snapshot_factory(a=4 * circumcircle_diameter,
n=10,
dimensions=2))
sim.seed = 175

# Initialize moments of inertia since original simulation was HPMC.
mass = area
# https://math.stackexchange.com/questions/2004798/moment-of-inertia-for-a-n-sided-regular-polygon # noqa
moment_inertia = (mass * circumcircle_diameter**2 / 6)
moment_inertia *= (1 + 2 * np.cos(np.pi / num_vertices)**2)

with sim.state.cpu_local_snapshot as snapshot:
snapshot.particles.mass[:] = mass
snapshot.particles.moment_inertia[:] = np.array([0, 0, moment_inertia])
# Not sure if this should be incircle or circumcircle;
# probably doesn't matter based on current usage, but may
# matter in the future for the potential if it's modified
# to actually use diameter.
snapshot.particles.diameter[:] = circumcircle_diameter
kT = 0.3
sim.state.thermalize_particle_momenta(hoomd.filter.All(), kT)

# Create box resize updater
packing_fraction = 0.4
final_area = area * sim.state.N_particles / packing_fraction
L_final = np.sqrt(final_area)
final_box = hoomd.Box.square(L_final)

n_compression_start = int(1e4)
n_compression_end = int(1e5)
n_compression_total = n_compression_end - n_compression_start

box_resize = hoomd.update.BoxResize(
box1=sim.state.box,
box2=final_box,
trigger=int(n_compression_total / 10000),
variant=hoomd.variant.Ramp(0, 1, n_compression_start,
n_compression_total),
filter=hoomd.filter.All())
sim.operations += box_resize

# Define forces and methods
r_cut_scale = 1.3
kernel_scale = (1 / np.cos(np.pi / num_vertices))
incircle_diameter = 2 * incircle_radius
r_cut_set = incircle_diameter * kernel_scale * r_cut_scale

alj = md.pair.aniso.ALJ(default_r_cut=r_cut_set, nlist=md.nlist.Cell(0.4))

alj.shape["A"] = {
"vertices": particle_vertices,
"faces": [],
"rounding_radii": 0
}

alpha = 0 # Make it WCA-only (no attraction)
eps_att = 1.0
alj.params[("A", "A")] = {
"epsilon": eps_att,
"sigma_i": incircle_diameter,
"sigma_j": incircle_diameter,
"alpha": alpha
}

nve = md.methods.NVE(filter=hoomd.filter.All())
integrator = md.Integrator(dt=1e-4,
forces=[alj],
methods=[nve],
integrate_rotational_dof=True)
sim.operations.integrator = integrator

# Compress box
sim.run(n_compression_end)

thermo = md.compute.ThermodynamicQuantities(hoomd.filter.All())
sim.operations += thermo

# Reset velocities after the compression, and equilibriate
sim.state.thermalize_particle_momenta(hoomd.filter.All(), kT)
sim.run(1000)

# run sim and get values back
w = hoomd.conftest.ManyListWriter([(thermo, 'potential_energy'),
(thermo, 'kinetic_energy'),
(integrator, 'linear_momentum')])
writer = hoomd.write.CustomWriter(action=w, trigger=1)
sim.operations.writers.append(writer)
sim.run(1000)
pe, ke, momentum = w.data
total_energies = np.array(pe) + np.array(ke)

# Ensure energy conservation up to the 3 digit per-particle.
npt.assert_allclose(total_energies,
total_energies[0],
atol=0.03 * sim.state.N_particles)

# Test momentum conservation.
p_magnitude = np.linalg.norm(momentum, axis=-1)
npt.assert_allclose(p_magnitude, p_magnitude[0], atol=1e-13)


def test_type_shapes(simulation_factory, two_particle_snapshot_factory):
alj = md.pair.aniso.ALJ(md.nlist.Cell(buffer=0.1))
sim = simulation_factory(two_particle_snapshot_factory(d=2.0))
Expand Down
Loading