<commit_before>
#ifndef PHOTOSPLINE_SPLINETABLE_H
#define PHOTOSPLINE_SPLINETABLE_H

#include <algorithm>
#include <cassert>
#include <memory>
#include <numeric>
#include <sstream>
#include <vector>

#include "photospline/bspline.h"
#include "photospline/detail/simd.h"

#include <string.h>
#include <fitsio.h>
#include <fitsio2.h>

#ifdef PHOTOSPLINE_INCLUDES_SPGLAM
#include "photospline/detail/splineutil.h"
#include "photospline/detail/glam.h"
#endif //PHOTOSPLINE_INCLUDES_SPGLAM

#include <iostream>

namespace photospline{
	
#ifdef PHOTOSPLINE_INCLUDES_SPGLAM
///A more user-friendly version of the C ndsparse
///Note that this does not mean entirely safe; as all of the internals are still
///exposed, leaving direct modification of the structure possible.
struct ndsparse : public ::ndsparse{
	size_t entriesInserted;
	///\brief 'Raw' constructor which leaves the object essentially uninitialized
	ndsparse(){
		ndim=0;
		rows=0;
		i=nullptr;
		ranges=nullptr;
		x=nullptr;
		entriesInserted=0;
	}
	///\brief Proper constructor which allocates memory
	///\param rows the number of nonzero entries for which space will be allocated
	///\param ndim the number of dimensions
	ndsparse(size_t rows, size_t ndim){
		if(!ndim)
			throw std::logic_error("Tried to allocate an ndsparse with dimension 0");
		if(!rows)
			throw std::logic_error("Tried to allocate an ndsparse with 0 entries");
		if(ndsparse_allocate(this,rows,ndim)!=0)
			throw std::bad_alloc();
		entriesInserted=0;
	}
	~ndsparse(){
		ndsparse_free(this);
	}
	///\brief A convenience interface for adding entries.
	///It is a bad idea to mix use of this interface with direct alterations to
	///the x and i arrays.
	void insertEntry(double value, unsigned int* indices){
		//assert(entriesInserted+1<rows && "Attempt to insert and entry into a full ndsparse");
		if(!(entriesInserted<rows))
			throw std::runtime_error("Attempt to insert an entry into a full ndsparse");
		x[entriesInserted]=value;
		for(size_t j=0; j<ndim; j++){
			i[j][entriesInserted] = indices[j];
			ranges[j] = std::max(ranges[j],i[j][entriesInserted]+1);
		}
		entriesInserted++;
	}
};
#endif //PHOTOSPLINE_INCLUDES_SPGLAM

namespace detail{
	///Work-around to replace passing 2D VLAs
	template<typename T>
	struct buffer2d{
		T* buf;
		size_t dim1;
		buffer2d(T* b, size_t d1):buf(b),dim1(d1){}
		
		T* operator[](size_t idx1){
			return(buf+dim1*idx1); }
		const T* operator[](size_t idx1) const{ return(buf+dim1*idx1); }
	};
	///A simple view into un-owned data.
	template<typename T>
	struct array_view{
	public:
		using value_type=const T;
		using reference=const T&;
		using const_reference=const T&;
		using pointer=const T*;
		using const_pointer=const T*;
		using iterator=const T*;
		using const_iterator=const T*;
		using difference_type=typename std::iterator_traits<iterator>::difference_type;
		using size_type=size_t;
	private:
		const_pointer d;
		size_type s;
	public:
		array_view():d(nullptr),s(0){}
		array_view(const_pointer d, size_type s):d(d),s(s){}
		void reset(const_pointer d, size_type s){
			this->d=d;
			this->s=s;
		}
		size_type size() const{ return(s); }
		size_type max_size() const{ return(s); }
		bool empty() const{ return(s==0); }
		const_pointer data() const{ return(d); }
		const_reference operator[](size_type i) const{ return(*(d+i)); }
		iterator begin() const{ return(d); }
		const_iterator cbegin() const{ return(d); }
		iterator end() const{ return(d+s); }
		const_iterator cend() const{ return(d+s); }
	};
}
	
template<typename Alloc = std::allocator<void> >
class splinetable{
public:
	typedef Alloc allocator_type;
	typedef std::allocator_traits<allocator_type> allocator_traits;
	
