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;
    }
 
    // Use last interval for x above range
    if (x_query >= tbl.x[n - 1]) {
        double dx = tbl.x[n - 1] - tbl.x[n - 2];
        return (dx > 0.0) ? (tbl.y[n - 1] - tbl.y[n - 2]) / dx : 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 */