Resolve masure_regression airfoils to polars with a pure-Julia Extra-Trees evaluator

CHANGELOG.md

@@ -4,6 +4,7 @@
 
 ### Added
 
+- `masure_regression` airfoils `{t, eta, kappa, delta, lambda, phi}` resolve to polars: `resolve_aero_geometry(...; ml_models_dir)` evaluates the Masure Extra-Trees regression (`masure_aero`, `load_masure_model`) in pure Julia, from models converted once with `scripts/export_masure_models.py`.
 - `plot_geometry`, `plot_distribution`, `plot_combined_analysis`, `plot_section_polars`,
   `plot_airfoil_fit` and `plot_airfoils` take `show_title=true`; with `false` the title
   is not drawn, and still names the saved file and window.

data/TUDELFT_V3_KITE/aero_geometry.yaml

@@ -75,7 +75,7 @@ wing_airfoils:
   #   - polars:
   #       polar_file_path: Path to polar CSV file (columns: alpha [rad], cl, cd, cm)
   #   - masure_regression:
-  #       t, eta, kappa, delta, lamba, phi: Regression parameters
+  #       t, eta, kappa, delta, lambda, phi: Regression parameters
   #   - inviscid:
   #       no further data is required 
   # ---------------------------------------------------------------

data/TUDELFT_V3_KITE/aero_geometry_coarse_discretisation.yaml

@@ -120,7 +120,7 @@ wing_airfoils:
   # ---------------------------------------------------------------
   # headers:
   #   - airfoil_id: integer, unique identifier for the airfoil
-  #   - type: one of [neuralfoil, breukels_regression, measure_regression, polars]
+  #   - type: one of [neuralfoil, breukels_regression, masure_regression, polars]
   #   - info_dict: dictionary with parameters depending on 'type'
   #
   # info_dict fields by type:
@@ -141,7 +141,7 @@ wing_airfoils:
   #           Dirty wind tunnel:   4–8
   #   - polars:
   #       polar_file_path: Path to polar CSV file (columns: alpha [rad], cl, cd, cm)
-  #   - measure_regression:
+  #   - masure_regression:
   #       t, eta, kappa, delta, lambda, phi: Regression parameters
   #   - inviscid:
   #       no further data is required

docs/src/airfoil_pipeline.md

@@ -124,6 +124,10 @@ keeping its own edge positions, and a warning lists the reuse. All floats are ro
 precision by the single [`write_yaml`](@ref VortexStepMethod.ObjAdapter.write_yaml) writer,
 so generated geometry files stay diff-friendly and consistent.
 
+## Parametric LEI airfoils: the masure regression
+
+A geometry YAML can also describe a leading-edge-inflatable airfoil by six shape parameters instead of a slice, as the `wing_airfoils` entry `[id, masure_regression, {t, eta, kappa, delta, lambda, phi}]`: tube diameter, chordwise and vertical position of maximum camber, trailing-edge reflex angle [deg], camber tension and leading-edge tension. [`resolve_aero_geometry`](@ref VortexStepMethod.ObjAdapter.resolve_aero_geometry) turns it into a `POLAR_VECTORS` table over the block's `alpha_range` with [`masure_aero`](@ref), which evaluates the Extra-Trees regression of K.R.G. Masure, trained on 2D RANS simulations, in pure Julia. The models exist for `reynolds` 1e6, 5e6 and 2e7. Download them from [Zenodo](https://doi.org/10.5281/zenodo.16925758), convert them once with `python scripts/export_masure_models.py ET_re*.pkl --out DIR` (needs numpy and scikit-learn), and pass `ml_models_dir=DIR`. The converted model gives the same coefficients as scikit-learn's `predict`. Keep `alpha_range` within the angles the models were trained on: outside them the trees return a constant.
+
 ## Why this is useful
 
 - **CAD in, solver-ready model out.** No hand-drawing airfoils or manually pairing them

docs/src/functions.md

@@ -41,6 +41,8 @@ deform_section
 analyze_section
 analyze_sweep
 neuralfoil_aero
+masure_aero
+load_masure_model
 deform_kulfan
 chord_residual
 chord_line

docs/src/private_functions.md

@@ -228,6 +228,15 @@ swish
 sigmoid
 ```
 
+### Masure regression
+```@docs
+MASURE_PARAMETERS
+MASURE_REYNOLDS
+forest_predict
+one_based
+scale_input
+```
+
 ### Polars and airfoil IO
 ```@docs
 create_2d_polars

docs/src/private_types.md

@@ -30,6 +30,8 @@ KulfanBasis
 LivePolarSettings
 LivePolars
 NeuralFoilModel
