Multimodel Example

This notebook shows a concrete multi_transit setup where two atmospheric regions are described in separate TauREx parameter files and combined into a single weighted transmission spectrum. The example mirrors a east-west style configuration: a hotter east side terminator, and a more cloudy west side terminator.

Data Note

The example below uses the opacity files prepared in the setup notebook together with checked-in parameter files under examples/parfiles. To keep the example runnable in a lightweight local setup, the regional models use absorption-only contributions.

The notebook changes into the project root before loading the main file so that the repo-relative paths inside the .par files resolve consistently.

from pathlib import Path
import os
import sys
import numpy as np

for candidate in (Path.cwd(), *Path.cwd().parents):
    if (candidate / '_shared.py').exists():
        sys.path.insert(0, str(candidate))
        break
    nested = candidate / 'examples' / 'notebooks'
    if (nested / '_shared.py').exists():
        sys.path.insert(0, str(nested))
        break

from _shared import PROJECT_ROOT, ensure_opacity_data

ensure_opacity_data(download=False)
os.chdir(PROJECT_ROOT)

project_root = Path.cwd()
parfile_dir = project_root / 'examples' / 'parfiles'
day_path = parfile_dir / 'multimodel_east.par'
night_path = parfile_dir / 'multimodel_west.par'
main_path = parfile_dir / 'multimodel_main.par'

print(f'Project root: {project_root.name}')
print(f'Day-side file: {day_path.relative_to(project_root)}')
print(f'Night-side file: {night_path.relative_to(project_root)}')
print(f'Main multimodel file: {main_path.relative_to(project_root)}')

Project root: taurex3
Day-side file: examples/parfiles/multimodel_east.par
Night-side file: examples/parfiles/multimodel_west.par
Main multimodel file: examples/parfiles/multimodel_main.par

PROJECT_ROOT

PosixPath('/home/quent/Software/taurex3')

How the parameter files are split

This example uses three checked-in .par files to create a forward model with limb asymetry.

• examples/parfiles/multimodel_main.par is the entry point. It defines the shared global paths, star, planet, a default atmosphere block, and the multi_transit model itself.

• examples/parfiles/multimodel_east.par defines the hotter, east terminator region.

• examples/parfiles/multimodel_weast.par defines the cooler, west terminator region.

Only the main file is printed below, because that is the file a user would normally start from. The regional files are linked through the parfiles entry inside that main file. More background on TauREx input files is in 08_parameter_files.ipynb and in the input-file docs at docs/user/taurex/inputfile.rst.

print(f'Main parameter file: {main_path.relative_to(project_root)}\n')
print(main_path.read_text())

print('\nCompanion regional files:')
print(f'  - {day_path.relative_to(project_root)}')
print(f'  - {night_path.relative_to(project_root)}')

Main parameter file: examples/parfiles/multimodel_main.par

[Global]
xsec_path = examples/tmp/xsec
cia_path = examples/tmp/cia

[Chemistry]
chemistry_type = taurex
fill_gases = H2,He
ratio = 0.1756

    [[H2O]]
    gas_type = constant
    mix_ratio = 1e-3

    [[CO2]]
    gas_type = constant
    mix_ratio = 1e-2

[Temperature]
profile_type = isothermal
T = 1800

[Pressure]
profile_type = Simple
atm_min_pressure = 1e-5
atm_max_pressure = 1e5
nlayers = 80

[Planet]
planet_type = Simple
planet_mass = 0.74
planet_radius = 1.38

[Star]
star_type = blackbody
temperature = 6117
radius = 1.16

[Model]
model_type = multi_transit
parfiles = examples/parfiles/multimodel_east.par,examples/parfiles/multimodel_west.par
fractions = 0.65,

[Observation]
observed_spectrum = examples/parfiles/multimodel_observation.dat

[Optimizer]
optimizer = nestle
num_live_points = 300
method = multi

[Fitting]

fr_1:fit = True
fr_1:bounds = 0.0, 1.0

planet_radius:fit = True
planet_radius:bounds = 1.0, 1.6

m1_H2O:fit = True
m1_H2O:bounds = 1e-12, 1e-1
m2_H2O:fit = True
m2_H2O:bounds = 1e-12, 1e-1

