mcpele
1.0.0
The Monte Carlo Python Energy Landscape Explorer
|
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 00342 Array<dtype> &operator/=(const Array<dtype> & rhs) 00343 { 00344 if (size() != rhs.size()) { 00345 throw std::runtime_error("operator/=: arrays must have the same size"); 00346 } 00347 const_iterator iter = rhs.begin(); 00348 for (dtype & val : *this) { 00349 val /= *iter; 00350 ++iter; 00351 } 00352 return *this; 00353 } 00354 00355 Array<dtype> &operator/=(const dtype &rhs) 00356 { 00357 for (dtype & val : (*this)) { 00358 val /= rhs; 00359 } 00360 return *this; 00361 } 00362 00363 /*SOME OTHER ARITHMETIC UTILITIES*/ 00364 00368 dtype sum() const 00369 { 00370 if (empty()) { 00371 throw std::runtime_error("array::sum(): array is empty, can't sum array elements"); 00372 } 00373 return std::accumulate(begin(), end(), dtype(0)); 00374 } 00375 00383 dtype prod() const 00384 { 00385 if (empty()) { 00386 throw std::runtime_error("array::prod(): array is empty, can't take product of array elements"); 00387 } 00388 return std::accumulate(begin(), end(), dtype(1), std::multiplies<dtype>()); 00389 } 00390 00394 Array<dtype> view(size_t ibegin, size_t iend) const 00395 { 00396 if (iend <= ibegin) { 00397 throw std::invalid_argument("iend must larger than ibegin"); 00398 } 00399 if (iend > size()) { 00400 throw std::invalid_argument("iend cannot be larger than array size"); 00401 } 00402 Array<dtype> newarray(*this); 00403 newarray._data += ibegin; 00404 newarray._size = iend - ibegin; 00405 return newarray; 00406 } 00407 00411 dtype get_max() const { return *std::max_element(begin(), end()); } 00412 dtype get_min() const { return *std::min_element(begin(), end()); } 00413 }; 00414 00415 00416 00417 // for array printing 00418 template<class dtype> 00419 inline std::ostream &operator<<(std::ostream &out, const Array<dtype> &a) 00420 { 00421 out << "[ "; 00422 for(size_t i = 0; i < a.size(); ++i) { 00423 if(i>0) out << ", "; 00424 out << a[i]; 00425 } 00426 out << " ]"; 00427 return out; 00428 } 00429 00433 inline double dot(Array<double> const &v1, Array<double> const &v2) 00434 { 00435 assert(v1.size() == v2.size()); 00436 return std::inner_product(v1.begin(), v1.end(), v2.begin(), double(0)); 00437 } 00438 00442 inline double norm(Array<double> const &v) 00443 { 00444 return sqrt(dot(v, v)); 00445 } 00446 00447 } // namespace pele 00448 00449 #endif // #ifndef PYGMIN_ARRAY_H