TableData.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: Apache-2.0
//
// Copyright 2026 Caleb Buahin
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

 
#ifndef OPENSWMM_ENGINE_TABLE_DATA_HPP
#define OPENSWMM_ENGINE_TABLE_DATA_HPP
 
#include "../core/FilePathPair.hpp"
#include "../core/StringCase.hpp"
#include <string>
#include <vector>
#include <unordered_map>
#include <cstdio>
#include <cstdint>
#include <cassert>
#include <algorithm>
#include <cmath>
 
namespace openswmm {
 
// ============================================================================
// Table types (matching legacy SWMM enums.h)
// ============================================================================

enum class TableType : int {
    TIMESERIES = 0,      
    CURVE_STORAGE = 1,   
    CURVE_DIVERSION = 2, 
    CURVE_RATING = 3,    
    CURVE_SHAPE = 4,     
    CURVE_CONTROL = 5,   
    CURVE_TIDAL = 6,     
    CURVE_PUMP1 = 7,     
    CURVE_PUMP2 = 8,     
    CURVE_PUMP3 = 9,     
    CURVE_PUMP4 = 10,    
    CURVE_PUMP5 = 11,    
    CURVE_WEIR = 12      
};
 
// ============================================================================
// TableCursor
// ============================================================================

struct TableCursor {
    int index = 0;      
    int direction = +1; 

    void reset() noexcept {
        index = 0;
        direction = +1;
    }
};
 
// ============================================================================
// Table
// ============================================================================

struct TableBlock {
    std::vector<double> x;          
    std::vector<double> y;          
    int num_cols = 1;               
    std::size_t file_row_start = 0; 

    std::size_t num_rows() const noexcept { return x.size(); }
    bool empty() const noexcept { return x.empty(); }
};
 
struct Table {
    std::string id; 
    TableType type; 

    std::string comment;
    std::vector<double> x; 
    std::vector<double> y; 
    TableCursor cursor;    
 
    // ---- Time-only (relative) timeseries rows ----
    // Legacy SWMM seeds every timeseries' lastDate with StartDate+StartTime
    // (input.c:170), so rows authored WITHOUT a date are elapsed times
    // anchored at the simulation start, until an explicit date re-anchors
    // the series. The parser counts those leading date-less rows here and
    // resolve_cross_references() adds options.start_date to exactly those
    // rows (recording the applied offset in rel_anchor so re-resolution is
    // idempotent and a later START_DATE change re-anchors by the delta).
    // InpWriter subtracts rel_anchor to emit the rows back in their
    // authored time-only form; rows at index >= n_relative are absolute
    // date/times and are written with explicit dates.
    int n_relative = 0;      
    double rel_anchor = 0.0; 
 
    // ---- File-backed time series support ----
    bool is_file_based = false;                      
    std::FILE* file_handle = nullptr;                
    FilePathPair file_path;                          
    TableBlock first_boundary;                       
    TableBlock last_boundary;                        
    TableBlock cache;                                
    double dx_min = 0.0;                             
    double x_min = 0.0;                              
    double x_max = 0.0;                              
    int num_cols = 1;                                
    std::size_t total_rows = 0;                      
    std::size_t num_cache_rows = 8192;               
    long data_start_offset = 0;                      
    std::vector<long> row_offsets;                   
    std::vector<std::string> column_ids;             
    std::unordered_map<std::string, int> column_map; 
 
    static constexpr std::size_t INDEX_STRIDE = 4096; 

    std::size_t size() const noexcept { return x.size(); }