m1_T:fit = True
m1_T:bounds = 1000, 2500
m2_T:fit = True
m2_T:bounds = 1000, 2500

m2_clouds_pressure:fit = True
m2_clouds_pressure:bounds = 1e3,1e-2


Companion regional files:
  - examples/parfiles/multimodel_east.par
  - examples/parfiles/multimodel_west.par

What TauREx is doing here

multimodel_main.par is the only file passed directly to ParameterParser. TauREx reads that main file first, sees model_type = multi_transit, and then follows the parfiles entry to load the east-side and west-side region files.

The main file supplies the shared context for the run: global data paths, star, planet, and a default atmospheric setup. The regional files then override the parts that differ between regions, here mainly the temperature and clouds for the west side of the terminator. The fractions entry in the main file sets how strongly each regional spectrum contributes to the final combined spectrum. For retrievals, provide only N-1 fraction values so TauREx can infer the last region automatically and keep the combined weights normalized.

from taurex.parameter import ParameterParser

parser = ParameterParser()
parser.read(str(main_path))
parser.setup_globals()
model = parser.generate_appropriate_model()
model.build()

internal_model = model._multimodel

print('Top-level model:', type(model).__name__)
print('Internal multimodel:', type(internal_model).__name__)
print('Number of regions:', len(internal_model._sub_models))
print('Fractions:', internal_model._fractions)
print('First fitting parameters:', list(model.fittingParameters)[:8])

Numba not installed, using numpy instead
Top-level model: MultiParameterTransitModel
Internal multimodel: MultiTransitModel
Number of regions: 2
Fractions: [0.65, 0.35]
First fitting parameters: ['fr_1', 'planet_mass', 'planet_radius', 'planet_distance', 'planet_sma', 'distance', 'm1_T', 'm2_T']

import matplotlib.pyplot as plt

combined_grid, combined_flux, _, _ = model.model()
combined_wavelength = 10000 / combined_grid

fig, ax = plt.subplots(figsize=(8, 4.5))
for index, sub_model in enumerate(internal_model._sub_models, start=1):
    region_grid, region_flux, _, _ = sub_model.model(wngrid=combined_grid, cutoff_grid=False)
    ax.plot(10000 / region_grid, region_flux, label=f'Region {index}', alpha=0.7)

ax.plot(combined_wavelength, combined_flux, color='black', linewidth=2.2, label='Weighted composite')
ax.set_xlabel('Wavelength [micron]')
ax.set_ylabel('Transit depth')
ax.set_title('Concrete multi_transit example')
ax.legend()
ax.grid(alpha=0.25)
plt.show()

Adding some binning to prepare a retrieval

This section is to prepare a future tutorial.

The full high-res spectrum is now binned down to the arbitrary \(\lambda = [1.2, 7.8]\) grid at \(\mathcal{R} = 20\) to get ready for a retrieval example with the parfiles.

The last cell saves the spectrum in a text file (it is already there, so saving is disabled)!

from taurex.binning import FluxBinner
from taurex.util.util import create_grid_res
np.random.seed(0)

grid = create_grid_res(20, 1.2, 7.8)
binner = FluxBinner(grid[:,0], grid[:,1])

observation_wl, observation_sp, _, _ = binner.bindown(combined_wavelength[::-1], combined_flux[::-1])
fake_error = 150 * 1e-6
observation_sp_rand = np.random.normal(observation_sp, fake_error) # add manually 150ppm error

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(observation_wl, observation_sp, color='black', label='Binned Spectrum')
ax.errorbar(observation_wl, observation_sp_rand, yerr=fake_error, fmt = '.', label='Simulated Observations')

ax.set_xscale('log')
ax.set_xlabel('Wavelength [micron]')
ax.set_ylabel('Transit depth')
ax.set_title('Simulated observation')
ax.set_xticks([2, 3, 4, 6], ['2', '3','4', '6'])
ax.legend()
ax.grid(alpha=0.25)
plt.show()

observation_path = parfile_dir / 'multimodel_observation.dat'

saved_spectrum = np.zeros((len(observation_wl), 3))
saved_spectrum[:,0] = observation_wl
saved_spectrum[:,1] = observation_sp_rand
saved_spectrum[:,2] = fake_error
#np.savetxt(observation_path, saved_spectrum) ## Disabled since file already there.