key.h

Go to the documentation of this file.

/*
 This file is part of MADNESS.

 Copyright (C) 2007,2010 Oak Ridge National Laboratory

 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 2 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

 For more information please contact:

 Robert J. Harrison
 Oak Ridge National Laboratory
 One Bethel Valley Road
 P.O. Box 2008, MS-6367

 email: harrisonrj@ornl.gov
 tel:   865-241-3937
 fax:   865-572-0680


 $Id$
 */

#ifndef MADNESS_MRA_KEY_H__INCLUDED
#define MADNESS_MRA_KEY_H__INCLUDED


#include <vector>
#include <madness/mra/power.h>
#include <madness/world/array.h>
#include <madness/world/binfsar.h>
#include <madness/world/worldhash.h>
#include <stdint.h>

namespace madness {

    //     // this has problems when nproc is a large-ish power of 2 such as 256
    //     // and leads to bad data distribution.
    //     static inline unsigned int sdbm(int n, const unsigned char* c, unsigned int sum=0) {
    //         for (int i=0; i<n; ++i) sum = c[i] + (sum << 6) + (sum << 16) - sum;
    //         return sum;
    //     }

    typedef int64_t Translation;
    typedef int Level;

    template<std::size_t NDIM>
    class KeyChildIterator;


    template<std::size_t NDIM>
    class Key {
        friend class KeyChildIterator<NDIM> ;
    private:
        Level n;
        Vector<Translation, NDIM> l;
        hashT hashval;


        // Helper function for operator <
        int
        encode(int dig) const {
            int retval = 0;
            for (std::size_t j = 0; j < NDIM; ++j) {
                // retval += ((l[j]/2^{n-1-dig}) mod 2) * 2^j
                retval += ((l[j] >> (n - 1 - dig)) % 2) << j;
            }
            return retval;
        }

        // Helper function for (Level, Translation) constructor
        Vector<Translation, NDIM>
        decode(Level level, Translation k) const {
            Vector<Translation, NDIM> L(0);
            int twotoD = power<static_cast<int>(NDIM)> ();
            int powr = 1, divisor = 2;
            for (Level i = 0; i < level; ++i) {
                Translation r = k % twotoD;
                for (int j = 0; j < NDIM; ++j) {
                    L[NDIM - j - 1] += (r % divisor) * powr;
                    r /= divisor;
                }
                k /= twotoD;
                powr *= 2;
            }
            return L;
        }
    public:
        Key() {
        }

        Key(Level n, const Vector<Translation, NDIM>& l) :
                n(n), l(l) {
            rehash();
        }

        Key(int n) :
                n(n), l(0) {
            rehash();
        }

        Key(Level n, Translation p) :
                n(n) {
            l = decode(n, p);
            rehash();
        }

        Key(const int n, const int l0) : n(n) {
            MADNESS_ASSERT(NDIM==1);
            l=Vector<Translation, NDIM>(l0);
            rehash();
        }

        Key(const int n, const int l0, const int l1, const int l2) : n(n) {
            MADNESS_ASSERT(NDIM==3);
            l=Vector<Translation, NDIM>(0);
            l[0]=l0; l[1]=l1; l[2]=l2;
            rehash();
        }

        static Key<NDIM>
        invalid() {
            return Key<NDIM> (-1);
        }

        bool
        is_invalid() const {
            return n == -1;
        }

        bool
        is_valid() const {
            return n != -1;
        }

        bool
        operator==(const Key& other) const {
            if (hashval != other.hashval)
                return false;
            if (n != other.n)
                return false;
            bool result = l == other.l;
            if (result && hashval != other.hashval) {
                print("!!  keys same but hash is different", hashval,
                      other.hashval, *this, other);
                MADNESS_EXCEPTION("Tell HQI not RJ3!",0);
            }
            return result;
        }

        bool
        operator!=(const Key& other) const {
            return !(*this == other);
        }

