From 9610bac375ce843b02cc2de0388b5a97ccbb8ab4 Mon Sep 17 00:00:00 2001
From: Jakob van Santen <jvansanten@gmail.com>
Date: Mon, 15 Apr 2024 15:49:55 +0200
Subject: [PATCH] Ignore partially-supported knots in spacing bounds (#43)

* Use only knots inside of support for spacing bound

It is perfectly legitimate for the first and last `order` knots to
be identical.  Exclude these from the minimum and maximum knot
spacing calculation to prevent elements of rmin_sep from becoming
infinite.

* Bump version
---
 CMakeLists.txt                            |  2 +-
 include/photospline/detail/bspline_eval.h |  2 +-
 include/photospline/detail/convolve.h     |  2 +-
 include/photospline/detail/fit.h          |  2 +-
 include/photospline/detail/fitsio.h       |  2 +-
 include/photospline/detail/permute.h      |  2 +-
 include/photospline/splinetable.h         | 12 ++++++------
 7 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7243eb4..4937ef9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,7 +1,7 @@
 cmake_minimum_required (VERSION 3.1.0 FATAL_ERROR)
 cmake_policy(VERSION 3.1.0)
 
-project (photospline VERSION 2.3.0 LANGUAGES C CXX)
+project (photospline VERSION 2.3.1 LANGUAGES C CXX)
 
 SET(CMAKE_CXX_STANDARD 11)
 SET(CMAKE_C_STANDARD 99)
diff --git a/include/photospline/detail/bspline_eval.h b/include/photospline/detail/bspline_eval.h
index 2b5138d..a3470cf 100644
--- a/include/photospline/detail/bspline_eval.h
+++ b/include/photospline/detail/bspline_eval.h
@@ -29,7 +29,7 @@ bool splinetable<Alloc>::searchcenters(const double* x, int* centers) const
 		}
 		
 		uint32_t min = order[i];
-		uint32_t max = nknots[i]-2;
+		uint32_t max = naxes[i];
 		double diff = x[i] - knots[i][min];
 		uint32_t hi = diff*rmin_sep[i];
 		uint32_t lo = diff*rmax_sep[i];
diff --git a/include/photospline/detail/convolve.h b/include/photospline/detail/convolve.h
index 6a5b597..af3b40c 100644
--- a/include/photospline/detail/convolve.h
+++ b/include/photospline/detail/convolve.h
@@ -146,7 +146,7 @@ void splinetable<Alloc>::convolve(const uint32_t dim, const double* conv_knots,
 	}
 
 	// Update knot separations
-	dknot_bounds();
+	fill_knot_spacing_bounds();
 	
 	/*
 	 * NB: A monotonic function remains monotonic after convolution
diff --git a/include/photospline/detail/fit.h b/include/photospline/detail/fit.h
index c4f282b..29cb203 100644
--- a/include/photospline/detail/fit.h
+++ b/include/photospline/detail/fit.h
@@ -104,7 +104,7 @@ void splinetable<Alloc>::fit(const ::ndsparse& data,
 		std::copy(knots[i].begin(),knots[i].end(),this->knots[i]);
 		dummy_knots[i]=&this->knots[i][0];
 	}
-	dknot_bounds();
+	fill_knot_spacing_bounds();
 	//same deal for the coordinates
 	std::unique_ptr<const double*[]> dummy_coords(new const double*[ndim]);
 	for(uint32_t i=0; i<ndim; i++)
diff --git a/include/photospline/detail/fitsio.h b/include/photospline/detail/fitsio.h
index fdbaaad..21e3422 100644
--- a/include/photospline/detail/fitsio.h
+++ b/include/photospline/detail/fitsio.h
@@ -382,7 +382,7 @@ bool splinetable<Alloc>::read_fits_core(fitsfile* fits, const std::string& fileP
 		}
 	}
 
-	dknot_bounds();
+	fill_knot_spacing_bounds();
 	
 	if(error!=0)
 		throw std::runtime_error("Error reading "+filePath+": Error "+std::to_string(error));
diff --git a/include/photospline/detail/permute.h b/include/photospline/detail/permute.h
index b1cf166..696d3af 100644
--- a/include/photospline/detail/permute.h
+++ b/include/photospline/detail/permute.h
@@ -84,7 +84,7 @@ void splinetable<Alloc>::permuteDimensions(const std::vector<size_t>& permutatio
 	std::copy(t_coefficients.get(),t_coefficients.get()+ncoeffs,coefficients);
 
 	// Update knot separations
-	dknot_bounds();
+	fill_knot_spacing_bounds();
 }
 	
 } //namespace photospline
diff --git a/include/photospline/splinetable.h b/include/photospline/splinetable.h
index 44c8e9e..2ee551a 100644
--- a/include/photospline/splinetable.h
+++ b/include/photospline/splinetable.h
@@ -277,16 +277,16 @@ class splinetable{
       lastKnots[nknots[inputDim]-1]=2*lastKnots[nknots[inputDim]-2]-lastKnots[nknots[inputDim]-3];
     }
 
-    rmin_sep=allocate<double>(ndim);
-    rmax_sep=allocate<double>(ndim);
-		dknot_bounds();
-
     //set naxes
     naxes=allocate<uint64_t>(ndim);
     for(unsigned int i=0; i<inputDim; i++)
       naxes[i] = tables.front()->get_ncoeffs(i);
     naxes[inputDim]=tables.size();
 
+    rmin_sep=allocate<double>(ndim);
+    rmax_sep=allocate<double>(ndim);
+    fill_knot_spacing_bounds();
+
     //copy coefficients
     unsigned long nCoeffs=std::accumulate(naxes, naxes+ndim, 1UL, std::multiplies<uint64_t>());
     unsigned long nInputCoeffs=std::accumulate(naxes, naxes+ndim-1, 1UL, std::multiplies<uint64_t>());
@@ -856,10 +856,10 @@ class splinetable{
 	///Write to a file
 	void write_fits_core(fitsfile*) const;
 
-	void dknot_bounds() {
+	void fill_knot_spacing_bounds() {
 		for (uint32_t i = 0; i < ndim; i++) {
 			uint32_t min = order[i];
-			uint32_t max = nknots[i]-2;
+			uint32_t max = naxes[i];
 			double mini = DBL_MAX;
 			double maxi = 0;
 			for (uint32_t j = min; j < max; j++) {