Parameters.h

#pragma once
#include <wmtk/OptimizerParameters.h>
#include <nlohmann/json.hpp>
#include <wmtk/Types.hpp>
#include <wmtk/components/simwild/expression_parser/Expression.hpp>
 
using ExpressionPtr = wmtk::components::simwild::expression_parser::ExpressionPtr;
 
namespace wmtk::components::topological_offset {
 
struct Parameters : public wmtk::OptimizerParameters
{
    ExpressionPtr offset_selection;
    std::set<std::string> offset_output_tag;
    std::set<std::string> protected_tags;
    bool offset_in;
    bool offset_out;
    double target_distance;
    double target_distance_rel;
    // Turn a non-converged run into a hard error instead of a warning. Off by default: a run that
    // misses the target is still a usable offset, and the warnings name the criterion that failed.
    // Integration tests set it true so a convergence regression fails rather than warns.
    bool throw_on_nonconvergence;
    bool optimize_offset;
    // Half-width of the envelope that contains every tag-region boundary during optimization.
    // Absolute; if < 0, computed from envelope_size_rel (relative to the bbox diagonal).
    double envelope_size;
    double envelope_size_rel;
 
    // ---- the smooth offset potential ----
    // Support radius of the potential, as a multiple of target_distance. Must be > 1: the offset
    // level set has to lie strictly inside the support, or the vertices on it get no gradient. A
    // band vertex that travels past the support is a hard error, not a silently frozen vertex.
    double offset_dhat_factor;
    double debug_manual_dhat;
    std::string offset_field; 
    // ---- the two convergence epsilons (was the single front_conv_rel) ----
    // Both are ABSOLUTE LENGTHS, resolved in init(); if < 0 each is computed from its _rel twin,
    // which is a fraction of the BOUNDING BOX DIAGONAL -- the same absolute/relative pair as
    // envelope_size / envelope_size_rel, and deliberately NOT a fraction of target_distance any
    // more, so changing the offset distance no longer silently changes the accuracy. init()
    // refuses either one above target_distance: an epsilon coarser than the offset it measures
    // cannot decide anything.
    //
    // THE ONE BAR. 3D measures a single quantity everywhere -- over a face's stencil, the RMS of
    // the distance to the level set along the field (OffsetPotential::relative_residual(), which
    // for the euclidean field is the relative error (Phi - c)/c) -- and compares it against
    // this. A vertex is placed when that same measure at the vertex alone is within it, which is
    // the order-0 stencil, so the vertex measure and the face measure are one measure at two
    // sample counts. Replaces vertex_conv / sag_conv, which split the two apart 2026-09-23. The
    // loop exits on the face measure alone (every face within the bar, nothing unmeasurable);
    // the vertex measure is a diagnostic since 2026-09-25. 2D exits on its chord measure, the
    // midpoint sag, which unlike the 3D stencil does not sample the chord's two ends.
    double front_conv;
    double front_conv_rel;
 
