File radialprojectionresult.h
File List > code_source > templat > include > TempLat > lattice > measuringtools > projectionhelpers > radialprojectionresult.h
Go to the documentation of this file
#ifndef TEMPLAT_LATTICE_MEASUREMENTS_PROJECTIONHELPERS_RADIALPROJECTIONRESULT_H
#define TEMPLAT_LATTICE_MEASUREMENTS_PROJECTIONHELPERS_RADIALPROJECTIONRESULT_H
/* This file is part of TempLat, available at https://cosmolattice.github.io/templat .
Copyright 2021-2026 The TempLat authors, see AUTHORS.md.
Released under the MIT license, see LICENSE.md. */
// File info: Main contributor(s): Wessel Valkenburg, Year: 2019
#include "TempLat/util/exception.h"
#include "TempLat/lattice/algebra/helpers/getfloattype.h"
#include "TempLat/lattice/measuringtools/projectionhelpers/radialprojectionsinglebinandvalue.h"
#include "TempLat/lattice/measuringtools/projectionhelpers/radialprojectionsinglequantity.h"
#include "TempLat/lattice/measuringtools/projectionhelpers/radialprojectionrebinner.h"
#include "TempLat/parallel/device.h"
namespace TempLat
{
MakeException(RadialProjectionResultSizeException);
MakeException(RadialProjectionResultFinalizationException);
template <typename T> class RadialProjectionResult : public std::vector<RadialProjectionSingleBinAndValue<T>>
{
public:
using floatType = typename GetFloatType<T>::type;
template <typename U> friend class RadialProjectionResult;
// Put public methods here. These should change very little over time.
RadialProjectionResult(size_t nBins, T deltakBins, bool pUseBinCentralValues = false, bool pIsInFourier = false)
: std::vector<RadialProjectionSingleBinAndValue<T>>(), finalizedOnce(false), mNBins(nBins), mDeltakBins(deltakBins), mValues(mNBins),
mBinBounds(mNBins), centralBinBounds(mNBins), mUseBinCentralValues(pUseBinCentralValues), mIsInFourier(pIsInFourier)
{
mMultiplicitiesDevice = DeviceView("RadialProjectionResult::mMultiplicitiesDevice", mNBins);
mMultiplicities = device::memory::createMirrorView(mMultiplicitiesDevice);
}
template <typename U>
RadialProjectionResult(const RadialProjectionResult<U> &other)
: std::vector<RadialProjectionSingleBinAndValue<T>>(), finalizedOnce(other.finalizedOnce), mNBins(other.mNBins),
mDeltakBins(static_cast<T>(other.mDeltakBins)), mValues(other.mNBins), mBinBounds(other.mNBins),
mUseBinCentralValues(other.mUseBinCentralValues), mIsInFourier(other.mIsInFourier)
{
this->reserve(other.size());
for (const auto &bv : other)
this->emplace_back(bv);
centralBinBounds.reserve(other.centralBinBounds.size());
for (const auto &b : other.centralBinBounds)
centralBinBounds.push_back(static_cast<T>(b));
mMultiplicitiesDevice = DeviceView("RadialProjectionResult::mMultiplicitiesDevice", mNBins);
mMultiplicities = device::memory::createMirrorView(mMultiplicitiesDevice);
}
RadialProjectionResult &rebin(device::Idx nBins, T customRange = -1)
{
std::vector<RadialProjectionSingleBinAndValue<T>>::operator=(
RadialProjectionRebinner<T>::rebin(*this, nBins, centralBinBounds, customRange));
return *this;
}
auto getNBins() { return mNBins; }
template <typename LL> RadialProjectionResult &rescale(LL rescaler)
{
// for (auto&& it : *this) {
for (size_t i = 0; i < this->size(); ++i) {
if (mUseBinCentralValues)
(*this)[i].getValue() *= rescaler(this->centralBinBounds[i]);
else
(*this)[i].getValue() *= rescaler((*this)[i].getBin().average);
}
return *this;
}
auto integrate(bool useCentralBin)
{
auto total = 0.;
auto kIR = (*this).getCentralBinBounds()[0];
if (useCentralBin) {
for (size_t i = 0; i < this->size(); ++i)
total += (*this)[i].getValue().average / centralBinBounds[i];
} else {
for (auto &&it : *this)
total += it.getValue().average * kIR / it.getBin().average;
}
return total * mDeltakBins;
}
RadialProjectionResult &rescaleBins(T scale)
{
for (auto &&it : *this) {
it.getBin() *= scale;
}
for (T &b : centralBinBounds) {
b *= scale;
}
return *this;
}
RadialProjectionResult &sumInsteadOfAverage()
{
floatType intMultiplicity = 0;
/*for (auto&& it : *this) {
intMultiplicity = 2 * it.getBin().multiplicity; // *2 is to get the full number of modes, against original
design it.getValue() *= intMultiplicity;
//it.getValue().sumInsteadOfAverage();
}*/
for (size_t i = 0; i < this->size(); ++i) {
intMultiplicity =
(*this)[i].getBin().multiplicity *
(mIsInFourier ? 2 : 1); // Multiply by 2 if in Fourier space (only half of the last coordinate is iterated
// over in this case because of reflection symmetry for real data).
(*this)[i].getValue() *= intMultiplicity;
}
return *this;
}
std::string toString(int verbosity = 0) const
{
if ((device::Idx)this->size() < 1) return "";
std::stringstream sstream;
sstream << this->front().getHeader(verbosity) << "\n";
for (auto &&it : *this) {
sstream << it.toString(verbosity) << "\n";
}
return sstream.str();
}
friend std::ostream &operator<<(std::ostream &ostream, const RadialProjectionResult &rpr)
{
ostream << rpr.toString() << "\n";
return ostream;
}
std::vector<T> &getCentralBinBounds() { return centralBinBounds; }
template <typename S> friend class RadialProjector;
template <typename S>
friend RadialProjectionResult<S> operator+(const RadialProjectionResult<S> &a, const RadialProjectionResult<S> &b);
DEVICE_INLINE_FUNCTION
void add_device(device::Idx i, const T &value, const T &position, const T &weight = (T)1) const
{
mValues.add_device(i, value, weight);
mBinBounds.add_device(i, position, weight);
device::atomic_add(&mMultiplicitiesDevice(i), weight);
}
RadialProjectionResult &finalize(MPICommReference comm)
{
if (finalizedOnce) throw RadialProjectionResultFinalizationException("Can only finalize once per instance.");
finalizedOnce = true;
pull();
comm.Allreduce(mMultiplicities, MPI_SUM);
mValues.finalize(comm);
mBinBounds.finalize(comm);
std::vector<T> updatedCentralBinBounds;
auto &mFullResult = *this;
mFullResult.clear();
for (size_t i = 0; i < mNBins; ++i) {
if (mMultiplicities[i] == 0) continue;
RadialProjectionSingleBinAndValue<T> next(mBinBounds.getFinal(i, mMultiplicities[i]),
mValues.getFinal(i, mMultiplicities[i]));
this->push_back(next);
updatedCentralBinBounds.push_back(centralBinBounds[i]);
}
centralBinBounds = updatedCentralBinBounds; //Needed to ensure the number of bins is consistent within the different members of this class.
return *this;
}
private:
/* Put all member variables and private methods here. These may change arbitrarily. */
bool finalizedOnce;
size_t mNBins;
T mDeltakBins;
RadialProjectionSingleQuantity<T> mValues, mBinBounds;
// Naive central values of the bin. Sized to mNBins at construction and zero-filled, so finalize() can
// always index it; RadialProjector overwrites it with the real values via setCentralBinBounds().
std::vector<T> centralBinBounds;
using DeviceView = device::memory::NDView<floatType, 1>;
using HostMirror = typename DeviceView::host_mirror_type;
HostMirror mMultiplicities;
DeviceView mMultiplicitiesDevice;
bool mUseBinCentralValues;
bool mIsInFourier;
void pull()
{
mValues.pull();
mBinBounds.pull();
device::memory::copyDeviceToHost(mMultiplicitiesDevice, mMultiplicities.data());
}
void push()
{
mValues.push();
mBinBounds.push();
device::memory::copyHostToDevice(mMultiplicities.data(), mMultiplicitiesDevice);
}
};
template <typename T>
RadialProjectionResult<T> operator+(const RadialProjectionResult<T> &a, const RadialProjectionResult<T> &b)
{
RadialProjectionResult<T> res(a);
for (size_t i = 0; i < a.size(); ++i) {
res[i].getValue().average += b[i].getValue().average;
}
return res;
}
template <typename T, class R> RadialProjectionResult<T> operator*(R &&func, const RadialProjectionResult<T> &obj)
{
RadialProjectionResult<T> res(obj);
return res.rescale(func);
}
template <typename T> RadialProjectionResult<T> operator*(double scale, const RadialProjectionResult<T> &obj)
{
auto func = [&](auto x) { return scale; };
return func * obj;
}
template <typename T> RadialProjectionResult<T> operator*(float scale, const RadialProjectionResult<T> &obj)
{
auto func = [&](auto x) { return scale; };
return func * obj;
}
template <typename T> RadialProjectionResult<T> operator*(int scale, const RadialProjectionResult<T> &obj)
{
auto func = [&](auto x) { return scale; };
return func * obj;
}
} // namespace TempLat
#endif