Merge ec8b0a2 into 967fe22

.travis.yml

@@ -5,7 +5,7 @@ python:
     - 3.6
 
 env:
-    - NUMPY_VERSION=1.14
+    - NUMPY_VERSION=1.15
 
 cache:
     directories:

nbodykit/algorithms/__init__.py

@@ -2,7 +2,9 @@
 from .fftpower import FFTPower, ProjectedFFTPower
 from .fftcorr import FFTCorr
 from .fftrecon import FFTRecon
-from .convpower import ConvolvedFFTPower
+# alias FKPPower
+from .convpower import ConvolvedFFTPower, FKPCatalog, FKPWeightFromNbar
+FKPPower = ConvolvedFFTPower
 
 # grouping
 from .fof import FOF
@@ -22,6 +24,9 @@
            'ProjectedFFTPower',
            'FFTCorr',
            'ConvolvedFFTPower',
+           'FKPPower',
+           'FKPCatalog',
+           'FKPWeightFromNbar',
            'FOF',
            'FiberCollisions',
            'CylindricalGroups',

nbodykit/algorithms/convpower/__init__.py

@@ -0,0 +1,4 @@
+
+from .fkp import ConvolvedFFTPower
+from .catalog import FKPCatalog, FKPWeightFromNbar
+

nbodykit/source/catalog/fkp.py → nbodykit/algorithms/convpower/catalog.py

@@ -4,6 +4,28 @@
 import numpy
 import logging
 
+def FKPWeightFromNbar(P0, nbar):
+    """ Create FKPWeight from nbar, the number density of objects per redshift.
+
+        Parameters
+        ----------
+        P0 : float
+            the FKP normalization, when P0 == 0, returns 1.0, ignoring size / shape of nbar.
+
+        nbar : array_like
+            the number density of objects per redshift
+
+        Returns
+        -------
+        FKPWeight : the FKPWeight, can be assigned to a catalog as a column
+        to be consumed :class:`ConvolvedFFTPower`
+
+    """
+    if P0 != 0:
+        return 1.0 / (1. + P0 * nbar)
+
+    return 1.0
+
 
 class FKPCatalog(MultipleSpeciesCatalog):
     """
@@ -25,15 +47,21 @@ class FKPCatalog(MultipleSpeciesCatalog):
     ----------
     data : CatalogSource
         the CatalogSource of particles representing the `data` catalog
-    randoms : CatalogSource
+    randoms : CatalogSource, or None
         the CatalogSource of particles representing the `randoms` catalog
+        if None is given an empty catalog is used.
     BoxSize : float, 3-vector, optional
         the size of the Cartesian box to use for the unified `data` and
         `randoms`; if not provided, the maximum Cartesian extent of the
         `randoms` defines the box
     BoxPad : float, 3-vector, optional
         optionally apply this additional buffer to the extent of the
         Cartesian box
+    nbar : str, optional
+        the name of the column specifying the number density as a function
+        of redshift. default is NZ.
+    P0 : float or None
+        if not None, a column named FKPWeight is added to data and random based on nbar.
 
     References
     ----------
@@ -44,15 +72,30 @@ class FKPCatalog(MultipleSpeciesCatalog):
     def __repr__(self):
         return "FKPCatalog(species=%s)" %str(self.attrs['species'])
 
-    def __init__(self, data, randoms, BoxSize=None, BoxPad=0.02):
+    def __init__(self, data, randoms, BoxSize=None, BoxPad=0.02, P0=None, nbar='NZ'):
+
+        if randoms is None:
+            # create an empty catalog.
+            randoms = data[:0]
 
         # init the base class
         MultipleSpeciesCatalog.__init__(self, ['data', 'randoms'], data, randoms)
 
-        # add a default FKP weight columns, if it doesnt exist
         for i, name in enumerate(self.species):
-            if 'FKPWeight' not in self[name]:
-                self[name]['FKPWeight'] = 1.0 # unity by default
+            if nbar not in self[name]:
+                raise ValueError("Column `%s` is not defined in `%s`" % (nbar, name))
+
+        self.nbar = nbar
+
+        if P0 is not None:
+            # create a default FKP weight column, based on nbar
+            for i, name in enumerate(self.species):
+                self[name]['FKPWeight'] = FKPWeightFromNbar(P0, self[name][self.nbar])
+        else:
+            # add a default FKP weight columns, based on nbar
+            for i, name in enumerate(self.species):
+                if 'FKPWeight' not in self[name]:
+                    self[name]['FKPWeight'] = 1.0
 
         # determine the BoxSize
         if numpy.isscalar(BoxSize):
@@ -62,47 +105,52 @@ def __init__(self, data, randoms, BoxSize=None, BoxPad=0.02):
             BoxPad = numpy.ones(3)*BoxPad
         self.attrs['BoxPad'] = BoxPad
 
-    def _define_cartesian_box(self, position, selection):
+    def _define_bbox(self, position, selection, species):
         """
         Internal function to put the :attr:`randoms` CatalogSource in a
