added an epsilon to peak_finder

fit/interpolation.h

@@ -276,21 +276,27 @@ T LinInterp<T>::operator()(T x) const
 // ----------------------------------------------------------------------------
 
 template<typename T>
-void find_peaks(std::size_t iLen, const T* px, const T* py, unsigned int iOrder,
+void find_peaks(std::size_t iLen, const T* _px, const T* _py, unsigned int iOrder,
 	std::vector<T>& vecMaximaX, std::vector<T>& vecMaximaSize, std::vector<T>& vecMaximaWidth,
-	std::size_t iNumSpline = 512)
+	T eps = tl::get_epsilon<T>(), std::size_t iNumSpline = 512)
 {
-	BSpline<T> spline(iLen, px, py, iOrder);
+	// allocate memory
+	std::unique_ptr<T, std::default_delete<T[]>> uptrMem{new T[6*iNumSpline]};
+	T *px = uptrMem.get() + 0*iNumSpline;
+	T *py = uptrMem.get() + 1*iNumSpline;
+	T *pSplineX = uptrMem.get() + 2*iNumSpline;
+	T *pSplineY = uptrMem.get() + 3*iNumSpline;
+	T *pSplineDiff = uptrMem.get() + 4*iNumSpline;
+	T *pSplineDiff2 = uptrMem.get() + 5*iNumSpline;
 
 
-	// allocate memory
-	std::unique_ptr<T, std::default_delete<T[]>> uptrMem{new T[4*iNumSpline]};
-	T *pSplineX = uptrMem.get();
-	T *pSplineY = uptrMem.get() + iNumSpline;
-	T *pSplineDiff = uptrMem.get() + 2*iNumSpline;
-	T *pSplineDiff2 = uptrMem.get() + 3*iNumSpline;
+	// sort input values
+	std::copy(_px, _px+iLen, px);
+	std::copy(_py, _py+iLen, py);
+	tl::sort_2<T*>(px, px+iLen, py);
 
 
+	BSpline<T> spline(iLen, px, py, iOrder);
 	const T* pdyMin = std::min_element(py, py+iLen);
 
 	for(std::size_t iSpline=0; iSpline<iNumSpline; ++iSpline)
@@ -304,10 +310,14 @@ void find_peaks(std::size_t iLen, const T* px, const T* py, unsigned int iOrder,
 
 	tl::diff(iNumSpline, pSplineX, pSplineY, pSplineDiff);
 	tl::diff(iNumSpline, pSplineX, pSplineDiff, pSplineDiff2);
-	std::vector<std::size_t> vecZeroes = tl::find_zeroes<T>(iNumSpline, pSplineDiff);
+	std::vector<std::size_t> vecZeroes = tl::find_zeroes<T>(iNumSpline, pSplineDiff, eps);
+
+	//std::ofstream ofDbg{"0.dat"};
+	//for(std::size_t i=0; i<iNumSpline; ++i)
+	//	ofDbg << pSplineX[i] << "\t" << pSplineY[i] << "\t" << pSplineDiff[i] << "\t" << pSplineDiff2[i] << std::endl;
 
 
-	for(std::size_t iZeroIdx = 0; iZeroIdx<vecZeroes.size(); ++iZeroIdx)
+	for(std::size_t iZeroIdx=0; iZeroIdx<vecZeroes.size(); ++iZeroIdx)
 	{
 		const std::size_t iZero = vecZeroes[iZeroIdx];
 
@@ -332,7 +342,7 @@ void find_peaks(std::size_t iLen, const T* px, const T* py, unsigned int iOrder,
 		if(iMinIdxLeft>=0)
 		{
 			dHeight += (pSplineY[iZero]-pSplineY[iMinIdxLeft]);
-			dWidth += fabs((pSplineX[iZero]-pSplineX[iMinIdxLeft]));
+			dWidth += std::abs((pSplineX[iZero]-pSplineX[iMinIdxLeft]));
 			dDiv += 1.;
 		}
 
@@ -340,7 +350,7 @@ void find_peaks(std::size_t iLen, const T* px, const T* py, unsigned int iOrder,
 		if(iMinIdxRight>=0)
 		{
 			dHeight += (pSplineY[iZero]-pSplineY[iMinIdxRight]);
-			dWidth += fabs((pSplineX[iZero]-pSplineX[iMinIdxRight]));
+			dWidth += std::abs((pSplineX[iZero]-pSplineX[iMinIdxRight]));
 			dDiv += 1.;
 		}
 

helper/misc.h

@@ -191,7 +191,7 @@ void sort_2(Iter begin1, Iter end1, Iter begin2)
 		pObj[i].vec.push_back(*(begin2+i));
 	}
 
-	std::sort(pObj, pObj+N, comp_fkt<T>);
+	std::stable_sort(pObj, pObj+N, comp_fkt<T>);
 	for(std::size_t i=0; i<N; ++i)
 	{
 		*(begin1+i) = pObj[i].vec[0];
@@ -219,7 +219,7 @@ void sort_3(Iter begin1, Iter end1, Iter begin2, Iter begin3)
 		pObj[i].vec.push_back(*(begin3+i));
 	}
 
-	std::sort(pObj, pObj+N, comp_fkt<T>);
+	std::stable_sort(pObj, pObj+N, comp_fkt<T>);
 	for(std::size_t i=0; i<N; ++i)
 	{
 		*(begin1+i) = pObj[i].vec[0];
@@ -243,7 +243,7 @@ t_cont<std::size_t> sorted_idx(const t_cont<t_val>& vec, t_func fkt)
 	for(std::size_t i=0; i<vec.size(); ++i)
 		vecIdx.push_back(i);
 
-	std::sort(vecIdx.begin(), vecIdx.end(),
+	std::stable_sort(vecIdx.begin(), vecIdx.end(),
 		[&vec, &fkt] (std::size_t iIdx0, std::size_t iIdx1) -> bool
 	{
 		return fkt(vec[iIdx0], vec[iIdx1]);

math/geo.h

@@ -1215,7 +1215,7 @@ public:
 
 
 template<typename T = double>
-std::vector<std::size_t> find_zeroes(std::size_t N, const T* pIn)
+std::vector<std::size_t> find_zeroes(std::size_t N, const T* pIn, T eps = tl::get_epsilon<T>())
 {
 	using t_vec = ublas::vector<T>;
 	std::vector<std::size_t> vecIndices;
@@ -1234,8 +1234,8 @@ std::vector<std::size_t> find_zeroes(std::size_t N, const T* pIn)
 		pos1[0] = 1.; pos1[1] = pIn[i+1];
 		Line<T> line(pos0, pos1-pos0);
 
-		T param;
-		if(!line.intersect(xaxis, param))
+		T param{};
+		if(!line.intersect(xaxis, param, eps))
 			continue;
 
 		t_vec posInters = line(param);