Multiple species¶

This example demonstrates how to combine multiple species in one go for the simulation.

First let's start with how you could combine data from multiple species together for a composite spectrum.

To do that, you could simply concatenate the databases for two or more species.

To ensure we retain the information which line corresponds with which species, we add a species column with this information.

For simplicity of this interactive example, the spectra are simulated based on the linestrength only, inherently assuming an equal amount for each selected species.

Furthermore, it assumes equilibration between rotational or vibrational states, such that a single $T_{rot}$ and $T_{vib}$ can be used to describe the spectrum.

In [4]:

from pathlib import Path
import Moose
import lmfit
import pandas as pd
import numpy as np
import panel as pn
import panel.widgets as pnw
from importlib import resources
import plotly.graph_objects as go

pn.extension('plotly')

species = Moose.database_files
params = lmfit.create_params(**Moose.default_params)
params['A'].vary = False 
params['b'].vary = False
params['wl_pad'].value = 200
params['T_vib'].value= 4000

def load_species_db(name, *args):
    df = Moose.query_DB(name,*args)
    df['Species'] = name
    return df.sort_values('air_wavelength')

dbs = pd.concat([load_species_db(n,(300,500)) for n in ['N2CB', 'C2_swan']])
x = np.linspace(300,500,2000)
y = Moose.model_for_fit(x,**{p:params[p].value for p in params}, sim_db=dbs)

from pathlib import Path import Moose import lmfit import pandas as pd import numpy as np import panel as pn import panel.widgets as pnw from importlib import resources import plotly.graph_objects as go pn.extension('plotly') species = Moose.database_files params = lmfit.create_params(**Moose.default_params) params['A'].vary = False params['b'].vary = False params['wl_pad'].value = 200 params['T_vib'].value= 4000 def load_species_db(name, *args): df = Moose.query_DB(name,*args) df['Species'] = name return df.sort_values('air_wavelength') dbs = pd.concat([load_species_db(n,(300,500)) for n in ['N2CB', 'C2_swan']]) x = np.linspace(300,500,2000) y = Moose.model_for_fit(x,**{p:params[p].value for p in params}, sim_db=dbs)

In [5]:

fig = go.Figure()
fig.add_scatter(x=x, y=y, name='Spectrum')
fig.update_xaxes(title_text='λ (nm)')
fig.update_yaxes(title_text='I (a.u.)')

select_db = pnw.MultiChoice(value=['N2CB','C2_swan',],options = species, name='Select transitions', max_items=4)
select_range = pnw.RangeSlider(name='λ (nm)', start=200, end=700, value=(300,500))
pop_fracs = [pnw.FloatInput(name=f'Pop. species {i}', value=1,start=0,end=100, step=0.1) for i in range(1,5)]
status_msg = pnw.TextInput(name='Log', disabled=True)

input_widgets = {p:pnw.FloatSlider(name=p, value=params[p].value, start=params[p].min, end=params[p].max, step=params[p].value/10) for p in params if params[p].vary is True}

def cb_update_spectrum(event):
    x = np.linspace(*select_range.value,2000)
    if event.new != []:
        dbs = pd.concat([load_species_db(n, select_range.value) for n in select_db.value])
        y = np.zeros(x.shape[0])
        for i, (n,g) in enumerate(dbs.groupby('Species')):
            y+=Moose.model_for_fit(x,**{k:v.value for k,v in input_widgets.items()}, **{p:params[p].value for p in params if p not in input_widgets.keys()}, sim_db=g)*pop_fracs[i].value
        with fig.batch_update():
            fig.data[0].x = x
            fig.data[0].y=y
    else:
        fig.data[0].x = x
        fig.data[0].y = np.zeros(x.shape[0])
        status_msg.value = "Please select at least one species"
    

for wid in list(input_widgets.values())+[select_range,select_db]+pop_fracs:
    if hasattr(wid,'value_throttled'):
        to_watch = 'value_throttled'
    else:
        to_watch = 'value'
    wid.param.watch(cb_update_spectrum,to_watch)