	typedef typename allocator_traits::template rebind_traits<uint32_t>::pointer uint32_t_ptr;
	typedef typename allocator_traits::template rebind_traits<uint64_t>::pointer uint64_t_ptr;
	typedef typename allocator_traits::template rebind_traits<float>::pointer float_ptr;
	typedef typename allocator_traits::template rebind_traits<double>::pointer double_ptr;
	typedef typename allocator_traits::template rebind_traits<double_ptr>::pointer double_ptr_ptr;
	typedef typename allocator_traits::template rebind_traits<char>::pointer char_ptr;
	typedef typename allocator_traits::template rebind_traits<char_ptr>::pointer char_ptr_ptr;
	typedef typename allocator_traits::template rebind_traits<char_ptr_ptr>::pointer char_ptr_ptr_ptr;
	
	///Construct an empty splinetable.
	///The resulting object is useful only for calling read_fits, read_fits_mem, or fit.
	explicit splinetable(allocator_type alloc=Alloc()):
	ndim(0),order(NULL),knots(NULL),nknots(NULL),extents(NULL),periods(NULL),
	coefficients(NULL),naxes(NULL),strides(NULL),naux(0),aux(NULL),allocator(alloc)
	{}
	
	///Construct a splinetable from serialized data previously stored in a FITS file.
	///\param filePath the path to the input file
	explicit splinetable(const std::string& filePath, allocator_type alloc=Alloc()):
	ndim(0),order(NULL),knots(NULL),nknots(NULL),extents(NULL),periods(NULL),
	coefficients(NULL),naxes(NULL),strides(NULL),naux(0),aux(NULL),allocator(alloc)
	{
		read_fits(filePath);
	}
	
	splinetable(splinetable&& other):
	ndim(other.ndim),order(std::move(other.order)),knots(std::move(other.knots)),
	nknots(other.nknots),extents(std::move(other.extents)),
	periods(std::move(other.periods)),coefficients(std::move(other.coefficients)),
	naxes(std::move(other.naxes)),strides(std::move(other.strides)),
	naux(other.naux),aux(std::move(other.aux)),
	allocator(std::move(other.allocator))
	{
		other.ndim=0;
		other.order=NULL;
		other.knots=NULL;
		other.nknots=NULL;
		other.extents=NULL;
		other.periods=NULL;
		other.coefficients=NULL;
		other.naxes=NULL;
		other.strides=NULL;
		other.naux=0;
		other.aux=NULL;
		other.allocator=Alloc();
	}
	
	~splinetable(){
		if(ndim){
			uint64_t ncoeffs=strides[0]*naxes[0];
			for(uint32_t i=0; i<ndim; i++)
				deallocate(knots[i]-order[i],nknots[i]+2*order[i]);
			deallocate(knots,ndim);
			deallocate(nknots,ndim);
			deallocate(order,ndim);
			if(extents){
				deallocate(extents[0],2*ndim);
				deallocate(extents,ndim);
			}
			if(periods)
				deallocate(periods,ndim);
			deallocate(coefficients,ncoeffs);
			deallocate(naxes,ndim);
			deallocate(strides,ndim);
			for(uint32_t i=0; i<naux; i++){
				deallocate(aux[i][0],strlen(&aux[i][0][0])+1);
				deallocate(aux[i][1],strlen(&aux[i][1][0])+1);
				deallocate(aux[i],2);
			}
			deallocate(aux,naux);
		}
	}
	
	splinetable& operator=(splinetable&& other){
		if(&other==this)
			return(*this);
		using std::swap;
		swap(ndim,other.ndim);
		swap(order,other.order);
		swap(knots,other.knots);
		swap(nknots,other.nknots);
		swap(extents,other.extents);
		swap(periods,other.periods);
		swap(coefficients,other.coefficients);
		swap(naxes,other.naxes);
		swap(strides,other.strides);
		swap(naux,other.naux);
		swap(aux,other.aux);
		swap(allocator,other.allocator);
		return(*this);
	}
	