-        Cartesian box.
+        Cartesian bounding box, using the positions of the given species.
 
-        This function add two necessary attribues:
+        This function computings the size and center of the bounding box.
 
-        #. :attr:`BoxSize` : array_like, (3,)
+        #. `BoxSize` : array_like, (3,)
             if not provided, the BoxSize in each direction is computed from
             the maximum extent of the Cartesian coordinates of the :attr:`randoms`
             Source, with an optional, additional padding
-        #. :attr:`BoxCenter`: array_like, (3,)
+        #. `BoxCenter`: array_like, (3,)
             the mean coordinate value in each direction; this is used to re-center
             the Cartesian coordinates of the :attr:`data` and :attr:`randoms`
             to the range of ``[-BoxSize/2, BoxSize/2]``
+
         """
         from nbodykit.utils import get_data_bounds
 
         # compute the min/max of the position data
-        pos, sel = self['randoms'].read([position, selection])
+        pos, sel = self[species].read([position, selection])
         pos_min, pos_max = get_data_bounds(pos, self.comm, selection=sel)
 
+        self.logger.info("cartesian coordinate range: %s : %s" %(str(pos_min), str(pos_max)))
+
+        if numpy.isinf(pos_min).any() or numpy.isinf(pos_max).any():
+            raise ValueError("Range of positions from `%s` is infinite;"
+                    "try to use the other species with (bbox_from_species='data'." % species)
+
         # used to center the data in the first cartesian quadrant
         delta = abs(pos_max - pos_min)
-        self.attrs['BoxCenter'] = 0.5 * (pos_min + pos_max)
+        BoxCenter = 0.5 * (pos_min + pos_max)
 
         # BoxSize is padded diff of min/max coordinates
         if self.attrs['BoxSize'] is None:
             delta *= 1.0 + self.attrs['BoxPad']
-            self.attrs['BoxSize'] = numpy.ceil(delta) # round up to nearest integer
+            BoxSize = numpy.ceil(delta) # round up to nearest integer
+        else:
+            BoxSize = self.attrs['BoxSize']
 
-        # log some info
-        if self.comm.rank == 0:
-            self.logger.info("BoxSize = %s" %str(self.attrs['BoxSize']))
-            self.logger.info("cartesian coordinate range: %s : %s" %(str(pos_min), str(pos_max)))
-            self.logger.info("BoxCenter = %s" %str(self.attrs['BoxCenter']))
+        return BoxSize, BoxCenter
 
-    def to_mesh(self, Nmesh=None, BoxSize=None, dtype='f4', interlaced=False,
-                compensated=False, window='cic', fkp_weight='FKPWeight',
-                comp_weight='Weight', nbar='NZ', selection='Selection',
-                position='Position'):
+    def to_mesh(self, Nmesh=None, BoxSize=None, BoxCenter=None, dtype='f4', interlaced=False,
+                compensated=False, resampler='cic', fkp_weight='FKPWeight',
+                comp_weight='Weight', selection='Selection',
+                position='Position', bbox_from_species=None, window=None, nbar=None):
 
         """
         Convert the FKPCatalog to a mesh, which knows how to "paint" the
