[yt-svn] commit/yt: 14 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Nov 16 13:27:24 PST 2015
14 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/55ba775f16f6/
Changeset: 55ba775f16f6
Branch: yt
User: chummels
Date: 2015-11-11 22:18:32+00:00
Summary: Correcting max_tau to be min_tau.
Affected #: 1 file
diff -r 25a0dbd51257bb3fd355338fd4594985c5963f77 -r 55ba775f16f670f5c5f7b667227ee16885ab1bea yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -203,7 +203,7 @@
Add continuum features to the spectrum.
"""
# Only add continuum features down to tau of 1.e-4.
- tau_min = 1.e-4
+ min_tau = 1.e-4
for continuum in self.continuum_list:
column_density = field_data[continuum['field_name']] * field_data['dl']
@@ -215,12 +215,12 @@
this_wavelength = delta_lambda + continuum['wavelength']
right_index = np.digitize(this_wavelength, self.lambda_bins).clip(0, self.n_lambda)
left_index = np.digitize((this_wavelength *
- np.power((tau_min * continuum['normalization'] /
+ np.power((min_tau * continuum['normalization'] /
column_density), (1. / continuum['index']))),
self.lambda_bins).clip(0, self.n_lambda)
valid_continuua = np.where(((column_density /
- continuum['normalization']) > tau_min) &
+ continuum['normalization']) > min_tau) &
(right_index - left_index > 1))[0]
pbar = get_pbar("Adding continuum feature - %s [%f A]: " % \
(continuum['label'], continuum['wavelength']),
@@ -240,8 +240,9 @@
"""
# Only make voigt profile for slice of spectrum that is 10 times the line width.
spectrum_bin_ratio = 5
- # Widen wavelength window until optical depth reaches a max value at the ends.
- max_tau = 0.001
+ # Widen wavelength window until optical depth reaches a min value at
+ # the ends to assure that the wings of a line have been fully resolved.
+ min_tau = 0.001
for line in parallel_objects(self.line_list, njobs=njobs):
column_density = field_data[line['field_name']] * field_data['dl']
@@ -297,7 +298,7 @@
lambda_bins=lbins[left_index[lixel]:right_index[lixel]])
# Widen wavelength window until optical depth reaches a max value at the ends.
- if (line_tau[0] < max_tau and line_tau[-1] < max_tau) or \
+ if (line_tau[0] < min_tau and line_tau[-1] < min_tau) or \
(left_index[lixel] <= 0 and right_index[lixel] >= self.n_lambda):
break
my_bin_ratio *= 2
https://bitbucket.org/yt_analysis/yt/commits/4278eff7a0ed/
Changeset: 4278eff7a0ed
Branch: yt
User: chummels
Date: 2015-11-12 23:38:56+00:00
Summary: Refactoring absorption spectrum to deal with unresolved lines
Affected #: 1 file
diff -r 55ba775f16f670f5c5f7b667227ee16885ab1bea -r 4278eff7a0edf7e474092313445e9ebc9efbad91 yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -238,95 +238,132 @@
"""
Add the absorption lines to the spectrum.
"""
- # Only make voigt profile for slice of spectrum that is 10 times the line width.
- spectrum_bin_ratio = 5
- # Widen wavelength window until optical depth reaches a min value at
- # the ends to assure that the wings of a line have been fully resolved.
+ # Widen wavelength window until optical depth falls below this tau
+ # value at the ends to assure that the wings of a line have been
+ # fully resolved.
min_tau = 0.001
+ # step through each ionic transition (e.g. HI, HII, MgII) specified
+ # and deposit the lines into the spectrum
for line in parallel_objects(self.line_list, njobs=njobs):
column_density = field_data[line['field_name']] * field_data['dl']
+
# redshift_eff field combines cosmological and velocity redshifts
+ # so delta_lambda gives the offset in angstroms from the rest frame
+ # wavelength to the observed wavelength of the transition
if use_peculiar_velocity:
delta_lambda = line['wavelength'] * field_data['redshift_eff']
else:
delta_lambda = line['wavelength'] * field_data['redshift']
+ # lambda_obs is central wavelength of line after redshift
+ lambda_obs = line['wavelength'] + delta_lambda
+ # thermal broadening b parameter
thermal_b = np.sqrt((2 * boltzmann_constant_cgs *
field_data['temperature']) /
line['atomic_mass'])
- center_bins = np.digitize((delta_lambda + line['wavelength']),
- self.lambda_bins)
- # ratio of line width to bin width
- width_ratio = ((line['wavelength'] + delta_lambda) * \
- thermal_b / speed_of_light_cgs / self.bin_width).in_units("").d
+ # the actual thermal width of the lines
+ thermal_width = (lambda_obs * thermal_b /
+ speed_of_light_cgs).convert_to_units("angstrom")
- if (width_ratio < 1.0).any():
- mylog.warn(("%d out of %d line components are unresolved, " +
- "consider increasing spectral resolution.") %
- ((width_ratio < 1.0).sum(), width_ratio.size))
+ if (thermal_width < self.bin_width).any():
+ mylog.warn(("%d out of %d line components are unresolved.") %
+ ((thermal_width < self.bin_width).sum(),
+ thermal_width.size))
- # do voigt profiles for a subset of the full spectrum
- left_index = (center_bins -
- spectrum_bin_ratio * width_ratio).astype(int).clip(0, self.n_lambda)
- right_index = (center_bins +
- spectrum_bin_ratio * width_ratio).astype(int).clip(0, self.n_lambda)
+ # Sanitize units for faster runtime of the tau_profile machinery.
+ lambda_0 = line['wavelength'].d # line's rest frame; angstroms
+ lambda_1 = lambda_obs.d # line's observed frame; angstroms
+ cdens = column_density.in_units("cm**-2").d # cm**-2
+ thermb = thermal_b.in_cgs().d # thermal b coefficient; cm / s
+ dlambda = delta_lambda.d # lambda offset; angstroms
+ vlos = field_data['velocity_los'].in_units("km/s").d # km/s
- # loop over all lines wider than the bin width
- valid_lines = np.where((width_ratio >= 1.0) &
- (right_index - left_index > 1))[0]
+ # When we actually deposit the voigt profile, sometimes we will
+ # have underresolved lines (ie lines with smaller widths than
+ # the spectral bin size). Here, we create virtual bins small
+ # enough in width to well resolve each line, deposit the voigt
+ # profile into them, then numerically integrate their tau values
+ # and sum them to redeposit them into the actual spectral bins.
+
+ # virtual bins will be the smaller of 2 options:
+ # --1/10th of the thermal width
+ # --the actual spectral bin width
+ # to start, the window over which the voigt profile will be
+ # deposited will be 50 v_bins, but this may grow if necessary
+ vbin_width = np.amin(zip(thermal_width / 10.,
+ self.bin_width *
+ np.ones(len(thermal_width))),
+ axis=1)
+ n_vbins = 50
+ vbin_window_width = n_vbins*vbin_width
+
+ valid_lines = np.arange(len(thermal_width))
pbar = get_pbar("Adding line - %s [%f A]: " % (line['label'], line['wavelength']),
- valid_lines.size)
+ thermal_width.size)
- # Sanitize units here
- column_density.convert_to_units("cm ** -2")
- lbins = self.lambda_bins.d # Angstroms
- lambda_0 = line['wavelength'].d # Angstroms
- v_doppler = thermal_b.in_cgs().d # cm / s
- cdens = column_density.d
- dlambda = delta_lambda.d # Angstroms
- vlos = field_data['velocity_los'].in_units("km/s").d
-
+ # for a given transition, step through each location in the
+ # observed spectrum where it occurs and deposit a voigt profile
for i, lixel in parallel_objects(enumerate(valid_lines), njobs=-1):
- my_bin_ratio = spectrum_bin_ratio
+ my_vbin_window_width = vbin_window_width[i]
+ my_n_vbins = n_vbins
+ my_vbin_width = vbin_width[i]
while True:
- lambda_bins, line_tau = \
+ vbins = \
+ np.linspace(lambda_1[i]-my_vbin_window_width/2.,
+ lambda_1[i]+my_vbin_window_width/2.,
+ my_n_vbins, endpoint=False)
+
+ vbins, vtau = \
tau_profile(
- lambda_0, line['f_value'], line['gamma'], v_doppler[lixel],
- cdens[lixel], delta_lambda=dlambda[lixel],
- lambda_bins=lbins[left_index[lixel]:right_index[lixel]])
+ lambda_0, line['f_value'], line['gamma'], thermb[i],
+ cdens[i], delta_lambda=dlambda[i],
+ lambda_bins=vbins)
- # Widen wavelength window until optical depth reaches a max value at the ends.
- if (line_tau[0] < min_tau and line_tau[-1] < min_tau) or \
- (left_index[lixel] <= 0 and right_index[lixel] >= self.n_lambda):
+ # If tau has not dropped below min tau threshold by the
+ # edges (ie the wings), then widen the wavelength
+ # window and repeat process.
+ if (vtau[0] < min_tau and vtau[-1] < min_tau):
break
- my_bin_ratio *= 2
- left_index[lixel] = (center_bins[lixel] -
- my_bin_ratio *
- width_ratio[lixel]).astype(int).clip(0, self.n_lambda)
- right_index[lixel] = (center_bins[lixel] +
- my_bin_ratio *
- width_ratio[lixel]).astype(int).clip(0, self.n_lambda)
+ my_vbin_window_width *= 2
+ my_n_vbins *= 2
- self.tau_field[left_index[lixel]:right_index[lixel]] += line_tau
- if save_line_list and line['label_threshold'] is not None and \
- cdens[lixel] >= line['label_threshold']:
- if use_peculiar_velocity:
- peculiar_velocity = vlos[lixel]
- else:
- peculiar_velocity = 0.0
- self.spectrum_line_list.append({'label': line['label'],
- 'wavelength': (lambda_0 + dlambda[lixel]),
- 'column_density': column_density[lixel],
- 'b_thermal': thermal_b[lixel],
- 'redshift': field_data['redshift'][lixel],
- 'v_pec': peculiar_velocity})
+ # run digitize to identify which vbins are deposited into which
+ # global lambda bins.
+ # shift global lambda bins over by half a bin width;
+ # this has the effect of assuring np.digitize will place
+ # the vbins in the closest bin center.
+ binned = np.digitize(vbins,
+ self.lambda_bins + (0.5 * self.bin_width))
+
+ # numerically integrate the virtual bins to calculate a
+ # virtual equivalent width; then sum the virtual equivalent
+ # widths and deposit into each spectral bin
+ vEW = vtau * my_vbin_width
+ EW = [vEW[binned == j].sum() for j in range(len(self.lambda_bins))]
+ EW = np.array(EW)/self.bin_width.d
+ self.tau_field += EW
+
+ #self.tau_field[left_index[lixel]:right_index[lixel]] += line_tau
+ #if save_line_list and line['label_threshold'] is not None and \
+ # cdens[lixel] >= line['label_threshold']:
+ # if use_peculiar_velocity:
+ # peculiar_velocity = vlos[lixel]
+ # else:
+ # peculiar_velocity = 0.0
+ # self.spectrum_line_list.append({'label': line['label'],
+ # 'wavelength': (lambda_0 + dlambda[lixel]),
+ # 'column_density': column_density[lixel],
+ # 'b_thermal': thermal_b[lixel],
+ # 'redshift': field_data['redshift'][lixel],
+ # 'v_pec': peculiar_velocity})
pbar.update(i)
pbar.finish()
del column_density, delta_lambda, thermal_b, \
- center_bins, width_ratio, left_index, right_index
+ thermal_width, vbins, vtau, vbin_width, n_vbins, \
+ vbin_window_width
comm = _get_comm(())
self.tau_field = comm.mpi_allreduce(self.tau_field, op="sum")
https://bitbucket.org/yt_analysis/yt/commits/3a1cce502200/
Changeset: 3a1cce502200
Branch: yt
User: chummels
Date: 2015-11-13 00:43:41+00:00
Summary: Speeding up absorption spectrum by decreasing size of window to apply virtual voigt profile.
Affected #: 1 file
diff -r 4278eff7a0edf7e474092313445e9ebc9efbad91 -r 3a1cce5022002ff3597e1edfd20cd44684cc6456 yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -257,6 +257,9 @@
delta_lambda = line['wavelength'] * field_data['redshift']
# lambda_obs is central wavelength of line after redshift
lambda_obs = line['wavelength'] + delta_lambda
+ # bin index in lambda_bins of central wavelength of line after z
+ center_index = np.digitize(lambda_obs, self.lambda_bins)
+
# thermal broadening b parameter
thermal_b = np.sqrt((2 * boltzmann_constant_cgs *
field_data['temperature']) /
@@ -291,11 +294,11 @@
# --the actual spectral bin width
# to start, the window over which the voigt profile will be
# deposited will be 50 v_bins, but this may grow if necessary
- vbin_width = np.amin(zip(thermal_width / 10.,
+ vbin_width = np.amin(zip(thermal_width / 100.,
self.bin_width *
np.ones(len(thermal_width))),
axis=1)
- n_vbins = 50
+ n_vbins = 500
vbin_window_width = n_vbins*vbin_width
valid_lines = np.arange(len(thermal_width))
@@ -329,21 +332,30 @@
my_vbin_window_width *= 2
my_n_vbins *= 2
+ # identify the extrema of the vbin_window so as to speed
+ # up searching over the entire lambda_bins array
+ bins_from_center = np.ceil((my_vbin_window_width/2.) / \
+ self.bin_width.d) + 1
+ left_index = (center_index[i] - bins_from_center).clip(0, self.n_lambda)
+ right_index = (center_index[i] + bins_from_center).clip(0, self.n_lambda)
+ window_width = right_index - left_index
+
# run digitize to identify which vbins are deposited into which
# global lambda bins.
# shift global lambda bins over by half a bin width;
# this has the effect of assuring np.digitize will place
# the vbins in the closest bin center.
binned = np.digitize(vbins,
- self.lambda_bins + (0.5 * self.bin_width))
+ self.lambda_bins[left_index:right_index] \
+ + (0.5 * self.bin_width))
# numerically integrate the virtual bins to calculate a
# virtual equivalent width; then sum the virtual equivalent
# widths and deposit into each spectral bin
vEW = vtau * my_vbin_width
- EW = [vEW[binned == j].sum() for j in range(len(self.lambda_bins))]
+ EW = [vEW[binned == j].sum() for j in np.arange(window_width)]
EW = np.array(EW)/self.bin_width.d
- self.tau_field += EW
+ self.tau_field[left_index:right_index] += EW
#self.tau_field[left_index[lixel]:right_index[lixel]] += line_tau
#if save_line_list and line['label_threshold'] is not None and \
https://bitbucket.org/yt_analysis/yt/commits/e3771adcaeea/
Changeset: e3771adcaeea
Branch: yt
User: chummels
Date: 2015-11-13 01:38:06+00:00
Summary: Fixing bug causing spikes in overresolved lines.
Affected #: 1 file
diff -r 3a1cce5022002ff3597e1edfd20cd44684cc6456 -r e3771adcaeeae311ee74a27cb633b3e2443cbb8c yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -289,16 +289,20 @@
# profile into them, then numerically integrate their tau values
# and sum them to redeposit them into the actual spectral bins.
- # virtual bins will be the smaller of 2 options:
- # --1/10th of the thermal width
- # --the actual spectral bin width
- # to start, the window over which the voigt profile will be
- # deposited will be 50 v_bins, but this may grow if necessary
- vbin_width = np.amin(zip(thermal_width / 100.,
- self.bin_width *
- np.ones(len(thermal_width))),
- axis=1)
- n_vbins = 500
+ # virtual bins (vbins) will be:
+ # 1) <= the bin_width; assures at least as good as spectral bins
+ # 2) <= 1/10th the thermal width; assures resolving voigt profiles
+ # 3) a bin width will be divisible by vbin_width times a power of
+ # 10; this will assure we don't get spikes in the deposited
+ # spectra from uneven numbers of vbins per bin
+ resolution = thermal_width / self.bin_width
+ vbin_width = self.bin_width / \
+ 10**(np.ceil(np.log10(10/resolution)).clip(0, np.inf))
+ vbin_width = vbin_width.in_units('angstrom').d
+
+ # the virtual window into which the line is deposited initially
+ # spans a region of 5 thermal_widths, but this may expand
+ n_vbins = np.ceil(5*thermal_width.d/vbin_width)
vbin_window_width = n_vbins*vbin_width
valid_lines = np.arange(len(thermal_width))
@@ -309,7 +313,7 @@
# observed spectrum where it occurs and deposit a voigt profile
for i, lixel in parallel_objects(enumerate(valid_lines), njobs=-1):
my_vbin_window_width = vbin_window_width[i]
- my_n_vbins = n_vbins
+ my_n_vbins = n_vbins[i]
my_vbin_width = vbin_width[i]
while True:
https://bitbucket.org/yt_analysis/yt/commits/9c4ef612cd4a/
Changeset: 9c4ef612cd4a
Branch: yt
User: chummels
Date: 2015-11-13 01:58:45+00:00
Summary: Cleaning up the changes i've made and added a few more comments to absorption spectrum.
Affected #: 1 file
diff -r e3771adcaeeae311ee74a27cb633b3e2443cbb8c -r 9c4ef612cd4ad5dff7aa296741627e2524245740 yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -269,11 +269,6 @@
thermal_width = (lambda_obs * thermal_b /
speed_of_light_cgs).convert_to_units("angstrom")
- if (thermal_width < self.bin_width).any():
- mylog.warn(("%d out of %d line components are unresolved.") %
- ((thermal_width < self.bin_width).sum(),
- thermal_width.size))
-
# Sanitize units for faster runtime of the tau_profile machinery.
lambda_0 = line['wavelength'].d # line's rest frame; angstroms
lambda_1 = lambda_obs.d # line's observed frame; angstroms
@@ -305,13 +300,18 @@
n_vbins = np.ceil(5*thermal_width.d/vbin_width)
vbin_window_width = n_vbins*vbin_width
+ if (thermal_width < self.bin_width).any():
+ mylog.warn(("%d out of %d line components are unresolved.") %
+ ((thermal_width < self.bin_width).sum(),
+ thermal_width.size))
+
valid_lines = np.arange(len(thermal_width))
pbar = get_pbar("Adding line - %s [%f A]: " % (line['label'], line['wavelength']),
thermal_width.size)
# for a given transition, step through each location in the
# observed spectrum where it occurs and deposit a voigt profile
- for i, lixel in parallel_objects(enumerate(valid_lines), njobs=-1):
+ for i in parallel_objects(valid_lines, njobs=-1):
my_vbin_window_width = vbin_window_width[i]
my_n_vbins = n_vbins[i]
my_vbin_width = vbin_width[i]
@@ -361,25 +361,25 @@
EW = np.array(EW)/self.bin_width.d
self.tau_field[left_index:right_index] += EW
- #self.tau_field[left_index[lixel]:right_index[lixel]] += line_tau
- #if save_line_list and line['label_threshold'] is not None and \
- # cdens[lixel] >= line['label_threshold']:
- # if use_peculiar_velocity:
- # peculiar_velocity = vlos[lixel]
- # else:
- # peculiar_velocity = 0.0
- # self.spectrum_line_list.append({'label': line['label'],
- # 'wavelength': (lambda_0 + dlambda[lixel]),
- # 'column_density': column_density[lixel],
- # 'b_thermal': thermal_b[lixel],
- # 'redshift': field_data['redshift'][lixel],
- # 'v_pec': peculiar_velocity})
+ if save_line_list and line['label_threshold'] is not None and \
+ cdens[i] >= line['label_threshold']:
+ if use_peculiar_velocity:
+ peculiar_velocity = vlos[i]
+ else:
+ peculiar_velocity = 0.0
+ self.spectrum_line_list.append({'label': line['label'],
+ 'wavelength': (lambda_0 + dlambda[i]),
+ 'column_density': column_density[i],
+ 'b_thermal': thermal_b[i],
+ 'redshift': field_data['redshift'][i],
+ 'v_pec': peculiar_velocity})
pbar.update(i)
pbar.finish()
- del column_density, delta_lambda, thermal_b, \
- thermal_width, vbins, vtau, vbin_width, n_vbins, \
- vbin_window_width
+ del column_density, delta_lambda, lambda_obs, center_index, \
+ thermal_b, thermal_width, lambda_1, cdens, thermb, dlambda, \
+ vlos, resolution, vbin_width, n_vbins, vbin_window_width, \
+ valid_lines, vbins, vtau, vEW
comm = _get_comm(())
self.tau_field = comm.mpi_allreduce(self.tau_field, op="sum")
https://bitbucket.org/yt_analysis/yt/commits/ca15d8b6d74c/
Changeset: ca15d8b6d74c
Branch: yt
User: chummels
Date: 2015-11-14 00:31:52+00:00
Summary: Adding equivalent width conservation test to AbsorptionSpectrum analysis module
Affected #: 1 file
diff -r 9c4ef612cd4ad5dff7aa296741627e2524245740 -r ca15d8b6d74c9868d4a1934e918001e9880b15f6 yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
@@ -12,7 +12,7 @@
import numpy as np
from yt.testing import \
- assert_allclose_units, requires_file, requires_module
+ assert_allclose_units, requires_file, requires_module, assert_almost_equal
from yt.analysis_modules.absorption_spectrum.absorption_line import \
voigt_old, voigt_scipy
from yt.analysis_modules.absorption_spectrum.api import AbsorptionSpectrum
@@ -22,10 +22,10 @@
import shutil
COSMO_PLUS = "enzo_cosmology_plus/AMRCosmology.enzo"
-
+COSMO_PLUS_SINGLE = "enzo_cosmology_plus/RD0009/RD0009"
@requires_file(COSMO_PLUS)
-def test_absorption_spectrum():
+def test_absorption_spectrum_cosmo():
"""
This test is simply following the description in the docs for how to
generate an absorption spectrum from a cosmological light ray for one
@@ -37,7 +37,7 @@
curdir = os.getcwd()
os.chdir(tmpdir)
- lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.1)
+ lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.03)
lr.make_light_ray(seed=1234567,
fields=['temperature', 'density', 'H_number_density'],
@@ -73,6 +73,105 @@
os.chdir(curdir)
shutil.rmtree(tmpdir)
+ at requires_file(COSMO_PLUS_SINGLE)
+def test_absorption_spectrum_non_cosmo():
+ """
+ This test is simply following the description in the docs for how to
+ generate an absorption spectrum from a non-cosmological light ray for one
+ of the sample datasets
+ """
+
+ # Set up in a temp dir
+ tmpdir = tempfile.mkdtemp()
+ curdir = os.getcwd()
+ os.chdir(tmpdir)
+
+ lr = LightRay(COSMO_PLUS_SINGLE)
+
+ ray_start = [0,0,0]
+ ray_end = [1,1,1]
+ lr.make_light_ray(start_position=ray_start, end_position=ray_end,
+ fields=['temperature', 'density', 'H_number_density'],
+ data_filename='lightray.h5')
+
+ sp = AbsorptionSpectrum(900.0, 1800.0, 10000)
+
+ my_label = 'HI Lya'
+ field = 'H_number_density'
+ wavelength = 1215.6700 # Angstromss
+ f_value = 4.164E-01
+ gamma = 6.265e+08
+ mass = 1.00794
+
+ sp.add_line(my_label, field, wavelength, f_value,
+ gamma, mass, label_threshold=1.e10)
+
+ my_label = 'HI Lya'
+ field = 'H_number_density'
+ wavelength = 912.323660 # Angstroms
+ normalization = 1.6e17
+ index = 3.0
+
+ sp.add_continuum(my_label, field, wavelength, normalization, index)
+
+ wavelength, flux = sp.make_spectrum('lightray.h5',
+ output_file='spectrum.txt',
+ line_list_file='lines.txt',
+ use_peculiar_velocity=True)
+
+ # clean up
+ os.chdir(curdir)
+ shutil.rmtree(tmpdir)
+
+ at requires_file(COSMO_PLUS_SINGLE)
+def test_equivalent_width_conserved():
+ """
+ This test assures that the equivalent width of the optical depth
+ is conserved regardless of the bin width employed in wavelength space.
+ Unresolved lines should still deposit optical depth into the spectrum.
+ """
+
+ # Set up in a temp dir
+ tmpdir = tempfile.mkdtemp()
+ curdir = os.getcwd()
+ os.chdir(tmpdir)
+
+ lr = LightRay(COSMO_PLUS_SINGLE)
+
+ ray_start = [0,0,0]
+ ray_end = [1,1,1]
+ lr.make_light_ray(start_position=ray_start, end_position=ray_end,
+ fields=['temperature', 'density', 'H_number_density'],
+ data_filename='lightray.h5')
+
+ my_label = 'HI Lya'
+ field = 'H_number_density'
+ wave = 1215.6700 # Angstromss
+ f_value = 4.164E-01
+ gamma = 6.265e+08
+ mass = 1.00794
+
+ lambda_min= 1200
+ lambda_max= 1300
+ lambda_bin_widths = [1e-3, 1e-2, 1e-1, 1e0, 1e1]
+ total_tau = []
+
+ for lambda_bin_width in lambda_bin_widths:
+ n_lambda = ((lambda_max - lambda_min)/ lambda_bin_width) + 1
+ sp = AbsorptionSpectrum(lambda_min=lambda_min, lambda_max=lambda_max,
+ n_lambda=n_lambda)
+ sp.add_line(my_label, field, wave, f_value, gamma, mass)
+ wavelength, flux = sp.make_spectrum('ray.h5')
+ total_tau.append((lambda_bin_width * sp.tau_field).sum())
+
+ # assure that the total tau values are all within 1e-5 of each other
+ for tau in total_tau:
+ assert_almost_equal(tau, total_tau[0], 5)
+
+ # clean up
+ os.chdir(curdir)
+ shutil.rmtree(tmpdir)
+
@requires_file(COSMO_PLUS)
@requires_module("astropy")
https://bitbucket.org/yt_analysis/yt/commits/b5fa21bb2db6/
Changeset: b5fa21bb2db6
Branch: yt
User: chummels
Date: 2015-11-14 06:28:17+00:00
Summary: Adding answer testing to AbsorptionSpectrum analysis module.
Affected #: 1 file
diff -r ca15d8b6d74c9868d4a1934e918001e9880b15f6 -r b5fa21bb2db6d533de2cdeb97d1634d25f86a43e yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
@@ -12,17 +12,24 @@
import numpy as np
from yt.testing import \
- assert_allclose_units, requires_file, requires_module, assert_almost_equal
+ assert_allclose_units, requires_file, requires_module, \
+ assert_almost_equal, assert_array_almost_equal
from yt.analysis_modules.absorption_spectrum.absorption_line import \
voigt_old, voigt_scipy
from yt.analysis_modules.absorption_spectrum.api import AbsorptionSpectrum
from yt.analysis_modules.cosmological_observation.api import LightRay
+from yt.config import ytcfg
import tempfile
import os
import shutil
+import h5py as h5
+
+test_dir = ytcfg.get("yt", "test_data_dir")
COSMO_PLUS = "enzo_cosmology_plus/AMRCosmology.enzo"
COSMO_PLUS_SINGLE = "enzo_cosmology_plus/RD0009/RD0009"
+HI_SPECTRUM = "absorption_spectrum_data/enzo_lyman_alpha_spec.h5"
+HI_SPECTRUM_FILE = os.path.join(test_dir, HI_SPECTRUM)
@requires_file(COSMO_PLUS)
def test_absorption_spectrum_cosmo():
@@ -74,6 +81,7 @@
shutil.rmtree(tmpdir)
@requires_file(COSMO_PLUS_SINGLE)
+ at requires_file(HI_SPECTRUM)
def test_absorption_spectrum_non_cosmo():
"""
This test is simply following the description in the docs for how to
@@ -94,7 +102,7 @@
fields=['temperature', 'density', 'H_number_density'],
data_filename='lightray.h5')
- sp = AbsorptionSpectrum(900.0, 1800.0, 10000)
+ sp = AbsorptionSpectrum(1200.0, 1300.0, 10001)
my_label = 'HI Lya'
field = 'H_number_density'
@@ -106,19 +114,22 @@
sp.add_line(my_label, field, wavelength, f_value,
gamma, mass, label_threshold=1.e10)
- my_label = 'HI Lya'
- field = 'H_number_density'
- wavelength = 912.323660 # Angstroms
- normalization = 1.6e17
- index = 3.0
-
- sp.add_continuum(my_label, field, wavelength, normalization, index)
-
wavelength, flux = sp.make_spectrum('lightray.h5',
- output_file='spectrum.txt',
+ output_file='spectrum.h5',
line_list_file='lines.txt',
use_peculiar_velocity=True)
+ # load just-generated hdf5 file of spectral data (for consistency)
+ f_new = h5.File('spectrum.h5', 'r')
+
+ # load standard data for comparison
+ f_old = h5.File(HI_SPECTRUM_FILE, 'r')
+
+ # compare between standard data and current data for each array saved
+ # (wavelength, flux, tau)
+ for key in f_old.keys():
+ assert_array_almost_equal(f_new[key].value, f_old[key].value, 10)
+
# clean up
os.chdir(curdir)
shutil.rmtree(tmpdir)
@@ -187,7 +198,7 @@
curdir = os.getcwd()
os.chdir(tmpdir)
- lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.1)
+ lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.03)
lr.make_light_ray(seed=1234567,
fields=['temperature', 'density', 'H_number_density'],
https://bitbucket.org/yt_analysis/yt/commits/d1bb9143fe95/
Changeset: d1bb9143fe95
Branch: yt
User: chummels
Date: 2015-11-14 19:17:07+00:00
Summary: Fixing bug in tests.
Affected #: 1 file
diff -r b5fa21bb2db6d533de2cdeb97d1634d25f86a43e -r d1bb9143fe95f0535b66d3559d5af5a33eed91c7 yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
@@ -172,7 +172,7 @@
sp = AbsorptionSpectrum(lambda_min=lambda_min, lambda_max=lambda_max,
n_lambda=n_lambda)
sp.add_line(my_label, field, wave, f_value, gamma, mass)
- wavelength, flux = sp.make_spectrum('ray.h5')
+ wavelength, flux = sp.make_spectrum('lightray.h5')
total_tau.append((lambda_bin_width * sp.tau_field).sum())
# assure that the total tau values are all within 1e-5 of each other
https://bitbucket.org/yt_analysis/yt/commits/67d46907ab6a/
Changeset: 67d46907ab6a
Branch: yt
User: chummels
Date: 2015-11-16 15:36:25+00:00
Summary: Changing message when AbsorptionSpectrum deposits unresolved lines.
Affected #: 1 file
diff -r d1bb9143fe95f0535b66d3559d5af5a33eed91c7 -r 67d46907ab6a7d8d85a8959a0123bd2c0adf1162 yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -301,12 +301,14 @@
vbin_window_width = n_vbins*vbin_width
if (thermal_width < self.bin_width).any():
- mylog.warn(("%d out of %d line components are unresolved.") %
+ mylog.info(("%d out of %d line components will be " + \
+ "deposited as unresolved lines.") %
((thermal_width < self.bin_width).sum(),
thermal_width.size))
valid_lines = np.arange(len(thermal_width))
- pbar = get_pbar("Adding line - %s [%f A]: " % (line['label'], line['wavelength']),
+ pbar = get_pbar("Adding line - %s [%f A]: " % \
+ (line['label'], line['wavelength']),
thermal_width.size)
# for a given transition, step through each location in the
https://bitbucket.org/yt_analysis/yt/commits/b0285115281d/
Changeset: b0285115281d
Branch: yt
User: chummels
Date: 2015-11-16 15:43:53+00:00
Summary: Adding import on demand to h5py import in AbsorptionSpectrum test
Affected #: 1 file
diff -r 67d46907ab6a7d8d85a8959a0123bd2c0adf1162 -r b0285115281d2e55f116c2c9679e2913500cf490 yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
@@ -22,7 +22,8 @@
import tempfile
import os
import shutil
-import h5py as h5
+from yt.utilities.on_demand_imports import \
+ _h5py as h5
test_dir = ytcfg.get("yt", "test_data_dir")
https://bitbucket.org/yt_analysis/yt/commits/042ef88323f3/
Changeset: 042ef88323f3
Branch: yt
User: chummels
Date: 2015-11-16 15:45:22+00:00
Summary: Setting min_tau values to match for AbsorptionSpectrum add_continua and add_line
Affected #: 1 file
diff -r b0285115281d2e55f116c2c9679e2913500cf490 -r 042ef88323f3085fb4e5f0d20247990e3005ab97 yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -203,7 +203,7 @@
Add continuum features to the spectrum.
"""
# Only add continuum features down to tau of 1.e-4.
- min_tau = 1.e-4
+ min_tau = 1.e-3
for continuum in self.continuum_list:
column_density = field_data[continuum['field_name']] * field_data['dl']
@@ -241,7 +241,7 @@
# Widen wavelength window until optical depth falls below this tau
# value at the ends to assure that the wings of a line have been
# fully resolved.
- min_tau = 0.001
+ min_tau = 1e-3
# step through each ionic transition (e.g. HI, HII, MgII) specified
# and deposit the lines into the spectrum
https://bitbucket.org/yt_analysis/yt/commits/2aff455a4323/
Changeset: 2aff455a4323
Branch: yt
User: chummels
Date: 2015-11-16 16:52:03+00:00
Summary: Making subgrid_resolution a keyword in make_spectrum() instead of assuming a value of 10.
Affected #: 1 file
diff -r 042ef88323f3085fb4e5f0d20247990e3005ab97 -r 2aff455a4323785a1bae533e38831d948915c59b yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -115,7 +115,8 @@
def make_spectrum(self, input_file, output_file="spectrum.h5",
line_list_file="lines.txt",
- use_peculiar_velocity=True, njobs="auto"):
+ use_peculiar_velocity=True,
+ subgrid_resolution=10, njobs="auto"):
"""
Make spectrum from ray data using the line list.
@@ -140,6 +141,17 @@
use_peculiar_velocity : optional, bool
if True, include line of sight velocity for shifting lines.
Default: True
+ subgrid_resolution : optional, int
+ When a line is being added that is unresolved (ie its thermal
+ width is less than the spectral bin width), the voigt profile of
+ the line is deposited into an array of virtual bins at higher
+ resolution. The optical depth from these virtual bins is integrated
+ and then added to the coarser spectral bin. The subgrid_resolution
+ value determines the ratio between the thermal width and the
+ bin width of the virtual bins. Increasing this value yields smaller
+ virtual bins, which increases accuracy, but is more expensive.
+ A value of 10 yields accuracy to the 4th significant digit.
+ Default: 10
njobs : optional, int or "auto"
the number of process groups into which the loop over
absorption lines will be divided. If set to -1, each
@@ -181,7 +193,9 @@
njobs = min(comm.size, len(self.line_list))
self._add_lines_to_spectrum(field_data, use_peculiar_velocity,
- line_list_file is not None, njobs=njobs)
+ line_list_file is not None,
+ subgrid_resolution=subgrid_resolution,
+ njobs=njobs)
self._add_continua_to_spectrum(field_data, use_peculiar_velocity)
self.flux_field = np.exp(-self.tau_field)
@@ -234,7 +248,7 @@
pbar.finish()
def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,
- save_line_list, njobs=-1):
+ save_line_list, subgrid_resolution=10, njobs=-1):
"""
Add the absorption lines to the spectrum.
"""
@@ -287,12 +301,13 @@
# virtual bins (vbins) will be:
# 1) <= the bin_width; assures at least as good as spectral bins
# 2) <= 1/10th the thermal width; assures resolving voigt profiles
+ # (actually 1/subgrid_resolution value, default is 1/10)
# 3) a bin width will be divisible by vbin_width times a power of
# 10; this will assure we don't get spikes in the deposited
# spectra from uneven numbers of vbins per bin
resolution = thermal_width / self.bin_width
vbin_width = self.bin_width / \
- 10**(np.ceil(np.log10(10/resolution)).clip(0, np.inf))
+ 10**(np.ceil(np.log10(subgrid_resolution/resolution)).clip(0, np.inf))
vbin_width = vbin_width.in_units('angstrom').d
# the virtual window into which the line is deposited initially
https://bitbucket.org/yt_analysis/yt/commits/de13fe315e19/
Changeset: de13fe315e19
Branch: yt
User: chummels
Date: 2015-11-16 17:06:05+00:00
Summary: Adding "answer" test to cosmological light ray absorption spectrum test
Affected #: 1 file
diff -r 2aff455a4323785a1bae533e38831d948915c59b -r de13fe315e19d583967c48f936c938fb3a32c1d8 yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
@@ -29,15 +29,16 @@
COSMO_PLUS = "enzo_cosmology_plus/AMRCosmology.enzo"
COSMO_PLUS_SINGLE = "enzo_cosmology_plus/RD0009/RD0009"
+HI_SPECTRUM_COSMO = "absorption_spectrum_data/enzo_lyman_alpha_cosmo_spec.h5"
+HI_SPECTRUM_COSMO_FILE = os.path.join(test_dir, HI_SPECTRUM_COSMO)
HI_SPECTRUM = "absorption_spectrum_data/enzo_lyman_alpha_spec.h5"
HI_SPECTRUM_FILE = os.path.join(test_dir, HI_SPECTRUM)
@requires_file(COSMO_PLUS)
+ at requires_file(HI_SPECTRUM_COSMO)
def test_absorption_spectrum_cosmo():
"""
- This test is simply following the description in the docs for how to
- generate an absorption spectrum from a cosmological light ray for one
- of the sample datasets
+ This test generates an absorption spectrum from a cosmological light ray
"""
# Set up in a temp dir
@@ -73,10 +74,21 @@
sp.add_continuum(my_label, field, wavelength, normalization, index)
wavelength, flux = sp.make_spectrum('lightray.h5',
- output_file='spectrum.txt',
+ output_file='spectrum.h5',
line_list_file='lines.txt',
use_peculiar_velocity=True)
+ # load just-generated hdf5 file of spectral data (for consistency)
+ f_new = h5.File('spectrum.h5', 'r')
+
+ # load standard data for comparison
+ f_old = h5.File(HI_SPECTRUM_COSMO_FILE, 'r')
+
+ # compare between standard data and current data for each array saved
+ # (wavelength, flux, tau)
+ for key in f_old.keys():
+ assert_array_almost_equal(f_new[key].value, f_old[key].value, 10)
+
# clean up
os.chdir(curdir)
shutil.rmtree(tmpdir)
@@ -85,9 +97,7 @@
@requires_file(HI_SPECTRUM)
def test_absorption_spectrum_non_cosmo():
"""
- This test is simply following the description in the docs for how to
- generate an absorption spectrum from a non-cosmological light ray for one
- of the sample datasets
+ This test generates an absorption spectrum from a non-cosmological light ray
"""
# Set up in a temp dir
@@ -138,8 +148,8 @@
@requires_file(COSMO_PLUS_SINGLE)
def test_equivalent_width_conserved():
"""
- This test assures that the equivalent width of the optical depth
- is conserved regardless of the bin width employed in wavelength space.
+ This tests that the equivalent width of the optical depth is conserved
+ regardless of the bin width employed in wavelength space.
Unresolved lines should still deposit optical depth into the spectrum.
"""
@@ -185,13 +195,11 @@
shutil.rmtree(tmpdir)
- at requires_file(COSMO_PLUS)
+ at requires_file(COSMO_PLUS_SINGLE)
@requires_module("astropy")
def test_absorption_spectrum_fits():
"""
- This test is simply following the description in the docs for how to
- generate an absorption spectrum from a cosmological light ray for one
- of the sample datasets. Outputs to fits file if astropy is installed.
+ This test generates an absorption spectrum and saves it as a fits file.
"""
# Set up in a temp dir
@@ -199,11 +207,12 @@
curdir = os.getcwd()
os.chdir(tmpdir)
- lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.03)
+ lr = LightRay(COSMO_PLUS_SINGLE)
- lr.make_light_ray(seed=1234567,
+ ray_start = [0,0,0]
+ ray_end = [1,1,1]
+ lr.make_light_ray(start_position=ray_start, end_position=ray_end,
fields=['temperature', 'density', 'H_number_density'],
- get_los_velocity=True,
data_filename='lightray.h5')
sp = AbsorptionSpectrum(900.0, 1800.0, 10000)
https://bitbucket.org/yt_analysis/yt/commits/90f900be7a36/
Changeset: 90f900be7a36
Branch: yt
User: ngoldbaum
Date: 2015-11-16 21:27:16+00:00
Summary: Merged in chummels/yt (pull request #1858)
Enabling AbsorptionSpectrum to deposit unresolved spectral lines
Affected #: 2 files
diff -r 50f82a593a983ca8ee997a0650fc03eddfdbb30b -r 90f900be7a36433fdd48941cae4bc91066ff5c76 yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -116,7 +116,8 @@
def make_spectrum(self, input_file, output_file="spectrum.h5",
line_list_file="lines.txt",
- use_peculiar_velocity=True, njobs="auto"):
+ use_peculiar_velocity=True,
+ subgrid_resolution=10, njobs="auto"):
"""
Make spectrum from ray data using the line list.
@@ -141,6 +142,17 @@
use_peculiar_velocity : optional, bool
if True, include line of sight velocity for shifting lines.
Default: True
+ subgrid_resolution : optional, int
+ When a line is being added that is unresolved (ie its thermal
+ width is less than the spectral bin width), the voigt profile of
+ the line is deposited into an array of virtual bins at higher
+ resolution. The optical depth from these virtual bins is integrated
+ and then added to the coarser spectral bin. The subgrid_resolution
+ value determines the ratio between the thermal width and the
+ bin width of the virtual bins. Increasing this value yields smaller
+ virtual bins, which increases accuracy, but is more expensive.
+ A value of 10 yields accuracy to the 4th significant digit.
+ Default: 10
njobs : optional, int or "auto"
the number of process groups into which the loop over
absorption lines will be divided. If set to -1, each
@@ -182,7 +194,9 @@
njobs = min(comm.size, len(self.line_list))
self._add_lines_to_spectrum(field_data, use_peculiar_velocity,
- line_list_file is not None, njobs=njobs)
+ line_list_file is not None,
+ subgrid_resolution=subgrid_resolution,
+ njobs=njobs)
self._add_continua_to_spectrum(field_data, use_peculiar_velocity)
self.flux_field = np.exp(-self.tau_field)
@@ -204,7 +218,7 @@
Add continuum features to the spectrum.
"""
# Only add continuum features down to tau of 1.e-4.
- tau_min = 1.e-4
+ min_tau = 1.e-3
for continuum in self.continuum_list:
column_density = field_data[continuum['field_name']] * field_data['dl']
@@ -216,12 +230,12 @@
this_wavelength = delta_lambda + continuum['wavelength']
right_index = np.digitize(this_wavelength, self.lambda_bins).clip(0, self.n_lambda)
left_index = np.digitize((this_wavelength *
- np.power((tau_min * continuum['normalization'] /
+ np.power((min_tau * continuum['normalization'] /
column_density), (1. / continuum['index']))),
self.lambda_bins).clip(0, self.n_lambda)
valid_continuua = np.where(((column_density /
- continuum['normalization']) > tau_min) &
+ continuum['normalization']) > min_tau) &
(right_index - left_index > 1))[0]
pbar = get_pbar("Adding continuum feature - %s [%f A]: " % \
(continuum['label'], continuum['wavelength']),
@@ -235,98 +249,155 @@
pbar.finish()
def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,
- save_line_list, njobs=-1):
+ save_line_list, subgrid_resolution=10, njobs=-1):
"""
Add the absorption lines to the spectrum.
"""
- # Only make voigt profile for slice of spectrum that is 10 times the line width.
- spectrum_bin_ratio = 5
- # Widen wavelength window until optical depth reaches a max value at the ends.
- max_tau = 0.001
+ # Widen wavelength window until optical depth falls below this tau
+ # value at the ends to assure that the wings of a line have been
+ # fully resolved.
+ min_tau = 1e-3
+ # step through each ionic transition (e.g. HI, HII, MgII) specified
+ # and deposit the lines into the spectrum
for line in parallel_objects(self.line_list, njobs=njobs):
column_density = field_data[line['field_name']] * field_data['dl']
+
# redshift_eff field combines cosmological and velocity redshifts
+ # so delta_lambda gives the offset in angstroms from the rest frame
+ # wavelength to the observed wavelength of the transition
if use_peculiar_velocity:
delta_lambda = line['wavelength'] * field_data['redshift_eff']
else:
delta_lambda = line['wavelength'] * field_data['redshift']
+ # lambda_obs is central wavelength of line after redshift
+ lambda_obs = line['wavelength'] + delta_lambda
+ # bin index in lambda_bins of central wavelength of line after z
+ center_index = np.digitize(lambda_obs, self.lambda_bins)
+
+ # thermal broadening b parameter
thermal_b = np.sqrt((2 * boltzmann_constant_cgs *
field_data['temperature']) /
line['atomic_mass'])
- center_bins = np.digitize((delta_lambda + line['wavelength']),
- self.lambda_bins)
- # ratio of line width to bin width
- width_ratio = ((line['wavelength'] + delta_lambda) * \
- thermal_b / speed_of_light_cgs / self.bin_width).in_units("").d
+ # the actual thermal width of the lines
+ thermal_width = (lambda_obs * thermal_b /
+ speed_of_light_cgs).convert_to_units("angstrom")
- if (width_ratio < 1.0).any():
- mylog.warn(("%d out of %d line components are unresolved, " +
- "consider increasing spectral resolution.") %
- ((width_ratio < 1.0).sum(), width_ratio.size))
+ # Sanitize units for faster runtime of the tau_profile machinery.
+ lambda_0 = line['wavelength'].d # line's rest frame; angstroms
+ lambda_1 = lambda_obs.d # line's observed frame; angstroms
+ cdens = column_density.in_units("cm**-2").d # cm**-2
+ thermb = thermal_b.in_cgs().d # thermal b coefficient; cm / s
+ dlambda = delta_lambda.d # lambda offset; angstroms
+ vlos = field_data['velocity_los'].in_units("km/s").d # km/s
- # do voigt profiles for a subset of the full spectrum
- left_index = (center_bins -
- spectrum_bin_ratio * width_ratio).astype(int).clip(0, self.n_lambda)
- right_index = (center_bins +
- spectrum_bin_ratio * width_ratio).astype(int).clip(0, self.n_lambda)
+ # When we actually deposit the voigt profile, sometimes we will
+ # have underresolved lines (ie lines with smaller widths than
+ # the spectral bin size). Here, we create virtual bins small
+ # enough in width to well resolve each line, deposit the voigt
+ # profile into them, then numerically integrate their tau values
+ # and sum them to redeposit them into the actual spectral bins.
- # loop over all lines wider than the bin width
- valid_lines = np.where((width_ratio >= 1.0) &
- (right_index - left_index > 1))[0]
- pbar = get_pbar("Adding line - %s [%f A]: " % (line['label'], line['wavelength']),
- valid_lines.size)
+ # virtual bins (vbins) will be:
+ # 1) <= the bin_width; assures at least as good as spectral bins
+ # 2) <= 1/10th the thermal width; assures resolving voigt profiles
+ # (actually 1/subgrid_resolution value, default is 1/10)
+ # 3) a bin width will be divisible by vbin_width times a power of
+ # 10; this will assure we don't get spikes in the deposited
+ # spectra from uneven numbers of vbins per bin
+ resolution = thermal_width / self.bin_width
+ vbin_width = self.bin_width / \
+ 10**(np.ceil(np.log10(subgrid_resolution/resolution)).clip(0, np.inf))
+ vbin_width = vbin_width.in_units('angstrom').d
- # Sanitize units here
- column_density.convert_to_units("cm ** -2")
- lbins = self.lambda_bins.d # Angstroms
- lambda_0 = line['wavelength'].d # Angstroms
- v_doppler = thermal_b.in_cgs().d # cm / s
- cdens = column_density.d
- dlambda = delta_lambda.d # Angstroms
- vlos = field_data['velocity_los'].in_units("km/s").d
+ # the virtual window into which the line is deposited initially
+ # spans a region of 5 thermal_widths, but this may expand
+ n_vbins = np.ceil(5*thermal_width.d/vbin_width)
+ vbin_window_width = n_vbins*vbin_width
- for i, lixel in parallel_objects(enumerate(valid_lines), njobs=-1):
- my_bin_ratio = spectrum_bin_ratio
+ if (thermal_width < self.bin_width).any():
+ mylog.info(("%d out of %d line components will be " + \
+ "deposited as unresolved lines.") %
+ ((thermal_width < self.bin_width).sum(),
+ thermal_width.size))
+
+ valid_lines = np.arange(len(thermal_width))
+ pbar = get_pbar("Adding line - %s [%f A]: " % \
+ (line['label'], line['wavelength']),
+ thermal_width.size)
+
+ # for a given transition, step through each location in the
+ # observed spectrum where it occurs and deposit a voigt profile
+ for i in parallel_objects(valid_lines, njobs=-1):
+ my_vbin_window_width = vbin_window_width[i]
+ my_n_vbins = n_vbins[i]
+ my_vbin_width = vbin_width[i]
while True:
- lambda_bins, line_tau = \
+ vbins = \
+ np.linspace(lambda_1[i]-my_vbin_window_width/2.,
+ lambda_1[i]+my_vbin_window_width/2.,
+ my_n_vbins, endpoint=False)
+
+ vbins, vtau = \
tau_profile(
- lambda_0, line['f_value'], line['gamma'], v_doppler[lixel],
- cdens[lixel], delta_lambda=dlambda[lixel],
- lambda_bins=lbins[left_index[lixel]:right_index[lixel]])
+ lambda_0, line['f_value'], line['gamma'], thermb[i],
+ cdens[i], delta_lambda=dlambda[i],
+ lambda_bins=vbins)
- # Widen wavelength window until optical depth reaches a max value at the ends.
- if (line_tau[0] < max_tau and line_tau[-1] < max_tau) or \
- (left_index[lixel] <= 0 and right_index[lixel] >= self.n_lambda):
+ # If tau has not dropped below min tau threshold by the
+ # edges (ie the wings), then widen the wavelength
+ # window and repeat process.
+ if (vtau[0] < min_tau and vtau[-1] < min_tau):
break
- my_bin_ratio *= 2
- left_index[lixel] = (center_bins[lixel] -
- my_bin_ratio *
- width_ratio[lixel]).astype(int).clip(0, self.n_lambda)
- right_index[lixel] = (center_bins[lixel] +
- my_bin_ratio *
- width_ratio[lixel]).astype(int).clip(0, self.n_lambda)
+ my_vbin_window_width *= 2
+ my_n_vbins *= 2
- self.tau_field[left_index[lixel]:right_index[lixel]] += line_tau
+ # identify the extrema of the vbin_window so as to speed
+ # up searching over the entire lambda_bins array
+ bins_from_center = np.ceil((my_vbin_window_width/2.) / \
+ self.bin_width.d) + 1
+ left_index = (center_index[i] - bins_from_center).clip(0, self.n_lambda)
+ right_index = (center_index[i] + bins_from_center).clip(0, self.n_lambda)
+ window_width = right_index - left_index
+
+ # run digitize to identify which vbins are deposited into which
+ # global lambda bins.
+ # shift global lambda bins over by half a bin width;
+ # this has the effect of assuring np.digitize will place
+ # the vbins in the closest bin center.
+ binned = np.digitize(vbins,
+ self.lambda_bins[left_index:right_index] \
+ + (0.5 * self.bin_width))
+
+ # numerically integrate the virtual bins to calculate a
+ # virtual equivalent width; then sum the virtual equivalent
+ # widths and deposit into each spectral bin
+ vEW = vtau * my_vbin_width
+ EW = [vEW[binned == j].sum() for j in np.arange(window_width)]
+ EW = np.array(EW)/self.bin_width.d
+ self.tau_field[left_index:right_index] += EW
+
if save_line_list and line['label_threshold'] is not None and \
- cdens[lixel] >= line['label_threshold']:
+ cdens[i] >= line['label_threshold']:
if use_peculiar_velocity:
- peculiar_velocity = vlos[lixel]
+ peculiar_velocity = vlos[i]
else:
peculiar_velocity = 0.0
self.spectrum_line_list.append({'label': line['label'],
- 'wavelength': (lambda_0 + dlambda[lixel]),
- 'column_density': column_density[lixel],
- 'b_thermal': thermal_b[lixel],
- 'redshift': field_data['redshift'][lixel],
+ 'wavelength': (lambda_0 + dlambda[i]),
+ 'column_density': column_density[i],
+ 'b_thermal': thermal_b[i],
+ 'redshift': field_data['redshift'][i],
'v_pec': peculiar_velocity})
pbar.update(i)
pbar.finish()
- del column_density, delta_lambda, thermal_b, \
- center_bins, width_ratio, left_index, right_index
+ del column_density, delta_lambda, lambda_obs, center_index, \
+ thermal_b, thermal_width, lambda_1, cdens, thermb, dlambda, \
+ vlos, resolution, vbin_width, n_vbins, vbin_window_width, \
+ valid_lines, vbins, vtau, vEW
comm = _get_comm(())
self.tau_field = comm.mpi_allreduce(self.tau_field, op="sum")
diff -r 50f82a593a983ca8ee997a0650fc03eddfdbb30b -r 90f900be7a36433fdd48941cae4bc91066ff5c76 yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py
@@ -12,24 +12,33 @@
import numpy as np
from yt.testing import \
- assert_allclose_units, requires_file, requires_module
+ assert_allclose_units, requires_file, requires_module, \
+ assert_almost_equal, assert_array_almost_equal
from yt.analysis_modules.absorption_spectrum.absorption_line import \
voigt_old, voigt_scipy
from yt.analysis_modules.absorption_spectrum.api import AbsorptionSpectrum
from yt.analysis_modules.cosmological_observation.api import LightRay
+from yt.config import ytcfg
import tempfile
import os
import shutil
+from yt.utilities.on_demand_imports import \
+ _h5py as h5
+
+test_dir = ytcfg.get("yt", "test_data_dir")
COSMO_PLUS = "enzo_cosmology_plus/AMRCosmology.enzo"
-
+COSMO_PLUS_SINGLE = "enzo_cosmology_plus/RD0009/RD0009"
+HI_SPECTRUM_COSMO = "absorption_spectrum_data/enzo_lyman_alpha_cosmo_spec.h5"
+HI_SPECTRUM_COSMO_FILE = os.path.join(test_dir, HI_SPECTRUM_COSMO)
+HI_SPECTRUM = "absorption_spectrum_data/enzo_lyman_alpha_spec.h5"
+HI_SPECTRUM_FILE = os.path.join(test_dir, HI_SPECTRUM)
@requires_file(COSMO_PLUS)
-def test_absorption_spectrum():
+ at requires_file(HI_SPECTRUM_COSMO)
+def test_absorption_spectrum_cosmo():
"""
- This test is simply following the description in the docs for how to
- generate an absorption spectrum from a cosmological light ray for one
- of the sample datasets
+ This test generates an absorption spectrum from a cosmological light ray
"""
# Set up in a temp dir
@@ -37,7 +46,7 @@
curdir = os.getcwd()
os.chdir(tmpdir)
- lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.1)
+ lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.03)
lr.make_light_ray(seed=1234567,
fields=['temperature', 'density', 'H_number_density'],
@@ -65,22 +74,30 @@
sp.add_continuum(my_label, field, wavelength, normalization, index)
wavelength, flux = sp.make_spectrum('lightray.h5',
- output_file='spectrum.txt',
+ output_file='spectrum.h5',
line_list_file='lines.txt',
use_peculiar_velocity=True)
+ # load just-generated hdf5 file of spectral data (for consistency)
+ f_new = h5.File('spectrum.h5', 'r')
+
+ # load standard data for comparison
+ f_old = h5.File(HI_SPECTRUM_COSMO_FILE, 'r')
+
+ # compare between standard data and current data for each array saved
+ # (wavelength, flux, tau)
+ for key in f_old.keys():
+ assert_array_almost_equal(f_new[key].value, f_old[key].value, 10)
+
# clean up
os.chdir(curdir)
shutil.rmtree(tmpdir)
-
- at requires_file(COSMO_PLUS)
- at requires_module("astropy")
-def test_absorption_spectrum_fits():
+ at requires_file(COSMO_PLUS_SINGLE)
+ at requires_file(HI_SPECTRUM)
+def test_absorption_spectrum_non_cosmo():
"""
- This test is simply following the description in the docs for how to
- generate an absorption spectrum from a cosmological light ray for one
- of the sample datasets. Outputs to fits file if astropy is installed.
+ This test generates an absorption spectrum from a non-cosmological light ray
"""
# Set up in a temp dir
@@ -88,11 +105,114 @@
curdir = os.getcwd()
os.chdir(tmpdir)
- lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.1)
+ lr = LightRay(COSMO_PLUS_SINGLE)
- lr.make_light_ray(seed=1234567,
+ ray_start = [0,0,0]
+ ray_end = [1,1,1]
+ lr.make_light_ray(start_position=ray_start, end_position=ray_end,
fields=['temperature', 'density', 'H_number_density'],
- get_los_velocity=True,
+ data_filename='lightray.h5')
+
+ sp = AbsorptionSpectrum(1200.0, 1300.0, 10001)
+
+ my_label = 'HI Lya'
+ field = 'H_number_density'
+ wavelength = 1215.6700 # Angstromss
+ f_value = 4.164E-01
+ gamma = 6.265e+08
+ mass = 1.00794
+
+ sp.add_line(my_label, field, wavelength, f_value,
+ gamma, mass, label_threshold=1.e10)
+
+ wavelength, flux = sp.make_spectrum('lightray.h5',
+ output_file='spectrum.h5',
+ line_list_file='lines.txt',
+ use_peculiar_velocity=True)
+
+ # load just-generated hdf5 file of spectral data (for consistency)
+ f_new = h5.File('spectrum.h5', 'r')
+
+ # load standard data for comparison
+ f_old = h5.File(HI_SPECTRUM_FILE, 'r')
+
+ # compare between standard data and current data for each array saved
+ # (wavelength, flux, tau)
+ for key in f_old.keys():
+ assert_array_almost_equal(f_new[key].value, f_old[key].value, 10)
+
+ # clean up
+ os.chdir(curdir)
+ shutil.rmtree(tmpdir)
+
+ at requires_file(COSMO_PLUS_SINGLE)
+def test_equivalent_width_conserved():
+ """
+ This tests that the equivalent width of the optical depth is conserved
+ regardless of the bin width employed in wavelength space.
+ Unresolved lines should still deposit optical depth into the spectrum.
+ """
+
+ # Set up in a temp dir
+ tmpdir = tempfile.mkdtemp()
+ curdir = os.getcwd()
+ os.chdir(tmpdir)
+
+ lr = LightRay(COSMO_PLUS_SINGLE)
+
+ ray_start = [0,0,0]
+ ray_end = [1,1,1]
+ lr.make_light_ray(start_position=ray_start, end_position=ray_end,
+ fields=['temperature', 'density', 'H_number_density'],
+ data_filename='lightray.h5')
+
+ my_label = 'HI Lya'
+ field = 'H_number_density'
+ wave = 1215.6700 # Angstromss
+ f_value = 4.164E-01
+ gamma = 6.265e+08
+ mass = 1.00794
+
+ lambda_min= 1200
+ lambda_max= 1300
+ lambda_bin_widths = [1e-3, 1e-2, 1e-1, 1e0, 1e1]
+ total_tau = []
+
+ for lambda_bin_width in lambda_bin_widths:
+ n_lambda = ((lambda_max - lambda_min)/ lambda_bin_width) + 1
+ sp = AbsorptionSpectrum(lambda_min=lambda_min, lambda_max=lambda_max,
+ n_lambda=n_lambda)
+ sp.add_line(my_label, field, wave, f_value, gamma, mass)
+ wavelength, flux = sp.make_spectrum('lightray.h5')
+ total_tau.append((lambda_bin_width * sp.tau_field).sum())
+
+ # assure that the total tau values are all within 1e-5 of each other
+ for tau in total_tau:
+ assert_almost_equal(tau, total_tau[0], 5)
+
+ # clean up
+ os.chdir(curdir)
+ shutil.rmtree(tmpdir)
+
+
+ at requires_file(COSMO_PLUS_SINGLE)
+ at requires_module("astropy")
+def test_absorption_spectrum_fits():
+ """
+ This test generates an absorption spectrum and saves it as a fits file.
+ """
+
+ # Set up in a temp dir
+ tmpdir = tempfile.mkdtemp()
+ curdir = os.getcwd()
+ os.chdir(tmpdir)
+
+ lr = LightRay(COSMO_PLUS_SINGLE)
+
+ ray_start = [0,0,0]
+ ray_end = [1,1,1]
+ lr.make_light_ray(start_position=ray_start, end_position=ray_end,
+ fields=['temperature', 'density', 'H_number_density'],
data_filename='lightray.h5')
sp = AbsorptionSpectrum(900.0, 1800.0, 10000)
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