quantumbilliards2d.cpp

Go to the documentation of this file.

/**
 * @file quantumbilliards2d.cpp
 * @brief Implements the quantumbilliards2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumbilliards2d.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 "quantumbilliards2d.h"
 
#include "core/quantum/gaussianpacket2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <vector>
 
namespace qpe {
namespace {
 
constexpr double kPi = 3.141592653589793238462643383279502884;
 
/**
 * @brief Performs fail.
 * @param reason The reason value used by the operation.
 * @param message The message value used by the operation.
 */
void fail(std::string *reason, const char *message) noexcept
{
    if (reason) *reason = message;
}
 
/**
 * @brief Performs near.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 * @param scale The scale value used by the operation.
 */
bool near(double a, double b, double scale = 1.0) noexcept
{
    return std::fabs(a - b) <= 1.0e-4 * std::max(1.0, std::fabs(scale));
}
 
/**
 * @brief Performs same grid.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 */
bool sameGrid(const ComplexField2D &a, const ComplexField2D &b) noexcept
{
    return !a.empty() && !b.empty() &&
           a.width() == b.width() && a.height() == b.height() &&
           std::fabs(a.dx() - b.dx()) <= 1.0e-9f &&
           std::fabs(a.dy() - b.dy()) <= 1.0e-9f;
}
 
/**
 * @brief Performs inside accessible.
 * @param g The g value used by the operation.
 * @param x X coordinate/value.
 * @param y Y coordinate/value.
 */
bool insideAccessible(const QuantumBilliardGeometry2D &g, double x, double y) noexcept
{
    if (!g.valid) return false;
    const double hx = 0.5 * g.innerWidth;
    const double hy = 0.5 * g.innerHeight;
    if (x < g.centerX - hx || x > g.centerX + hx ||
        y < g.centerY - hy || y > g.centerY + hy) {
        return false;
    }
    if (g.hasObstacle) {
        const double dx = x - g.obstacleCenterX;
        const double dy = y - g.obstacleCenterY;
        if (dx * dx + dy * dy <= g.obstacleRadius * g.obstacleRadius) return false;
    }
    return true;
}
 
} // namespace
 
