Fix a number of bugs/corner cases

hyperion/conf/conf_files.py

@@ -512,7 +512,7 @@ def set_convergence(self, convergence, percentile=100., absolute=0., relative=0.
             The relative threshold below which the ratio in the percentile
             value has to be for convergence.
         '''
-        self.check_convergence = True
+        self.check_convergence = convergence
         self.convergence_percentile = percentile
         self.convergence_absolute = absolute
         self.convergence_relative = relative

hyperion/conf/tests/test_run_conf.py

@@ -22,3 +22,12 @@ def test_propagation_check_frequency_invalid2(value):
     with pytest.raises(ValueError) as exc:
         r.set_propagation_check_frequency(value)
     assert exc.value.args[0] == "frequency should be between 0 and 1"
+
+
+def test_set_convergence_flag():
+    r = RunConf()
+    r.set_convergence(True, percentile=99., absolute=2., relative=1.02)
+    assert r.check_convergence is True
+    r = RunConf()
+    r.set_convergence(False)
+    assert r.check_convergence is False

hyperion/densities/alpha_disk.py

@@ -448,6 +448,9 @@ def density(self, grid):
         # Find density factor
         rho *= self.rho_0
 
+        if np.sum(rho * grid.volumes) == 0. and self.mass > 0:
+            raise Exception("Discretized disk mass is zero, suggesting that the grid is too coarse")
+
         norm = self.mass / np.sum(rho * grid.volumes)
 
         logger.info("Disk density is being re-scaled by a factor of %.2f to give the correct mass." % norm)

hyperion/densities/ambient_medium.py

@@ -132,8 +132,8 @@ def dust(self, value):
 
     def _check_all_set(self):
 
-        if self.density is None:
-            raise Exception("density is not set")
+        if self.rho is None:
+            raise Exception("rho is not set")
         if self.rmin is None:
             raise Exception("rmin is not set")
         if self.rmax is None:

hyperion/densities/tests/test_densities.py

@@ -835,3 +835,56 @@ def test_ambient_densities_4():
 
     expected = np.repeat(2., 9)
     assert_array_almost_equal_nulp((p1.density(g) + c.density(g) + a.density(g))[0, 0, :], expected, 10)
+
+
+def test_ulrich_midplane_cumulative_density_inner():
+    # Regression test for a factor of 2 overestimate of the midplane
+    # cumulative density inside the centrifugal radius
+    e = UlrichEnvelope()
+    e.star = Star()
+    e.star.mass = 0.5
+    e.star.radius = 0.5
+    e.rho_0 = 1.
+    e.rc = 1.
+    e.rmin = 0.01
+    e.rmax = 10.
+    r = np.logspace(np.log10(e.rmin), np.log10(0.9 * e.rc), 100000)
+    rho_mid = e.rho_0 / np.sqrt(r / e.rc) / (1. - r / e.rc) / 2.
+    expected = np.sum(0.5 * (rho_mid[1:] + rho_mid[:-1]) * np.diff(r))
+    actual = e.midplane_cumulative_density(np.array([r[-1]]))[0]
+    assert abs(actual - expected) / expected < 1.e-3
+
+
+def test_ambient_medium_missing_rho():
+    a = AmbientMedium()
+    a.rmin = 1.
+    a.rmax = 10.
+    with pytest.raises(Exception) as exc:
+        a._check_all_set()
+    assert exc.value.args[0] == "rho is not set"
+
+
+def test_alpha_disk_zero_discretized_mass():
+    # A grid too coarse to resolve the disk gives a zero discretized mass,
+    # which used to silently produce NaN densities
+    e = AlphaDisk()
+    e.star = Star()
+    e.star.mass = 0.5
+    e.star.radius = 0.5
+    e.rmin = 1.
+    e.rmax = 2.
+    e.mass = 1.
+    e.r_0 = 5.
+    e.h_0 = 1.
+    e.p = -1.
+    e.beta = 1.25
+    e.lvisc = G * (1.5 - 2. / np.sqrt(2.) + 0.5)
+
+    r = np.array([0., 0.5, 100.])
+    t = np.linspace(0., np.pi, 3)
+    p = np.linspace(0., 2. * np.pi, 3)
+    g = SphericalPolarGrid(r, t, p)
+
+    with pytest.raises(Exception) as exc:
+        e.density(g)
+    assert "zero" in exc.value.args[0]

hyperion/densities/ulrich_envelope.py

@@ -442,7 +442,7 @@ def midplane_cumulative_density(self, r):
             # entries are overwritten just below, so the warnings can be
             # ignored.
             with np.errstate(invalid='ignore'):
-                rho[:] = (self.rho_0 * self.rc
+                rho[:] = (0.5 * self.rho_0 * self.rc
                           * (np.log((np.sqrt(gamma_1) + 1) / (1. - np.sqrt(gamma_1)))
                              - np.log((np.sqrt(gamma_0) + 1) / (1. - np.sqrt(gamma_0)))))
 

hyperion/dust/dust_type.py

@@ -747,7 +747,7 @@ def __init__(self, model):
             if np.any(np.isnan(values)):
                 logger.warning("NaN values found inside MieX %s file - interpolating" % quantity)
                 invalid = np.isnan(values)
-                values[invalid] = interp1d_fast_loglog(wav[~invalid], values[~invalid], wav[invalid])
+                values[invalid] = interp1d_fast_loglog(wav[~invalid][::-1], values[~invalid][::-1], wav[invalid])
                 if np.any(np.isnan(values)):
                     raise Exception("Did not manage to fix NaN values in MieX %s" % quantity)
 
@@ -808,7 +808,7 @@ def __init__(self, model):
                 if np.any(np.isnan(values[:, i])):
                     logger.warning("NaN values found inside MieX %s file - interpolating" % quantity)
                     invalid = np.isnan(values[:, i])
-                    values[:, i][invalid] = interp1d_fast_loglog(wav[~invalid], values[:, i][~invalid], wav[invalid])
+                    values[:, i][invalid] = interp1d_fast_loglog(wav[~invalid][::-1], values[:, i][~invalid][::-1], wav[invalid])
                     if np.any(np.isnan(values[:, i])):
                         raise Exception("Did not manage to fix NaN values in MieX %s" % quantity)
 

hyperion/filter/filter.py

@@ -136,7 +136,7 @@ def from_hdf5_group(cls, group, name):
 
         self.alpha = group[name].attrs['alpha']
         self._beta = group[name].attrs['beta']
-        self.central_spectral_coords = group[name].attrs['nu0'] * u.Hz
+        self.central_spectral_coord = group[name].attrs['nu0'] * u.Hz
 
         return self
 

hyperion/filter/tests/test_filter.py

@@ -104,3 +104,5 @@ def test_roundtrip():
                                f2.spectral_coord.to(u.Hz).value)
     np.testing.assert_allclose(f.transmission.to(u.percent).value,
                                f2.transmission.to(u.percent).value)
+    np.testing.assert_allclose(f.central_spectral_coord.to(u.Hz).value,
+                               f2.central_spectral_coord.to(u.Hz).value)

hyperion/grid/amr_grid.py

@@ -1,6 +1,7 @@
 from __future__ import print_function, division
 
 import os
+import posixpath
 import struct
 import hashlib
 from copy import deepcopy
@@ -509,13 +510,13 @@ def __setitem__(self, item, value):
                     grid_ref.quantities[item] = deepcopy(grid.quantities[value.viewed_quantity])
         elif isinstance(value, h5py.ExternalLink):
             filename = value.filename
-            base_path = os.path.dirname(value.path)
-            array_name = os.path.basename(value.path)
+            base_path = posixpath.dirname(value.path)
+            array_name = posixpath.basename(value.path)
             for ilevel, level_ref in enumerate(self.levels):
                 level_path = 'level_%05i' % (ilevel + 1)
                 for igrid, grid_ref in enumerate(level_ref.grids):
-                    grid_path = 'grid_%05i' % (ilevel + 1)
-                    grid_ref.quantities[item] = h5py.ExternalLink(filename, os.path.join(base_path, level_path, grid_path, array_name))
+                    grid_path = 'grid_%05i' % (igrid + 1)
+                    grid_ref.quantities[item] = h5py.ExternalLink(filename, posixpath.join(base_path, level_path, grid_path, array_name))
         elif value == []:
             for level in self.levels:
                 for grid in level.grids:

hyperion/grid/tests/test_views.py

@@ -253,3 +253,22 @@ def test_dustless_array(self, grid_type):
         g2['density1'] = g1['density1']
         g2['density2'] = g1['density1'][0]
         assert g2.n_dust == 1
+
+
+def test_amr_external_link_paths():
+    # Regression test for links that used the level index instead of the
+    # grid index in the per-grid path
+    import h5py
+    amr = AMRGrid()
+    level = amr.add_level()
+    for i in range(2):
+        grid = level.add_grid()
+        grid.xmin, grid.xmax = -1., 1.
+        grid.ymin, grid.ymax = -1., 1.
+        grid.zmin, grid.zmax = -1., 1.
+        grid.nx, grid.ny, grid.nz = 4, 4, 4
+    amr['density'] = h5py.ExternalLink('other_model.rtin', 'Grid/Quantities/density')
+    link_1 = amr.levels[0].grids[0].quantities['density']
+    link_2 = amr.levels[0].grids[1].quantities['density']
+    assert link_1.path.endswith('level_00001/grid_00001/density')
+    assert link_2.path.endswith('level_00001/grid_00002/density')

hyperion/grid/tests/test_voronoi.py

@@ -134,3 +134,16 @@ def test_init_sparse(tmpdir):
     g3 = VoronoiGrid(g2)
     g3.write(f3)
 
+
+
+def test_setitem_from_view():
+    # Regression test for VoronoiGrid.__setitem__, which used to check a
+    # nonexistent 'refined' attribute when assigning a view to an empty grid
+    g1 = VoronoiGrid([1., 2., 3.], [3., 4., 5.], [2., 3., 4.])
+    g1.quantities['density'] = np.array([1., 2., 3.])
+    g2 = VoronoiGrid()
+    g2['density'] = g1['density']
+    np.testing.assert_allclose(g2['density'].array, g1['density'].array)
+    assert g2.xmin == g1.xmin and g2.xmax == g1.xmax
+    assert g2.ymin == g1.ymin and g2.ymax == g1.ymax
+    assert g2.zmin == g1.zmin and g2.zmax == g1.zmax

hyperion/grid/voronoi_grid.py

@@ -545,9 +545,12 @@ def __getitem__(self, item):
 
     def __setitem__(self, item, value):
         if isinstance(value, VoronoiGridView):
-            if self.refined is None:
+            if self._x is None:
                 logger.warning("No geometry in target grid - copying from original grid")
-                self.set_points(value.x, value.y, value.z)
+                self.set_points(value.x, value.y, value.z,
+                                xmin=value.xmin, xmax=value.xmax,
+                                ymin=value.ymin, ymax=value.ymax,
+                                zmin=value.zmin, zmax=value.zmax)
             self.quantities[item] = deepcopy(value.quantities[value.viewed_quantity])
         elif isinstance(value, h5py.ExternalLink):
             self.quantities[item] = value

hyperion/grid/voronoi_helpers.py

@@ -197,12 +197,12 @@ def plot(self):
             # Drop the empty vertices coordinates, signalled by NaN.
             arr = v[~np.isnan(v)]
             assert(len(arr) % 3 == 0)
-            tmp = np.split(arr, len(arr) / 3)
+            tmp = np.split(arr, len(arr) // 3)
             # Append the vertices.
             points = points + tmp
             # Append the cell description.
             cells = cells + \
-                [len(tmp)] + range(cur_cell_idx, cur_cell_idx + len(tmp))
+                [len(tmp)] + list(range(cur_cell_idx, cur_cell_idx + len(tmp)))
             cur_cell_idx += len(tmp)
             # Append the offset info.
             offset.append(cur_offset)

hyperion/importers/orion.py

@@ -58,9 +58,9 @@ def read_data(self, filename, offset, quantity_indices, verbose=False):
         p8 = header.find('))', p7)
         bytes = [int(x) for x in header[p7:p8].split()]
 
-        if bytes == range(1, n_bytes + 1):
+        if bytes == list(range(1, n_bytes + 1)):
             endian = '>'
-        elif bytes == range(n_bytes, 0, -1):
+        elif bytes == list(range(n_bytes, 0, -1)):
             endian = '<'
         else:
             raise Exception("Unexpected bytes: %s" % str(bytes))
@@ -72,7 +72,7 @@ def read_data(self, filename, offset, quantity_indices, verbose=False):
         for quantity in quantity_indices:
 
             f.seek(pos + quantity_indices[quantity] * n_bytes * gridsize)
-            array = np.fromstring(f.read(n_bytes * gridsize),
+            array = np.frombuffer(f.read(n_bytes * gridsize),
                                   dtype='%sf%i' % (endian, n_bytes))
             self.quantities[quantity] = array.reshape(self.nz, self.ny, self.nx)
 

hyperion/model/model.py

@@ -125,8 +125,8 @@ def _read_monochromatic(self, group):
         self._monochromatic = str2bool(group.attrs['monochromatic'])
         if self._monochromatic:
             self._frequencies = np.array(group['frequencies']['nu'])
-            if 'energy_threshold' in group.attrs:
-                self._monochromatic_energy_threshold = group.attrs['energy_threshold']
+            if 'monochromatic_energy_threshold' in group.attrs:
+                self._monochromatic_energy_threshold = group.attrs['monochromatic_energy_threshold']
             else:
                 self._monochromatic_energy_threshold = 1.e-10
 
@@ -472,8 +472,8 @@ def use_image_config(self, filename):
             g_image = f['/Input/Output/']
 
         # Read in binned images
-        if 'n_theta' in g_image['Binned']:
-            self.binned_output = BinnedImageConf.read(g_image['Binned'])
+        if 'group_00001' in g_image['Binned']:
+            self.binned_output = BinnedImageConf.read(g_image['Binned']['group_00001'])
 
         # Read in peeled images
         for peeled in g_image['Peeled']:

hyperion/model/model_output.py

@@ -340,6 +340,8 @@ def get_sed(self, stokes='I', group=0, technique='peeled',
         if technique == 'peeled':
             n_groups = len(self.file['Peeled'])
             if (group < 0 and -group <= n_groups) or (group >= 0 and group < n_groups):
+                if group < 0:
+                    group = n_groups + group
                 g = self.file['Peeled/group_%05i' % (group + 1)]
             else:
                 raise ValueError('File only contains %i image/SED group(s)' % n_groups)
@@ -662,6 +664,8 @@ def get_image(self, stokes='I', group=0, technique='peeled',
         if technique == 'peeled':
             n_groups = len(self.file['Peeled'])
             if (group < 0 and -group <= n_groups) or (group >= 0 and group < n_groups):
+                if group < 0:
+                    group = n_groups + group
                 g = self.file['Peeled/group_%05i' % (group + 1)]
             else:
                 raise ValueError('File only contains %i image/SED group(s)' % n_groups)
@@ -958,6 +962,8 @@ def get_available_components(self, iteration=-1):
         # If iteration is last one, find iteration number
         if iteration == -1:
             iteration = find_last_iteration(self.file)
+        else:
+            iteration = iteration + 1  # Python value is zero based
 
         # Return components
         components = list(self.file['iteration_%05i' % iteration].keys())

hyperion/model/tests/test_get_quantities.py

@@ -18,3 +18,22 @@ def test_no_initial(tmpdir):
     mo = m.run(tmpdir.join(random_id()).strpath)
     g = mo.get_quantities()
     assert 'density' in g
+
+
+@pytest.mark.requires_hyperion_binaries
+def test_get_available_components_zero_based(tmpdir):
+    # Regression test for get_available_components using 1-based iteration
+    # numbering while get_quantities is zero-based
+    m = Model()
+    m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.])
+    m.add_density_grid(np.array([[[1.]]]), get_test_dust())
+    source = m.add_point_source()
+    source.luminosity = 1.
+    source.temperature = 6000.
+    m.set_n_initial_iterations(1)
+    m.set_n_photons(initial=1, imaging=0)
+    m.write(tmpdir.join(random_id()).strpath)
+    mo = m.run(tmpdir.join(random_id()).strpath)
+    components = mo.get_available_components(iteration=0)
+    assert 'specific_energy' in components
+    assert mo.get_available_components() == components

hyperion/model/tests/test_model.py

@@ -5,6 +5,7 @@
 import shutil
 from copy import deepcopy
 
+import h5py
 import numpy as np
 import pytest
 
@@ -396,3 +397,39 @@ def test_binned_forced_first_interaction(tmpdir):
     with pytest.raises(Exception) as exc:
         m.write(tmpdir.join(random_id()).strpath)
     assert exc.value.args[0] == "can't use binned images with forced first interaction - use set_forced_first_interaction(False) to disable"
+
+
+def test_monochromatic_energy_threshold_round_trip(tmpdir):
+    # Regression test for the threshold being read back with the wrong
+    # attribute name and silently reset to the default
+    m = Model()
+    m.set_monochromatic(True, wavelengths=[1., 2.], energy_threshold=1.e-5)
+    f = h5py.File(tmpdir.join(random_id()).strpath, 'w')
+    m._write_monochromatic(f)
+    m2 = Model()
+    m2._read_monochromatic(f)
+    f.close()
+    assert m2._monochromatic_energy_threshold == 1.e-5
+
+
+def test_use_image_config_binned(tmpdir):
+    # Regression test for binned image configuration never being restored
+    m = Model()
+    m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.])
+    m.add_density_grid(np.array([[[1.]]]), get_test_dust())
+    s = m.add_point_source()
+    s.luminosity = 1.
+    s.temperature = 6000.
+    i = m.add_binned_images(sed=True, image=False)
+    i.set_wavelength_range(5, 1, 10)
+    i.set_viewing_bins(3, 2)
+    m.set_forced_first_interaction(False)
+    m.set_n_photons(initial=1, imaging=1)
+    filename = tmpdir.join(random_id()).strpath
+    m.write(filename)
+
+    m2 = Model()
+    m2.use_image_config(filename)
+    assert m2.binned_output is not None
+    assert m2.binned_output.n_theta == 3
+    assert m2.binned_output.n_phi == 2