BmnAlignPriors.cxx

Go to the documentation of this file.

/*
    BM@N alignment routine
    BM@N experiment at NICA complex, JINR, 2025
 
    Department: Math & Soft Group of HEP lab
    Author:     Claude (Opus 5)
    Supervisor: Igor Polev, polev@jinr.int
 
    BmnAlignPriors class implementation
*/
 
#include "BmnAlignPriors.h"
 
#include "TDecompSVD.h"
#include "TString.h"
#include "fairlogger/Logger.h"
 
#include <algorithm>
#include <cmath>
#include <utility>
#include <vector>
 
namespace
{
// relative to the norm of a constraint row: below it the row adds nothing
constexpr Double_t ROW_RANK_TOL{1.0E-9};
// cm or rad per unit row norm: a dropped row asking for more contradicts
constexpr Double_t CONTRADICTION_TOL{1.0E-9};
}   // namespace
 
/*
 *  Fold the prior knowledge (presigmas) into the reduced system of normal
 *  equations C' dA = b' produced by the Millepede dimension reduction.
 *
 *  A presigma is a pseudo-measurement of a single global parameter:
 *
 *      "parameter k equals reference_k, with uncertainty sigma_k"
 *
 *  Its residual is (reference_k - A_k), its derivative with respect to A_k is
 *  exactly 1 (and zero for every other global and local parameter), and its
 *  weight is 1/sigma_k^2. A measurement with residual r, weight w and
 *  derivative row g contributes  C += w*g*g^T,  b += w*g*r  - and since g is
 *  the unit vector e_k, this collapses onto a single diagonal element and a
 *  single RHS entry.
 *
 *  Two details are essential here:
 *
 *  1. Each iteration solves for the step dA while A accumulates, so the
 *     residual must reference the *accumulated* correction of the previous
 *     iteration - exactly as the constraint RHS does in IterateAlignment().
 *     With A = A_acc + dA the residual becomes (reference - A_acc).
 *     CAUTION: A_acc must be stated in the basis the system is solved in. For
 *     an element that is its own share, while the accumulated solution holds
 *     its effective correction a_acc + J P_acc. Measuring a leaf prior against
 *     the effective value makes the composites drift: the data cannot see the
 *     level split, so every iteration moves P by about the effective value.
 *
 *  2. The weight is stated in the same units as the hit weights of
 *     ProcessTracks, which are plain 1/sigma^2. A presigma of sigma therefore
 *     acts precisely like one extra hit measuring that parameter with
 *     resolution sigma.
 *
 *  A single pseudo-measurement is weak against millions of real hits - and
 *  that is the intent: the prior is negligible for well determined parameters
 *  and dominates only where the data says (almost) nothing, i.e. along the
 *  weak / chi2-invariant modes it is meant to lift.
 */
void BmnAlignPriors::Apply(TVectorD& vB, TMatrixD& mC, const BmnAlignHierarchy& hierarchy, const TVectorD* accA) const
{
    if (!UsePresigmas())
        return;
 
    const TVectorD zero(hierarchy.AllParams());
    const TVectorD& acc{accA ? *accA : zero};
 
    // the ancestors' part of each element's accumulated correction, J P_acc
    TVectorD vJP(hierarchy.LeafParams());
    if (hierarchy.CompositeParams())
        vJP = hierarchy.J() * acc.GetSub(hierarchy.CompositeParamBase(), hierarchy.AllParams() - 1);
 
    for (const ParamRef& p : hierarchy.Layout().All()) {
        const Double_t s{SigmaAt(p.full)};
        if (s <= 0.0)
            continue;
        const Double_t w{1.0 / (s * s)}, own{p.IsComposite() ? acc(p.full) : acc(p.full) - vJP[p.full]};
        mC(p.full, p.full) += w;                 // Hessian : + 1/sigma^2
        vB[p.full] += w * (ref[p.full] - own);   // gradient: + (ref - acc_own)/sigma^2
    }
}
 
