Changeset 19156

Timestamp:

02/25/15 16:51:50 (10 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixes covertree by computing maxdist on the fly

File:

• issm/trunk-jpl/src/c/classes/kriging/Observations.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/classes/kriging/Observations.cpp

@@ -155 +155 @@
 
         /*Intermediaries*/
-    int          maxdepth = 20;
          IssmPDouble  minspacing,mintrimming,maxtrimming;
 
@@ -166 +165 @@
         options->Get(&minspacing,"minspacing",0.01);
         if(minspacing<=0) _error_("minspacing must > 0");
+
+        /*Get maximum distance between 2 points
+         *  maxDist should be the maximum distance that any two points
+         *  can have between each other. IE p.distance(q) < maxDist for all
+         *  p,q that you will ever try to insert. The cover tree may be invalid
+         *  if an inaccurate maxDist is given.*/
+        IssmPDouble xmin = x[0];
+        IssmPDouble xmax = x[0];
+        IssmPDouble ymin = y[0];
+        IssmPDouble ymax = y[0];
+        for(int i=1;i<n;i++){
+                if(x[i]<xmin) xmin=x[i];
+                if(x[i]>xmax) xmax=x[i];
+                if(y[i]<ymin) ymin=y[i];
+                if(y[i]>ymax) ymax=y[i];
+        }
+        IssmPDouble maxDist = sqrt(pow(xmax-xmin,2)+pow(ymax-ymin,2));
+        IssmPDouble base    = 2.;
+        int         maxdepth = ceilf(log(maxDist)/log(base));
 
          _printf0_("Generating covertree with a maximum depth " <<  maxdepth <<"... ");
@@ -534 +552 @@
         /*First: Create semivariogram matrix for observations*/
         for(i=0;i<n_obs;i++){
+                //printf("%g %g ==> %g\n",x[i],y[i],sqrt(pow(x[i]-x_interp,2)+pow(y[i]-y_interp,2)));
                 for(j=0;j<=i;j++){
                         A[i*(n_obs+1)+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);