bool QuantumBilliards2D::identify(const PotentialScene2D &scene,
                                  QuantumBilliardGeometry2D &geometry,
                                  std::string *reason) noexcept
{
    geometry = QuantumBilliardGeometry2D();
    if (reason) reason->clear();
 
    std::vector<const PotentialObject2D *> vertical;
    std::vector<const PotentialObject2D *> horizontal;
    const PotentialObject2D *circle = nullptr;
 
    for (const PotentialObject2D &object : scene.objects()) {
        if (!object.enabled) continue;
        if (object.type == PotentialObjectType2D::Rectangle && object.value > 0.0f) {
            if (object.height > object.width) vertical.push_back(&object);
            else if (object.width > object.height) horizontal.push_back(&object);
            else {
                fail(reason, "Billiard walls must be clearly vertical or horizontal rectangles.");
                return false;
            }
        } else if (object.type == PotentialObjectType2D::Circle && object.value > 0.0f) {
            if (circle != nullptr) {
                fail(reason, "Billiard analysis supports at most one circular scatterer.");
                return false;
            }
            circle = &object;
        } else {
            fail(reason, "Billiard scene expects four positive Rectangle walls and optionally one positive Circle.");
            return false;
        }
    }
 
    if (vertical.size() != 2u || horizontal.size() != 2u) {
        fail(reason, "Billiard scene needs exactly two vertical and two horizontal walls.");
        return false;
    }
 
    std::sort(vertical.begin(), vertical.end(),
              [](const PotentialObject2D *a, const PotentialObject2D *b) {
                  return a->centerX < b->centerX;
              });
    std::sort(horizontal.begin(), horizontal.end(),
              [](const PotentialObject2D *a, const PotentialObject2D *b) {
                  return a->centerY < b->centerY;
              });
 
    const PotentialObject2D &left = *vertical[0];
    const PotentialObject2D &right = *vertical[1];
    const PotentialObject2D &bottom = *horizontal[0];
    const PotentialObject2D &top = *horizontal[1];
 
    const double leftInner = static_cast<double>(left.centerX) + 0.5 * left.width;
    const double rightInner = static_cast<double>(right.centerX) - 0.5 * right.width;
    const double bottomInner = static_cast<double>(bottom.centerY) + 0.5 * bottom.height;
    const double topInner = static_cast<double>(top.centerY) - 0.5 * top.height;
    const double width = rightInner - leftInner;
    const double height = topInner - bottomInner;
    if (!(width > 0.0) || !(height > 0.0)) {
        fail(reason, "Billiard walls overlap or leave no interior.");
        return false;
    }
 
    const double wallValue = left.value;
    if (!near(right.value, wallValue, wallValue) ||
        !near(bottom.value, wallValue, wallValue) ||
        !near(top.value, wallValue, wallValue)) {
        fail(reason, "All billiard walls must use the same potential height.");
        return false;
    }
 
    const double centerX = 0.5 * (leftInner + rightInner);
    const double centerY = 0.5 * (bottomInner + topInner);
    if (!near(left.centerY, centerY, height) || !near(right.centerY, centerY, height) ||
        !near(bottom.centerX, centerX, width) || !near(top.centerX, centerX, width)) {
        fail(reason, "Billiard walls are not centred on one rectangular interior.");
        return false;
    }
 
    geometry.valid = true;
    geometry.centerX = centerX;
    geometry.centerY = centerY;
    geometry.innerWidth = width;
    geometry.innerHeight = height;
    geometry.wallPotential = wallValue;
 
    if (circle != nullptr) {
        const double radius = static_cast<double>(circle->radius);
        if (!(radius > 0.0)) {
            fail(reason, "Circular billiard scatterer must have positive radius.");
            geometry = QuantumBilliardGeometry2D();
            return false;
        }
        const double hx = 0.5 * width;
        const double hy = 0.5 * height;
        if (circle->centerX - radius <= centerX - hx ||
            circle->centerX + radius >= centerX + hx ||
            circle->centerY - radius <= centerY - hy ||
            circle->centerY + radius >= centerY + hy) {
            fail(reason, "Circular scatterer must lie completely inside the billiard.");
            geometry = QuantumBilliardGeometry2D();
            return false;
        }
        geometry.hasObstacle = true;
        geometry.obstacleCenterX = circle->centerX;
        geometry.obstacleCenterY = circle->centerY;
        geometry.obstacleRadius = radius;
        geometry.obstaclePotential = circle->value;
    }
 
    return true;
}
 
double QuantumBilliards2D::accessibleArea(const QuantumBilliardGeometry2D &geometry) noexcept
{
    if (!geometry.valid || !(geometry.innerWidth > 0.0) || !(geometry.innerHeight > 0.0)) return 0.0;
    double area = geometry.innerWidth * geometry.innerHeight;
    if (geometry.hasObstacle && geometry.obstacleRadius > 0.0) {
        area -= kPi * geometry.obstacleRadius * geometry.obstacleRadius;
    }
    return std::max(0.0, area);
}
 
double QuantumBilliards2D::fidelity(const ComplexField2D &reference,
                                    const ComplexField2D &field) noexcept
{
    if (!sameGrid(reference, field)) return 0.0;
    std::complex<double> overlap(0.0, 0.0);
    double na = 0.0;
    double nb = 0.0;
    for (std::size_t i = 0u; i < field.size(); ++i) {
        const std::complex<double> a(reference.data()[i].real(), reference.data()[i].imag());
        const std::complex<double> b(field.data()[i].real(), field.data()[i].imag());
        overlap += std::conj(a) * b;
        na += std::norm(a);
        nb += std::norm(b);
    }
    if (!(na > 0.0) || !(nb > 0.0)) return 0.0;
    return std::norm(overlap) / (na * nb);
}
 
