quantumanderson2d.cpp

Go to the documentation of this file.

/**
 * @file quantumanderson2d.cpp
 * @brief Implements the quantumanderson2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumanderson2d.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 "core/quantum/quantumanderson2d.h"
 
#include "core/quantum/gaussianpacket2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <limits>
 
namespace qpe {
namespace {
 
const double kPi = 3.14159265358979323846;
 
/**
 * @brief Sets reason.
 * @param reason The reason value used by the operation.
 * @param message The message value used by the operation.
 */
void setReason(std::string *reason, const char *message)
{
    if (reason) *reason = message;
}
 
/**
 * @brief Performs finite.
 * @param v The v value used by the operation.
 */
bool finite(double v) noexcept
{
    return std::isfinite(v);
}
 
struct Row
{
    double y = 0.0;
    std::vector<const PotentialObject2D *> objects;
};
 
/**
 * @brief Performs std dev.
 * @param values The values value used by the operation.
 */
double stdDev(const std::vector<double> &values)
{
    if (values.empty()) return 0.0;
    double mean = 0.0;
    for (double v : values) mean += v;
    mean /= static_cast<double>(values.size());
    double var = 0.0;
    for (double v : values) {
        const double d = v - mean;
        var += d * d;
    }
    return std::sqrt(var / static_cast<double>(values.size()));
}
 
} // namespace
 
bool QuantumAnderson2D::validateConfig(const QuantumAndersonInitialState2DConfig &config,
                                       std::string *reason)
{
    if (!finite(config.centerX) || !finite(config.cleanLaneY) ||
        !finite(config.disorderLaneY) || !finite(config.sigma) ||
        !finite(config.relativePhaseDegrees)) {
        setReason(reason, "QuantumAnderson2D: initial-state values must be finite");
        return false;
    }
    if (!(config.sigma > 0.0)) {
        setReason(reason, "QuantumAnderson2D: sigma must be positive");
        return false;
    }
    if (std::abs(config.cleanLaneY - config.disorderLaneY) < 4.0 * config.sigma) {
        setReason(reason, "QuantumAnderson2D: the two lanes are too close for a clean comparison");
        return false;
    }
    return true;
}
 
bool QuantumAnderson2D::initializeDualLane(ComplexField2D &field,
                                           const QuantumAndersonInitialState2DConfig &config,
                                           std::string *reason)
{
    if (field.empty()) {
        setReason(reason, "QuantumAnderson2D: field must not be empty");
        return false;
    }
    if (!validateConfig(config, reason)) return false;
 
    const double inv4Sigma2 = 1.0 / (4.0 * config.sigma * config.sigma);
    const double phase = config.relativePhaseDegrees * kPi / 180.0;
    const std::complex<double> disorderPhase(std::cos(phase), std::sin(phase));
 
    for (std::size_t y = 0; y < field.height(); ++y) {
        const double py = GaussianPacket2D::yCoordinate(field, y);
        for (std::size_t x = 0; x < field.width(); ++x) {
            const double px = GaussianPacket2D::xCoordinate(field, x);
            const double dx = px - config.centerX;
            const double dyClean = py - config.cleanLaneY;
            const double dyDisorder = py - config.disorderLaneY;
            const double a = std::exp(-(dx * dx + dyClean * dyClean) * inv4Sigma2);
            const double b = std::exp(-(dx * dx + dyDisorder * dyDisorder) * inv4Sigma2);
            const std::complex<double> value = std::complex<double>(a, 0.0) + disorderPhase * b;
            field(x, y) = ComplexField2D::Value(static_cast<float>(value.real()),
                                                static_cast<float>(value.imag()));
        }
    }
 
    if (!field.normalize()) {
        setReason(reason, "QuantumAnderson2D: normalization failed");
        return false;
    }
    return true;
}
 