	///Estimate the memory needed to load an existing spline.
	///This function is intended for users using specialized allocators for which
	///it is useful to know the total memory required before constructing the
	///allocator and splinetable. If a convolution will be applied to the spline
	///after it is loaded, doing so will increase its memory requirements, which
	///this function can also take into account.
	///\param filePath the path to the input FITS file
	///\param n_convolution_knots the number of knots possessed by the other
	///       spline with which one dimension of this spline will be convolved.
	///       A value of one is equivalent to no convolution.
	///\param convolution_dimension the dimension of this spline which will have
	///       a convolution applied.
	static size_t estimateMemory(const std::string& filePath,
	                             uint32_t n_convolution_knots = 1,
	                             uint32_t convolution_dimension = 0);
	
	///Read from a FITS file
	///\param path the path to the input file
	bool read_fits(const std::string& path);
	
	///Read from a FITS 'file' in a memory buffer
	///\param the input data buffer
	///\param buffer_size the length of the input buffer
	bool read_fits_mem(void* buffer, size_t buffer_size);
	
	///Write to a FITS file
	///\param path the path to the output file
	void write_fits(const std::string& path) const;
	
	///Write to a FITS memory 'file'.<commit_msg>Add an equality operator for splinetable<commit_after>
#ifndef PHOTOSPLINE_SPLINETABLE_H
#define PHOTOSPLINE_SPLINETABLE_H

#include <algorithm>
#include <cassert>
#include <memory>
#include <numeric>
#include <sstream>
#include <vector>

#include "photospline/bspline.h"
#include "photospline/detail/simd.h"

#include <string.h>
#include <fitsio.h>
#include <fitsio2.h>

#ifdef PHOTOSPLINE_INCLUDES_SPGLAM
#include "photospline/detail/splineutil.h"
#include "photospline/detail/glam.h"
#endif //PHOTOSPLINE_INCLUDES_SPGLAM

#include <iostream>

namespace photospline{
	
#ifdef PHOTOSPLINE_INCLUDES_SPGLAM
///A more user-friendly version of the C ndsparse
///Note that this does not mean entirely safe; as all of the internals are still
///exposed, leaving direct modification of the structure possible.
struct ndsparse : public ::ndsparse{
	size_t entriesInserted;
	///\brief 'Raw' constructor which leaves the object essentially uninitialized
	ndsparse(){
		ndim=0;
		rows=0;
		i=nullptr;
		ranges=nullptr;
		x=nullptr;
		entriesInserted=0;
	}
	///\brief Proper constructor which allocates memory
	///\param rows the number of nonzero entries for which space will be allocated
	///\param ndim the number of dimensions
	ndsparse(size_t rows, size_t ndim){
		if(!ndim)
			throw std::logic_error("Tried to allocate an ndsparse with dimension 0");
		if(!rows)
			throw std::logic_error("Tried to allocate an ndsparse with 0 entries");
		if(ndsparse_allocate(this,rows,ndim)!=0)
			throw std::bad_alloc();
		entriesInserted=0;
	}
	~ndsparse(){
		ndsparse_free(this);
	}
	///\brief A convenience interface for adding entries.
	///It is a bad idea to mix use of this interface with direct alterations to
	///the x and i arrays.
	void insertEntry(double value, unsigned int* indices){
		//assert(entriesInserted+1<rows && "Attempt to insert and entry into a full ndsparse");
		if(!(entriesInserted<rows))
			throw std::runtime_error("Attempt to insert an entry into a full ndsparse");
		x[entriesInserted]=value;
		for(size_t j=0; j<ndim; j++){
			i[j][entriesInserted] = indices[j];
			ranges[j] = std::max(ranges[j],i[j][entriesInserted]+1);
		}
		entriesInserted++;
	}
};
#endif //PHOTOSPLINE_INCLUDES_SPGLAM

namespace detail{
	///Work-around to replace passing 2D VLAs
	template<typename T>
	struct buffer2d{
		T* buf;
		size_t dim1;
		buffer2d(T* b, size_t d1):buf(b),dim1(d1){}
		