plot_widget = pn.pane.Plotly(fig, sizing_mode='stretch_width', height=350)
layout = pn.Column(plot_widget,pn.FlexBox(select_db, select_range, status_msg), pn.FlexBox(*input_widgets.values()), pn.Row(*pop_fracs))
display(layout)

fig = go.Figure() fig.add_scatter(x=x, y=y, name='Spectrum') fig.update_xaxes(title_text='λ (nm)') fig.update_yaxes(title_text='I (a.u.)') select_db = pnw.MultiChoice(value=['N2CB','C2_swan',],options = species, name='Select transitions', max_items=4) select_range = pnw.RangeSlider(name='λ (nm)', start=200, end=700, value=(300,500)) pop_fracs = [pnw.FloatInput(name=f'Pop. species {i}', value=1,start=0,end=100, step=0.1) for i in range(1,5)] status_msg = pnw.TextInput(name='Log', disabled=True) input_widgets = {p:pnw.FloatSlider(name=p, value=params[p].value, start=params[p].min, end=params[p].max, step=params[p].value/10) for p in params if params[p].vary is True} def cb_update_spectrum(event): x = np.linspace(*select_range.value,2000) if event.new != []: dbs = pd.concat([load_species_db(n, select_range.value) for n in select_db.value]) y = np.zeros(x.shape[0]) for i, (n,g) in enumerate(dbs.groupby('Species')): y+=Moose.model_for_fit(x,**{k:v.value for k,v in input_widgets.items()}, **{p:params[p].value for p in params if p not in input_widgets.keys()}, sim_db=g)*pop_fracs[i].value with fig.batch_update(): fig.data[0].x = x fig.data[0].y=y else: fig.data[0].x = x fig.data[0].y = np.zeros(x.shape[0]) status_msg.value = "Please select at least one species" for wid in list(input_widgets.values())+[select_range,select_db]+pop_fracs: if hasattr(wid,'value_throttled'): to_watch = 'value_throttled' else: to_watch = 'value' wid.param.watch(cb_update_spectrum,to_watch) plot_widget = pn.pane.Plotly(fig, sizing_mode='stretch_width', height=350) layout = pn.Column(plot_widget,pn.FlexBox(select_db, select_range, status_msg), pn.FlexBox(*input_widgets.values()), pn.Row(*pop_fracs)) display(layout)

Model and fit the combined spectrum of two species¶

Alternatively, the provided sample file contains emission of the CN violet system and the N$_2$ second positive system.

In the case of the former there are multiple overlapping vibrational bands, but the latter only contains one vibrational band, for the (0,2) transition.

Thus to model the spectrum we allow for a vibrational and rotational temperature for CN, while N$_2$ only is assigned a rotational temperature.

In [6]:

def model_two_species(
    x: np.array,
    mu: float = 0,
    sigma: float = 0.01,
    gamma: float = 0.01,
    T_rot_1: float = 300,
    T_rot_2: float = 300,
    T_vib: float = 300,
    ratio: float = 1,
    A: float = 1,
    b: float = 0,
    resolution: int = 100,
    wl_pad: float = 10,
    sim_db1=None,
    sim_db2=None,
):
    """A model function that combines two line-by-line dataframes for two species.
    The respective contributions between the two species are scaled according to the `ratio` argument"""
    sticks = np.concatenate(
        [
            Moose.create_stick_spectrum(300, T_rot_1, df_db=sim_db1),
            Moose.create_stick_spectrum(T_vib, T_rot_2, pop=ratio, df_db=sim_db2),
        ]
    )
    equid = Moose.equidistant_mesh(sticks, wl_pad, resolution)
    v = Moose.apply_voigt(equid, sigma, gamma)
    matched = Moose.match_spectra((x - mu).reshape(-1, 1), v)
    matched[:, 1] = (matched[:, 1] - matched[:, 1].min()) / (matched[:, 1].max() - matched[:, 1].min())

    return A * matched[:, 1] + b


