"""Independent attenuation and single-scattering validation. Keep this file in scripts/ and the case files in ../data/attenuation/. """ from pathlib import Path import argparse import json import less import numpy as np DATA = Path(__file__).resolve().parents[1] / "data" / "attenuation" BACKENDS = ("optix", "vulkan", "embree") def base_scene(config, backend, zenith): scene = less.Scene(backend=backend) scene.sampler = less.Sampler(config.get("sampling_method", "independent")) scene.size = tuple(config["size_m"]) scene.repetitive = config["periodic_wrap_limit"] scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0)) scene.illumination = less.Illumination( source=less.Sun(zenith=zenith, azimuth=90, irradiance=config["normal_irradiance_W_m2_nm"]), atmosphere=less.NoAtmosphere()) return scene def volume_scene(config, backend, zenith, tau, model): scene = base_scene(config, backend, zenith) region = less.Object("layer") region.add_component("volume", less.TurbidBoundary(str(DATA / "layer.obj"))) height = config["layer_height_m"] if model == "medium": prop = less.HomogeneousMedium(absorption=tau / height, scattering=0, temperature=None) else: prop = less.OpticalVegetation( leaf=less.Lambertian(reflectance=0, transmittance=0), leaf_area_density=2 * tau / height, leaf_angle_distribution="spherical") region.set_property("volume", prop) scene.add(region, positions=[[0, 0, 0]]) return scene def leaf_positions(config, seed): settings = config["discrete"] rng = np.random.default_rng(seed) positions = rng.uniform(size=(settings["leaf_count"], 3)) positions[:, :2] *= np.asarray(config["size_m"]) positions[:, 2] = config["layer_bottom_m"] + positions[:, 2] * config["layer_height_m"] # Include the portions of leaves that cross the periodic cell edges. # Ray wrapping alone does not create these neighbouring geometry instances. half_side = settings["leaf_side_m"] / 2 width, depth = config["size_m"] tiled = np.concatenate([positions + [dx, dy, 0] for dx in (-width, 0, width) for dy in (-depth, 0, depth)]) intersects = ((tiled[:, 0] > -half_side) & (tiled[:, 0] < width + half_side) & (tiled[:, 1] > -half_side) & (tiled[:, 1] < depth + half_side)) return tiled[intersects] def discrete_scene(config, backend, zenith, seed): scene = base_scene(config, backend, zenith) leaf = less.Object("leaves", mesh=str(DATA / "leaf.obj")) leaf.set_property(less.Lambertian(reflectance=0, transmittance=0)) scene.add(leaf, positions=leaf_positions(config, seed)) return scene def measure(scene, config, photons, backend): area = float(np.prod(config["size_m"])) product = scene.simulate(less.ForwardRadiationFieldSensor( bands=[config["wavelength_nm"]], mode="grid", grid=(1, 1), layer_heights=[0, .1, .35, .6, .85, 1.1, 1.2], photons_per_m2=int(photons / area))) identity = product.metadata["execution"] assert identity["backend"] == backend and identity["native"] is True assert identity["workflow"] == "optical_forward" assert identity["executor"] == f"{backend}.slang.optical_forward" incident = float(product.incident_energy[0]) transmission = float(product.terrain_absorbed[0]) / incident layers = product.absorption_grid[:, 0, 0, 0].astype(float) / incident absorption = float(layers.sum()) escaped = float(product.escaped_energy[0]) / incident truncated = float(product.truncated_energy[0]) / incident return { "transmission": transmission, "absorption": absorption, "layer_absorption": layers.tolist(), "escaped": escaped, "truncated": truncated, "incident_W_nm": incident, "closure_error": abs(transmission + absorption + escaped + truncated - 1), "execution": {key: identity[key] for key in ("backend", "workflow", "executor", "native", "random_seed")}, } def run(config): assert config["size_m"] == [2, 2] assert config["layer_bottom_m"] == .1 and config["layer_height_m"] == 1 assert config["discrete"]["leaf_side_m"] == .2 thresholds = config["thresholds"] records = [] for backend in BACKENDS: for model in ("medium", "turbid_spherical"): for zenith in config["zenith_degrees"]: cosine = float(np.cos(np.deg2rad(zenith))) for tau in config["vertical_optical_depths"]: scene = volume_scene(config, backend, zenith, tau, model) reference = float(np.exp(-tau / cosine)) depths = np.array([1, 1, .75, .5, .25, 0, 0]) profile = np.diff(np.exp(-tau * depths / cosine)) previous = None for photons in config["photon_budgets"]: values = measure(scene, config, photons, backend) error = max(abs(values["transmission"] - reference), abs(values["absorption"] - (1 - reference)), float(np.max(np.abs(values["layer_absorption"] - profile)))) current = np.r_[values["transmission"], values["absorption"], values["layer_absorption"]] change = None if previous is None else float(np.max(abs(current - previous))) passed = (error <= thresholds["fraction_absolute_error"] and values["closure_error"] <= thresholds["closure_absolute_error"] and values["escaped"] <= thresholds["escaped_fraction"] and values["truncated"] <= thresholds["truncated_fraction"] and (change is None or change <= thresholds["budget_change_absolute"])) records.append(dict(backend=backend, model=model, zenith=zenith, tau=tau, photons=photons, reference=reference, reference_layers=profile.tolist(), error=error, budget_change=change, passed=passed, **values)) previous = current print(backend, model, zenith, "completed", flush=True) paired = [] final = [r for r in records if r["photons"] == config["photon_budgets"][-1]] for row in final: if row["model"] != "medium": continue match = next(r for r in final if r["backend"] == row["backend"] and r["model"] == "turbid_spherical" and r["zenith"] == row["zenith"] and r["tau"] == row["tau"]) difference = max(abs(row[key] - match[key]) for key in ("transmission", "absorption")) paired.append(dict(backend=row["backend"], zenith=row["zenith"], tau=row["tau"], difference=difference, passed=difference <= thresholds["paired_representation_absolute"])) settings = config["discrete"] probability = settings["leaf_side_m"] ** 2 / np.prod(config["size_m"]) discrete_reference = float((1 - probability) ** settings["leaf_count"]) discrete = [] for backend in BACKENDS: for zenith in config["zenith_degrees"]: samples = {budget: [] for budget in settings["photon_budgets"]} controls = [] for seed in settings["geometry_seeds"]: scene = discrete_scene(config, backend, zenith, seed) for photons in settings["photon_budgets"]: values = measure(scene, config, photons, backend) samples[photons].append(values["transmission"]) controls.append(values) means = {budget: float(np.mean(values)) for budget, values in samples.items()} final_samples = samples[settings["photon_budgets"][-1]] sem = float(np.std(final_samples, ddof=1) / np.sqrt(len(final_samples))) error = abs(means[settings["photon_budgets"][-1]] - discrete_reference) change = abs(means[settings["photon_budgets"][-1]] - means[settings["photon_budgets"][0]]) maximum = lambda key: max(v[key] for v in controls) limit = thresholds["fraction_absolute_error"] passed = (error <= limit and error <= thresholds["ensemble_standard_errors"] * sem and change <= thresholds["budget_change_absolute"] and maximum("closure_error") <= thresholds["closure_absolute_error"] and maximum("escaped") <= thresholds["escaped_fraction"] and maximum("truncated") <= thresholds["truncated_fraction"]) discrete.append(dict(backend=backend, zenith=zenith, reference=discrete_reference, budget_means=means, final_scene_transmissions=final_samples, standard_error=sem, error=error, acceptance_limit=limit, budget_change=change, closure_error=maximum("closure_error"), escaped=maximum("escaped"), truncated=maximum("truncated"), execution=controls[-1]["execution"], passed=passed)) print(backend, "discrete", zenith, "error", error, "SEM", sem, flush=True) return dict(less_version=less.__version__, input=config, records=records, paired=paired, discrete=discrete, passed=all(r["passed"] for r in final + paired + discrete)) def scattering_cases(config): settings = config["scattering"] cases = [] def add(sun, signed, tau, albedo, g, cross=False): relative = (settings["cross_plane_view"]["relative_azimuth"] if cross else 180 if signed < 0 else 0) cases.append(dict(sun_zenith=sun, view_zenith=abs(signed), relative_azimuth=relative, signed_view=None if cross else signed, tau=tau, albedo=albedo, g=g)) for sun in settings["sun_zenith_degrees"]: for g in settings["phase_g"]: for signed in settings["signed_view_degrees"]: add(sun, signed, settings["sweep_tau"], settings["sweep_single_scattering_albedo"], g) add(sun, settings["cross_plane_view"]["zenith"], settings["sweep_tau"], settings["sweep_single_scattering_albedo"], g, cross=True) for tau in settings["control_optical_depths"]: for signed in settings["control_signed_views"]: add(settings["control_sun_zenith"], signed, tau, settings["sweep_single_scattering_albedo"], 0) for albedo in settings["control_single_scattering_albedos"]: for signed in settings["control_signed_views"]: add(settings["control_sun_zenith"], signed, settings["sweep_tau"], albedo, 0) return cases def scattering_reference(config, case): """Exact first-order radiance of an infinite plane-parallel HG layer.""" sun, view, delta = np.deg2rad([ case["sun_zenith"], case["view_zenith"], case["relative_azimuth"]]) mu0, mu = np.cos(sun), np.cos(view) cosine = -(mu0*mu+np.sin(sun)*np.sin(view)*np.cos(delta)) g = case["g"] phase = (1-g*g)/(4*np.pi*(1+g*g-2*g*cosine)**1.5) layer = -np.expm1(-case["tau"]*(1/mu0+1/mu)) brf = np.pi*case["albedo"]*phase/(mu0+mu)*layer radiance = config["normal_irradiance_W_m2_nm"]*mu0/np.pi*brf return dict(cos_scattering=float(cosine), phase_sr_inverse=float(phase), radiance=float(radiance), brf=float(brf)) def run_scattering(config): from scipy.special import roots_legendre from scipy.stats import t settings = config["scattering"] cases = scattering_cases(config) groups = settings["seed_groups"] critical = float(t.ppf(1-settings["family_error_rate"]/(2*len(cases)*len(BACKENDS)), groups-1)) nodes, weights = roots_legendre(128) phase_checks = [] for g in settings["phase_g"]: values = (1-g*g)/(4*np.pi*(1+g*g-2*g*nodes)**1.5) norm = float(2*np.pi*np.sum(weights*values)) mean = float(2*np.pi*np.sum(weights*nodes*values)) assert abs(norm-1) <= settings["thresholds"]["phase_reference_normalization_error"] assert abs(mean-g) <= settings["thresholds"]["phase_reference_mean_cosine_error"] phase_checks.append(dict(g=g, normalization=norm, mean_cosine=mean)) records, convergence, executions = [], [], {} limits = settings["thresholds"] for backend in BACKENDS: scene = base_scene(config, backend, settings["sun_zenith_degrees"][0]) layer = less.Object("scattering_layer") layer.add_component("volume", less.TurbidBoundary(str(DATA/"layer.obj"))) scene.add(layer, positions=[(0, 0, 0)]) accumulators = {n: dict(abs_errors=[], rel_errors=[], normalized_groups=[], variances=[]) for n in settings["qualities"]} for index, case in enumerate(cases): reference = scattering_reference(config, case) scene.illumination = less.Illumination( source=less.Sun(zenith=case["sun_zenith"], azimuth=settings["sun_azimuth_degrees"], irradiance=config["normal_irradiance_W_m2_nm"]), atmosphere=less.NoAtmosphere()) extinction = case["tau"]/config["layer_height_m"] layer.set_property("volume", less.HomogeneousMedium( absorption=extinction*(1-case["albedo"]), scattering=extinction*case["albedo"], g=case["g"], temperature=None)) previous = None for quality in settings["qualities"]: sensor = less.OpticalImager(less.Orthographic( image_size=tuple(settings["image_size"]), extent=settings["extent_m"], center=(1, 1), view_zenith=case["view_zenith"], view_azimuth=(settings["sun_azimuth_degrees"]+case["relative_azimuth"])%360), bands=[config["wavelength_nm"]], quality=quality, outputs=("radiance", "scattering_order"), max_output_order=1) image = scene.simulate(sensor) identity = image.metadata["execution"] assert identity["backend"] == backend and identity["native"] is True assert identity["workflow"] == "optical_backward" assert identity["executor"] == f"{backend}.slang.optical_backward" executions[backend] = {key: identity[key] for key in ("backend", "workflow", "executor", "native")} np.testing.assert_allclose(image.scattering_orders[0], 0, atol=1e-10) samples = np.asarray(image.scattering_orders[1], float).reshape(groups, -1).mean(axis=1) factor = np.pi/(config["normal_irradiance_W_m2_nm"] *np.cos(np.deg2rad(case["sun_zenith"]))) values = samples*factor mean = float(values.mean()) sem = float(values.std(ddof=1)/np.sqrt(groups)) absolute = abs(mean-reference["brf"]) relative = absolute/reference["brf"] accumulator = accumulators[quality] accumulator["abs_errors"].append(absolute) accumulator["rel_errors"].append(relative) accumulator["normalized_groups"].extend((values/reference["brf"]-1).tolist()) accumulator["variances"].append(float(np.var(values/reference["brf"], ddof=1))) change = None if previous is None else abs(mean-previous) if quality == settings["qualities"][-1]: interval = [mean-critical*sem, mean+critical*sem] in_interval = interval[0] <= reference["brf"] <= interval[1] passed = (absolute <= limits["brf_absolute_error"] and relative <= limits["radiance_relative_error"] and change <= limits["last_budget_brf_change"] and in_interval) records.append(dict(backend=backend, **case, reference=reference, brf=mean, radiance=float(samples.mean()), sem_brf=sem, group_brf=values.tolist(), confidence_interval_brf=interval, absolute_error=absolute, relative_error=relative, budget_change=change, passed=bool(passed))) previous = mean print(backend, index+1, "/", len(cases), case["g"], case["sun_zenith"], case["signed_view"], case["relative_azimuth"], absolute, relative, passed, flush=True) for quality, accumulator in accumulators.items(): convergence.append(dict(backend=backend, quality=quality, max_absolute_error=max(accumulator["abs_errors"]), max_relative_error=max(accumulator["rel_errors"]), normalized_group_rms=float(np.sqrt(np.mean(np.square(accumulator["normalized_groups"])))), normalized_group_variance=float(np.mean(accumulator["variances"])))) trend_passed = all( next(r for r in convergence if r["backend"] == b and r["quality"] == settings["qualities"][-1])[key] < next(r for r in convergence if r["backend"] == b and r["quality"] == settings["qualities"][0])[key] for b in BACKENDS for key in ("normalized_group_rms", "normalized_group_variance")) return dict(records=records, convergence=convergence, executions=executions, phase_reference_checks=phase_checks, critical_t=critical, simultaneous_comparisons=len(cases)*len(BACKENDS), trend_passed=bool(trend_passed), passed=all(r["passed"] for r in records) and bool(trend_passed)) def scattering_figure(result, output): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"svg.fonttype": "none", "font.size": 10}) colors = dict(optix="#0072B2", vulkan="#D55E00", embree="#009E73") markers = dict(optix="o", vulkan="x", embree="+") config, data = result["input"], result["scattering"] settings = config["scattering"] fig, axes = plt.subplots(2, 2, figsize=(11, 7), layout="constrained") for sun, ax in zip(settings["sun_zenith_degrees"], axes[0]): for g, style in zip(settings["phase_g"], ("-", "--", ":")): angles = np.linspace(-60, 60, 181) references = [scattering_reference(config, dict(sun_zenith=sun, view_zenith=abs(a), relative_azimuth=180 if a<0 else 0, tau=settings["sweep_tau"], albedo=settings["sweep_single_scattering_albedo"], g=g))["brf"] for a in angles] ax.plot(angles, references, color="black", ls=style, label=f"Analytic g={g:g}") for backend in BACKENDS: rows = [r for r in data["records"] if r["backend"]==backend and r["sun_zenith"]==sun and r["g"]==g and r["signed_view"] is not None and r["tau"]==settings["sweep_tau"] and r["albedo"]==settings["sweep_single_scattering_albedo"]] ax.plot([r["signed_view"] for r in rows], [r["brf"] for r in rows], ls="none", marker=markers[backend], color=colors[backend], label=backend if g==settings["phase_g"][0] else None) ax.set(title=f"Sun zenith {sun}°; optical depth 1; albedo 0.7", xlabel="Signed view zenith (degrees)", ylabel="First-order BRF") ax.grid(alpha=.2) ax.legend(fontsize=8, frameon=False) for backend in BACKENDS: rows = [r for r in data["convergence"] if r["backend"]==backend] axes[1, 0].plot([r["quality"] for r in rows], [r["max_absolute_error"] for r in rows], marker=markers[backend], color=colors[backend], label=backend) axes[1, 1].plot([r["quality"] for r in rows], [r["normalized_group_rms"] for r in rows], marker=markers[backend], color=colors[backend], label=backend) axes[1, 0].axhline(settings["thresholds"]["brf_absolute_error"], color="k", ls="--") axes[1, 0].set(xscale="log", yscale="log", xlabel="Samples per pixel", ylabel="Maximum absolute BRF error", title="All parameter and view combinations") axes[1, 1].set(xscale="log", yscale="log", xlabel="Samples per pixel", ylabel="Normalized group RMS error", title="Eight independent pixel groups") for ax in axes[1]: ax.legend(frameon=False) ax.grid(alpha=.2) fig.savefig(output/"scattering.svg") fig.savefig(output/"scattering.png", dpi=140) plt.close(fig) def figures(result, output): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"svg.fonttype": "none", "font.size": 10}) colors = dict(optix="#0072B2", vulkan="#D55E00", embree="#009E73") config = result["input"] fig, axes = plt.subplots(1, 3, figsize=(12, 3.7), layout="constrained") final = [r for r in result["records"] if r["photons"] == config["photon_budgets"][-1]] for zenith, ax in zip(config["zenith_degrees"], axes): tau = np.linspace(0, 5, 200) ax.plot(tau, np.exp(-tau / np.cos(np.deg2rad(zenith))), "k-", label="Analytic") for backend in BACKENDS: for model, marker in (("medium", "o"), ("turbid_spherical", "+")): rows = [r for r in final if r["backend"] == backend and r["model"] == model and r["zenith"] == zenith] ax.plot([r["tau"] for r in rows], [r["transmission"] for r in rows], marker, color=colors[backend], fillstyle="none", label=f"{backend}: {model.split('_')[0]}") ax.set(xlabel="Vertical optical depth", ylabel="Transmission", title=f"Sun zenith {zenith} deg") ax.grid(alpha=.2) axes[0].legend(fontsize=7) fig.savefig(output / "transmission.svg"); plt.close(fig) fig, axes = plt.subplots(1, 2, figsize=(10, 3.7), layout="constrained") for backend, color in colors.items(): rows = [r for r in result["records"] if r["backend"] == backend] errors = [max(r["error"] for r in rows if r["photons"] == budget) for budget in config["photon_budgets"]] axes[0].plot(config["photon_budgets"], errors, "o-", color=color, label=backend) rows = [r for r in result["discrete"] if r["backend"] == backend] axes[1].errorbar([r["zenith"] for r in rows], [float(list(r["budget_means"].values())[-1]) for r in rows], yerr=[4 * r["standard_error"] for r in rows], fmt="o", capsize=4, color=color, label=backend) axes[0].axhline(.005, ls="--", color="#444444", label="Error limit") axes[0].set(xscale="log", xlabel="Photon budget", ylabel="Max absolute fraction error", title="Medium and turbid budget refinement") reference = result["discrete"][0]["reference"] axes[1].axhline(reference, color="black", label="Finite-leaf analytic mean") axes[1].axhline(np.exp(-1), color="#555555", ls="--", label="Infinitesimal-leaf limit") axes[1].set(xlabel="Sun zenith (deg)", ylabel="Mean transmission (+/- 4 SEM)", title="12 independent horizontal-leaf scenes") for ax in axes: ax.grid(alpha=.2); ax.legend(fontsize=8) fig.savefig(output / "sampling.svg"); plt.close(fig) fig, axes = plt.subplots(1, 2, figsize=(10, 3.7), layout="constrained") selected = next(r for r in final if r["tau"] == 1 and r["zenith"] == 45) heights = [.225, .475, .725, .975] axes[0].plot(selected["reference_layers"][1:5], heights, "k-", label="Analytic") for backend, color in colors.items(): for model, marker in (("medium", "o"), ("turbid_spherical", "+")): row = next(r for r in final if r["backend"] == backend and r["model"] == model and r["tau"] == 1 and r["zenith"] == 45) axes[0].plot(row["layer_absorption"][1:5], heights, marker, color=color, fillstyle="none", label=f"{backend}: {model.split('_')[0]}") axes[0].set(xlabel="Absorbed fraction per 0.25 m layer", ylabel="Height (m)", title="Depth profile; optical depth 1; sun 45 deg") zenith = np.linspace(0, 75, 200) transmission = np.exp(-1 / np.cos(np.deg2rad(zenith))) axes[1].plot(zenith, transmission, "k-", label="Analytic ground absorption") axes[1].plot(zenith, 1 - transmission, "k--", label="Analytic layer absorption") for backend, color in colors.items(): rows = [r for r in final if r["backend"] == backend and r["model"] == "medium" and r["tau"] == 1] for quantity, marker in (("transmission", "o"), ("absorption", "+")): axes[1].plot([r["zenith"] for r in rows], [r[quantity] for r in rows], marker, color=color, fillstyle="none", label=backend if quantity == "transmission" else None) axes[1].set(xlabel="Sun zenith (deg)", ylabel="Fraction of incident flux", title="Energy partition; optical depth 1") for ax in axes: ax.grid(alpha=.2); ax.legend(fontsize=7) fig.savefig(output / "absorption.svg"); plt.close(fig) from mpl_toolkits.mplot3d.art3d import Poly3DCollection fig = plt.figure(figsize=(7, 4.5), layout="constrained") ax = fig.add_subplot(projection="3d", computed_zorder=False) centers = leaf_positions(config, config["discrete"]["geometry_seeds"][0]) offsets = np.array([[-.1, -.1, 0], [.1, -.1, 0], [.1, .1, 0], [-.1, .1, 0]]) ax.add_collection3d(Poly3DCollection(centers[:, None, :] + offsets, facecolors="#344d3a", edgecolors="#1c2c20", linewidths=.3, zorder=2)) for z in (.1, 1.1): ax.plot([0, 2, 2, 0, 0], [0, 0, 2, 2, 0], z, color="#6086a4", zorder=1) for x, y in ((0, 0), (0, 2), (2, 0), (2, 2)): ax.plot([x, x], [y, y], [.1, 1.1], color="#6086a4", zorder=1) ax.quiver(.8, .8, 1.7, 0, 0, -.45, color="#ce641c", linewidth=2, zorder=5) ax.set(xlim=(0, 2), ylim=(0, 2), zlim=(0, 1.8), xlabel="X / m", ylabel="Y / m", zlabel="Z / m", title="Horizontal black leaves; LAI = 1; seed = 101") ax.view_init(elev=24, azim=-65) fig.savefig(output / "scene.svg"); plt.close(fig) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--output", type=Path, default=Path("build/attenuation")) parser.add_argument("--scattering-only", action="store_true") parser.add_argument("--plot-only", type=Path) args = parser.parse_args() args.output.mkdir(parents=True, exist_ok=True) if args.plot_only: result = json.loads(args.plot_only.read_text(encoding="utf-8")) scattering_figure(result, args.output) return config = json.loads((DATA / "input.json").read_text(encoding="utf-8")) if args.scattering_only: result = json.loads((DATA/"results.json").read_text(encoding="utf-8")) assert {k:v for k,v in result["input"].items() if k!="scattering"} == { k:v for k,v in config.items() if k!="scattering"} else: result = run(config) absorption_passed = result.get("absorption_passed", result["passed"]) result["input"] = config result["scattering"] = run_scattering(config) result["absorption_passed"] = bool(absorption_passed) result["passed"] = bool(absorption_passed and result["scattering"]["passed"]) (args.output / "results.json").write_text( json.dumps(result, indent=2, allow_nan=False), encoding="utf-8") if not args.scattering_only: figures(result, args.output) scattering_figure(result, args.output) raise SystemExit(0 if result["passed"] else 1) if __name__ == "__main__": main()