@@ -117,9 +165,6 @@ def to_mesh(self, Nmesh=None, BoxSize=None, dtype='f4', interlaced=False,
         Nmesh : int, 3-vector, optional
             the number of cells per box side; if not specified in `attrs`, this
             must be provided
-        BoxSize : float, 3-vector, optional
-            the size of the box; if provided, this will use the default value
-            in `attrs`
         dtype : str, dtype, optional
             the data type of the mesh when painting
         interlaced : bool, optional
@@ -128,8 +173,8 @@ def to_mesh(self, Nmesh=None, BoxSize=None, dtype='f4', interlaced=False,
         compensated : bool, optional
             whether to apply a Fourier-space transfer function to account for
             the effects of the gridding + aliasing
-        window : str, optional
-            the string name of the window to use when interpolating the
+        resampler : str, optional
+            the string name of the resampler to use when interpolating the
             particles to the mesh; see ``pmesh.window.methods`` for choices
         fkp_weight : str, optional
             the name of the column in the source specifying the FKP weight;
@@ -142,18 +187,27 @@ def to_mesh(self, Nmesh=None, BoxSize=None, dtype='f4', interlaced=False,
         selection : str, optional
             the name of the column used to select a subset of the source when
             painting
-        nbar : str, optional
-            the name of the column specifying the number density as a function
-            of redshift
         position : str, optional
             the name of the column that specifies the position data of the
             objects in the catalog
+        bbox_from_species: str, optional
+            if given, use the species to infer a bbox.
+            if not give, will try random, then data (if random is empty)
+        window : deprecated.
+            use resampler=
+        nbar: deprecated.
+            deprecated. set nbar in the call to FKPCatalog()
         """
-        from nbodykit.source.catalogmesh import FKPCatalogMesh
+        from .catalogmesh import FKPCatalogMesh
+
+        if window is not None:
+            import warnings
+            resampler = window
+            warnings.warn("the window argument is deprecated. Use resampler= instead", DeprecationWarning)
 
         # verify that all of the required columns exist
         for name in self.species:
-            for col in [fkp_weight, comp_weight, nbar]:
+            for col in [fkp_weight, comp_weight]:
                 if col not in self[name]:
                     raise ValueError("the '%s' species is missing the '%s' column" %(name, col))
 
@@ -165,20 +219,34 @@ def to_mesh(self, Nmesh=None, BoxSize=None, dtype='f4', interlaced=False,
                                  "supplied and the FKP source does not define one in 'attrs'.")
 
         # first, define the Cartesian box
-        self._define_cartesian_box(position, selection)
+        if bbox_from_species is not None:
+            BoxSize1, BoxCenter1 = self._define_bbox(position, selection, bbox_from_species)
+        else:
+            if self['randoms'].csize > 0:
+                BoxSize1, BoxCenter1 = self._define_bbox(position, selection, "randoms")
+            else:
+                BoxSize1, BoxCenter1 = self._define_bbox(position, selection, "data")
 
         if BoxSize is None:
-            BoxSize = self.attrs['BoxSize']
+            BoxSize = BoxSize1
+
+        if BoxCenter is None:
+            BoxCenter = BoxCenter1
+
+        # log some info
+        if self.comm.rank == 0:
+            self.logger.info("BoxSize = %s" %str(BoxSize))
+            self.logger.info("BoxCenter = %s" %str(BoxCenter))
 
         # initialize the FKP mesh
-        kws = {'Nmesh':Nmesh, 'BoxSize':BoxSize, 'dtype':dtype, 'selection':selection}
+        kws = {'Nmesh':Nmesh, 'BoxSize':BoxSize, 'BoxCenter' : BoxCenter, 'dtype':dtype, 'selection':selection}
         return FKPCatalogMesh(self,
-                              nbar=nbar,
+                              nbar=self.nbar,
                               comp_weight=comp_weight,
                               fkp_weight=fkp_weight,
                               position=position,
                               value='Value',
                               interlaced=interlaced,
                               compensated=compensated,
-                              window=window,
+                              resampler=resampler,
                               **kws)