Changeset 14048

Timestamp:

11/28/12 21:04:08 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: working on minimum curvature interpolation module

Location:

issm/trunk-jpl/src/c

Files:

• Container/Observations.cpp (modified) (1 diff)
• Container/Observations.h (modified) (1 diff)
• modules/Krigingx/Krigingx.cpp (modified) (2 diffs)
• modules/Krigingx/Krigingx.h (modified) (1 diff)

issm/trunk-jpl/src/c/Container/Observations.cpp

@@ -424 +424 @@
         *pprediction = obs;
 }/*}}}*/
+/*FUNCTION Observations::InterpolationV4{{{*/
+void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){
+
+        /*Intermediaries*/
+        int         i,j,n_obs;
+        IssmPDouble prediction,h;
+        IssmPDouble *x       = NULL;
+        IssmPDouble *y       = NULL;
+        IssmPDouble *obs     = NULL;
+        IssmPDouble *Green   = NULL;
+        IssmPDouble *weights = NULL;
+        IssmPDouble *d       = NULL;
+
+        /*Some checks*/
+        _assert_(maxdata>0);
+        _assert_(pprediction);
+
+        /*If radius is not provided or is 0, return all observations*/
+        if(radius==0) radius=this->quadtree->root->length;
+
+        /*Get list of observations for current point*/
+        this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
+
+        /*If we have less observations than mindata, return UNDEF*/
+        if(n_obs<mindata){
+                prediction = UNDEF; 
+        }
+        else{
+
+                /*Allocate intermediary matrix and vectors*/
+                Green = xNew<IssmPDouble>(n_obs*n_obs);
+                d     = xNew<IssmPDouble>(n_obs);
+
+                /*First: distance vector*/
+                for(i=0;i<n_obs;i++){
+                        d[i] = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
+                }
+
+                /*Build Green function matrix*/
+                for(i=0;i<n_obs;i++){
+                        for(j=0;j<=i;j++){
+                                h = sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) );
+                                if(h>0){
+                                        Green[j*n_obs+i] = h*h*(log(h)-1.);
+                                }
+                                else{
+                                        Green[j*n_obs+i] = 0.;
+                                }
+                                Green[i*n_obs+j] = Green[j*n_obs+i];
+                        }
+                }
+
+                /*Compute weights*/
+#if _HAVE_GSL_
+                SolverxSeq(&weights,Green,obs,n_obs); // Green^-1 obs
+#else
+                _error_("GSL is required");
+#endif
+
+                /*ok now we can interpolate with these weights*/
+
+                _error_("Not implemented yet");
+        }
+
+        /*clean-up*/
+        *pprediction = prediction;
+        xDelete<IssmPDouble>(x);
+        xDelete<IssmPDouble>(y);
+        xDelete<IssmPDouble>(obs);
+}/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/
 void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){

issm/trunk-jpl/src/c/Container/Observations.h

@@ -28 +28 @@
                 void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
                 void InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power);
+                void InterpolationV4(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata);
                 void InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram);
                 void InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);

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

@@ -117 +117 @@
                 xDelete<double>(gate.percent);
         }
+        else if(strcmp(output,"v4")==0){ //Inverse distance weighting
+                /*initialize thread parameters: */
+                gate.n_interp     = n_interp;
+                gate.x_interp     = x_interp;
+                gate.y_interp     = y_interp;
+                gate.radius       = radius;
+                gate.mindata      = mindata;
+                gate.maxdata      = maxdata;
+                gate.variogram    = variogram;
+                gate.observations = observations;
+                gate.predictions  = predictions;
+                gate.error        = error;
+                gate.percent      = xNewZeroInit<double>(num);
+
+                /*launch the thread manager with Krigingxt as a core: */
+                LaunchThread(v4t,(void*)&gate,num);
+                _printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
+                xDelete<double>(gate.percent);
+        }
         else if(strcmp(output,"prediction")==0){
 
@@ -292 +311 @@
         return NULL;
 }/*}}}*/
+/*FUNCTION v4t{{{*/
+void* v4t(void* vpthread_handle){
+
+        /*gate variables :*/
+        KrigingxThreadStruct *gate        = NULL;
+        pthread_handle       *handle      = NULL;
+        int my_thread;
+        int num_threads;
+        int i0,i1;
+
+        /*recover handle and gate: */
+        handle      = (pthread_handle*)vpthread_handle;
+        gate        = (KrigingxThreadStruct*)handle->gate;
+        my_thread   = handle->id;
+        num_threads = handle->num;
+
+        /*recover parameters :*/
+        int           n_interp     = gate->n_interp;
+        double       *x_interp     = gate->x_interp;
+        double       *y_interp     = gate->y_interp;
+        double        radius       = gate->radius;
+        int           mindata      = gate->mindata;
+        int           maxdata      = gate->maxdata;
+        Variogram    *variogram    = gate->variogram;
+        Observations *observations = gate->observations;
+        double       *predictions  = gate->predictions;
+        double       *error        = gate->error;
+        double       *percent      = gate->percent;
+
+        /*Intermediaries*/
+        double localpercent;
+
+        /*partition loop across threads: */
+        PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
+        for(int idx=i0;idx<i1;idx++){
+
+                /*Print info*/
+                percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
+                localpercent=percent[0];
+                for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
+                if(my_thread==0) _printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");
+
+                observations->InterpolationV4(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata);
+        }
+        return NULL;
+}/*}}}*/
 
 void ProcessVariogram(Variogram **pvariogram,Options* options){/*{{{*/

issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h

@@ -35 +35 @@
 void* NearestNeighbort(void*);
 void* idwt(void*);
+void* v4t(void*);
 #endif /* _KRIGINGX_H */