SpatialOps
FieldComparisons.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2014 The University of Utah
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining a copy
5  * of this software and associated documentation files (the "Software"), to
6  * deal in the Software without restriction, including without limitation the
7  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
8  * sell copies of the Software, and to permit persons to whom the Software is
9  * furnished to do so, subject to the following conditions:
10  *
11  * The above copyright notice and this permission notice shall be included in
12  * all copies or substantial portions of the Software.
13  *
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
19  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
20  * IN THE SOFTWARE.
21  */
22 #ifndef SpatialOps_FieldComparisons_h
23 #define SpatialOps_FieldComparisons_h
24 
25 #include <spatialops/Nebo.h>
26 #include <spatialops/structured/SpatialFieldStore.h>
27 #include <sstream>
28 #include <cmath>
29 
30 // Boost includes //
31 #include <boost/math/special_functions/next.hpp>
32 
33 #define FIELDCOMPARISONS_ABS_ERROR_CONST .000001
34 
42 namespace SpatialOps{
43 
44 template<typename FieldT>
46 
66 template<typename FieldT>
67 bool field_not_equal(const FieldT& f1, const FieldT& f2, double error=0.0) {
68  double error_abs = error ? nebo_norm(f1)*error*FIELDCOMPARISONS_ABS_ERROR_CONST : 0;
69  return !field_equal(f1, f2, error, error_abs);
70 }
71 
87 template<typename FieldT>
88 bool field_not_equal(const FieldT& f1, const FieldT& f2, double error, const double error_abs) {
89  return !field_equal(f1, f2, error, error_abs);
90 }
91 //------------------------------------------------------------------
92 
112 template<typename FieldT>
113 bool field_equal(const FieldT& f1, const FieldT& f2, double error=0.0)
114 {
115  double error_abs = error ? nebo_norm(f1)*error*FIELDCOMPARISONS_ABS_ERROR_CONST : 0;
116  return field_equal(f1, f2, error, error_abs);
117 }
118 
135 template<typename FieldT>
136 bool field_equal(const FieldT& f1, const FieldT& f2, double error, const double errorAbs)
137 {
138  const MemoryWindow& w1 = f1.window_with_ghost();
139  const MemoryWindow& w2 = f2.window_with_ghost();
140 
141  if( w1 != w2 ){
142  throw( std::runtime_error( "Attempted comparison between fields of unequal size." ) );
143  }
144 
145  error = std::abs(error);
146  const bool exactComparison = error == 0.0;
147 
148  SpatFldPtr<FieldT> tmp1, tmp2;
149  typename FieldT::const_iterator if1, iend, if2;
150 
153 
154  //do comparison
155  for( ; if2 != iend; ++if1, ++if2 ){
156  const typename FieldT::value_type combined = *if1 - *if2;
157  if( std::isnan( combined ) || std::isinf( combined ) ) return false;
158  if( exactComparison ){
159  if( *if1 != *if2 ) return false;
160  }
161  else {
162  const double denom = std::abs(*if1) + errorAbs;
163  if( std::abs(combined)/denom > error ) return false;
164  }
165  }
166  return true;
167 }
168 //------------------------------------------------------------------
169 
183 template<typename FieldT>
184 bool field_not_equal_abs(const FieldT& f1, const FieldT& f2, double error=0.0) {
185  return !field_equal_abs(f1, f2, error);
186 }
187 //------------------------------------------------------------------
188 
200 template<typename FieldT>
201 bool field_equal_abs( const FieldT& f1, const FieldT& f2, double error=0.0 )
202 {
203  const MemoryWindow w1 = f1.window_with_ghost();
204  const MemoryWindow w2 = f2.window_with_ghost();
205 
206  if( w1 != w2 ){
207  throw( std::runtime_error( "Attempted comparison between fields of unequal size." ) );
208  }
209 
210  error = std::abs(error);
211  const bool exactComparison = (error == 0.0);
212 
213  SpatFldPtr<FieldT> tmp1, tmp2;
214  typename FieldT::const_iterator if1, iend, if2;
215 
216  FieldComparisonHelper<FieldT>::init_iterator( f1, tmp1, if1, iend );
217  FieldComparisonHelper<FieldT>::init_iterator( f2, tmp2, if2, iend );
218 
219  // do comparison at each point
220  for( ; if2 != iend; ++if1, ++if2 ){
221  if( std::isnan( *if1 + *if2 ) ) return false;
222  if( exactComparison ){
223  if(*if1 != *if2) return false;
224  }
225  else{
226  if( std::abs(*if1 - *if2) > error ) return false;
227  }
228  }
229  return true;
230 }
231 //------------------------------------------------------------------
232 
246 template<typename FieldT>
247 bool field_not_equal_ulp(const FieldT& f1, const FieldT& f2, const unsigned int ulps) {
248  return !field_equal_ulp(f1, f2, ulps);
249 }
250 //------------------------------------------------------------------
251 
269 template<typename FieldT>
270 bool field_equal_ulp(const FieldT& f1, const FieldT& f2, const unsigned int ulps)
271 {
272  const MemoryWindow& w1 = f1.window_with_ghost();
273  const MemoryWindow& w2 = f2.window_with_ghost();
274 
275  if(w1 != w2) {
276  throw( std::runtime_error( "Attempted comparison between fields of unequal size." ) );
277  }
278 
279  const bool exactComparison = ulps == 0;
280  SpatFldPtr<FieldT> tmp1, tmp2;
281  typename FieldT::const_iterator if1, iend, if2;
282 
283  FieldComparisonHelper<FieldT>::init_iterator( f1, tmp1, if1, iend );
284  FieldComparisonHelper<FieldT>::init_iterator( f2, tmp2, if2, iend );
285 
286  //do comparison
287  for( ; if2 != iend; ++if1, ++if2 ){
288  if( std::isnan( *if1 + *if2 ) ) return false;
289  if( exactComparison ){
290  if( boost::math::float_distance(*if1, *if2) != 0) return false;
291  }
292  else {
293  if (std::abs(boost::math::float_distance(*if1, *if2)) > ulps) return false;
294  }
295  }
296  return true;
297 }
298 
304 
323 template<typename FieldT>
324 bool field_not_equal(const double d, const FieldT& f1, double error=0.0) {
325  double error_abs = error ? nebo_norm(f1)*error*FIELDCOMPARISONS_ABS_ERROR_CONST : 0;
326  return !field_equal(d, f1, error, error_abs);
327 }
328 
344 template<typename FieldT>
345 bool field_not_equal(const double d, const FieldT& f1, double error, const double error_abs) {
346  return !field_equal(d, f1, error, error_abs);
347 }
348 //------------------------------------------------------------------
349 
367 template<typename FieldT>
368 bool field_equal(const double d, const FieldT& f1, double error=0.0)
369 {
370  double error_abs = error ? nebo_norm(f1)*error*FIELDCOMPARISONS_ABS_ERROR_CONST : 0;
371  return field_equal(d, f1, error, error_abs);
372 }
373 
391 template<typename FieldT>
392 bool field_equal(const double d, const FieldT& f1, double error, const double errorAbs)
393 {
394  error = std::abs(error);
395  const bool exactComparison = error == 0.0;
396 
397  SpatFldPtr<FieldT> tmp;
398  typename FieldT::const_iterator if1, iend;
399 
401 
402  //do comparison
403  const double denom = std::abs(d) + errorAbs;
404  for( ; if1 != iend; ++if1 ){
405  const typename FieldT::value_type combined = d - *if1;
406  if( std::isnan( combined ) || std::isinf( combined ) ) return false;
407  if( exactComparison ){
408  if( *if1 != d ) return false;
409  }
410  else{
411  if( std::abs(combined)/denom > error ) return false;
412  }
413  }
414  return true;
415 }
416 //------------------------------------------------------------------
417 
432 template<typename FieldT>
433 bool field_not_equal_abs(const double d, const FieldT& f1, const double error=0.0) {
434  return !field_equal_abs(d, f1, error);
435 }
436 //------------------------------------------------------------------
437 
451 template<typename FieldT>
452 bool field_equal_abs(const double d, const FieldT& f1, double error=0.0)
453 {
454  error = std::abs(error);
455  const bool exactComparison = error == 0.0;
456 
457  SpatFldPtr<FieldT> tmp;
458  typename FieldT::const_iterator if1, iend;
459 
461 
462  //do comparison
463  for( ; if1 != iend; ++if1 ){
464  if( std::isnan( d + *if1 ) ) return false;
465  if( exactComparison ){
466  if( *if1 != d ) return false;
467  }
468  else{
469  if( std::abs(d - *if1) > error ) return false;
470  }
471  }
472  return true;
473 }
474 //------------------------------------------------------------------
475 
491 template<typename FieldT>
492 bool field_not_equal_ulp(const double d, const FieldT& f1, const unsigned int ulps) {
493  return !field_equal_ulp(d, f1, ulps);
494 }
495 //------------------------------------------------------------------
496 
516 template<typename FieldT>
517 bool field_equal_ulp(const double d, const FieldT& f1, const unsigned int ulps)
518 {
519  bool exactComparison = ulps == 0;
520 
521  SpatFldPtr<FieldT> tmp;
522  typename FieldT::const_iterator if1, iend;
523 
525 
526  //do comparison
527  for( ; if1 != iend; ++if1 ){
528  if( std::isnan( *if1 ) ) return false;
529  if( exactComparison ){
530  if( boost::math::float_distance(d, *if1) != 0 ) return false;
531  }
532  else{
533  if( std::abs(boost::math::float_distance(d, *if1)) > ulps ) return false;
534  }
535  }
536  return true;
537 }
538 
550 template<typename FieldT>
552 {
570  private:
571  inline static void
572  init_iterator( const FieldT& field,
573  SpatFldPtr<FieldT>& fcopy,
574  typename FieldT::const_iterator& ibegin,
575  typename FieldT::const_iterator& iend )
576  {
577  //we need to transfer the field to local ram to iterate
578  if( IS_CPU_INDEX( field.active_device_index() ) ){
579  ibegin = field.begin();
580  iend = field.end();
581  }
582 # ifdef ENABLE_CUDA
583  else if( IS_GPU_INDEX(field.active_device_index()) ){
584  fcopy = SpatialFieldStore::get<FieldT>(field);
585  fcopy->set_stream(field.get_stream());
586  *fcopy <<= field;
587  fcopy->add_device(CPU_INDEX);
588  //
589  // jcs hack to get things working. There are two issues here:
590  //
591  // 1. in order to get begin() to use the const version, we need to have
592  // a const field, which "fcopy" is not.
593  //
594  // 2. For the non-const field to work with a non-const version of begin()
595  // we need to first set the active field location to the CPU. This
596  // Seems to cause problems with the memory pool due to a bug that
597  // has not yet been resolved.
598  //
599  const FieldT* const fcopyA = &(*fcopy);
600  ibegin = fcopyA->begin();
601  iend = fcopyA->end();
602  return;
603  }
604 # endif
605  else{
606  throw std::runtime_error("Unrecognized field location (FieldComparisons.h FieldCopmarisonHelper::init_iterator()");
607  }
608  }
609 
610  friend bool field_equal<FieldT>(const FieldT& f1, const FieldT& f2, double error, const double error_abs);
611  friend bool field_equal_abs<FieldT>(const FieldT& f1, const FieldT& f2, double error);
612  friend bool field_equal_ulp<FieldT>(const FieldT& f1, const FieldT& f2, const unsigned int ulps);
613 
614  friend bool field_equal<FieldT>(const double d, const FieldT& f1, double error, const double error_abs);
615  friend bool field_equal_abs<FieldT>(const double d, const FieldT& f1, double error);
616  friend bool field_equal_ulp<FieldT>(const double d, const FieldT& f1, const unsigned int ulps);
617 };
618 
619 } // namespace SpatialOps
620 
621 #endif //SpatialOps_FieldComparisons_h
bool field_equal_ulp(const FieldT &f1, const FieldT &f2, const unsigned int ulps)
Returns if f1 is element-wise equal to f2 within a certain number of ulps.
bool field_not_equal_abs(const FieldT &f1, const FieldT &f2, double error=0.0)
Returns if f1 is element-wise not equal to f2 within a certain absolute tolerance.
Provides tools to index into a sub-block of memory.
Definition: MemoryWindow.h:63
bool field_equal_abs(const FieldT &f1, const FieldT &f2, double error=0.0)
Returns if f1 is element-wise equal to f2 within a certain absolute tolerance.
bool field_not_equal_ulp(const FieldT &f1, const FieldT &f2, const unsigned int ulps)
Returns if f1 is element-wise not equal to f2 within a certain number of ulps.
static helper class for field comparison functions
bool field_not_equal(const FieldT &f1, const FieldT &f2, double error=0.0)
Returns if f1 is element-wise not equal to f2 within a certain relative tolerance.This function simply calls field_equal and negates it. return !field_equal(f1, f2, error, error_abs); error_abs is defined as default to be the L2 norm of f1 multiplied by error and FIELDCOMPARISONS_ABS_ERROR_CONST.
bool field_equal(const double d, const FieldT &f1, double error, const double errorAbs)
Returns if f1 is element-wise equal to the scalar value d within a certain relative tolerance...
Wrapper for pointers to SpatialField objects. Provides reference counting.