    bool empty() const noexcept { return x.empty(); }
};
 
// ============================================================================
// Cursor-optimized lookup
// ============================================================================

inline double table_lookup_cursor(Table& tbl, double x_query) noexcept {
    const int n = static_cast<int>(tbl.x.size());
    if (n == 0)
        return 0.0;
    if (n == 1)
        return tbl.y[0];
 
    // Clamp below first entry
    if (x_query <= tbl.x[0]) {
        tbl.cursor.index = 0;
        tbl.cursor.direction = +1;
        return tbl.y[0];
    }
 
    // Clamp above last entry — strictly above only: legacy table_lookup
    // (src/legacy/engine/table.c:400) interpolates the last segment when
    // x == x[n-1] (its scan condition is `x <= x2`), so equality must fall
    // through to the interpolation below for bit parity.
    if (x_query > tbl.x[n - 1]) {
        tbl.cursor.index = n - 1;
        tbl.cursor.direction = -1;
        return tbl.y[n - 1];
    }
 
    // Start from cursor position and seek in the most likely direction
    int idx = std::clamp(tbl.cursor.index, 0, n - 2);
 
    // Forward seek
    while (idx < n - 1 && tbl.x[idx + 1] < x_query) {
        ++idx;
        tbl.cursor.direction = +1;
    }
 
    // Backward seek — `>=` so a query landing exactly on an interior knot
    // selects the segment ENDING at that knot, matching legacy table_lookup's
    // first-match linear scan (`if (x <= x2) return table_interpolate(...)`,
    // src/legacy/engine/table.c:426).
    while (idx > 0 && tbl.x[idx] >= x_query) {
        --idx;
        tbl.cursor.direction = -1;
    }
 
    tbl.cursor.index = idx;
 
    // PARITY: op-for-op transliteration of legacy table_interpolate
    // (src/legacy/engine/table.c:54-67): multiply THEN divide —
    //   y1 + (x - x1) * (y2 - y1) / dx
    // (the previous t=(x-x1)/dx; y1+t*(y2-y1) form rounds differently by
    // 1 ULP), and the legacy degenerate-segment guard |dx| < 1e-20 → mean.
    const double x1 = tbl.x[idx], y1 = tbl.y[idx];
    const double x2 = tbl.x[idx + 1], y2 = tbl.y[idx + 1];
    const double dx = x2 - x1;
    if (std::fabs(dx) < 1.0e-20)
        return (y1 + y2) / 2.;
    return y1 + (x_query - x1) * (y2 - y1) / dx;
}

inline double table_tseries_lookup_cursor(Table& tbl, double x_query) noexcept {
    const int n = static_cast<int>(tbl.x.size());
    if (n == 0)
        return 0.0;
 
    // Before first entry or after last entry → 0 (series not active)
    if (x_query <= tbl.x[0]) {
        tbl.cursor.index = 0;
        tbl.cursor.direction = +1;
        return (x_query < tbl.x[0]) ? 0.0 : tbl.y[0];
    }
    if (x_query >= tbl.x[n - 1]) {
        tbl.cursor.index = n - 1;
        tbl.cursor.direction = -1;
        return (x_query > tbl.x[n - 1]) ? 0.0 : tbl.y[n - 1];
    }
 
    // Delegate to cursor-based linear interpolation for the in-range case
    return table_lookup_cursor(tbl, x_query);
}

inline double table_step_cursor(Table& tbl, double x_query) noexcept {
    const int n = static_cast<int>(tbl.x.size());
    if (n == 0)
        return 0.0;
    if (n == 1)
        return tbl.y[0];
 
    // Before first entry → 0 (no rain before first recorded value)
    if (x_query < tbl.x[0]) {
        tbl.cursor.index = 0;
        tbl.cursor.direction = +1;
        return 0.0;
    }
 
    // At or past last entry → return last value.
    // The caller (Gage.cpp) handles the rain interval cutoff and returns 0
    // after entry_time + rainInterval. This allows the last entry's value to
    // be used for its full recording interval before going to zero.
    if (x_query >= tbl.x[n - 1]) {
        tbl.cursor.index = n - 1;
        tbl.cursor.direction = -1;
        return tbl.y[n - 1];
    }
 
    // Seek to the interval containing x_query: x[idx] <= x_query < x[idx+1]
    int idx = std::clamp(tbl.cursor.index, 0, n - 2);
 
    while (idx < n - 1 && tbl.x[idx + 1] <= x_query) {
        ++idx;
        tbl.cursor.direction = +1;
    }
    while (idx > 0 && tbl.x[idx] > x_query) {
        --idx;
        tbl.cursor.direction = -1;
    }
 
    tbl.cursor.index = idx;
    return tbl.y[idx];
}
 
// ============================================================================
// Storage volume by trapezoidal integration of area curve
// ============================================================================

inline double table_getStorageVolume(Table& tbl, double depth) noexcept {
    const int n = static_cast<int>(tbl.x.size());
    if (n == 0 || depth <= 0.0)
        return 0.0;
 
    double x1 = tbl.x[0];
    double a1 = tbl.y[0];
 
    // Target below first entry — triangular approximation
    if (depth <= x1) {
        if (x1 < 1.0e-6)
            return 0.0;
        return (a1 / x1) * depth * depth / 2.0;
    }
 
    // Traverse entries using end-area (trapezoidal) method
    double v = 0.0;
    double dx = 0.0, dy = 0.0;
    for (int i = 1; i < n; ++i) {
        double x2 = tbl.x[i];
        double a2 = tbl.y[i];
        if (x2 >= depth) {
            // Bracketed — interpolate area at target depth.
            // PARITY: legacy table_getStorageVolume (table.c:624) calls
            // table_interpolate (table.c:54): multiply THEN divide,
            // |dx| < 1e-20 degenerate guard → mean of end areas.
            const double dxi = x2 - x1;
            double a;
            if (std::fabs(dxi) < 1.0e-20)
                a = (a1 + a2) / 2.;
            else
                a = a1 + (depth - x1) * (a2 - a1) / dxi;
            return v + (a1 + a) / 2.0 * (depth - x1);
        }
        dx = x2 - x1;
        dy = a2 - a1;
        v += (a1 + a2) / 2.0 * dx;
        x1 = x2;
        a1 = a2;
    }
 
    // Extrapolate beyond last entry
    if (dx > 1.0e-6) {
        double s = dy / dx;
        double a = a1 + s * (depth - x1);
        if (a < 0.0) {
            v -= a1 * a1 / s / 2.0;
        } else {
            v += (a1 + a) / 2.0 * (depth - x1);
        }
    }
    return v;
}
 
// ============================================================================
// Inverse volume-to-depth for storage curves
// ============================================================================

inline double table_getStorageDepth(Table& tbl, double volume) noexcept {
    const int n = static_cast<int>(tbl.x.size());
    if (n == 0 || volume <= 0.0)
        return 0.0;
 
    double x1 = tbl.x[0];
    double a1 = tbl.y[0];
    double v_accum = 0.0;
 
    for (int i = 1; i < n; ++i) {
        double x2 = tbl.x[i];
        double a2 = tbl.y[i];
        double dv = (a1 + a2) / 2.0 * (x2 - x1);
 
        if (v_accum + dv >= volume) {
            // Solve the quadratic: target volume falls in this interval.
            // Area varies linearly: a(d) = a1 + s*(d - x1), s = (a2-a1)/(x2-x1)
            // Volume:  dv = a1*(d-x1) + s*(d-x1)^2/2
            // Rearrange to standard form and apply quadratic formula.
            double dx = x2 - x1;
            double dv_need = volume - v_accum;
            if (dx < 1.0e-10)
                return x1;
 
            double s = (a2 - a1) / dx;
            double depth;
            if (std::fabs(s) < 1.0e-10) {
                // Rectangular slice: dv = a1 * dd
                depth = (a1 > 0.0) ? x1 + dv_need / a1 : x1;
            } else {
                // Quadratic: s/2 * dd^2 + a1 * dd - dv_need = 0
                double disc = a1 * a1 + 2.0 * s * dv_need;
                if (disc < 0.0)
                    disc = 0.0;
                depth = x1 + (-a1 + std::sqrt(disc)) / s;
            }
            return std::min(depth, x2);
        }
 
        v_accum += dv;
        x1 = x2;
        a1 = a2;
    }
 
    // Volume exceeds curve range — extrapolate linearly using last slope
    if (tbl.x.size() >= 2) {
        int last = n - 1;
        double dx = tbl.x[last] - tbl.x[last - 1];
        double da = tbl.y[last] - tbl.y[last - 1];
        double s = (dx > 1.0e-10) ? da / dx : 0.0;
        double dv_need = volume - v_accum;
        double a1_ext = tbl.y[last];
        double depth;
        if (std::fabs(s) < 1.0e-10) {
            depth = (a1_ext > 0.0) ? tbl.x[last] + dv_need / a1_ext : tbl.x[last];
        } else {
            double disc = a1_ext * a1_ext + 2.0 * s * dv_need;
            if (disc < 0.0)
                disc = 0.0;
            depth = tbl.x[last] + (-a1_ext + std::sqrt(disc)) / s;
        }
        return depth;
    }
    return tbl.x[n - 1];
}
 
// ============================================================================
// Generic inverse lookup (y → x)
// ============================================================================

inline double table_inverseLookup(const Table& tbl, double y_query) noexcept {
    const int n = static_cast<int>(tbl.x.size());
    if (n == 0)
        return 0.0;
    if (n == 1)
        return tbl.x[0];
 
    if (y_query <= tbl.y[0])
        return tbl.x[0];
    if (y_query >= tbl.y[n - 1])
        return tbl.x[n - 1];
 
    for (int i = 1; i < n; ++i) {
        if (tbl.y[i] >= y_query) {
            double dy = tbl.y[i] - tbl.y[i - 1];
            if (dy <= 0.0)
                return tbl.x[i - 1];
            double t = (y_query - tbl.y[i - 1]) / dy;
            return tbl.x[i - 1] + t * (tbl.x[i] - tbl.x[i - 1]);
        }
    }
    return tbl.x[n - 1];
}
 
// ============================================================================
// Extrapolating lookup (extends beyond table bounds)
// ============================================================================

inline double table_lookupEx(const Table& tbl, double x_query) noexcept {
    // PARITY: op-for-op transliteration of legacy table_lookupEx
    // (src/legacy/engine/table.c:469-505):
    //   - below first entry: x/x1*y1 when x1 > 0 (line through the origin),
    //     else y1 — NOT first-interval slope extrapolation;
    //   - in range: table_interpolate (multiply THEN divide, |dx| < 1e-20
    //     degenerate guard → mean);
    //   - above last entry: extrapolate with the slope `s` of the last
    //     segment seen during the walk (skipping x2 == x1 segments),
    //     clamped at s >= 0.
    const int n = static_cast<int>(tbl.x.size());
    if (n == 0)
        return 0.0;
 
    double x1 = tbl.x[0];
    double y1 = tbl.y[0];
    double s = 0.0;
    if (x_query <= x1) {
        if (x1 > 0.0)
            return x_query / x1 * y1;
        else
            return y1;
    }
    for (int i = 1; i < n; ++i) {
        const double x2 = tbl.x[i];
        const double y2 = tbl.y[i];
        if (x2 != x1)
            s = (y2 - y1) / (x2 - x1);
        if (x_query <= x2) {
            const double dx = x2 - x1;
            if (std::fabs(dx) < 1.0e-20)
                return (y1 + y2) / 2.;
            return y1 + (x_query - x1) * (y2 - y1) / dx;
        }
        x1 = x2;
        y1 = y2;
    }
    if (s < 0.0)
        s = 0.0;
    return y1 + s * (x_query - x1);
}
 
// ============================================================================
// Step-function interval lookup (first entry > x)
// ============================================================================

inline double table_intervalLookup(const Table& tbl, double x_query) noexcept {
    const int n = static_cast<int>(tbl.x.size());
    if (n == 0)
        return 0.0;
 
    for (int i = 0; i < n; ++i) {
        if (tbl.x[i] > x_query)
            return tbl.y[i];
    }
    return tbl.y[n - 1];
}
 
// ============================================================================
// Slope at a point
// ============================================================================

inline double table_getSlope(const Table& tbl, double x_query) noexcept {
    const int n = static_cast<int>(tbl.x.size());
    if (n < 2)
        return 0.0;
 
    // Use first interval for x below range
    if (x_query <= tbl.x[0]) {
        double dx = tbl.x[1] - tbl.x[0];
        return (dx > 0.0) ? (tbl.y[1] - tbl.y[0]) / dx : 0.0;
    }
 
    // Above the last entry the slope is 0: legacy table_getSlope's scan
    // runs off the end with x1 = x2 = the last entry (dx == 0 → 0.0). A
    // head-based pump driven past its curve therefore contributes NO dqdh
    // to its end nodes' surcharge denominator (session81-small-pump-fm-
    // model's 5002 at 142 ft of head: the last segment's -15.6 GPM/ft
    // here kept node 16's Picard update damped where legacy let it run).
    // Exactly AT the last entry legacy breaks on `x <= x2` with the last
    // segment bracketed, so that case falls through to the loop.
    if (x_query > tbl.x[n - 1])
        return 0.0;
 
    for (int i = 1; i < n; ++i) {
        if (tbl.x[i] >= x_query) {
            double dx = tbl.x[i] - tbl.x[i - 1];
            return (dx > 0.0) ? (tbl.y[i] - tbl.y[i - 1]) / dx : 0.0;
        }
    }
    return 0.0;
}
 
// ============================================================================
// Maximum y in non-decreasing portion
// ============================================================================

inline double table_getMaxY(const Table& tbl) noexcept {
    const int n = static_cast<int>(tbl.x.size());
    if (n == 0)
        return 0.0;
    double ymax = tbl.y[0];
    for (int i = 1; i < n; ++i) {
        if (tbl.y[i] < ymax)
            break;
        ymax = tbl.y[i];
    }
    return ymax;
}
 
// ============================================================================
// TableData — SoA collection of all tables
// ============================================================================

struct TableData {
    std::vector<Table> tables; 

