Changeset 16989

Timestamp:

12/03/13 08:33:17 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing some stuff in Kriging

Location:

issm/trunk-jpl/src/c

Files:

• bamg/Mesh.cpp (modified) (1 diff)
• classes/kriging/Observations.cpp (modified) (3 diffs)
• classes/kriging/Quadtree.cpp (modified) (1 diff)
• modules/Krigingx/pKrigingx.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/bamg/Mesh.cpp

@@ -4101 +4101 @@
         detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);
 
-        while(t->det>0) { 
+        while(t->det>0){
 
                 /*Increase counter*/

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

@@ -202 +202 @@
 
         /*If radius is not provided or is 0, return all observations*/
-        if(radius==0) radius=this->quadtree->root->length;
+        if(radius==0.) radius=this->quadtree->root->length;
 
         /*Compute radius square*/
@@ -214 +214 @@
         }
         nobs = 0;
-        for (i=0;i<tempnobs;i++){
+        for(i=0;i<tempnobs;i++){
                 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(tempindices[i]));
                 h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
 
                 if(nobs==maxdata && h2>radius2) continue;
-                if(nobs<=maxdata){
+                if(nobs<maxdata){
                         indices[nobs]   = tempindices[i];
                         dists[nobs]     = h2;
@@ -256 +256 @@
 
                 /*Loop over all observations and fill in x, y and obs*/
-                for (i=0;i<nobs;i++){
+                for(i=0;i<nobs;i++){
                         observation=dynamic_cast<Observation*>(this->GetObjectByOffset(indices[i]));
                         observation->WriteXYObs(&x[i],&y[i],&obs[i]);

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

@@ -237 +237 @@
 
                 /*We found the closest observation, now average observation (do not change xi and yi to avoid round off errors*/
-                box->obs[index]->x = (box->obs[index]->weight*box->obs[index]->x + x)/(box->obs[index]->weight+1);
-                box->obs[index]->y = (box->obs[index]->weight*box->obs[index]->y + y)/(box->obs[index]->weight+1);
-                box->obs[index]->xi= int((box->obs[index]->weight*double(box->obs[index]->xi) + double(xi))/(box->obs[index]->weight+1));
-                box->obs[index]->yi= int((box->obs[index]->weight*double(box->obs[index]->yi) + double(yi))/(box->obs[index]->weight+1));
-                box->obs[index]->value   = (box->obs[index]->weight*box->obs[index]->value + value)/(box->obs[index]->weight+1);
-                box->obs[index]->weight += 1;
+                box->obs[index]->x = (box->obs[index]->weight*box->obs[index]->x + x)/(box->obs[index]->weight+1.);
+                box->obs[index]->y = (box->obs[index]->weight*box->obs[index]->y + y)/(box->obs[index]->weight+1.);
+                box->obs[index]->xi= int((box->obs[index]->weight*double(box->obs[index]->xi) + double(xi))/(box->obs[index]->weight+1.));
+                box->obs[index]->yi= int((box->obs[index]->weight*double(box->obs[index]->yi) + double(yi))/(box->obs[index]->weight+1.));
+                box->obs[index]->value   = (box->obs[index]->weight*box->obs[index]->value + value)/(box->obs[index]->weight+1.);
+                box->obs[index]->weight += 1.;
         }
         else{

issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp

@@ -22 +22 @@
         /*Intermediaries*/
         int           mindata,maxdata;
+        double        dmindata,dmaxdata;
         double        radius;
         char         *output       = NULL;
@@ -39 +40 @@
         ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); start=ISSM_MPI_Wtime();
         options->Get(&radius,"searchradius",0.);
-        options->Get(&mindata,"mindata",1);
-        options->Get(&maxdata,"maxdata",50);
+
+        options->Get(&dmindata,"mindata",1.);  mindata=int(dmindata);//FIXME (Options come as double but we want to retrive integers)
+        options->Get(&dmaxdata,"maxdata",50.); maxdata=int(dmaxdata);//FIXME (Options come as double but we want to retrive integers)
 
         /*Process observation dataset*/