[yt-svn] commit/yt: xarthisius: Merged in brittonsmith/yt/yt-3.0 (pull request #1086)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Jul 29 13:13:24 PDT 2014
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/b2a39372ea82/
Changeset: b2a39372ea82
Branch: yt-3.0
User: xarthisius
Date: 2014-07-29 22:13:15
Summary: Merged in brittonsmith/yt/yt-3.0 (pull request #1086)
Updating Halo Mass Function
Affected #: 8 files
diff -r db184f28bbac06e15f81292ec306ed0b7f102a69 -r b2a39372ea82bfa7a22b78af7fddfe6137073a8c doc/source/analyzing/analysis_modules/halo_catalogs.rst
--- a/doc/source/analyzing/analysis_modules/halo_catalogs.rst
+++ b/doc/source/analyzing/analysis_modules/halo_catalogs.rst
@@ -1,3 +1,4 @@
+.. _halo_catalog:
Creating Halo Catalogs
======================
diff -r db184f28bbac06e15f81292ec306ed0b7f102a69 -r b2a39372ea82bfa7a22b78af7fddfe6137073a8c doc/source/analyzing/analysis_modules/halo_mass_function.rst
--- a/doc/source/analyzing/analysis_modules/halo_mass_function.rst
+++ b/doc/source/analyzing/analysis_modules/halo_mass_function.rst
@@ -2,29 +2,34 @@
Halo Mass Function
==================
-.. sectionauthor:: Stephen Skory <sskory at physics.ucsd.edu>
-.. versionadded:: 1.6
The Halo Mass Function extension is capable of outputting the halo mass function
-for a collection haloes (input), and/or an analytical fit over a given mass range
+for a collection halos (input), and/or an analytical fit over a given mass range
for a set of specified cosmological parameters.
-
This extension is based on code generously provided by Brian O'Shea.
General Overview
----------------
-In order to run this extension on a dataset, the haloes need to be located
-(using HOP, FOF or Parallel HOP, see :ref:`halo_finding`),
-and their virial masses determined using the
-HaloProfiler.
-Please see the step-by-step how-to which puts these steps together
-(:ref:`hmf_howto`).
-If an optional analytical fit is desired, the correct initial
-cosmological parameters will need to be input as well. These initial parameters
-are not stored in an Enzo dataset, so they must be set by hand.
-An analytical fit can be found without referencing a particular dataset or
-set of haloes, but all the cosmological parameters need to be set by hand.
+A halo mass function can be created for the halos identified in a cosmological
+simulation, as well as analytic fits using any arbitrary set of cosmological
+paramters. In order to create a mass function for simulated halos, they must
+first be identified (using HOP, FOF, or Rockstar, see
+:ref:`halo_catalog`) and loaded as a halo dataset object. The distribution of
+halo masses will then be found, and can be compared to the analytic prediction
+at the same redshift and using the same cosmological parameters as were used
+in the simulation. Care should be taken in this regard, as the analytic fit
+requires the specification of cosmological parameters that are not necessarily
+stored in the halo or simulation datasets, and must be specified by the user.
+Efforts have been made to set reasonable defaults for these parameters, but
+setting them to identically match those used in the simulation will produce a
+much better comparison.
+
+Analytic halo mass functions can also be created without a halo dataset by
+providing either a simulation dataset or specifying cosmological parameters by
+hand. yt includes 5 analytic fits for the halo mass function which can be
+selected.
+
Analytical Fits
---------------
@@ -45,115 +50,181 @@
The Tinker fit is for the :math:`\Delta=300` fits given in the paper, which
appears to fit HOP threshold=80.0 fairly well.
-Analyze Simulated Haloes
-------------------------
-If an analytical fit is not needed, it is simple to analyze a set of
-haloes. The ``halo_file`` needs to be specified, and
-``fitting_function`` does not need to be specified.
-``num_sigma_bins`` is how many bins the halo masses are sorted into.
-The default is 360. ``mass_column`` is the zero-indexed column of the
-``halo_file`` file that contains the halo masses. The default is 5, which
-corresponds to the sixth column of data in the file.
+Basic Halo Mass Function Creation
+---------------------------------
+
+The simplest way to create a halo mass function object is to simply pass it no
+arguments and let it use the default cosmological parameters.
+
+.. code-block:: python
+
+ from yt.analysis_modules.halo_mass_function.api import *
+
+ hmf = HaloMassFcn()
+
+This will create a HaloMassFcn object off of which arrays holding the information
+about the analytic mass function hang. Creating the halo mass function for a set
+of simulated halos requires only the loaded halo dataset to be passed as an
+argument. This also creates the analytic mass function using all parameters that
+can be extracted from the halo dataset, at the same redshift, spanning a similar
+range of halo masses.
.. code-block:: python
from yt.mods import *
from yt.analysis_modules.halo_mass_function.api import *
- ds = load("data0030")
- hmf = HaloMassFcn(ds, halo_file="FilteredQuantities.out", num_sigma_bins=200,
- mass_column=5)
-Attached to ``hmf`` is the convenience function ``write_out``, which saves
-the halo mass function to a text file. By default, both the halo analysis (``haloes``) and
-fit (``fit``) are written to (different) text files, but they can be turned on or off
-explicitly. ``prefix`` sets the name used for the file(s). The haloes file
-is named ``prefix-haloes.dat``, and the fit file ``prefix-fit.dat``.
-Continued from above, invoking this command:
+ my_halos = load("rockstar_halos/halos_0.0.bin")
+ hmf = HaloMassFcn(halos_ds=my_halos)
-.. code-block:: python
-
- hmf.write_out(prefix='hmf', fit=False, haloes=True)
-
-will save the haloes data to a file named ``hmf-haloes.dat``. The contents
-of the ``-haloes.dat`` file is three columns:
-
- 1. log10 of mass (Msolar, NOT Msolar/h) for this bin.
- 2. mass (Msolar/h) for this bin.
- 3. cumulative number density of halos (per Mpc^3, NOT h^3/Mpc^3) in this bin.
-
-Analytical Halo Mass Function Fit
----------------------------------
-
-When an analytical fit is desired, in nearly all cases several cosmological
-parameters will need to be specified by hand. These parameters are not
-stored with Enzo datasets. In the case where both the haloes and an analytical
-fit are desired, the analysis is instantiated as below.
-``sigma8input``, ``primordial_index`` and ``omega_baryon0`` should be set to
-the same values as
-``PowerSpectrumSigma8``, ``PowerSpectrumPrimordialIndex`` and
-``CosmologyOmegaBaryonNow`` from the
-`inits <http://lca.ucsd.edu/projects/enzo/wiki/UserGuide/RunningInits>`_
-parameter file used to set up the simulation.
-``fitting_function`` is set to values 1 through 4 from the list of available
-fits above.
+A simulation dataset can be passed along with additonal cosmological parameters
+to create an analytic mass function.
.. code-block:: python
from yt.mods import *
from yt.analysis_modules.halo_mass_function.api import *
- ds = load("data0030")
- hmf = HaloMassFcn(ds, halo_file="FilteredQuantities.out",
- sigma8input=0.9, primordial_index=1., omega_baryon0=0.06,
- fitting_function=4)
- hmf.write_out(prefix='hmf')
-Both the ``-haloes.dat`` and ``-fit.dat`` files are written to disk.
-The contents of the ``-fit.dat`` file is four columns:
+ my_ds = load("RD0027/RedshiftOutput0027")
+ hmf = HaloMassFcn(simulation_ds=my_ds, omega_baryon0=0.05, primordial_index=0.96,
+ sigma8 = 0.8, log_mass_min=5, log_mass_max=9)
- 1. log10 of mass (Msolar, NOT Msolar/h) for this bin.
- 2. mass (Msolar/h) for this bin.
- 3. (dn/dM)*dM (differential number density of halos, per Mpc^3 (NOT h^3/Mpc^3) in this bin.
- 4. cumulative number density of halos (per Mpc^3, NOT h^3/Mpc^3) in this bin.
-
-Below is an example of the output for both the haloes and the (Warren)
-analytical fit, for three datasets. The black lines are the calculated
-halo mass functions, and the blue lines the analytical fit set by initial
-conditions. This simulation shows typical behavior, in that there are too
-few small haloes compared to the fit due to lack of mass and gravity resolution
-for small haloes. But at higher mass ranges, the simulated haloes are quite close
-to the analytical fit.
-
-.. image:: _images/halo_mass_function.png
- :width: 350
- :height: 400
-
-The analytical fit can be found without referencing a particular dataset. In this
-case, all the various cosmological parameters need to be specified by hand.
-``omega_matter0`` is the fraction of universe that is made up of matter
-(baryons and dark matter). ``omega_lambda0`` is the fractional proportion due
-to dark energy. In a flat universe, ``omega_matter0`` + ``omega_lambda0`` = 1.
-``this_redshift`` is the redshift for which you wish to generate a fit.
-``log_mass_min`` and ``log_mass_max`` are the logarithmic ends of the mass range for which
-you wish to calculate the fit.
+The analytic mass function can be created for a set of arbitrary cosmological
+parameters without any dataset being passed as an argument.
.. code-block:: python
from yt.mods import *
from yt.analysis_modules.halo_mass_function.api import *
- hmf = HaloMassFcn(None, omega_matter0=0.3, omega_lambda0=0.7,
- omega_baryon0=0.06, hubble0=.7, this_redshift=0., log_mass_min=8.,
- log_mass_max=13., sigma8input=0.9, primordial_index=1.,
- fitting_function=1)
- hmf.write_out(prefix="hmf-press-schechter", fit=True, haloes=False)
-It is possible to access the output of the halo mass function without saving
-to disk. The content is stored in arrays hanging off the ``HaloMassFcn``
-object:
+ hmf = HaloMassFcn(omega_baryon0=0.05, omega_matter0=0.27,
+ omega_lambda0=0.73, hubble0=0.7, this_redshift=10,
+ log_mass_min=5, log_mass_max=9, fitting_function=5)
- * ``hmf.logmassarray`` for log10 of mass bin.
- * ``hmf.massarray`` for mass bin.
- * ``hmf.dn_M_z`` for (dn/dM)*dM (analytical fit).
- * ``hmf.nofmz_cum`` for cumulative number density of halos (analytical fit).
- * ``hmf.dis`` for cumulative number density of halos (from provided halo
- halo information).
+Keyword Arguments
+-----------------
+
+ * **simulation_ds** (*Simulation dataset object*)
+ The loaded simulation dataset, used to set cosmological paramters.
+ Default : None.
+
+ * **halos_ds** (*Halo dataset object*)
+ The halos from a simulation to be used for creation of the
+ halo mass function in the simulation.
+ Default : None.
+
+ * **make_analytic** (*bool*)
+ Whether or not to calculate the analytic mass function to go with
+ the simulated halo mass function. Automatically set to true if a
+ simulation dataset is provided.
+ Default : True.
+
+ * **omega_matter0** (*float*)
+ The fraction of the universe made up of matter (dark and baryonic).
+ Default : 0.2726.
+
+ * **omega_lambda0** (*float*)
+ The fraction of the universe made up of dark energy.
+ Default : 0.7274.
+
+ * **omega_baryon0** (*float*)
+ The fraction of the universe made up of baryonic matter. This is not
+ always stored in the datset and should be checked by hand.
+ Default : 0.0456.
+
+ * **hubble0** (*float*)
+ The expansion rate of the universe in units of 100 km/s/Mpc.
+ Default : 0.704.
+
+ * **sigma8** (*float*)
+ The amplitude of the linear power spectrum at z=0 as specified by
+ the rms amplitude of mass-fluctuations in a top-hat sphere of radius
+ 8 Mpc/h. This is not always stored in the datset and should be
+ checked by hand.
+ Default : 0.86.
+
+ * **primoridal_index** (*float*)
+ This is the index of the mass power spectrum before modification by
+ the transfer function. A value of 1 corresponds to the scale-free
+ primordial spectrum. This is not always stored in the datset and
+ should be checked by hand.
+ Default : 1.0.
+
+ * **this_redshift** (*float*)
+ The current redshift.
+ Default : 0.
+
+ * **log_mass_min** (*float*)
+ The log10 of the mass of the minimum of the halo mass range. This is
+ set automatically by the range of halo masses if a simulated halo
+ dataset is provided. If a halo dataset if not provided and no value
+ is specified, it will be set to 5. Units: M_solar
+ Default : None.
+
+ * **log_mass_max** (*float*)
+ The log10 of the mass of the maximum of the halo mass range. This is
+ set automatically by the range of halo masses if a simulated halo
+ dataset is provided. If a halo dataset if not provided and no value
+ is specified, it will be set to 16. Units: M_solar
+ Default : None.
+
+ * **num_sigma_bins** (*float*)
+ The number of bins (points) to use for the calculation of the
+ analytic mass function.
+ Default : 360.
+
+ * **fitting_function** (*int*)
+ Which fitting function to use. 1 = Press-Schechter, 2 = Jenkins,
+ 3 = Sheth-Tormen, 4 = Warren, 5 = Tinker
+ Default : 4.
+
+
+Outputs
+-------
+
+A HaloMassFnc object has several arrays hanging off of it containing the
+ * **masses_sim**: Halo masses from simulated halos. Units: M_solar
+
+ * **n_cumulative_sim**: Number density of halos with mass greater than the
+ corresponding mass in masses_sim. Units: comoving Mpc^-3
+
+ * **masses_analytic**: Masses used for the generation of the analytic mass
+ function. Units: M_solar
+
+ * **n_cumulative_analytic**: Number density of halos with mass greater then
+ the corresponding mass in masses_analytic. Units: comoving Mpc^-3
+
+ * **dndM_dM_analytic**: Differential number density of halos, (dn/dM)*dM.
+
+After the mass function has been created for both simulated halos and the
+corresponding analytic fits, they can be plotted though something along the
+lines of
+
+.. code-block:: python
+
+ from yt.mods import *
+ from yt.analysis_modules.halo_mass_function.api import *
+ import matplotlib.pyplot as plt
+
+ my_halos = load("rockstar_halos/halos_0.0.bin")
+ hmf = HaloMassFcn(halos_ds=my_halos)
+
+ plt.loglog(hmf.masses_sim, hmf.n_cumulative_sim)
+ plt.loglog(hmf.masses_analytic, hmf.n_cumulative_analytic)
+
+Attached to ``hmf`` is the convenience function ``write_out``, which saves the
+halo mass function to a text file. (continued from above)
+.. code-block:: python
+
+ hmf.write_out(prefix='hmf', analytic=True, simulated=True)
+
+This writes the files `hmf-analytic.dat' with columns
+ * mass [Msun]
+ * cumulative number density of halos [comoving Mpc^-3]
+ * (dn/dM)*dM (differential number density of halos) [comoving Mpc^-3
+
+and the file `hmf-simulated.dat' with columns
+
+ * mass [Msun]
+ * cumulative number density of halos [comoving Mpc^-3]
diff -r db184f28bbac06e15f81292ec306ed0b7f102a69 -r b2a39372ea82bfa7a22b78af7fddfe6137073a8c doc/source/analyzing/analysis_modules/hmf_howto.rst
--- a/doc/source/analyzing/analysis_modules/hmf_howto.rst
+++ /dev/null
@@ -1,143 +0,0 @@
-.. _hmf_howto:
-
-Halo Mass Function: Start to Finish
-===================================
-
-This how-to steps through the three simple steps it takes to find the halo
-mass function for an ``enzo`` dataset. The first step is to find the haloes,
-second to find the virial mass (gas + dark matter) contained in each halo using
-the halo profiler, and
-finally build a halo mass function of the haloes and an analytical fit
-for comparison.
-
-Halo Finding
-------------
-
-The first step is find the haloes in the simulation. There are a number of ways
-to do this explained in detail in :ref:`halo_finding`.
-You are encouraged to read about the differences between the halo finders and
-experiment with different settings and functions.
-For the purposes of
-the halo mass function, Friends of Friends (FOF) is probably not the best choice.
-Therefore, below
-is a simple example of how to run HOP on a dataset. This example will also
-write out a text file with the halo particulars, which will be used in the next
-step.
-
-.. code-block:: python
-
- from yt.mods import *
- ds = load("data0001")
- halo_list = HaloFinder(ds)
- halo_list.write_out("HopAnalysis.out")
-
-The only important columns of data in the text file ``HopAnalysis.out``
-are the halo number, center of
-mass, and maximum radius. These are the only values that the next step requires.
-
-Halo Profiling
---------------
-
-The halo profiler is a powerful tool that can analyze
-haloes in many ways. It is beneficial to read its documentation to become
-familiar with it before using it.
-For this exercise, only the virial mass of each
-halo is important. This script below will take the output from the previous step
-and find the virial mass for each halo, and save it to a text file.
-
-.. code-block:: python
-
- import yt.analysis_modules.halo_profiler.api as HP
- hp = HP.HaloProfiler("data0001", halo_list_file='HopAnalysis.out')
- hp.add_halo_filter(HP.VirialFilter,must_be_virialized=True,
- overdensity_field='ActualOverdensity',
- virial_overdensity=200,
- virial_filters=[['TotalMassMsun','>=','1e8']],
- virial_quantities=['TotalMassMsun','RadiusMpc'])
- hp.make_profiles(filename="VirialHaloes.out")
-
-This script limits the output to virialized haloes with mass greater than or
-equal to 1e8 solar masses. If you run into problems, try pre-filtering problem
-haloes.
-
-Halo Mass Function
-------------------
-
-The halo mass function extension (see :ref:`halo_mass_function`) reads in the
-contents of ``VirialHaloes.out`` and can output files which can be used
-to make a plot. In this script below, the results from the previous step are
-read in, analyzed, and at the same time an analytical fit is calculated, and
-finally the results are written out to two files. The plot found from the haloes
-is saved to a file named ``hmf-haloes.dat`` and the fit to ``hmf-fit.dat``.
-For these files, the columns that are most likely to be needed for plotting are
-the first (halo mass bin), and the fourth and third (cumulative number density
-of haloes) for the ``-fit.dat`` and ``-haloes.dat`` files. See
-the full halo mass function documentation for more details on the contents of
-the files.
-
-.. code-block:: python
-
- from yt.mods import *
- from yt.analysis_modules.halo_mass_function.api import *
- ds = load("data0001")
- hmf = HaloMassFcn(ds, halo_file="VirialHaloes.out",
- sigma8input=0.9, primordial_index=1., omega_baryon0=0.06,
- fitting_function=4, mass_column=5, num_sigma_bins=200)
- hmf.write_out(prefix='hmf')
-
-Inside the call to ``HaloMassFcn`` there are several hand-set parameters that
-*must* be correct for the analytical fit to be correct. The three cosmological
-parameters (``sigma8input``, ``primordial_index`` and ``omega_baryon0``) are
-not stored with ``Enzo`` datasets, so they must be found from the initial
-conditions of the simulation. ``mass_column`` is set to 5 because that is the
-zero-indexed ordinal of the mass column in the ``VirialHaloes.out`` file.
-``num_sigma_bins`` controls how many mass bins the haloes are dropped into,
-and ``fitting_function`` controls which analytical function to use.
-
-Putting it All Together
------------------------
-
-It is not necessary to run each step separately from the others. This example
-below will run all steps at once.
-
-.. code-block:: python
-
- from yt.mods import *
- import yt.analysis_modules.halo_profiler.api as HP
- from yt.analysis_modules.halo_mass_function.api import *
-
- # If desired, start loop here.
- ds = load("data0001")
-
- halo_list = HaloFinder(ds)
- halo_list.write_out("HopAnalysis.out")
-
- hp = HP.HaloProfiler("data0001", halo_list_file='HopAnalysis.out')
- hp.add_halo_filter(HP.VirialFilter,must_be_virialized=True,
- overdensity_field='ActualOverdensity',
- virial_overdensity=200,
- virial_filters=[['TotalMassMsun','>=','1e8']],
- virial_quantities=['TotalMassMsun','RadiusMpc'])
- hp.make_profiles(filename="VirialHaloes.out")
-
- hmf = HaloMassFcn(ds, halo_file="VirialHaloes.out",
- sigma8input=0.9, primordial_index=1., omega_baryon0=0.06,
- fitting_function=4, mass_column=5, num_sigma_bins=200)
- hmf.write_out(prefix='hmf')
- # End loop here.
-
-The script above will work in parallel which can reduce runtimes substantially.
-If this analysis is to be run on a sequence of datasets, the section that needs
-to be inside the loop is shown bracketed by comments. Be careful how the
-output files are named as to not over-write output from previous loop cycles.
-
-Plotting
---------
-
-When plotting the output, be careful about the units of the output for the
-halo mass function. The figure shown in the documentation (on this page:
-:ref:`halo_mass_function`) has the number density of haloes per (h^-1 Mpc)^3,
-which is different than the output of the halo mass extension (which is
-haloes per (Mpc)^3). To get the same units as the figure for the ``-fit.dat``
-and ``-haloes.dat`` files, divide the fourth and third column by the comoving
-volume cubed of the simulation, respectively when plotting.
diff -r db184f28bbac06e15f81292ec306ed0b7f102a69 -r b2a39372ea82bfa7a22b78af7fddfe6137073a8c yt/analysis_modules/halo_analysis/halo_finding_methods.py
--- a/yt/analysis_modules/halo_analysis/halo_finding_methods.py
+++ b/yt/analysis_modules/halo_analysis/halo_finding_methods.py
@@ -21,6 +21,7 @@
HaloCatalogDataset
from yt.frontends.stream.data_structures import \
load_particles
+from yt.units.dimensions import length
from yt.utilities.operator_registry import \
OperatorRegistry
@@ -136,5 +137,13 @@
attr_val = getattr(data_ds, attr)
setattr(particle_ds, attr, attr_val)
particle_ds.current_time = particle_ds.current_time.in_cgs()
+
+ particle_ds.unit_registry.modify("h", particle_ds.hubble_constant)
+ # Comoving lengths
+ for my_unit in ["m", "pc", "AU", "au"]:
+ new_unit = "%scm" % my_unit
+ particle_ds.unit_registry.add(new_unit, particle_ds.unit_registry.lut[my_unit][0] /
+ (1 + particle_ds.current_redshift),
+ length, "\\rm{%s}/(1+z)" % my_unit)
return particle_ds
diff -r db184f28bbac06e15f81292ec306ed0b7f102a69 -r b2a39372ea82bfa7a22b78af7fddfe6137073a8c yt/analysis_modules/halo_mass_function/halo_mass_function.py
--- a/yt/analysis_modules/halo_mass_function/halo_mass_function.py
+++ b/yt/analysis_modules/halo_mass_function/halo_mass_function.py
@@ -1,5 +1,5 @@
"""
-halo_mass_function - Halo Mass Function and supporting functions.
+Halo Mass Function and supporting functions.
@@ -16,198 +16,314 @@
import numpy as np
import math, time
-from yt.funcs import *
-from yt.utilities.parallel_tools.parallel_analysis_interface import \
- ParallelDummy, \
- ParallelAnalysisInterface, \
- parallel_blocking_call
-from yt.utilities.physical_constants import \
- cm_per_mpc, \
- mass_sun_cgs
+from yt.funcs import mylog
+from yt.units.yt_array import \
+ YTArray, \
+ YTQuantity
from yt.utilities.physical_ratios import \
rho_crit_g_cm3_h2
-class HaloMassFcn(ParallelAnalysisInterface):
+class HaloMassFcn():
+ r"""
+ Initalize a HaloMassFcn object to analyze the distribution of halos as
+ a function of mass. A mass function can be created for a set of
+ simulated halos, an analytic fit to can be created for a redshift and
+ set of cosmological parameters, or both can be created.
+
+ Provided with a halo dataset object, this will make a the mass function
+ for simulated halos. Prodiving a simulation dataset will set as many
+ of the cosmological parameters as possible for the creation of the
+ analytic mass function.
+
+ The HaloMassFcn object has arrays hanging off of it containing the mass
+ function information.
+
+ masses_sim : Array
+ Halo masses from simulated halos. Units: M_solar.
+ n_cumulative_sim : Array
+ Number density of halos with mass greater than the corresponding
+ mass in masses_sim (simulated). Units: comoving Mpc^-3
+ masses_analytic : Array
+ Masses used for the generation of the analytic mass function, Units:
+ M_solar.
+ n_cumulative_analytic : Array
+ Number density of halos with mass greater then the corresponding
+ mass in masses_analytic (analytic). Units: comoving Mpc^-3
+ dndM_dM_analytic : Array
+ Differential number density of halos, (dn/dM)*dM (analytic).
+
+ The HaloMassFcn object also has a convenience function write_out() that
+ will write out the data to disk.
+
+ Creating a HaloMassFcn object with no arguments will produce an analytic
+ mass function at redshift = 0 using default cosmolocigal values.
+
+ Parameters
+ ----------
+ simulation_ds : Simulation dataset object
+ The loaded simulation dataset, used to set cosmological paramters.
+ Default : None.
+ halos_ds : Halo dataset object
+ The halos from a simulation to be used for creation of the
+ halo mass function in the simulation.
+ Default : None.
+ make_analytic : bool
+ Whether or not to calculate the analytic mass function to go with
+ the simulated halo mass function. Automatically set to true if a
+ simulation dataset is provided.
+ Default : True.
+ omega_matter0 : float
+ The fraction of the universe made up of matter (dark and baryonic).
+ Default : 0.2726.
+ omega_lambda0 : float
+ The fraction of the universe made up of dark energy.
+ Default : 0.7274.
+ omega_baryon0 : float
+ The fraction of the universe made up of baryonic matter. This is not
+ always stored in the datset and should be checked by hand.
+ Default : 0.0456.
+ hubble0 : float
+ The expansion rate of the universe in units of 100 km/s/Mpc.
+ Default : 0.704.
+ sigma8 : float
+ The amplitude of the linear power spectrum at z=0 as specified by
+ the rms amplitude of mass-fluctuations in a top-hat sphere of radius
+ 8 Mpc/h. This is not always stored in the datset and should be
+ checked by hand.
+ Default : 0.86.
+ primoridal_index : float
+ This is the index of the mass power spectrum before modification by
+ the transfer function. A value of 1 corresponds to the scale-free
+ primordial spectrum. This is not always stored in the datset and
+ should be checked by hand.
+ Default : 1.0.
+ this_redshift : float
+ The current redshift.
+ Default : 0.
+ log_mass_min : float
+ The log10 of the mass of the minimum of the halo mass range. This is
+ set automatically by the range of halo masses if a simulated halo
+ dataset is provided. If a halo dataset if not provided and no value
+ is specified, it will be set to 5. Units: M_solar
+ Default : None.
+ log_mass_max : float
+ The log10 of the mass of the maximum of the halo mass range. This is
+ set automatically by the range of halo masses if a simulated halo
+ dataset is provided. If a halo dataset if not provided and no value
+ is specified, it will be set to 16. Units: M_solar
+ Default : None.
+ num_sigma_bins : float
+ The number of bins (points) to use for the calculation of the
+ analytic mass function.
+ Default : 360.
+ fitting_function : int
+ Which fitting function to use. 1 = Press-Schechter, 2 = Jenkins,
+ 3 = Sheth-Tormen, 4 = Warren, 5 = Tinker
+ Default : 4.
+
+ Examples
+ --------
+
+ This creates the halo mass function for a halo dataset from a simulation
+ and the analytic mass function at the same redshift as the dataset,
+ using as many cosmological parameters as can be pulled from the dataset.
+
+ >>> halos_ds = load("rockstar_halos/halo_0.0.bin")
+ >>> hmf = HaloMassFcn(halos_ds=halos_ds)
+ >>> plt.loglog(hmf.masses_sim, hmf.n_cumulative_sim)
+ >>> plt.loglog(hmf.masses_analytic, hmf.n_cumulative_analytic)
+ >>> plt.savefig("mass_function.png")
+
+ This creates only the analytic halo mass function for a simulation
+ dataset, with default values for cosmological paramters not stored in
+ the dataset.
+
+ >>> ds = load("enzo_tiny_cosmology/DD0046/DD0046")
+ >>> hmf = HaloMassFcn(simulation_ds=ds)
+ >>> plt.loglog(hmf.masses_analytic, hmf.n_cumulative_analytic)
+ >>> plt.savefig("mass_function.png")
+
+ This creates the analytic mass function for an arbitrary set of
+ cosmological parameters, with neither a simulation nor halo dataset.
+
+ >>> hmf = HaloMassFcn(omega_baryon0=0.05, omega_matter0=0.27,
+ omega_lambda0=0.73, hubble0=0.7, this_redshift=10,
+ log_mass_min=5, log_mass_max=9)
+ >>> plt.loglog(hmf.masses_analytic, hmf.n_cumulative_analytic)
+ >>> plt.savefig("mass_function.png")
"""
- Initalize a HaloMassFcn object to analyze the distribution of haloes
- as a function of mass.
- :param halo_file (str): The filename of the output of the Halo Profiler.
- Default=None.
- :param omega_matter0 (float): The fraction of the universe made up of
- matter (dark and baryonic). Default=None.
- :param omega_lambda0 (float): The fraction of the universe made up of
- dark energy. Default=None.
- :param omega_baryon0 (float): The fraction of the universe made up of
- ordinary baryonic matter. This should match the value
- used to create the initial conditions, using 'inits'. This is
- *not* stored in the enzo datset so it must be checked by hand.
- Default=0.05.
- :param hubble0 (float): The expansion rate of the universe in units of
- 100 km/s/Mpc. Default=None.
- :param sigma8input (float): The amplitude of the linear power
- spectrum at z=0 as specified by the rms amplitude of mass-fluctuations
- in a top-hat sphere of radius 8 Mpc/h. This should match the value
- used to create the initial conditions, using 'inits'. This is
- *not* stored in the enzo datset so it must be checked by hand.
- Default=0.86.
- :param primoridal_index (float): This is the index of the mass power
- spectrum before modification by the transfer function. A value of 1
- corresponds to the scale-free primordial spectrum. This should match
- the value used to make the initial conditions using 'inits'. This is
- *not* stored in the enzo datset so it must be checked by hand.
- Default=1.0.
- :param this_redshift (float): The current redshift. Default=None.
- :param log_mass_min (float): The log10 of the mass of the minimum of the
- halo mass range. Default=None.
- :param log_mass_max (float): The log10 of the mass of the maximum of the
- halo mass range. Default=None.
- :param num_sigma_bins (float): The number of bins (points) to use for
- the calculations and generated fit. Default=360.
- :param fitting_function (int): Which fitting function to use.
- 1 = Press-schechter, 2 = Jenkins, 3 = Sheth-Tormen, 4 = Warren fit
- 5 = Tinker
- Default=4.
- :param mass_column (int): The column of halo_file that contains the
- masses of the haloes. Default=4.
- """
- def __init__(self, ds, halo_file=None, omega_matter0=None, omega_lambda0=None,
- omega_baryon0=0.05, hubble0=None, sigma8input=0.86, primordial_index=1.0,
- this_redshift=None, log_mass_min=None, log_mass_max=None, num_sigma_bins=360,
- fitting_function=4, mass_column=5):
- ParallelAnalysisInterface.__init__(self)
- self.ds = ds
- self.halo_file = halo_file
+ def __init__(self, simulation_ds=None, halos_ds=None, make_analytic=True,
+ omega_matter0=0.2726, omega_lambda0=0.7274, omega_baryon0=0.0456, hubble0=0.704,
+ sigma8=0.86, primordial_index=1.0, this_redshift=0, log_mass_min=None,
+ log_mass_max=None, num_sigma_bins=360, fitting_function=4):
+ self.simulation_ds = simulation_ds
+ self.halos_ds = halos_ds
self.omega_matter0 = omega_matter0
self.omega_lambda0 = omega_lambda0
self.omega_baryon0 = omega_baryon0
self.hubble0 = hubble0
- self.sigma8input = sigma8input
+ self.sigma8 = sigma8
self.primordial_index = primordial_index
self.this_redshift = this_redshift
self.log_mass_min = log_mass_min
self.log_mass_max = log_mass_max
self.num_sigma_bins = num_sigma_bins
self.fitting_function = fitting_function
- self.mass_column = mass_column
-
- # Determine the run mode.
- if halo_file is None:
- # We are hand-picking our various cosmological parameters
- self.mode = 'single'
- else:
- # Make the fit using the same cosmological parameters as the dataset.
- self.mode = 'haloes'
- self.omega_matter0 = self.ds.omega_matter
- self.omega_lambda0 = self.ds.omega_lambda
- self.hubble0 = self.ds.hubble_constant
- self.this_redshift = self.ds.current_redshift
- self.read_haloes()
- if self.log_mass_min == None:
- self.log_mass_min = math.log10(min(self.haloes))
- if self.log_mass_max == None:
- self.log_mass_max = math.log10(max(self.haloes))
+ self.make_analytic = make_analytic
+ self.make_simulated = False
+ """
+ If we want to make an analytic mass function, grab what we can from either the
+ halo file or the data set, and make sure that the user supplied everything else
+ that is needed.
+ """
+ # If we don't have any datasets, make the analytic function with user values
+ if simulation_ds is None and halos_ds is None:
+ # Set a reasonable mass min and max if none were provided
+ if log_mass_min is None:
+ self.log_mass_min = 5
+ if log_mass_max is None:
+ self.log_mass_max = 16
+ # If we're making the analytic function...
+ if self.make_analytic == True:
+ # Try to set cosmological parameters from the simulation dataset
+ if simulation_ds is not None:
+ self.omega_matter0 = self.simulation_ds.omega_matter
+ self.omega_lambda0 = self.simulation_ds.omega_lambda
+ self.hubble0 = self.simulation_ds.hubble_constant
+ self.this_redshift = self.simulation_ds.current_redshift
+ # Set a reasonable mass min and max if none were provided
+ if log_mass_min is None:
+ self.log_mass_min = 5
+ if log_mass_max is None:
+ self.log_mass_max = 16
+ # If we have a halo dataset but not a simulation dataset, use that instead
+ if simulation_ds is None and halos_ds is not None:
+ self.omega_matter0 = self.halos_ds.omega_matter
+ self.omega_lambda0 = self.halos_ds.omega_lambda
+ self.hubble0 = self.halos_ds.hubble_constant
+ self.this_redshift = self.halos_ds.current_redshift
+ # If the user didn't specify mass min and max, set them from the halos
+ if log_mass_min is None:
+ self.set_mass_from_halos("min_mass")
+ if log_mass_max is None:
+ self.set_mass_from_halos("max_mass")
+ # Do the calculations.
+ self.sigmaM()
+ self.dndm()
+ # Return the mass array in M_solar rather than M_solar/h
+ self.masses_analytic = YTArray(self.masses_analytic/self.hubble0, "Msun")
+ # The halo arrays will already have yt units, but the analytic forms do
+ # not. If a dataset has been provided, use that to give them units. At the
+ # same time, convert to comoving (Mpc)^-3
+ if simulation_ds is not None:
+ self.n_cumulative_analytic = simulation_ds.arr(self.n_cumulative_analytic,
+ "(Mpccm)**(-3)")
+ elif halos_ds is not None:
+ self.n_cumulative_analytic = halos_ds.arr(self.n_cumulative_analytic,
+ "(Mpccm)**(-3)")
+ else:
+ from yt.units.unit_registry import UnitRegistry
+ from yt.units.dimensions import length
+ hmf_registry = UnitRegistry()
+ for my_unit in ["m", "pc", "AU", "au"]:
+ new_unit = "%scm" % my_unit
+ hmf_registry.add(new_unit,
+ hmf_registry.lut[my_unit][0] /
+ (1 + self.this_redshift),
+ length, "\\rm{%s}/(1+z)" % my_unit)
+ self.n_cumulative_analytic = YTArray(self.n_cumulative_analytic,
+ "(Mpccm)**(-3)",
+ registry=hmf_registry)
- # Input error check.
- if self.mode == 'single':
- if omega_matter0 == None or omega_lambda0 == None or \
- hubble0 == None or this_redshift == None or log_mass_min == None or\
- log_mass_max == None:
- mylog.error("All of these parameters need to be set:")
- mylog.error("[omega_matter0, omega_lambda0, \
- hubble0, this_redshift, log_mass_min, log_mass_max]")
- mylog.error("[%s,%s,%s,%s,%s,%s]" % (omega_matter0,\
- omega_lambda0, hubble0, this_redshift,\
- log_mass_min, log_mass_max))
- return None
-
- # Poke the user to make sure they're doing it right.
- mylog.info(
+
"""
- Please make sure these are the correct values! They are
- not stored in enzo datasets, so must be entered by hand.
- sigma8input=%f primordial_index=%f omega_baryon0=%f
- """ % (self.sigma8input, self.primordial_index, self.omega_baryon0))
-
- # Do the calculations.
- self.sigmaM()
- self.dndm()
-
- if self.mode == 'haloes':
- self.bin_haloes()
+ If a halo file has been supplied, make a mass function for the simulated halos.
+ """
+ if halos_ds is not None:
+ # Used to check if a simulated halo mass funciton exists to write out
+ self.make_simulated=True
+ # Calculate the simulated halo mass function
+ self.create_sim_hmf()
- def write_out(self, prefix='HMF', fit=True, haloes=True):
+ """
+ If we're making an analytic fit and have a halo dataset, but don't have log_mass_min
+ or log_mass_max from the user, set it from the range of halo masses.
+ """
+ def set_mass_from_halos(self, which_limit):
+ data_source = self.halos_ds.all_data()
+ if which_limit is "min_mass":
+ self.log_mass_min = \
+ int(np.log10(np.amin(data_source["particle_mass"].in_units("Msun"))))
+ if which_limit is "max_mass":
+ self.log_mass_max = \
+ int(np.log10(np.amax(data_source["particle_mass"].in_units("Msun"))))+1
+
+ """
+ Here's where we create the halo mass functions from simulated halos
+ """
+ def create_sim_hmf(self):
+ data_source = self.halos_ds.all_data()
+ # We're going to use indices to count the number of halos above a given mass
+ masses_sim = np.sort(data_source["particle_mass"].in_units("Msun"))
+ # Determine the size of the simulation volume in comoving Mpc**3
+ sim_volume = self.halos_ds.domain_width.in_units('Mpccm').prod()
+ n_cumulative_sim = np.arange(len(masses_sim),0,-1)
+ # We don't want repeated halo masses, and the unique indices tell us which values
+ # correspond to distinct halo masses.
+ self.masses_sim, unique_indices = np.unique(masses_sim, return_index=True)
+ # Now make this an actual number density of halos as a function of mass.
+ self.n_cumulative_sim = n_cumulative_sim[unique_indices]/sim_volume
+ # masses_sim and n_cumulative_sim are now set, but remember that the log10 quantities
+ # are what is usually plotted for a halo mass function.
+
+ def write_out(self, prefix='HMF', analytic=True, simulated=True):
"""
Writes out the halo mass functions to file(s) with prefix *prefix*.
"""
- # First the fit file.
- if fit:
- fitname = prefix + '-fit.dat'
- fp = self.comm.write_on_root(fitname)
- line = \
- """#Columns:
-#1. log10 of mass (Msolar, NOT Msolar/h)
-#2. mass (Msolar/h)
-#3. (dn/dM)*dM (differential number density of haloes, per Mpc^3 (NOT h^3/Mpc^3)
-#4. cumulative number density of halos (per Mpc^3, NOT h^3/Mpc^3)
-"""
- fp.write(line)
- for i in xrange(self.logmassarray.size - 1):
- line = "%e\t%e\t%e\t%e\n" % (self.logmassarray[i], self.massarray[i],
- self.dn_M_z[i], self.nofmz_cum[i])
+ # First the analytic file, check that analytic fit exists and was requested
+ if analytic:
+ if self.make_analytic:
+ fitname = prefix + '-analytic.dat'
+ fp = open(fitname, "w")
+ line = \
+ "#Columns:\n" + \
+ "#1. mass (M_solar)\n" + \
+ "#2. cumulative number density of halos [comoving Mpc^-3]\n" + \
+ "#3. (dn/dM)*dM (differential number density of halos) [comoving Mpc^-3]\n"
fp.write(line)
- fp.close()
- if self.mode == 'haloes' and haloes:
- haloname = prefix + '-haloes.dat'
- fp = self.comm.write_on_root(haloname)
- line = \
- """#Columns:
-#1. log10 of mass (Msolar, NOT Msolar/h)
-#2. mass (Msolar/h)
-#3. cumulative number density of haloes (per Mpc^3, NOT h^3/Mpc^3)
-"""
- fp.write(line)
- for i in xrange(self.logmassarray.size - 1):
- line = "%e\t%e\t%e\n" % (self.logmassarray[i], self.massarray[i],
- self.dis[i])
+ for i in xrange(self.masses_analytic.size - 1):
+ line = "%e\t%e\t%e\n" % (self.masses_analytic[i],
+ self.n_cumulative_analytic[i],
+ self.dndM_dM_analytic[i])
+ fp.write(line)
+ fp.close()
+ # If the analytic halo mass function wasn't created, warn the user
+ else:
+ mylog.warning("The analytic halo mass function was not created and cannot be " +
+ "written out! Specify its creation with " +
+ "HaloMassFcn(make_analytic=True, other_args) when creating the " +
+ "HaloMassFcn object.")
+ # Write out the simulated mass fucntion if it exists and was requested
+ if simulated:
+ if self.make_simulated:
+ haloname = prefix + '-simulated.dat'
+ fp = open(haloname, "w")
+ line = \
+ "#Columns:\n" + \
+ "#1. mass [Msun]\n" + \
+ "#2. cumulative number density of halos [comoving Mpc^-3]\n"
fp.write(line)
- fp.close()
-
- def read_haloes(self):
- """
- Read in the virial masses of the haloes.
- """
- mylog.info("Reading halo masses from %s" % self.halo_file)
- f = open(self.halo_file,'r')
- line = f.readline()
- if line == "":
- self.haloes = np.array([])
- return
- while line[0] == '#':
- line = f.readline()
- self.haloes = []
- while line:
- line = line.split()
- mass = float(line[self.mass_column])
- if mass > 0:
- self.haloes.append(float(line[self.mass_column]))
- line = f.readline()
- f.close()
- self.haloes = np.array(self.haloes)
-
- def bin_haloes(self):
- """
- With the list of virial masses, find the halo mass function.
- """
- bins = np.logspace(self.log_mass_min,
- self.log_mass_max,self.num_sigma_bins)
- avgs = (bins[1:]+bins[:-1])/2.
- dis, bins = np.histogram(self.haloes,bins)
- # add right to left
- for i,b in enumerate(dis):
- dis[self.num_sigma_bins-i-3] += dis[self.num_sigma_bins-i-2]
- if i == (self.num_sigma_bins - 3): break
-
- self.dis = dis / (self.ds.domain_width * self.ds.units["mpccm"]).prod()
+ for i in xrange(self.masses_sim.size - 1):
+ line = "%e\t%e\n" % (self.masses_sim[i],
+ self.n_cumulative_sim[i])
+ fp.write(line)
+ fp.close()
+ # If the simulated halo mass function wasn't created, warn the user
+ else:
+ mylog.warning("The simulated halo mass function was not created and cannot " +
+ "be written out! Specify its creation by providing a loaded " +
+ "halo dataset with HaloMassFcn(ds_halos=loaded_halo_dataset, " +
+ "other_args) when creating the HaloMassFcn object.")
def sigmaM(self):
"""
@@ -223,10 +339,8 @@
Outputs: four columns of data containing the following information:
- 1) log mass (Msolar)
- 2) mass (Msolar/h)
- 3) Radius (comoving Mpc/h)
- 4) sigma (normalized) using Msun/h as the input
+ 1) mass (Msolar/h)
+ 2) sigma (normalized) using Msun/h as the input
The arrays output are used later.
"""
@@ -239,24 +353,21 @@
mylog.error("You should probably fix your cosmology parameters!")
# output arrays
- # 1) log10 of mass (Msolar, NOT Msolar/h)
- self.Rarray = np.empty(self.num_sigma_bins,dtype='float64')
- # 2) mass (Msolar/h)
- self.logmassarray = np.empty(self.num_sigma_bins, dtype='float64')
- # 3) spatial scale corresponding to that radius (Mpc/h)
- self.massarray = np.empty(self.num_sigma_bins, dtype='float64')
- # 4) sigma(M, z=0, where mass is in Msun/h)
+ # 1) mass (M_solar/h), changed to M_solar/h at output
+ self.masses_analytic = np.empty(self.num_sigma_bins, dtype='float64')
+ # 2) sigma(M, z=0, where mass is in Msun/h)
self.sigmaarray = np.empty(self.num_sigma_bins, dtype='float64')
# get sigma_8 normalization
R = 8.0; # in units of Mpc/h (comoving)
sigma8_unnorm = math.sqrt(self.sigma_squared_of_R(R));
- sigma_normalization = self.sigma8input / sigma8_unnorm;
+ sigma_normalization = self.sigma8 / sigma8_unnorm;
# rho0 in units of h^2 Msolar/Mpc^3
- rho0 = self.omega_matter0 * \
- rho_crit_g_cm3_h2 * cm_per_mpc**3 / mass_sun_cgs
+ rho0 = YTQuantity(self.omega_matter0 * rho_crit_g_cm3_h2 * self.hubble0**2,
+ 'g/cm**3').in_units('Msun/Mpc**3')
+ rho0 = rho0.value.item()
# spacing in mass of our sigma calculation
dm = (float(self.log_mass_max) - self.log_mass_min)/self.num_sigma_bins;
@@ -280,32 +391,31 @@
R = thisradius; # h^-1 Mpc (comoving)
- self.Rarray[i] = thisradius; # h^-1 Mpc (comoving)
- self.logmassarray[i] = thislogmass; # Msun (NOT Msun/h)
- self.massarray[i] = thismass; # Msun/h
+ self.masses_analytic[i] = thismass; # Msun/h
# get normalized sigma(R)
self.sigmaarray[i] = math.sqrt(self.sigma_squared_of_R(R)) * sigma_normalization;
# All done!
def dndm(self):
-
# constants - set these before calling any functions!
# rho0 in units of h^2 Msolar/Mpc^3
- rho0 = self.omega_matter0 * \
- rho_crit_g_cm3_h2 * cm_per_mpc**3 / mass_sun_cgs
+ rho0 = YTQuantity(self.omega_matter0 * rho_crit_g_cm3_h2 * self.hubble0**2,
+ 'g/cm**3').in_units('Msun/Mpc**3')
+ rho0 = rho0.value.item()
+
self.delta_c0 = 1.69; # critical density for turnaround (Press-Schechter)
- nofmz_cum = 0.0; # keep track of cumulative number density
+ n_cumulative_analytic = 0.0; # keep track of cumulative number density
# Loop over masses, going BACKWARD, and calculate dn/dm as well as the
# cumulative mass function.
# output arrays
# 5) (dn/dM)*dM (differential number density of halos, per Mpc^3 (NOT h^3/Mpc^3)
- self.dn_M_z = np.empty(self.num_sigma_bins, dtype='float64')
+ self.dndM_dM_analytic = np.empty(self.num_sigma_bins, dtype='float64')
# 6) cumulative number density of halos (per Mpc^3, NOT h^3/Mpc^3)
- self.nofmz_cum = np.zeros(self.num_sigma_bins, dtype='float64')
+ self.n_cumulative_analytic = np.zeros(self.num_sigma_bins, dtype='float64')
for j in xrange(self.num_sigma_bins - 1):
i = (self.num_sigma_bins - 2) - j
@@ -313,25 +423,25 @@
thissigma = self.sigmaof_M_z(i, self.this_redshift);
nextsigma = self.sigmaof_M_z(i+1, self.this_redshift);
- # calc dsigmadm - has units of h (since massarray has units of h^-1)
- dsigmadm = (nextsigma-thissigma) / (self.massarray[i+1] - self.massarray[i]);
+ # calc dsigmadm - has units of h (since masses_analytic has units of h^-1)
+ dsigmadm = (nextsigma-thissigma) / (self.masses_analytic[i+1] - self.masses_analytic[i]);
# calculate dn(M,z) (dn/dM * dM)
# this has units of h^3 since rho0 has units of h^2, dsigmadm
- # has units of h, and massarray has units of h^-1
- dn_M_z = -1.0 / thissigma * dsigmadm * rho0 / self.massarray[i] * \
- self.multiplicityfunction(thissigma)*(self.massarray[i+1] - self.massarray[i]);
+ # has units of h, and masses_analytic has units of h^-1
+ dndM_dM_analytic = -1.0 / thissigma * dsigmadm * rho0 / self.masses_analytic[i] * \
+ self.multiplicityfunction(thissigma)*(self.masses_analytic[i+1] - self.masses_analytic[i]);
# scale by h^3 to get rid of all factors of h
- dn_M_z *= math.pow(self.hubble0, 3.0);
+ dndM_dM_analytic *= math.pow(self.hubble0, 3.0);
# keep track of cumulative number density
- if dn_M_z > 1.0e-20:
- nofmz_cum += dn_M_z;
+ if dndM_dM_analytic > 1.0e-20:
+ n_cumulative_analytic += dndM_dM_analytic;
# Store this.
- self.nofmz_cum[i] = nofmz_cum
- self.dn_M_z[i] = dn_M_z
+ self.n_cumulative_analytic[i] = n_cumulative_analytic
+ self.dndM_dM_analytic[i] = dndM_dM_analytic
def sigma_squared_of_R(self, R):
diff -r db184f28bbac06e15f81292ec306ed0b7f102a69 -r b2a39372ea82bfa7a22b78af7fddfe6137073a8c yt/frontends/enzo/tests/test_outputs.py
--- a/yt/frontends/enzo/tests/test_outputs.py
+++ b/yt/frontends/enzo/tests/test_outputs.py
@@ -18,7 +18,9 @@
requires_ds, \
small_patch_amr, \
big_patch_amr, \
- data_dir_load
+ data_dir_load, \
+ AnalyticHaloMassFunctionTest, \
+ SimulatedHaloMassFunctionTest
from yt.frontends.enzo.api import EnzoDataset
_fields = ("temperature", "density", "velocity_magnitude",
@@ -69,6 +71,19 @@
test_galaxy0030.__name__ = test.description
yield test
+enzotiny = "enzo_tiny_cosmology/DD0046/DD0046"
+ at requires_ds(enzotiny)
+def test_simulated_halo_mass_function():
+ ds = data_dir_load(enzotiny)
+ for finder in ["fof", "hop"]:
+ yield SimulatedHaloMassFunctionTest(ds, finder)
+
+ at requires_ds(enzotiny)
+def test_analytic_halo_mass_function():
+ ds = data_dir_load(enzotiny)
+ for fit in range(1, 6):
+ yield AnalyticHaloMassFunctionTest(ds, fit)
+
ecp = "enzo_cosmology_plus/DD0046/DD0046"
@requires_ds(ecp, big_data=True)
def test_ecp():
diff -r db184f28bbac06e15f81292ec306ed0b7f102a69 -r b2a39372ea82bfa7a22b78af7fddfe6137073a8c yt/utilities/answer_testing/framework.py
--- a/yt/utilities/answer_testing/framework.py
+++ b/yt/utilities/answer_testing/framework.py
@@ -576,6 +576,55 @@
for newc, oldc in zip(new_result["children"], old_result["children"]):
assert(newp == oldp)
+class SimulatedHaloMassFunctionTest(AnswerTestingTest):
+ _type_name = "SimulatedHaloMassFunction"
+ _attrs = ("finder",)
+
+ def __init__(self, ds_fn, finder):
+ super(SimulatedHaloMassFunctionTest, self).__init__(ds_fn)
+ self.finder = finder
+
+ def run(self):
+ from yt.analysis_modules.halo_analysis.api import HaloCatalog
+ from yt.analysis_modules.halo_mass_function.api import HaloMassFcn
+ hc = HaloCatalog(data_ds=self.ds, finder_method=self.finder)
+ hc.create()
+
+ hmf = HaloMassFcn(halos_ds=hc.halos_ds)
+ result = np.empty((2, hmf.masses_sim.size))
+ result[0] = hmf.masses_sim.d
+ result[1] = hmf.n_cumulative_sim.d
+ return result
+
+ def compare(self, new_result, old_result):
+ err_msg = ("Simulated halo mass functions not equation for " +
+ "%s halo finder.") % self.finder
+ assert_equal(new_result, old_result,
+ err_msg=err_msg, verbose=True)
+
+class AnalyticHaloMassFunctionTest(AnswerTestingTest):
+ _type_name = "AnalyticHaloMassFunction"
+ _attrs = ("fitting_function",)
+
+ def __init__(self, ds_fn, fitting_function):
+ super(AnalyticHaloMassFunctionTest, self).__init__(ds_fn)
+ self.fitting_function = fitting_function
+
+ def run(self):
+ from yt.analysis_modules.halo_mass_function.api import HaloMassFcn
+ hmf = HaloMassFcn(simulation_ds=self.ds,
+ fitting_function=self.fitting_function)
+ result = np.empty((2, hmf.masses_analytic.size))
+ result[0] = hmf.masses_analytic.d
+ result[1] = hmf.n_cumulative_analytic.d
+ return result
+
+ def compare(self, new_result, old_result):
+ err_msg = ("Analytic halo mass functions not equation for " +
+ "fitting function %d.") % self.fitting_function
+ assert_equal(new_result, old_result,
+ err_msg=err_msg, verbose=True)
+
def compare_image_lists(new_result, old_result, decimals):
fns = ['old.png', 'new.png']
num_images = len(old_result)
@@ -677,7 +726,7 @@
return comp_imgs
def compare(self, new_result, old_result):
compare_image_lists(new_result, old_result, self.decimals)
-
+
def requires_ds(ds_fn, big_data = False, file_check = False):
def ffalse(func):
diff -r db184f28bbac06e15f81292ec306ed0b7f102a69 -r b2a39372ea82bfa7a22b78af7fddfe6137073a8c yt/utilities/physical_ratios.py
--- a/yt/utilities/physical_ratios.py
+++ b/yt/utilities/physical_ratios.py
@@ -97,6 +97,8 @@
# flux
jansky_cgs = 1.0e-23
# Cosmological constants
+# Calculated with H = 100 km/s/Mpc, value given in units of h^2 g cm^-3
+# Multiply by h^2 to get the critical density in units of g cm^-3
rho_crit_g_cm3_h2 = 1.8788e-29
primordial_H_mass_fraction = 0.76
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