1 /* 2 This file is part of MADNESS. 3 4 Copyright (C) 2007,2010 Oak Ridge National Laboratory 5 6 This program is free software; you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation; either version 2 of the License, or 9 (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with this program; if not, write to the Free Software 18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 19 20 For more information please contact: 21 22 Robert J. Harrison 23 Oak Ridge National Laboratory 24 One Bethel Valley Road 25 P.O. Box 2008, MS-6367 26 27 email: harrisonrj@ornl.gov 28 tel: 865-241-3937 29 fax: 865-572-0680 30 */ 31 32 #ifndef MADNESS_MRA_KEY_H__INCLUDED 33 #define MADNESS_MRA_KEY_H__INCLUDED 34 35 /// \file key.h 36 /// \brief Multidimension Key for MRA tree and associated iterators 37 38 #include <vector> 39 #include <madness/mra/power.h> 40 #include <madness/world/vector.h> 41 #include <madness/world/binary_fstream_archive.h> 42 #include <madness/world/worldhash.h> 43 #include <stdint.h> 44 45 namespace madness { 46 47 // // this has problems when nproc is a large-ish power of 2 such as 256 48 // // and leads to bad data distribution. 49 // static inline unsigned int sdbm(int n, const unsigned char* c, unsigned int sum=0) { 50 // for (int i=0; i<n; ++i) sum = c[i] + (sum << 6) + (sum << 16) - sum; 51 // return sum; 52 // } 53 54 typedef int64_t Translation; 55 typedef int Level; 56 57 template<std::size_t NDIM> 58 class KeyChildIterator; 59 60 /// Key is the index for a node of the 2^NDIM-tree 61 62 /// See KeyChildIterator for facile generation of children, 63 /// and foreach_child(parent,op) for facile applicaiton of operators 64 /// to child keys. 65 template<std::size_t NDIM> 66 class Key { 67 friend class KeyChildIterator<NDIM> ; 68 private: 69 Level n; 70 Vector<Translation, NDIM> l; 71 hashT hashval; 72 73 74 public: 75 /// Default constructor makes an \em uninitialized key Key()76 Key() {} 77 78 /// Constructor with given n, l Key(Level n,const Vector<Translation,NDIM> & l)79 Key(Level n, const Vector<Translation, NDIM>& l) : n(n), l(l) 80 { 81 rehash(); 82 } 83 84 /// Constructor with given n and l=0 Key(int n)85 Key(int n) : n(n), l(0) 86 { 87 rehash(); 88 } 89 90 // /// easy constructor ... UGH !!!!!!!!!!!!!!!!!!!!!! 91 // Key(const int n, const int l0, const int l1, const int l2) : n(n) { 92 // MADNESS_ASSERT(NDIM==3); 93 // l=Vector<Translation, NDIM>(0); 94 // l[0]=l0; l[1]=l1; l[2]=l2; 95 // rehash(); 96 // } 97 98 /// Returns an invalid key invalid()99 static Key<NDIM> invalid() { 100 return Key<NDIM> (-1); 101 } 102 103 /// Checks if a key is invalid is_invalid()104 bool is_invalid() const { 105 return n == -1; 106 } 107 108 /// Checks if a key is valid is_valid()109 bool is_valid() const { 110 return n != -1; 111 } 112 113 /// Equality test 114 bool operator==(const Key& other) const { 115 if (hashval != other.hashval) return false; 116 if (n != other.n) return false; 117 for (unsigned int i=0; i<NDIM; i++) 118 if (l[i] != other.l[i]) return false; 119 return true; // everything is equal 120 } 121 122 bool operator!=(const Key& other) const { 123 return !(*this == other); 124 } 125 126 /// Comparison operator less than to enable storage in STL map 127 bool operator<(const Key& other) const { 128 if (hashval < other.hashval) return true; 129 if (hashval > other.hashval) return false; 130 131 if (n < other.n) return true; 132 if (n > other.n) return false; 133 134 for (unsigned int i=0; i<NDIM; i++) { 135 if (l[i] < other.l[i]) return true; 136 if (l[i] > other.l[i]) return false; 137 } 138 139 return false; // everything is equal 140 } 141 142 inline hashT hash()143 hash() const { 144 return hashval; 145 } 146 147 // template<typename Archive> 148 // inline void 149 // serialize(Archive& ar) { 150 // ar & archive::wrap((unsigned char*) this, sizeof(*this)); 151 // } 152 153 Level level()154 level() const { 155 return n; 156 } 157 158 const Vector<Translation, NDIM>& translation()159 translation() const { 160 return l; 161 } 162 163 uint64_t distsq()164 distsq() const { 165 uint64_t dist = 0; 166 for (std::size_t d = 0; d < NDIM; ++d) { 167 dist += l[d] * l[d]; 168 } 169 return dist; 170 } 171 172 /// Returns the key of the parent 173 174 /// Default is the immediate parent (generation=1). To get 175 /// the grandparent use generation=2, and similarly for 176 /// great-grandparents. 177 /// 178 /// !! If there is no such parent it quietly returns the 179 /// closest match (which may be self if this is the top of the 180 /// tree). 181 Key 182 parent(int generation = 1) const { 183 Vector<Translation, NDIM> pl; 184 if (generation > n) 185 generation = n; 186 for (std::size_t i = 0; i < NDIM; ++i) 187 pl[i] = l[i] >> generation; 188 return Key(n - generation, pl); 189 } 190 191 192 bool is_child_of(const Key & key)193 is_child_of(const Key& key) const { 194 if (this->n < key.n) { 195 return false; // I can't be child of something lower on the tree 196 } 197 else if (this->n == key.n) { 198 return (*this == key); // I am child of myself 199 } 200 else { 201 Level dn = this->n - key.n; 202 Key mama = this->parent(dn); 203 return (mama == key); 204 } 205 } 206 207 208 bool is_parent_of(const Key & key)209 is_parent_of(const Key& key) const { 210 return (key.is_child_of(*this)); 211 } 212 213 /// Assuming keys are at the same level, returns true if displaced by no more than 1 in any direction 214 215 /// Assumes key and this are at the same level 216 bool is_neighbor_of(const Key & key,const std::vector<bool> & bperiodic)217 is_neighbor_of(const Key& key, const std::vector<bool>& bperiodic) const { 218 Translation dist = 0; 219 Translation TWON1 = (Translation(1)<<n) - 1; 220 for (std::size_t i=0; i<NDIM; ++i) 221 { 222 Translation ll = std::abs(l[i] - key.l[i]); 223 if (bperiodic[i] && ll==TWON1) ll=1; 224 dist = std::max(dist, ll); 225 } 226 return (dist <= 1); 227 } 228 229 /// given a displacement, generate a neighbor key; ignore boundary conditions and disp's level 230 231 /// @param[in] disp the displacement 232 /// @return a new key neighbor(const Key<NDIM> & disp)233 Key neighbor(const Key<NDIM>& disp) const { 234 Vector<Translation,NDIM> l = this->translation()+disp.translation(); 235 return Key(this->level(),l); 236 } 237 238 239 /// check if this MultiIndex contains point x, disregarding these two dimensions thisKeyContains(const Vector<double,NDIM> & x,const unsigned int & dim0,const unsigned int & dim1)240 bool thisKeyContains(const Vector<double,NDIM>& x, const unsigned int& dim0, 241 const unsigned int& dim1) const { 242 243 // it's sufficient if one single dimension is out 244 bool contains=true; 245 const double twotoN = std::pow(2.0,double(n)); 246 MADNESS_ASSERT(dim0<NDIM and dim1<NDIM); 247 248 for (unsigned int i=0; i<NDIM; i++ ) { 249 250 // check bounds 251 MADNESS_ASSERT((x[i]>=0.0) and (x[i]<=1.0)); 252 253 // leave these two dimensions out 254 if ((i==dim0) or (i==dim1)) continue; 255 256 const int ll=int (x[i]*twotoN); 257 if (not (l[i]==ll)) contains=false; 258 } 259 return contains; 260 } 261 262 /// break key into two low-dimensional keys 263 template<std::size_t LDIM, std::size_t KDIM> break_apart(Key<LDIM> & key1,Key<KDIM> & key2)264 void break_apart(Key<LDIM>& key1, Key<KDIM>& key2) const { 265 266 // if LDIM==NDIM the 2nd key will be constructed empty 267 MADNESS_ASSERT((LDIM+KDIM==NDIM) or (LDIM==NDIM)); 268 Vector<Translation, LDIM> l1; 269 Vector<Translation, KDIM> l2; 270 for (int i=0; i<static_cast<int>(LDIM); ++i) { 271 l1[i]=l[i]; 272 } 273 for (size_t i=LDIM; i<NDIM; ++i) { 274 l2[i-LDIM]=l[i]; 275 } 276 key1=Key<LDIM>(n,l1); 277 key2=Key<KDIM>(n,l2); 278 } 279 280 /// merge with other key (ie concatenate), use level of rhs, not of this 281 template<std::size_t LDIM> merge_with(const Key<LDIM> & rhs)282 Key<NDIM+LDIM> merge_with(const Key<LDIM>& rhs) const { 283 Vector<Translation,NDIM+LDIM> t; 284 for (int i=0; i<static_cast<int>(NDIM); ++i) t[i] =this->l[i]; 285 for (int i=0; i<static_cast<int>(LDIM); ++i) t[NDIM+i]=rhs.translation()[i]; 286 return Key<NDIM+LDIM>(rhs.level(),t); 287 } 288 289 /// return if the other key is pointing in the same direction and is farther out 290 291 /// unlike in distsq() the direction is taken into account, and other must be 292 /// longer than this in each dimension 293 /// @param[in] other a key 294 /// @return if other is farther out is_farther_out_than(const Key<NDIM> & other)295 bool is_farther_out_than(const Key<NDIM>& other) const { 296 for (size_t i=0; i<NDIM; ++i) { 297 if ((other.translation()[i]>0) and (other.translation()[i]>l[i])) return false; 298 if ((other.translation()[i]<0) and (other.translation()[i]<l[i])) return false; 299 } 300 return true; 301 } 302 303 304 /// Recomputes hashval ... presently only done when reading from external storage 305 void rehash()306 rehash() { 307 //hashval = sdbm(sizeof(n)+sizeof(l), (unsigned char*)(&n)); 308 // default hash is still best 309 310 hashval = hash_value(l); 311 hash_combine(hashval, n); 312 } 313 }; 314 315 template<std::size_t NDIM> 316 std::ostream& 317 operator<<(std::ostream& s, const Key<NDIM>& key) { 318 s << "(" << key.level() << "," << key.translation() << ")"; 319 return s; 320 } 321 322 /// given a source and a target, return the displacement in translation 323 324 /// @param[in] source the source key 325 /// @param[in] target the target key 326 /// @return disp such that target = source + disp 327 template<size_t NDIM> displacement(const Key<NDIM> & source,const Key<NDIM> & target)328 Key<NDIM> displacement(const Key<NDIM>& source, const Key<NDIM>& target) { 329 MADNESS_ASSERT(source.level()==target.level()); 330 const Vector<Translation,NDIM> l = target.translation()-source.translation(); 331 return Key<NDIM>(source.level(),l); 332 } 333 334 335 336 /// Iterates in lexical order thru all children of a key 337 338 /// Example usage: 339 /// \code 340 /// for (KeyChildIterator<NDIM> it(key); it; ++it) print(it.key()); 341 /// \endcode 342 template<std::size_t NDIM> 343 class KeyChildIterator { 344 Key<NDIM> parent; 345 Key<NDIM> child; 346 Vector<Translation, NDIM> p; 347 bool finished; 348 349 public: KeyChildIterator()350 KeyChildIterator() : 351 p(0), finished(true) { 352 } 353 KeyChildIterator(const Key<NDIM> & parent)354 KeyChildIterator(const Key<NDIM>& parent) : 355 parent(parent), child(parent.n + 1, parent.l * 2), p(0), 356 finished(false) { 357 } 358 359 /// Pre-increment of an iterator (i.e., ++it) 360 KeyChildIterator& 361 operator++() { 362 if (finished) 363 return *this; 364 std::size_t i; 365 for (i = 0; i < NDIM; ++i) { 366 if (p[i] == 0) { 367 ++(p[i]); 368 ++(child.l[i]); 369 for (std::size_t j = 0; j < i; ++j) { 370 --(p[j]); 371 --(child.l[j]); 372 } 373 break; 374 } 375 } 376 finished = (i == NDIM); 377 child.rehash(); 378 return *this; 379 } 380 381 /// True if iterator is not at end 382 operator bool() const { 383 return !finished; 384 } 385 386 template<typename Archive> 387 inline void serialize(Archive & ar)388 serialize(Archive& ar) { 389 ar & archive::wrap((unsigned char*) this, sizeof(*this)); 390 } 391 392 /// Returns the key of the child 393 inline const Key<NDIM>& key()394 key() const { 395 return child; 396 } 397 }; 398 399 /// Applies op(key) to each child key of parent 400 template<std::size_t NDIM, typename opT> 401 inline void foreach_child(const Key<NDIM> & parent,opT & op)402 foreach_child(const Key<NDIM>& parent, opT& op) { 403 for (KeyChildIterator<NDIM> 404 it(parent); it; ++it) 405 op(it.key()); 406 } 407 408 /// Applies member function of obj to each child key of parent 409 template<std::size_t NDIM, typename objT> 410 inline void foreach_child(const Key<NDIM> & parent,objT * obj,void (objT::* memfun)(const Key<NDIM> &))411 foreach_child(const Key<NDIM>& parent, objT* obj, void 412 (objT::*memfun)(const Key<NDIM>&)) { 413 for (KeyChildIterator<NDIM> 414 it(parent); it; ++it) 415 (obj ->* memfun)(it.key()); 416 } 417 418 namespace archive { 419 420 // For efficiency serialize opaque so is just one memcpy, but 421 // when reading from external storage rehash() so that we 422 // can read data even if hash algorithm/function has changed. 423 424 template <class Archive, std::size_t NDIM> 425 struct ArchiveLoadImpl< Archive, Key<NDIM> > { 426 static void load(const Archive& ar, Key<NDIM>& t) { 427 ar & archive::wrap((unsigned char*) &t, sizeof(t)); 428 } 429 }; 430 431 template <std::size_t NDIM> 432 struct ArchiveLoadImpl< BinaryFstreamInputArchive, Key<NDIM> > { 433 static void load(const BinaryFstreamInputArchive& ar, Key<NDIM>& t) { 434 ar & archive::wrap((unsigned char*) &t, sizeof(t)); 435 t.rehash(); // <<<<<<<<<< This is the point 436 } 437 }; 438 439 template <class Archive, std::size_t NDIM> 440 struct ArchiveStoreImpl< Archive, Key<NDIM> > { 441 static void store(const Archive& ar, const Key<NDIM>& t) { 442 ar & archive::wrap((unsigned char*) &t, sizeof(t)); 443 } 444 }; 445 } 446 447 } 448 449 #endif // MADNESS_MRA_KEY_H__INCLUDED 450 451