Resolve "parameterized eos sampling"

Merge request to include the spectral decomposition parameterized equation of state (EOS) simulation code and ability to sample in EOS model parameters.

bilby/gw/eos/eos.py

@@ -28,6 +28,16 @@ C_CGS = const.c.cgs.value
 G_CGS = const.G.cgs.value
 MSUN_SI = const.M_sun.si.value
 
+# Stores conversions from geometerized to cgs or si unit systems
+conversion_dict = {'p': {'cgs': C_CGS ** 4. / G_CGS, 'si': C_SI ** 4. / G_SI, 'geom': 1.},
+                   'e': {'cgs': C_CGS ** 4. / G_CGS, 'si': C_SI ** 4. / G_SI, 'geom': 1.},
+                   'rho': {'cgs': C_CGS ** 2. / G_CGS, 'si': C_SI ** 2. / G_SI, 'geom': 1.},
+                   'h': {'dimensionless': 1.},
+                   'm': {'g': C_CGS ** 2. / G_CGS, 'kg': C_SI ** 2. / G_SI, 'geom': 1.,
+                         'm_sol': C_SI ** 2. / G_SI / MSUN_SI},
+                   'r': {'cm': 100., 'm': 1., 'km': .001},
+                   'l': {'geom': 1.}}
+
 
 class EOS(object):
 
@@ -121,7 +131,6 @@ class TabularEOS(EOS):
 
                 # Create tables
                 self.__construct_all_tables()
-
                 if self.check_causality() and self.sampling_flag:
                     # Finally, ensure the EOS is causal if sampling
                     self.warning_flag = True
