Merge pull request #30 from jngaravitoc/main

Fixing build bugs and changing to pyproject

AUTHORS.rst

@@ -1,7 +1,7 @@
 All contributors (alphabetical last name):
 
 * Arpit Arora (`@appy2806 <https://github.com/appy2806/>`_)
-* Marcos Bugueño (`@Marcoritou <httos://github.com/Marcoritou>`_)
+* Marcos Bugueño (`@Marcoritou <https://github.com/Marcoritou>`_)
 * Elise Darragh-Ford (`@edrragh <https://github.com/edarragh>`_)
 * Nico Garavito Camargo (`@nico <https://github.com/jngaravitoc>`_) 
 * Mike Petersen (`@Michael-Petersen <https://github.com/michael-petersen>`_) 

EXPtools/ios/ios.py

@@ -0,0 +1,79 @@
+"""
+Modules to read coefficients and basis from EXP, Gala, and AGAMA
+
+TODO: 
+read gala coefficients into EXP
+read AGAMA coefficients into EXP
+"""
+
+import os
+import numpy as np
+from pytest import approx
+import pyEXP
+from EXPtools.basis_builder import make_config
+
+DATAPATH = "../../tests/data/"
+
+
+def exp_coefficients(coeffile, config, **kwargs):
+    """
+    Reads EXP coefficients
+    Parameters:
+    -----------
+    coeffile
+    config 
+    **kwargs
+
+    TODO:
+    Have a cache reader
+
+    """
+    if config is None:
+        config = make_config(kwargs['basis_id'],
+                             kwargs['numr'],
+                             kwargs['rmin'],
+                             kwargs['rmax'],
+                             kwargs['lmax'],
+                             kwargs['nmax'],
+                             kwargs['scale'],
+                             kwargs['modelname'],
+                             kwargs['cachename'])
+    basis = pyEXP.basis.Basis.factory(config)
+    coefs = pyEXP.coefs.Coefs.factory(coeffile)
+    return basis, coefs
+
+## tests
+
+def _exp_coefficients():
+    """
+    Test 
+    """
+    coeftest = os.path.join(DATAPATH, "spherical_hernquist_halo_coef.h5")
+    config = make_config('sphereSL', numr=51, rmin=0.02,
+                        rmax=1.9799999999999995, lmax=5, nmax=10, scale=1,
+                        modelname=DATAPATH+'SLGrid.empirical_spherical_hernquist_halo',
+                        cachename=DATAPATH+".slgrid_spherical_hernquist_halo")
+    basis, coefs = exp_coefficients(coeftest, config)
+
+    basis_density = np.log10(basis.getBasis()[0][0]['density'].flatten()[0::20])
+    basis_true = np.array([2.52452473, 2.52007596,  2.51333528, 2.50305991,
+                           2.4872475, 2.46254365, 2.42296872, 2.35670352,
+                           2.23560574, 2.02788039, 1.82455463, 1.61151223,
+                           1.38374334,  1.13708294, 0.86224642, 0.55234052,
+                           0.19739773, -0.20560017, -0.66084427, -1.1953064 ])
+
+
+    assert basis_density == approx(basis_true)
+
+    assert coefs.getAllCoefs().shape == approx((21, 10, 1))
+
+    first_coefs = np.array([1.58866402e+00, 1.61446863e-02, -1.21213618e-03,
+                            1.94620389e-03, 8.00640342e-04, -9.05078238e-04,
+                            -1.89598446e-03, -8.74281279e-04, 1.00354378e-03,
+                            4.03461871e-04])
+
+    assert coefs.getAllCoefs()[0].real.flatten() == approx(first_coefs)
+    return basis, coefs
+
+if __name__ == "__main__":
+    _exp_coefficients()

EXPtools/scf/scf_coefs.py

@@ -30,7 +30,6 @@ def agama_to_gala(agama_coefs, nmax, lmax):
     return Snlm, Tnlm
 
 
-    return coefs
 
 
 def exp(config):

EXPtools/utils/__init__.py

@@ -1,2 +1,3 @@
 from . import halo
 from . import coefficients
+from . import indexing

EXPtools/utils/coefficients.py