    double front_conv_frac() const { return front_conv / std::max(target_distance, 1e-16); }
    // 2D ONLY since the 3D criteria were unified. Which convergence test gates a 2D run's vertex
    // placement. F is the vertex's front objective, g its gradient, H its Gauss-Newton Hessian,
    // n its move direction; all four compare against front_conv. See front_vertex_conv_ratio().
    //   "step_size_rel" (the default): the remaining 1-D Newton step, |n.g| / (n^T H n), against
    //     front_conv.
    //   "decrement": the Newton decrement, half of (n.g)^2 / (n^T H n), against
    //     front_conv_frac() x F.
    //   "gradient_norm_rel": |n.g| against front_conv_frac() x the reference gradient, measured
    //     once on the band as constructed.
    //   "residual_error": not a stationarity measure at all -- the field's own residual at the
    //     vertex as a length (OffsetPotential::residual_length(), so |d - target_distance| for
    //     the euclidean field and the distance to the level set along the field for the smooth
    //     one), against front_conv.
    //     No objective is built and n does not enter.
    //
    // 3D DOES NOT READ THIS. Its one measure is the stencil RMS of the relative error, which is
    // not a stationarity test and has no variants -- see TopoOffsetTetMesh::face_conv_ratio().
    std::string front_conv_criterion;
    // The front is placed by a one-dimensional solve along its field normal
    // n = grad Phi / |grad Phi| -- same objective, solver and accept test, restricted to the line
    // x0 + s n -- instead of a free solve. Where a vertex sits along the front carries no offset
    // information, and in the free solve that tangential motion made fronts slide and fold where
    // two of them meet.
    bool front_normal_projection = true;
    bool front_alignment_energy = true; 
    bool sizing_collapse_min = true;
    bool deform_others = true;
    int max_rounds = 40;
    // Sampling density of THE measure -- in 3D both the criterion's and the energy's -- and of
    // the residual diagnostics that share the lattice. Order k puts these points on a face:
    //   0 -> 3, the CORNERS alone, so the measure is exactly the three vertices' placement error;
    //   k >= 1 -> the vertices of the triangle subdivided k-1 times by 4-way midpoint refinement
    //   plus the centroid of each of its 4^(k-1) sub-triangles, i.e. 4, 10, 31, 109, ...
    // The corners are IN the stencil, unlike the strictly interior lattice this replaces, because
    // the quantity measured is a distance to the level set rather than an interpolation error and
    // so is not identically zero there. Raising it costs a Phi value and gradient per sample, in
    // the ops guard's hot path as well as in energy_criterion() and every smoothing solve.
    // 2D reads this key for its diagnostics only; its chord test is still the MIDPOINT.
    // See TopoOffsetTetMesh::for_each_face_sample, TopoOffsetTriMesh::offset_edge_samples.
    int stencil_order;
    std::string front_measure;
    bool sorted_marching;
    bool sphere_trace_initialization;
    double sphere_trace_target_rel_tol; 
    bool experimental_consistent_construction_split = true;
    bool experimental_aggresive_refine = true;
    bool experimental_nonoverlapping_gates = false;
    std::string output_path; // no extension
    bool save_vtu;
 
    // Samples per side of the grid the smooth offset potential is written on, beside the result,
    // for the viewer. 0 disables it. The whole domain in 2D, one plane through the box in 3D.
    int phi_grid_resolution;
 
    int num_threads; // number of threads for parallel execution (smoothing, collapse). 0 = serial
    int max_iterations;
    // The operation passes' offset envelope half-width: the leash the front is kept inside while
    // the operation passes run, the same tube every turn, rebuilt after every smoothing pass; see
    // rebuild_offset_envelope(). Absolute-or-relative exactly as envelope_size / envelope_size_rel
    // and against the same reference, the BOUNDING BOX DIAGONAL -- it is a distance in space, and
    // tying it to target_distance made every change of the offset distance a silent change of the
    // leash as well. Also feeds the derived sizing floor (min_edge_length_rel < 0).
    double offset_envelope; 
    double offset_envelope_rel;
 
    // l_min from the paper: the shortest edge the sizing field may ask for, given as a multiple of
    // target_distance rather than of the bounding box because that is the scale the offset has;
    // min_edge_length is derived from min_edge_length_rel in init() when negative. A floor on
    // refinement, so raising it makes the result coarser. When not given, it falls back to the
    // offset envelope eps, following TetWild: a surface pinned only to within eps cannot buy
    // fidelity from shorter edges, so this is a runaway rail, not a resolution setting.
    double min_edge_length;
    double min_edge_length_rel;
 
    // ---- sizing field ----
    // bounds for VertexAttributes::m_sizing_scalar
    double min_sizing_scalar;
    double max_sizing_scalar;
    // gradation cap: neighboring vertices' sizing scalars may differ by at most this factor,
    // enforced by propagating the refinement outward (monotone, only ever lowers a
    // neighbor's scalar). <= 1 disables gradation entirely. Only used by "ring" mode.
    double sizing_gradation;
    std::string sizing_gradation_mode;
    bool adaptive_smoothing;
    int adaptive_smoothing_max_passes; 
    double adaptive_smoothing_stall_rel; 
    double adaptive_smoothing_step_rel; 
    bool pre_smooth;
 
    VectorXd box_min;
    VectorXd box_max;
 
    Parameters() = default;
 
