positionmeasurement2d.cpp

Go to the documentation of this file.

/**
 * @file positionmeasurement2d.cpp
 * @brief Implements the positionmeasurement2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/positionmeasurement2d.cpp. This file is part of Quantum Phenomena Explorer (QPE),
 * maintained by Officina Turini. See LICENSE_IT.txt / LICENSE.txt for the
 * QPE freeware source-available terms and THIRD_PARTY_NOTICES.txt for
 * third-party components.
 */
 
#include "positionmeasurement2d.h"
 
#include "core/math/complexfield2d.h"
 
#include <cmath>
#include <random>
 
namespace qpe {
namespace {
 
std::vector<PositionMeasurement2DPoint> sampleWithEngine(
    const ComplexField2D &field,
    std::size_t shots,
    std::mt19937 &engine)
{
    std::vector<PositionMeasurement2DPoint> result;
    if (field.empty() || shots == 0u) {
        return result;
    }
 
    std::vector<double> weights(field.size(), 0.0);
    double total = 0.0;
    const ComplexField2D::Value *values = field.data();
    for (std::size_t i = 0; i < field.size(); ++i) {
        const double re = static_cast<double>(values[i].real());
        const double im = static_cast<double>(values[i].imag());
        const double p = re * re + im * im;
        if (std::isfinite(p) && p > 0.0) {
            weights[i] = p;
            total += p;
        }
    }
 
    if (!(total > 0.0) || !std::isfinite(total)) {
        return result;
    }
 
    std::discrete_distribution<std::size_t> distribution(weights.begin(), weights.end());
    result.reserve(shots);
 
    const double cx = 0.5 * static_cast<double>(field.width() - 1u);
    const double cy = 0.5 * static_cast<double>(field.height() - 1u);
    const double dx = static_cast<double>(field.dx());
    const double dy = static_cast<double>(field.dy());
 
    for (std::size_t shot = 0; shot < shots; ++shot) {
        const std::size_t index = distribution(engine);
        const std::size_t xIndex = index % field.width();
        const std::size_t yIndex = index / field.width();
 
        PositionMeasurement2DPoint point;
        point.x = static_cast<float>((static_cast<double>(xIndex) - cx) * dx);
        point.y = static_cast<float>((static_cast<double>(yIndex) - cy) * dy);
        result.push_back(point);
    }
 
    return result;
}
 
} // namespace
 
std::vector<PositionMeasurement2DPoint> PositionMeasurement2D::sample(
    const ComplexField2D &field,
    std::size_t shots)
{
    std::random_device rd;
    std::seed_seq seed{
        rd(), rd(), rd(), rd(),
        static_cast<unsigned int>(field.width()),
        static_cast<unsigned int>(field.height()),
        static_cast<unsigned int>(shots)
    };
    std::mt19937 engine(seed);
    return sampleWithEngine(field, shots, engine);
}
 
std::vector<PositionMeasurement2DPoint> PositionMeasurement2D::sampleSeeded(
    const ComplexField2D &field,
    std::size_t shots,
    std::uint32_t seed)
{
    std::mt19937 engine(seed);
    return sampleWithEngine(field, shots, engine);
}
 
} // namespace qpe