Cmake repo update

.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

.gitignore

@@ -50,3 +50,4 @@ modules.order
 Module.symvers
 Mkfile.old
 dkms.conf
+build/

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)

polyfit.c

@@ -14,6 +14,7 @@
  */
 
 #include "polyfit.h"
+#include <math.h>
 
 /*============================================================================*/
 /* 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,31 +604,11 @@ 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;
     }
   }
 
   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;
-}