lambda_lanczos_util.hpp

Go to the documentation of this file.

#ifndef LAMBDA_LANCZOS_UTIL_H_
#define LAMBDA_LANCZOS_UTIL_H_
 
#include <vector>
#include <complex>
#include <limits>
#include <cmath>
#include <numeric>
#include <cassert>
#include <sstream>
#include <map>
#include <functional>
 
 
namespace lambda_lanczos { namespace util {
 
template <size_t I, typename container_type>
class TupleViewIterator {
private:
  typename container_type::const_iterator iter;
 
public:
  TupleViewIterator(const typename container_type::const_iterator& iter) : iter(iter) {}
 
  const typename std::tuple_element<I, typename container_type::value_type>::type& operator*() const {
    return std::get<I>(*iter);
  }
 
  bool operator==(const TupleViewIterator& obj) const {
    return this->iter == obj.iter;
  }
 
  bool operator!=(const TupleViewIterator& obj) const {
    return this->iter != obj.iter;
  }
 
  TupleViewIterator& operator++() {
    this->iter++;
    return *this;
  }
 
  TupleViewIterator operator++(int dummy) {
    auto before = *this;
    this->iter++;
    return before;
  }
};
 
template <typename map_type>
using MapValueIterator = TupleViewIterator<1, map_type>;
 
template <typename map_type>
class MapValueIterable {
private:
  const typename map_type::const_iterator itr_cbegin;
  const typename map_type::const_iterator itr_cend;
 
public:
  MapValueIterable(const map_type& map) : itr_cbegin(map.cbegin()), itr_cend(map.cend()) {}
 
  MapValueIterator<map_type> cbegin() const {
    return itr_cbegin;
  }
 
  MapValueIterator<map_type> cend() const {
    return itr_cend;
  }
};
 
 
template <typename T>
struct realTypeMap {
  typedef T type;
};
 
template <typename T>
struct realTypeMap<std::complex<T>> {
  typedef T type;
};
 
template <typename T>
using real_t = typename realTypeMap<T>::type;
 
 
template <typename T>
struct TypedConjugate {
  static T invoke(const T& val) {
    return val;
  }
};
 
template <typename T>
struct TypedConjugate<std::complex<T>> {
  static std::complex<T> invoke(const std::complex<T> val) {
    return std::conj(val);
  }
};
 
 
template <typename T>
inline T typed_conj(const T& val) {
  return TypedConjugate<T>::invoke(val);
}
 
 
template <typename T>
inline T inner_prod(const std::vector<T>& v1, const std::vector<T>& v2) {
  assert(v1.size() == v2.size());
  return std::inner_product(std::begin(v1), std::end(v1),
                            std::begin(v2), T(),
                            [](const T& v, const T& u) -> T { return v+u; },
                            [](const T& a, const T& b) -> T { return typed_conj(a)*b; });
  // T() means zero value of type T
}
 
template <typename T>
inline real_t<T> norm(const std::vector<T>& vec) {
  return std::sqrt(std::real(inner_prod(vec, vec)));
  // The norm of any complex vector <v|v> is real by definition.
}
 
template <typename T1, typename T2>
inline void scalar_mul(T1 a, std::vector<T2>& vec) {
  for(auto& elem : vec) {
    elem *= a;
  }
}
 
template <typename T>
inline void normalize(std::vector<T>& vec) {
  scalar_mul(T(1)/norm(vec), vec);
}
 
 
template <typename T>
inline real_t<T> l1_norm(const std::vector<T>& vec) {
  real_t<T> norm = real_t<T>(); // Zero initialization
 
  for(const T& element : vec) {
    norm += std::abs(element);
  }
 
  return norm;
}
 
 
template <typename ForwardIterator, typename T>
inline void schmidt_orth(std::vector<T>& uorth,
                         ForwardIterator first,
                         ForwardIterator last) {
  const auto n = uorth.size();
 
  for(auto iter = first; iter != last; ++iter) {
    const auto& uk = *iter;
    T innprod = util::inner_prod(uk, uorth);
 
    for(size_t i = 0; i < n; ++i) {
      uorth[i] -= innprod * uk[i];
    }
  }
}
 
 
template <typename T>
void initAsIdentity(std::vector<std::vector<T>>& a, size_t n) {
  a.resize(n);
  for(size_t i = 0; i < n; ++i) {
    a[i].resize(n);
    std::fill(a[i].begin(), a[i].end(), T());
    a[i][i] = 1.0;
  }
}
 
 
template <typename T>
inline void sort_eigenpairs(std::vector<real_t<T>>& eigenvalues,
                            std::vector<std::vector<T>>& eigenvectors,
                            bool sort_eigenvector,
                            const std::function<bool (real_t<T>, real_t<T>)> predicate = std::less<real_t<T>>()) {
  std::vector<std::pair<real_t<T>, size_t>> ev_index_pairs;
  ev_index_pairs.reserve(eigenvalues.size());
  for(size_t i = 0; i < eigenvalues.size(); ++i) {
    ev_index_pairs.emplace_back(eigenvalues[i], i);
  }
 
  std::sort(ev_index_pairs.begin(),
            ev_index_pairs.end(),
            [&predicate](const auto& x, const auto& y ) {
              return predicate(x.first, y.first);
            });
 
  std::vector<real_t<T>> eigenvalues_new;
  eigenvalues_new.reserve(eigenvalues.size());
  for(const auto& ev_index : ev_index_pairs) {
    size_t k = ev_index.second;
    eigenvalues_new.emplace_back(eigenvalues[k]);
  }
  eigenvalues = std::move(eigenvalues_new);
 
  if(sort_eigenvector) {
    std::vector<std::vector<T>> eigenvectors_new;
    eigenvectors_new.reserve(eigenvalues.size());
    for(const auto& ev_index : ev_index_pairs) {
      size_t k = ev_index.second;
      eigenvectors_new.push_back(std::move(eigenvectors[k]));
    }
    eigenvectors = std::move(eigenvectors_new);
  }
}
 
 
template <typename T>
inline constexpr int sig_decimal_digit() {
  return (int)(std::numeric_limits<T>::digits *
               log10(std::numeric_limits<T>::radix));
}
 
 
template <typename T>
inline constexpr T minimum_effective_decimal() {
  return pow(10, -sig_decimal_digit<T>());
}
 
 
template <typename T>
T sgn(T val) {
  if(val >= 0) {
    return (T)1;
  } else {
    return (T)(-1);
  }
}
 
 
template <typename T>
std::string vectorToString(const std::vector<T>& vec, std::string delimiter = " ") {
  std::stringstream ss;
 
  for(const auto& elem : vec) {
    ss << elem << delimiter;
  }
 
  /* Remove the last space */
  std::string result = ss.str();
  if(!result.empty()) {
    result.pop_back();
  }
 
  return result;
}
 
}} /* namespace lambda_lanczos::util */
 
#endif /* LAMBDA_LANCZOS_UTIL_H_ */