		T* operator[](size_t idx1){
			return(buf+dim1*idx1); }
		const T* operator[](size_t idx1) const{ return(buf+dim1*idx1); }
	};
	///A simple view into un-owned data.
	template<typename T>
	struct array_view{
	public:
		using value_type=const T;
		using reference=const T&;
		using const_reference=const T&;
		using pointer=const T*;
		using const_pointer=const T*;
		using iterator=const T*;
		using const_iterator=const T*;
		using difference_type=typename std::iterator_traits<iterator>::difference_type;
		using size_type=size_t;
	private:
		const_pointer d;
		size_type s;
	public:
		array_view():d(nullptr),s(0){}
		array_view(const_pointer d, size_type s):d(d),s(s){}
		void reset(const_pointer d, size_type s){
			this->d=d;
			this->s=s;
		}
		size_type size() const{ return(s); }
		size_type max_size() const{ return(s); }
		bool empty() const{ return(s==0); }
		const_pointer data() const{ return(d); }
		const_reference operator[](size_type i) const{ return(*(d+i)); }
		iterator begin() const{ return(d); }
		const_iterator cbegin() const{ return(d); }
		iterator end() const{ return(d+s); }
		const_iterator cend() const{ return(d+s); }
	};
}
	
template<typename Alloc = std::allocator<void> >
class splinetable{
public:
	typedef Alloc allocator_type;
	typedef std::allocator_traits<allocator_type> allocator_traits;
	
	typedef typename allocator_traits::template rebind_traits<uint32_t>::pointer uint32_t_ptr;
	typedef typename allocator_traits::template rebind_traits<uint64_t>::pointer uint64_t_ptr;
	typedef typename allocator_traits::template rebind_traits<float>::pointer float_ptr;
	typedef typename allocator_traits::template rebind_traits<double>::pointer double_ptr;
	typedef typename allocator_traits::template rebind_traits<double_ptr>::pointer double_ptr_ptr;
	typedef typename allocator_traits::template rebind_traits<char>::pointer char_ptr;
	typedef typename allocator_traits::template rebind_traits<char_ptr>::pointer char_ptr_ptr;
	typedef typename allocator_traits::template rebind_traits<char_ptr_ptr>::pointer char_ptr_ptr_ptr;
	
	///Construct an empty splinetable.
	///The resulting object is useful only for calling read_fits, read_fits_mem, or fit.
	explicit splinetable(allocator_type alloc=Alloc()):
	ndim(0),order(NULL),knots(NULL),nknots(NULL),extents(NULL),periods(NULL),
	coefficients(NULL),naxes(NULL),strides(NULL),naux(0),aux(NULL),allocator(alloc)
	{}
	
	///Construct a splinetable from serialized data previously stored in a FITS file.
	///\param filePath the path to the input file
	explicit splinetable(const std::string& filePath, allocator_type alloc=Alloc()):
	ndim(0),order(NULL),knots(NULL),nknots(NULL),extents(NULL),periods(NULL),
	coefficients(NULL),naxes(NULL),strides(NULL),naux(0),aux(NULL),allocator(alloc)
	{
		read_fits(filePath);
	}
	
	splinetable(splinetable&& other):
	ndim(other.ndim),order(std::move(other.order)),knots(std::move(other.knots)),
	nknots(other.nknots),extents(std::move(other.extents)),
	periods(std::move(other.periods)),coefficients(std::move(other.coefficients)),
	naxes(std::move(other.naxes)),strides(std::move(other.strides)),
	naux(other.naux),aux(std::move(other.aux)),
	allocator(std::move(other.allocator))
	{
		other.ndim=0;
		other.order=NULL;
		other.knots=NULL;
		other.nknots=NULL;
		other.extents=NULL;
		other.periods=NULL;
		other.coefficients=NULL;
		other.naxes=NULL;
		other.strides=NULL;
		other.naux=0;
		other.aux=NULL;
		other.allocator=Alloc();
	}
	
