BALL 1.5.0
regularData1D.h
Go to the documentation of this file.
1// -*- Mode: C++; tab-width: 2; -*-
2// vi: set ts=2:
3//
4
5#ifndef BALL_DATATYPE_REGULARDATA1D_H
6#define BALL_DATATYPE_REGULARDATA1D_H
7
8#ifndef BALL_COMMON_H
9# include <BALL/common.h>
10#endif
11
12#ifndef BALL_SYSTEM_FILE_H
13# include <BALL/SYSTEM/file.h>
14#endif
15
16#ifndef BALL_SYSTEM_BINARYFILEADAPTOR_H
18#endif
19
20#include <vector>
21#include <iostream>
22#include <fstream>
23#include <iterator>
24#include <algorithm>
25#include <cmath>
26
27namespace BALL
28{
40 template <typename ValueType>
42 {
43 public:
44
46
47
50
51
54 typedef std::vector<ValueType> VectorType;
56 typedef double CoordinateType;
58 typedef typename std::vector<ValueType>::iterator Iterator;
60 typedef typename std::vector<ValueType>::const_iterator ConstIterator;
62
63 // STL compatibility types
64 //
65 typedef ValueType value_type;
66 typedef typename std::vector<ValueType>::iterator iterator;
67 typedef typename std::vector<ValueType>::const_iterator const_iterator;
68 typedef typename std::vector<ValueType>::reference reference;
69 typedef typename std::vector<ValueType>::const_reference const_reference;
70 typedef typename std::vector<ValueType>::pointer pointer;
71 typedef typename std::vector<ValueType>::difference_type difference_type;
72 typedef typename std::vector<ValueType>::size_type size_type;
73
77
80
85
89 TRegularData1D(const CoordinateType& origin, const CoordinateType& dimension, const CoordinateType& spacing);
90
95
99 TRegularData1D(const VectorType& data, const CoordinateType& origin = 0.0, const CoordinateType& dimension = 1.0);
100
102 virtual ~TRegularData1D();
103
105 virtual void clear();
107
108
112
117 TRegularData1D& operator = (const TRegularData1D<ValueType>& data);
118
123 TRegularData1D& operator = (const VectorType& data);
124
126
131 bool operator == (const TRegularData1D& data) const;
132
134 BALL_INLINE bool operator != (const TRegularData1D& data) const { return !this->operator == (data); }
135
137 BALL_INLINE bool empty() const { return data_.empty(); }
138
140 bool isInside(const CoordinateType& x) const;
142
147 BALL_INLINE ConstIterator begin() const { return data_.begin(); }
149 BALL_INLINE ConstIterator end() const { return data_.end(); }
151 BALL_INLINE Iterator begin() { return data_.begin(); }
153 BALL_INLINE Iterator end() { return data_.end(); }
155
159
160 // STL compatibility
161 BALL_INLINE size_type size() const { return data_.size(); }
162 BALL_INLINE size_type max_size() const { return data_.max_size(); }
163 BALL_INLINE void swap(TRegularData1D<ValueType>& data) { std::swap(*this, data); }
164
169 const ValueType& getData(const IndexType& index) const;
170
175 ValueType& getData(const IndexType& index);
176
181 const ValueType& operator [] (const IndexType& index) const { return data_[index]; }
182
187 ValueType& operator [] (const IndexType& index) { return data_[index]; }
188
197 ValueType operator () (const CoordinateType& x) const;
198
205 ValueType getInterpolatedValue(const CoordinateType& x) const;
206
213 void getEnclosingIndices(const CoordinateType& x, Position& lower, Position& upper) const;
214
219 void getEnclosingValues(const CoordinateType& x, ValueType& lower, ValueType& upper) const;
220
226
234
242
249 const ValueType& getClosestValue(const CoordinateType& x) const;
250
257 ValueType& getClosestValue(const CoordinateType& x);
258
260 BALL_INLINE IndexType getSize() const { return (IndexType)data_.size(); }
261
266 BALL_INLINE const CoordinateType& getOrigin() const { return origin_; }
267
273
276 BALL_INLINE void setOrigin(const CoordinateType& origin) { origin_ = origin; }
277
284
290 BALL_INLINE void setDimension(const CoordinateType& dimension) { dimension_ = dimension; }
291
304 void resize(const IndexType& size);
305
315 void rescale(const IndexType& new_size);
316
320 ValueType calculateMean() const;
321
325 ValueType calculateSD() const;
326
330 void binaryWrite(const String& filename) const;
331
335 void binaryRead(const String& filename);
337
338
339 protected:
342
345
348
351
353 typedef struct { ValueType bt[1024]; } BlockValueType;
354 };
355
359
360 template <typename ValueType>
362 : origin_(0.0),
363 dimension_(0.0),
364 spacing_(1.0),
365 data_()
366 {
367 }
368
369 template <typename ValueType>
371 {
372 }
373
374 template <typename ValueType>
376 : origin_(data.origin_),
377 dimension_(data.dimension_),
378 spacing_(data.spacing_),
379 data_()
380 {
381 // Try to catch allocation errors and rethrow them as OutOfMemory
382 try
383 {
384 data_ = data.data_;
385 }
386 catch (std::bad_alloc&)
387 {
388 throw Exception::OutOfMemory(__FILE__, __LINE__, data.size() * sizeof(ValueType));
389 }
390 }
391
392 template <typename ValueType>
394 (const typename TRegularData1D<ValueType>::CoordinateType& origin,
395 const typename TRegularData1D<ValueType>::CoordinateType& dimension,
396 const typename TRegularData1D<ValueType>::CoordinateType& spacing)
397 : origin_(origin),
398 dimension_(dimension),
399 spacing_(spacing),
400 data_()
401 {
402 // Determine the size of the vector
404
405 // Try to catch allocation errors and rethrow them as OutOfMemory
406 try
407 {
408 data_.resize(size);
409 }
410 catch (std::bad_alloc&)
411 {
412 throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
413 }
414 }
415
416 template <typename ValueType>
418 (const typename TRegularData1D<ValueType>::VectorType& data,
419 const typename TRegularData1D<ValueType>::CoordinateType& origin,
420 const typename TRegularData1D<ValueType>::CoordinateType& dimension)
421 : origin_(origin),
422 dimension_(dimension),
423 spacing_(dimension / ((double)data.size()-1)),
424 data_()
425 {
426 // Try to catch allocation errors and rethrow them as OutOfMemory
427 try
428 {
429 data_ = data;
430 }
431 catch (std::bad_alloc&)
432 {
433 throw Exception::OutOfMemory(__FILE__, __LINE__, data.size() * sizeof(ValueType));
434 }
435 }
436
437
438 template <class ValueType>
440 (const typename TRegularData1D<ValueType>::IndexType& size)
441 : origin_(0.0),
442 dimension_(1.0),
443 data_()
444 {
445 // Compute the grid spacing
446 spacing_ = dimension_ / (double)(size - 1);
447
448 try
449 {
450 data_.resize(size);
451 }
452 catch (std::bad_alloc&)
453 {
454 data_.resize(0);
455 throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
456 }
457 }
458
459 template <typename ValueType>
461 {
462 // iterate over the data and reset all values to their default
463 // boundaries and vector size remain unchanged
464 static ValueType default_value = ValueType();
465 std::fill(data_.begin(), data_.end(), default_value);
466 }
467
468 template <typename ValueType>
470 {
471 // copy all members...
472 origin_ = rhs.origin_;
473 dimension_ = rhs.dimension_;
474 spacing_ = rhs.spacing_;
475 try
476 {
477 data_ = rhs.data_;
478 }
479 catch (std::bad_alloc&)
480 {
481 data_.resize(0);
482 throw Exception::OutOfMemory(__FILE__, __LINE__, rhs.size() * sizeof(ValueType));
483 }
484
485 return *this;
486 }
487
488 template <typename ValueType>
490 {
491 // Copy the data. The boundaries remain unchanged.
492 try
493 {
494 data_ = rhs;
495 }
496 catch (std::bad_alloc&)
497 {
498 data_.resize(0);
499 throw Exception::OutOfMemory(__FILE__, __LINE__, rhs.size() * sizeof(ValueType));
500 }
501
502 return *this;
503 }
504
505 template <typename ValueType>
507 {
508 return (origin_ == data.origin_
509 && dimension_ == data.dimension_
510 && data_ == data.data_);
511 }
512
513 template <class ValueType>
516 {
517 return ((r >= origin_) && (r <= (origin_ + dimension_)));
518 }
519
520 template <typename ValueType>
522 const ValueType& TRegularData1D<ValueType>::getData(const IndexType& index) const
523 {
524 if (index >= data_.size())
525 {
526 throw Exception::OutOfGrid(__FILE__, __LINE__);
527 }
528 return data_[index];
529 }
530
531 template <typename ValueType>
534 {
535 if (index >= data_.size())
536 {
537 throw Exception::OutOfGrid(__FILE__, __LINE__);
538 }
539 return data_[index];
540 }
541
542 template <typename ValueType>
545 Position& lower, Position& upper) const
546 {
547 if (!isInside(x) || (data_.size() < 2))
548 {
549 throw Exception::OutOfGrid(__FILE__, __LINE__);
550 }
551 lower = (Position)std::floor((x - origin_) / spacing_);
552 if (lower == data_.size() - 1)
553 {
554 // If we are on the right most data point, we cannot interpolate to the right!
555 lower = data_.size() - 2;
556 }
557 upper = lower + 1;
558 }
559
560 template <typename ValueType>
563 ValueType& lower, ValueType& upper) const
564 {
565 Position lower_index;
566 Position upper_index;
567 getEnclosingIndices(x, lower_index, upper_index);
568 lower = data_[lower_index];
569 upper = data_[upper_index];
570 }
571
572 template <typename ValueType>
575 {
576 if (!isInside(x))
577 {
578 throw Exception::OutOfGrid(__FILE__, __LINE__);
579 }
580 return operator () (x);
581 }
582
583 template <typename ValueType>
586 (const typename TRegularData1D<ValueType>::IndexType& index) const
587 {
588 if ((index >= data_.size()) || (data_.size() == 0))
589 {
590 throw Exception::OutOfGrid(__FILE__, __LINE__);
591 }
592
593 return (CoordinateType)(origin_ + (double)index / ((double)data_.size()-1) * dimension_);
594 }
595
596 template <typename ValueType>
599 {
600 if ((x < origin_) || (x > (origin_ + dimension_)))
601 {
602 throw Exception::OutOfGrid(__FILE__, __LINE__);
603 }
604
605 return (IndexType)(size_type)std::floor((x - origin_) / spacing_ + 0.5);
606 }
607
608 template <typename ValueType>
611 {
612 if ((x < origin_) || (x > (origin_ + dimension_)))
613 {
614 throw Exception::OutOfGrid(__FILE__, __LINE__);
615 }
616
617 return (IndexType)(size_type)std::floor((x - origin_) / spacing_);
618 }
619
620 template <typename ValueType>
623 {
624 if ((x < origin_) || (x > (origin_ + dimension_)))
625 {
626 throw Exception::OutOfGrid(__FILE__, __LINE__);
627 }
628
629 // Round to the closest data point.
630 size_type index = (size_type)std::floor((x - origin_) / spacing_ + 0.5);
631 return data_[index];
632 }
633
634 template <typename ValueType>
637 {
638 if ((x < origin_) || (x > (origin_ + dimension_)))
639 {
640 throw Exception::OutOfGrid(__FILE__, __LINE__);
641 }
642
643 // Round to the closest data point.
644 size_type index = (size_type)std::floor((x - origin_) / spacing_ + 0.5);
645 return data_[index];
646 }
647
648 template <typename ValueType>
651 {
652 IndexType data_points = this->getSize();
653 ValueType mean = 0;
654 for (IndexType i = 0; i < data_points; i++)
655 {
656 mean += data_[i];
657 }
658 mean /= data_points;
659 return mean;
660 }
661
662 template <typename ValueType>
665 {
666 IndexType data_points = this->getSize();
667 ValueType stddev = 0;
668 ValueType mean = this->calculateMean();
669 for (IndexType i = 0; i < data_points; i++)
670 {
671 stddev += (pow(data_[i]-mean,2));
672 }
673 stddev /= (data_points-1);
674 stddev = sqrt(stddev);
675 return stddev;
676 }
677
678 template <typename ValueType>
681 {
682 size_type left_index = (size_type)std::floor((x - origin_) / spacing_);
683 if (left_index == data_.size() - 1)
684 {
685 // If we are on the right most data point, we cannot interpolate to the right!
686 return data_[data_.size() - 1];
687 }
688
689 // Interpolate between the point to the left and the point to the right.
690 double d = 1.0 - (((x - origin_) - (double)left_index * spacing_) / spacing_);
691 return data_[left_index] * d + (1.0 - d) * data_[left_index + 1];
692 }
693
694 template <typename ValueType>
696 (const typename TRegularData1D<ValueType>::IndexType& new_size)
697 {
698 // Rescale dimension to the new size.
699 if (data_.size() > 0)
700 {
701 dimension_ *= (double)new_size / (double)data_.size();
702 }
703
704 // Try to resize the vactor and rethrow any bad_allocs.
705 try
706 {
707 data_.resize(new_size);
708 }
709 catch (std::bad_alloc&)
710 {
711 // The resulting vector is empty and thus well-defined.
712 data_.resize(0);
713 throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * sizeof(ValueType));
714 }
715 }
716
717 template <typename ValueType>
719 (const typename TRegularData1D<ValueType>::IndexType& new_size)
720 {
721 // if the new and the old size coincide: done.
722 if (new_size == (IndexType)data_.size())
723 {
724 return;
725 }
726
727 // Catch any bad_allocs throw by vector::resize
728 try
729 {
730 // if the data set is empty...
731 if (data_.size() == 0)
732 {
733 // ...there's nothing to do: a resize was requested
734 data_.resize(new_size);
735 return;
736 }
737
738 // if the data set contains only a single value,
739 // we fill everything with this value
740 if ((data_.size() == 1) && (new_size > 1))
741 {
742 ValueType old_value = data_[0];
743 data_.resize(new_size);
744 for (IndexType i = 1; i < new_size; i++)
745 {
746 data_[i] = old_value;
747 }
748
749 return;
750 }
751
752 // that's the default case: use linear interpolation
753 // to determine the values at the new positions
754 VectorType new_data(new_size);
755 CoordinateType factor1 = (CoordinateType)data_.size() / (CoordinateType)new_size;
756 CoordinateType factor2 = (CoordinateType)(data_.size() - 1) / (new_size - 1);
757
758 for (Size i = 0; i < new_size; i++)
759 {
760 // determine the interval of the old data set we are currently in
761 // ([old_idx, old_idx + 1])
762 IndexType old_idx = (IndexType)((CoordinateType)i * factor1);
763
764 // consider numerical inaccuracies...
765 if (old_idx >= (data_.size() - 1))
766 {
767 old_idx = data_.size() - 2;
768 }
769 CoordinateType factor3 = (CoordinateType)i * factor2 - (CoordinateType)old_idx;
770 new_data[i] = data_[old_idx] * (1 - factor3) + factor3 * data_[old_idx + 1];
771 }
772
773 // assign the new data
774 data_ = new_data;
775 }
776 catch (std::bad_alloc&)
777 {
778 // Make sure we are in a well-defined state.
779 data_.resize(0);
780 throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * sizeof(ValueType));
781 }
782 }
783
787 template <typename ValueType>
788 std::ostream& operator << (std::ostream& os, const TRegularData1D<ValueType>& data)
789 {
790 // Write the grid origin, dimension, and number of grid points
791 os << data.getOrigin() << std::endl
792 << data.getOrigin() + data.getDimension() << std::endl
793 << data.getSize() - 1 << std::endl;
794
795 // Write the array contents.
796 std::copy(data.begin(), data.end(), std::ostream_iterator<ValueType>(os, "\n"));
797 return os;
798 }
799
801 template <typename ValueType>
802 std::istream& operator >> (std::istream& is, TRegularData1D<ValueType>& grid)
803 {
807
808 is >> origin;
809 is >> dimension;
810 is >> size;
811
812 dimension -= origin;
813 size++;
814
815 grid.resize(size);
816 grid.setOrigin(origin);
817 grid.setDimension(dimension);
818
819 std::copy(std::istream_iterator<ValueType>(is),
820 std::istream_iterator<ValueType>(),
821 grid.begin());
822 // std::copy_n(std::istream_iterator<ValueType>(is), grid.size(), grid.begin());
823
824 return is;
825 }
826
827 template <typename ValueType>
829 {
830 File outfile(filename, std::ios::out|std::ios::binary);
831 if (!outfile.isValid())
832 {
833 throw Exception::FileNotFound(__FILE__, __LINE__, filename);
834 }
835
837 BinaryFileAdaptor<ValueType> adapt_single;
838
839 // write all information we need to recreate the grid
840 BinaryFileAdaptor<CoordinateType> adapt_coordinate;
841 BinaryFileAdaptor<Size> adapt_size;
842
843 adapt_size.setData(data_.size());
844 outfile << adapt_size;
845
846 adapt_coordinate.setData(origin_);
847 outfile << adapt_coordinate;
848
849 adapt_coordinate.setData(dimension_);
850 outfile << adapt_coordinate;
851
852 adapt_coordinate.setData(spacing_);
853 outfile << adapt_coordinate;
854
855 // we slide a window of size 1024 over our data
856 Index window_pos = 0;
857 while (((int)data_.size() - (1024 + window_pos)) >= 0 )
858 {
859 adapt_block.setData(*(BlockValueType*)&(data_[window_pos]));
860 outfile << adapt_block;
861 window_pos += 1024;
862 }
863
864 // now we have to write the remaining data one by one
865 for (Size i = window_pos; i < data_.size(); i++)
866 {
867 adapt_single.setData(data_[i]);
868 outfile << adapt_single;
869 }
870
871 // that's it. I hope...
872 outfile.close();
873 }
874
875 template <typename ValueType>
877 {
878 File infile(filename, std::ios::in|std::ios::binary);
879 if (!infile.isValid()) throw Exception::FileNotFound(__FILE__, __LINE__, filename);
880
882 BinaryFileAdaptor<ValueType> adapt_single;
883
884 // read all information we need to recreate the grid
885 BinaryFileAdaptor<CoordinateType> adapt_coordinate;
886 BinaryFileAdaptor<Size> adapt_size;
887
888 infile >> adapt_size;
889 Size new_size = adapt_size.getData();
890
891 infile >> adapt_coordinate;
892 origin_ = adapt_coordinate.getData();
893
894 infile >> adapt_coordinate;
895 dimension_ = adapt_coordinate.getData();
896
897 infile >> adapt_coordinate;
898 spacing_ = adapt_coordinate.getData();
899
900 data_.resize(new_size);
901
902 // we slide a window of size 1024 over our data
903 Index window_pos = 0;
904
905 while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
906 {
907 infile >> adapt_block;
908 *(BlockValueType*)(&(data_[window_pos])) = adapt_block.getData();
909 /*
910 for (Size i=0; i<1024; i++)
911 {
912 data_[i+window_pos] = adapt_block.getData().bt[i];
913 }
914 */
915 window_pos+=1024;
916 }
917
918 // now we have to read the remaining data one by one
919 for (Size i=window_pos; i<data_.size(); i++)
920 {
921 infile >> adapt_single;
922 data_[i] = adapt_single.getData();
923 }
924
925 // that's it. I hope...
926 infile.close();
927 }
928} // namespace BALL
929
930#endif // BALL_DATATYPE_REGULARDATA1D_H
#define BALL_CREATE(name)
Definition: create.h:62
#define BALL_INLINE
Definition: debug.h:15
BALL_EXPORT std::ostream & operator<<(std::ostream &os, const Exception::GeneralException &e)
STL namespace.
Definition: constants.h:13
TRegularData1D< float > RegularData1D
BALL_SIZE_TYPE Position
std::istream & operator>>(std::istream &is, TRegularData1D< ValueType > &grid)
Input operator.
std::ifstream infile
Definition: classTest.h:27
long floor(const T &t)
Definition: MATHS/common.h:284
CoordinateType dimension_
The dimension (length)
const ValueType & getClosestValue(const CoordinateType &x) const
BALL_INLINE size_type max_size() const
ValueType getInterpolatedValue(const CoordinateType &x) const
BALL_INLINE size_type size() const
void rescale(const IndexType &new_size)
IndexType getLowerIndex(const CoordinateType &x) const
CoordinateType origin_
The origin of the data set.
BALL_INLINE const CoordinateType & getOrigin() const
std::vector< ValueType > VectorType
The type containing an STL vector of the corresponding ValueType.
Definition: regularData1D.h:54
BALL_INLINE const CoordinateType & getDimension() const
virtual void clear()
Clear the contents.
std::vector< ValueType >::pointer pointer
Definition: regularData1D.h:70
IndexType getClosestIndex(const CoordinateType &x) const
void getEnclosingValues(const CoordinateType &x, ValueType &lower, ValueType &upper) const
BALL_INLINE void setDimension(const CoordinateType &dimension)
std::vector< ValueType >::iterator iterator
Definition: regularData1D.h:66
void resize(const IndexType &size)
bool operator==(const TRegularData1D &data) const
Equality operator.
const ValueType & getData(const IndexType &index) const
BALL_INLINE void setOrigin(const CoordinateType &origin)
BALL_INLINE IndexType getSize() const
Return the number of points in the data set.
std::vector< ValueType >::iterator Iterator
A mutable iterator.
Definition: regularData1D.h:58
double CoordinateType
The coordinate type.
Definition: regularData1D.h:56
ValueType calculateSD() const
std::vector< ValueType >::difference_type difference_type
Definition: regularData1D.h:71
std::vector< ValueType >::size_type size_type
Definition: regularData1D.h:72
std::vector< ValueType >::const_iterator const_iterator
Definition: regularData1D.h:67
BALL_INLINE void swap(TRegularData1D< ValueType > &data)
ValueType calculateMean() const
virtual ~TRegularData1D()
Destructor.
TRegularData1D()
Default constructor.
Position IndexType
The IndexType.
Definition: regularData1D.h:52
BALL_INLINE Iterator end()
BALL_INLINE const CoordinateType & getSpacing() const
TRegularData1D & operator=(const TRegularData1D< ValueType > &data)
VectorType data_
The data.
CoordinateType getCoordinates(const IndexType &index) const
std::vector< ValueType >::reference reference
Definition: regularData1D.h:68
ValueType operator()(const CoordinateType &x) const
BALL_INLINE ConstIterator begin() const
ValueType & getClosestValue(const CoordinateType &x)
BALL_INLINE bool empty() const
Empty predicate.
BALL_INLINE Iterator begin()
CoordinateType spacing_
The spacing.
std::vector< ValueType >::const_iterator ConstIterator
A constant iterator.
Definition: regularData1D.h:60
void binaryWrite(const String &filename) const
void binaryRead(const String &filename)
const ValueType & operator[](const IndexType &index) const
ValueType & getData(const IndexType &index)
std::vector< ValueType >::const_reference const_reference
Definition: regularData1D.h:69
bool isInside(const CoordinateType &x) const
Test whether a point is inside the grid.
void getEnclosingIndices(const CoordinateType &x, Position &lower, Position &upper) const
BALL_INLINE ConstIterator end() const
The block data type for reading and writing binary data.
const T & getData() const
void setData(const T &data)
bool isValid() const
void close()