pele
Python energy landscape explorer
/home/js850/projects/pele/source/pele/array.h
00001 #ifndef PYGMIN_ARRAY_H
00002 #define PYGMIN_ARRAY_H
00003 
00004 #include <algorithm>
00005 #include <cassert>
00006 #include <cmath>
00007 #include <iostream>
00008 #include <memory>
00009 #include <numeric>
00010 #include <stdexcept>
00011 #include <vector>
00012 
00013 namespace pele {
00014 
00019 template<typename dtype>
00020 class _ArrayMemory {
00021     std::vector<dtype> _vector;
00022     dtype *_data; 
00025     size_t _size;  
00027 public:
00028     _ArrayMemory()
00029         : _vector(),
00030           _data(_vector.data()),
00031           _size(_vector.size())
00032     {}
00033 
00034     _ArrayMemory(size_t size)
00035         : _vector(size),
00036           _data(_vector.data()),
00037           _size(_vector.size())
00038     {}
00039 
00040     _ArrayMemory(size_t size, dtype const & val)
00041         : _vector(size, val),
00042           _data(_vector.data()),
00043           _size(_vector.size())
00044     {}
00045 
00049     _ArrayMemory(dtype * data, size_t size)
00050         : _vector(),
00051           _data(data),
00052           _size(size)
00053     {}
00054 
00058     inline size_t size() const { return _size; }
00059 
00063     inline dtype *data() { return _data; }
00064     inline dtype const *data() const { return _data; }
00065 };
00066 
00067 template<typename dtype>
00068 class Array
00069 {
00070     std::shared_ptr<_ArrayMemory<dtype> > _memory;
00071     dtype * _data; 
00074     size_t _size;   
00075 public:
00081     Array()
00082         : _memory(new _ArrayMemory<dtype>()),
00083           _data(_memory->data()),
00084           _size(_memory->size())
00085     {}
00086 
00090     Array(size_t size)
00091         : _memory(new _ArrayMemory<dtype>(size)),
00092           _data(_memory->data()),
00093           _size(_memory->size())
00094     {}
00095 
00099     Array(size_t size, dtype const & val)
00100         : _memory(new _ArrayMemory<dtype>(size, val)),
00101           _data(_memory->data()),
00102           _size(_memory->size())
00103 
00104     {}
00105 
00109     Array(dtype *data, size_t size)
00110         : _memory(new _ArrayMemory<dtype>(data, size)),
00111           _data(_memory->data()),
00112           _size(_memory->size())
00113     {}
00114 
00118     Array(std::vector<dtype> &x)
00119         : _memory(new _ArrayMemory<dtype>(x.data(), x.size())),
00120           _data(_memory->data()),
00121           _size(_memory->size())
00122     {}
00123 
00140     /*
00141      * wrap another array
00142      */
00143     inline void wrap(Array<dtype> x)
00144     {
00145         _memory = x._memory;
00146         _data = x._data;
00147         _size = x._size;
00148     }
00149 
00153     inline dtype *data() { return _data; }
00154     inline dtype const *data() const { return _data; }
00155 
00159     inline size_t size() const { return _size; }
00160 
00164     inline dtype &operator[](const size_t i) { return data()[i]; }
00165     inline dtype const &operator[](const size_t i) const { return data()[i]; }
00166 
00170     typedef dtype * iterator;
00171     typedef dtype const * const_iterator;
00172     inline iterator begin() { return data(); }
00173     inline iterator end() { return data() + size(); }
00174     inline const_iterator begin() const { return data(); }
00175     inline const_iterator end() const { return data() + size(); }
00176 
00177     /*
00178      * Assignment operator: WRAP the data
00179      *
00180      * This is commented because it just duplicates the default assignment operator
00181      *
00182     Array<dtype> &operator=(const Array<dtype> & rhs){
00183         _memory = rhs._memory;
00184         _data = rhs._data;
00185         _size = rhs._size;
00186     }
00187     */
00188 
00189 
00195     inline bool operator==(Array<dtype> const rhs) const
00196     {
00197         return data() == rhs.data() and size() == rhs.size();
00198     }
00199     inline bool operator!=(Array<dtype> const rhs) const
00200     {
00201         return !operator==(rhs);
00202     }
00203 
00204 
00210     Array<dtype> &assign(const Array<dtype> & rhs)
00211     {
00212         //check for self assignment
00213         if ((*this) != rhs) {
00214             if (size() != rhs.size()) {
00215                 throw std::runtime_error("arrays must have the same size during assignment");
00216             }
00217             std::copy(rhs.begin(), rhs.end(), begin());
00218         }
00219         return *this;
00220     }
00221 
00225     Array<dtype> &assign(dtype const & d)
00226     {
00227         std::fill(begin(), end(), d);
00228         return *this;
00229     }
00230 
00231 //    /**
00232 //     * assign each element of the array to be
00233 //     */
00234 //    Array<dtype> &assign(dtype const * const d) {
00235 //        std::copy(rhs.begin(), rhs.end(), begin());
00236 //        for (size_t i = 0; i < size(); ++i) {
00237 //
00238 //        }
00239 //        return *this;
00240 //    }
00241 
00245     Array<dtype> copy() const
00246     {
00247         Array<dtype> newarray(size());
00248         newarray.assign(*this);
00249         return newarray;
00250     }
00251 
00255     inline void free()
00256     {
00257         _memory = std::make_shared<_ArrayMemory<dtype> >();
00258         _data = _memory->data();
00259         _size = _memory->size();
00260     }
00261 
00265     inline bool empty() const
00266     {
00267         return size() == 0;
00268     }
00269 
00270     inline long int reference_count() const
00271     {
00272         return _memory.use_count();
00273     }
00274 
00278     Array<dtype> &operator+=(const Array<dtype> & rhs)
00279     {
00280         if (size() != rhs.size()) {
00281             throw std::runtime_error("operator+=: arrays must have the same size");
00282         }
00283         const_iterator iter = rhs.begin();
00284         for (dtype & val : *this) {
00285             val += *iter;
00286             ++iter;
00287         }
00288         return *this;
00289     }
00290 
00291     Array<dtype> &operator+=(const dtype &rhs)
00292     {
00293         for (dtype & val : (*this)) {
00294             val += rhs;
00295         }
00296         return *this;
00297     }
00298 
00299     Array<dtype> &operator-=(const Array<dtype> & rhs)
00300     {
00301         if (size() != rhs.size()) {
00302             throw std::runtime_error("operator-=: arrays must have the same size");
00303         }
00304         const_iterator iter = rhs.begin();
00305         for (dtype & val : *this) {
00306             val -= *iter;
00307             ++iter;
00308         }
00309         return *this;
00310     }
00311 
00312     Array<dtype> &operator-=(const dtype &rhs)
00313     {
00314         for (dtype & val : (*this)) {
00315             val -= rhs;
00316         }
00317         return *this;
00318     }
00319 
00320     Array<dtype> &operator*=(const Array<dtype> & rhs)
00321     {
00322         if (size() != rhs.size()) {
00323             throw std::runtime_error("operator*=: arrays must have the same size");
00324         }
00325         const_iterator iter = rhs.begin();
00326         for (dtype & val : *this) {
00327             val *= *iter;
00328             ++iter;
00329         }
00330         return *this;
00331     }
00332 
00333     Array<dtype> &operator*=(const dtype &rhs)
00334     {
00335         for (dtype & val : (*this)) {
00336             val *= rhs;
00337         }
00338         return *this;
00339     }
00340     
00341     Array<dtype> operator*(const dtype rhs)
00342     {
00343         Array<dtype> result = this->copy();
00344         return (result *= rhs).copy();
00345     }
00346 
00347     Array<dtype> &operator/=(const Array<dtype> & rhs)
00348     {
00349         if (size() != rhs.size()) {
00350             throw std::runtime_error("operator/=: arrays must have the same size");
00351         }
00352         const_iterator iter = rhs.begin();
00353         for (dtype & val : *this) {
00354             val /= *iter;
00355             ++iter;
00356         }
00357         return *this;
00358     }
00359 
00360     Array<dtype> &operator/=(const  dtype &rhs)
00361     {
00362         for (dtype & val : (*this)) {
00363             val /= rhs;
00364         }
00365         return *this;
00366     }
00367 
00368 /*SOME OTHER ARITHMETIC UTILITIES*/
00369 
00373     dtype sum() const
00374     {
00375         if (empty()) {
00376             throw std::runtime_error("array::sum(): array is empty, can't sum array elements");
00377         }
00378         return std::accumulate(begin(), end(), dtype(0));
00379     }
00380 
00388     dtype prod() const
00389     {
00390         if (empty()) {
00391             throw std::runtime_error("array::prod(): array is empty, can't take product of array elements");
00392         }
00393         return std::accumulate(begin(), end(), dtype(1), std::multiplies<dtype>());
00394     }
00395 
00399     Array<dtype> view(size_t ibegin, size_t iend) const
00400     {
00401         if (iend <= ibegin) {
00402             throw std::invalid_argument("iend must larger than ibegin");
00403         }
00404         if (iend > size()) {
00405             throw std::invalid_argument("iend cannot be larger than array size");
00406         }
00407         Array<dtype> newarray(*this);
00408         newarray._data += ibegin;
00409         newarray._size = iend - ibegin;
00410         return newarray;
00411     }
00412     
00416     dtype get_max() const { return *std::max_element(begin(), end()); }
00417     dtype get_min() const { return *std::min_element(begin(), end()); }
00418 };
00419 
00420 
00421 
00422 // for array printing
00423 template<class dtype>
00424 inline std::ostream &operator<<(std::ostream &out, const Array<dtype> &a)
00425 {
00426     out << "[ ";
00427     for(size_t i = 0; i < a.size(); ++i) {
00428         if(i>0) out << ", ";
00429         out << a[i];
00430     }
00431     out << " ]";
00432     return out;
00433 }
00434 
00438 inline double dot(Array<double> const &v1, Array<double> const &v2)
00439 {
00440   assert(v1.size() == v2.size());
00441   return std::inner_product(v1.begin(), v1.end(), v2.begin(), double(0));
00442 }
00443 
00447 inline double norm(Array<double> const &v)
00448 {
00449   return sqrt(dot(v, v));
00450 }
00451 
00452 template<class T, class U>
00453 Array<T> operator*(const U rhs, const Array<T>& lhs)
00454 {
00455     Array<T> result = lhs.copy();
00456     return (result *= rhs).copy();
00457 }
00458 
00459 } // namespace pele
00460 
00461 #endif // #ifndef PYGMIN_ARRAY_H
 All Classes Namespaces Functions Variables Typedefs