db1 = Moose.query_DB("N2CB", (350, 400))
db2 = Moose.query_DB("CNBX", (350, 400))

x = np.linspace(370, 390, 1000)
y = model_two_species(x, 0, 0.1, 0.1, 300, 300, 300, 1, sim_db1=db1[db1.Dv==0], sim_db2=db2)

def model_two_species( x: np.array, mu: float = 0, sigma: float = 0.01, gamma: float = 0.01, T_rot_1: float = 300, T_rot_2: float = 300, T_vib: float = 300, ratio: float = 1, A: float = 1, b: float = 0, resolution: int = 100, wl_pad: float = 10, sim_db1=None, sim_db2=None, ): """A model function that combines two line-by-line dataframes for two species. The respective contributions between the two species are scaled according to the `ratio` argument""" sticks = np.concatenate( [ Moose.create_stick_spectrum(300, T_rot_1, df_db=sim_db1), Moose.create_stick_spectrum(T_vib, T_rot_2, pop=ratio, df_db=sim_db2), ] ) equid = Moose.equidistant_mesh(sticks, wl_pad, resolution) v = Moose.apply_voigt(equid, sigma, gamma) matched = Moose.match_spectra((x - mu).reshape(-1, 1), v) matched[:, 1] = (matched[:, 1] - matched[:, 1].min()) / (matched[:, 1].max() - matched[:, 1].min()) return A * matched[:, 1] + b db1 = Moose.query_DB("N2CB", (350, 400)) db2 = Moose.query_DB("CNBX", (350, 400)) x = np.linspace(370, 390, 1000) y = model_two_species(x, 0, 0.1, 0.1, 300, 300, 300, 1, sim_db1=db1[db1.Dv==0], sim_db2=db2)

In [7]:

df = pd.read_csv("CN/370-400nm_9ns_calib.txt").dropna()

model = lmfit.Model(model_two_species, sim_db1=db1, sim_db2=db2)
fit_params = lmfit.create_params(
    **{k: v for k, v in Moose.default_params.items() if k not in ["T_rot"]}
)
fit_params.add("T_rot_1", value=1000, min=300, max=1e4)
fit_params.add("T_rot_2", value=6000, min=300, max=1e4)
fit_params.add("ratio", value=1, min=0, max=100)


fig_fit = go.Figure()
fig_fit.add_scatter(x=df["Wavelength"], y=df["I"], name="Spectrum")
fig_fit.add_scatter(
    x=df["Wavelength"],
    y=model.eval(x=df["Wavelength"].values, params=fit_params),
    name="Fit",
)
fig_fit.add_scatter(x=df["Wavelength"], y=np.zeros(df.shape[0]), name="Residual")

fig_fit.update_xaxes(title_text="λ (nm)")
fig_fit.update_yaxes(title_text="I (a.u.)")

btn_fit = pnw.Button(name="Fit", button_type="primary")
widgets = {
    p: pnw.FloatInput(name=p, value=fit_params[p].value)
    for p in fit_params
    if p not in ["A", "b"]
}


def iter_cb(params, iter, resid, *args, **kwargs):
    for p in params:
        if p not in ["A", "b"]:
            widgets[p].value = params[p].value
    if iter % 10 == 0:
        with fig_fit.batch_update():
            fig_fit.data[2].y = resid
            # fig_fit.data[1].y = model_two_species(**params.valuesdict(),**kwargs, sim_db1=db1,sim_db)
            fig_fit.data[1].y = model.eval(params=params, x=kwargs.get("x"))


def cb_fit(event):
    btn_fit.button_type = "warning"
    result = model.fit(
        df["I"].values, x=df["Wavelength"].values, params=fit_params, iter_cb=iter_cb
    )
    btn_fit.button_type = "primary"
    iter_cb(result.params, -10, result.residual, x=df["Wavelength"].values)


btn_fit.on_click(cb_fit)

layout = pn.Column(
    pn.pane.Plotly(fig_fit, sizing_mode="stretch_width"),
    btn_fit,
    pn.FlexBox(*widgets.values()),
)
layout

