amap2_fit_gamma.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_MAM_AMAP2_FIT_GAMMA_H
#define LINE_API_MAM_AMAP2_FIT_GAMMA_H
 
/**
 * @file
 * @ingroup api_mam
 * AMAP(2) fit of three moments and the autocorrelation decay rate
 * (matlab/lib/m3a/m3a/amap2/amap2_fit_gamma.m).
 *
 * A unit SCV short-circuits to a Poisson process, since the AMAP(2) canonical
 * forms are degenerate there. Otherwise every exact solution is enumerated
 * with amap2_fitall_gamma, normalized, and the first is returned.
 *
 * When no exact solution exists the characteristics are relaxed into the
 * feasible region with amap2_adjust_gamma at the reference default (method 3:
 * priority M2 > M3 > GAMMA, closed form through aph2_adjust and the GAMMA
 * bounds) and refitted. As in the reference, those approximate solutions are
 * not passed through map_normalize. Only when the refit is empty as well does
 * the Poisson process matching the mean come back, flagged through
 * Amap2FitGammaResult::poisson_fallback.
 *
 * Gated on transcendental arithmetic through amap2_fitall_gamma.
 */
 
#include <vector>
 
#include "line/api/mam/amap2_adjust_gamma.h"
#include "line/api/mam/amap2_fitall_gamma.h"
#include "line/api/mam/map_moment.h"
#include "line/api/mam/map_transform.h"
#include "line/num/number.h"
 
namespace line {
namespace mam {
 
/** Result of amap2_fit_gamma. */
template <class T>
struct Amap2FitGammaResult {
    Map<T> amap;                ///< the selected fit
    std::vector<Map<T>> amaps;  ///< every solution found (of the adjusted characteristics when adjusted)
    bool adjusted;              ///< true when the fit is of characteristics relaxed by amap2_adjust_gamma
    bool poisson_fallback;      ///< true when SCV == 1, or no AMAP(2) exists even after adjustment
};
 
/**
 * Fit an AMAP(2) to (M1, M2, M3, GAMMA). cvtol is the tolerance on
 * |M2 - 2 M1^2| below which the Poisson short-circuit fires (MATLAB 1e-6).
 */
template <class T>
Amap2FitGammaResult<T> amap2_fit_gamma(const T& M1, const T& M2, const T& M3, const T& GAMMA,
                                       const T& cvtol) {
    static_assert(num_traits<T>::has_transcendental,
                  "amap2_fit_gamma requires transcendental arithmetic");
    const T two = num_traits<T>::from_int(2);
    Amap2FitGammaResult<T> r;
    r.adjusted = false;
    r.poisson_fallback = false;
 
    if (num_abs(T(M2 - two * M1 * M1)) < cvtol) {
        r.amap = map_exponential_mean(M1);
        r.amaps.push_back(r.amap);
        r.poisson_fallback = true;
        return r;
    }
 
    std::vector<Map<T>> all = amap2_fitall_gamma(M1, M2, M3, GAMMA);
    for (std::size_t j = 0; j < all.size(); ++j) all[j] = map_normalize(all[j]);
    r.amaps = all;
 
    if (r.amaps.empty()) {
        // approximate fit of the nearest feasible characteristics (reference default method 3)
        const Amap2AdjustGammaResult<T> adj = amap2_adjust_gamma(M1, M2, M3, GAMMA);
        r.amaps = amap2_fitall_gamma(M1, adj.M2a, adj.M3a, adj.GAMMAa);
        r.adjusted = true;
    }
 
    if (r.amaps.empty()) {
        r.amap = map_exponential_mean(M1);
        r.amaps.push_back(r.amap);
        r.poisson_fallback = true;
        return r;
    }
    r.amap = r.amaps.front();
    return r;
}
 
/** amap2_fit_gamma with the MATLAB default cvtol = 1e-6. */
template <class T>
Amap2FitGammaResult<T> amap2_fit_gamma(const T& M1, const T& M2, const T& M3, const T& GAMMA) {
    return amap2_fit_gamma(M1, M2, M3, GAMMA, T(num_traits<T>::from_double(1e-6)));
}
 
}  // namespace mam
}  // namespace line
 
#endif  // LINE_API_MAM_AMAP2_FIT_GAMMA_H