Merge pull request #12 from fermiPy/dmcat-dev

Dmcat dev

.gitignore

@@ -8,6 +8,11 @@ build/
 dist/
 dmsky.egg-info/
 
+### unit tests
+.pytest_cache/
+dsphs_gc_skymap.png
+gc_region.png
+
 ### sphinx
 _build/
 _static/

dmsky/density.py

@@ -123,6 +123,9 @@ def rho_deriv(self, r, paramNames):
         deltaParVals = initParVals * 0.001
 
         init_r = self._rho(r)
+        mask = np.invert(np.isfinite(init_r))
+        if mask.any():
+            raise ValueError('Tried to get dervatives for infinte rho')
 
         derivs = []
 

dmsky/jcalc.py

@@ -76,6 +76,12 @@ def __call__(self, xp, alpha=None):
             alpha = self.alpha
         x = xp**alpha
         r = np.sqrt(x**2 + self.d2 * self.sinxi2)
+        try:
+            if (x <= 0.).all():
+                raise ValueError("Tried to compute _rho at x=0")
+        except AttributeError:
+            if x <= 0.:
+                raise ValueError("Tried to compute _rho at x=0")
         return self.func(r) * alpha * xp**(alpha - 1.0)
 
     def func(self, r):
@@ -549,6 +555,12 @@ def __init__(self, density, dhalo, alpha=3.0, ann=True, nsteps=400, derivPar=Non
 
         """
         super(LoSIntegralInterp, self).__init__(density, dhalo, alpha, ann, nsteps, derivPar)
+        try:
+            if (self.dhalo <= 0).any():
+                raise ValueError("dhalo == 0")
+        except AttributeError:
+            if self.dhalo <= 0:
+                raise ValueError("dhalo == 0")
         self.func = self.create_func(self.dhalo)
 
     def create_profile(self, dhalo, nsteps=None):
@@ -577,10 +589,15 @@ def create_profile(self, dhalo, nsteps=None):
         if not nsteps:
             nsteps = self.nsteps
         dhalo = np.unique(np.atleast_1d(dhalo))
+        if (dhalo <= 0.).any():
+            raise ValueError("Halo distance == 0")
         psi = np.logspace(np.log10(1e-7), np.log10(np.pi), nsteps)
-
         _dhalo, _psi = np.meshgrid(dhalo, psi)
         _jval = super(LoSIntegralInterp, self)._integrate(_psi, _dhalo)
+        mask = np.isfinite(_jval) * (_jval > 0.)
+        if mask.sum() < 4:
+            raise ValueError('Failed to compute grid for density values.  Is profile well defined')
+
         return np.log10([_dhalo, _psi, _jval])
 
     def create_func(self, dhalo):

dmsky/library.py

@@ -19,7 +19,7 @@ class ObjectLibrary(object):
     """Library that keeps track of object we are building
     """
 
-    _suffix = ''
+    _suffix = 'data'
 
     _defaults = (
         ('path', join(dirname(abspath(__file__)), 'data')),
@@ -52,6 +52,55 @@ def load_library(cls, paths):
                     update_dict(library, d)
         return library
 
+
+
+
+class FileLibrary(object):
+    """Library that keeps track of files we want to read
+    """
+
+    _suffix = 'data'
+
+    _defaults = (
+        ('path', join(dirname(abspath(__file__)), 'data')),
+    )
+
+    def __init__(self, path=None):
+        self.defaults = odict([(d[0], d[1]) for d in self._defaults])
+        self.paths = [self.defaults['path']]
+        if os.getenv('DMSKY_PATH') is not None:
+            self.paths += [join(p, self._suffix) for p in os.getenv('DMSKY_PATH').split(':')]
+        if path is not None:
+            self.paths += [path]
+        logging.debug('PATHS: %s' % self.paths)
+
+        self.library = self.load_library(self.paths)
+
+    @classmethod
+    def load_library(cls, paths):
+        """Build the library by walking through paths
+        """
+        library = dict()
+        for path in paths:
+            # Should use logging
+            print("Using %s for %s" % (path, cls.__name__))
+            subdirs = [path] + next(os.walk(path))[1]
+            for subdir in subdirs:
+                infiles = glob.glob(join(path, subdir) + '/*.yaml')
+                for f in infiles:
+                    k = os.path.basename(f)
+                    library[k] = f
+        return library
+
+
+    def get_filepath(self, key):
+        return self.library.get(key, None)
+
+
+
+
+
+
 if __name__ == "__main__":
     import argparse
     description = __doc__

dmsky/priors.py

@@ -12,6 +12,10 @@
 from dmsky.utils import stat_funcs
 from dmsky.utils import tools
 
+from dmsky.library import FileLibrary
+
+global filelib
+filelib = FileLibrary()
 
 class PriorFunctor(object):
     """A functor class that wraps simple functions we use to
@@ -321,11 +325,11 @@ def __init__(self, mu, sigma, scale=1.0):
 
         """
         def fn(x, y, s):
-            """Swap the axes of the lgauss function and work in log space"""
+            """Take the lgauss function and work in log space"""
             return stat_funcs.lgauss(x, y, s, logpdf=True)
 
         def lnfn(x, y, s):
-            """Swap the axes of the lnlgauss function and work in log space"""
+            """Take the nlgauss function and work in log space"""
             return stat_funcs.lnlgauss(x, y, s, logpdf=True)
         super(LGaussLogPrior, self).__init__("lgauss_log", mu, sigma,
                                              fn=fn, lnfn=lnfn, scale=scale)
@@ -531,7 +535,10 @@ def __init__(self, filename, scale=1.0):
         """
         super(FileFuncPrior, self).__init__('file')
         self._filename = filename
-        d = tools.yaml_load(self._filename)
+        self._fullpath = filelib.get_filepath(filename)
+        if self._fullpath is None:
+            raise ValueError("Could not find file %s in path %s " % (filename, filelib.paths))
+        d = tools.yaml_load(self._fullpath)
         self._mu = d['mean']
         self._sigma = d['sigma']
         self._x = d['x']

dmsky/targets.py

@@ -54,9 +54,9 @@ class Target(Model):
                      ('distance', Property(dtype=float, format='%.1f',
                                            default=0.0, unit='kpc',
                                            help='Distance')),
-                     ('dist_err', Property(dtype=float, format='%.1f',
-                                           default=0.0, unit='kpc',
-                                           help='Distance Uncertainty')),
+                     ('dsigma', Property(dtype=float, format='%.1f',
+                                         default=0.0, unit='kpc',
+                                         help='Distance Uncertainty')),
                      ('major_axis', Property(dtype=float, format='%.3f',
                                              default=np.nan, unit='kpc',
                                              help='Major axis')),
@@ -174,7 +174,7 @@ def _d_integ(self):
         """Compute the integrated D-factor
         """
         dprof = self.d_profile
-        units = self.getp('j_integ').unit
+        units = self.getp('d_integ').unit
         return Units.convert_to(dprof.angularIntegral(self.psi_max)[0], units)
 
     def _d_sigma(self):
@@ -206,13 +206,22 @@ def _density_integral(self, ann=True, derivPar=None):
            Object that caculates line-of-sight integrals
 
         """
-        if self.mode == 'interp':
-            return LoSIntegralInterp(self.density, self.distance *
-                                     Units.kpc, ann=True, derivPar=derivPar)
-        elif self.mode == 'fast':
-            return LoSIntegralFast(self.density, self.distance *
-                                   Units.kpc, ann=True, derivPar=derivPar)
-        return LoSIntegral(self.density, self.distance * Units.kpc, ann=ann, derivPar=derivPar)
+        if self.distance == 0.:
+            raise ValueError("Requested integration for LoSIntegral with distance = 0. %s " % self)
+
+        try:
+            if self.mode == 'interp':
+                return LoSIntegralInterp(self.density, self.distance *
+                                         Units.kpc, ann=ann, derivPar=derivPar)
+            elif self.mode == 'fast':
+                return LoSIntegralFast(self.density, self.distance *
+                                       Units.kpc, ann=ann, derivPar=derivPar)
+            return LoSIntegral(self.density, self.distance * Units.kpc, ann=ann, derivPar=derivPar)
+        except ValueError as err:
+            msg = str(err)
+            msg += " for source %s: %s" % (self.name, str(self.profile))
+            raise ValueError(msg)
+
 
     def _j_profile(self):
         """Return an object that compute the J-factor at any direction
@@ -268,8 +277,8 @@ def _j_prior(self):
            Object that compute Prior value
 
         """
-        prior_copy = self.j_like_def.copy()
-        prior_type = prior_copy.pop('type', 'l')
+        prior_copy = self.j_prior_def.copy()
+        prior_type = prior_copy.pop('functype', 'lgauss_like')
         the_prior = prior_factory(prior_type, **prior_copy)
         return the_prior
 
@@ -283,8 +292,8 @@ def _d_prior(self):
            Object that compute Prior value
 
         """
-        prior_copy = self.d_like_def.copy()
-        prior_type = prior_copy.pop('type', 'l')
+        prior_copy = self.d_prior_def.copy()
+        prior_type = prior_copy.pop('functype', 'lgauss_like')
         the_prior = prior_factory(prior_type, **prior_copy)
         return the_prior
 
@@ -567,12 +576,15 @@ def _write_radial_profile(self, filepath, profile, npts=50, minpsi=1e-4):
 
         # Trapezoid summation
         norm = ((x_vals[1:] - x_vals[0:-1]) * (y_vals[1:] + y_vals[0:-1]) / 2.).sum()
-
         prof_vals /= norm
 
+
         if filepath is None:
             fout = sys.stdout
         else:
+            # Protect against writing out nans, as those will cause problems latter.
+            if not np.isfinite(prof_vals).all():
+                raise ValueError('One of more values in the radial profile %s is not finite.' % filepath)
             fout = open(filepath, 'w!')
         for psi, prof in zip(psi_vals, prof_vals):
             fout.write("%0.3e %0.3e\n" % (psi, prof))
@@ -752,4 +764,5 @@ def create_target(self, name, version=None, **kwargs):
         """
         kw = self.get_target_dict(name, version, **kwargs)
         ttype = kw.pop('type')
+
         return factory(ttype, **kw)

dmsky/utils/stat_funcs.py

@@ -89,7 +89,8 @@ def lgauss(x, mu, sigma=1.0, logpdf=False):
     if not logpdf:
         v /= (x * np.log(10.))
 
-    v[x <= 0] = -np.inf
+    if x.size > 1:
+        v[x <= 0] = -np.inf
 
     return v
 

setup.py

@@ -62,7 +62,7 @@ def find_package_data(datadir=None):
         'numpy >= 1.9.0',
         'scipy >= 0.14.0',
         'healpy >= 1.6.0',
-        'pyyaml >= 3.10',
+        'pyyaml >= 3.12',
         'pymodeler >= 0.1.0',
     ],
     packages=find_packages(),

tests/test_jcalc.py

@@ -71,9 +71,9 @@ def test_profiles_2d():
         fast = LoSIntegralFast(dp,dhalo)
         itrp = LoSIntegralInterp(dp,dhalo)
 
-        slow_jval = slow(_psi,_dhalo)
-        fast_jval = fast(_psi,_dhalo)
-        itrp_jval = itrp(_psi,_dhalo)
+        slow_jval = slow(_psi,_dhalo).clip(1e-99, np.inf)
+        fast_jval = fast(_psi,_dhalo).clip(1e-99, np.inf)
+        itrp_jval = itrp(_psi,_dhalo).clip(1e-99, np.inf)
 
         fast_vs_slow = (fast_jval-slow_jval)/slow_jval
         itrp_vs_slow = (itrp_jval-slow_jval)/slow_jval