Variance Report

Variance Report#

Overview

In this tutorial, we introduce variance report in Eradiate.

Prerequisites

Ability to run 1D simulations and visualize the results (see eradiate_quickstart).

What you will learn

How to turn on variance report for a simulation.
How to access and visualize the variance reported by a simulation.

In this tutorial, we will show you how to use Eradiate to retrieve the variance of the results of your simulations. Let’s start by running a simple 1D simulation and viewing its result.

[2]:

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import eradiate
from eradiate import unit_registry as ureg
from eradiate import scenes as ertsc

eradiate.set_mode("mono")

[3]:

def run_exp(spp=10000, integrator=None):
    """
    Generate and run an atmosphere experiment over a hemisphere plane cut.
    """

    if integrator is None:
        integrator = ertsc.integrators.VolPathIntegrator()

    exp = eradiate.experiments.AtmosphereExperiment(
        surface={
            "type":"lambertian",
            "reflectance":0.5
        },
        illumination={
            "type": "directional",
            "zenith": 30,
            "azimuth": 0,
        },
        atmosphere={
            "type": "molecular",
            "has_absorption":False
        },
        measures= ertsc.measure.MultiDistantMeasure.hplane(
            zeniths=np.arange(-85,86,5),
            azimuth=0.,
            spp = spp,
            srf={"type": "multi_delta", "wavelengths": 550*ureg.nm},
        ),
        integrator=integrator,
    )

    return eradiate.run(exp,spp)

result = run_exp(spp=10000)

result

By default, variance is not available in the results. Variance report can be activated by setting the moment parameter of any integrator to True:

[10]:

integrator = ertsc.integrators.PiecewiseVolPathIntegrator(moment=True)

We can now run the simulations again, with our new integrator.

[11]:

var_result = run_exp(spp=10000, integrator=integrator)
var_result

Note the new radiance_var data variable, which is the radiance variance of our simulation. We can visualize it by plotting two standard deviations around the resulting radiance.

[19]:

# Get the standard deviation from the variance.
std = np.sqrt(var_result.radiance_var)

[23]:

fig, ax = plt.subplots(1,1,figsize=(1*4.5,4), layout="constrained")

upper = (var_result.radiance + 2*std).squeeze()
lower = (var_result.radiance - 2*std).squeeze()
vza = var_result.vza.squeeze()

ax.plot(var_result.vza, var_result.radiance.squeeze().values, linestyle=":", marker=".")
ax.set_xlabel("viewing zenith angle [deg]")
ax.ticklabel_format(axis='y', scilimits=[-1,1])
ax.fill_between(vza, upper, lower, alpha=0.2, color="C0")
ax.set_ylabel("radiance")
plt.show()
plt.close()

../../_images/tutorials_howto_variance_report_11_0.png

We can repeat this process a final time with a higher sample count and see the variance go down.

[9]:

var_result = run_exp(spp=100000, integrator = integrator)

std = np.sqrt(var_result.radiance_var)

fig, ax = plt.subplots(1,1,figsize=(1*4.5,4), layout="constrained")

upper = (var_result.radiance + 2*std).squeeze()
lower = (var_result.radiance - 2*std).squeeze()
vza = var_result.vza.squeeze()

ax.plot(var_result.vza, var_result.radiance.squeeze().values, linestyle=":", marker=".")
ax.set_xlabel("viewing zenith angle [deg]")
ax.ticklabel_format(axis='y', scilimits=[-1,1])
ax.fill_between(vza, upper, lower, alpha=0.2, color="C0")
ax.set_ylabel("radiance")
plt.show()
plt.close()

../../_images/tutorials_howto_variance_report_13_0.png

Final words#

In this tutorial, you have learned how to turn on variance report and access and visualize radiance variance. Note that variance is also available for Stokes vectors when running in polarization mode (see polarization).

Variance Report

Contents

Variance Report#

Final words#