bool QuantumAnderson2D::identifyLattice(const PotentialScene2D &scene,
                                        QuantumAndersonLattice2D &lattice,
                                        std::string *reason)
{
    lattice = QuantumAndersonLattice2D();
    std::vector<const PotentialObject2D *> gaussians;
    for (const PotentialObject2D &object : scene.objects()) {
        if (!object.enabled) continue;
        if (object.type != PotentialObjectType2D::Gaussian) {
            setReason(reason, "QuantumAnderson2D: guided comparison expects Gaussian lattice sites only");
            return false;
        }
        gaussians.push_back(&object);
    }
    if (gaussians.size() < 14u || (gaussians.size() % 2u) != 0u) {
        setReason(reason, "QuantumAnderson2D: two equal Gaussian rows are required");
        return false;
    }
 
    std::sort(gaussians.begin(), gaussians.end(), [](const PotentialObject2D *a,
                                                      const PotentialObject2D *b) {
        if (a->centerY != b->centerY) return a->centerY < b->centerY;
        return a->centerX < b->centerX;
    });
 
    const double yLow = gaussians.front()->centerY;
    const double yHigh = gaussians.back()->centerY;
    if (std::abs(yHigh - yLow) < 1.0e-6) {
        setReason(reason, "QuantumAnderson2D: two distinct lattice rows were not found");
        return false;
    }
 
    Row low; low.y = yLow;
    Row high; high.y = yHigh;
    for (const PotentialObject2D *object : gaussians) {
        if (std::abs(object->centerY - yLow) < 1.0e-5) low.objects.push_back(object);
        else if (std::abs(object->centerY - yHigh) < 1.0e-5) high.objects.push_back(object);
        else {
            setReason(reason, "QuantumAnderson2D: more than two Gaussian rows were found");
            return false;
        }
    }
    if (low.objects.size() != high.objects.size() || low.objects.size() < 7u) {
        setReason(reason, "QuantumAnderson2D: lattice rows must contain the same number of sites");
        return false;
    }
 
    auto sortX = [](Row &row) {
        std::sort(row.objects.begin(), row.objects.end(), [](const PotentialObject2D *a,
                                                             const PotentialObject2D *b) {
            return a->centerX < b->centerX;
        });
    };
    sortX(low); sortX(high);
 
    const double spacing = low.objects[1]->centerX - low.objects[0]->centerX;
    if (!(spacing > 0.0)) {
        setReason(reason, "QuantumAnderson2D: invalid lattice spacing");
        return false;
    }
 
    std::vector<double> lowDepths;
    std::vector<double> highDepths;
    for (std::size_t i = 0; i < low.objects.size(); ++i) {
        const PotentialObject2D &a = *low.objects[i];
        const PotentialObject2D &b = *high.objects[i];
        if (std::abs(a.centerX - b.centerX) > 1.0e-4 ||
            std::abs((a.centerX - low.objects[0]->centerX) - static_cast<double>(i) * spacing) > 2.0e-3 ||
            std::abs((b.centerX - high.objects[0]->centerX) - static_cast<double>(i) * spacing) > 2.0e-3) {
            setReason(reason, "QuantumAnderson2D: the two rows do not share one regular lattice");
            return false;
        }
        if (!(a.sigmaX > 0.0f && a.sigmaY > 0.0f && b.sigmaX > 0.0f && b.sigmaY > 0.0f)) {
            setReason(reason, "QuantumAnderson2D: Gaussian site widths must be positive");
            return false;
        }
        lowDepths.push_back(static_cast<double>(a.value));
        highDepths.push_back(static_cast<double>(b.value));
    }
 
    const double lowStd = stdDev(lowDepths);
    const double highStd = stdDev(highDepths);
    const Row *clean = lowStd <= highStd ? &low : &high;
    const Row *disorder = lowStd <= highStd ? &high : &low;
    const std::vector<double> &cleanDepths = lowStd <= highStd ? lowDepths : highDepths;
    const std::vector<double> &disorderDepths = lowStd <= highStd ? highDepths : lowDepths;
 
    double cleanMean = 0.0, disorderMean = 0.0;
    for (double v : cleanDepths) cleanMean += v;
    for (double v : disorderDepths) disorderMean += v;
    cleanMean /= static_cast<double>(cleanDepths.size());
    disorderMean /= static_cast<double>(disorderDepths.size());
 
    if (stdDev(cleanDepths) > 1.0e-3 || stdDev(disorderDepths) < 1.0e-2) {
        setReason(reason, "QuantumAnderson2D: one uniform row and one disordered row are required");
        return false;
    }
 
    lattice.valid = true;
    lattice.cleanLaneY = clean->y;
    lattice.disorderLaneY = disorder->y;
    lattice.spacing = spacing;
    lattice.siteSigmaX = clean->objects.front()->sigmaX;
    lattice.siteSigmaY = clean->objects.front()->sigmaY;
    lattice.cleanDepth = cleanMean;
    lattice.disorderDepthMean = disorderMean;
    lattice.disorderDepthStdDev = stdDev(disorderDepths);
    lattice.cleanDepths = cleanDepths;
    lattice.disorderDepths = disorderDepths;
    lattice.siteX.reserve(clean->objects.size());
    for (const PotentialObject2D *object : clean->objects) lattice.siteX.push_back(object->centerX);
    return true;
}
 
