Ice Sheet System Model  4.18
Code documentation
Krigingx.cpp
Go to the documentation of this file.
1 
5 #include "./Krigingx.h"
6 #include "../../shared/shared.h"
7 #include "../../toolkits/toolkits.h"
8 #include "../../classes/classes.h"
9 #include "../modules.h"
10 #ifdef _HAVE_GSL_
11 #include <gsl/gsl_linalg.h>
12 #endif
13 int Krigingx(double** ppredictions,double **perror,double* obs_x, double* obs_y, double* obs_list, int obs_length,double* x_interp,double* y_interp,int n_interp,Options* options){/*{{{*/
14 
15  /*output*/
16  double *predictions = NULL;
17  double *error = NULL;
18 
19  /*Intermediaries*/
20  int mindata,maxdata;
21  double dmindata,dmaxdata,dnumthreads; //FIXME (Options come as double but we want to retrive integers)
22  double radius;
23  char *output = NULL;
24  Variogram *variogram = NULL;
25  Observations *observations = NULL;
26 
27  /*threading: */
28  KrigingxThreadStruct gate;
29  int num = _NUMTHREADS_;
30 
31  /*Get Variogram from Options*/
32  ProcessVariogram(&variogram,options);
33  options->Get(&radius,"searchradius",0.);
34  options->Get(&dmindata,"mindata",1.); mindata=int(dmindata);//FIXME (Options come as double but we want to retrive integers)
35  options->Get(&dmaxdata,"maxdata",50.); maxdata=int(dmaxdata);//FIXME (Options come as double but we want to retrive integers)
36  options->Get(&dnumthreads,"numthreads",double(num)); num=int(dnumthreads);//FIXME (Options come as double but we want to retrive integers)
37 
38  /*Process observation dataset*/
39  observations=new Observations(obs_list,obs_x,obs_y,obs_length,options);
40 
41  /*Allocate output*/
42  predictions =xNewZeroInit<double>(n_interp);
43  error =xNewZeroInit<double>(n_interp);
44 
45  /*Get output*/
46  options->Get(&output,"output",(char*)"prediction");
47 
48  if(strcmp(output,"quadtree")==0){
49  observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
50  }
51  else if(strcmp(output,"variomap")==0){
52  observations->Variomap(predictions,x_interp,n_interp);
53  }
54  else if(strcmp(output,"distance")==0){
55  /*initialize thread parameters: */
56  gate.n_interp = n_interp;
57  gate.x_interp = x_interp;
58  gate.y_interp = y_interp;
59  gate.radius = radius;
60  gate.mindata = mindata;
61  gate.maxdata = maxdata;
62  gate.variogram = variogram;
63  gate.observations = observations;
64  gate.predictions = predictions;
65  gate.error = error;
66  gate.numdone = xNewZeroInit<int>(num);
67 
68  /*launch the thread manager with Krigingxt as a core: */
69  LaunchThread(Distancest,(void*)&gate,num);
70  xDelete<int>(gate.numdone);
71  }
72  else if(strcmp(output,"delaunay")==0){
73 
74  #ifdef _HAVE_BAMG_
75  int nobs,nel;
76  double *x = NULL;
77  double *y = NULL;
78  double *data = NULL;
79  int *index = NULL;
80 
81  observations->ObservationList(&x,&y,&data,&nobs);
82 
83  _printf_("Generation Delaunay Triangulation\n");
84  BamgTriangulatex(&index,&nel,x,y,nobs);
85 
86  _printf_("Interpolating\n");
87  xDelete<double>(predictions);
88  InterpFromMeshToMesh2dx(&predictions,index,x,y,nobs,nel,data,nobs,1,x_interp,y_interp,n_interp,options);
89  xDelete<double>(x);
90  xDelete<double>(y);
91  xDelete<double>(data);
92  xDelete<int>(index);
93  #else
94  _error_("you did not compile ISSM with bamg");
95  #endif
96  }
97  else if(strcmp(output,"nearestneighbor")==0){
98  /*initialize thread parameters: */
99  gate.n_interp = n_interp;
100  gate.x_interp = x_interp;
101  gate.y_interp = y_interp;
102  gate.radius = radius;
103  gate.mindata = mindata;
104  gate.maxdata = maxdata;
105  gate.variogram = variogram;
106  gate.observations = observations;
107  gate.predictions = predictions;
108  gate.error = error;
109  gate.numdone = xNewZeroInit<int>(num);
110 
111  /*launch the thread manager with Krigingxt as a core: */
112  LaunchThread(NearestNeighbort,(void*)&gate,num);
113  _printf_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"% \n");
114  xDelete<int>(gate.numdone);
115  }
116  else if(strcmp(output,"idw")==0){ //Inverse distance weighting
117  double power;
118  options->Get(&power,"power",2.);
119  /*initialize thread parameters: */
120  gate.n_interp = n_interp;
121  gate.x_interp = x_interp;
122  gate.y_interp = y_interp;
123  gate.radius = radius;
124  gate.mindata = mindata;
125  gate.maxdata = maxdata;
126  gate.variogram = variogram;
127  gate.observations = observations;
128  gate.predictions = predictions;
129  gate.error = error;
130  gate.numdone = xNewZeroInit<int>(num);
131  gate.power = power;
132 
133  /*launch the thread manager with Krigingxt as a core: */
134  LaunchThread(idwt,(void*)&gate,num);
135  _printf_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"% \n");
136  xDelete<int>(gate.numdone);
137  }
138  else if(strcmp(output,"v4")==0){ //Inverse distance weighting
139 #if !defined(_HAVE_GSL_)
140  _error_("GSL is required for v4 interpolation");
141 #endif
142  /*initialize thread parameters: */
143  gate.n_interp = n_interp;
144  gate.x_interp = x_interp;
145  gate.y_interp = y_interp;
146  gate.radius = radius;
147  gate.mindata = mindata;
148  gate.maxdata = maxdata;
149  gate.variogram = variogram;
150  gate.observations = observations;
151  gate.predictions = predictions;
152  gate.error = error;
153  gate.numdone = xNewZeroInit<int>(num);
154 
155  /*launch the thread manager with Krigingxt as a core: */
156  LaunchThread(v4t,(void*)&gate,num);
157  _printf_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"% \n");
158  xDelete<int>(gate.numdone);
159  }
160  else if(strcmp(output,"prediction")==0){
161 #if !defined(_HAVE_GSL_)
162  _error_("GSL is required for v4 interpolation");
163 #endif
164 
165  /*initialize thread parameters: */
166  gate.n_interp = n_interp;
167  gate.x_interp = x_interp;
168  gate.y_interp = y_interp;
169  gate.radius = radius;
170  gate.mindata = mindata;
171  gate.maxdata = maxdata;
172  gate.variogram = variogram;
173  gate.observations = observations;
174  gate.predictions = predictions;
175  gate.error = error;
176  gate.numdone = xNewZeroInit<int>(num);
177 
178  /*launch the thread manager with Krigingxt as a core: */
179  LaunchThread(Krigingxt,(void*)&gate,num);
180  _printf_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"% \n");
181  xDelete<int>(gate.numdone);
182  }
183  else{
184  _error_("output '" << output << "' not supported yet");
185  }
186 
187  /*clean-up and Assign output pointer*/
188  delete variogram;
189  delete observations;
190  xDelete<char>(output);
191  *ppredictions = predictions;
192  *perror = error;
193  return 1;
194 }/*}}}*/
195 void* Krigingxt(void* vpthread_handle){/*{{{*/
196 
197  /*gate variables :*/
198  KrigingxThreadStruct *gate = NULL;
199  pthread_handle *handle = NULL;
200  int my_thread;
201  int num_threads;
202  int i0,i1;
203 
204  /*recover handle and gate: */
205  handle = (pthread_handle*)vpthread_handle;
206  gate = (KrigingxThreadStruct*)handle->gate;
207  my_thread = handle->id;
208  num_threads = handle->num;
209 
210  /*recover parameters :*/
211  int n_interp = gate->n_interp;
212  double *x_interp = gate->x_interp;
213  double *y_interp = gate->y_interp;
214  double radius = gate->radius;
215  int mindata = gate->mindata;
216  int maxdata = gate->maxdata;
217  Variogram *variogram = gate->variogram;
218  Observations *observations = gate->observations;
219  double *predictions = gate->predictions;
220  double *error = gate->error;
221  int *numdone = gate->numdone;
222 
223  /*partition loop across threads: */
224  PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
225  for(int idx=i0;idx<i1;idx++){
226 
227  /*Print info*/
228  numdone[my_thread]=idx-i0;
229  if(my_thread==0){
230  int alldone=numdone[0];
231  for(int i=1;i<num_threads;i++) alldone+=numdone[i];
232  _printf_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"% ");
233  }
234 
235  /*Kriging interpolation*/
236  observations->InterpolationKriging(&predictions[idx],&error[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,variogram);
237  }
238 
239  return NULL;
240 }/*}}}*/
241 void* NearestNeighbort(void* vpthread_handle){/*{{{*/
242 
243  /*gate variables :*/
244  KrigingxThreadStruct *gate = NULL;
245  pthread_handle *handle = NULL;
246  int my_thread;
247  int num_threads;
248  int i0,i1;
249 
250  /*recover handle and gate: */
251  handle = (pthread_handle*)vpthread_handle;
252  gate = (KrigingxThreadStruct*)handle->gate;
253  my_thread = handle->id;
254  num_threads = handle->num;
255 
256  /*recover parameters :*/
257  int n_interp = gate->n_interp;
258  double *x_interp = gate->x_interp;
259  double *y_interp = gate->y_interp;
260  double radius = gate->radius;
261  int mindata = gate->mindata;
262  int maxdata = gate->maxdata;
263  Variogram *variogram = gate->variogram;
264  Observations *observations = gate->observations;
265  double *predictions = gate->predictions;
266  double *error = gate->error;
267  int *numdone = gate->numdone;
268 
269  /*partition loop across threads: */
270  PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
271  for(int idx=i0;idx<i1;idx++){
272 
273  /*Print info*/
274  numdone[my_thread]=idx-i0;
275  if(my_thread==0){
276  int alldone=numdone[0];
277  for(int i=1;i<num_threads;i++) alldone+=numdone[i];
278  _printf_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"% ");
279  }
280 
281  observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
282  }
283 
284  return NULL;
285 }/*}}}*/
286 void* idwt(void* vpthread_handle){/*{{{*/
287 
288  /*gate variables :*/
289  KrigingxThreadStruct *gate = NULL;
290  pthread_handle *handle = NULL;
291  int my_thread;
292  int num_threads;
293  int i0,i1;
294 
295  /*recover handle and gate: */
296  handle = (pthread_handle*)vpthread_handle;
297  gate = (KrigingxThreadStruct*)handle->gate;
298  my_thread = handle->id;
299  num_threads = handle->num;
300 
301  /*recover parameters :*/
302  int n_interp = gate->n_interp;
303  double *x_interp = gate->x_interp;
304  double *y_interp = gate->y_interp;
305  double radius = gate->radius;
306  int mindata = gate->mindata;
307  int maxdata = gate->maxdata;
308  Variogram *variogram = gate->variogram;
309  Observations *observations = gate->observations;
310  double *predictions = gate->predictions;
311  double *error = gate->error;
312  int *numdone = gate->numdone;
313  double power = gate->power;
314 
315  /*partition loop across threads: */
316  PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
317  for(int idx=i0;idx<i1;idx++){
318 
319  /*Print info*/
320  numdone[my_thread]=idx-i0;
321  if(my_thread==0){
322  int alldone=numdone[0];
323  for(int i=1;i<num_threads;i++) alldone+=numdone[i];
324  _printf_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"% ");
325  }
326 
327  observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
328  }
329  return NULL;
330 }/*}}}*/
331 void* v4t(void* vpthread_handle){/*{{{*/
332 
333  /*gate variables :*/
334  KrigingxThreadStruct *gate = NULL;
335  pthread_handle *handle = NULL;
336  int my_thread;
337  int num_threads;
338  int i0,i1;
339 
340  /*recover handle and gate: */
341  handle = (pthread_handle*)vpthread_handle;
342  gate = (KrigingxThreadStruct*)handle->gate;
343  my_thread = handle->id;
344  num_threads = handle->num;
345 
346  /*recover parameters :*/
347  int n_interp = gate->n_interp;
348  double *x_interp = gate->x_interp;
349  double *y_interp = gate->y_interp;
350  double radius = gate->radius;
351  int mindata = gate->mindata;
352  int maxdata = gate->maxdata;
353  Variogram *variogram = gate->variogram;
354  Observations *observations = gate->observations;
355  double *predictions = gate->predictions;
356  double *error = gate->error;
357  int *numdone = gate->numdone;
358 
359  /*partition loop across threads: */
360  PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
361  for(int idx=i0;idx<i1;idx++){
362 
363  /*Print info*/
364  numdone[my_thread]=idx-i0;
365  if(my_thread==0){
366  int alldone=numdone[0];
367  for(int i=1;i<num_threads;i++) alldone+=numdone[i];
368  _printf_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"% ");
369  }
370 
371  observations->InterpolationV4(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata);
372  }
373  return NULL;
374 }/*}}}*/
375 void* Distancest(void* vpthread_handle){/*{{{*/
376 
377  /*gate variables :*/
378  KrigingxThreadStruct *gate = NULL;
379  pthread_handle *handle = NULL;
380  int my_thread;
381  int num_threads;
382  int i0,i1;
383 
384  /*recover handle and gate: */
385  handle = (pthread_handle*)vpthread_handle;
386  gate = (KrigingxThreadStruct*)handle->gate;
387  my_thread = handle->id;
388  num_threads = handle->num;
389 
390  /*recover parameters :*/
391  int n_interp = gate->n_interp;
392  double *x_interp = gate->x_interp;
393  double *y_interp = gate->y_interp;
394  double radius = gate->radius;
395  int mindata = gate->mindata;
396  int maxdata = gate->maxdata;
397  Variogram *variogram = gate->variogram;
398  Observations *observations = gate->observations;
399  double *predictions = gate->predictions;
400  double *error = gate->error;
401  int *numdone = gate->numdone;
402 
403  /*partition loop across threads: */
404  PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
405  observations->Distances(&predictions[i0],&x_interp[i0],&y_interp[i0],i1-i0,radius);
406  return NULL;
407 }/*}}}*/
408 
409 void ProcessVariogram(Variogram **pvariogram,Options* options){/*{{{*/
410 
411  /*Intermediaries*/
412  Variogram* variogram = NULL;
413  char *model = NULL;
414 
415  if(options->GetOption("model")){
416  options->Get(&model,"model");
417  if (strcmp(model,"gaussian")==0) variogram = new GaussianVariogram(options);
418  else if(strcmp(model,"exponential")==0) variogram = new ExponentialVariogram(options);
419  else if(strcmp(model,"spherical")==0) variogram = new SphericalVariogram(options);
420  else if(strcmp(model,"power")==0) variogram = new PowerVariogram(options);
421  else _error_("variogram " << model << " not supported yet (list of supported variogram: gaussian, exponential, spherical and power)");
422  }
423  else variogram = new GaussianVariogram(options);
424 
425  /*Assign output pointer*/
426  xDelete<char>(model);
427  *pvariogram = variogram;
428 }/*}}}*/
KrigingxThreadStruct::observations
Observations * observations
Definition: Krigingx.h:28
KrigingxThreadStruct::variogram
Variogram * variogram
Definition: Krigingx.h:27
Options
Definition: Options.h:9
pthread_handle::id
int id
Definition: issm_threads.h:12
v4t
void * v4t(void *vpthread_handle)
Definition: Krigingx.cpp:331
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
LaunchThread
void LaunchThread(void *function(void *), void *gate, int num_threads)
Definition: LaunchThread.cpp:25
pthread_handle::num
int num
Definition: issm_threads.h:13
GaussianVariogram
Definition: GaussianVariogram.h:12
KrigingxThreadStruct
Definition: Krigingx.h:20
KrigingxThreadStruct::error
double * error
Definition: Krigingx.h:30
KrigingxThreadStruct::mindata
int mindata
Definition: Krigingx.h:25
KrigingxThreadStruct::y_interp
double * y_interp
Definition: Krigingx.h:23
KrigingxThreadStruct::n_interp
int n_interp
Definition: Krigingx.h:21
Krigingx.h
PowerVariogram
Definition: PowerVariogram.h:11
Observations::QuadtreeColoring
void QuadtreeColoring(IssmDouble *A, IssmDouble *x, IssmDouble *y, int n)
Definition: Observations.cpp:692
BamgTriangulatex
int BamgTriangulatex(int **pindex, int *pnels, double *x, double *y, int nods)
Definition: BamgTriangulatex.cpp:14
Observations::Variomap
void Variomap(IssmDouble *gamma, IssmDouble *x, int n)
Definition: Observations.cpp:704
pthread_handle
Definition: issm_threads.h:10
PartitionRange
void PartitionRange(int *pi0, int *pi1, int num_el, int num_threads, int my_thread)
Definition: PartitionRange.cpp:13
KrigingxThreadStruct::numdone
int * numdone
Definition: Krigingx.h:31
Krigingxt
void * Krigingxt(void *vpthread_handle)
Definition: Krigingx.cpp:195
Options::Get
void Get(OptionType *pvalue, const char *name)
Definition: Options.h:21
KrigingxThreadStruct::radius
double radius
Definition: Krigingx.h:24
Observations::ObservationList
void ObservationList(IssmDouble **px, IssmDouble **py, IssmDouble **pobs, int *pnobs)
Definition: Observations.cpp:316
ProcessVariogram
void ProcessVariogram(Variogram **pvariogram, Options *options)
Definition: Krigingx.cpp:409
idwt
void * idwt(void *vpthread_handle)
Definition: Krigingx.cpp:286
Krigingx
int Krigingx(double **ppredictions, double **perror, double *obs_x, double *obs_y, double *obs_list, int obs_length, double *x_interp, double *y_interp, int n_interp, Options *options)
Definition: Krigingx.cpp:13
Distancest
void * Distancest(void *vpthread_handle)
Definition: Krigingx.cpp:375
Options::GetOption
Option * GetOption(const char *name)
Definition: Options.cpp:67
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
KrigingxThreadStruct::predictions
double * predictions
Definition: Krigingx.h:29
NearestNeighbort
void * NearestNeighbort(void *vpthread_handle)
Definition: Krigingx.cpp:241
SphericalVariogram
Definition: SphericalVariogram.h:11
Variogram
Definition: Variogram.h:10
ExponentialVariogram
Definition: ExponentialVariogram.h:11
KrigingxThreadStruct::power
double power
Definition: Krigingx.h:32
pthread_handle::gate
void * gate
Definition: issm_threads.h:11
Observations
Declaration of Observations class.
Definition: Observations.h:16
KrigingxThreadStruct::x_interp
double * x_interp
Definition: Krigingx.h:22
InterpFromMeshToMesh2dx
int InterpFromMeshToMesh2dx(double **pdata_interp, int *index_data, double *x_data, double *y_data, int nods_data, int nels_data, double *data, int M_data, int N_data, double *x_interp, double *y_interp, int N_interp, Options *options)
Definition: InterpFromMeshToMesh2dx.cpp:16
KrigingxThreadStruct::maxdata
int maxdata
Definition: Krigingx.h:26
Generated on Thu Jul 2 2020 08:09:19 for Ice Sheet System Model by   doxygen 1.8.19