+MasureModel
+ExtraTreesForest
 NeuralFoilResult
 NeuralFoilWorkspace
 ContourPressureScratch

scripts/export_masure_models.py

@@ -0,0 +1,125 @@
+"""Export the masure-regression scikit-learn models to the .npz files AirfoilAero reads.
+
+    python export_masure_models.py ET_re1e6.pkl ET_re5e6.pkl ET_re2e7.pkl --out DIR
+    python export_masure_models.py --fixture DIR
+
+The first form converts the pickles from https://doi.org/10.5281/zenodo.16925758 into
+`DIR/ET_re<Re>.npz`; pass `DIR` as `ml_models_dir` to `resolve_aero_geometry`. The
+second trains a small model of the same layout on synthetic data and writes it with
+the inputs and `predict` outputs the Julia test compares against. Needs numpy and
+scikit-learn.
+
+An exported file holds `input_mean` and `input_scale` (the StandardScaler) and, for
+output k = 0, 1, 2 (CD, CL, CM), the nodes of all its trees concatenated:
+`output{k}_roots`, `_left`, `_right`, `_feature`, `_threshold`, `_value`. Node and
+feature indices are 0-based, and `_left` is -1 on a leaf.
+"""
+
+import argparse
+import pickle
+from pathlib import Path
+
+import numpy as np
+from sklearn.ensemble import ExtraTreesRegressor
+from sklearn.multioutput import MultiOutputRegressor
+from sklearn.pipeline import Pipeline
+from sklearn.preprocessing import StandardScaler
+
+
+def forest_arrays(forest):
+    """Concatenate the nodes of every tree in `forest`, children indexed globally."""
+    roots, left, right, feature, threshold, value = [], [], [], [], [], []
+    offset = 0
+    for estimator in forest.estimators_:
+        tree = estimator.tree_
+        if tree.n_outputs != 1:
+            raise ValueError(f"expected single-output trees, got {tree.n_outputs}")
+        is_leaf = tree.children_left == -1
+        roots.append(offset)
+        left.append(np.where(is_leaf, -1, tree.children_left + offset))
+        right.append(np.where(is_leaf, -1, tree.children_right + offset))
+        feature.append(tree.feature)
+        threshold.append(tree.threshold)
+        value.append(tree.value[:, 0, 0])
+        offset += tree.node_count
+    return {
+        "roots": np.asarray(roots, dtype=np.int32),
+        "left": np.concatenate(left).astype(np.int32),
+        "right": np.concatenate(right).astype(np.int32),
+        "feature": np.concatenate(feature).astype(np.int32),
+        "threshold": np.concatenate(threshold).astype(np.float64),
+        "value": np.concatenate(value).astype(np.float64),
+    }
+
+
+def export_model(model, npz_path):
+    """Write a StandardScaler -> MultiOutputRegressor(ExtraTreesRegressor) pipeline."""
+    steps = [step for _, step in model.steps] if isinstance(model, Pipeline) else []
+    if not (len(steps) == 2 and isinstance(steps[0], StandardScaler)
+            and isinstance(steps[1], MultiOutputRegressor)):
+        raise ValueError(f"unexpected model layout: {model!r}")
+    scaler, multi_output = steps
+    if len(multi_output.estimators_) != 3:
+        raise ValueError(f"expected 3 outputs, got {len(multi_output.estimators_)}")
+    arrays = {"input_mean": scaler.mean_, "input_scale": scaler.scale_}
+    for k, forest in enumerate(multi_output.estimators_):
+        if not isinstance(forest, ExtraTreesRegressor):
+            raise ValueError(f"output {k} is a {type(forest).__name__}")
+        for name, array in forest_arrays(forest).items():
+            arrays[f"output{k}_{name}"] = array
+    np.savez(npz_path, **arrays)
+
+
+def threshold_rows(model, row):
+    """Copies of `row` moved onto each output's first root split, where float32 matters."""
+    scaler, multi_output = (step for _, step in model.steps)
+    rows = []
+    for forest in multi_output.estimators_:
+        tree = forest.estimators_[0].tree_
+        feature, threshold = tree.feature[0], tree.threshold[0]
+        for step in (-1e-9, 0.0, 1e-9):
+            moved = row.copy()
+            moved[feature] = threshold * (1 + step) * scaler.scale_[feature] + \
+                scaler.mean_[feature]
+            rows.append(moved)
+    return np.array(rows)
+
+
+def write_fixture(out_dir):
+    """Train a small pipeline of the published layout and write it with its predictions."""
+    rng = np.random.default_rng(42)
+    low = np.array([0.05, 0.1, 0.0, -10.0, 0.1, 0.1, -10.0])
+    high = np.array([0.12, 0.6, 0.15, 5.0, 0.4, 0.9, 30.0])
+    X = rng.uniform(low, high, size=(400, 7))
+    alpha = np.deg2rad(X[:, 6])
+    y = np.column_stack([0.02 + 0.3 * alpha**2 + X[:, 0],
+                         2 * np.pi * alpha + 5 * X[:, 2],
+                         -0.1 - X[:, 2] + 0.01 * X[:, 3]])
+    model = Pipeline([
+        ("scale", StandardScaler()),
+        ("model", MultiOutputRegressor(ExtraTreesRegressor(
+            n_estimators=5, max_depth=6, max_features="log2", random_state=42))),
+    ]).fit(X, y)
+    out_dir.mkdir(parents=True, exist_ok=True)
+    export_model(model, out_dir / "ET_re1e6.npz")
+    X_test = np.vstack([rng.uniform(low, high, size=(20, 7)), threshold_rows(model, X[0])])
+    np.savez(out_dir / "reference.npz", X=X_test, Y=model.predict(X_test))
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
+    parser.add_argument("pickles", nargs="*", type=Path)
+    parser.add_argument("--out", type=Path, default=Path("."))
+    parser.add_argument("--fixture", type=Path)
+    args = parser.parse_args()
+    if args.fixture is not None:
+        write_fixture(args.fixture)
+    for pickle_path in args.pickles:
+        with open(pickle_path, "rb") as file:
+            model = pickle.load(file)
+        args.out.mkdir(parents=True, exist_ok=True)
+        export_model(model, args.out / pickle_path.with_suffix(".npz").name)
+
+
+if __name__ == "__main__":
+    main()

