[yt-svn] commit/yt: xarthisius: Merged in jzuhone/yt (pull request #1802)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Oct 20 08:39:34 PDT 2015
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/35741b799d57/
Changeset: 35741b799d57
Branch: yt
User: xarthisius
Date: 2015-10-20 15:39:21+00:00
Summary: Merged in jzuhone/yt (pull request #1802)
[bugfix] Fixing bugs and big speedups for the spectral models in photon_simulator
Affected #: 7 files
diff -r 96e7ac3e7175758b3c0b2a780d3ef0920adaaaa9 -r 35741b799d57ac96134c4f869cfa92e8c4f9e2be .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -10,6 +10,7 @@
yt/analysis_modules/halo_finding/rockstar/rockstar_groupies.c
yt/analysis_modules/halo_finding/rockstar/rockstar_interface.c
yt/analysis_modules/ppv_cube/ppv_utils.c
+yt/analysis_modules/photon_simulator/utils.c
yt/frontends/ramses/_ramses_reader.cpp
yt/frontends/sph/smoothing_kernel.c
yt/geometry/fake_octree.c
diff -r 96e7ac3e7175758b3c0b2a780d3ef0920adaaaa9 -r 35741b799d57ac96134c4f869cfa92e8c4f9e2be doc/source/analyzing/analysis_modules/photon_simulator.rst
--- a/doc/source/analyzing/analysis_modules/photon_simulator.rst
+++ b/doc/source/analyzing/analysis_modules/photon_simulator.rst
@@ -132,7 +132,7 @@
.. code:: python
- apec_model = TableApecModel("atomdb_v2.0.2",
+ apec_model = TableApecModel("$SPECTRAL_DATA/spectral",
0.01, 20.0, 20000,
thermal_broad=False,
apec_vers="2.0.2")
diff -r 96e7ac3e7175758b3c0b2a780d3ef0920adaaaa9 -r 35741b799d57ac96134c4f869cfa92e8c4f9e2be yt/analysis_modules/photon_simulator/photon_simulator.py
--- a/yt/analysis_modules/photon_simulator/photon_simulator.py
+++ b/yt/analysis_modules/photon_simulator/photon_simulator.py
@@ -577,7 +577,7 @@
"responses are normalized properly. If not, you may "+
"get inconsistent results.")
f.close()
- Aratio = eff_area.max()/self.parameters["FiducialArea"]
+ Aratio = eff_area.max()/self.parameters["FiducialArea"].v
else:
mylog.info("Using constant effective area.")
Aratio = parse_value(area_new, "cm**2")/self.parameters["FiducialArea"]
diff -r 96e7ac3e7175758b3c0b2a780d3ef0920adaaaa9 -r 35741b799d57ac96134c4f869cfa92e8c4f9e2be yt/analysis_modules/photon_simulator/setup.py
--- a/yt/analysis_modules/photon_simulator/setup.py
+++ b/yt/analysis_modules/photon_simulator/setup.py
@@ -4,6 +4,8 @@
def configuration(parent_package='', top_path=None):
from numpy.distutils.misc_util import Configuration
config = Configuration('photon_simulator', parent_package, top_path)
+ config.add_extension("utils",
+ ["yt/analysis_modules/photon_simulator/utils.pyx"])
config.add_subpackage("tests")
config.make_config_py() # installs __config__.py
#config.make_svn_version_py()
diff -r 96e7ac3e7175758b3c0b2a780d3ef0920adaaaa9 -r 35741b799d57ac96134c4f869cfa92e8c4f9e2be yt/analysis_modules/photon_simulator/spectral_models.py
--- a/yt/analysis_modules/photon_simulator/spectral_models.py
+++ b/yt/analysis_modules/photon_simulator/spectral_models.py
@@ -17,8 +17,10 @@
from yt.funcs import mylog
from yt.units.yt_array import YTArray, YTQuantity
-from yt.utilities.on_demand_imports import _astropy, _scipy
-from yt.utilities.physical_constants import hcgs, clight, erg_per_keV, amu_cgs
+from yt.utilities.on_demand_imports import _astropy
+from yt.utilities.physical_constants import hcgs, clight
+from yt.utilities.physical_ratios import erg_per_keV, amu_grams
+from yt.analysis_modules.photon_simulator.utils import broaden_lines
hc = (hcgs*clight).in_units("keV*angstrom").v
cl = clight.v
@@ -30,7 +32,7 @@
self.emin = YTQuantity(emin, "keV")
self.emax = YTQuantity(emax, "keV")
self.nchan = nchan
- self.ebins = np.linspace(self.emin, self.emax, nchan+1)
+ self.ebins = YTArray(np.linspace(self.emin, self.emax, nchan+1), "keV")
self.de = np.diff(self.ebins)
self.emid = 0.5*(self.ebins[1:]+self.ebins[:-1])
@@ -63,7 +65,7 @@
--------
>>> mekal_model = XSpecThermalModel("mekal", 0.05, 50.0, 1000)
"""
- def __init__(self, model_name, emin, emax, nchan,
+ def __init__(self, model_name, emin, emax, nchan,
thermal_broad=False, settings=None):
self.model_name = model_name
self.thermal_broad = thermal_broad
@@ -133,7 +135,7 @@
--------
>>> abs_model = XSpecAbsorbModel("wabs", 0.1)
"""
- def __init__(self, model_name, nH, emin=0.01, emax=50.0,
+ def __init__(self, model_name, nH, emin=0.01, emax=50.0,
nchan=100000, settings=None):
self.model_name = model_name
self.nH = nH
@@ -191,8 +193,8 @@
Examples
--------
- >>> apec_model = TableApecModel("/Users/jzuhone/Data/atomdb_v2.0.2/", 0.05, 50.0,
- ... 1000, thermal_broad=True)
+ >>> apec_model = TableApecModel("$SPECTRAL_DATA/spectral/", 0.05, 50.0,
+ ... 1000, apec_vers="3.0", thermal_broad=True)
"""
def __init__(self, apec_root, emin, emax, nchan,
apec_vers="2.0.2", thermal_broad=False):
@@ -202,7 +204,7 @@
self.apec_prefix+"_coco.fits")
self.linefile = os.path.join(self.apec_root,
self.apec_prefix+"_line.fits")
- super(TableApecModel, self).__init__(emin, emax, nchan)
+ super(TableApecModel, self).__init__(emin, emax, nchan)
self.wvbins = hc/self.ebins[::-1].d
# H, He, and trace elements
self.cosmic_elem = [1,2,3,4,5,9,11,15,17,19,21,22,23,24,25,27,29,30]
@@ -237,47 +239,47 @@
self.maxlam = self.wvbins.max()
self.scale_factor = 1.0/(1.+zobs)
- def _make_spectrum(self, element, tindex):
+ def _make_spectrum(self, kT, element, tindex):
tmpspec = np.zeros(self.nchan)
- i = np.where((self.line_handle[tindex].data.field('element') == element) &
- (self.line_handle[tindex].data.field('lambda') > self.minlam) &
- (self.line_handle[tindex].data.field('lambda') < self.maxlam))[0]
+ line_data = self.line_handle[tindex].data
+ coco_data = self.coco_handle[tindex].data
- vec = np.zeros(self.nchan)
- E0 = hc/self.line_handle[tindex].data.field('lambda')[i]*self.scale_factor
- amp = self.line_handle[tindex].data.field('epsilon')[i]
+ i = np.where((line_data.field('element') == element) &
+ (line_data.field('lambda') > self.minlam) &
+ (line_data.field('lambda') < self.maxlam))[0]
+
+ E0 = hc/line_data.field('lambda')[i].astype("float64")*self.scale_factor
+ amp = line_data.field('epsilon')[i].astype("float64")
ebins = self.ebins.d
+ de = self.de.d
+ emid = self.emid.d
if self.thermal_broad:
- vec = np.zeros(self.nchan)
- sigma = E0*np.sqrt(self.Tvals[tindex]*erg_per_keV/(self.A[element]*amu_cgs))/cl
- for E, sig, a in zip(E0, sigma, amp):
- cdf = _scipy.stats.norm(E,sig).cdf(ebins)
- vec += np.diff(cdf)*a
+ sigma = E0*np.sqrt(2.*kT*erg_per_keV/(self.A[element]*amu_grams))/cl
+ vec = broaden_lines(E0, sigma, amp, emid)*de
else:
- ie = np.searchsorted(ebins, E0, side='right')-1
- for i, a in zip(ie, amp): vec[i] += a
+ vec = np.histogram(E0, ebins, weights=amp)[0]
tmpspec += vec
- ind = np.where((self.coco_handle[tindex].data.field('Z') == element) &
- (self.coco_handle[tindex].data.field('rmJ') == 0))[0]
+ ind = np.where((coco_data.field('Z') == element) &
+ (coco_data.field('rmJ') == 0))[0]
if len(ind) == 0:
return tmpspec
else:
ind = ind[0]
- n_cont = self.coco_handle[tindex].data.field('N_Cont')[ind]
- e_cont = self.coco_handle[tindex].data.field('E_Cont')[ind][:n_cont]*self.scale_factor
- continuum = self.coco_handle[tindex].data.field('Continuum')[ind][:n_cont]
+ n_cont = coco_data.field('N_Cont')[ind]
+ e_cont = coco_data.field('E_Cont')[ind][:n_cont]
+ continuum = coco_data.field('Continuum')[ind][:n_cont]
- tmpspec += np.interp(self.emid.d, e_cont, continuum)*self.de.d
+ tmpspec += np.interp(emid, e_cont*self.scale_factor, continuum)*de/self.scale_factor
- n_pseudo = self.coco_handle[tindex].data.field('N_Pseudo')[ind]
- e_pseudo = self.coco_handle[tindex].data.field('E_Pseudo')[ind][:n_pseudo]*self.scale_factor
- pseudo = self.coco_handle[tindex].data.field('Pseudo')[ind][:n_pseudo]
+ n_pseudo = coco_data.field('N_Pseudo')[ind]
+ e_pseudo = coco_data.field('E_Pseudo')[ind][:n_pseudo]
+ pseudo = coco_data.field('Pseudo')[ind][:n_pseudo]
- tmpspec += np.interp(self.emid.d, e_pseudo, pseudo)*self.de.d
+ tmpspec += np.interp(emid, e_pseudo*self.scale_factor, pseudo)*de/self.scale_factor
return tmpspec
@@ -295,12 +297,12 @@
dT = (kT-self.Tvals[tindex])/self.dTvals[tindex]
# First do H,He, and trace elements
for elem in self.cosmic_elem:
- cspec_l += self._make_spectrum(elem, tindex+2)
- cspec_r += self._make_spectrum(elem, tindex+3)
+ cspec_l += self._make_spectrum(kT, elem, tindex+2)
+ cspec_r += self._make_spectrum(kT, elem, tindex+3)
# Next do the metals
for elem in self.metal_elem:
- mspec_l += self._make_spectrum(elem, tindex+2)
- mspec_r += self._make_spectrum(elem, tindex+3)
+ mspec_l += self._make_spectrum(kT, elem, tindex+2)
+ mspec_r += self._make_spectrum(kT, elem, tindex+3)
cosmic_spec = YTArray(cspec_l*(1.-dT)+cspec_r*dT, "cm**3/s")
metal_spec = YTArray(mspec_l*(1.-dT)+mspec_r*dT, "cm**3/s")
return cosmic_spec, metal_spec
diff -r 96e7ac3e7175758b3c0b2a780d3ef0920adaaaa9 -r 35741b799d57ac96134c4f869cfa92e8c4f9e2be yt/analysis_modules/photon_simulator/tests/test_spectra.py
--- /dev/null
+++ b/yt/analysis_modules/photon_simulator/tests/test_spectra.py
@@ -0,0 +1,44 @@
+from yt.analysis_modules.photon_simulator.api import \
+ TableApecModel, XSpecThermalModel
+import numpy as np
+from yt.testing import requires_module, fake_random_ds
+from yt.utilities.answer_testing.framework import \
+ GenericArrayTest, data_dir_load
+from yt.config import ytcfg
+
+def setup():
+ ytcfg["yt", "__withintesting"] = "True"
+
+test_data_dir = ytcfg.get("yt", "test_data_dir")
+
+ds = fake_random_ds(64)
+
+ at requires_module("xspec")
+ at requires_module("astropy")
+def test_apec():
+
+ settings = {"APECROOT":test_data_dir+"/xray_data/apec_v2.0.2"}
+ xmod = XSpecThermalModel("apec", 0.1, 10.0, 10000, thermal_broad=True,
+ settings=settings)
+ xmod.prepare_spectrum(0.2)
+
+ xcspec, xmspec = xmod.get_spectrum(6.0)
+ spec1 = xcspec+0.3*xmspec
+
+ amod = TableApecModel(test_data_dir+"/xray_data", 0.1, 10.0,
+ 10000, thermal_broad=True)
+ amod.prepare_spectrum(0.2)
+
+ acspec, amspec = amod.get_spectrum(6.0)
+ spec2 = acspec+0.3*amspec
+
+ def spec1_test():
+ return spec1.v
+ def spec2_test():
+ return spec2.v
+
+ for test in [GenericArrayTest(ds, spec1_test),
+ GenericArrayTest(ds, spec2_test)]:
+ test_apec.__name__ = test.description
+ yield test
+
diff -r 96e7ac3e7175758b3c0b2a780d3ef0920adaaaa9 -r 35741b799d57ac96134c4f869cfa92e8c4f9e2be yt/analysis_modules/photon_simulator/utils.pyx
--- /dev/null
+++ b/yt/analysis_modules/photon_simulator/utils.pyx
@@ -0,0 +1,31 @@
+import numpy as np
+cimport numpy as np
+cimport cython
+from libc.math cimport exp
+
+cdef double gfac = 1.0/np.sqrt(np.pi)
+
+ at cython.cdivision(True)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+def broaden_lines(np.ndarray[np.float64_t, ndim=1] E0,
+ np.ndarray[np.float64_t, ndim=1] sigma,
+ np.ndarray[np.float64_t, ndim=1] amp,
+ np.ndarray[np.float64_t, ndim=1] E):
+
+ cdef int i, j, n
+ cdef double x, isigma, iamp
+ cdef np.ndarray[np.float64_t, ndim=1] lines
+
+ n = E0.shape[0]
+ m = E.shape[0]
+ lines = np.zeros(m)
+
+ for i in range(n):
+ isigma = 1.0/sigma[i]
+ iamp = gfac*amp[i]*isigma
+ for j in range(m):
+ x = (E[j]-E0[i])*isigma
+ lines[j] += iamp*exp(-x*x)
+
+ return lines
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list