Ffi2

nrcatalogtools/__init__.py

@@ -2,6 +2,7 @@
 nr-catalog-tools is a toolkit for interfacing with gravitational-wave
 catalogs generated via Numerical Relativity simulations
 """
+
 from __future__ import absolute_import
 
 from . import lvc, maya, metadata, rit, sxs, utils, waveform

nrcatalogtools/catalog.py

@@ -106,12 +106,18 @@ def get(self, sim_name, quantity="waveform"):
                         " {}".format(self.psi4_url_from_simname(sim_name))
                     )
                 self.download_psi4_data(sim_name)
-            try:
-                return waveform.WaveformModes.load_from_h5(filepath, metadata=metadata)
-            except OSError:
+
+            from nrcatalogtools.rit import RITCatalog
+
+            if isinstance(self, RITCatalog):
                 return waveform.WaveformModes.load_from_targz(
                     filepath, metadata=metadata
                 )
+
+            else:
+                # MAYA catalog
+                return waveform.WaveformModes.load_from_h5(filepath, metadata=metadata)
+
         else:
             raise IOError(
                 f"Cannot provide quantity: {quantity}. Only supported options are [waveform, psi4]"

nrcatalogtools/maya.py

@@ -181,19 +181,19 @@ def _add_paths_to_metadata(self):
         if any([col not in existing_cols for col in new_cols]):
             for sim_name in metadata_dict:
                 if "metadata_location" not in existing_cols:
-                    metadata_dict[sim_name][
-                        "metadata_location"
-                    ] = self.metadata_filepath_from_simname(sim_name)
+                    metadata_dict[sim_name]["metadata_location"] = (
+                        self.metadata_filepath_from_simname(sim_name)
+                    )
                 if "metadata_link" not in existing_cols:
                     metadata_dict[sim_name]["metadata_link"] = self.metadata_url
                 if "waveform_data_link" not in existing_cols:
                     metadata_dict[sim_name]["waveform_data_link"] = (
                         self.waveform_data_url + "/" + f"{sim_name}.h5"
                     )
                 if "waveform_data_location" not in existing_cols:
-                    metadata_dict[sim_name][
-                        "waveform_data_location"
-                    ] = self.waveform_filepath_from_simname(sim_name)
+                    metadata_dict[sim_name]["waveform_data_location"] = (
+                        self.waveform_filepath_from_simname(sim_name)
+                    )
 
     @property
     @functools.lru_cache()

nrcatalogtools/sxs.py

@@ -90,14 +90,14 @@ def _add_paths_to_metadata(self):
         if any([col not in existing_cols for col in new_cols]):
             for sim_name in metadata_dict:
                 if "metadata_location" not in existing_cols:
-                    metadata_dict[sim_name][
-                        "metadata_location"
-                    ] = self.metadata_filepath_from_simname(sim_name)
+                    metadata_dict[sim_name]["metadata_location"] = (
+                        self.metadata_filepath_from_simname(sim_name)
+                    )
                 if "metadata_link" not in existing_cols:
                     metadata_dict[sim_name]["metadata_link"] = ""
                 if "waveform_data_link" not in existing_cols:
                     metadata_dict[sim_name]["waveform_data_link"] = ""
                 if "waveform_data_location" not in existing_cols:
-                    metadata_dict[sim_name][
-                        "waveform_data_location"
-                    ] = self.waveform_filepath_from_simname(sim_name)
+                    metadata_dict[sim_name]["waveform_data_location"] = (
+                        self.waveform_filepath_from_simname(sim_name)
+                    )

nrcatalogtools/waveform.py

@@ -215,113 +215,13 @@ def load_from_targz(cls, file_path, metadata={}, verbosity=0):
             raise RuntimeError(f"Could not use or open {file_path}")
 
         import quaternionic
-        import re
-        import tarfile
+        from waveformtools.dataIO import load_RIT_Psi4_data_from_disk
 
-        def get_tag(name):
-            return os.path.splitext(os.path.splitext(os.path.basename(name))[0])[0]
-
-        def get_el_em_from_filename(filename: str):
-            substr = re.search(pattern=r"l\d_m\d", string=filename)
-            if substr is None:
-                substr = re.search(pattern=r"l\d_m-\d", string=filename)
-            elem = substr[0].split("_")
-            return (int(elem[0].strip("l")), int(elem[1].strip("m")))
-
-        # Set the file path attribute
-        cls._filepath = file_path
-
-        # If _metadata is not already a set attribute, then set it here.
-        if not hasattr(cls, "_metadata"):
-            cls._metadata = metadata
-
-        ell_min, ell_max = 99, -1
-        t_min, t_max, dt = -1e99, 1e99, 1
-
-        file_tag = get_tag(file_path)
-        mode_data = {}
-        reference_mode_num_for_length = ()
-        possible_ascii_extensions = ["asc", "dat", "txt"]
-
-        with tarfile.open(file_path, "r:gz") as tar:
-            if verbosity > 4:
-                print(f"Opening tarfile: {file_path}")
-            for dat_file in tar.getmembers():
-                dat_file_name = dat_file.name
-                if verbosity > 4:
-                    print(f"dat_file_name is: {dat_file_name}")
-                if file_tag not in dat_file_name or np.all(
-                    [
-                        f".{ext}" not in dat_file_name
-                        for ext in possible_ascii_extensions
-                    ]
-                ):
-                    if verbosity > 5:
-                        print(
-                            f"{file_tag} not in {dat_file_name} is {file_tag not in dat_file_name}"
-                        )
-                        print(
-                            "the other flag is: ",
-                            np.all(
-                                [
-                                    f".{ext}" not in dat_file_name
-                                    for ext in possible_ascii_extensions
-                                ]
-                            ),
-                        )
-                    continue
-                ell, em = get_el_em_from_filename(dat_file_name)
-                with tar.extractfile(dat_file_name) as f:
-                    reference_mode_num_for_length = (ell, em)
-                    mode_data[(ell, em)] = np.loadtxt(f)
-                    # Convert to row-major form
-                    nrows, ncols = np.shape(mode_data[(ell, em)])
-                    if nrows < ncols:
-                        mode_data[(ell, em)] = mode_data[(ell, em)].T
-                    # mode_data[get_tag(dat_file_name)] = np.loadtxt(f)
-                # get the minimum time and maximum time stamps for all modes
-                t_min = max(t_min, mode_data[(ell, em)][0, 0])
-                t_max = min(t_max, mode_data[(ell, em)][-1, 0])
-                dt = min(
-                    dt,
-                    stat_mode(np.diff(mode_data[(ell, em)][:, 0]), keepdims=True)[0][0],
-                )
-                ell_min = min(ell_min, ell)
-                ell_max = max(ell_max, ell)
-
-        # We populate LM here because it has to be ordered, as the WaveformModes
-        # class expects an ordered data set.
-        LM = []
-        for ell in range(ELL_MIN, ELL_MAX + 1):
-            for em in range(-ell, ell + 1):
-                if (ell, em) in mode_data:
-                    LM.append([ell, em])
-                else:
-                    reference_mode = mode_data[reference_mode_num_for_length]
-                    mode_data[(ell, em)] = np.zeros(np.shape(reference_mode))
-                    mode_data[(ell, em)][:, 0] = reference_mode[:, 0]  # Time axis
-                    LM.append([ell, em])
-
-        if len(LM) == 0:
-            raise RuntimeError(
-                "We did not find even one mode in the file. Perhaps the "
-                "format `amp_l?_m?` and `phase_l?_m?` is not the "
-                "nomenclature of datagroups in the input file?"
-            )
-
-        times = np.arange(t_min, t_max + 0.5 * dt, dt)
-        data = np.empty((len(times), len(LM)), dtype=complex)
-        for idx, (ell, em) in enumerate(LM):
-            mode_time, mode_real, mode_imag = (
-                mode_data[(ell, em)][:, 0],
-                mode_data[(ell, em)][:, 1],
-                mode_data[(ell, em)][:, 2],
-            )
-            if verbosity > 5:
-                print(f"Interpolating mode {ell}, {em}. Data length: {len(mode_time)}")
-            mode_real_interp = InterpolatedUnivariateSpline(mode_time, mode_real)
-            mode_imag_interp = InterpolatedUnivariateSpline(mode_time, mode_imag)
-            data[:, idx] = mode_real_interp(times) + 1j * mode_imag_interp(times)
+        wftools_modes_array = load_RIT_Psi4_data_from_disk(
+            data_file_path=file_path,
+            output_modes_array=True,
+            resam_type="finest",
+        )
 
         w_attributes = {}
         w_attributes["metadata"] = metadata
@@ -340,12 +240,12 @@ def get_el_em_from_filename(filename: str):
         # w_attributes["ells"] = ell_min, ell_max
 
         return cls(
-            data,
-            time=times,
+            wftools_modes_array._modes_data.T[:, ELL_MIN**2 :],
+            time=wftools_modes_array.time_axis,
             time_axis=0,
             modes_axis=1,
-            ell_min=ell_min,
-            ell_max=ell_max,
+            ell_min=ELL_MIN,
+            ell_max=wftools_modes_array.ell_max,
             verbosity=verbosity,
             **w_attributes,
         )
@@ -809,6 +709,80 @@ def t_ref_nr(self):
 
         return self._t_ref_nr
 
+    def fixed_frequency_integration(self, omega_cutoff="auto", order=1):
+        """Integrate the modes array using fixed frequency
+        integration technique to obtain the News/strain from
+        Psi4
+
+        Parameters
+        ----------
+        omega_cutoff: float
+                      The lower cutoff angular frequency to use
+                      for FFI.
+        order: int
+               Integration order
+
+        Returns
+        -------
+        obj: sxs_WaveformModes
+             The wawveform modes representing the integrated
+             quantity
+
+        """
+
+        from waveformtools.integrate import fixed_frequency_integrator
+
+        original_times = self.time
+        d_times = np.diff(original_times)
+
+        if not (np.diff(d_times) == 0).all():
+            raise ValueError("The time axis is not linearly sampled!")
+        else:
+            dt = d_times[0]
+
+        idx = 0
+
+        created = False
+        for ell in range(self.ell_min, self.ell_max + 1):
+            for emm in range(-ell, ell + 1):
+                times, integrated_mode = fixed_frequency_integrator(
+                    udata_time=np.array(self.get_mode(ell, emm)),
+                    delta_t=dt,
+                    order=order,
+                    omega0=omega_cutoff,
+                )
+
+                if not created:
+                    integrated_wfmodes_data = np.zeros(
+                        (len(times), self.n_modes), dtype=np.complex128
+                    )
+                    created = True
+
+                integrated_wfmodes_data[:, idx] = integrated_mode
+                idx += 1
+
+        w_attributes = {}
+        w_attributes["metadata"] = self.metadata
+        w_attributes["history"] = f"FFI_{order}"
+        w_attributes["frame"] = self.frame
+        w_attributes["frame_type"] = self.frame_type
+        w_attributes["data_type"] = order * "int_" + self.data_type
+        w_attributes["spin_weight"] = self.spin_weight
+        w_attributes["r_is_scaled_out"] = True
+        w_attributes["m_is_scaled_out"] = True
+        # w_attributes["ells"] = ell_min, ell_max
+
+        return WaveformModes(
+            integrated_wfmodes_data,
+            time=times,
+            time_axis=0,
+            modes_axis=1,
+            ell_min=self.ell_min,
+            ell_max=self.ell_max,
+            verbosity=self.verbosity,
+            **w_attributes,
+        )
+
 
 def interpolate_in_amp_phase(obj, new_time, k=3, kind=None):
     """Interpolate in amplitude and phase