.. _program_listing_file_src_array_latvec.hpp:

Program Listing for File array_latvec.hpp
=========================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_array_latvec.hpp>` (``src/array_latvec.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   /* This file is part of brille.
   
   Copyright © 2019,2020 Greg Tucker <greg.tucker@stfc.ac.uk>
   
   brille is free software: you can redistribute it and/or modify it under the
   terms of the GNU Affero General Public License as published by the Free
   Software Foundation, either version 3 of the License, or (at your option)
   any later version.
   
   brille 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 Affero General Public License for more details.
   
   You should have received a copy of the GNU Affero General Public License
   along with brille. If not, see <https://www.gnu.org/licenses/>.            */
   
   #ifndef BRILLE_LATVEC_CLASS_H_
   #define BRILLE_LATVEC_CLASS_H_
   
   #include <typeinfo> // for std::bad_cast
   #include <exception>
   #include "lattice.hpp"
   #include "array2.hpp"
   
   template<class T>
   using bArray = brille::Array2<T>;
   
   namespace brille {
   
   
   class LatVec{};
   template<class T> class LDVec;
   template<class T> class LQVec;
   
   
   template<class... T> struct ArrayTraits{
     static constexpr bool array = false;
     static constexpr bool latvec = false;
   };
   #ifndef DOXYGEN_SHOULD_SKIP_THIS
   template<class T> struct ArrayTraits<bArray<T>>{
     static constexpr bool array = true;
     static constexpr bool latvec = false;
   };
   template<class T> struct ArrayTraits<LDVec<T>>{
     static constexpr bool array = true;
     static constexpr bool latvec = true;
   };
   template<class T> struct ArrayTraits<LQVec<T>>{
     static constexpr bool array = true;
     static constexpr bool latvec = true;
   };
   #endif
   
   
   template<class T, template<class> class A>
   inline constexpr bool isArray = ArrayTraits<A<T>>::array;
   template<class T, template<class> class A>
   inline constexpr bool isLatVec = ArrayTraits<A<T>>::latvec;
   template<class T, template<class> class A>
   inline constexpr bool isBareArray = isArray<T,A> && !isLatVec<T,A>;
   
   template<class T, template<class> class A, class R, template<class> class B>
   inline constexpr bool bothArrays = isArray<T,A> && isArray<R,B>;
   template<class T, template<class> class A, class R, template<class> class B>
   inline constexpr bool bareArrays = isBareArray<T,A> && isBareArray<R,B>;
   template<class T, template<class> class A, class R, template<class> class B>
   inline constexpr bool bothLatVecs = isLatVec<T,A> && isLatVec<R,B>;
   
   template<class T>
   bool equal_shapes(const std::vector<T>& a, const std::vector<T>&b){
     bool ok = a.size() == b.size();
     if (ok) ok = std::equal(a.begin(), a.end(), b.begin());
     return ok;
   }
   template<class T, size_t Nel>
   bool equal_shapes(const std::array<T,Nel>& a, const std::array<T,Nel>&b){
     return std::equal(a.begin(), a.end(), b.begin());
   }
   
   #define LVEC_SCALAR_INPLACE_OP_DEF(L,X) \
   L <T>&\
   operator X (const T& v){\
     for (auto x: this->subItr()) this->_data[this->s2l_d(x)] X v;\
     return *this;\
   }
   #define LVEC_ARRAY_INPLACE_OP_DEF(L,X) \
   template<template<class> class A>\
   std::enable_if_t<isArray<T,A>, L <T>>&\
   operator X (const A<T>& b){\
     auto itr = this->broadcastItr(b.shape());\
     if (!equal_shapes(itr.shape(), this->shape())) throw std::runtime_error("Incompatible shaped Array for binary operation");\
     for (auto [ox, ax, bx]: itr) this->_data[this->s2l_d(ax)] X b[bx];\
     return *this;\
   }
   
   #define LVEC_INIT_INT(N,L,X) N(const L &lat, const X &n)\
   : bArray<T>(static_cast<brille::ind_t>(n),3), lattice(lat)\
   {}
   
   template<class T>
   class LDVec: public LatVec, public bArray<T>{
     Direct lattice;
   public:
     // default constructor for zero three-vectors:
     explicit LDVec(const Direct& lat=Direct())
     : bArray<T>(0,3), lattice(lat)
     {
     }
     LVEC_INIT_INT(LDVec,Direct,int)
     LVEC_INIT_INT(LDVec,Direct,long)
     LVEC_INIT_INT(LDVec,Direct,long long)
     LVEC_INIT_INT(LDVec,Direct,unsigned)
     LVEC_INIT_INT(LDVec,Direct,unsigned long)
     LVEC_INIT_INT(LDVec,Direct,unsigned long long)
     template<typename... Args> LDVec(const Direct& lat, Args... args)
     : bArray<T>(args...), lattice(lat)
     {
       this->check_array();
     }
     template<class R>
     LDVec(const LDVec<R>& other)
     : bArray<T>(other.get_hkl()), lattice(other.get_lattice())
     {
     }
   
   
     Direct get_lattice() const { return lattice; }
     template<class R> bool samelattice(const LDVec<R> *vec) const { return lattice.issame(vec->get_lattice()); }
     template<class R> bool samelattice(const LQVec<R> *)    const { return false; }
     template<class R> bool starlattice(const LDVec<R> *)    const { return false; }
     template<class R> bool starlattice(const LQVec<R> *vec) const { return lattice.isstar(vec->get_lattice()); }
     template<class R> bool samelattice(const LDVec<R> &vec) const { return lattice.issame(vec.get_lattice()); }
     template<class R> bool samelattice(const LQVec<R> &)    const { return false; }
     template<class R> bool starlattice(const LDVec<R> &)    const { return false; }
     template<class R> bool starlattice(const LQVec<R> &vec) const { return lattice.isstar(vec.get_lattice()); }
     // extract overloads preserving lattice information
     template<typename... A> LDVec<T> view(A... args) const;
     template<typename... A> LDVec<T> extract(A... args) const;
     bArray<T> get_hkl() const;
     bArray<double> get_xyz() const;
     LQVec<double> star() const;
   
     double dot(const size_t i, const size_t j) const;
     double norm(const size_t i) const { return sqrt(this->dot(i,i)); }
     LDVec<double> cross(const size_t i, const size_t j) const;
   
     // LDVec<T>& operator-=(const LDVec<T>& av);
     // LDVec<T>& operator+=(const LDVec<T>& av);
   
     LVEC_SCALAR_INPLACE_OP_DEF(LDVec,+=)
     LVEC_SCALAR_INPLACE_OP_DEF(LDVec,-=)
     LVEC_SCALAR_INPLACE_OP_DEF(LDVec,*=)
     LVEC_SCALAR_INPLACE_OP_DEF(LDVec,/=)
   
     LVEC_ARRAY_INPLACE_OP_DEF(LDVec,+=)
     LVEC_ARRAY_INPLACE_OP_DEF(LDVec,-=)
     LVEC_ARRAY_INPLACE_OP_DEF(LDVec,*=)
     LVEC_ARRAY_INPLACE_OP_DEF(LDVec,/=)
   
   
     template<class R>
     void binary_operation_check(const LDVec<R>& b) const{
       assert(this->samelattice(b));
     }
     template<class R>
     void binary_operation_check(const LQVec<R>& b) const{
       assert(this->starlattice(b));
     }
     template<class R, template<class> class A,
     typename=typename std::enable_if<!std::is_base_of<LatVec,A<R>>::value>::type>
     void binary_operation_check(const A<R>&) const {}
   
     template<class R>
     bool is(const LDVec<R>& that){
       return (this->samelattice(that) && this->bArray<T>::is(that));
     }
     LDVec<int> round() const {
       return LDVec<int>(this->lattice, this->bArray<T>::round());
     }
     LDVec<int> floor() const {
       return LDVec<int>(this->lattice, this->bArray<T>::floor());
     }
     LDVec<int> ceil() const {
       return LDVec<int>(this->lattice, this->bArray<T>::ceil());
     }
   
     LDVec<T> decouple() {
       return LDVec<T>(lattice, this->bArray<T>::decouple());
     }
   protected:
     void check_array();
   };
   
   
   template<class T>
   class LQVec:  public LatVec, public bArray<T>{
     Reciprocal lattice;
   public:
     // default constructor for zero three-vectors:
     explicit LQVec(const Reciprocal& lat=Reciprocal())
     : bArray<T>(0,3), lattice(lat)
     {
     }
     LVEC_INIT_INT(LQVec,Reciprocal,int)
     LVEC_INIT_INT(LQVec,Reciprocal,long)
     LVEC_INIT_INT(LQVec,Reciprocal,long long)
     LVEC_INIT_INT(LQVec,Reciprocal,unsigned)
     LVEC_INIT_INT(LQVec,Reciprocal,unsigned long)
     LVEC_INIT_INT(LQVec,Reciprocal,unsigned long long)
     template<typename... Args>
     LQVec(const Reciprocal& lat, Args... args)
     : bArray<T>(args...), lattice(lat)
     {
       this->check_array();
     }
   
     template<class R>
     LQVec(const LQVec<R>& other)
     : bArray<T>(other.get_hkl()), lattice(other.get_lattice())
     {
     }
   
     // LQVec(const LQVec<T>& other, const bool m, const bool d): bArray<T>(other.get_hkl(), m, d), lattice{other.get_lattice()} {}
     // LQVec(const LQVec<T>& other, const bool m=false, const bool d=false): bArray<T>(other.get_hkl(),m,d), lattice(other.get_lattice()) {}
   
     // //! Type conversion assignment:
     // template<class R> LQVec<T>& operator=(const LQVec<R>& other){
     //   this->lattice = other.get_lattice();
     //   *this = bArray<T>(other.get_hkl());
     //   return *this;
     // }
     // LQVec<T>& operator=(const LQVec<T>& other){
     //   if (this != &other){ // do nothing if called by, e.g., a = a;
     //     this->lattice = other.get_lattice();
     //     this->bArray<T>::operator=(other);
     //   }
     //   return *this;
     // }
     // LQVec<T>& operator=(LQVec<T>&& other) noexcept {
     //
     // }
   
     Reciprocal get_lattice() const { return lattice; }
     template<class R> bool samelattice(const LQVec<R> *vec) const { return lattice.issame(vec->get_lattice()); }
     template<class R> bool samelattice(const LDVec<R> *)    const { return false; }
     template<class R> bool starlattice(const LQVec<R> *)    const { return false; }
     template<class R> bool starlattice(const LDVec<R> *vec) const { return lattice.isstar(vec->get_lattice()); }
     template<class R> bool samelattice(const LQVec<R> &vec) const { return lattice.issame(vec.get_lattice()); }
     template<class R> bool samelattice(const LDVec<R> &)    const { return false; }
     template<class R> bool starlattice(const LQVec<R> &)    const { return false; }
     template<class R> bool starlattice(const LDVec<R> &vec) const { return lattice.isstar(vec.get_lattice()); }
     // extract overloads preserving lattice information
     template<typename... A> LQVec<T> view(A... args) const;
     template<typename... A> LQVec<T> extract(A... args) const;
     bArray<T> get_hkl() const;
     bArray<double> get_xyz() const;
     LDVec<double> star() const;
   
     double dot(const size_t i, const size_t j) const;
     double norm(const size_t i) const { return sqrt(this->dot(i,i)); }
     LQVec<double> cross(const size_t i, const size_t j) const;
   
     LVEC_SCALAR_INPLACE_OP_DEF(LQVec,+=)
     LVEC_SCALAR_INPLACE_OP_DEF(LQVec,-=)
     LVEC_SCALAR_INPLACE_OP_DEF(LQVec,*=)
     LVEC_SCALAR_INPLACE_OP_DEF(LQVec,/=)
   
     LVEC_ARRAY_INPLACE_OP_DEF(LQVec,+=)
     LVEC_ARRAY_INPLACE_OP_DEF(LQVec,-=)
     LVEC_ARRAY_INPLACE_OP_DEF(LQVec,*=)
     LVEC_ARRAY_INPLACE_OP_DEF(LQVec,/=)
     // LQVec<T> operator -();
   
     template<class R>
     void binary_operation_check(const LQVec<R>& b) const{
       assert(this->samelattice(b));
     }
     template<class R>
     void binary_operation_check(const LDVec<R>& b) const{
       assert(this->starlattice(b));
     }
     template<class R, template<class> class A,
     typename=typename std::enable_if<!std::is_base_of<LatVec,A<R>>::value>::type>
     void binary_operation_check(const A<R>&) const {}
   
     template<class R>
     bool is(const LQVec<R>& that) const {
       return (this->samelattice(that) && this->bArray<T>::is(that));
     }
     LQVec<int> round() const {
       return LQVec<int>(this->lattice, this->bArray<T>::round());
     }
     LQVec<int> floor() const {
       return LQVec<int>(this->lattice, this->bArray<T>::floor());
     }
     LQVec<int> ceil() const {
       return LQVec<int>(this->lattice, this->bArray<T>::ceil());
     }
     //
     static LQVec<T> from_invA(const Reciprocal& lat, const bArray<T>& vec, const int=1){
       assert(vec.is_row_ordered() && vec.is_contiguous() && vec.size(vec.ndim()-1)==3);
       // initialise an empy array for the relative lattice unit components
       bArray<T> rlu(vec.shape(), vec.stride());
       auto vsh = vec.shape();
       vsh.back() = 0;
       // grab the transformation matric from the lattice
       std::vector<double> fromxyz = lat.get_inverse_xyz_transform();
       for (auto i: vec.subItr(vsh))
         brille::utils::multiply_matrix_vector(rlu.ptr(i), fromxyz.data(), vec.ptr(i));
       // finally make the vector(s) with lattice and components:
       return LQVec<T>(lat, rlu);
     }
   
     LQVec<T> decouple() {
       return LQVec<T>(lattice, this->bArray<T>::decouple());
     }
   protected:
     void check_array();
   };
   
   #undef LVEC_SCALAR_INPLACE_OP_DEF
   #undef LVEC_ARRAY_INPLACE_OP_DEF
   #undef LVEC_INIT_INT
   
   #ifndef DOXYGEN_SHOULD_SKIP_THIS
   // extend the lattice traits structs
   template<class T> struct LatticeTraits<LDVec<T>>{
     using type = Direct;
     using star = Reciprocal;
   };
   template<class T> struct LatticeTraits<LQVec<T>>{
     using type = Reciprocal;
     using star = Direct;
   };
   #endif
   
   template<class, class> struct LatVecTraits{
     using type = void; //< LDVec<R> or LQVec<R>
     using star = void; //< LQVec<R> or LDVec<R>
   };
   #ifndef DOXYGEN_SHOULD_SKIP_THIS
   template <class R, class S> struct LatVecTraits<LDVec<R>,S>{
     using type = LDVec<S>;
     using star = LQVec<S>;
   };
   template <class R, class S> struct LatVecTraits<LQVec<R>,S>{
     using type = LQVec<S>;
     using star = LDVec<S>;
   };
   template <class S> struct LatVecTraits<Direct,S>{
     using type = LDVec<S>;
     using star = LQVec<S>;
   };
   template <class S> struct LatVecTraits<Reciprocal,S>{
     using type = LQVec<S>;
     using star = LDVec<S>;
   };
   #endif
   
   template<class... T> using LatVecType = typename LatVecTraits<T...>::type;
   template<class... T> using LatVecStar = typename LatVecTraits<T...>::star;
   
   #include "array_latvec.tpp"
   #include "array_ldvec.tpp"
   #include "array_lqvec.tpp"
   
   }
   
   #endif