Bool_t BmnAlignPriors::Check(const BmnAlignHierarchy& hierarchy)
{
    /*
     *  Constraint rows are kept only if they add something to the rows before
     *  them, tested by modified Gram-Schmidt on the promoted rows. A repeated
     *  row makes the bordered matrix singular, and a config may repeat itself.
     *  A row dropped so must ask for the value the kept rows already give; if
     *  it does not, the constraints contradict each other.
     *
     *  The test does not depend on the iteration. The promotion is c -> c S
     *  with S = [ I J ; 0 I ] invertible, so a promoted row is a combination
     *  of other promoted rows exactly when the original row is the same
     *  combination of the original rows. The right-hand side the aligner uses,
     *  rhs - C A_acc, then reduces to that combination of the rhs from the file:
     *  the A_acc terms cancel.
     */
    const TMatrixD rows{hierarchy.PromoteConstraints(constraints)};
    std::vector<TVectorD> basis;   // orthonormal span of the rows kept so far
    VecDouble_t basisRhs;          // their right-hand sides, transformed alike
    Int_t contradictions{0};
    fKept.clear();
    for (Int_t j = 0; j < ConstCount(); j++) {
        const TVectorD row{TMatrixDRow_const(rows, j)};
        TVectorD res{row};
        Double_t resRhs{rhs[j]};
        for (std::size_t m = 0; m < basis.size(); m++) {
            const Double_t d{Dot(basis[m], res)};
            Add(res, -d, basis[m]);
            resRhs -= d * basisRhs[m];
        }
        const Double_t rowNorm{std::sqrt(row.Norm2Sqr())}, resNorm{std::sqrt(res.Norm2Sqr())};
        if (resNorm > ROW_RANK_TOL * rowNorm) {
            basis.push_back((1.0 / resNorm) * res);
            basisRhs.push_back(resRhs / resNorm);
            fKept.push_back(j);
            continue;
        }
        const Double_t scale{rowNorm > 0.0 ? rowNorm : 1.0};   // an all-zero row asks 0 = rhs
        if (std::fabs(resRhs) <= CONTRADICTION_TOL * scale)
            continue;
        LOG(error) << Form(" <!> BmnAlignPriors::Check: constraint row #%d asks for a value %.2e off what the"
                           " constraints before it give (or has an all-zero LHS with a non-zero RHS)",
                           j + 1, resRhs / scale);
        contradictions++;
    }
    if (contradictions) {
        LOG(error) << " <!> BmnAlignPriors::Check: " << contradictions << " contradicting constraint(s)";
        return kFALSE;
    }
    fKeptRows.ResizeTo(KeptCount(), fAllParams);
    for (Int_t i = 0; i < KeptCount(); i++)
        TMatrixDRow(fKeptRows, i) = TMatrixDRow_const(rows, fKept[i]);
    LOG_IF(info, ConstCount()) << "    " << KeptCount() << " constraint(s) kept, " << ConstCount() - KeptCount()
                               << " dropped as empty or repeated";
 
    /*
     *  Composite parameters are redundant by construction: T = [ I | J ] has a
     *  null space { (-J c, c) }, one dimension per composite parameter, and the
     *  data says exactly nothing along it. Something else has to determine
     *  those directions, and only three things can:
     *
     *    - a prior on leaf parameter l         sees the direction as (J c)_l
     *    - a prior on composite parameter a    sees it as c_a
     *    - a constraint                        sees it as c_direct^T c
     *
     *  The constraint case is the one worth spelling out: its row in the
     *  hierarchical basis is [ c_leaf ; J^T c_leaf + c_direct ], which on
     *  (-J c, c) evaluates to c_direct^T c - the leaf part cancels exactly. So
     *  a constraint written on leaf parameters is *blind* to the composite
     *  split, however many of them there are, because it constrains the
     *  effective corrections A = T P and T v = 0 means A does not move. Only
     *  coefficients placed directly on a composite (NodeName) count here.
     *
     *  Stacking those rows, the hierarchy is determined if they have full
     *  column rank. Only the presence of a prior matters, not how tight it is:
     *  any finite sigma lifts the exact degeneracy.
     */
    const Int_t nComp{hierarchy.CompositeParams()};
    if (!nComp)
        return kTRUE;   // flat parametrization leaves nothing extra undetermined
 
    const Int_t nFull{hierarchy.AllParams()};
    // sigma stays empty when no presigmas file was given, hence SigmaAt()
    TMatrixD probe(nFull + ConstCount(), nComp);
    probe.Zero();
    for (const ParamRef& p : hierarchy.Layout().All()) {
        if (SigmaAt(p.full) <= 0.0)
            continue;
        if (p.IsComposite())
            probe(p.full, p.comp) = 1.0;   // prior placed on the node itself
        else                               // leaf prior, seen through J
            TMatrixDRow(probe, p.full) = TMatrixDRow_const(hierarchy.J(), p.full);
    }
    if (ConstCount())
        probe.SetSub(nFull, 0, constraints.GetSub(0, ConstCount() - 1, hierarchy.CompositeParamBase(), nFull - 1));
 
    // An all-zero probe is rank 0 by inspection, and SVD does not survive it.
    // Worth its own message: it means nothing at all was configured, rather
    // than what was configured being insufficient.
    if (!probe.NonZeros()) {
        LOG(error) << " <!> none of the " << nComp
                   << " composite parameter(s) is"
                      " determined: neither a presigma nor a composite-addressed"
                      " constraint was given, so the split between the levels rests"
                      " on nothing at all";
        return kFALSE;
    }
    TDecompSVD svd(probe);
    if (!svd.Decompose()) {
        LOG(error) << " <!> BmnAlignPriors::Check: SVD failed, solvability not verified";
        return kFALSE;
    }
    const TVectorD& sig{svd.GetSig()};                  // TDecompSVD returns them in decreasing order
    constexpr static const Double_t RANK_TOL{1.0E-9};   // relative to the largest singular value
    Int_t rank{0};
    while (rank < sig.GetNrows() && sig[rank] > sig[0] * RANK_TOL)
        rank++;
 
    if (rank == nComp) {
        LOG(info) << "    hierarchy check: all " << nComp
                  << " composite parameter(s) determined by priors and constraints";
        return kTRUE;
    }
 
    // Report the undetermined directions themselves rather than the parameters
    // lacking a direct prior: leaf priors pin composite parameters indirectly
    // through J, so the two lists differ and only the former is actionable.
    // Singular values come back in decreasing order, so the trailing right
    // singular vectors span the null space.
    const TMatrixD& V{svd.GetV()};
    LOG(error) << " <!> " << nComp - rank << " of " << nComp
               << " composite parameter direction(s) left undetermined -"
                  " the fit would be singular along them:";
    for (Int_t d = rank; d < nComp; d++) {
        std::vector<std::pair<Double_t, Int_t>> weight;
        weight.reserve(nComp);
        for (Int_t a = 0; a < nComp; a++)
            weight.emplace_back(std::fabs(V(a, d)), a);
        std::sort(weight.rbegin(), weight.rend());
        TString descr;
        for (Int_t t = 0; t < nComp && weight[t].first > 0.1; t++)
            descr += Form("%s%+.2f x [%s]", descr.IsNull() ? "" : "  ", V(weight[t].second, d),
                          hierarchy.Layout().Composites()[weight[t].second].Label().Data());
        LOG(error) << "       " << descr;
    }
    LOG(error) << "     give the composites involved a presigma (NodeName entry) or a"
                  " NodeName constraint - leaf-addressed constraints cannot fix a"
                  " composite split, and a leaf presigma only reaches the combinations"
                  " that move its own element";
    return kFALSE;
}