BLI2D.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab - May 30, 2009
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ May 30, 2009 - CA
14  Was first added in v2.5.1.
15  @ Mar 08, 2010 - CA
16  Fix `converting to `int' from `double'' warnings by using TMath::FloorInt
17  in all lines similar to `int ix = (xmax-xmin)/dx'.
18  @ July 29, 2011 - AM
19  Added BLI2DNonUnifGrid.
20 
21 */
22 //____________________________________________________________________________
23 
24 #include <TMath.h>
25 
26 #include <cassert>
27 #include <limits>
28 
31 
32 using namespace genie;
33 
35 
36 //___________________________________________________________________________
38 {
39 
40 }
41 //___________________________________________________________________________
43 {
44  if (fX) { delete [] fX; }
45  if (fY) { delete [] fY; }
46  if (fZ) { delete [] fZ; }
47 }
48 //___________________________________________________________________________
49 int BLI2DGrid::IdxZ(int ix, int iy) const
50 {
51  return ix*fNY+iy;
52 }
53 //___________________________________________________________________________
54 //___________________________________________________________________________
55 //___________________________________________________________________________
57 
59 {
60  this->Init();
61 }
62 //___________________________________________________________________________
64  int nx, double xmin, double xmax, int ny, double ymin, double ymax)
65 {
66  this->Init(nx, xmin, xmax, ny, ymin, ymax);
67 }
68 //___________________________________________________________________________
70  int nx, int ny, double *x, double *y, double *z)
71 {
72  double xmin = x[0];
73  double xmax = x[nx-1];
74  double ymin = y[0];
75  double ymax = y[nx-1];
76 
77  this->Init(nx, xmin, xmax, ny, ymin, ymax);
78 
79  for(int ix=0; ix<nx; ix++) {
80  for(int iy=0; iy<ny; iy++) {
81  this->AddPoint(x[ix], y[iy], z[this->IdxZ(ix,iy)]);
82  }
83  }
84 }
85 //___________________________________________________________________________
86 bool BLI2DUnifGrid::AddPoint(double x, double y, double z)
87 {
88  int ix = TMath::FloorNint( (x - fXmin + fDX/2) / fDX );
89  int iy = TMath::FloorNint( (y - fYmin + fDY/2) / fDY );
90  int iz = this->IdxZ(ix,iy);
91 
92  fZ[iz] = z;
93 
94  fZmin = TMath::Min(z, fZmin);
95  fZmax = TMath::Max(z, fZmax);
96 
97  LOG("BLI2DUnifGrid", pDEBUG)
98  << "Added x = " << x << " (ix = " << ix << ")"
99  << " y = " << y << " (iy = " << iy << ") -> "
100  << " z = " << z << " (iz = " << iz << ")";
101 
102  return true;
103 }
104 //___________________________________________________________________________
105 double BLI2DUnifGrid::Evaluate(double x, double y) const
106 {
107  if(x < fXmin || x > fXmax) return 0.;
108  if(y < fYmin || y > fYmax) return 0.;
109 
110  int ix_lo = TMath::FloorNint( (x - fXmin) / fDX );
111  int iy_lo = TMath::FloorNint( (y - fYmin) / fDY );
112  int ix_hi = ix_lo + 1;
113  int iy_hi = iy_lo + 1;
114 
115  double x1 = fX[ix_lo];
116  double x2 = fX[ix_hi];
117  double y1 = fY[iy_lo];
118  double y2 = fY[iy_hi];
119 
120  double z11 = fZ[ this->IdxZ(ix_lo,iy_lo) ];
121  double z21 = fZ[ this->IdxZ(ix_hi,iy_lo) ];
122  double z12 = fZ[ this->IdxZ(ix_lo,iy_hi) ];
123  double z22 = fZ[ this->IdxZ(ix_hi,iy_hi) ];
124 
125  double z1 = z11 * (x2-x)/(x2-x1) + z21 * (x-x1)/(x2-x1);
126  double z2 = z12 * (x2-x)/(x2-x1) + z22 * (x-x1)/(x2-x1);
127  double z = z1 * (y2-y)/(y2-y1) + z2 * (y-y1)/(y2-y1);
128 
129 /*
130  LOG("BLI2DUnifGrid", pDEBUG) << "x = " << x << " -> nearby nodes: " << x1 << ", " << x2;
131  LOG("BLI2DUnifGrid", pDEBUG) << "y = " << x << " -> nearby nodes: " << y1 << ", " << y2;
132  LOG("BLI2DUnifGrid", pDEBUG) << "z11 := z(" << this->IdxZ(ix_lo,iy_lo) << ") = " << z11;
133  LOG("BLI2DUnifGrid", pDEBUG) << "z21 := z(" << this->IdxZ(ix_hi,iy_lo) << ") = " << z21;
134  LOG("BLI2DUnifGrid", pDEBUG) << "z12 := z(" << this->IdxZ(ix_lo,iy_hi) << ") = " << z12;
135  LOG("BLI2DUnifGrid", pDEBUG) << "z22 := z(" << this->IdxZ(ix_hi,iy_hi) << ") = " << z22;
136  LOG("BLI2DUnifGrid", pDEBUG) << "z1 = " << z1 << ", z2 = " << z2;
137  LOG("BLI2DUnifGrid", pDEBUG) << "interpolated z(x,y) = " << z;
138 */
139 
140  return z;
141 }
142 //___________________________________________________________________________
144  int nx, double xmin, double xmax, int ny, double ymin, double ymax)
145 {
146  fNX = 0;
147  fNY = 0;
148  fNZ = 0;
149  fXmin = 0.;
150  fXmax = 0.;
151  fYmin = 0.;
152  fYmax = 0.;
155  fDX = 0.;
156  fDY = 0.;
157  fX = 0;
158  fY = 0;
159  fZ = 0;
160 
161  if(nx>1 && ny>1) {
162  fNX = nx;
163  fNY = ny;
164  fNZ = nx * ny;
165 
166  fXmin = xmin;
167  fXmax = xmax;
168  fYmin = ymin;
169  fYmax = ymax;
172 
173  fDX = (xmax-xmin)/(nx-1);
174  fDY = (ymax-ymin)/(ny-1);
175 
176  fX = new double[fNX];
177  fY = new double[fNY];
178  fZ = new double[fNZ];
179 
180  for(int i=0; i<fNX; i++) { fX[i] = xmin + i*fDX; }
181  for(int i=0; i<fNY; i++) { fY[i] = ymin + i*fDY; }
182  for(int i=0; i<fNZ; i++) { fZ[i] = 0.; }
183  }
184 }
185 //___________________________________________________________________________
186 //___________________________________________________________________________
187 //___________________________________________________________________________
189 
191 {
192  this->Init();
193 }
194 //___________________________________________________________________________
196  int nx, double xmin, double xmax, int ny, double ymin, double ymax)
197 {
198  this->Init(nx, xmin, xmax, ny, ymin, ymax);
199 }
200 //___________________________________________________________________________
202  int nx, int ny, double *x, double *y, double *z)
203 {
204  double xmin = x[0];
205  double xmax = x[nx-1];
206  double ymin = y[0];
207  double ymax = y[ny-1];
208 
209  this->Init(nx, xmin, xmax, ny, ymin, ymax);
210 
211  for(int ix=0; ix<nx; ix++) {
212  for(int iy=0; iy<ny; iy++) {
213  this->AddPoint(x[ix], y[iy], z[this->IdxZ(ix,iy)]);
214  }
215  }
216 }
217 //___________________________________________________________________________
218 bool BLI2DNonUnifGrid::AddPoint(double x, double y, double z)
219 {
220 
221  // check the x,y values' existence before moving anything
222  // if they do, use them
223  // if they don't, add them in the appropriate place
224  //
225  int xidx = -1;
226  int yidx = -1;
227  bool changex = true;
228  bool changey = true;
229  // first, check and see if the x,y values exist
230  // NOTE: == should not be used with double values:
231  // instead, a tolerance of ~5 significant digits is checked
232  for(int i=0;i<fNFillX;i++)
233  {
234  if (!(TMath::Abs(x-fX[i])<=.5*.0000001*(TMath::Abs(x)+TMath::Abs(fX[i]))) && x<fX[i])
235  {
236  // missed one, make sure there aren't too many x and move everything up
237  if (fNFillX==fNX)
238  {
239  LOG("BLI2DNonUnifGrid", pWARN) << "Warning: too many x values, cannot add x = "<<x;
240  return false;
241  }
242  else
243  {
244  xidx=i;
245  changex=true;
246  break;
247  }
248  }
249  else
250  {
251  if (TMath::Abs(x-fX[i])<=.5*.0000001*(TMath::Abs(x)+TMath::Abs(fX[i]))) {
252  xidx=i;
253  LOG("BLI2DNonUnifGrid", pDEBUG) << "x value found at index "<<i;
254  changex=false; break;}
255  changex=true;
256  }
257  }
258  if (changex && xidx<0) xidx=fNFillX;
259  if (changex) LOG("BLI2DNonUnifGrid", pDEBUG) << "Adding x value at index "<<xidx;
260 
261  for(int i=0;i<fNFillY;i++)
262  {
263  if (!(TMath::Abs(y-fY[i])<=.5*.0000001*(TMath::Abs(y)+TMath::Abs(fY[i]))) && y<fY[i])
264  {
265  if (fNFillY==fNY)
266  {
267  LOG("BLI2DNonUnifGrid", pWARN) << "Warning: too many y values, cannot add y = "<<y;
268  return false;
269  }
270  else
271  {
272  yidx=i;
273  changey=true;
274  break;
275  }
276  }
277  else
278  {
279  if (TMath::Abs(y-fY[i])<=.5*.0000001*(TMath::Abs(y)+TMath::Abs(fY[i]))) {
280  yidx=i;
281  LOG("BLI2DNonUnifGrid", pDEBUG) << "y value found at index "<<i;
282  changey=false; break;}
283  changey=true;
284  }
285  }
286  if (changey && yidx<0) yidx=fNFillY;
287  if (changey) LOG("BLI2DNonUnifGrid", pDEBUG) << "Adding y value at index "<<yidx;
288 
289  // make new entries if needed
290  if (changex && xidx>=0)
291  {
292  for (int j=0;j<fNFillX-xidx;j++)
293  {
294  fX[fNFillX-j]=fX[fNFillX-j-1];
295  for (int k=0;k<fNFillY;k++)
296  {
297  fZ[ this->IdxZ(fNFillX-j,k) ]
298  = fZ[ this->IdxZ(fNFillX-j-1,k) ];
299  }
300  }
301  fX[xidx]=x;
302  fNFillX++;
303  fXmin = TMath::Min(x,fXmin);
304  fXmax = TMath::Max(x,fXmax);
305  }
306  if (changey && yidx>=0)
307  {
308  for (int j=0;j<fNFillY-yidx;j++)
309  {
310  fY[fNFillY-j]=fY[fNFillY-j-1];
311  for (int k=0;k<fNFillX;k++)
312  {
313  fZ[ this->IdxZ(k,fNFillY-j) ]
314  = fZ[ this->IdxZ(k,fNFillY-j-1) ];
315  }
316  }
317  fY[yidx]=y;
318  fNFillY++;
319  fYmin = TMath::Min(y,fYmin);
320  fYmax = TMath::Max(y,fYmax);
321  }
322 
323  int iz = this->IdxZ(xidx,yidx);
324 
325  fZ[iz] = z;
326 
327  fZmin = TMath::Min(z, fZmin);
328  fZmax = TMath::Max(z, fZmax);
329 
330  LOG("BLI2DNonUnifGrid", pDEBUG)
331  << "Added x = " << x << " (ix = " << xidx << ")"
332  << " y = " << y << " (iy = " << yidx << ") -> "
333  << " z = " << z << " (iz = " << iz << ")";
334 
335  return true;
336 }
337 //___________________________________________________________________________
338 double BLI2DNonUnifGrid::Evaluate(double x, double y) const
339 {
340  double evalx=TMath::Min(x,fXmax);
341  evalx=TMath::Max(evalx,fXmin);
342  double evaly=TMath::Min(y,fYmax);
343  evaly=TMath::Max(evaly,fYmin);
344 
345  int ix_lo = -2;
346  int iy_lo = -2;
347  for (int i=0;i<fNFillX;i++)
348  {
349  if (evalx<=fX[fNFillX-1-i]) ix_lo=fNFillX-2-i;
350  else break;
351  }
352  for (int i=0;i<fNFillY;i++)
353  {
354  if (evaly<=fY[fNFillY-1-i]) iy_lo=fNFillY-2-i;
355  else break;
356  }
357  int ix_hi = ix_lo + 1;
358  int iy_hi = iy_lo + 1;
359 
360  // in case x = xmin
361  if (ix_lo==-1) {ix_lo++; ix_hi++;}
362  if (iy_lo==-1) {iy_lo++; iy_hi++;}
363  // in case x = xmax
364  if (ix_lo==-2) {ix_lo=fNFillX-2; ix_hi=fNFillX-1;}
365  if (iy_lo==-2) {iy_lo=fNFillY-2; iy_hi=fNFillY-1;}
366  // if an error occurs
367  if (ix_lo<0 || iy_lo<0 ) return 0.;
368  if (ix_hi>fNFillX || iy_hi>fNFillY) return 0.;
369 
370  double x1 = fX[ix_lo];
371  double x2 = fX[ix_hi];
372  double y1 = fY[iy_lo];
373  double y2 = fY[iy_hi];
374 
375  double z11 = fZ[ this->IdxZ(ix_lo,iy_lo) ];
376  double z21 = fZ[ this->IdxZ(ix_hi,iy_lo) ];
377  double z12 = fZ[ this->IdxZ(ix_lo,iy_hi) ];
378  double z22 = fZ[ this->IdxZ(ix_hi,iy_hi) ];
379 
380  double z1 = z11 * (x2-evalx)/(x2-x1) + z21 * (evalx-x1)/(x2-x1);
381  double z2 = z12 * (x2-evalx)/(x2-x1) + z22 * (evalx-x1)/(x2-x1);
382  double z = z1 * (y2-evaly)/(y2-y1) + z2 * (evaly-y1)/(y2-y1);
383 
384  /*
385  LOG("BLI2DNonUnifGrid", pINFO) << "x = " << evalx << " -> nearby nodes: " << x1 << ", " << x2;
386  LOG("BLI2DNonUnifGrid", pINFO) << "y = " << evaly << " -> nearby nodes: " << y1 << ", " << y2;
387  LOG("BLI2DNonUnifGrid", pINFO) << "xmin = " << fXmin << ", xmax = " << fXmax;
388  LOG("BLI2DNonUnifGrid", pINFO) << "ymin = " << fYmin << ", ymax = " << fYmax;
389  LOG("BLI2DNonUnifGrid", pINFO) << "z11 := z(" << this->IdxZ(ix_lo,iy_lo) << ") = " << z11;
390  LOG("BLI2DNonUnifGrid", pINFO) << "z21 := z(" << this->IdxZ(ix_hi,iy_lo) << ") = " << z21;
391  LOG("BLI2DNonUnifGrid", pINFO) << "z12 := z(" << this->IdxZ(ix_lo,iy_hi) << ") = " << z12;
392  LOG("BLI2DNonUnifGrid", pINFO) << "z22 := z(" << this->IdxZ(ix_hi,iy_hi) << ") = " << z22;
393  LOG("BLI2DNonUnifGrid", pINFO) << "z1 = " << z1 << ", z2 = " << z2;
394  LOG("BLI2DNonUnifGrid", pINFO) << "interpolated z(x,y) = " << z;
395  */
396 
397  return z;
398 }
399 //___________________________________________________________________________
401  int nx, double xmin, double xmax, int ny, double ymin, double ymax)
402 {
403  fNX = 0;
404  fNY = 0;
405  fNZ = 0;
406  fNFillX= 0;
407  fNFillY= 0;
408  fXmin = 0.;
409  fXmax = 0.;
410  fYmin = 0.;
411  fYmax = 0.;
414  fX = 0;
415  fY = 0;
416  fZ = 0;
417 
418  if(nx>1 && ny>1) {
419  fNX = nx;
420  fNY = ny;
421  fNZ = nx * ny;
422 
423  fXmin = xmin;
424  fXmax = xmax;
425  fYmin = ymin;
426  fYmax = ymax;
429 
430  fX = new double[fNX];
431  fY = new double[fNY];
432  fZ = new double[fNZ];
433 
434  for(int i=0; i<fNX; i++) { fX[i] = 0.; }
435  for(int i=0; i<fNY; i++) { fY[i] = 0.; }
436  for(int i=0; i<fNZ; i++) { fZ[i] = 0.; }
437  }
438 }
439 //___________________________________________________________________________
double fXmin
Definition: BLI2D.h:63
double fYmax
Definition: BLI2D.h:66
double * fY
Definition: BLI2D.h:59
Bilinear interpolation of 2D functions on a regular grid.
Definition: BLI2D.h:76
bool AddPoint(double x, double y, double z)
Definition: BLI2D.cxx:218
double fZmax
Definition: BLI2D.h:68
std::map< std::string, double > xmax
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
double * fX
Definition: BLI2D.h:58
double Evaluate(double x, double y) const
Definition: BLI2D.cxx:338
double fDX
Definition: BLI2D.h:61
void Init(int nx=0, double xmin=0, double xmax=0, int ny=0, double ymin=0, double ymax=0)
Definition: BLI2D.cxx:143
Double_t ymax
Definition: plot.C:25
bool AddPoint(double x, double y, double z)
Definition: BLI2D.cxx:86
double fDY
Definition: BLI2D.h:62
virtual void Init(int nx, double xmin, double xmax, int ny, double ymin, double ymax)=0
void Init(int nx=0, double xmin=0, double xmax=0, int ny=0, double ymin=0, double ymax=0)
Definition: BLI2D.cxx:400
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
virtual bool AddPoint(double x, double y, double z)=0
const double j
Definition: BetheBloch.cxx:29
ClassImp(BLI2DGrid) BLI2DGrid
Definition: BLI2D.cxx:34
z
Definition: test.py:28
#define pWARN
Definition: Messenger.h:61
double Evaluate(double x, double y) const
Definition: BLI2D.cxx:105
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
int IdxZ(int ix, int iy) const
Definition: BLI2D.cxx:49
double fYmin
Definition: BLI2D.h:65
virtual ~BLI2DGrid()
Definition: BLI2D.cxx:42
Double_t ymin
Definition: plot.C:24
double fXmax
Definition: BLI2D.h:64
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
double fZmin
Definition: BLI2D.h:67
double * fZ
Definition: BLI2D.h:60
#define pDEBUG
Definition: Messenger.h:64