    std::unordered_map<std::string, std::vector<int>, CiHash, CiEqual> by_name;

    std::size_t count() const noexcept { return tables.size(); }

    Table& operator[](int idx) { return tables[static_cast<std::size_t>(idx)]; }
    const Table& operator[](int idx) const { return tables[static_cast<std::size_t>(idx)]; }

    int add(const std::string& id, TableType type) {
        Table t;
        t.id = id;
        t.type = type;
        tables.push_back(std::move(t));
        const int idx = static_cast<int>(tables.size()) - 1;
        by_name[id].push_back(idx);
        return idx;
    }

    void rebuild_index() {
        by_name.clear();
        by_name.reserve(tables.size());
        for (std::size_t i = 0; i < tables.size(); ++i)
            by_name[tables[i].id].push_back(static_cast<int>(i));
    }

    int find_by_kind(std::string_view name, bool want_timeseries) const noexcept {
        const auto it = by_name.find(name);
        if (it == by_name.end())
            return -1;
        for (const int i : it->second) {
            const bool is_ts = tables[static_cast<std::size_t>(i)].type == TableType::TIMESERIES;
            if (is_ts == want_timeseries)
                return i;
        }
        return -1;
    }

    void reset_cursors() noexcept {
        for (auto& t : tables)
            t.cursor.reset();
    }
};
 
// ============================================================================
// Table validation
// ============================================================================
 
struct TableValidation {
    bool valid = true;
    std::vector<std::string> errors;
    std::vector<std::string> warnings;
};

TableValidation validate_table(Table& tbl);
 
// ============================================================================
// File-backed table API
// ============================================================================

bool table_open_file(Table& tbl, std::size_t boundary_rows = 128);
std::size_t table_load_cache(Table& tbl, std::size_t start_row);
 
// ============================================================================
// Multicolumn lookup API
// ============================================================================

double table_lookup_column(Table& tbl, int col_idx, double x_query);
double table_step_column(Table& tbl, int col_idx, double x_query);
 
} /* namespace openswmm */
 
#endif /* OPENSWMM_ENGINE_TABLE_DATA_HPP */