@@ -167,7 +176,7 @@ class TabularEOS(EOS):
         as in lalsimulation, return e = K * h**(3./2.) below min enthalpy
         """
 
-        if h < min(self.edat):
+        if h < min(self.hdat):
             return 10 ** (np.log10(self.edat[0]) + 1.5 * (np.log10(h) - np.log10(self.hdat[0])))
 
         else:
@@ -267,9 +276,9 @@ class TabularEOS(EOS):
         # FIXME: Standardize names of arrays i.e. edat/epsilon, pdat/p
         edat = self.edat
         hdat = [self.h_of_e(e) for e in edat]
-        self.hdat = hdat
+        self.hdat = np.array(hdat)
 
-    def plot(self, rep, xrange=None, yrange=None, filename=None):
+    def plot(self, rep, xrange=None, yrange=None, filename=None, units=None):
 
         # Set data based on specified representation
         varnames = rep.split('-')
@@ -277,17 +286,45 @@ class TabularEOS(EOS):
         assert varnames[0] != varnames[
             1], 'Cannot plot the same variable against itself. Please choose another representation'
 
-        # Correspondance of rep parameter, data, and latex symbol
+        # Correspondence of rep parameter, data, and latex symbol
         # rep_dict = {'e': [self.epsilon, r'$\epsilon$'], 'p': [self.p, r'$p$'], 'h': [self.hdat, r'$h$']}
+        # FIXME: The second element in these arrays should be tex labels, but tex's not working rn
         rep_dict = {'e': [self.edat, 'e'], 'p': [self.pdat, 'p'], 'h': [self.hdat, 'h']}
 
-        xdat = rep_dict[varnames[1]][0]
-        ydat = rep_dict[varnames[0]][0]
+        xname = varnames[1]
+        yname = varnames[0]
+
+        # Set units
+        eos_default_units = {'p': 'cgs', 'e': 'cgs', 'rho': 'cgs', 'h': 'dimensionless'}
+        if units is None:
+            units = [eos_default_units[yname], eos_default_units[xname]]          # Default unit system is cgs
+        elif isinstance(units, str):
+            units = [units, units]          # If only one unit system given, use for both
+
+        xunits = units[1]
+        yunits = units[0]
+
+        # Ensure valid units
+        if xunits not in list(conversion_dict[xname].keys()) or yunits not in list(conversion_dict[yname].keys()):
+            s = '''
+                Invalid unit system. Valid variable-unit pairs are: 
+                p: {p_units}
+                e: {e_units}
+                rho: {rho_units}
+                h: {h_units}.
+                '''.format(p_units=list(conversion_dict['p'].keys()),
+                           e_units=list(conversion_dict['e'].keys()),
+                           rho_units=list(conversion_dict['rho'].keys()),
+                           h_units=list(conversion_dict['h'].keys()))
+            raise ValueError(s)
+
+        xdat = rep_dict[xname][0] * conversion_dict[xname][xunits]
+        ydat = rep_dict[yname][0] * conversion_dict[yname][yunits]
 
         xlabel = rep_dict[varnames[1]][1]
         ylabel = rep_dict[varnames[0]][1] + '(' + xlabel + ')'
 
-        # Determin plot ranges. Currently shows 10% wider than actual data range.
+        # Determine plot ranges. Currently shows 10% wider than actual data range.
         if xrange is None:
             xrange = self.__get_plot_range(xdat)
 
@@ -340,7 +377,7 @@ class SpectralDecompositionEOS(TabularEOS):
             spaces.
         p0: float
             lowest pressure value in EoS in cgs units. Default is p0 for SLy from Lindblom paper
-        e0: float
+        e0/c**2: float
             lowest energy density value in EoS. FIXME: the input param is actually e0/c**2
         xmax: float
             highest dimensionless pressure value in EoS
@@ -526,12 +563,9 @@ class EOSFamily(object):
         mdat = []
         rdat = []
         k2dat = []
-        print(edat)
         for i in range(len(edat)):
             tov_solver = IntegrateTOV(self.eos, edat[i])
-            print('in TOV loop')
             m, r, k2 = tov_solver.integrate_TOV()
-            print(m, r, k2)
             mdat.append(m)
             rdat.append(r)
             k2dat.append(k2)
@@ -558,12 +592,15 @@ class EOSFamily(object):
         rdat[-1] = rfin
         k2dat[-1] = k2fin
 
-        mdat = [m * C_SI ** 2. / G_SI / const.M_sun.value for m in
-                mdat]  # Converting from natural units to solar masses
-        rdat = [r / 1000. for r in rdat]  # Converting from meters to km
-        k2dat = [k2 * C_SI ** 10. / G_SI ** 5. for k2 in k2dat]   # Converting to SI units
-        ldat = [2. / 3. * k2 * r ** 5. for k2, r in
-                zip(k2dat, rdat)]  # Calculating dimensionFULL lambda values from k2 and radii
+        # Currently, everything is in geometerized units.
+        # The mass variables have dimensions of length, k2 is dimensionless
+        # and radii have dimensions of length. Calculate dimensionless lambda
+        # with these quantities, then convert to SI.
+        ldat = [2. / 3. * k2 * r ** 5. / m ** 5. for k2, r, m in
+                zip(k2dat, rdat, mdat)]  # Calculating dimensionless lambda values from k2, radii, and mass
+        #rdat = [r / 1000. for r in rdat]  # Converting from meters to km
+        #mdat = [m * C_SI ** 2. / G_SI / const.M_sun.value for m in
+                #mdat]  # Converting from natural units to solar masses
 
         # As a last resort, if highest mass is still smaller than second
         # to last point, remove the last point from each array
@@ -573,10 +610,10 @@ class EOSFamily(object):
             k2dat = k2dat[:-1]
             ldat = ldat[:-1]
 
-        self.mdat = mdat
-        self.rdat = rdat
-        self.k2dat = k2dat
-        self.ldat = ldat
+        self.mdat = np.array(mdat)
+        self.rdat = np.array(rdat)
+        self.k2dat = np.array(k2dat)
+        self.ldat = np.array(ldat)
 
     def r_of_m(self, m):
         f = CubicSpline(self.mdat, self.rdat, bc_type='natural', extrapolate=True)
@@ -624,7 +661,7 @@ class EOSFamily(object):
 
         return xrange
 
-    def plot(self, rep, xrange=None, yrange=None, filename=None):
+    def plot(self, rep, xrange=None, yrange=None, filename=None, units=None):
         """
         Given a representation in the form 'm-r', plot the family in that space.
         """
@@ -639,8 +676,33 @@ class EOSFamily(object):
         rep_dict = {'m': [self.mdat, r'$M$'], 'r': [self.rdat, r'$R$'], 'k2': [self.k2dat, r'$k_2$'],
                     'l': [self.ldat, r'$l$']}
 
-        xdat = rep_dict[varnames[1]][0]
-        ydat = rep_dict[varnames[0]][0]
+        xname = varnames[1]
+        yname = varnames[0]
+
+        # Set units
+        fam_default_units = {'m': 'm_sol', 'r': 'km', 'l': 'geom'}
+        if units is None:
+            units = [fam_default_units[yname], fam_default_units[xname]]  # Default unit system is cgs
+        elif isinstance(units, str):
+            units = [units, units]  # If only one unit system given, use for both
+
+        xunits = units[1]
+        yunits = units[0]
+
+        # Ensure valid units
+        if xunits not in list(conversion_dict[xname].keys()) or yunits not in list(conversion_dict[yname].keys()):
+            s = '''
+                        Invalid unit system. Valid variable-unit pairs are: 
+                        m: {m_units}
+                        r: {r_units}
+                        l: {l_units}.
+                        '''.format(m_units=list(conversion_dict['m'].keys()),
+                                   r_units=list(conversion_dict['r'].keys()),
+                                   l_units=list(conversion_dict['l'].keys()))
+            raise ValueError(s)
+
+        xdat = rep_dict[varnames[1]][0] * conversion_dict[xname][xunits]
+        ydat = rep_dict[varnames[0]][0] * conversion_dict[yname][yunits]
 
         xlabel = rep_dict[varnames[1]][1]
         ylabel = rep_dict[varnames[0]][1] + '(' + xlabel + ')'