src/airfoil_aero/AirfoilAero.jl

@@ -16,6 +16,7 @@ include("kulfan.jl")
 include("deform.jl")
 include("shrink_wrap.jl")
 include("neuralfoil.jl")
+include("masure_regression.jl")
 include("poly.jl")
 include("airfoil_solvers/common.jl")
 include("airfoil_solvers/xfoil_solver.jl")
@@ -33,6 +34,7 @@ export ShrinkWrap, shrink_wrap
 export fit_kulfan_parameters, kulfan_to_coordinates
 export NeuralFoilModel, NeuralFoilResult, load_neuralfoil_model
 export neuralfoil_aero, neuralfoil_section
+export MasureModel, load_masure_model, masure_aero
 export KulfanBasis, deform_kulfan, control_point_deflection
 export chord_residual, chord_line
 export LivePolarSettings, LivePolars, panel_kulfan_parameters

src/airfoil_aero/airfoil_io.jl

@@ -16,15 +16,10 @@ end
 """
     write_polar(filepath, result::NeuralFoilResult)
 
-Write a NeuralFoil result to a `POLAR_VECTORS` table (`alpha, Cd, Cs, Cl, Cm`), CSV or
-Arrow as the suffix of `filepath` says (see [`write_node_rows`](@ref
-VortexStepMethod.write_node_rows)).
+Write a NeuralFoil result to a `POLAR_VECTORS` table.
 """
-function write_polar(filepath::String, result::NeuralFoilResult)
-    values = [result.CD zero(result.CD) result.CL result.CM]
-    return write_node_rows(filepath, deg2rad.(result.alpha), nothing, values;
-                           columns=["Cd", "Cs", "Cl", "Cm"])
-end
+write_polar(filepath::String, result::NeuralFoilResult) =
+    write_polar(filepath, deg2rad.(result.alpha), result.CL, result.CD, result.CM)
 
 """
     generate_polar_from_coordinates(x, y, output_path; Re, alpha_range=-180:1:180,
@@ -69,34 +64,43 @@ function generate_polar_from_coordinates(x::Vector, y::Vector, output_path::Stri
 end
 
 """
-    resolve_airfoil(type, info, out_dir, id; Re, alpha_range, table_format=:csv)
-        -> (new_type, new_info)
+    resolve_airfoil(type, info, out_dir, id; Re, alpha_range, table_format=:csv,
+                    ml_models_dir=nothing) -> (new_type, new_info)
 
 Resolve one awesIO `wing_airfoils` entry to a core-loadable form. `breukels_regression`
 `(t, kappa)` → `poly` coeffs (via [`lei_poly_coeffs`](@ref)); `neuralfoil`
