positioncollapse2d.cpp

Go to the documentation of this file.

/**
 * @file positioncollapse2d.cpp
 * @brief Implements the positioncollapse2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/positioncollapse2d.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 "positioncollapse2d.h"
 
#include "core/math/complexfield2d.h"
#include "core/quantum/gaussianpacket2d.h"
#include "core/quantum/positionmeasurement2d.h"
 
#include <cmath>
#include <random>
#include <vector>
 
namespace qpe {
namespace {
 
/**
 * @brief Performs collapse with seed.
 * @param field Field data to read or modify.
 * @param resolution The resolution value used by the operation.
 * @param seed The seed value used by the operation.
 */
PositionCollapse2DResult collapseWithSeed(ComplexField2D &field,
                                          double resolution,
                                          std::uint32_t seed)
{
    PositionCollapse2DResult result;
    result.resolution = resolution;
    result.normBefore = field.norm2();
 
    if (field.empty() || !(resolution > 0.0) || !std::isfinite(resolution) ||
        !(result.normBefore > 0.0) || !std::isfinite(result.normBefore)) {
        return result;
    }
 
    // For a Gaussian position POVM, the outcome distribution is the Born
    // density convolved with a Gaussian of standard deviation 'resolution'.
    // Sampling x~|psi|^2 and then adding Gaussian readout noise is an exact
    // Monte-Carlo draw from that convolution.
    const std::vector<PositionMeasurement2DPoint> latent =
        PositionMeasurement2D::sampleSeeded(field, 1u, seed);
    if (latent.empty()) {
        return result;
    }
 
    std::mt19937 engine(seed ^ 0x9e3779b9u);
    std::normal_distribution<double> noise(0.0, resolution);
    result.outcome.x = static_cast<float>(static_cast<double>(latent.front().x) + noise(engine));
    result.outcome.y = static_cast<float>(static_cast<double>(latent.front().y) + noise(engine));
 
    const double invFourSigma2 = 1.0 / (4.0 * resolution * resolution);
    for (std::size_t y = 0u; y < field.height(); ++y) {
        const double py = GaussianPacket2D::yCoordinate(field, y);
        const double dy = py - static_cast<double>(result.outcome.y);
        for (std::size_t x = 0u; x < field.width(); ++x) {
            const double px = GaussianPacket2D::xCoordinate(field, x);
            const double dx = px - static_cast<double>(result.outcome.x);
            const double exponent = -(dx * dx + dy * dy) * invFourSigma2;
            const float weight = static_cast<float>(std::exp(exponent));
            field(x, y) *= weight;
        }
    }
 
    if (!field.normalize()) {
        return result;
    }
 
    result.normAfter = field.norm2();
    result.valid = result.normAfter > 0.0 && std::isfinite(result.normAfter);
    return result;
}
 
} // namespace
 
PositionCollapse2DResult PositionCollapse2D::collapse(ComplexField2D &field,
                                                       double resolution)
{
    std::random_device rd;
    std::seed_seq seedSequence{rd(), rd(), rd(), rd(),
                               static_cast<unsigned int>(field.width()),
                               static_cast<unsigned int>(field.height())};
    std::uint32_t seed = 0u;
    seedSequence.generate(&seed, &seed + 1);
    return collapseWithSeed(field, resolution, seed);
}
 
PositionCollapse2DResult PositionCollapse2D::collapseSeeded(ComplexField2D &field,
                                                             double resolution,
                                                             std::uint32_t seed)
{
    return collapseWithSeed(field, resolution, seed);
}
 
} // namespace qpe