	~splinetable(){
		if(ndim){
			uint64_t ncoeffs=strides[0]*naxes[0];
			for(uint32_t i=0; i<ndim; i++)
				deallocate(knots[i]-order[i],nknots[i]+2*order[i]);
			deallocate(knots,ndim);
			deallocate(nknots,ndim);
			deallocate(order,ndim);
			if(extents){
				deallocate(extents[0],2*ndim);
				deallocate(extents,ndim);
			}
			if(periods)
				deallocate(periods,ndim);
			deallocate(coefficients,ncoeffs);
			deallocate(naxes,ndim);
			deallocate(strides,ndim);
			for(uint32_t i=0; i<naux; i++){
				deallocate(aux[i][0],strlen(&aux[i][0][0])+1);
				deallocate(aux[i][1],strlen(&aux[i][1][0])+1);
				deallocate(aux[i],2);
			}
			deallocate(aux,naux);
		}
	}
	
	splinetable& operator=(splinetable&& other){
		if(&other==this)
			return(*this);
		using std::swap;
		swap(ndim,other.ndim);
		swap(order,other.order);
		swap(knots,other.knots);
		swap(nknots,other.nknots);
		swap(extents,other.extents);
		swap(periods,other.periods);
		swap(coefficients,other.coefficients);
		swap(naxes,other.naxes);
		swap(strides,other.strides);
		swap(naux,other.naux);
		swap(aux,other.aux);
		swap(allocator,other.allocator);
		return(*this);
	}
	
	///Compare splines for equality.
	///Splines are considered equal if evaluation would return the same result
	///at all possible coordinate values. Auxiliary FITS header entries are not
	///considered.
	bool operator==(const splinetable& other) const{
		if (ndim != other.ndim)
			return false;
		if (!std::equal(order,order+ndim,other.order))
			return false;
		if (!std::equal(naxes,naxes+ndim,other.naxes))
			return false;
		if (!std::equal(nknots,nknots+ndim,other.nknots))
			return false;
		for (uint32_t i=0; i<ndim; i++)
			if (!std::equal(knots[i],knots[i]+nknots[i],other.knots[i]))
				return false;
		if (get_ncoeffs() != other.get_ncoeffs())
			return false;
		if (!std::equal(coefficients,coefficients+get_ncoeffs(),other.coefficients))
			return false;
	}

	///Estimate the memory needed to load an existing spline.
	///This function is intended for users using specialized allocators for which
	///it is useful to know the total memory required before constructing the
	///allocator and splinetable. If a convolution will be applied to the spline
	///after it is loaded, doing so will increase its memory requirements, which
	///this function can also take into account.
	///\param filePath the path to the input FITS file
	///\param n_convolution_knots the number of knots possessed by the other
	///       spline with which one dimension of this spline will be convolved.
	///       A value of one is equivalent to no convolution.
	///\param convolution_dimension the dimension of this spline which will have
	///       a convolution applied.
	static size_t estimateMemory(const std::string& filePath,
	                             uint32_t n_convolution_knots = 1,
	                             uint32_t convolution_dimension = 0);
	
	///Read from a FITS file
	///\param path the path to the input file
	bool read_fits(const std::string& path);
	
	///Read from a FITS 'file' in a memory buffer
	///\param the input data buffer
	///\param buffer_size the length of the input buffer
	bool read_fits_mem(void* buffer, size_t buffer_size);
	
	///Write to a FITS file
	///\param path the path to the output file
	void write_fits(const std::string& path) const;
	
	///Write to a FITS memory 'file'.
{'fraction_non_alphanumeric': 0.13137110063071766, 'fraction_numerical': 0.007898176032729132, 'mean_word_length': 5.0109289617486334, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2866320', 'n_tokens_mistral': 6733, 'n_tokens_neox': 6087, 'n_words': 1473}