Fix for close line test. (#1565)

* Fix for close line test.  Appears that we don't actually want to look for close lines as much as check for numerical errors when computing the comving fredquency.  The test now occurs in calculate_distance_line.  The packet specs have been updated to remove the close line flags and the tests for the close lines have been removed

* Removed close_line parameter initialization for r_packet

* prevented test_montecarlo_main_loop from running comparison to C_tests when generating references

* Updated the value of the CLOSE_LINE_THRESHOLD to the double precision limit

* Removed statement claiming we are running C close lines (cause we are not anymore)

* Empty commit to rerun tests

* Updated montecarlo numba tests to do a full test of a tardis run and compare output estimators, energies, and frequencies.  This test should now compare a PR to the last stable tardis version

* useless commit to get the tests running again

* another empty commit for tests

* updated test so that it no longer import unecessary modules

* Cleaned up base test a tiny bit, really should be noted that this is a full test and isn't directly testing mainloop on its own.

tardis/montecarlo/montecarlo_numba/base.py

@@ -151,7 +151,6 @@ def montecarlo_main_loop(
     virt_packet_last_line_interaction_in_id = []
     virt_packet_last_line_interaction_out_id = []
 
-    print("Running post-merge numba montecarlo (with C close lines)!")
     for i in prange(len(output_nus)):
         if montecarlo_configuration.single_packet_seed != -1:
             seed = packet_seeds[montecarlo_configuration.single_packet_seed]
@@ -166,7 +165,6 @@ def montecarlo_main_loop(
             packet_collection.packets_input_energy[i],
             seed,
             i,
-            0,
         )
         vpacket_collection = VPacketCollection(
             r_packet.index,
@@ -202,7 +200,7 @@ def montecarlo_main_loop(
         ).astype(np.int64)
         # if we're only in a single-packet mode
         # if montecarlo_configuration.single_packet_seed == -1:
-        #     break
+        #    break
         for j, idx in enumerate(v_packets_idx):
             if (vpackets_nu[j] < spectrum_frequency[0]) or (
                 vpackets_nu[j] > spectrum_frequency[-1]

tardis/montecarlo/montecarlo_numba/calculate_distances.py

@@ -11,6 +11,7 @@
     MISS_DISTANCE,
     SIGMA_THOMSON,
     ENABLE_FULL_RELATIVITY,
+    CLOSE_LINE_THRESHOLD
 )
 
 from tardis.montecarlo.montecarlo_numba.utils import MonteCarloException
@@ -91,9 +92,8 @@ def calculate_distance_line(
     nu_diff = comov_nu - nu_line
 
     # for numerical reasons, if line is too close, we set the distance to 0.
-    if r_packet.is_close_line:
+    if abs(nu_diff/nu) < CLOSE_LINE_THRESHOLD:
         nu_diff = 0.0
-        r_packet.is_close_line = False
 
     if nu_diff >= 0:
         distance = (nu_diff / nu) * C_SPEED_OF_LIGHT * time_explosion

tardis/montecarlo/montecarlo_numba/formal_integral.py

@@ -109,7 +109,7 @@ def numba_formal_integral(model, plasma, iT, inu, inu_size, att_S_ul, Jred_lu, J
                     nu_ends, 
                     side='right'
             )
-            # loop over all interactions
+            # loop over all interactions 
             for i in range(size_z - 1):
                 escat_op = electron_density[int(shell_id[i])] * SIGMA_THOMSON
                 nu_end = nu_ends[i]
@@ -236,7 +236,7 @@ def check(self, raises=True):
 
         The function returns False if the configuration conflicts with the
         required settings. If raises evaluates to True, then a
-        IntegrationError is raised instead
+        IntegrationError is raised instead 
         """
 
         def raise_or_return(message):

tardis/montecarlo/montecarlo_numba/interaction.py

@@ -13,7 +13,6 @@
     angle_aberration_CMF_to_LF,
 )
 from tardis.montecarlo.montecarlo_numba.r_packet import (
-    test_for_close_line,
     InteractionType,
 )
 from tardis.montecarlo.montecarlo_numba.utils import get_random_mu
@@ -112,10 +111,6 @@ def line_emission(r_packet, emission_line_id, time_explosion, numba_plasma):
     r_packet.next_line_id = emission_line_id + 1
     nu_line = numba_plasma.line_list_nu[emission_line_id]
 
-    if emission_line_id != (len(numba_plasma.line_list_nu) - 1):
-        test_for_close_line(
-            r_packet, emission_line_id + 1, nu_line, numba_plasma
-        )
 
     if montecarlo_configuration.full_relativity:
         r_packet.mu = angle_aberration_CMF_to_LF(

tardis/montecarlo/montecarlo_numba/numba_config.py

@@ -1,7 +1,7 @@
 from tardis import constants as const
 
 SIGMA_THOMSON = const.sigma_T.to("cm^2").value
-CLOSE_LINE_THRESHOLD = 1e-7
+CLOSE_LINE_THRESHOLD = 1e-14
 C_SPEED_OF_LIGHT = const.c.to("cm/s").value
 MISS_DISTANCE = 1e99
 

tardis/montecarlo/montecarlo_numba/r_packet.py

@@ -53,7 +53,6 @@ class PacketStatus(IntEnum):
     ("status", int64),
     ("seed", int64),
     ("index", int64),
-    ("is_close_line", boolean),
     ("last_interaction_type", int64),
     ("last_interaction_in_nu", float64),
     ("last_line_interaction_in_id", int64),
@@ -63,7 +62,7 @@ class PacketStatus(IntEnum):
 
 @jitclass(rpacket_spec)
 class RPacket(object):
-    def __init__(self, r, mu, nu, energy, seed, index=0, is_close_line=False):
+    def __init__(self, r, mu, nu, energy, seed, index=0):
         self.r = r
         self.mu = mu
         self.nu = nu
@@ -72,7 +71,6 @@ def __init__(self, r, mu, nu, energy, seed, index=0, is_close_line=False):
         self.status = PacketStatus.IN_PROCESS
         self.seed = seed
         self.index = index
-        self.is_close_line = is_close_line
         self.last_interaction_type = -1
         self.last_interaction_in_nu = 0.0
         self.last_line_interaction_in_id = -1
@@ -141,7 +139,6 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators):
     cur_line_id = start_line_id  # initializing varibale for Numba
     # - do not remove
     last_line_id = len(numba_plasma.line_list_nu) - 1
-
     for cur_line_id in range(start_line_id, len(numba_plasma.line_list_nu)):
 
         # Going through the lines
@@ -218,10 +215,6 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators):
             distance = distance_trace
             break
 
-        if not is_last_line:
-            test_for_close_line(
-                r_packet, cur_line_id + 1, nu_line, numba_plasma
-            )
 
         # Recalculating distance_electron using tau_event -
         # tau_trace_line_combined
@@ -322,9 +315,3 @@ def move_packet_across_shell_boundary(packet, delta_shell, no_of_shells):
     else:
         packet.current_shell_id = next_shell_id
 
-
-@njit(**njit_dict_no_parallel)
-def test_for_close_line(r_packet, line_id, nu_line, numba_plasma):
-    r_packet.is_close_line = abs(
-        numba_plasma.line_list_nu[line_id] - nu_line
-    ) < (nu_line * CLOSE_LINE_THRESHOLD)

tardis/montecarlo/montecarlo_numba/tests/conftest.py

@@ -118,7 +118,6 @@ def packet(verysimple_packet_collection):
         energy=verysimple_packet_collection.packets_input_energy[0],
         seed=1963,
         index=0,
-        is_close_line=0,
     )
 
 
@@ -131,7 +130,6 @@ def static_packet():
         energy=0.9,
         seed=1963,
         index=0,
-        is_close_line=0,
     )
 
 

tardis/montecarlo/montecarlo_numba/tests/test_base.py

@@ -2,171 +2,70 @@
 import pandas as pd
 import os
 import numpy.testing as npt
-import numpy as np
 from copy import deepcopy
-from astropy import units as u
 
 from tardis.montecarlo import (
     montecarlo_configuration as montecarlo_configuration,
 )
-import tardis.montecarlo.montecarlo_numba.base as base
 from tardis.simulation import Simulation
-import tardis.montecarlo.montecarlo_numba.r_packet as r_packet
-import tardis.montecarlo.montecarlo_numba.single_packet_loop as spl
-from tardis.montecarlo.montecarlo_numba.numba_interface import (
-    PacketCollection,
-    VPacketCollection,
-    NumbaModel,
-    numba_plasma_initialize,
-    Estimators,
-    configuration_initialize,
-)
 
 
 @pytest.mark.xfail(reason="To be implemented")
 def test_montecarlo_radial1d():
     assert False
 
-
 def test_montecarlo_main_loop(
     config_verysimple,
     atomic_dataset,
     tardis_ref_path,
     tmpdir,
     set_seed_fixture,
     random_call_fixture,
+    request
 ):
 
     montecarlo_configuration.LEGACY_MODE_ENABLED = True
-
-    # Load C data from refdata
-    C_fname = os.path.join(tardis_ref_path, "montecarlo_1e5_compare_data.h5")
-    expected_nu = pd.read_hdf(
-        C_fname, key="/simulation/runner/output_nu"
-    ).values
-    expected_energy = pd.read_hdf(
-        C_fname, key="/simulation/runner/output_energy"
-    ).values
-    expected_nu_bar_estimator = pd.read_hdf(
-        C_fname, key="/simulation/runner/nu_bar_estimator"
-    ).values
-    expected_j_estimator = pd.read_hdf(
-        C_fname, key="/simulation/runner/j_estimator"
-    ).values
-
     # Setup model config from verysimple
     atomic_data = deepcopy(atomic_dataset)
     config_verysimple.montecarlo.last_no_of_packets = 1e5
     config_verysimple.montecarlo.no_of_virtual_packets = 0
     config_verysimple.montecarlo.iterations = 1
     config_verysimple.montecarlo.single_packet_seed = 0
+    config_verysimple.plasma.line_interaction_type = 'macroatom'
     del config_verysimple["config_dirname"]
 
     sim = Simulation.from_config(config_verysimple, atom_data=atomic_data)
+    sim.run()
 
-    # Init model
-    numba_plasma = numba_plasma_initialize(
-        sim.plasma, line_interaction_type="macroatom"
-    )
-
-    runner = sim.runner
-    model = sim.model
-
-    runner._initialize_geometry_arrays(model)
-    runner._initialize_estimator_arrays(numba_plasma.tau_sobolev.shape)
-    runner._initialize_packets(
-        model.t_inner.value, 100000, 0, model.r_inner[0].value
-    )
-
-    # Init parameters
-    montecarlo_configuration.v_packet_spawn_start_frequency = (
-        runner.virtual_spectrum_spawn_range.end.to(
-            u.Hz, equivalencies=u.spectral()
-        ).value
-    )
-    montecarlo_configuration.v_packet_spawn_end_frequency = (
-        runner.virtual_spectrum_spawn_range.start.to(
-            u.Hz, equivalencies=u.spectral()
-        ).value
-    )
-    montecarlo_configuration.temporary_v_packet_bins = 20000
-    montecarlo_configuration.full_relativity = runner.enable_full_relativity
-    montecarlo_configuration.single_packet_seed = 0
-
-    # Init packet collection from runner
-    packet_collection = PacketCollection(
-        runner.input_r,
-        runner.input_nu,
-        runner.input_mu,
-        runner.input_energy,
-        runner._output_nu,
-        runner._output_energy,
-    )
-
-    # Init model from runner
-    numba_model = NumbaModel(
-        runner.r_inner_cgs,
-        runner.r_outer_cgs,
-        model.time_explosion.to("s").value,
-    )
-
-    # Init estimators from runner
-    estimators = Estimators(
-        runner.j_estimator,
-        runner.nu_bar_estimator,
-        runner.j_blue_estimator,
-        runner.Edotlu_estimator,
-    )
-
-    # Empty vpacket collection
-    vpacket_collection = VPacketCollection(
-        0, np.array([0, 0], dtype=np.float64), 0, np.inf, 0, 0
-    )
-
-    # output arrays
-    output_nus = np.empty_like(packet_collection.packets_output_nu)
-    output_energies = np.empty_like(packet_collection.packets_output_nu)
-
-    # IMPORTANT: seeds RNG state within JIT
-    seed = 23111963
-    set_seed_fixture(seed)
-    for i in range(len(packet_collection.packets_input_nu)):
-        # Generate packet
-        packet = r_packet.RPacket(
-            packet_collection.packets_input_radius[i],
-            packet_collection.packets_input_mu[i],
-            packet_collection.packets_input_nu[i],
-            packet_collection.packets_input_energy[i],
-            seed,
-            i,
-            0,
-        )
 
-        # Loop packet
-        spl.single_packet_loop(
-            packet, numba_model, numba_plasma, estimators, vpacket_collection
-        )
-        output_nus[i] = packet.nu
-        if packet.status == r_packet.PacketStatus.REABSORBED:
-            output_energies[i] = -packet.energy
-        elif packet.status == r_packet.PacketStatus.EMITTED:
-            output_energies[i] = packet.energy
-
-        # RNG to match C
-        random_call_fixture()
+    compare_fname = os.path.join(tardis_ref_path, "montecarlo_1e5_compare_data.h5")
+    if request.config.getoption("--generate-reference"):
+        sim.to_hdf(compare_fname, overwrite=True)
+
+    # Load compare data from refdata
+    expected_nu = pd.read_hdf(
+        compare_fname, key="/simulation/runner/output_nu"
+    ).values
+    expected_energy = pd.read_hdf(
+        compare_fname, key="/simulation/runner/output_energy"
+    ).values
+    expected_nu_bar_estimator = pd.read_hdf(
+        compare_fname, key="/simulation/runner/nu_bar_estimator"
+    ).values
+    expected_j_estimator = pd.read_hdf(
+        compare_fname, key="/simulation/runner/j_estimator"
+    ).values
 
-    packet_collection.packets_output_energy[:] = output_energies[:]
-    packet_collection.packets_output_nu[:] = output_nus[:]
 
-    actual_energy = packet_collection.packets_output_energy
-    actual_nu = packet_collection.packets_output_nu
-    actual_nu_bar_estimator = estimators.nu_bar_estimator
-    actual_j_estimator = estimators.j_estimator
+    actual_energy = sim.runner.output_energy
+    actual_nu = sim.runner.output_nu
+    actual_nu_bar_estimator = sim.runner.nu_bar_estimator
+    actual_j_estimator = sim.runner.j_estimator
 
     # Compare
     npt.assert_allclose(
         actual_nu_bar_estimator, expected_nu_bar_estimator, rtol=1e-13
     )
     npt.assert_allclose(actual_j_estimator, expected_j_estimator, rtol=1e-13)
-    npt.assert_allclose(actual_energy, expected_energy, rtol=1e-13)
-    npt.assert_allclose(actual_nu, expected_nu, rtol=1e-13)
+    npt.assert_allclose(actual_energy.value, expected_energy, rtol=1e-13)
+    npt.assert_allclose(actual_nu.value, expected_nu, rtol=1e-13)

tardis/montecarlo/montecarlo_numba/tests/test_packet.py

@@ -402,16 +402,3 @@ def test_move_packet_across_shell_boundary_increment(
     assert packet.current_shell_id == current_shell_id + delta_shell
 
 
-@pytest.mark.parametrize(
-    ["line_id", "nu_line", "expected"],
-    [(5495, 1629252823683562.5, True), (3000, 0, False)],
-)
-def test_test_for_close_line(
-    packet, line_id, nu_line, verysimple_numba_plasma, expected
-):
-
-    r_packet.test_for_close_line(
-        packet, line_id, nu_line, verysimple_numba_plasma
-    )
-
-    assert packet.is_close_line == expected

tardis/montecarlo/montecarlo_numba/tests/test_vpacket.py

@@ -29,7 +29,6 @@ def v_packet():
         current_shell_id=0,
         next_line_id=0,
         index=0,
-        is_close_line=0,
     )
 
 
@@ -87,7 +86,6 @@ def test_trace_vpacket(
     npt.assert_almost_equal(v_packet.energy, 0.0)
     npt.assert_almost_equal(v_packet.mu, 0.8309726858508629)
     assert v_packet.next_line_id == 2773
-    assert v_packet.is_close_line == False
     assert v_packet.current_shell_id == 1
 
 
@@ -121,7 +119,6 @@ def broken_packet():
         mu=0.4916053094346575,
         energy=2.474533071386993e-07,
         index=3,
-        is_close_line=True,
         current_shell_id=0,
         next_line_id=5495,
     )