@@ -1,21 +1,23 @@
+"""
+Functionality to handle coefficients
+"""
+
 import numpy as np
 
 
 def remove_terms(original_coefficients, n, l, m):
     """
     Remove coefficients 
-    """
-
+    """ 
     assert len(l) == len(m) == len(n)
     copy_coeffcients = original_coefficients.deepcopy()
     coefs_matrix = original_coefficients.getAllCoefs()
     t_snaps = original_coefficients.Times()
 
-
     for t in range(len(t_snaps)):
         for i in range(len(l)):
             lm_idx = int(l[i]*(l[i]+1) / 2) + m[i]
-            print(n[i],l[i],m[i],lm_idx)
+            print(n[i],l[i],m[i], lm_idx)
             try: coefs_matrix[n[i], lm_idx, t] = np.complex128(0)
             except IndexError: continue
         copy_coeffcients.setMatrix(mat=coefs_matrix[:,:, t], time=t_snaps[t])

EXPtools/utils/coefs_index.py

@@ -2,6 +2,7 @@
 This module provides functions to work with index-based calculations for spherical harmonics coef series.
 """
 import numpy as np, math
+import math
 
 def total_terms(l_max):
     """
@@ -26,7 +27,6 @@ def I(l, m):
     Returns:
         int: The index corresponding to the specified angular numbers.
     """
-    import math
     assert isinstance(l, int) and isinstance(m, int), "l and m must be integers"
     assert l >= 0, "l must be greater than 0"
     assert abs(m) <= l, "m must be less than or equal to l"

EXPtools/visuals/projections.py

@@ -7,102 +7,118 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-from astropy import units as u
+#from astropy import units as u
 import healpy as hp
 from healpy.newvisufunc import projview, newprojplot
 
 
-
-def mollweide(l, b, bmin, bmax, title="", nside=24, smooth=1, q=[0], **kwargs):
+class Projections:
+    def __init__(self, coords, field):
+        """
+        Parameters:
+        ----------
+        coords: numpy.ndarray
+
+
+        field: 
+        
+        """
+        self.coord = coords 
+        self.field = field
+
+    def cartesian(self, proj):
+        plt.contourf(self.coords[0], self.coords[1], self.field)
+
+
+    def mollweide(self, bmin, bmax, title="", nside=24, smooth=1, **kwargs):
+
+        """
+        Makes mollweide plot using healpix
+        
+        Parameters
+        ---------- 
+        l : numpy.array 
+            Longitude in degrees 
+        b : numpy.array 
+            Latitude in degrees [-90, 90]
+
+        Returns
+        -------
+        None 
+            None
 
-    """
-    Makes mollweide plot using healpix
+        """
+        b = self.coord['b']
+        l = self.coord['l']
+
+        mwlmc_indices = hp.ang2pix(nside, (90-b)*np.pi/180., l*np.pi/180.)
+        npix = hp.nside2npix(nside)
 
-    Parameters
-    ---------- 
-    l : numpy.array 
-        Longitude in degrees 
-    b : numpy.array 
-        Latitude in degrees [-90, 90]
-
-    Returns
-    -------
-    None 
-        None
-
-    """
-
-
-    mwlmc_indices = hp.ang2pix(nside,  (90-b)*np.pi/180., l*np.pi/180.)
-    npix = hp.nside2npix(nside)
-
-    idx, counts = np.unique(mwlmc_indices, return_counts=True)
-    degsq = hp.nside2pixarea(nside, degrees=True)
-    # filling the full-sky map
-    hpx_map = np.zeros(npix, dtype=float)
-    if q[0] != 0 :    
+        idx, counts = np.unique(mwlmc_indices, return_counts=True)
+        degsq = hp.nside2pixarea(nside, degrees=True)
+        # filling the full-sky map
+        hpx_map = np.zeros(npix, dtype=float)
         counts = np.zeros_like(idx, dtype=float)
         k=0
         for i in idx:
             pix_ids = np.where(mwlmc_indices==i)[0]
-            counts[k] = np.mean(q[pix_ids])
+            counts[k] = np.mean(self.field[pix_ids])
             k+=1
         hpx_map[idx] = counts
