// velocity_curve.cpp — see velocity_curve.h. Pure eval + editing/clamp/inverse map; no host types.

#include "core/instrument/engine/velocity_curve.h"

#include <algorithm>
#include <cmath>
#include <utility>

namespace reasampler::instrument::engine {

namespace {

double clampVelocity(double v) { return std::clamp(v, kVelMin, kVelMax); }
double clampAmp(double a) { return std::clamp(a, kAmpMin, kAmpMax); }

// X spans the width for [0,127]; Y spans (height-1) rows for amp [0,1] with amp 1 at the TOP
// (pixel y increases downward, so this axis is inverted relative to amp).
double velPerPixel(const VelocityCurve::Box& box) {
    const int w = std::max(0, box.width);
    if (w <= 0) return 0.0;
    return (kVelMax - kVelMin) / static_cast<double>(w);
}
double ampPerPixel(const VelocityCurve::Box& box) {
    const int h = std::max(0, box.height);
    if (h <= 1) return 0.0;
    return (kAmpMax - kAmpMin) / static_cast<double>(h - 1);
}
int velToX(const VelocityCurve::Box& box, double velocity) {
    const int w = std::max(0, box.width);
    if (w <= 0) return box.left;
    const double frac = (clampVelocity(velocity) - kVelMin) / (kVelMax - kVelMin);
    return box.left + static_cast<int>(frac * static_cast<double>(w) + 0.5);
}
int ampToY(const VelocityCurve::Box& box, double amp) {
    const int h = std::max(0, box.height);
    if (h <= 1) return box.top;
    const double frac = (clampAmp(amp) - kAmpMin) / (kAmpMax - kAmpMin);
    return box.top + static_cast<int>((1.0 - frac) * static_cast<double>(h - 1) + 0.5);
}

}  // namespace

VelocityCurve VelocityCurve::flat() {
    VelocityCurve c;
    c.points_ = {{kVelMin, kAmpMax}, {kVelMax, kAmpMax}};
    return c;
}

VelocityCurve VelocityCurve::linear() {
    VelocityCurve c;
    c.points_ = {{kVelMin, kAmpMin}, {kVelMax, kAmpMax}};
    return c;
}

VelocityCurve VelocityCurve::fromPoints(std::vector<VelocityPoint> pts) {
    // Stable sort so coincident-X points keep their wire order (eval stays well-defined for
    // duplicate-X knots).
    for (VelocityPoint& p : pts) {
        p.velocity = clampVelocity(p.velocity);
        p.amp = clampAmp(p.amp);
    }
    std::stable_sort(pts.begin(), pts.end(),
                     [](const VelocityPoint& a, const VelocityPoint& b) {
                         return a.velocity < b.velocity;
                     });
    if (pts.size() < 2) return flat();
    if (pts.front().velocity > kVelMin) {
        pts.insert(pts.begin(), VelocityPoint{kVelMin, pts.front().amp});
    } else {
        pts.front().velocity = kVelMin;
    }
    if (pts.back().velocity < kVelMax) {
        pts.push_back(VelocityPoint{kVelMax, pts.back().amp});
    } else {
        pts.back().velocity = kVelMax;
    }
    VelocityCurve c;
    c.points_ = std::move(pts);
    return c;
}

namespace {

// Fritsch-Carlson monotone-cubic tangent: a sign change (or flat) neighbour is a local extremum,
// so the tangent pins to 0 to avoid overshoot; otherwise the weighted-harmonic-mean tangent,
// which for collinear knots (dPrev==dNext) reduces exactly to the shared secant — this is what
// makes the spline reproduce a straight line to ~1e-15 for linear()-style input.
double fritschCarlsonTangent(double dPrev, double dNext, double spanPrev, double spanNext) {
    if (dPrev * dNext <= 0.0) return 0.0;
    const double w1 = 2.0 * spanNext + spanPrev;
    const double w2 = spanNext + 2.0 * spanPrev;
    return (w1 + w2) / (w1 / dPrev + w2 / dNext);
}

}  // namespace

double VelocityCurve::eval(double velocity) const {
    if (points_.empty()) return kAmpMax;
    if (points_.size() == 1) return clampAmp(points_[0].amp);
    const double v = clampVelocity(velocity);
    if (v <= points_.front().velocity) return clampAmp(points_.front().amp);
    if (v >= points_.back().velocity) return clampAmp(points_.back().amp);
    for (std::size_t i = 0; i + 1 < points_.size(); ++i) {
        const VelocityPoint& a = points_[i];
        const VelocityPoint& b = points_[i + 1];
        if (v >= a.velocity && v <= b.velocity) {
            const double span = b.velocity - a.velocity;
            // Coincident-X neighbours (a step): zero-width segment, no interior to blend.
            if (span <= 0.0) return clampAmp(b.amp);

            // Monotone cubic Hermite (Fritsch-Carlson): provably stays within [a.amp, b.amp]
            // between the two knots (no overshoot), reproducing a straight line for collinear input.
            const double d = (b.amp - a.amp) / span;

            double mA = d;
            if (i > 0) {
                const VelocityPoint& prev = points_[i - 1];
                const double spanPrev = a.velocity - prev.velocity;
                if (spanPrev > 0.0) {
                    const double dPrev = (a.amp - prev.amp) / spanPrev;
                    mA = fritschCarlsonTangent(dPrev, d, spanPrev, span);
                } else {
                    mA = 0.0;
                }
            }
            double mB = d;
            if (i + 2 < points_.size()) {
                const VelocityPoint& next = points_[i + 2];
                const double spanNext = next.velocity - b.velocity;
                if (spanNext > 0.0) {
                    const double dNext = (next.amp - b.amp) / spanNext;
                    mB = fritschCarlsonTangent(d, dNext, span, spanNext);
                } else {
                    mB = 0.0;
                }
            }

            const double t = (v - a.velocity) / span;
            const double t2 = t * t;
            const double t3 = t2 * t;
            const double h00 = 2.0 * t3 - 3.0 * t2 + 1.0;
            const double h10 = t3 - 2.0 * t2 + t;
            const double h01 = -2.0 * t3 + 3.0 * t2;
            const double h11 = t3 - t2;
            const double y = h00 * a.amp + h10 * span * mA + h01 * b.amp + h11 * span * mB;
            return clampAmp(y);
        }
    }
    return clampAmp(points_.back().amp);  // unreachable (v is between the endpoints)
}

std::size_t VelocityCurve::addPoint(double velocity, double amp) {
    const VelocityPoint p{clampVelocity(velocity), clampAmp(amp)};
    // First index strictly greater, so a duplicate-X point lands immediately after the existing one.
    std::size_t i = 0;
    while (i < points_.size() && points_[i].velocity <= p.velocity) ++i;
    points_.insert(points_.begin() + static_cast<std::ptrdiff_t>(i), p);
    return i;
}

VelocityPoint VelocityCurve::movePoint(std::size_t index, double velocity, double amp) {
    if (index >= points_.size()) return VelocityPoint{};  // no-op (out of range)
    const bool isFirst = (index == 0);
    const bool isLast  = (index + 1 == points_.size());

    double newAmp = clampAmp(amp);
    double newVel;
    if (isFirst) {
        newVel = kVelMin;
    } else if (isLast) {
        newVel = kVelMax;
    } else {
        const double lo = points_[index - 1].velocity;
        const double hi = points_[index + 1].velocity;
        newVel = std::clamp(clampVelocity(velocity), lo, hi);
    }
    points_[index] = VelocityPoint{newVel, newAmp};
    return points_[index];
}

bool VelocityCurve::deletePoint(std::size_t index) {
    if (index >= points_.size()) return false;
    if (index == 0 || index + 1 == points_.size()) return false;  // endpoints are not deletable
    points_.erase(points_.begin() + static_cast<std::ptrdiff_t>(index));
    return true;
}

VelocityCurve::CurvePixel VelocityCurve::pixelFromPoint(const Box& box, const VelocityPoint& p) {
    return CurvePixel{velToX(box, p.velocity), ampToY(box, p.amp)};
}

VelocityPoint VelocityCurve::pointFromPixel(const Box& box, int x, int y) {
    // Exact inverse of velToX/ampToY (within one pixel); degenerate dims collapse the same way.
    VelocityPoint p;
    const int w = std::max(0, box.width);
    const int h = std::max(0, box.height);
    p.velocity = (w <= 0)
        ? kVelMin
        : clampVelocity(kVelMin + static_cast<double>(x - box.left) / static_cast<double>(w) *
                                      (kVelMax - kVelMin));
    p.amp = (h <= 1)
        ? kAmpMax
        : clampAmp(kAmpMax - static_cast<double>(y - box.top) / static_cast<double>(h - 1) *
                                 (kAmpMax - kAmpMin));
    return p;
}

int VelocityCurve::pointAtPixel(const Box& box, int x, int y) const {
    for (std::size_t i = 0; i < points_.size(); ++i) {
        const int px = velToX(box, points_[i].velocity);
        const int py = ampToY(box, points_[i].amp);
        if (std::abs(x - px) <= kCurveNodeGrabRadius && std::abs(y - py) <= kCurveNodeGrabRadius) {
            return static_cast<int>(i);
        }
    }
    return -1;
}

VelocityCurve VelocityCurve::resolvePointDrag(const VelocityCurve& grabCurve, std::size_t index,
                                              const Box& box, int dxPixels, int dyPixels) {
    VelocityCurve out = grabCurve;
    if (index >= out.points_.size()) return out;  // out of range -> no motion
    const double velPerPx = velPerPixel(box);
    const double ampPerPx = ampPerPixel(box);
    if (velPerPx <= 0.0 || ampPerPx <= 0.0) return out;  // degenerate box -> no motion

    const VelocityPoint& grab = grabCurve.points_[index];
    const double newVel = grab.velocity + static_cast<double>(dxPixels) * velPerPx;
    // Y increases downward but amp increases upward, so a downward drag (positive dy) LOWERS amp.
    const double newAmp = grab.amp - static_cast<double>(dyPixels) * ampPerPx;
    out.movePoint(index, newVel, newAmp);  // applies box + neighbour-X + endpoint-pin clamps
    return out;
}

bool VelocityCurve::equals(const VelocityCurve& other, double eps) const {
    if (points_.size() != other.points_.size()) return false;
    for (std::size_t i = 0; i < points_.size(); ++i) {
        if (std::fabs(points_[i].velocity - other.points_[i].velocity) > eps) return false;
        if (std::fabs(points_[i].amp - other.points_[i].amp) > eps) return false;
    }
    return true;
}

}  // namespace reasampler::instrument::engine