Ice Sheet System Model  4.18
Code documentation
Functions
pKrigingx.cpp File Reference
#include "./Krigingx.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../../classes/classes.h"
#include "../../shared/io/io.h"

Go to the source code of this file.

Functions

int pKrigingx (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)
 
void ProcessVariogram2 (Variogram **pvariogram, Options *options)
 

Function Documentation

◆ pKrigingx()

int pKrigingx ( 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 at line 11 of file pKrigingx.cpp.

11  {/*{{{*/
12 
13 #ifdef _HAVE_MPI_
14  int num_procs;
15  int my_rank;
16 
17  /*output*/
18  double *predictions = NULL;
19  double *error = NULL;
20 
21  /*Intermediaries*/
22  int mindata,maxdata;
23  double dmindata,dmaxdata;
24  double radius;
25  char *output = NULL;
26  Variogram *variogram = NULL;
27  Observations *observations = NULL;
28 
29  /*timing*/
30  double start, finish;
31  double start_core, finish_core;
32  double start_init, finish_init;
33 
34  /*Get my_rank: */
35  my_rank=IssmComm::GetRank();
36  num_procs=IssmComm::GetSize();
37 
38  /*Get some Options*/
39  ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); start=ISSM_MPI_Wtime();
40  options->Get(&radius,"searchradius",0.);
41 
42  options->Get(&dmindata,"mindata",1.); mindata=int(dmindata);//FIXME (Options come as double but we want to retrive integers)
43  options->Get(&dmaxdata,"maxdata",50.); maxdata=int(dmaxdata);//FIXME (Options come as double but we want to retrive integers)
44 
45  /*Process observation dataset*/
46  ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); start_init=ISSM_MPI_Wtime();
47  observations=new Observations(obs_list,obs_x,obs_y,obs_length,options);
48  ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); finish_init=ISSM_MPI_Wtime();
49 
50  /*Allocate output*/
51  predictions =xNewZeroInit<double>(n_interp);
52  error =xNewZeroInit<double>(n_interp);
53 
54  /*Get output*/
55  options->Get(&output,"output",(char*)"prediction");
56 
57  ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); start_core=ISSM_MPI_Wtime( );
58  if(strcmp(output,"quadtree")==0){
59  observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
60  }
61  else if(strcmp(output,"variomap")==0){
62  observations->Variomap(predictions,x_interp,n_interp);
63  }
64  else if(strcmp(output,"prediction")==0){
65 
66  /*Process Variogram*/
67  ProcessVariogram2(&variogram,options);
68 
69  /*partition loop across threads: */
70  for(int idx=my_rank;idx<n_interp;idx+=num_procs){
71  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<double(idx)/double(n_interp)*100.<<"% \n");
72  observations->InterpolationKriging(&predictions[idx],&error[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,variogram);
73  }
74  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<100.<<"% \n");
75 
76  double *sumpredictions =xNew<double>(n_interp);
77  double *sumerror =xNew<double>(n_interp);
78  ISSM_MPI_Allreduce(predictions,sumpredictions,n_interp,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
79  ISSM_MPI_Allreduce(error,sumerror,n_interp,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
80  xDelete<double>(error); error=sumerror;
81  xDelete<double>(predictions); predictions=sumpredictions;
82  }
83  else if(strcmp(output,"v4")==0){
84 
85  /*partition loop across threads: */
86  for(int idx=my_rank;idx<n_interp;idx+=num_procs){
87  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<double(idx)/double(n_interp)*100.<<"% \n");
88  observations->InterpolationV4(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata);
89  }
90  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<100.<<"% \n");
91 
92  double *sumpredictions =xNew<double>(n_interp);
93  ISSM_MPI_Allreduce(predictions,sumpredictions,n_interp,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
94  xDelete<double>(predictions); predictions=sumpredictions;
95  }
96  else if(strcmp(output,"nearestneighbor")==0){
97 
98  /*partition loop across threads: */
99  for(int idx=my_rank;idx<n_interp;idx+=num_procs){
100  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<double(idx)/double(n_interp)*100.<<"% \n");
101  observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
102  }
103  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<100.<<"% \n");
104 
105  double *sumpredictions =xNew<double>(n_interp);
106  ISSM_MPI_Allreduce(predictions,sumpredictions,n_interp,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
107  xDelete<double>(predictions); predictions=sumpredictions;
108  }
109  else if(strcmp(output,"distance")==0){
110 
111  /*partition loop across threads: */
112  for(int idx=my_rank;idx<n_interp;idx+=num_procs){
113  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<double(idx)/double(n_interp)*100.<<"% \n");
114  observations->Distances(&predictions[idx],&x_interp[idx],&y_interp[idx],1,radius);
115  }
116  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<100.<<"% \n");
117 
118  double *sumpredictions =xNew<double>(n_interp);
119  ISSM_MPI_Allreduce(predictions,sumpredictions,n_interp,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
120  xDelete<double>(predictions); predictions=sumpredictions;
121  }
122  else if(strcmp(output,"idw")==0){
123  double power;
124  options->Get(&power,"power",2.);
125 
126  /*partition loop across threads: */
127  for(int idx=my_rank;idx<n_interp;idx+=num_procs){
128  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<double(idx)/double(n_interp)*100.<<"% \n");
129  observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
130  }
131  _printf0_(" interpolation progress: "<<fixed<<setw(6)<<setprecision(4)<<100.<<"% \n");
132 
133  double *sumpredictions =xNew<double>(n_interp);
134  ISSM_MPI_Allreduce(predictions,sumpredictions,n_interp,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
135  xDelete<double>(predictions); predictions=sumpredictions;
136  }
137  else{
138  _error_("output '" << output << "' not supported yet");
139  }
140  ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); finish_core=ISSM_MPI_Wtime( );
141 
142  /*clean-up and Assign output pointer*/
143  delete variogram;
144  delete observations;
145  xDelete<char>(output);
146  *ppredictions = predictions;
147  *perror = error;
148 
149  ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); finish=ISSM_MPI_Wtime( );
150  _printf0_("\n " << setw(34) << left << "Observation fitering elapsed time: " << finish_init-start_init << " seconds \n\n");
151  _printf0_(" " << setw(34) << left << "Kriging prediction elapsed time: " << finish_core-start_core << " seconds \n\n");
152  _printf0_("\n " << "Total elapsed time " << int((finish-start)/3600) << " hrs " << int(int(finish-start)%3600/60) << " min " << int(finish-start)%60 << " sec\n\n\n");
153  return 1;
154 #else
155  _error_("MPI not available");
156 #endif
157 }/*}}}*/