    Parameters(const nlohmann::json& json_params)
    {
        for (const std::string& tag : json_params["offset_output_tags"]) {
            if (tag == "ambient") {
                logger().warn(
                    "'ambient' tag cannot be given explicitly to offset_output_tags, ignoring. To "
                    "set offset to 'ambient', pass offset_output_tags=[].");
                continue;
            }
            offset_output_tag.insert(tag);
        }
        for (const std::string& tag : json_params["protected_tags"]) {
            if (tag == "ambient") {
                logger().warn("'ambient' tag cannot be protected, ignoring.");
                continue;
            }
            protected_tags.insert(tag);
        }
        offset_in = json_params["offset_in"];
        offset_out = json_params["offset_out"];
        target_distance = json_params["target_distance"];
        target_distance_rel = json_params["target_distance_rel"];
        throw_on_nonconvergence = json_params["throw_on_nonconvergence"];
        optimize_offset = json_params["optimize_offset"];
        envelope_size = json_params["envelope_size"];
        envelope_size_rel = json_params["envelope_size_rel"];
        offset_dhat_factor = json_params["offset_dhat_factor"];
        debug_manual_dhat = json_params["DEBUG_manual_dhat"];
        offset_field = json_params["offset_field"];
        front_conv = json_params["front_conv"];
        front_conv_rel = json_params["front_conv_rel"];
        front_conv_criterion = json_params["front_conv_criterion"];
        stencil_order = json_params["stencil_order"];
        front_measure = json_params["front_measure"];
 
        sorted_marching = json_params["sorted_marching"];
        sphere_trace_initialization = json_params["sphere_trace_initialization"];
        sphere_trace_target_rel_tol = json_params["sphere_trace_target_rel_tol"];
        experimental_consistent_construction_split =
            json_params["EXPERIMENTAL_consistent_construction_split"];
        experimental_aggresive_refine = json_params["EXPERIMENTAL_aggresive_refine"];
        experimental_nonoverlapping_gates = json_params["EXPERIMENTAL_nonoverlapping_gates"];
        output_path = json_params["output"];
        save_vtu = json_params["save_vtu"];
        phi_grid_resolution = json_params["phi_grid_resolution"];
 
        num_threads = json_params["num_threads"];
        max_iterations = json_params["max_iterations"];
        offset_envelope = json_params["offset_envelope"];
        offset_envelope_rel = json_params["offset_envelope_rel"];
 
        min_edge_length = json_params["min_edge_length"];
        min_edge_length_rel = json_params["min_edge_length_rel"];
 
        min_sizing_scalar = json_params["min_sizing_scalar"];
        max_sizing_scalar = json_params["max_sizing_scalar"];
        sizing_gradation = json_params["sizing_gradation"];
        sizing_gradation_mode = json_params["sizing_gradation_mode"];
        adaptive_smoothing = json_params["adaptive_smoothing"];
        adaptive_smoothing_max_passes = json_params["adaptive_smoothing_max_passes"];
        adaptive_smoothing_stall_rel = json_params["adaptive_smoothing_stall_rel"];
        adaptive_smoothing_step_rel = json_params["adaptive_smoothing_step_rel"];
        pre_smooth = json_params["pre_smooth"];
 
        // ---- inherited from wmtk::OptimizerParameters ----
        debug_output = json_params["DEBUG_output"];
        lr = json_params["length_rel"];
        l = json_params["length"];
        stop_energy = json_params["stop_energy"];
        num_smoothing_passes = json_params["num_smoothing_passes"];
        interleaved_smoothing = json_params["interleaved_smoothing"];
        interleaved_smoothing_passes = json_params["interleaved_smoothing_passes"];
        split_high_valence_threshold = json_params["split_high_valence_threshold"];
        // skip_good_regions is deliberately not exposed: it would restrict a smoothing pass to
        // cells still far from stop_energy, but the smoother is what places the offset boundary,
        // so a well-shaped yet badly-placed patch is exactly what must not be skipped.
        // Every key of the coarsening group is copied here, not just the on/off switch: declaring
        // a key in the spec only makes jse inject its default into the json, so a key nothing
        // copies into this struct silently keeps whatever OptimizerParameters holds.
        coarsen_pass = json_params["coarsen_pass"];
        coarsen_unbounded = json_params["coarsen_unbounded"];
        coarsen_local_smoothing_passes = json_params["coarsen_local_smoothing_passes"];
        coarsen_smooth_ring = json_params["coarsen_smooth_ring"];
        coarsen_global_smoothing_passes = json_params["coarsen_global_smoothing_passes"];
        coarsen_max_rounds = json_params["coarsen_max_rounds"];
        stuck_refine_stall_eps = json_params["stuck_refine_stall_eps"];
        stuck_refine_cooldown = json_params["stuck_refine_cooldown"];
        stuck_refine_num_worst = json_params["stuck_refine_num_worst"];
        stuck_refine_rings = json_params["stuck_refine_rings"];
        stuck_refine_factor = json_params["stuck_refine_factor"];
        stuck_refine_min_scalar = json_params["stuck_refine_min_scalar"];
        stuck_refine_gradation = json_params["stuck_refine_gradation"];
        stuck_refine_force_split = json_params["stuck_refine_force_split"];
        front_normal_projection = json_params["front_normal_projection"];
        front_alignment_energy = json_params["front_alignment_energy"];
        sizing_collapse_min = json_params["sizing_collapse_min"];
        deform_others = json_params["deform_others"];
        max_rounds = json_params["max_rounds"];
        w_amips = json_params["w_amips"];
        smoothing_mode = json_params["smoothing_mode"];
        project_line_search_steps = json_params["project_line_search_steps"];
        project_line_search_nested_steps = json_params["project_line_search_nested_steps"];
        smooth_quality_veto = json_params["smooth_quality_veto"];
        w_envelope = 1. - w_amips;
        perform_sanity_checks = json_params["perform_sanity_checks"];
    }
 