QuantumAndersonLaneSample2D QuantumAnderson2D::measureLane(
    const ComplexField2D &field,
    const QuantumAndersonLattice2D &lattice,
    bool cleanLane)
{
    QuantumAndersonLaneSample2D out;
    if (field.empty() || !lattice.valid || lattice.siteX.empty()) return out;
 
    const double splitY = 0.5 * (lattice.cleanLaneY + lattice.disorderLaneY);
    const bool cleanAbove = lattice.cleanLaneY > lattice.disorderLaneY;
    const double cellArea = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
    out.marginalX.assign(field.width(), 0.0f);
    out.profileDx = static_cast<double>(field.dx());
    out.profileXMin = GaussianPacket2D::xCoordinate(field, 0u);
 
    for (std::size_t y = 0; y < field.height(); ++y) {
        const double py = GaussianPacket2D::yCoordinate(field, y);
        const bool above = py >= splitY;
        if ((cleanLane && above != cleanAbove) || (!cleanLane && above == cleanAbove)) continue;
        for (std::size_t x = 0; x < field.width(); ++x) {
            const ComplexField2D::Value &v = field(x, y);
            const double p = (static_cast<double>(v.real()) * v.real() +
                              static_cast<double>(v.imag()) * v.imag()) * cellArea;
            out.probability += p;
            out.marginalX[x] += static_cast<float>(p / static_cast<double>(field.dx()));
        }
    }
    if (!(out.probability > 0.0)) return out;
 
    double weightedX = 0.0;
    for (std::size_t x = 0; x < field.width(); ++x) {
        const double px = GaussianPacket2D::xCoordinate(field, x);
        weightedX += static_cast<double>(out.marginalX[x]) * static_cast<double>(field.dx()) * px;
        out.marginalX[x] = static_cast<float>(static_cast<double>(out.marginalX[x]) / out.probability);
    }
    out.meanX = weightedX / out.probability;
 
    double variance = 0.0;
    for (std::size_t x = 0; x < field.width(); ++x) {
        const double px = GaussianPacket2D::xCoordinate(field, x);
        const double p = static_cast<double>(out.marginalX[x]) * static_cast<double>(field.dx());
        const double d = px - out.meanX;
        variance += p * d * d;
    }
    out.sigmaX = std::sqrt(std::max(0.0, variance));
 
    std::vector<double> siteProbabilities(lattice.siteX.size(), 0.0);
    const double halfCell = 0.5 * lattice.spacing;
    for (std::size_t x = 0; x < field.width(); ++x) {
        const double px = GaussianPacket2D::xCoordinate(field, x);
        std::size_t nearest = 0u;
        double best = std::numeric_limits<double>::max();
        for (std::size_t i = 0; i < lattice.siteX.size(); ++i) {
            const double d = std::abs(px - lattice.siteX[i]);
            if (d < best) { best = d; nearest = i; }
        }
        if (best <= halfCell) {
            siteProbabilities[nearest] += static_cast<double>(out.marginalX[x]) * static_cast<double>(field.dx());
        }
    }
    double sumSquares = 0.0;
    for (double p : siteProbabilities) sumSquares += p * p;
    if (sumSquares > 0.0) out.siteParticipation = 1.0 / sumSquares;
    out.centralSiteProbability = siteProbabilities[siteProbabilities.size() / 2u];
 
    if (!cleanLane) {
        double sx = 0.0, sy = 0.0, sxx = 0.0, sxy = 0.0;
        int count = 0;
        const double centreX = lattice.siteX[lattice.siteX.size() / 2u];
        for (std::size_t i = 0; i < siteProbabilities.size(); ++i) {
            const double r = std::abs(lattice.siteX[i] - centreX);
            const double p = siteProbabilities[i];
            if (r < 0.75 * lattice.spacing || p < 1.0e-5) continue;
            const double yy = std::log(p);
            sx += r; sy += yy; sxx += r * r; sxy += r * yy; ++count;
        }
        const double denom = static_cast<double>(count) * sxx - sx * sx;
        if (count >= 4 && std::abs(denom) > 1.0e-12) {
            const double slope = (static_cast<double>(count) * sxy - sx * sy) / denom;
            if (slope < -1.0e-6) {
                out.localizationLength = -2.0 / slope; // P ~ exp(-2 |x| / xi)
                out.localizationLengthValid = finite(out.localizationLength);
            }
        }
    }
    return out;
}
 
