[yt-svn] commit/yt: 14 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Oct 20 08:39:33 PDT 2015
14 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/b12a659d76f5/
Changeset: b12a659d76f5
Branch: yt
User: jzuhone
Date: 2015-10-14 14:43:36+00:00
Summary: Changing this in the example to something a bit less confusing
Affected #: 1 file
diff -r ca754f897f2f44a4c2a9415b3aed977a0d1a4df0 -r b12a659d76f541bf5863e0fabd50e4055d452ee2 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("/Users/jzuhone/spectral",
0.01, 20.0, 20000,
thermal_broad=False,
apec_vers="2.0.2")
https://bitbucket.org/yt_analysis/yt/commits/d5fb9f08fcea/
Changeset: d5fb9f08fcea
Branch: yt
User: jzuhone
Date: 2015-10-14 14:45:00+00:00
Summary: Simplifications to TableApecModel and massive speedups when setting thermal_broad = True
Affected #: 3 files
diff -r b12a659d76f541bf5863e0fabd50e4055d452ee2 -r d5fb9f08fceaf6f4c3453f087145c7a911e71545 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,9 @@
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"],
+ libraries=["m"])
config.add_subpackage("tests")
config.make_config_py() # installs __config__.py
#config.make_svn_version_py()
diff -r b12a659d76f541bf5863e0fabd50e4055d452ee2 -r d5fb9f08fceaf6f4c3453f087145c7a911e71545 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 compute_lines
hc = (hcgs*clight).in_units("keV*angstrom").v
cl = clight.v
@@ -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("/Users/jzuhone/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]
@@ -241,43 +243,43 @@
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(self.Tvals[tindex]*erg_per_keV/(self.A[element]*amu_grams))/cl
+ vec = compute_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]*self.scale_factor
+ 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, continuum)*de
- 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]*self.scale_factor
+ 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, pseudo)*de
return tmpspec
@@ -296,7 +298,7 @@
# 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_r += self._make_spectrum(elem, tindex+3)
# Next do the metals
for elem in self.metal_elem:
mspec_l += self._make_spectrum(elem, tindex+2)
diff -r b12a659d76f541bf5863e0fabd50e4055d452ee2 -r d5fb9f08fceaf6f4c3453f087145c7a911e71545 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
+
+cdef extern from "math.h":
+ double exp(double x) nogil
+
+cdef double gfac = 1.0/np.sqrt(2*np.pi)
+
+ at cython.cdivision(True)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+def compute_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
+ 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):
+ for j in range(m):
+ x = (E[j]-E0[i])/sigma[i]
+ lines[j] += amp[i]*gfac*exp(-0.5*x*x)/sigma[i]
+
+ return lines
https://bitbucket.org/yt_analysis/yt/commits/d1f05ff2cbde/
Changeset: d1f05ff2cbde
Branch: yt
User: jzuhone
Date: 2015-10-14 17:27:58+00:00
Summary: Frigging bugs
Affected #: 1 file
diff -r d5fb9f08fceaf6f4c3453f087145c7a911e71545 -r d1f05ff2cbde1ac9e6b3568058e15c41b5cc4c52 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
@@ -65,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
@@ -270,13 +270,13 @@
ind = ind[0]
n_cont = coco_data.field('N_Cont')[ind]
- e_cont = coco_data.field('E_Cont')[ind][:n_cont]*self.scale_factor
+ e_cont = coco_data.field('E_Cont')[ind][:n_cont]
continuum = coco_data.field('Continuum')[ind][:n_cont]
tmpspec += np.interp(emid, e_cont, continuum)*de
n_pseudo = coco_data.field('N_Pseudo')[ind]
- e_pseudo = coco_data.field('E_Pseudo')[ind][:n_pseudo]*self.scale_factor
+ e_pseudo = coco_data.field('E_Pseudo')[ind][:n_pseudo]
pseudo = coco_data.field('Pseudo')[ind][:n_pseudo]
tmpspec += np.interp(emid, e_pseudo, pseudo)*de
https://bitbucket.org/yt_analysis/yt/commits/aaa8d5cfeeb7/
Changeset: aaa8d5cfeeb7
Branch: yt
User: jzuhone
Date: 2015-10-14 20:00:55+00:00
Summary: Bugfixes
Affected #: 1 file
diff -r d1f05ff2cbde1ac9e6b3568058e15c41b5cc4c52 -r aaa8d5cfeeb7050df0e91d98cbc72aef41c69a08 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
@@ -239,7 +239,7 @@
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)
@@ -250,16 +250,16 @@
(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
+ E0 = hc/line_data.field('lambda')[i].astype("float64")
amp = line_data.field('epsilon')[i].astype("float64")
ebins = self.ebins.d
de = self.de.d
emid = self.emid.d
if self.thermal_broad:
- sigma = E0*np.sqrt(self.Tvals[tindex]*erg_per_keV/(self.A[element]*amu_grams))/cl
- vec = compute_lines(E0, sigma, amp, emid)*de
+ sigma = E0*np.sqrt(2.*kT*erg_per_keV/(self.A[element]*amu_grams))/cl
+ vec = compute_lines(E0*self.scale_factor, sigma, amp, emid)*de
else:
- vec = np.histogram(E0, ebins, weights=amp)[0]
+ vec = np.histogram(E0*self.scale_factor, ebins, weights=amp)[0]*self.scale_factor
tmpspec += vec
ind = np.where((coco_data.field('Z') == element) &
@@ -273,13 +273,15 @@
e_cont = coco_data.field('E_Cont')[ind][:n_cont]
continuum = coco_data.field('Continuum')[ind][:n_cont]
- tmpspec += np.interp(emid, e_cont, continuum)*de
+ tmpspec += np.interp(emid, e_cont*self.scale_factor, continuum)*de
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(emid, e_pseudo, pseudo)*de
+ tmpspec += np.interp(emid, e_pseudo*self.scale_factor, pseudo)*de
+
+ tmpspec /= self.scale_factor
return tmpspec
@@ -297,8 +299,8 @@
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)
https://bitbucket.org/yt_analysis/yt/commits/f201f33dca63/
Changeset: f201f33dca63
Branch: yt
User: jzuhone
Date: 2015-10-14 20:33:11+00:00
Summary: More bugfixes
Affected #: 2 files
diff -r aaa8d5cfeeb7050df0e91d98cbc72aef41c69a08 -r f201f33dca63323c70f703b075d773d9fa337a07 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
@@ -32,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])
@@ -250,17 +250,17 @@
(line_data.field('lambda') > self.minlam) &
(line_data.field('lambda') < self.maxlam))[0]
- E0 = hc/line_data.field('lambda')[i].astype("float64")
+ 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:
sigma = E0*np.sqrt(2.*kT*erg_per_keV/(self.A[element]*amu_grams))/cl
- vec = compute_lines(E0*self.scale_factor, sigma, amp, emid)*de
+ vec = compute_lines(E0, sigma, amp, emid)*de
else:
- vec = np.histogram(E0*self.scale_factor, ebins, weights=amp)[0]*self.scale_factor
- tmpspec += vec
+ vec = np.histogram(E0, ebins, weights=amp)[0]
+ tmpspec += vec*self.scale_factor
ind = np.where((coco_data.field('Z') == element) &
(coco_data.field('rmJ') == 0))[0]
@@ -303,8 +303,8 @@
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 aaa8d5cfeeb7050df0e91d98cbc72aef41c69a08 -r f201f33dca63323c70f703b075d773d9fa337a07 yt/analysis_modules/photon_simulator/utils.pyx
--- a/yt/analysis_modules/photon_simulator/utils.pyx
+++ b/yt/analysis_modules/photon_simulator/utils.pyx
@@ -5,7 +5,7 @@
cdef extern from "math.h":
double exp(double x) nogil
-cdef double gfac = 1.0/np.sqrt(2*np.pi)
+cdef double gfac = 1.0/np.sqrt(np.pi)
@cython.cdivision(True)
@cython.boundscheck(False)
@@ -26,6 +26,6 @@
for i in range(n):
for j in range(m):
x = (E[j]-E0[i])/sigma[i]
- lines[j] += amp[i]*gfac*exp(-0.5*x*x)/sigma[i]
+ lines[j] += amp[i]*gfac*exp(-x*x)/sigma[i]
return lines
https://bitbucket.org/yt_analysis/yt/commits/96c9fd7f28fd/
Changeset: 96c9fd7f28fd
Branch: yt
User: jzuhone
Date: 2015-10-14 21:14:53+00:00
Summary: Minor nits
Affected #: 2 files
diff -r f201f33dca63323c70f703b075d773d9fa337a07 -r 96c9fd7f28fde33ae7b7da052c87359d3ff363ae 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
@@ -20,7 +20,7 @@
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 compute_lines
+from yt.analysis_modules.photon_simulator.utils import broaden_lines
hc = (hcgs*clight).in_units("keV*angstrom").v
cl = clight.v
@@ -257,10 +257,10 @@
emid = self.emid.d
if self.thermal_broad:
sigma = E0*np.sqrt(2.*kT*erg_per_keV/(self.A[element]*amu_grams))/cl
- vec = compute_lines(E0, sigma, amp, emid)*de
+ vec = broaden_lines(E0, sigma, amp, emid)*de
else:
vec = np.histogram(E0, ebins, weights=amp)[0]
- tmpspec += vec*self.scale_factor
+ tmpspec += vec
ind = np.where((coco_data.field('Z') == element) &
(coco_data.field('rmJ') == 0))[0]
@@ -273,15 +273,13 @@
e_cont = coco_data.field('E_Cont')[ind][:n_cont]
continuum = coco_data.field('Continuum')[ind][:n_cont]
- tmpspec += np.interp(emid, e_cont*self.scale_factor, continuum)*de
+ tmpspec += np.interp(emid, e_cont*self.scale_factor, continuum)*de/self.scale_factor
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(emid, e_pseudo*self.scale_factor, pseudo)*de
-
- tmpspec /= self.scale_factor
+ tmpspec += np.interp(emid, e_pseudo*self.scale_factor, pseudo)*de/self.scale_factor
return tmpspec
diff -r f201f33dca63323c70f703b075d773d9fa337a07 -r 96c9fd7f28fde33ae7b7da052c87359d3ff363ae yt/analysis_modules/photon_simulator/utils.pyx
--- a/yt/analysis_modules/photon_simulator/utils.pyx
+++ b/yt/analysis_modules/photon_simulator/utils.pyx
@@ -10,7 +10,7 @@
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
-def compute_lines(np.ndarray[np.float64_t, ndim=1] E0,
+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):
https://bitbucket.org/yt_analysis/yt/commits/aa9897c51d64/
Changeset: aa9897c51d64
Branch: yt
User: jzuhone
Date: 2015-10-14 23:38:36+00:00
Summary: Adding APEC test
Affected #: 1 file
diff -r 96c9fd7f28fde33ae7b7da052c87359d3ff363ae -r aa9897c51d64d17e329ef3fcc251896954bb4e89 yt/analysis_modules/photon_simulator/tests/test_apec.py
--- /dev/null
+++ b/yt/analysis_modules/photon_simulator/tests/test_apec.py
@@ -0,0 +1,43 @@
+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
+from yt.config import ytcfg
+
+def setup():
+ from yt.config import ytcfg
+ ytcfg["yt", "__withintesting"] = "True"
+
+test_data_dir = ytcfg.get("yt", "test_data_dir")
+
+ at requires_module("xspec")
+ at requires_module("astropy")
+def test_apec():
+ ds = fake_random_ds(64)
+ 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
+ def spec2_test():
+ return spec2
+
+ for test in [GenericArrayTest(ds, spec1_test),
+ GenericArrayTest(ds, spec2_test)]:
+ test_apec.__name__ = test.description
+ yield test
+
https://bitbucket.org/yt_analysis/yt/commits/76543f1c0925/
Changeset: 76543f1c0925
Branch: yt
User: jzuhone
Date: 2015-10-15 14:18:53+00:00
Summary: An answer test for spectra
Affected #: 2 files
diff -r aa9897c51d64d17e329ef3fcc251896954bb4e89 -r 76543f1c092507574edc972b8ec739819a28a2e3 yt/analysis_modules/photon_simulator/tests/test_apec.py
--- a/yt/analysis_modules/photon_simulator/tests/test_apec.py
+++ /dev/null
@@ -1,43 +0,0 @@
-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
-from yt.config import ytcfg
-
-def setup():
- from yt.config import ytcfg
- ytcfg["yt", "__withintesting"] = "True"
-
-test_data_dir = ytcfg.get("yt", "test_data_dir")
-
- at requires_module("xspec")
- at requires_module("astropy")
-def test_apec():
- ds = fake_random_ds(64)
- 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
- def spec2_test():
- return spec2
-
- for test in [GenericArrayTest(ds, spec1_test),
- GenericArrayTest(ds, spec2_test)]:
- test_apec.__name__ = test.description
- yield test
-
diff -r aa9897c51d64d17e329ef3fcc251896954bb4e89 -r 76543f1c092507574edc972b8ec739819a28a2e3 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
+
https://bitbucket.org/yt_analysis/yt/commits/2d286c0b273b/
Changeset: 2d286c0b273b
Branch: yt
User: jzuhone
Date: 2015-10-15 15:20:36+00:00
Summary: Slightly more efficient (thanks Kacper!)
Affected #: 1 file
diff -r 76543f1c092507574edc972b8ec739819a28a2e3 -r 2d286c0b273b1022a861f83f629e384387d4a8d5 yt/analysis_modules/photon_simulator/utils.pyx
--- a/yt/analysis_modules/photon_simulator/utils.pyx
+++ b/yt/analysis_modules/photon_simulator/utils.pyx
@@ -16,7 +16,7 @@
np.ndarray[np.float64_t, ndim=1] E):
cdef int i, j, n
- cdef double x
+ cdef double x, isigma, iamp
cdef np.ndarray[np.float64_t, ndim=1] lines
n = E0.shape[0]
@@ -24,8 +24,10 @@
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])/sigma[i]
- lines[j] += amp[i]*gfac*exp(-x*x)/sigma[i]
+ x = (E[j]-E0[i])*isigma
+ lines[j] += iamp*exp(-x*x)
return lines
https://bitbucket.org/yt_analysis/yt/commits/e6f8446d910c/
Changeset: e6f8446d910c
Branch: yt
User: jzuhone
Date: 2015-10-15 17:11:59+00:00
Summary: Make sure this ends up dimensionless
Affected #: 1 file
diff -r 2d286c0b273b1022a861f83f629e384387d4a8d5 -r e6f8446d910cd1e6cc39b1250f48850b54329ea0 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"]
https://bitbucket.org/yt_analysis/yt/commits/6a5e5975b742/
Changeset: 6a5e5975b742
Branch: yt
User: jzuhone
Date: 2015-10-16 13:31:21+00:00
Summary: Making this less specific
Affected #: 2 files
diff -r e6f8446d910cd1e6cc39b1250f48850b54329ea0 -r 6a5e5975b7424e0e511959e92cbf2f714e2b45ab 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("/Users/jzuhone/spectral",
+ apec_model = TableApecModel("$SPECTRAL_DATA/spectral",
0.01, 20.0, 20000,
thermal_broad=False,
apec_vers="2.0.2")
diff -r e6f8446d910cd1e6cc39b1250f48850b54329ea0 -r 6a5e5975b7424e0e511959e92cbf2f714e2b45ab 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
@@ -193,7 +193,7 @@
Examples
--------
- >>> apec_model = TableApecModel("/Users/jzuhone/Data/spectral/", 0.05, 50.0,
+ >>> 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,
https://bitbucket.org/yt_analysis/yt/commits/7ba025fe57b0/
Changeset: 7ba025fe57b0
Branch: yt
User: jzuhone
Date: 2015-10-16 13:32:03+00:00
Summary: Don't need to see this
Affected #: 1 file
diff -r 6a5e5975b7424e0e511959e92cbf2f714e2b45ab -r 7ba025fe57b0da71a126b91c96ffc28a8176a1f9 .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
https://bitbucket.org/yt_analysis/yt/commits/5438e77cc0c9/
Changeset: 5438e77cc0c9
Branch: yt
User: jzuhone
Date: 2015-10-16 14:02:25+00:00
Summary: Import the exp function this way
Affected #: 2 files
diff -r 7ba025fe57b0da71a126b91c96ffc28a8176a1f9 -r 5438e77cc0c94ab2038bbc78462911de445f6b2e yt/analysis_modules/photon_simulator/setup.py
--- a/yt/analysis_modules/photon_simulator/setup.py
+++ b/yt/analysis_modules/photon_simulator/setup.py
@@ -5,8 +5,7 @@
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"],
- libraries=["m"])
+ ["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 7ba025fe57b0da71a126b91c96ffc28a8176a1f9 -r 5438e77cc0c94ab2038bbc78462911de445f6b2e yt/analysis_modules/photon_simulator/utils.pyx
--- a/yt/analysis_modules/photon_simulator/utils.pyx
+++ b/yt/analysis_modules/photon_simulator/utils.pyx
@@ -1,9 +1,7 @@
import numpy as np
cimport numpy as np
cimport cython
-
-cdef extern from "math.h":
- double exp(double x) nogil
+from libc.math cimport exp
cdef double gfac = 1.0/np.sqrt(np.pi)
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