    void init(const VectorXd& min_, const VectorXd& max_)
    {
        box_min = min_;
        box_max = max_;
 
        // Not a user knob: a topological offset preserves the topology of the region it wraps, so
        // the shared collapse must always apply the substructure link condition, or a collapse
        // across a thin band pinches the two sides together and the region stops being manifold.
        // Set here because tetwild and simwild leave the flag off.
        preserve_topology = true;
 
        // Fills diag_l, l/lr and splitting_l2 / collapsing_l2. It also derives eps from epsr,
        // which the offset never reads: its envelope tolerance is m_envelope_eps, set on the mesh.
        init_lengths_from_diagonal((max_ - min_).norm());
 
        // The split gate moved out to twice the collapse gate, so no half a split produces is a
        // collapse candidate: splitting_l2 = (8/5 l)^2 against collapsing_l2 = (4/5 l)^2, which
        // is left where it is. The overlapping TetWild gates (4/3 and 4/5), still the default, let
        // the passes trade the same edges every turn -- see the declaration for the measurement.
        if (experimental_nonoverlapping_gates) {
            splitting_l2 = l * l * (64 / 25.);
        }
 
        if (target_distance > 0) {
            target_distance_rel = target_distance / diag_l;
        } else {
            target_distance = target_distance_rel * diag_l;
        }
 
        // An ordinary relative length: it bounds how far a region boundary may drift in space, so
        // the bounding box is the right reference.
        if (envelope_size > 0) {
            envelope_size_rel = envelope_size / diag_l;
        } else {
            envelope_size = envelope_size_rel * diag_l;
        }
 
        // The convergence epsilon, the same absolute-or-relative pair as the envelope and
        // against the same reference. It is a length in space, so the bounding box diagonal is
        // the reference, not target_distance: tying the accuracy to the offset distance made
        // every change of target_distance a silent change of accuracy as well.
        if (front_conv > 0) {
            front_conv_rel = front_conv / diag_l;
        } else {
            front_conv = front_conv_rel * diag_l;
        }
 
        // An epsilon coarser than the offset it measures decides nothing: every front vertex is
        // "placed" and every face "resolved" from the first turn, whatever the offset looks like.
        // Checked on the resolved ABSOLUTE values, so it catches the mistake whichever of the two
        // forms the config used to state it.
        if (front_conv > target_distance) {
            log_and_throw_error(
                "front_conv {} must be <= target_distance {}: the convergence epsilon cannot be "
                "coarser than the offset distance it measures, or every front face reads as "
                "resolved from the first turn",
                front_conv,
                target_distance);
        }
 
        // The operation leash, the same absolute-or-relative pair as the envelope and the
        // convergence epsilon, against the same reference.
        if (offset_envelope > 0) {
            offset_envelope_rel = offset_envelope / diag_l;
        } else {
            offset_envelope = offset_envelope_rel * diag_l;
        }
 
        // l_min is relative to the offset distance rather than the bounding box: it is the offset
        // that has to be resolved. See the declaration.
        if (min_edge_length_rel < 0) {
            // The envelope eps expressed as a multiple of target_distance, which is what this
            // wants. offset_envelope_rel is a fraction of the BBOX DIAGONAL since 2026-09-24, so
            // the conversion goes through the resolved absolute rather than being the identity it
            // used to be -- the derived floor is unchanged in model units either way.
            min_edge_length_rel =
                std::max(offset_envelope / std::max(target_distance, 1e-16), 1e-12);
        }
        if (min_edge_length < 0) {
            min_edge_length = min_edge_length_rel * target_distance;
        } else {
            min_edge_length_rel = min_edge_length / std::max(target_distance, 1e-16);
        }
    }
};
} // namespace wmtk::components::topological_offset