double QuantumAnderson2D::freeGaussianSigma(double time, double sigma0,
                                             double mass, double hbar) noexcept
{
    if (!(sigma0 > 0.0) || !(mass > 0.0) || !(hbar > 0.0) || !finite(time)) return 0.0;
    const double q = hbar * time / (2.0 * mass * sigma0 * sigma0);
    return sigma0 * std::sqrt(1.0 + q * q);
}
 
QuantumAndersonHistory2D::QuantumAndersonHistory2D(std::size_t maxSamples)
    : maxSamples_(std::max<std::size_t>(2u, maxSamples))
{
}
 
void QuantumAndersonHistory2D::clear()
{
    configured_ = false;
    config_ = QuantumAndersonInitialState2DConfig();
    lattice_ = QuantumAndersonLattice2D();
    samples_.clear();
}
 
void QuantumAndersonHistory2D::configure(const QuantumAndersonInitialState2DConfig &config,
                                         const QuantumAndersonLattice2D &lattice)
{
    clear();
    if (!lattice.valid || !QuantumAnderson2D::validateConfig(config, nullptr)) return;
    config_ = config;
    lattice_ = lattice;
    configured_ = true;
}
 
bool QuantumAndersonHistory2D::configured() const noexcept { return configured_; }
bool QuantumAndersonHistory2D::empty() const noexcept { return samples_.empty(); }
 
void QuantumAndersonHistory2D::append(double time, const ComplexField2D &field)
{
    if (!configured_ || field.empty()) return;
    QuantumAndersonSample2D sample;
    sample.time = time;
    sample.clean = QuantumAnderson2D::measureLane(field, lattice_, true);
    sample.disorder = QuantumAnderson2D::measureLane(field, lattice_, false);
    sample.norm = sample.clean.probability + sample.disorder.probability;
    samples_.push_back(sample);
    trim();
}
 
const QuantumAndersonInitialState2DConfig &QuantumAndersonHistory2D::config() const noexcept { return config_; }
const QuantumAndersonLattice2D &QuantumAndersonHistory2D::lattice() const noexcept { return lattice_; }
const std::vector<QuantumAndersonSample2D> &QuantumAndersonHistory2D::samples() const noexcept { return samples_; }
 
void QuantumAndersonHistory2D::trim()
{
    if (samples_.size() <= maxSamples_) return;
    const std::size_t removeCount = samples_.size() - maxSamples_;
    samples_.erase(samples_.begin(), samples_.begin() + static_cast<std::ptrdiff_t>(removeCount));
}
 
} // namespace qpe