        bool
        operator<(const Key& other) const {
            if (*this == other)
                return false; // I am not less than self
            Level nmin;
            bool retval = false;

            if (this->n > other.n) {
                nmin = other.n;
                retval = true;
            }
            else {
                nmin = this->n;
            }

            for (Level i = 0; i < nmin; ++i) {
                int tthis = this->encode(i), tother = other.encode(i);
                if (tthis != tother) {
                    return (tthis < tother);
                }
            }
            return retval;
        }

        inline hashT
        hash() const {
            return hashval;
        }

        // template<typename Archive>
        // inline void
        // serialize(Archive& ar) {
        //     ar & archive::wrap((unsigned char*) this, sizeof(*this));
        // }

        Level
        level() const {
            return n;
        }

        const Vector<Translation, NDIM>&
        translation() const {
            return l;
        }

        uint64_t
        distsq() const {
            uint64_t dist = 0;
            for (std::size_t d = 0; d < NDIM; ++d) {
                dist += l[d] * l[d];
            }
            return dist;
        }


        Key
        parent(int generation = 1) const {
            Vector<Translation, NDIM> pl;
            if (generation > n)
                generation = n;
            for (std::size_t i = 0; i < NDIM; ++i)
                pl[i] = l[i] >> generation;
            return Key(n - generation, pl);
        }


        bool
        is_child_of(const Key& key) const {
            if (this->n < key.n) {
                return false; // I can't be child of something lower on the tree
            }
            else if (this->n == key.n) {
                return (*this == key); // I am child of myself
            }
            else {
                Level dn = this->n - key.n;
                Key mama = this->parent(dn);
                return (mama == key);
            }
        }


        bool
        is_parent_of(const Key& key) const {
            return (key.is_child_of(*this));
        }


        bool
        is_neighbor_of(const Key& key, const std::vector<bool>& bperiodic) const {
          Translation dist = 0;
          Translation TWON1 = (Translation(1)<<n) - 1;
                for (std::size_t i=0; i<NDIM; ++i)
                {
                  Translation ll = std::abs(l[i] - key.l[i]);
                  if (bperiodic[i] && ll==TWON1) ll=1;
                  dist = std::max(dist, ll);
                }
                return (dist <= 1);
        }


        Key neighbor(const Key<NDIM>& disp) const {
            Vector<Translation,NDIM> l = this->translation()+disp.translation();
            return Key(this->level(),l);
        }


        bool thisKeyContains(const Vector<double,NDIM>& x, const unsigned int& dim0,
                        const unsigned int& dim1) const {

                // it's sufficient if one single dimension is out
                bool contains=true;
                const double twotoN = std::pow(2.0,double(n));
                MADNESS_ASSERT(dim0<NDIM and dim1<NDIM);

                for (unsigned int i=0; i<NDIM; i++ ) {

                        // check bounds
                        MADNESS_ASSERT((x[i]>=0.0) and (x[i]<=1.0));

                        // leave these two dimensions out
                        if ((i==dim0) or (i==dim1)) continue;

                        const int ll=int (x[i]*twotoN);
                        if (not (l[i]==ll)) contains=false;
                }
                return contains;
        }

        template<std::size_t LDIM, std::size_t KDIM>
        void break_apart(Key<LDIM>& key1, Key<KDIM>& key2) const {

            // if LDIM==NDIM the 2nd key will be constructed empty
            MADNESS_ASSERT((LDIM+KDIM==NDIM) or (LDIM==NDIM));
            Vector<Translation, LDIM> l1;
            Vector<Translation, KDIM> l2;
            for (int i=0; i<static_cast<int>(LDIM); ++i) {
                l1[i]=l[i];
            }
            for (size_t i=LDIM; i<NDIM; ++i) {
                l2[i-LDIM]=l[i];
            }
            key1=Key<LDIM>(n,l1);
            key2=Key<KDIM>(n,l2);
        }

        template<std::size_t LDIM>
        Key<NDIM+LDIM> merge_with(const Key<LDIM>& rhs) const {
            Vector<Translation,NDIM+LDIM> t;
            for (int i=0; i<static_cast<int>(NDIM); ++i) t[i]     =this->l[i];
            for (int i=0; i<static_cast<int>(LDIM); ++i) t[NDIM+i]=rhs.translation()[i];
            return Key<NDIM+LDIM>(rhs.level(),t);
        }


