Changeset 42319 for branches/eam_branches/psLib.20230123/src/math/psSpline.c

Timestamp:

Jan 23, 2023, 7:10:15 AM (3 years ago)

Author:

eugene

Message:

workin on spline implementation

File:

• branches/eam_branches/psLib.20230123/src/math/psSpline.c (modified) (9 diffs)

branches/eam_branches/psLib.20230123/src/math/psSpline.c

@@ -1 +1 @@
 /** @file psSpline.c
-*
-*  @brief Contains basic spline allocation, deallocation, fitting,
-*         and evaluation routines.
-*
-*  This file contains the routines that allocate, free, and evaluate splines.
-*
-*  @version $Revision: 1.158 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2007-03-14 02:36:28 $
-*
-*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
+    re-written by EAM 2023.01.23
 */
 
@@ -16 +7 @@
 #endif
 
-/*****************************************************************************/
-/*  INCLUDE FILES                                                            */
-/*****************************************************************************/
 #include <stdio.h>
 #include <stdbool.h>
@@ -36 +24 @@
 #include "psMathUtils.h"
 
-/*****************************************************************************/
-/* DEFINE STATEMENTS                                                         */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* TYPE DEFINITIONS                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* GLOBAL VARIABLES                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* FILE STATIC VARIABLES                                                     */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* FUNCTION IMPLEMENTATION - LOCAL                                           */
-/*****************************************************************************/
-
-static void spline1DFree(psSpline1D *tmpSpline)
-{
-    if (tmpSpline == NULL) {
-        return;
-    }
+static void psSpline1DFree(psSpline1D *tmpSpline)
+{
+    if (tmpSpline == NULL) return;
 
     if (tmpSpline->spline != NULL) {
@@ -78 +44 @@
 }
 
-
-static void PS_POLY1D_PRINT(
-    psPolynomial1D *poly)
-{
-    printf("-------------- PS_POLY1D_PRINT() --------------\n");
-    printf("poly->nX is %d\n", poly->nX);
-    for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) {
-        printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]);
-    }
-}
-
 static void PS_PRINT_SPLINE2(psSpline1D *mySpline)
 {
@@ -94 +49 @@
     printf("mySpline->n is %d\n", mySpline->n);
     for (psS32 i = 0 ; i < mySpline->n ; i++) {
-        PS_POLY1D_PRINT(mySpline->spline[i]);
-    }
-    PS_VECTOR_PRINT_F32(mySpline->knots);
+        // print the knots
+    }
 }
 
 /*****************************************************************************
-CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a
-tabulated function at n points, this routine calculates the second
-derivatives of the interpolating cubic splines at those n points.
-
-The first and second derivatives at the endpoints were previously undefined in
-the SDR.  From bugzilla #???, they are required to be 0.0, implementing natural
-splines.
-
-Endpoints are defined by
-    PS_LEFT_SPLINE_DERIV
-    PS_RIGHT_SPLINE_DERIV
-
-This routine assumes that vectors x and y are of the appropriate types/sizes
-(F32).
-
-XXX: use recycled vectors for internal data.
-XXX: do an F64 version?
+CalculateSecondDerivs(): Given a set of x,y vectors corresponding to the spline knots, this
+routine calculates the second derivatives of the interpolating cubic splines at those n points.
+
+The boundary conditions may be:
+* first derivative specified OR
+* second derivatives = 0
+
+One or the other condition may be true for each of the upper and lower bounds
+
+NOTE: vectors must be F32
+
  *****************************************************************************/
-#define PS_LEFT_SPLINE_DERIV 0.0
-#define PS_RIGHT_SPLINE_DERIV 0.0
+
 static psF32 *calculateSecondDerivs(
-    const psVector* x,                  ///< Ordinates
-    const psVector* y)                  ///< Coordinates
+    const psVector* x, ///< Ordinates
+    const psVector* y, ///< Coordinates
+    psF32 dyLower; // if not NAN, lower-bound 1st derivative is defined
+    psF32 dyUpper; // if not NAN, lower-bound 1st derivative is defined
+    )                  
 {
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
@@ -190 +140 @@
 
 /*****************************************************************************
-psVectorFitSpline1D(): given a set of x/y vectors, this routine generates the
+psSpline1DFitVector(): given a set of x/y vectors, this routine generates the
 linear or cublic splines which satisfy those data points.
 
@@ -209 +159 @@
 XXX: What types must be supported?
  *****************************************************************************/
-psSpline1D *psVectorFitSpline1D(
+psSpline1D *psSpline1DFitVector(
     const psVector* x,                  ///< Ordinates.
     const psVector* y)                  ///< Coordinates.
@@ -361 +311 @@
 psF64?
  *****************************************************************************/
-float psSpline1DEval(
-    const psSpline1D *spline,
+float psSpline1DEval_Old(
+    const psSpline1D_Old *spline,
     float x)
 {
@@ -397 +347 @@
 }
 
+// XXX EAM : changing implementation to use yKnot, d2yKnot instead of polynomials
+float psSpline1DEval_New(
+    const psSpline1D *spline,
+    float x)
+{
+    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(spline, NAN);
+    PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN);
+    PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NAN);
+
+    psS32 n = spline->n;
+    if ((x < spline->knots->data.F32[0]) || (x > spline->knots->data.F32[spline->knots->n-1])) {
+        // If x is outside the range of spline->knots, generate a warning
+        // message, then return the left, or right, endpoint.
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f) (%d).",
+                 x, spline->knots->data.F32[0], spline->knots->data.F32[n-1], n);
+
+        psS32 binNum = (x < spline->knots->data.F32[0]) ? 0 : n-1;
+        psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
+        return(psPolynomial1DEval(spline->spline[binNum], x));
+    }
+
+    psScalar tmpScalar;
+    tmpScalar.type.type = PS_TYPE_F32;
+    tmpScalar.data.F32 = x;
+    psVectorBinaryDisectResult result;
+    psS32 binNum = psVectorBinaryDisect(&result, spline->knots, &tmpScalar);
+    if (result != PS_BINARY_DISECT_PASS) {
+        psError(PS_ERR_UNKNOWN, false, "Could not perform bin dissection on spline->knots.\n");
+        return(NAN);
+    }
+
+    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
+    return(psPolynomial1DEval(spline->spline[binNum], x));
+}
+
 
 /*****************************************************************************