From 9df2340837fc308b865aa92f117e0e681bc63f6d Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem@mithril.(none)>
Date: Sun, 8 Mar 2009 22:05:15 -0700
Subject: [PATCH] Better error reporting

---
 src/dinterpolate.hpp |  4 +++-
 src/dinterpolate.tcc | 13 ++++++++-----
 2 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/src/dinterpolate.hpp b/src/dinterpolate.hpp
index 2777764..5109ce9 100644
--- a/src/dinterpolate.hpp
+++ b/src/dinterpolate.hpp
@@ -45,6 +45,7 @@ namespace CosmoTool {
      */
     DelaunayInterpolate(CoordType *positions, IType *values, uint32_t *simplex_list,
 			uint32_t numPoints, uint32_t numSimplex)
+	throw (InvalidArgumentException)
     {
       this->positions = positions;
       this->values = values;
@@ -69,7 +70,8 @@ namespace CosmoTool {
       gsl_eigen_symmv_free(eigen_work);
     }
 
-    void buildPreweight();
+    void buildPreweight()
+	throw (InvalidArgumentException);
     void buildQuickAccess();
     void buildHyperplane(const PType *v, CoordType& hyper);
 
diff --git a/src/dinterpolate.tcc b/src/dinterpolate.tcc
index 810e1e7..d545788 100644
--- a/src/dinterpolate.tcc
+++ b/src/dinterpolate.tcc
@@ -70,6 +70,7 @@ namespace CosmoTool {
 
   template<typename PType, typename IType, int N>
   void DelaunayInterpolate<PType,IType,N>::buildPreweight()
+	throw(InvalidArgumentException)
   {
     double preweight[N*N];
     double preweight_inverse[N*N];
@@ -97,6 +98,8 @@ namespace CosmoTool {
 	int signum;
 
 	gsl_linalg_LU_decomp(&M.matrix, p, &signum);
+	if (fabs(gsl_linalg_LU_det(&M.matrix, signum)) < 1e-10)
+	  throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar.");
 	gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix);
 
 	for (int j = 0; j < N*N; j++)
@@ -239,7 +242,7 @@ namespace CosmoTool {
 
 	// The point does not belong to any simplex.
 	if (i != N_ngb)
-	  throw InvalidArgumentException();
+	  throw InvalidArgumentException("the given point does not belong to any simplex");
 
 	N_ngb *= 2;
 	cell_Ngb = new QuickCell *[N_ngb];
@@ -271,13 +274,13 @@ namespace CosmoTool {
 	for (int j = 0; j < N; j++)
 	  weight[i] += preweight[(i-1)*N+j]*(c[j]-p0[j]);
 
-	assert(weight[i] >= 0);
-	assert(weight[i] <= 1);
+	assert(weight[i] > -1e-7);
+	assert(weight[i] < 1+1e-7);
 	sum_weight += weight[i];
       }   
     weight[0] = 1-sum_weight;
-    assert(weight[0] >= 0);
-    assert(weight[0] <= 1);
+    assert(weight[0] > -1e-7);
+    assert(weight[0] < (1+1e-7));
     
     // We compute the final value by weighing the value at the N+1
     // points by the proper weight.