QuantumBilliardMetrics2D QuantumBilliards2D::measure(const ComplexField2D &field,
                                                       const QuantumBilliardGeometry2D &geometry,
                                                       const ComplexField2D *reference) noexcept
{
    QuantumBilliardMetrics2D out;
    if (field.empty() || !geometry.valid) return out;
 
    const double cellArea = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
    double interiorNorm = 0.0;
    double densitySquareIntegral = 0.0;
    double shellProbability = 0.0;
    const double shellWidth = geometry.hasObstacle
        ? std::max(0.45, 2.5 * std::max(static_cast<double>(field.dx()), static_cast<double>(field.dy())))
        : 0.0;
    const double shellOuter = geometry.obstacleRadius + shellWidth;
    const double shellOuter2 = shellOuter * shellOuter;
    const double shellInner2 = geometry.obstacleRadius * geometry.obstacleRadius;
 
    for (std::size_t y = 0u; y < field.height(); ++y) {
        const double py = GaussianPacket2D::yCoordinate(field, y);
        for (std::size_t x = 0u; x < field.width(); ++x) {
            const double px = GaussianPacket2D::xCoordinate(field, x);
            const ComplexField2D::Value v = field(x, y);
            const double density = static_cast<double>(v.real()) * v.real() +
                                   static_cast<double>(v.imag()) * v.imag();
            const double probability = density * cellArea;
            if (insideAccessible(geometry, px, py)) {
                interiorNorm += probability;
                densitySquareIntegral += density * density * cellArea;
            }
            if (geometry.hasObstacle) {
                const double dx = px - geometry.obstacleCenterX;
                const double dy = py - geometry.obstacleCenterY;
                const double r2 = dx * dx + dy * dy;
                if (r2 > shellInner2 && r2 <= shellOuter2) shellProbability += probability;
            }
        }
    }
 
    out.accessibleProbability = interiorNorm;
    if (interiorNorm > 1.0e-20 && densitySquareIntegral > 1.0e-30) {
        // Renormalize density to the accessible probability before evaluating IPR.
        out.effectiveArea = interiorNorm * interiorNorm / densitySquareIntegral;
        const double area = accessibleArea(geometry);
        if (area > 0.0) out.occupiedFraction = std::max(0.0, std::min(1.0, out.effectiveArea / area));
    }
    out.obstacleShellProbability = shellProbability;
    if (reference != nullptr) out.fidelity = fidelity(*reference, field);
    return out;
}
 
QuantumBilliardHistory2D::QuantumBilliardHistory2D(std::size_t maximumSamples)
    : maximumSamples_(std::max<std::size_t>(1u, maximumSamples))
{
}
 
void QuantumBilliardHistory2D::clear() noexcept
{
    reference_.clear();
    samples_.clear();
    meanProbabilityDensity_.clear();
    accumulatedFrames_ = 0u;
}
 
bool QuantumBilliardHistory2D::captureReference(const ComplexField2D &field,
                                                 const QuantumBilliardGeometry2D &geometry)
{
    clear();
    if (field.empty() || !geometry.valid) return false;
    reference_ = field;
    meanProbabilityDensity_.assign(field.size(), 0.0f);
    return true;
}
 
void QuantumBilliardHistory2D::setMaximumSamples(std::size_t maximumSamples)
{
    maximumSamples_ = std::max<std::size_t>(1u, maximumSamples);
    trimToLimit();
}
 
bool QuantumBilliardHistory2D::append(double time,
                                      const ComplexField2D &field,
                                      const QuantumBilliardGeometry2D &geometry)
{
    if (!hasReference() || !sameGrid(reference_, field) || !geometry.valid) return false;
 
    const QuantumBilliardMetrics2D metrics = QuantumBilliards2D::measure(field, geometry, &reference_);
    QuantumBilliardSample2D sample;
    sample.time = time;
    sample.effectiveArea = metrics.effectiveArea;
    sample.occupiedFraction = metrics.occupiedFraction;
    sample.obstacleShellProbability = metrics.obstacleShellProbability;
    sample.fidelity = metrics.fidelity;
    samples_.push_back(sample);
 
    ++accumulatedFrames_;
    const double n = static_cast<double>(accumulatedFrames_);
    for (std::size_t i = 0u; i < field.size(); ++i) {
        const ComplexField2D::Value v = field.data()[i];
        const double density = static_cast<double>(v.real()) * v.real() +
                               static_cast<double>(v.imag()) * v.imag();
        const double previous = meanProbabilityDensity_[i];
        meanProbabilityDensity_[i] = static_cast<float>(previous + (density - previous) / n);
    }
 
    trimToLimit();
    return true;
}
 
void QuantumBilliardHistory2D::trimToLimit()
{
    if (samples_.size() <= maximumSamples_) return;
    const std::size_t removeCount = samples_.size() - maximumSamples_;
    samples_.erase(samples_.begin(), samples_.begin() + static_cast<std::ptrdiff_t>(removeCount));
}
 
} // namespace qpe