df = pd.read_csv("CN/370-400nm_9ns_calib.txt").dropna() model = lmfit.Model(model_two_species, sim_db1=db1, sim_db2=db2) fit_params = lmfit.create_params( **{k: v for k, v in Moose.default_params.items() if k not in ["T_rot"]} ) fit_params.add("T_rot_1", value=1000, min=300, max=1e4) fit_params.add("T_rot_2", value=6000, min=300, max=1e4) fit_params.add("ratio", value=1, min=0, max=100) fig_fit = go.Figure() fig_fit.add_scatter(x=df["Wavelength"], y=df["I"], name="Spectrum") fig_fit.add_scatter( x=df["Wavelength"], y=model.eval(x=df["Wavelength"].values, params=fit_params), name="Fit", ) fig_fit.add_scatter(x=df["Wavelength"], y=np.zeros(df.shape[0]), name="Residual") fig_fit.update_xaxes(title_text="λ (nm)") fig_fit.update_yaxes(title_text="I (a.u.)") btn_fit = pnw.Button(name="Fit", button_type="primary") widgets = { p: pnw.FloatInput(name=p, value=fit_params[p].value) for p in fit_params if p not in ["A", "b"] } def iter_cb(params, iter, resid, *args, **kwargs): for p in params: if p not in ["A", "b"]: widgets[p].value = params[p].value if iter % 10 == 0: with fig_fit.batch_update(): fig_fit.data[2].y = resid # fig_fit.data[1].y = model_two_species(**params.valuesdict(),**kwargs, sim_db1=db1,sim_db) fig_fit.data[1].y = model.eval(params=params, x=kwargs.get("x")) def cb_fit(event): btn_fit.button_type = "warning" result = model.fit( df["I"].values, x=df["Wavelength"].values, params=fit_params, iter_cb=iter_cb ) btn_fit.button_type = "primary" iter_cb(result.params, -10, result.residual, x=df["Wavelength"].values) btn_fit.on_click(cb_fit) layout = pn.Column( pn.pane.Plotly(fig_fit, sizing_mode="stretch_width"), btn_fit, pn.FlexBox(*widgets.values()), ) layout

Out[7]:

Using the 'multi_species_objective' with lmfit.¶

Moose also includes a convenience function to simplify fitting multiple species together, using the lmfit.minimize function, or lmfit.Minimizer class.

This allows you to define an optimization by providing line-by-line data for each species as keyword arguments to the multi_species_objective.

For each of these species (named species), you then add a fraction_{species} parameter and optionally T_vib_{species} and T_rot_{species} parameters, to a lmfit.Parameters instance.

In [8]:

# import the objective, load databases and define parameters for the problem
from Moose.lmfit import multi_species_objective

db_n2 = Moose.query_DB("N2CB", (350, 400))
db_cn = Moose.query_DB("CNBX", (350, 400))

fit_params = lmfit.create_params(**Moose.default_params)

fit_params.add("fraction_N2", value=0.5, min=0, max=1)
fit_params.add("fraction_CN", value=0.5, min=0, max=1)
# For N2 the parameters T_rot and T_vib will be used, for CN we add separate ones.
fit_params.add("T_rot_CN", value=1000, min=300, max=1e4)
fit_params.add("T_vib_CN", value=6000, min=300, max=1e4)

# import the objective, load databases and define parameters for the problem from Moose.lmfit import multi_species_objective db_n2 = Moose.query_DB("N2CB", (350, 400)) db_cn = Moose.query_DB("CNBX", (350, 400)) fit_params = lmfit.create_params(**Moose.default_params) fit_params.add("fraction_N2", value=0.5, min=0, max=1) fit_params.add("fraction_CN", value=0.5, min=0, max=1) # For N2 the parameters T_rot and T_vib will be used, for CN we add separate ones. fit_params.add("T_rot_CN", value=1000, min=300, max=1e4) fit_params.add("T_vib_CN", value=6000, min=300, max=1e4)

In [9]:

# Example of using `multi_species_objective`, as initial guess based on the `fit_params` 
y_init = multi_species_objective(fit_params, x=df["Wavelength"].to_numpy(), N2=db_n2, CN=db_cn, normalize=True)

# Example of using `multi_species_objective`, as initial guess based on the `fit_params` y_init = multi_species_objective(fit_params, x=df["Wavelength"].to_numpy(), N2=db_n2, CN=db_cn, normalize=True)

In [10]:

fig_fit = go.Figure()
fig_fit.add_scatter(x=df["Wavelength"], y=df["I"], name="Spectrum")
fig_fit.add_scatter(
    x=df["Wavelength"],
    y=y_init,
    name="Fit",
)
fig_fit.add_scatter(x=df["Wavelength"], y=np.zeros(df.shape[0]), name="Residual")

fig_fit.update_xaxes(title_text="λ (nm)")
fig_fit.update_yaxes(title_text="I (a.u.)")

btn_fit = pnw.Button(name="Fit", button_type="primary")
widgets = {p: pnw.FloatInput(name=p, value=fit_params[p].value) for p in fit_params if p not in ["A", "b"]}


def iter_cb(params, iter, resid, *args, **kwargs):
    if iter % 10 == 0:
        for p in params:
            if p not in ["A", "b"]:
                widgets[p].value = params[p].value
        with fig_fit.batch_update():
            fig_fit.data[2].y = resid
            fig_fit.data[1].y = multi_species_objective(params, x=kwargs['x'],normalize=True, CN=db_cn,N2=db_n2)


def cb_fit(event):
    btn_fit.button_type = "warning"
    
    # The below function performs the fitting. 
    result = lmfit.minimize(
        multi_species_objective,
        fit_params,
        kws={"x": df["Wavelength"].to_numpy(), "y": df["I"].to_numpy(), "normalize": True, "CN": db_cn, "N2": db_n2},
        iter_cb = iter_cb
    )
    btn_fit.button_type = "primary"
    iter_cb(result.params, -10, result.residual, x=df["Wavelength"].values)


btn_fit.on_click(cb_fit)

layout = pn.Column(
    pn.pane.Plotly(fig_fit, sizing_mode="stretch_width"),
    btn_fit,
    pn.FlexBox(*widgets.values()),
)
layout

fig_fit = go.Figure() fig_fit.add_scatter(x=df["Wavelength"], y=df["I"], name="Spectrum") fig_fit.add_scatter( x=df["Wavelength"], y=y_init, name="Fit", ) fig_fit.add_scatter(x=df["Wavelength"], y=np.zeros(df.shape[0]), name="Residual") fig_fit.update_xaxes(title_text="λ (nm)") fig_fit.update_yaxes(title_text="I (a.u.)") btn_fit = pnw.Button(name="Fit", button_type="primary") widgets = {p: pnw.FloatInput(name=p, value=fit_params[p].value) for p in fit_params if p not in ["A", "b"]} def iter_cb(params, iter, resid, *args, **kwargs): if iter % 10 == 0: for p in params: if p not in ["A", "b"]: widgets[p].value = params[p].value with fig_fit.batch_update(): fig_fit.data[2].y = resid fig_fit.data[1].y = multi_species_objective(params, x=kwargs['x'],normalize=True, CN=db_cn,N2=db_n2) def cb_fit(event): btn_fit.button_type = "warning" # The below function performs the fitting. result = lmfit.minimize( multi_species_objective, fit_params, kws={"x": df["Wavelength"].to_numpy(), "y": df["I"].to_numpy(), "normalize": True, "CN": db_cn, "N2": db_n2}, iter_cb = iter_cb ) btn_fit.button_type = "primary" iter_cb(result.params, -10, result.residual, x=df["Wavelength"].values) btn_fit.on_click(cb_fit) layout = pn.Column( pn.pane.Plotly(fig_fit, sizing_mode="stretch_width"), btn_fit, pn.FlexBox(*widgets.values()), ) layout

Out[10]: