Skip to content

Commit

Permalink
Implement EmissionAbsorptionMismatchedModel (#10)
Browse files Browse the repository at this point in the history
* rework matched/mismatched models
  • Loading branch information
tvwenger authored Oct 15, 2024
1 parent cbda584 commit abcaab1
Show file tree
Hide file tree
Showing 15 changed files with 4,125 additions and 28 deletions.
7 changes: 3 additions & 4 deletions .readthedocs.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,9 @@ build:
os: ubuntu-22.04
tools:
python: "3.12"
# You can also specify other tool versions:
# nodejs: "19"
# rust: "1.64"
# golang: "1.19"
jobs:
pre_build:
- sphinx-apidoc -f --templatedir=docs/source/_templates/ -o docs/source/ caribou_hi/

# Build documentation in the "docs/" directory with Sphinx
sphinx:
Expand Down
17 changes: 13 additions & 4 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,8 @@ Read below to get started, and check out the tutorials and guides here: https://
- [`EmissionModel`](#emissionmodel)
- [`AbsorptionModel`](#absorptionmodel)
- [`EmissionAbsorptionModel`](#emissionabsorptionmodel)
- [`EmissionAbsorptionFFModel`](#emissionabsorptionffmodel)
- [`EmissionAbsorptionMatchedModel`](#emissionabsorptionmatchedmodel)
- [`EmissionAbsorptionMismatchedModel`](#emissionabsorptionmismatchedmodel)
- [`ordered`](#ordered)
- [Syntax \& Examples](#syntax--examples)
- [Issues and Contributing](#issues-and-contributing)
Expand Down Expand Up @@ -125,13 +126,21 @@ The models provided by `caribou_hi` are implemented in the [`bayes_spec`](https:
| `rms_emission` | Emission spectrum rms noise | `K` | ${\rm rms}_{T} \sim {\rm HalfNormal}(\sigma=p)$ | `1.0` |
| `rms_absorption` | Optical depth spectrum rms noise | `K` | ${\rm rms}_{\tau} \sim {\rm HalfNormal}(\sigma=p)$ | `0.01` |

## `EmissionAbsorptionFFModel`
## `EmissionAbsorptionMatchedModel`

Finally, `EmissionAbsorptionFFModel` is like `EmissionAbsorptionModel`, except it allows for the possibility of beam dilution in the emission spectrum. That is, the expected brightness temperature contribution from a cloud is modified by a parameter, `filling_factor`, which takes values between zero and one.
`EmissionAbsorptionMatchedModel` is like `EmissionAbsorptionModel`, except it allows for the possibility of beam dilution in the emission spectrum. That is, the expected brightness temperature contribution from a cloud is modified by a parameter, `filling_factor`, which takes values between zero and one. Use this model when the emission and absorption data have matching beam sizes.

| Cloud Parameter<br>`variable` | Parameter | Units | Prior, where<br>($p_0, p_1, \dots$) = `prior_{variable}` | Default<br>`prior_{variable}` |
| :---------------------------- | :------------- | :---- | :------------------------------------------------------- | :---------------------------- |
| `filling_factor` | Filling Factor | `` | $f \sim {\rm Uniform}(0, 1)$ | `` |
| `filling_factor` | Filling Factor | `` | $f \sim {\rm Uniform}(0, 1)$ | `` |

## `EmissionAbsorptionMismatchedModel`

`EmissionAbsorptionMismatchedModel` is like `EmissionAbsorptionMatchedModel`, except it also allows for the possibility that a cloud detected in emission is not seen in absorption, for example, due to a mis-match between the emission and absorption beam size. The absorption optical depth of each cloud is modified by a parameter, `absorption_weight`, which takes values between zero and one.

| Cloud Parameter<br>`variable` | Parameter | Units | Prior, where<br>($p_0, p_1, \dots$) = `prior_{variable}` | Default<br>`prior_{variable}` |
| :---------------------------- | :---------------- | :---- | :------------------------------------------------------- | :---------------------------- |
| `absorption_weight` | Absorption weight | `` | $w_\tau \sim {\rm Uniform}(0, 1)$ | `` |

## `ordered`

Expand Down
6 changes: 4 additions & 2 deletions caribou_hi/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,15 @@
"EmissionModel",
"AbsorptionModel",
"EmissionAbsorptionModel",
"EmissionAbsorptionFFModel",
"EmissionAbsorptionMatchedModel",
"EmissionAbsorptionMismatchedModel",
]

from caribou_hi.emission_model import EmissionModel
from caribou_hi.absorption_model import AbsorptionModel
from caribou_hi.emission_absorption_model import EmissionAbsorptionModel
from caribou_hi.emission_absorption_ff_model import EmissionAbsorptionFFModel
from caribou_hi.emission_absorption_matched_model import EmissionAbsorptionMatchedModel
from caribou_hi.emission_absorption_mismatched_model import EmissionAbsorptionMismatchedModel

from . import _version

Expand Down
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
"""
emission_absorption_ff_model.py
EmissionAbsorptionFFModel definition
emission_absorption_matched_model.py
EmissionAbsorptionMatchedModel definition
Copyright(C) 2024 by
Trey V. Wenger; [email protected]
Expand Down Expand Up @@ -28,11 +28,11 @@
from caribou_hi import physics


class EmissionAbsorptionFFModel(HIModel):
"""Definition of the EmissionAbsorptionFFModel model. SpecData keys must be "emission" and "absorption"."""
class EmissionAbsorptionMatchedModel(HIModel):
"""Definition of the EmissionAbsorptionMatchedModel model. SpecData keys must be "emission" and "absorption"."""

def __init__(self, *args, bg_temp: float = 3.77, **kwargs):
"""Initialize a new EmissionAbsorptionFFModel instance
"""Initialize a new EmissionAbsorptionMatchedModel instance
Parameters
----------
Expand All @@ -48,7 +48,7 @@ def __init__(self, *args, bg_temp: float = 3.77, **kwargs):
# Define TeX representation of each parameter
self.var_name_map.update(
{
"filling_factor": r"f",
"filling_factor": r"$f$",
}
)

Expand Down
111 changes: 111 additions & 0 deletions caribou_hi/emission_absorption_mismatched_model.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
"""
emission_absorption_mismatched_model.py
EmissionAbsorptionMismatchedModel definition
Copyright(C) 2024 by
Trey V. Wenger; [email protected]
GNU General Public License v3 (GNU GPLv3)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published
by the Free Software Foundation, either version 3 of the License,
or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""

import pymc as pm
import pytensor.tensor as pt

from caribou_hi.hi_model import HIModel
from caribou_hi import physics


class EmissionAbsorptionMismatchedModel(HIModel):
"""Definition of the EmissionAbsorptionMismatchedModel model. SpecData keys must be "emission" and "absorption"."""

def __init__(self, *args, bg_temp: float = 3.77, **kwargs):
"""Initialize a new EmissionAbsorptionMismatchedModel instance
Parameters
----------
bg_temp : float, optional
Assumed background temperature (K), by default 3.77
"""
# Initialize HIModel
super().__init__(*args, **kwargs)

# Save inputs
self.bg_temp = bg_temp

# Define TeX representation of each parameter
self.var_name_map.update(
{
"filling_factor": r"$f$",
"absorption_weight": r"$w_\tau$",
}
)

def add_priors(self, *args, **kwargs):
"""Add priors and deterministics to the model"""
super().add_priors(*args, **kwargs)

with self.model:
# Filling factor
_ = pm.Uniform("filling_factor", lower=0.0, upper=1.0, dims="cloud")

# Absorption weight
_ = pm.Uniform("absorption_weight", lower=0.0, upper=1.0, dims="cloud")

def add_likelihood(self):
"""Add likelihood to the model. SpecData key must be "emission"."""
# Predict optical depth spectrum (shape: spectral, clouds)
absorption_optical_depth = self.model["absorption_weight"] * physics.calc_optical_depth(
self.data["absorption"].spectral,
self.model["velocity"],
10.0 ** self.model["log10_NHI"],
self.model["tspin"],
self.model["fwhm"],
)
emission_optical_depth = physics.calc_optical_depth(
self.data["emission"].spectral,
self.model["velocity"],
10.0 ** self.model["log10_NHI"],
self.model["tspin"],
self.model["fwhm"],
)

# Sum over clouds
predicted_absorption = 1.0 - pt.exp(-absorption_optical_depth.sum(axis=1))

# Evaluate radiative transfer
predicted_emission = physics.radiative_transfer(
emission_optical_depth, self.model["tspin"], self.model["filling_factor"], self.bg_temp
)

# Add baseline models
baseline_models = self.predict_baseline()
predicted_absorption = predicted_absorption + baseline_models["absorption"]
predicted_emission = predicted_emission + baseline_models["emission"]

with self.model:
# Evaluate likelihood
_ = pm.Normal(
"absorption",
mu=predicted_absorption,
sigma=self.data["absorption"].noise,
observed=self.data["absorption"].brightness,
)
_ = pm.Normal(
"emission",
mu=predicted_emission,
sigma=self.data["emission"].noise,
observed=self.data["emission"].brightness,
)
9 changes: 9 additions & 0 deletions docs/source/_templates/module.rst.jinja
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
{%- if show_headings %}
{{- [basename] | join(' ') | e | heading }}

{% endif -%}
.. automodule:: {{ qualname }}
{%- for option in automodule_options %}
:{{ option }}:
{%- endfor %}

57 changes: 57 additions & 0 deletions docs/source/_templates/package.rst.jinja
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
{%- macro automodule(modname, options) -%}
.. automodule:: {{ modname }}
{%- for option in options %}
:{{ option }}:
{%- endfor %}
{%- endmacro %}

{%- macro toctree(docnames) -%}
.. toctree::
:maxdepth: {{ maxdepth }}
{% for docname in docnames %}
{{ docname }}
{%- endfor %}
{%- endmacro %}

{%- if is_namespace %}
{{- [pkgname] | join(" ") | e | heading }}
{% else %}
{{- [pkgname] | join(" ") | e | heading }}
{% endif %}

{%- if is_namespace %}
.. py:module:: {{ pkgname }}
{% endif %}

{%- if modulefirst and not is_namespace %}
{{ automodule(pkgname, automodule_options) }}
{% endif %}

{%- if subpackages %}
Subpackages
-----------

{{ toctree(subpackages) }}
{% endif %}

{%- if submodules %}
Submodules
----------
{% if separatemodules %}
{{ toctree(submodules) }}
{% else %}
{%- for submodule in submodules %}
{% if show_headings %}
{{- [submodule] | join(" ") | e | heading(2) }}
{% endif %}
{{ automodule(submodule, automodule_options) }}
{% endfor %}
{%- endif %}
{%- endif %}

{%- if not modulefirst and not is_namespace %}
Module contents
---------------

{{ automodule(pkgname, automodule_options) }}
{% endif %}
8 changes: 8 additions & 0 deletions docs/source/_templates/toc.rst.jinja
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
{{ header | heading }}

.. toctree::
:maxdepth: {{ maxdepth }}
{% for docname in docnames %}
{{ docname }}
{%- endfor %}

3 changes: 2 additions & 1 deletion docs/source/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ Installation
notebooks/emission_model
notebooks/absorption_model
notebooks/emission_absorption_model
notebooks/emission_absorption_ff_model
notebooks/emission_absorption_matched_model
notebooks/emission_absorption_mismatched_model
notebooks/optimization

.. toctree::
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,11 @@
"id": "2a6f87e3-0d3c-4da4-90c6-eb9e00bd1e29",
"metadata": {},
"source": [
"# `EmissionAbsorptionFFModel` Tutorial\n",
"# `EmissionAbsorptionMatchedModel` Tutorial\n",
"\n",
"Trey V. Wenger (c) October 2024\n",
"\n",
"Here we demonstrate the basic features of the `EmissionAbsorptionFFModel` model. `EmissionAbsorptionFFModel` models both 21cm emission and absorption observations allowing for a \"filling factor\" term to account for beam dilution effects in the emission spectra."
"Here we demonstrate the basic features of the `EmissionAbsorptionMatchedModel` model. `EmissionAbsorptionMatchedModel` models both 21cm emission and absorption observations allowing for a \"filling factor\" term to account for beam dilution effects in the emission spectra."
]
},
{
Expand Down Expand Up @@ -154,12 +154,12 @@
},
"outputs": [],
"source": [
"from caribou_hi import EmissionAbsorptionFFModel\n",
"from caribou_hi import EmissionAbsorptionMatchedModel\n",
"\n",
"# Initialize and define the model\n",
"n_clouds = 3\n",
"baseline_degree = 0\n",
"model = EmissionAbsorptionFFModel(\n",
"model = EmissionAbsorptionMatchedModel(\n",
" dummy_data,\n",
" n_clouds=n_clouds,\n",
" baseline_degree=baseline_degree,\n",
Expand Down Expand Up @@ -290,7 +290,7 @@
"outputs": [],
"source": [
"# Initialize and define the model\n",
"model = EmissionAbsorptionFFModel(\n",
"model = EmissionAbsorptionMatchedModel(\n",
" data,\n",
" n_clouds=n_clouds,\n",
" baseline_degree=baseline_degree,\n",
Expand Down Expand Up @@ -814,7 +814,7 @@
],
"source": [
"# Plot model graph\n",
"model.graph().render(\"emission_absorption_ff_model\", format=\"png\")\n",
"model.graph().render(\"emission_absorption_matched_model\", format=\"png\")\n",
"model.graph()"
]
},
Expand Down Expand Up @@ -3874,7 +3874,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"display_name": "bayes_spec-dev",
"language": "python",
"name": "python3"
},
Expand Down
Loading

0 comments on commit abcaab1

Please sign in to comment.