-`(dat_file_path, …)` → a `polars` table `{id}.{table_format}` (`:csv` or `:arrow`, via
-[`generate_polar_from_dat`](@ref)) written under `out_dir`; `polars`/`poly`/`inviscid`
-pass through. `masure_regression` is not yet supported. `info` file paths should
-already be absolute.
+`(dat_file_path, …)` and `masure_regression` `(t, eta, kappa, delta, lambda, phi)` → a
+`polars` table `{id}.{table_format}` (`:csv` or `:arrow`) written under `out_dir`, the
+latter from the model in `ml_models_dir` (see [`load_masure_model`](@ref));
+`polars`/`poly`/`inviscid` pass through. `info` file paths should already be absolute.
 """
 function resolve_airfoil(type::AbstractString, info::AbstractDict, out_dir, id;
-                         Re, alpha_range, table_format::Symbol=:csv)
+                         Re, alpha_range, table_format::Symbol=:csv,
+                         ml_models_dir=nothing)
+    polar = joinpath(out_dir, "$(id).$(table_format)")
     if type == "breukels_regression"
         cl, cd, cm = lei_poly_coeffs(Float64(info["t"]), Float64(info["kappa"]))
         return "poly", Dict{String,Any}("cl_coeffs" => cl, "cd_coeffs" => cd,
                                         "cm_coeffs" => cm)
     elseif type == "neuralfoil"
-        polar = joinpath(out_dir, "$(id).$(table_format)")
         solver = NeuralFoilSolver(
             model_size=String(get(info, "model_size", "large")),
             n_crit=Float64(get(info, "n_crit", 9.0)))
         generate_polar_from_dat(String(info["dat_file_path"]), polar; Re,
             alpha_range=collect(alpha_range), solver)
         return "polars", Dict{String,Any}("polar_file_path" => polar)
+    elseif type == "masure_regression"
+        ml_models_dir === nothing &&
+            error("masure_regression airfoil $id needs ml_models_dir.")
+        absent = [name for name in MASURE_PARAMETERS if !haskey(info, name)]
+        isempty(absent) || error("masure_regression airfoil $id lacks " *
+            "$(join(absent, ", ")); it takes $(join(MASURE_PARAMETERS, ", ")).")
+        alpha = collect(Float64, alpha_range)
+        cl, cd, cm = masure_aero(load_masure_model(Re, ml_models_dir), info, alpha)
+        write_polar(polar, deg2rad.(alpha), cl, cd, cm)
+        return "polars", Dict{String,Any}("polar_file_path" => polar)
     elseif type in ("polars", "poly", "inviscid")
         return String(type), info
-    elseif type == "masure_regression"
-        error("masure_regression airfoils are not yet supported by AirfoilAero.")
     else
         error("Unknown airfoil type: $type")
     end

src/airfoil_aero/airfoil_solvers/common.jl

@@ -31,22 +31,29 @@ struct SectionSolution
     cf::Vector{Float64}
 end
 
+"""
+    write_polar(filepath, alpha, cl, cd, cm)
+
+Write a polar to a `POLAR_VECTORS` table (`alpha, Cd, Cs, Cl, Cm`), CSV or Arrow as the
+suffix of `filepath` says. `alpha` is in radians and is written in degrees.
+"""
+function write_polar(filepath::String, alpha, cl, cd, cm)
+    return write_node_rows(filepath, alpha, nothing, [cd zero(cd) cl cm];
+                           columns=["Cd", "Cs", "Cl", "Cm"])
+end
+
 """
     write_polar(filepath, sols::Vector{SectionSolution})
 
 Write a solver sweep (from any [`AbstractAirfoilSolver`](@ref)) to a `POLAR_VECTORS`
-table (`alpha, Cd, Cs, Cl, Cm`; alpha in degrees), CSV or Arrow as the suffix of
-`filepath` says. Non-converged angles (`NaN`) are skipped, so this works for both
-NeuralFoil and XFoil sweeps.
+table, skipping non-converged angles (`NaN`), so this works for both NeuralFoil and
+XFoil sweeps.
 """
 function write_polar(filepath::String, sols::Vector{SectionSolution})
     converged = filter(sol -> !isnan(sol.cl), sols)
-    alpha = [sol.alpha for sol in converged]
-    cd = [sol.cd for sol in converged]
-    cl = [sol.cl for sol in converged]
-    cm = [sol.cm for sol in converged]
-    return write_node_rows(filepath, alpha, nothing, [cd zero(cd) cl cm];
-                           columns=["Cd", "Cs", "Cl", "Cm"])
+    return write_polar(filepath, [sol.alpha for sol in converged],
+                       [sol.cl for sol in converged], [sol.cd for sol in converged],
+                       [sol.cm for sol in converged])
 end
 
 """