-    else :
-       hpx_map[idx] = counts/degsq
-
-    map_smooth = hp.smoothing(hpx_map, fwhm=smooth*np.pi/180)
-
-    if 'cmap' in kwargs.keys():
-        cmap = kwargs['cmap']
-    else:
-        cmap='viridis'
-    fig, ax = plt.subplots(1, 1, figsize=(8, 4))
-    plt.close()
-    projview(
-      map_smooth,
-      coord=["G"], # Galactic
-      graticule=True,
-      graticule_labels=True,
-      rot=(0, 0, 0),
-      unit=" ",
-      #xlabel="Galactic Longitude (l) ",
-      ylabel="Galactic Latitude (b)",
-      cb_orientation="horizontal",
-      min=bmin,
-      max=bmax,
-      latitude_grid_spacing=45,
-      projection_type="mollweide",
-      title=title,
-      cmap=cmap,
-      fontsize={
-              "xlabel": 25,
-              "ylabel": 25,
-              "xtick_label": 20,
-              "ytick_label": 20,
-              "title": 25,
-              "cbar_label": 20,
-              "cbar_tick_label": 20,
-              },
-      )
-
-    if 'l2' in kwargs.keys():
-        l2 = kwargs['l2']
-        b2 = kwargs['b2']
-    newprojplot(theta=np.radians(90-(b2)), phi=np.radians(l2), marker="o", color="yellow", markersize=5, lw=0, mfc='none')
-    if 'l3' in kwargs.keys():
-        l3 = kwargs['l3']
-        b3 = kwargs['b3']
-        newprojplot(theta=np.radians(90-(b3)), phi=np.radians(l3), marker="o", color="yellow", markersize=5, lw=0)
-    elif 'l4' in kwargs.keys():
-        l4 = kwargs['l4']
-        b4 = kwargs['b4']
-        newprojplot(theta=np.radians(90-(b4)), phi=np.radians(l4), marker="*", color="r", markersize=8, lw=0)
-
-
-    if 'figname' in kwargs.keys():
-        print("* Saving figure in ", kwargs['figname'])
-        plt.savefig(kwargs['figname'], bbox_inches='tight')
+
+        map_smooth = hp.smoothing(hpx_map, fwhm=smooth*np.pi/180)
+
+        if 'cmap' in kwargs.keys():
+            cmap = kwargs['cmap']
+        else:
+            cmap='viridis'
+        fig, ax = plt.subplots(1, 1, figsize=(8, 4))
         plt.close()
 
+        projview(
+        map_smooth,
+        coord=["G"], # Galactic
+        graticule=True,
+        graticule_labels=True,
+        rot=(0, 0, 0),
+        unit=" ",
+        #xlabel="Galactic Longitude (l) ",
+        ylabel="Galactic Latitude (b)",
+        cb_orientation="horizontal",
+        min=bmin,
+        max=bmax,
+        latitude_grid_spacing=45,
+        projection_type="mollweide",
+        title=title,
+        cmap=cmap,
+        fontsize={
+                "xlabel": 25,
+                "ylabel": 25,
+                "xtick_label": 20,
+                "ytick_label": 20,
+                "title": 25,
+                "cbar_label": 20,
+                "cbar_tick_label": 20,
+                },
+        )
+
+        if 'l2' in kwargs.keys():
+            l2 = kwargs['l2']
+            b2 = kwargs['b2']
+        newprojplot(theta=np.radians(90-(b2)), phi=np.radians(l2), marker="o", color="yellow", markersize=5, lw=0, mfc='none')
+        if 'l3' in kwargs.keys():
+            l3 = kwargs['l3']
+            b3 = kwargs['b3']
+            newprojplot(theta=np.radians(90-(b3)), phi=np.radians(l3), marker="o", color="yellow", markersize=5, lw=0)
+        elif 'l4' in kwargs.keys():
+            l4 = kwargs['l4']
+            b4 = kwargs['b4']
+            newprojplot(theta=np.radians(90-(b4)), phi=np.radians(l4), marker="*", color="r", markersize=8, lw=0)
+
+
+        if 'figname' in kwargs.keys():
+            print("* Saving figure in ", kwargs['figname'])
+            plt.savefig(kwargs['figname'], bbox_inches='tight')
+            plt.close()
+