        bool is_farther_out_than(const Key<NDIM>& other) const {
                for (size_t i=0; i<NDIM; ++i) {
                        if ((other.translation()[i]>0) and (other.translation()[i]>l[i])) return false;
                        if ((other.translation()[i]<0) and (other.translation()[i]<l[i])) return false;
                }
                return true;
        }


        void
        rehash() {
            //hashval = sdbm(sizeof(n)+sizeof(l), (unsigned char*)(&n));
            // default hash is still best
            hashval = hash_value(l);
            hash_combine(hashval, n);
        }
    };

    template<std::size_t NDIM>
    std::ostream&
    operator<<(std::ostream& s, const Key<NDIM>& key) {
        s << "(" << key.level() << "," << key.translation() << ")";
        return s;
    }


    template<size_t NDIM>
    Key<NDIM> displacement(const Key<NDIM>& source, const Key<NDIM>& target) {
        MADNESS_ASSERT(source.level()==target.level());
        const Vector<Translation,NDIM> l = target.translation()-source.translation();
        return Key<NDIM>(source.level(),l);
    }




    template<std::size_t NDIM>
    class KeyChildIterator {
        Key<NDIM> parent;
        Key<NDIM> child;
        Vector<Translation, NDIM> p;
        bool finished;

    public:
        KeyChildIterator() :
                p(0), finished(true) {
        }

        KeyChildIterator(const Key<NDIM>& parent) :
                parent(parent), child(parent.n + 1, parent.l * 2), p(0),
                finished(false) {
        }

        KeyChildIterator&
        operator++() {
            if (finished)
                return *this;
            std::size_t i;
            for (i = 0; i < NDIM; ++i) {
                if (p[i] == 0) {
                    ++(p[i]);
                    ++(child.l[i]);
                    for (std::size_t j = 0; j < i; ++j) {
                        --(p[j]);
                        --(child.l[j]);
                    }
                    break;
                }
            }
            finished = (i == NDIM);
            child.rehash();
            return *this;
        }

        operator bool() const {
            return !finished;
        }

        template<typename Archive>
        inline void
        serialize(Archive& ar) {
            ar & archive::wrap((unsigned char*) this, sizeof(*this));
        }

        inline const Key<NDIM>&
        key() const {
            return child;
        }
    };

    template<std::size_t NDIM, typename opT>
    inline void
    foreach_child(const Key<NDIM>& parent, opT& op) {
        for (KeyChildIterator<NDIM>
                it(parent); it; ++it)
            op(it.key());
    }

    template<std::size_t NDIM, typename objT>
    inline void
    foreach_child(const Key<NDIM>& parent, objT* obj, void
                  (objT::*memfun)(const Key<NDIM>&)) {
        for (KeyChildIterator<NDIM>
                it(parent); it; ++it)
            (obj ->* memfun)(it.key());
    }

    namespace archive {

        // For efficiency serialize opaque so is just one memcpy, but
        // when reading from external storage rehash() so that we
        // can read data even if hash algorithm/function has changed.

        template <class Archive, std::size_t NDIM>
        struct ArchiveLoadImpl< Archive, Key<NDIM> > {
            static void load(const Archive& ar, Key<NDIM>& t) {
                ar & archive::wrap((unsigned char*) &t, sizeof(t));
            }
        };

        template <std::size_t NDIM>
        struct ArchiveLoadImpl< BinaryFstreamInputArchive, Key<NDIM> > {
            static void load(const BinaryFstreamInputArchive& ar, Key<NDIM>& t) {
                ar & archive::wrap((unsigned char*) &t, sizeof(t));
                t.rehash(); // <<<<<<<<<< This is the point
            }
        };

        template <class Archive, std::size_t NDIM>
        struct ArchiveStoreImpl< Archive, Key<NDIM> > {
            static void store(const Archive& ar, const Key<NDIM>& t) {
                ar & archive::wrap((unsigned char*) &t, sizeof(t));
            }
        };
    }

}

#endif // MADNESS_MRA_KEY_H__INCLUDED