◆ ProcessVariogram2()

void ProcessVariogram2 ( Variogram **  pvariogram,
Options *  options 
)

Definition at line 158 of file pKrigingx.cpp.

158  {/*{{{*/
159 
160  /*Intermediaries*/
161  Variogram* variogram = NULL;
162  char *model = NULL;
163 
164  if(options->GetOption("model")){
165  options->Get(&model,"model");
166  if (strcmp(model,"gaussian")==0) variogram = new GaussianVariogram(options);
167  else if(strcmp(model,"exponential")==0) variogram = new ExponentialVariogram(options);
168  else if(strcmp(model,"spherical")==0) variogram = new SphericalVariogram(options);
169  else if(strcmp(model,"power")==0) variogram = new PowerVariogram(options);
170  else _error_("variogram " << model << " not supported yet (list of supported variogram: gaussian, exponential, spherical and power)");
171  }
172  else variogram = new GaussianVariogram(options);
173 
174  /*Assign output pointer*/
175  xDelete<char>(model);
176  *pvariogram = variogram;
177 }/*}}}*/
Observations::Distances
void Distances(IssmPDouble *distances, IssmPDouble *x, IssmPDouble *y, int n, IssmPDouble radius)
Definition: Observations.cpp:302
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
ISSM_MPI_Allreduce
int ISSM_MPI_Allreduce(void *sendbuf, void *recvbuf, int count, ISSM_MPI_Datatype datatype, ISSM_MPI_Op op, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:116
ISSM_MPI_SUM
#define ISSM_MPI_SUM
Definition: issmmpi.h:134
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
GaussianVariogram
Definition: GaussianVariogram.h:12
ProcessVariogram2
void ProcessVariogram2(Variogram **pvariogram, Options *options)
Definition: pKrigingx.cpp:158
ISSM_MPI_COMM_WORLD
#define ISSM_MPI_COMM_WORLD
Definition: issmmpi.h:137
Observations::InterpolationIDW
void InterpolationIDW(IssmDouble *pprediction, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int mindata, int maxdata, IssmDouble power)
Definition: Observations.cpp:481
PowerVariogram
Definition: PowerVariogram.h:11
Observations::QuadtreeColoring
void QuadtreeColoring(IssmDouble *A, IssmDouble *x, IssmDouble *y, int n)
Definition: Observations.cpp:692
Observations::Variomap
void Variomap(IssmDouble *gamma, IssmDouble *x, int n)
Definition: Observations.cpp:704
IssmComm::GetSize
static int GetSize(void)
Definition: IssmComm.cpp:46
ISSM_MPI_PDOUBLE
#define ISSM_MPI_PDOUBLE
Definition: issmmpi.h:126
Options::Get
void Get(OptionType *pvalue, const char *name)
Definition: Options.h:21
IssmComm::GetRank
static int GetRank(void)
Definition: IssmComm.cpp:34
Options::GetOption
Option * GetOption(const char *name)
Definition: Options.cpp:67
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Observations::InterpolationKriging
void InterpolationKriging(IssmDouble *pprediction, IssmDouble *perror, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int mindata, int maxdata, Variogram *variogram)
Definition: Observations.cpp:526
SphericalVariogram
Definition: SphericalVariogram.h:11
ISSM_MPI_Wtime
double ISSM_MPI_Wtime(void)
Definition: issmmpi.cpp:511
Variogram
Definition: Variogram.h:10
ExponentialVariogram
Definition: ExponentialVariogram.h:11
ISSM_MPI_Barrier
int ISSM_MPI_Barrier(ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:148
Observations::InterpolationNearestNeighbor
void InterpolationNearestNeighbor(IssmDouble *pprediction, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius)
Definition: Observations.cpp:599
Observations::InterpolationV4
void InterpolationV4(IssmDouble *pprediction, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int mindata, int maxdata)
Definition: Observations.cpp:610
Observations
Declaration of Observations class.
Definition: Observations.h:16
Generated on Thu Jul 2 2020 08:09:24 for Ice Sheet System Model by   doxygen 1.8.19