validation.py - master - plasma/compton_ff_yield

validation.py
94 строки3 KB
Sergey Rykovanov
Initial commit
13 июл 2026, 19:42
13 июл 2026, 19:4265c7569
from dataclasses import replace
from pathlib import Path

import numpy as np
import pandas as pd
from numpy.polynomial.hermite import hermgauss
from model import yield_mean, parameters


def export_gh_cases(case, path, nodes=5, betas=(0, .5, .8, .95)):
    x, w = hermgauss(nodes)

    def offsets(sigma):
        return (np.array([0.]), np.array([1.])) if sigma == 0 else (np.sqrt(2)*sigma*x, w/np.sqrt(np.pi))

    dx, wx = offsets(case.sigma_J_perp_m)
    dy, wy = offsets(case.sigma_J_perp_m)
    dz, wz = offsets(case.sigma_JZ_m)
    dt, wt = offsets(case.sigma_Jt_s)
    rows = []
    for beta_ff in betas:
        current = replace(case, beta_ff=beta_ff)
        params = parameters(current)
        analytic_yield = yield_mean(current)
        for x0, ax in zip(dx, wx):
            for y0, ay in zip(dy, wy):
                for z0, az in zip(dz, wz):
                    for t0, at in zip(dt, wt):
                        rows.append({
                            "beta_ff": beta_ff,
                            "delta_x_m": x0,
                            "delta_y_m": y0,
                            "delta_Z_m": z0,
                            "delta_tau_m": t0 * 299792458.0,
                            "delta_t_s": t0,
                            "weight": ax * ay * az * at,
                            "electron_energy_MeV": current.electron_energy_MeV,
                            "charge_pC": current.charge_pC,
                            "epsn_m_rad": current.epsn_m_rad,
                            "beta_star_m": current.beta_star_m,
                            "sigma_le_m": current.sigma_le_m,
                            "laser_energy_J": current.laser_energy_J,
                            "wavelength_m": current.wavelength_m,
                            "a0_peak": current.a0_peak,
                            "sigma_p0_m": current.sigma_p0_m,
                            "sigma_lp_m": params["sigma_lp"],
                            "Z_R_m": params["rayleigh"],
                            "scattering_model": current.scattering_model,
                            "analytic_ensemble_yield": analytic_yield,
                        })
    frame = pd.DataFrame(rows)
    Path(path).parent.mkdir(parents=True, exist_ok=True)
    frame.to_csv(path, index=False)
    return frame


def compare_external_ensemble(exported_cases, external_csv):
    external = pd.read_csv(external_csv)
    required = {
        "beta_ff",
        "delta_x_m",
        "delta_y_m",
        "delta_Z_m",
        "delta_t_s",
        "yield",
    }
    missing = required - set(external.columns)
    if missing:
        raise ValueError(f"External CSV misses columns: {sorted(missing)}")

    keys = ["beta_ff", "delta_x_m", "delta_y_m", "delta_Z_m", "delta_t_s"]
    merged = exported_cases.merge(
        external[keys + ["yield"]],
        on=keys,
        how="inner",
        suffixes=("", "_external"),
    )
    if len(merged) != len(exported_cases):
        raise ValueError("External CSV does not contain all exported cases.")

    grouped = (
        merged.assign(weighted_yield=merged["weight"] * merged["yield"])
        .groupby("beta_ff")
        .agg(
            external_ensemble_yield=("weighted_yield", "sum"),
            analytic_ensemble_yield=("analytic_ensemble_yield", "first"),
            weight_sum=("weight", "sum"),
        )
        .reset_index()
    )
    grouped["relative_difference"] = (
        grouped["external_ensemble_yield"] / grouped["analytic_ensemble_yield"] - 1.0
    )
    return grouped