diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..c210347 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,36 @@ +name: CI + +on: + push: + branches: ["**"] + pull_request: + branches: ["**"] + +jobs: + build-and-test: + runs-on: ubuntu-24.04 + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Install build tools + run: | + sudo apt-get update -qq + sudo apt-get install -y -qq gcc g++ cmake make + + - name: Configure (CMake + FetchContent GTest) + run: cmake -B build -DCMAKE_BUILD_TYPE=Release + + - name: Build + run: cmake --build build --parallel "$(nproc)" + + - name: Run tests + run: build/tests/test_polyfit --gtest_output=xml:results.xml + + - name: Publish test results + if: always() + uses: mikepenz/action-junit-report@v4 + with: + report_paths: results.xml + check_name: Google Test Results diff --git a/.gitignore b/.gitignore index c6127b3..3e6adbc 100644 --- a/.gitignore +++ b/.gitignore @@ -50,3 +50,4 @@ modules.order Module.symvers Mkfile.old dkms.conf +build/ diff --git a/CMakeLists.txt b/CMakeLists.txt index 5ce154a..40875ac 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,16 +1,13 @@ -cmake_minimum_required(VERSION 3.10) -project(polyfit_demo C) +cmake_minimum_required(VERSION 3.14) +project(polyfit C CXX) set(CMAKE_C_STANDARD 11) +set(CMAKE_CXX_STANDARD 17) -# Source files -set(SOURCES demo.c polyfit.c) +# Build polyfit as a static library so both the demo and tests can link it +add_library(polyfit STATIC polyfit.c) +target_include_directories(polyfit PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) +target_link_libraries(polyfit PUBLIC m) -# Executable name -set(TARGET demo) - -# Add the executable -add_executable(${TARGET} ${SOURCES}) - -# Link the math library -target_link_libraries(${TARGET} m) +enable_testing() +add_subdirectory(tests) diff --git a/polyfit.c b/polyfit.c index 5a7e031..066a958 100644 --- a/polyfit.c +++ b/polyfit.c @@ -14,6 +14,7 @@ */ #include "polyfit.h" +#include /*============================================================================*/ /* PRIVATE FUNCTION DECLARATIONS */ @@ -26,16 +27,6 @@ static polyfit_error_t allocate_matrix(float*** matrix, int32_t rows, static void free_matrix(float** matrix, int32_t rows); static polyfit_error_t validate_input_arrays(const float* x, const float* y, int32_t num_points); -static bool is_matrix_singular(float** A, int32_t n); - -/*============================================================================*/ -/* GLOBAL VARIABLES */ -/*============================================================================*/ - -static const polyfit_config_t default_config = { - .absolute_threshold = POLYFIT_ABSOLUTE_THRESHOLD, - .relative_threshold = POLYFIT_RELATIVE_THRESHOLD, - .enable_pivot_check = true}; /*============================================================================*/ /* PUBLIC FUNCTION IMPLEMENTATIONS */ @@ -124,13 +115,6 @@ polyfit_error_t polyfit_least_squares(const float* x, const float* y, } } - // Check for singular matrix - if (is_matrix_singular(A, degree + 1)) { - free_matrix(A, degree + 1); - free(B); - return POLYFIT_ERROR_SINGULAR_MATRIX; - } - // Solve the system error = gaussian_elimination(A, B, result_poly->coefficients, degree + 1); @@ -163,16 +147,6 @@ polyfit_error_t polyfit_evaluate(const Polynomial* poly, float x, *result = *result * x + poly->coefficients[i]; } - // Apply thresholding - float max_magnitude; - polyfit_get_max_coefficient_magnitude(poly, &max_magnitude); - - if (polyfit_is_nearly_zero(*result, default_config.absolute_threshold) || - polyfit_is_nearly_zero( - *result, default_config.relative_threshold * max_magnitude)) { - *result = 0.0f; - } - return POLYFIT_SUCCESS; } @@ -267,16 +241,14 @@ polyfit_error_t polyfit_eval_at(const float* x, const float* y, } // Fit the polynomial - Polynomial* poly = polyfit(x, y, num_points, degree, NULL); + polyfit_error_t error; + Polynomial* poly = polyfit(x, y, num_points, degree, &error); if (poly == NULL) { - // Determine which error occurred by trying again with error output - polyfit_error_t error; - polyfit(x, y, num_points, degree, &error); return error; } // Evaluate at the specified point - polyfit_error_t error = polyfit_evaluate(poly, eval_x, result); + error = polyfit_evaluate(poly, eval_x, result); // Clean up polyfit_free(poly); @@ -306,6 +278,164 @@ polyfit_error_t polyfit_get_coefficients(const Polynomial* poly, float* coeffs, return POLYFIT_SUCCESS; } +/*============================================================================*/ +/* FIT QUALITY AND AUTO-DEGREE IMPLEMENTATIONS */ +/*============================================================================*/ + +polyfit_error_t polyfit_compute_residuals(const Polynomial* poly, + const float* x, const float* y, + int32_t num_points, + float* residuals) { + if (poly == NULL || x == NULL || y == NULL || residuals == NULL) { + return POLYFIT_ERROR_NULL_POINTER; + } + + if (!polyfit_is_valid(poly)) { + return POLYFIT_ERROR_INVALID_INPUT; + } + + if (num_points <= 0) { + return POLYFIT_ERROR_INVALID_INPUT; + } + + for (int32_t i = 0; i < num_points; i++) { + float y_hat; + polyfit_error_t err = polyfit_evaluate(poly, x[i], &y_hat); + if (err != POLYFIT_SUCCESS) { + return err; + } + residuals[i] = y_hat - y[i]; + } + + return POLYFIT_SUCCESS; +} + +polyfit_error_t polyfit_r_squared(const Polynomial* poly, const float* x, + const float* y, int32_t num_points, + float* r_squared) { + if (poly == NULL || x == NULL || y == NULL || r_squared == NULL) { + return POLYFIT_ERROR_NULL_POINTER; + } + + if (!polyfit_is_valid(poly)) { + return POLYFIT_ERROR_INVALID_INPUT; + } + + if (num_points <= 0) { + return POLYFIT_ERROR_INVALID_INPUT; + } + + // Compute mean of y + float y_mean = 0.0f; + for (int32_t i = 0; i < num_points; i++) { + y_mean += y[i]; + } + y_mean /= (float)num_points; + + // Compute total and residual sums of squares + float ss_tot = 0.0f; + float ss_res = 0.0f; + + for (int32_t i = 0; i < num_points; i++) { + float y_hat; + polyfit_error_t err = polyfit_evaluate(poly, x[i], &y_hat); + if (err != POLYFIT_SUCCESS) { + return err; + } + float diff_res = y_hat - y[i]; + float diff_tot = y[i] - y_mean; + ss_res += diff_res * diff_res; + ss_tot += diff_tot * diff_tot; + } + + if (ss_tot < 1e-20f) { + *r_squared = (ss_res < 1e-20f) ? 1.0f : 0.0f; + return POLYFIT_SUCCESS; + } + + *r_squared = 1.0f - (ss_res / ss_tot); + + return POLYFIT_SUCCESS; +} + +Polynomial* polyfit_best_degree(const float* x, const float* y, + int32_t num_points, int32_t max_degree, + int32_t* best_degree, + polyfit_error_t* error) { + if (x == NULL || y == NULL || best_degree == NULL) { + if (error != NULL) { + *error = POLYFIT_ERROR_NULL_POINTER; + } + return NULL; + } + + if (max_degree < 1 || max_degree > POLYFIT_MAX_DEGREE) { + if (error != NULL) { + *error = POLYFIT_ERROR_INVALID_DEGREE; + } + return NULL; + } + + if (num_points <= max_degree) { + if (error != NULL) { + *error = POLYFIT_ERROR_INSUFFICIENT_POINTS; + } + return NULL; + } + + Polynomial* best_poly = NULL; + float best_bic = 0.0f; + float fn = (float)num_points; + + for (int32_t d = 1; d <= max_degree; d++) { + polyfit_error_t local_err; + Polynomial* poly = polyfit(x, y, num_points, d, &local_err); + if (poly == NULL) { + continue; + } + + // Compute residual sum of squares + float ss_res = 0.0f; + for (int32_t i = 0; i < num_points; i++) { + float y_hat; + if (polyfit_evaluate(poly, x[i], &y_hat) == POLYFIT_SUCCESS) { + float diff = y_hat - y[i]; + ss_res += diff * diff; + } + } + + // BIC = n * ln(RSS/n) + (d+1) * ln(n) + float bic; + if (ss_res < 1e-30f) { + bic = -1e30f; // Near-perfect fit + } else { + bic = fn * logf(ss_res / fn) + (float)(d + 1) * logf(fn); + } + + if (best_poly == NULL || bic < best_bic) { + best_bic = bic; + *best_degree = d; + polyfit_free(best_poly); + best_poly = poly; + } else { + polyfit_free(poly); + } + } + + if (best_poly == NULL) { + if (error != NULL) { + *error = POLYFIT_ERROR_SINGULAR_MATRIX; + } + return NULL; + } + + if (error != NULL) { + *error = POLYFIT_SUCCESS; + } + + return best_poly; +} + /*============================================================================*/ /* UTILITY FUNCTION IMPLEMENTATIONS */ /*============================================================================*/ @@ -474,7 +604,7 @@ static polyfit_error_t validate_input_arrays(const float* x, const float* y, // Check for NaN or infinite values for (int32_t i = 0; i < num_points; i++) { - if (x[i] != x[i] || y[i] != y[i]) { // NaN check + if (x[i] - x[i] != 0.0f || y[i] - y[i] != 0.0f) { // NaN and Inf check return POLYFIT_ERROR_INVALID_INPUT; } } @@ -482,23 +612,3 @@ static polyfit_error_t validate_input_arrays(const float* x, const float* y, return POLYFIT_SUCCESS; } -static bool is_matrix_singular(float** A, int32_t n) { - const float determinant_threshold = 1e-12f; - - // Simple check: if any diagonal element is too small after partial pivoting - for (int32_t i = 0; i < n; i++) { - float max_in_column = 0.0f; - for (int32_t j = i; j < n; j++) { - float abs_val = polyfit_fabs(A[j][i]); - if (abs_val > max_in_column) { - max_in_column = abs_val; - } - } - - if (max_in_column < determinant_threshold) { - return true; - } - } - - return false; -} \ No newline at end of file diff --git a/polyfit.h b/polyfit.h index 24b95c9..2944b83 100644 --- a/polyfit.h +++ b/polyfit.h @@ -210,6 +210,67 @@ polyfit_error_t polyfit_eval_at(const float* x, const float* y, polyfit_error_t polyfit_get_coefficients(const Polynomial* poly, float* coeffs, int32_t size); +/*============================================================================*/ +/* FIT QUALITY AND AUTO-DEGREE FUNCTIONS */ +/*============================================================================*/ + +/** + * @brief Compute residuals between polynomial predictions and actual values + * @param poly Pointer to fitted Polynomial (must not be NULL) + * @param x Array of x values (must not be NULL) + * @param y Array of y values (must not be NULL) + * @param num_points Number of data points + * @param residuals Output array of size >= num_points; stores (y_hat[i] - y[i]) + * @return Error code indicating success or failure + */ +polyfit_error_t polyfit_compute_residuals(const Polynomial* poly, + const float* x, const float* y, + int32_t num_points, + float* residuals); + +/** + * @brief Compute R-squared goodness-of-fit metric + * @param poly Pointer to fitted Polynomial (must not be NULL) + * @param x Array of x values (must not be NULL) + * @param y Array of y values (must not be NULL) + * @param num_points Number of data points + * @param r_squared Output pointer for R² value (0 to 1, higher is better) + * @return Error code indicating success or failure + */ +polyfit_error_t polyfit_r_squared(const Polynomial* poly, const float* x, + const float* y, int32_t num_points, + float* r_squared); + +/** + * @brief Automatically select the best polynomial degree using BIC + * + * Tries degrees 1 to max_degree, fits each polynomial, and selects the degree + * minimising the Bayesian Information Criterion: + * BIC = n * ln(RSS/n) + (d+1) * ln(n) + * + * BIC penalises complexity, preventing overfitting without cross-validation. + * + * @param x Array of x values (must not be NULL) + * @param y Array of y values (must not be NULL) + * @param num_points Number of data points (must be > max_degree) + * @param max_degree Maximum degree to evaluate (1 to POLYFIT_MAX_DEGREE) + * @param best_degree Output pointer for selected degree (must not be NULL) + * @param error Optional pointer to store error code (can be NULL) + * @return Pointer to the best-fit Polynomial, or NULL on failure + * @note Caller is responsible for freeing with polyfit_free() + * + * @example + * int32_t deg; + * Polynomial *poly = polyfit_best_degree(x, y, n, 5, °, NULL); + * if (poly) { + * printf("Best degree: %d\n", deg); + * polyfit_free(poly); + * } + */ +Polynomial* polyfit_best_degree(const float* x, const float* y, + int32_t num_points, int32_t max_degree, + int32_t* best_degree, polyfit_error_t* error); + /*============================================================================*/ /* UTILITY FUNCTIONS */ /*============================================================================*/ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt new file mode 100644 index 0000000..56c87ce --- /dev/null +++ b/tests/CMakeLists.txt @@ -0,0 +1,17 @@ +include(FetchContent) + +FetchContent_Declare( + googletest + URL https://github.com/google/googletest/archive/refs/tags/v1.14.0.tar.gz + DOWNLOAD_EXTRACT_TIMESTAMP TRUE +) + +# Prevent overriding the parent project's compiler/linker settings on Windows +set(gtest_force_shared_crt ON CACHE BOOL "" FORCE) +FetchContent_MakeAvailable(googletest) + +add_executable(test_polyfit test_polyfit.cpp) +target_link_libraries(test_polyfit PRIVATE polyfit GTest::gtest_main) + +include(GoogleTest) +gtest_discover_tests(test_polyfit) diff --git a/tests/test_polyfit.cpp b/tests/test_polyfit.cpp new file mode 100644 index 0000000..dd79372 --- /dev/null +++ b/tests/test_polyfit.cpp @@ -0,0 +1,508 @@ +extern "C" { +#include "polyfit.h" +} + +#include +#include + +/*============================================================================*/ +/* SHARED TEST DATA */ +/*============================================================================*/ + +// y = 2x + 1 over 5 points +static const float kLinX[] = {0.0f, 1.0f, 2.0f, 3.0f, 4.0f}; +static const float kLinY[] = {1.0f, 3.0f, 5.0f, 7.0f, 9.0f}; +static const int kLinN = 5; + +// y = x^2 over 9 points +static const float kQuadX[] = {-3.0f, -2.0f, -1.0f, 0.0f, 1.0f, + 2.0f, 3.0f, 4.0f, 5.0f}; +static const float kQuadY[] = { 9.0f, 4.0f, 1.0f, 0.0f, 1.0f, + 4.0f, 9.0f, 16.0f, 25.0f}; +static const int kQuadN = 9; + +/*============================================================================*/ +/* INIT / FREE */ +/*============================================================================*/ + +TEST(PolyfitInit, ValidDegrees) { + Polynomial *p = polyfit_init(0); + ASSERT_NE(p, nullptr); + EXPECT_EQ(p->degree, 0); + EXPECT_TRUE(p->is_valid); + polyfit_free(p); + + p = polyfit_init(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(p->degree, 1); + polyfit_free(p); + + p = polyfit_init(POLYFIT_MAX_DEGREE); + ASSERT_NE(p, nullptr); + EXPECT_EQ(p->degree, POLYFIT_MAX_DEGREE); + polyfit_free(p); +} + +TEST(PolyfitInit, NegativeDegreeReturnsNull) { + EXPECT_EQ(polyfit_init(-1), nullptr); +} + +TEST(PolyfitInit, ExceedsMaxDegreeReturnsNull) { + EXPECT_EQ(polyfit_init(POLYFIT_MAX_DEGREE + 1), nullptr); +} + +TEST(PolyfitFree, NullIsSafe) { + EXPECT_NO_THROW(polyfit_free(nullptr)); +} + +/*============================================================================*/ +/* IS VALID */ +/*============================================================================*/ + +TEST(PolyfitIsValid, Null) { + EXPECT_FALSE(polyfit_is_valid(nullptr)); +} + +TEST(PolyfitIsValid, ValidPolyAfterInit) { + Polynomial *p = polyfit_init(2); + ASSERT_NE(p, nullptr); + EXPECT_TRUE(polyfit_is_valid(p)); + polyfit_free(p); +} + +TEST(PolyfitIsValid, NullCoefficientsIsInvalid) { + Polynomial p; + p.coefficients = nullptr; + p.degree = 1; + p.is_valid = true; + EXPECT_FALSE(polyfit_is_valid(&p)); +} + +/*============================================================================*/ +/* ERROR STRINGS */ +/*============================================================================*/ + +TEST(PolyfitErrorString, AllCodesReturnNonNull) { + EXPECT_STREQ(polyfit_error_string(POLYFIT_SUCCESS), "Success"); + EXPECT_NE(polyfit_error_string(POLYFIT_ERROR_NULL_POINTER), nullptr); + EXPECT_NE(polyfit_error_string(POLYFIT_ERROR_INVALID_DEGREE), nullptr); + EXPECT_NE(polyfit_error_string(POLYFIT_ERROR_MEMORY_ALLOC), nullptr); + EXPECT_NE(polyfit_error_string(POLYFIT_ERROR_SINGULAR_MATRIX), nullptr); + EXPECT_NE(polyfit_error_string(POLYFIT_ERROR_INSUFFICIENT_POINTS), nullptr); + EXPECT_NE(polyfit_error_string(POLYFIT_ERROR_INVALID_INPUT), nullptr); +} + +TEST(PolyfitErrorString, UnknownCodeReturnsNonNull) { + EXPECT_NE(polyfit_error_string((polyfit_error_t)999), nullptr); +} + +/*============================================================================*/ +/* LEAST SQUARES CORE */ +/*============================================================================*/ + +TEST(PolyfitLeastSquares, LinearFit) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_least_squares(kLinX, kLinY, kLinN, 1, p), POLYFIT_SUCCESS); + // y = 1 + 2x → coeffs[0]=1, coeffs[1]=2 + EXPECT_NEAR(p->coefficients[0], 1.0f, 1e-4f); + EXPECT_NEAR(p->coefficients[1], 2.0f, 1e-4f); + polyfit_free(p); +} + +TEST(PolyfitLeastSquares, QuadraticFit) { + Polynomial *p = polyfit_init(2); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_least_squares(kQuadX, kQuadY, kQuadN, 2, p), POLYFIT_SUCCESS); + // y = x^2 → coeffs[0]=0, coeffs[1]=0, coeffs[2]=1 + EXPECT_NEAR(p->coefficients[0], 0.0f, 1e-4f); + EXPECT_NEAR(p->coefficients[1], 0.0f, 1e-4f); + EXPECT_NEAR(p->coefficients[2], 1.0f, 1e-4f); + polyfit_free(p); +} + +TEST(PolyfitLeastSquares, NullXReturnsError) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_least_squares(nullptr, kLinY, kLinN, 1, p), + POLYFIT_ERROR_NULL_POINTER); + polyfit_free(p); +} + +TEST(PolyfitLeastSquares, NullYReturnsError) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_least_squares(kLinX, nullptr, kLinN, 1, p), + POLYFIT_ERROR_NULL_POINTER); + polyfit_free(p); +} + +TEST(PolyfitLeastSquares, NullResultReturnsError) { + EXPECT_EQ(polyfit_least_squares(kLinX, kLinY, kLinN, 1, nullptr), + POLYFIT_ERROR_NULL_POINTER); +} + +TEST(PolyfitLeastSquares, InsufficientPoints) { + Polynomial *p = polyfit_init(2); + ASSERT_NE(p, nullptr); + // num_points must be > degree; 2 points for degree 2 is not enough + EXPECT_EQ(polyfit_least_squares(kLinX, kLinY, 2, 2, p), + POLYFIT_ERROR_INSUFFICIENT_POINTS); + polyfit_free(p); +} + +TEST(PolyfitLeastSquares, NegativeDegreeReturnsError) { + Polynomial *p = polyfit_init(0); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_least_squares(kLinX, kLinY, kLinN, -1, p), + POLYFIT_ERROR_INVALID_DEGREE); + polyfit_free(p); +} + +TEST(PolyfitLeastSquares, ExceedsMaxDegreeReturnsError) { + Polynomial *p = polyfit_init(0); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_least_squares(kLinX, kLinY, kLinN, POLYFIT_MAX_DEGREE + 1, p), + POLYFIT_ERROR_INVALID_DEGREE); + polyfit_free(p); +} + +TEST(PolyfitLeastSquares, NaNInputRejected) { + float xi[] = {1.0f, 2.0f, 3.0f}; + float yi[] = {0.0f / 0.0f, 2.0f, 3.0f}; // NaN + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_least_squares(xi, yi, 3, 1, p), POLYFIT_ERROR_INVALID_INPUT); + polyfit_free(p); +} + +TEST(PolyfitLeastSquares, InfInputRejected) { + float xi[] = {1.0f, 2.0f, 3.0f}; + float yi[] = {1.0f / 0.0f, 2.0f, 3.0f}; // +Inf + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_least_squares(xi, yi, 3, 1, p), POLYFIT_ERROR_INVALID_INPUT); + polyfit_free(p); +} + +/*============================================================================*/ +/* EVALUATE */ +/*============================================================================*/ + +TEST(PolyfitEvaluate, LinearAtKnownPoints) { + Polynomial *p = polyfit(kLinX, kLinY, kLinN, 1, nullptr); + ASSERT_NE(p, nullptr); + + float result; + EXPECT_EQ(polyfit_evaluate(p, 0.0f, &result), POLYFIT_SUCCESS); + EXPECT_NEAR(result, 1.0f, 1e-4f); + EXPECT_EQ(polyfit_evaluate(p, 4.0f, &result), POLYFIT_SUCCESS); + EXPECT_NEAR(result, 9.0f, 1e-4f); + + polyfit_free(p); +} + +TEST(PolyfitEvaluate, NullPolyReturnsError) { + float result; + EXPECT_EQ(polyfit_evaluate(nullptr, 1.0f, &result), POLYFIT_ERROR_NULL_POINTER); +} + +TEST(PolyfitEvaluate, NullResultReturnsError) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_evaluate(p, 1.0f, nullptr), POLYFIT_ERROR_NULL_POINTER); + polyfit_free(p); +} + +TEST(PolyfitEvaluate, SmallLegitimateValueNotSnappedToZero) { + // y = 5e-7 * x — values well below the old 1e-6 absolute snap threshold + float xs[] = {0.0f, 1.0f, 2.0f, 3.0f, 4.0f}; + float ys[] = {0.0f, 5e-7f, 10e-7f, 15e-7f, 20e-7f}; + Polynomial *p = polyfit(xs, ys, 5, 1, nullptr); + ASSERT_NE(p, nullptr); + + float result; + EXPECT_EQ(polyfit_evaluate(p, 1.0f, &result), POLYFIT_SUCCESS); + EXPECT_NEAR(result, 5e-7f, 1e-8f); // must NOT be zero + + polyfit_free(p); +} + +/*============================================================================*/ +/* CONVENIENCE FUNCTIONS */ +/*============================================================================*/ + +TEST(PolyfitConvenience, PolyfitLinearFit) { + polyfit_error_t err; + Polynomial *p = polyfit(kLinX, kLinY, kLinN, 1, &err); + ASSERT_NE(p, nullptr); + EXPECT_EQ(err, POLYFIT_SUCCESS); + EXPECT_NEAR(p->coefficients[0], 1.0f, 1e-4f); + EXPECT_NEAR(p->coefficients[1], 2.0f, 1e-4f); + polyfit_free(p); +} + +TEST(PolyfitConvenience, PolyfitReturnsNullOnNullInput) { + polyfit_error_t err; + Polynomial *p = polyfit(nullptr, kLinY, kLinN, 1, &err); + EXPECT_EQ(p, nullptr); + EXPECT_EQ(err, POLYFIT_ERROR_NULL_POINTER); +} + +TEST(PolyfitConvenience, EvalAtQuadratic) { + float result; + EXPECT_EQ(polyfit_eval_at(kQuadX, kQuadY, kQuadN, 2, 4.0f, &result), + POLYFIT_SUCCESS); + EXPECT_NEAR(result, 16.0f, 1e-3f); +} + +TEST(PolyfitConvenience, EvalAtNullResultReturnsError) { + EXPECT_EQ(polyfit_eval_at(kQuadX, kQuadY, kQuadN, 2, 4.0f, nullptr), + POLYFIT_ERROR_NULL_POINTER); +} + +TEST(PolyfitConvenience, GetCoefficients) { + Polynomial *p = polyfit(kLinX, kLinY, kLinN, 1, nullptr); + ASSERT_NE(p, nullptr); + + float coeffs[2]; + EXPECT_EQ(polyfit_get_coefficients(p, coeffs, 2), POLYFIT_SUCCESS); + EXPECT_NEAR(coeffs[0], 1.0f, 1e-4f); + EXPECT_NEAR(coeffs[1], 2.0f, 1e-4f); + + polyfit_free(p); +} + +TEST(PolyfitConvenience, GetCoefficientsBufferTooSmall) { + Polynomial *p = polyfit(kLinX, kLinY, kLinN, 1, nullptr); + ASSERT_NE(p, nullptr); + + float coeffs[1]; // degree=1 needs 2 elements + EXPECT_EQ(polyfit_get_coefficients(p, coeffs, 1), POLYFIT_ERROR_INVALID_INPUT); + + polyfit_free(p); +} + +/*============================================================================*/ +/* RESIDUALS */ +/*============================================================================*/ + +TEST(PolyfitResiduals, PerfectFitNearZero) { + Polynomial *p = polyfit(kQuadX, kQuadY, kQuadN, 2, nullptr); + ASSERT_NE(p, nullptr); + + float residuals[9]; + EXPECT_EQ(polyfit_compute_residuals(p, kQuadX, kQuadY, kQuadN, residuals), + POLYFIT_SUCCESS); + for (int i = 0; i < kQuadN; i++) { + EXPECT_NEAR(residuals[i], 0.0f, 1e-3f) << "residuals[" << i << "]"; + } + + polyfit_free(p); +} + +TEST(PolyfitResiduals, NullPolyReturnsError) { + float residuals[5]; + EXPECT_EQ(polyfit_compute_residuals(nullptr, kLinX, kLinY, kLinN, residuals), + POLYFIT_ERROR_NULL_POINTER); +} + +TEST(PolyfitResiduals, NullXReturnsError) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + float residuals[5]; + EXPECT_EQ(polyfit_compute_residuals(p, nullptr, kLinY, kLinN, residuals), + POLYFIT_ERROR_NULL_POINTER); + polyfit_free(p); +} + +TEST(PolyfitResiduals, NullYReturnsError) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + float residuals[5]; + EXPECT_EQ(polyfit_compute_residuals(p, kLinX, nullptr, kLinN, residuals), + POLYFIT_ERROR_NULL_POINTER); + polyfit_free(p); +} + +TEST(PolyfitResiduals, NullResidualsOutputReturnsError) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_compute_residuals(p, kLinX, kLinY, kLinN, nullptr), + POLYFIT_ERROR_NULL_POINTER); + polyfit_free(p); +} + +/*============================================================================*/ +/* R-SQUARED */ +/*============================================================================*/ + +TEST(PolyfitRSquared, PerfectFitIsOne) { + Polynomial *p = polyfit(kQuadX, kQuadY, kQuadN, 2, nullptr); + ASSERT_NE(p, nullptr); + + float r2; + EXPECT_EQ(polyfit_r_squared(p, kQuadX, kQuadY, kQuadN, &r2), POLYFIT_SUCCESS); + EXPECT_NEAR(r2, 1.0f, 1e-5f); + + polyfit_free(p); +} + +TEST(PolyfitRSquared, WrongDegreeGivesLowerScore) { + // Fit degree-1 to quadratic data — R² should be well below 1 + Polynomial *p = polyfit(kQuadX, kQuadY, kQuadN, 1, nullptr); + ASSERT_NE(p, nullptr); + + float r2; + EXPECT_EQ(polyfit_r_squared(p, kQuadX, kQuadY, kQuadN, &r2), POLYFIT_SUCCESS); + EXPECT_LT(r2, 0.99f); + + polyfit_free(p); +} + +TEST(PolyfitRSquared, NullPolyReturnsError) { + float r2; + EXPECT_EQ(polyfit_r_squared(nullptr, kLinX, kLinY, kLinN, &r2), + POLYFIT_ERROR_NULL_POINTER); +} + +TEST(PolyfitRSquared, NullXReturnsError) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + float r2; + EXPECT_EQ(polyfit_r_squared(p, nullptr, kLinY, kLinN, &r2), + POLYFIT_ERROR_NULL_POINTER); + polyfit_free(p); +} + +TEST(PolyfitRSquared, NullYReturnsError) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + float r2; + EXPECT_EQ(polyfit_r_squared(p, kLinX, nullptr, kLinN, &r2), + POLYFIT_ERROR_NULL_POINTER); + polyfit_free(p); +} + +TEST(PolyfitRSquared, NullOutputReturnsError) { + Polynomial *p = polyfit_init(1); + ASSERT_NE(p, nullptr); + EXPECT_EQ(polyfit_r_squared(p, kLinX, kLinY, kLinN, nullptr), + POLYFIT_ERROR_NULL_POINTER); + polyfit_free(p); +} + +/*============================================================================*/ +/* BEST DEGREE AUTO-SELECTION */ +/*============================================================================*/ + +TEST(PolyfitBestDegree, SelectsLinearForLinearData) { + int32_t deg; + Polynomial *p = polyfit_best_degree(kLinX, kLinY, kLinN, 3, °, nullptr); + ASSERT_NE(p, nullptr); + EXPECT_EQ(deg, 1); + polyfit_free(p); +} + +TEST(PolyfitBestDegree, SelectsQuadraticForQuadraticData) { + int32_t deg; + Polynomial *p = polyfit_best_degree(kQuadX, kQuadY, kQuadN, 5, °, nullptr); + ASSERT_NE(p, nullptr); + EXPECT_EQ(deg, 2); + polyfit_free(p); +} + +TEST(PolyfitBestDegree, ReturnedPolyEvaluatesCorrectly) { + int32_t deg; + Polynomial *p = polyfit_best_degree(kQuadX, kQuadY, kQuadN, 4, °, nullptr); + ASSERT_NE(p, nullptr); + + float result; + EXPECT_EQ(polyfit_evaluate(p, 4.0f, &result), POLYFIT_SUCCESS); + EXPECT_NEAR(result, 16.0f, 1e-2f); + + polyfit_free(p); +} + +TEST(PolyfitBestDegree, NullXReturnsError) { + int32_t deg; + polyfit_error_t err; + EXPECT_EQ(polyfit_best_degree(nullptr, kQuadY, kQuadN, 4, °, &err), nullptr); + EXPECT_EQ(err, POLYFIT_ERROR_NULL_POINTER); +} + +TEST(PolyfitBestDegree, NullYReturnsError) { + int32_t deg; + polyfit_error_t err; + EXPECT_EQ(polyfit_best_degree(kQuadX, nullptr, kQuadN, 4, °, &err), nullptr); + EXPECT_EQ(err, POLYFIT_ERROR_NULL_POINTER); +} + +TEST(PolyfitBestDegree, NullBestDegreeOutputReturnsError) { + polyfit_error_t err; + EXPECT_EQ(polyfit_best_degree(kQuadX, kQuadY, kQuadN, 4, nullptr, &err), nullptr); + EXPECT_EQ(err, POLYFIT_ERROR_NULL_POINTER); +} + +TEST(PolyfitBestDegree, InsufficientPoints) { + // 3 points but max_degree=5 requires at least 6 points + float xi[] = {0.0f, 1.0f, 2.0f}; + float yi[] = {1.0f, 2.0f, 4.0f}; + int32_t deg; + polyfit_error_t err; + Polynomial *p = polyfit_best_degree(xi, yi, 3, 5, °, &err); + EXPECT_EQ(p, nullptr); + EXPECT_EQ(err, POLYFIT_ERROR_INSUFFICIENT_POINTS); +} + +/*============================================================================*/ +/* UTILITY FUNCTIONS */ +/*============================================================================*/ + +TEST(PolyfitPow, ZeroExponent) { + EXPECT_FLOAT_EQ(polyfit_pow(5.0f, 0), 1.0f); + EXPECT_FLOAT_EQ(polyfit_pow(0.0f, 0), 1.0f); // 0^0 = 1 by convention +} + +TEST(PolyfitPow, ZeroBase) { + EXPECT_FLOAT_EQ(polyfit_pow(0.0f, 3), 0.0f); +} + +TEST(PolyfitPow, OneBase) { + EXPECT_FLOAT_EQ(polyfit_pow(1.0f, 99), 1.0f); +} + +TEST(PolyfitPow, NegativeOneBase) { + EXPECT_FLOAT_EQ(polyfit_pow(-1.0f, 2), 1.0f); // even exponent + EXPECT_FLOAT_EQ(polyfit_pow(-1.0f, 3), -1.0f); // odd exponent +} + +TEST(PolyfitPow, PositiveExponents) { + EXPECT_NEAR(polyfit_pow(2.0f, 3), 8.0f, 1e-6f); + EXPECT_NEAR(polyfit_pow(3.0f, 4), 81.0f, 1e-5f); +} + +TEST(PolyfitPow, NegativeExponent) { + EXPECT_NEAR(polyfit_pow(2.0f, -2), 0.25f, 1e-7f); +} + +TEST(PolyfitFabs, Positive) { EXPECT_FLOAT_EQ(polyfit_fabs(3.14f), 3.14f); } +TEST(PolyfitFabs, Negative) { EXPECT_FLOAT_EQ(polyfit_fabs(-3.14f), 3.14f); } +TEST(PolyfitFabs, Zero) { EXPECT_FLOAT_EQ(polyfit_fabs(0.0f), 0.0f); } +TEST(PolyfitFabs, NaN) { + float nan_val = 0.0f / 0.0f; + EXPECT_TRUE(std::isnan(polyfit_fabs(nan_val))); +} + +TEST(PolyfitIsNearlyZero, BelowThresholdIsTrue) { + EXPECT_TRUE(polyfit_is_nearly_zero(1e-8f, 1e-6f)); +} + +TEST(PolyfitIsNearlyZero, AboveThresholdIsFalse) { + EXPECT_FALSE(polyfit_is_nearly_zero(1e-5f, 1e-6f)); +} + +TEST(PolyfitIsNearlyZero, ExactlyAtThresholdIsFalse) { + // strictly less than threshold required + EXPECT_FALSE(polyfit_is_nearly_zero(1e-6f, 1e-6f)); +}