Unverified Commit 53afb8dd authored by Peter Doak's avatar Peter Doak Committed by GitHub
Browse files

Merge pull request #270 from PDoakORNL/function_real_conversion_fix

This looks good in testing and fixes a serious bug.
parents a65290d9 877b058b
Loading
Loading
Loading
Loading
+48 −12
Original line number Diff line number Diff line
// Copyright (C) 2018 ETH Zurich
// Copyright (C) 2018 UT-Battelle, LLC
// Copyright (C) 2022 ETH Zurich
// Copyright (C) 2022 UT-Battelle, LLC
// All rights reserved.
//
// See LICENSE for terms of usage.
// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications.
//
// Author: Giovanni Balduzzi (gbalduzz@gitp.phys.ethz.ch)
//         Peter Doak (doakpw@ornl.gov)
//
// This file provides utility methods to convert functions from real to complex and vice versa.

@@ -19,6 +20,7 @@
#include "dca/function/function.hpp"

namespace dca {
  enum class ImagCheck { IGNORE, WARN, FAIL };
namespace func {
namespace util {
// dca::func::util::
@@ -34,23 +36,57 @@ auto complex(const function<Scalartype, Dmn>& f) {
  return f_complex;
}

// Returns a real valued function that is equal to the real part of f.
// If check_imaginary = true, the method checks whether the imaginary part of f is zero and throws a
// std::logic_error if this is not the case.
template <typename Scalartype, typename Dmn>
auto real(const function<std::complex<Scalartype>, Dmn>& f, const bool check_imaginary = false) {
namespace detail {
 /** implementation of function complex to real conversion function
  *  This allows good code reuse with no loss or low loss of performance
  *  in the ignore and warn cases.
  */
template <typename Scalartype, typename Dmn, ImagCheck IC>
auto real(const function<std::complex<Scalartype>, Dmn>& f) {
  function<Scalartype, Dmn> f_real;

  auto checkForImaginary = [](const std::complex<Scalartype> s) -> bool {
    return (std::abs(s.imag()) > (1000 * std::numeric_limits<Scalartype>::epsilon()));
  };
  auto writeCheckFail = [](const int i, const std::complex<Scalartype> s) -> std::string {
    std::ostringstream fail_message;
    fail_message << "Element " << i << " of function conversion complex -> real "
                 << " has a non negligible imaginary component of " << std::abs(s.imag()) << '\n';
    return fail_message.str();
  };
  bool have_not_warned = true;
  for (int i = 0; i < f_real.size(); ++i) {
    if (check_imaginary && std::abs(f(i).imag()) > 1000 * std::numeric_limits<Scalartype>::epsilon()) {
      std::ostringstream err_msg;
      throw(std::logic_error(err_msg.str()));
    if constexpr (IC == ImagCheck::WARN) {
      if (have_not_warned) {
        if (checkForImaginary(f(i)))
          std::cerr << "WARNING: " << writeCheckFail(i, f(i));
      }
    }
    else if constexpr (IC == ImagCheck::FAIL) {
      if (checkForImaginary(f(i)))
        throw(std::out_of_range(writeCheckFail(i, f(i))));
    }
    f_real(i) = f(i).real();
  }

  return f_real;
}
}  // namespace detail

// Returns a real valued function that is equal to the real part of f.
// If IC ==
//   ImagCheck::FAIL    -> throw exception if imaginary part of f is zero
//   ImagCheck::WARN    -> write one warning per conversion to std::err continue without further
//   checks ImagCheck::IGNORE  -> do not check for imaginary elements in f at all
template <typename Scalartype, typename Dmn>
auto real(const function<std::complex<Scalartype>, Dmn>& f, ImagCheck ic = ImagCheck::IGNORE) {
  switch (ic) {
    case ImagCheck::IGNORE:
      return detail::real<Scalartype, Dmn, ImagCheck::IGNORE>(f);
    case ImagCheck::WARN:
      return detail::real<Scalartype, Dmn, ImagCheck::WARN>(f);
    case ImagCheck::FAIL:
      return detail::real<Scalartype, Dmn, ImagCheck::FAIL>(f);
  }
}

}  // namespace util
}  // namespace func
+6 −3
Original line number Diff line number Diff line
@@ -283,11 +283,14 @@ bool ADIOS2Writer<CT>::execute(const std::string& name,
  for (int i = 0; i < n; ++i) {
    sizes[i] = value[i].size();
    nTotal += sizes[i];
    if (verbose_) {
      std::cout << sizes[i];
      if (i < n - 1) {
        std::cout << ", ";
      }
    }
  }

  if (verbose_) {
    std::cout << ")\n";
  }
+3 −5
Original line number Diff line number Diff line
@@ -22,16 +22,14 @@ namespace io {
 *  dummy at 0 debug automatic initialization of int to 0
 *  causes a bug. That is not a good code smell.
 */
enum class IOType : int {
  JSON = 0,
  HDF5,
  ADIOS2
};
enum class IOType : int { JSON = 0, HDF5, ADIOS2, UNKNOWN };

IOType stringToIOType(const std::string& name);

std::string toString(const IOType type);

IOType extensionToIOType(const std::string& file_name);

}  // namespace io
}  // namespace dca

+3 −3
Original line number Diff line number Diff line
// Copyright (C) 2009-2016 ETH Zurich
// Copyright (C) 2007?-2016 Center for Nanophase Materials Sciences, ORNL
// Copyright (C) 2022 ETH Zurich
// Copyright (C) 2022 UT-Battelle, LLC
// All rights reserved.
//
// See LICENSE.txt for terms of usage.
@@ -68,7 +68,7 @@ public:
      func::function<Real, func::dmn_variadic<NuDmn, NuDmn, RDmn, LastDmn>>& f_output) {
    func::function<Complex, func::dmn_variadic<NuDmn, NuDmn, RDmn, LastDmn>> f_out_cmplx;
    execute(f_input, f_out_cmplx);
    f_output = std::move(func::util::real(f_out_cmplx, true));
    f_output = std::move(func::util::real(f_out_cmplx, ImagCheck::WARN));
  }

  // Default to old implementation when no band is present.
+3 −3
Original line number Diff line number Diff line
// Copyright (C) 2009-2016 ETH Zurich
// Copyright (C) 2007?-2016 Center for Nanophase Materials Sciences, ORNL
// Copyright (C) 2022 ETH Zurich
// Copyright (C) 2022 UT-Battelle, LLC
// All rights reserved.
//
// See LICENSE.txt for terms of usage.
@@ -72,7 +72,7 @@ public:
      func::function<Real, func::dmn_variadic<NuDmn, NuDmn, KDmn, LastDmn>>& f_output) {
    func::function<Complex, func::dmn_variadic<NuDmn, NuDmn, KDmn, LastDmn>> f_out_cmplx;
    execute(f_input, f_out_cmplx);
    f_output = std::move(func::util::real(f_out_cmplx, false));
    f_output = std::move(func::util::real(f_out_cmplx, ImagCheck::WARN));
  }

  // Default to old implementation when no band is present.
Loading