// Part of BNC, a utility for retrieving decoding and
// converting GNSS data streams from NTRIP broadcasters.
//
// Copyright (C) 2007
// German Federal Agency for Cartography and Geodesy (BKG)
// http://www.bkg.bund.de
// Czech Technical University Prague, Department of Geodesy
// http://www.fsv.cvut.cz
//
// Email: euref-ip@bkg.bund.de
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation, version 2.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

/* -------------------------------------------------------------------------
 * BKG NTRIP Client
 * -------------------------------------------------------------------------
 *
 * Class:      t_ssrQc
 *
 * Purpose:    Quality control of decoded SSR streams
 *
 * Created:    05-Oct-2026
 *
 * Changes:
 *
 * -----------------------------------------------------------------------*/

#include <algorithm>
#include <cmath>
#include <iomanip>

#include "ssrQc.h"
#include "bnccore.h"
#include "bncutils.h"
#include "bncconst.h"
#include "bncsettings.h"
#include "combination/bnccomb.h"

using namespace std;

namespace {
  const double MAX_CLK_JUMP   = 0.5;  // [m] clock change beyond its rate between consecutive epochs
  const double EPH_WARMUP     = 120.0; // [s] no ephemeris checks at the start
  const double CMP_SAMPL      = 30.0;  // [s] sampling of the stream comparison
  const double CMP_DELAY      = 60.0;  // [s] epochs compared when data this much newer arrived
  const double HIST_LENGTH    = 300.0; // [s] corrections kept for the comparison
  const double DEV_MIN        = 0.05;  // [m] deviating satellites: RMS above max(DEV_MIN,
  const double DEV_FACTOR     = 3.0;   //     DEV_FACTOR x median RMS of the system)

  double tKey(const bncTime& tt) {
    return tt.gpsw() * 604800.0 + tt.gpssec();
  }
  const char*  TYPE_NAME[]    = {"ORBIT", "CLOCK", "CODE_BIAS", "PHASE_BIAS"};

  double median(QVector<double> vv) {
    if (vv.isEmpty()) {
      return 0.0;
    }
    std::sort(vv.begin(), vv.end());
    int nn = vv.size();
    return (nn % 2) ? vv[nn/2] : 0.5 * (vv[nn/2-1] + vv[nn/2]);
  }

  bncTime currentTime() {
    int    week = 0;
    double sec  = 0.0;
    currentGPSWeeks(week, sec);
    return bncTime(week, sec);
  }

  // satellite list "G01 G05(12)", counts shown if not all epochs concerned
  QString satList(const QMap<QString, unsigned>& sats, unsigned maxNum = 20) {
    QString str;
    unsigned num = 0;
    QMapIterator<QString, unsigned> it(sats);
    while (it.hasNext()) {
      it.next();
      if (++num > maxNum) {
        str += QString(" ... (%1 satellites)").arg(sats.size());
        break;
      }
      str += " " + it.key();
      if (it.value() > 0) {
        str += QString("(%1)").arg(it.value());
      }
    }
    return str.isEmpty() ? " -" : str;
  }
}

// Constructor
////////////////////////////////////////////////////////////////////////////
t_ssrQc::t_ssrQc(const QString& dirName, double interval, const QStringList& streams,
                 const QString& reference) : _ephUser(true) {
  _dirName   = dirName;
  _reference = reference.isEmpty() ? QString("CONSENSUS") : reference;
  _interval = interval > 0.0 ? interval : 3600.0;
  // "ALL" or a list of streams; streams with a leading '!' are excluded
  // (from the reports and from the consensus)
  bool all = false;
  for (const QString& sta : streams) {
    if (sta.startsWith('!')) {
      _excluded.insert(sta.mid(1));
    }
    else if (sta.toUpper() == "ALL") {
      all = true;
    }
    else {
      _streams.insert(sta);
    }
  }
  if (all) {
    _streams.clear();
  }

  connect(this, SIGNAL(newMessage(QByteArray,bool)), BNC_CORE, SLOT(slotMessage(const QByteArray,bool)));
  connect(BNC_CORE, SIGNAL(newOrbCorrections(QList<t_orbCorr>)),
          this, SLOT(slotNewOrbCorrections(QList<t_orbCorr>)), Qt::DirectConnection);
  connect(BNC_CORE, SIGNAL(newClkCorrections(QList<t_clkCorr>)),
          this, SLOT(slotNewClkCorrections(QList<t_clkCorr>)), Qt::DirectConnection);
  connect(BNC_CORE, SIGNAL(newCodeBiases(QList<t_satCodeBias>)),
          this, SLOT(slotNewCodeBiases(QList<t_satCodeBias>)), Qt::DirectConnection);
  connect(BNC_CORE, SIGNAL(newPhaseBiases(QList<t_satPhaseBias>)),
          this, SLOT(slotNewPhaseBiases(QList<t_satPhaseBias>)), Qt::DirectConnection);
  // combined corrections (stream INTERNAL)
  if (bncComb::exists()) {
    connect(BNC_CMB, SIGNAL(newOrbCorrections(QList<t_orbCorr>)),
            this, SLOT(slotNewOrbCorrections(QList<t_orbCorr>)), Qt::DirectConnection);
    connect(BNC_CMB, SIGNAL(newClkCorrections(QList<t_clkCorr>)),
            this, SLOT(slotNewClkCorrections(QList<t_clkCorr>)), Qt::DirectConnection);
  }
}

// Stream selected for the QC
////////////////////////////////////////////////////////////////////////////
bool t_ssrQc::selected(const string& staID) const {
  QString sta = QString::fromStdString(staID);
  return !_excluded.contains(sta) && (_streams.isEmpty() || _streams.contains(sta));
}

// Analysis center of a stream: the 3 characters of the AC in mountpoints
// like SSRA00CNE1, otherwise the stream name itself
////////////////////////////////////////////////////////////////////////////
QString t_ssrQc::acOf(const QString& staID) {
  static const QRegularExpression re("^SSR[AC][0-9]{2}([A-Z0-9]{3})[0-9]$");
  QRegularExpressionMatch mm = re.match(staID);
  return mm.hasMatch() ? mm.captured(1) : staID;
}

// Stream part of the consensus: not excluded, and not a combination product
// (INTERNAL, IGS combinations), as these are derived from the other streams
////////////////////////////////////////////////////////////////////////////
bool t_ssrQc::votes(const QString& staID) const {
  return !_excluded.contains(staID) && staID != "INTERNAL" && acOf(staID) != "IGS";
}

// Daily QC file of a stream: <staID>_S_<YYYYDDD>0000_01D_QC.txt
////////////////////////////////////////////////////////////////////////////
ofstream& t_ssrQc::streamOut(const QString& staID, const bncTime& time) {
  unsigned year, month, day;
  time.civil_date(year, month, day);
  QDate date(year, month, day);
  QString name = QDir(_dirName).filePath(QString("%1_S_%2%3%4_01D_QC.txt").arg(staID)
                 .arg(year).arg(date.dayOfYear(), 3, 10, QChar('0')).arg("0000"));
  if (_outNames.value(staID) != name) {
    delete _outs.value(staID, 0);
    bool exists = QFile::exists(name);
    ofstream* out = new ofstream(name.toLocal8Bit().data(), ios_base::out | ios_base::app);
    if (!exists) {
      *out << "SSR QC of stream " << staID.toStdString() << ", report interval " << _interval << " s" << endl
           << "Checks: epochs and update intervals, gaps (interval > 2 x max(declared, median)), latency,"
           << " completeness w.r.t. healthy broadcast ephemerides," << endl
           << "IOD matching, code and phase bias signals, code biases of the reference signals, clock jumps > "
           << MAX_CLK_JUMP << " m, largest orbit corrections" << endl;
    }
    _outs[staID]     = out;
    _outNames[staID] = name;
  }
  return *_outs[staID];
}

// Destructor
////////////////////////////////////////////////////////////////////////////
t_ssrQc::~t_ssrQc() {
  QMutexLocker locker(&_mutex);
  writeReport(true);
  for (ofstream* out : _outs) {
    delete out;
  }
}

// Key of the epoch statistics
////////////////////////////////////////////////////////////////////////////
QString t_ssrQc::key(const string& staID, char sys, e_type type) {
  return QString("%1 %2 %3").arg(staID.c_str()).arg(sys).arg(int(type));
}

// Reference signals (IGS conventions, as in the combination)
////////////////////////////////////////////////////////////////////////////
bool t_ssrQc::refSignals(char sys, QString& sig1, QString& sig2,
                         t_frequency::type& frq1, t_frequency::type& frq2) {
  switch (sys) {
    case 'G': sig1 = "1W"; sig2 = "2W"; frq1 = t_frequency::G1; frq2 = t_frequency::G2; return true;
    case 'R': sig1 = "1P"; sig2 = "2P"; frq1 = t_frequency::R1; frq2 = t_frequency::R2; return true;
    case 'E': sig1 = "1C"; sig2 = "5Q"; frq1 = t_frequency::E1; frq2 = t_frequency::E5; return true;
    case 'C': sig1 = "2I"; sig2 = "6I"; frq1 = t_frequency::C2; frq2 = t_frequency::C6; return true;
    case 'J': sig1 = "1C"; sig2 = "2L"; frq1 = t_frequency::J1; frq2 = t_frequency::J2; return true;
  }
  return false;
}

// Start a new report interval if the epoch is beyond the current one
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::checkInterval(const bncTime& time) {
  if (!_intervalStart.valid()) {
    double sec = floor(time.gpssec() / _interval) * _interval;
    _intervalStart = bncTime(time.gpsw(), sec);
  }
  // stream comparison of the epochs already complete for all streams
  if (!_nextCmpEpoch.valid()) {
    double sec = ceil(time.gpssec() / CMP_SAMPL) * CMP_SAMPL;
    _nextCmpEpoch = bncTime(time.gpsw(), sec);
  }
  while (time - _nextCmpEpoch >= CMP_DELAY) {
    compareEpoch(_nextCmpEpoch);
    _nextCmpEpoch = _nextCmpEpoch + CMP_SAMPL;
  }
  if (time - _intervalStart >= _interval) {
    writeReport(false);
    double sec = floor(time.gpssec() / _interval) * _interval;
    _intervalStart = bncTime(time.gpsw(), sec);
  }
  if (!_lastTime.valid() || _lastTime < time) {
    _lastTime = time;
  }
}

// Epoch statistics
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::newEpoch(const string& staID, char sys, e_type type, const bncTime& time,
                       unsigned updateInt, const QString& prn) {
  QString sta = QString::fromStdString(staID);
  if (!_staStart.contains(sta)) {
    _staStart[sta] = time;
  }
  if (!_staLast.contains(sta) || _staLast[sta] < time) {
    _staLast[sta] = time;
  }
  t_epoStat& stat = _epoStat[key(staID, sys, type)];
  if (updateInt < unsigned(ssrUpdateInt.size())) {
    stat.nominal = ssrUpdateInt[updateInt];
  }
  stat.sats.insert(prn);
  if (!stat.last.valid() || stat.last < time) {
    if (stat.last.valid()) {
      stat.intervals.push_back(time - stat.last);
    }
    stat.last = time;
    stat.numEpo += 1;
    stat.latencies.push_back(currentTime() - time);
  }
  else if (time < stat.last) {
    stat.numOutOfOrder += 1;
  }
}

// Broadcast ephemeris the correction refers to
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::checkEph(t_satStat& stat, const QString& prnInt, unsigned iod, const bncTime& time,
                       const string& staID) {
  // not before all broadcast ephemerides may have been received
  if (time - _staStart[QString::fromStdString(staID)] < EPH_WARMUP) {
    return;
  }
  t_eph* ephLast = _ephUser.ephLast(prnInt);
  t_eph* ephPrev = _ephUser.ephPrev(prnInt);
  t_eph* eph = 0;
  if      (ephLast && ephLast->IOD() == iod) eph = ephLast;
  else if (ephPrev && ephPrev->IOD() == iod) eph = ephPrev;
  if (!ephLast) {
    stat.noEph += 1;
  }
  else if (!eph) {
    stat.iodUnknown += 1;
  }
  else if (eph->checkState() == t_eph::unhealthy || eph->checkState() == t_eph::bad) {
    stat.ephUnhealthy += 1;
  }
}

// Orbit corrections
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::slotNewOrbCorrections(QList<t_orbCorr> orbCorrections) {
  QMutexLocker locker(&_mutex);
  for (const t_orbCorr& corr : orbCorrections) {
    if (!_excluded.contains(QString::fromStdString(corr._staID)) || corr._staID == _reference.toStdString()) {
      storeHist(&corr, 0);
    }
    if (!selected(corr._staID)) continue;
    QString prn = QString::fromStdString(corr._prn.toString());
    checkInterval(corr._time);
    newEpoch(corr._staID, corr._prn.system(), ORB, corr._time, corr._updateInt, prn);
    t_satStat& stat = _satStat[QString::fromStdString(corr._staID) + ' ' + prn];
    stat.numOrb += 1;
    checkEph(stat, QString::fromStdString(corr._prn.toInternalString()), corr._iod, corr._time, corr._staID);
    double radial = fabs(corr._xr[0]);
    double alongCross = max(fabs(corr._xr[1]), fabs(corr._xr[2]));
    stat.maxRadial     = max(stat.maxRadial, radial);
    stat.maxAlongCross = max(stat.maxAlongCross, alongCross);
  }
}

// Clock corrections
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::slotNewClkCorrections(QList<t_clkCorr> clkCorrections) {
  QMutexLocker locker(&_mutex);
  for (const t_clkCorr& corr : clkCorrections) {
    if (!_excluded.contains(QString::fromStdString(corr._staID)) || corr._staID == _reference.toStdString()) {
      storeHist(0, &corr);
    }
    if (!selected(corr._staID)) continue;
    QString prn = QString::fromStdString(corr._prn.toString());
    checkInterval(corr._time);
    newEpoch(corr._staID, corr._prn.system(), CLK, corr._time, corr._updateInt, prn);
    t_satStat& stat = _satStat[QString::fromStdString(corr._staID) + ' ' + prn];
    stat.numClk += 1;
    checkEph(stat, QString::fromStdString(corr._prn.toInternalString()), corr._iod, corr._time, corr._staID);

    // jump: change between consecutive epochs (same IOD) beyond the clock rate
    double clk = corr._dClk * t_CST::c;
    if (stat.lastClkTime.valid() && stat.lastIod == corr._iod) {
      double dt = corr._time - stat.lastClkTime;
      int nominal = (corr._updateInt < unsigned(ssrUpdateInt.size())) ? ssrUpdateInt[corr._updateInt] : 0;
      if (dt > 0.0 && dt <= max(2.0 * nominal, 60.0)) {
        double jump = fabs(clk - stat.lastClk - corr._dotDClk * t_CST::c * dt);
        stat.maxJump = max(stat.maxJump, jump);
        if (jump > MAX_CLK_JUMP) {
          stat.clkJumps += 1;
        }
      }
    }
    stat.lastClkTime = corr._time;
    stat.lastIod     = corr._iod;
    stat.lastClk     = clk;
  }
}

// Code biases
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::slotNewCodeBiases(QList<t_satCodeBias> codeBiases) {
  QMutexLocker locker(&_mutex);
  for (const t_satCodeBias& bias : codeBiases) {
    bool sel = selected(bias._staID);
    if (_excluded.contains(QString::fromStdString(bias._staID)) && bias._staID != _reference.toStdString()) continue;
    char sys = bias._prn.system();
    QString prn = QString::fromStdString(bias._prn.toString());
    if (sel) {
      checkInterval(bias._time);
      newEpoch(bias._staID, sys, CBIAS, bias._time, bias._updateInt, prn);
    }
    QString staSys = QString("%1 %2").arg(bias._staID.c_str()).arg(sys);
    QMap<QString, double> values;
    for (const t_frqCodeBias& frqBias : bias._bias) {
      QString sig = QString::fromStdString(frqBias._rnxType2ch);
      _codeSignals[staSys].insert(sig);
      values[sig] = frqBias._value;
    }
    t_satStat& stat = _satStat[QString::fromStdString(bias._staID) + ' ' + prn];
    QString sig1, sig2;
    t_frequency::type frq1, frq2;
    if (refSignals(sys, sig1, sig2, frq1, frq2)) {
      stat.missingRefSig.clear();
      if (!values.contains(sig1)) stat.missingRefSig.insert("C" + sig1);
      if (!values.contains(sig2)) stat.missingRefSig.insert("C" + sig2);
      stat.refBiasOK = stat.missingRefSig.isEmpty();
      if (stat.refBiasOK) {
        double channel = 0.0;
        if (sys == 'R') {
          t_eph* eph = _ephUser.ephLast(QString::fromStdString(bias._prn.toInternalString()));
          if (!eph) {
            stat.refBiasOK = false;
            continue;
          }
          channel = eph->slotNum();
        }
        double f1 = t_CST::freq(frq1, channel);
        double f2 = t_CST::freq(frq2, channel);
        stat.biasIF = (f1 * f1 * values[sig1] - f2 * f2 * values[sig2]) / (f1 * f1 - f2 * f2);
      }
    }
  }
}

// Phase biases
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::slotNewPhaseBiases(QList<t_satPhaseBias> phaseBiases) {
  QMutexLocker locker(&_mutex);
  for (const t_satPhaseBias& bias : phaseBiases) {
    if (!selected(bias._staID)) continue;
    checkInterval(bias._time);
    newEpoch(bias._staID, bias._prn.system(), PBIAS, bias._time, bias._updateInt,
             QString::fromStdString(bias._prn.toString()));
    QString staSys = QString("%1 %2").arg(bias._staID.c_str()).arg(bias._prn.system());
    for (const t_frqPhaseBias& frqBias : bias._bias) {
      _phaseSignals[staSys].insert(QString::fromStdString(frqBias._rnxType2ch));
    }
  }
}

// Reference point of the orbits of a stream: APC (SSRA...), CoM (SSRC...);
// the combined stream INTERNAL refers to the APC
////////////////////////////////////////////////////////////////////////////
char t_ssrQc::refPoint(const QString& staID) const {
  if (staID == "INTERNAL") {
    return 'A';  // the combined corrections always refer to the APC
  }
  const QString& sta = staID;
  if (sta.startsWith("SSRA")) return 'A';
  if (sta.startsWith("SSRC")) return 'C';
  return '?';
}

// Keep the recent orbit and clock corrections of a stream
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::storeHist(const t_orbCorr* orb, const t_clkCorr* clk) {
  if (orb) {
    QMap<double, t_orbCorr>& hist = _orbHist[QString::fromStdString(orb->_staID)]
                                            [QString::fromStdString(orb->_prn.toInternalString())];
    hist[tKey(orb->_time)] = *orb;
    while (!hist.isEmpty() && hist.firstKey() < tKey(orb->_time) - HIST_LENGTH) {
      hist.erase(hist.begin());
    }
  }
  if (clk) {
    QMap<double, t_clkCorr>& hist = _clkHist[QString::fromStdString(clk->_staID)]
                                            [QString::fromStdString(clk->_prn.toInternalString())];
    hist[tKey(clk->_time)] = *clk;
    while (!hist.isEmpty() && hist.firstKey() < tKey(clk->_time) - HIST_LENGTH) {
      hist.erase(hist.begin());
    }
  }
}

// Compare the orbits and clocks of all streams at epoch tt
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::compareEpoch(const bncTime& tt) {
  double tk = tKey(tt);

  // Satellite positions and clocks of all streams (latest corrections <= tt)
  QMap<QString, QMap<QString, t_satCrd> > crd;   // stream, PRN (external)
  QMapIterator<QString, QMap<QString, QMap<double, t_clkCorr> > > itSta(_clkHist);
  while (itSta.hasNext()) {
    itSta.next();
    const QString& staID = itSta.key();
    QMapIterator<QString, QMap<double, t_clkCorr> > itPrn(itSta.value());
    while (itPrn.hasNext()) {
      itPrn.next();
      const QString& prnInt = itPrn.key();
      QMap<double, t_clkCorr>::const_iterator itC = itPrn.value().upperBound(tk);
      if (itC == itPrn.value().constBegin()) {
        continue;
      }
      const t_clkCorr& clk = (--itC).value();
      const QMap<double, t_orbCorr>& orbHist = _orbHist[staID][prnInt];
      QMap<double, t_orbCorr>::const_iterator itO = orbHist.upperBound(tk);
      if (itO == orbHist.constBegin()) {
        continue;
      }
      const t_orbCorr& orb = (--itO).value();
      t_eph* eph = 0;
      t_eph* ephLast = _ephUser.ephLast(prnInt);
      t_eph* ephPrev = _ephUser.ephPrev(prnInt);
      if      (ephLast && ephLast->IOD() == orb._iod) eph = ephLast;
      else if (ephPrev && ephPrev->IOD() == orb._iod) eph = ephPrev;
      if (!eph || clk._iod != orb._iod) {
        continue;
      }
      eph->setOrbCorr(&orb);
      eph->setClkCorr(&clk);
      ColumnVector xc(6), vv(3);
      if (eph->getCrd(tt, xc, vv, true) != success) {
        continue;
      }
      t_satCrd& sc = crd[staID][QString::fromStdString(clk._prn.toString())];
      sc.xx  = xc.Rows(1,3);
      sc.vv  = vv;
      sc.clk = xc(4) * t_CST::c;
      sc.refPoint = refPoint(staID);
      // clock referred to the reference signals (as in the combination):
      // + B_IF of the stream's code biases; the combination (INTERNAL)
      // refers to them already
      sc.clkOK = (staID == "INTERNAL");
      QString satKey = staID + ' ' + QString::fromStdString(clk._prn.toString());
      if (!sc.clkOK && _satStat.contains(satKey) && _satStat[satKey].refBiasOK) {
        sc.clk  += _satStat[satKey].biasIF;
        sc.clkOK = true;
      }
    }
  }

  // Satellites per system
  QMap<char, QSet<QString> > satsSys;
  QMapIterator<QString, QMap<QString, t_satCrd> > itCrd(crd);
  while (itCrd.hasNext()) {
    itCrd.next();
    for (const QString& prn : itCrd.value().keys()) {
      satsSys[prn[0].toLatin1()].insert(prn);
    }
  }

  QMapIterator<char, QSet<QString> > itSys(satsSys);
  while (itSys.hasNext()) {
    itSys.next();
    const QSet<QString>& sats = itSys.value();
    QStringList staList;
    for (const QString& staID : crd.keys()) {
      for (const QString& prn : sats) {
        if (crd[staID].contains(prn)) {
          staList << staID;
          break;
        }
      }
    }

    // Consensus: one value per AC (median of its streams), so that an AC
    // with several streams (formats, APC/CoM) has one vote; the streams of
    // the AC compared are left out (leave-one-AC-out)
    auto acClocks = [&](const QString& prn, const QMap<QString, double>& off, bool useOff) {
      QMap<QString, QVector<double> > perAC;
      for (const QString& oo : staList) {
        if (!votes(oo) || !crd[oo].contains(prn) || !crd[oo][prn].clkOK) continue;
        if (useOff && !off.contains(oo)) continue;
        perAC[acOf(oo)] << crd[oo][prn].clk - (useOff ? off[oo] : 0.0);
      }
      QMap<QString, double> res;
      for (const QString& ac : perAC.keys()) res[ac] = median(perAC[ac]);
      return res;
    };

    // Clock datum: offset of each stream w.r.t. the per-satellite consensus
    // of all ACs (or w.r.t. the reference stream)
    QMap<QString, double> offset;
    QMap<QString, double> noOff;
    QMap<QString, QMap<QString, double> > rawAC;  // PRN, AC
    if (_reference == "CONSENSUS") {
      for (const QString& prn : sats) rawAC[prn] = acClocks(prn, noOff, false);
    }
    for (const QString& staID : staList) {
      QVector<double> diffs;
      for (const QString& prn : sats) {
        if (!crd[staID].contains(prn) || !crd[staID][prn].clkOK) continue;
        if (_reference == "CONSENSUS") {
          if (rawAC[prn].size() >= 3) diffs << crd[staID][prn].clk - median(rawAC[prn].values().toVector());
        }
        else if (crd.value(_reference).contains(prn) && crd[_reference][prn].clkOK) {
          diffs << crd[staID][prn].clk - crd[_reference][prn].clk;
        }
      }
      if (!diffs.isEmpty()) {
        offset[staID] = median(diffs);
      }
    }
    QMap<QString, QMap<QString, double> > alignedAC;  // PRN, AC
    if (_reference == "CONSENSUS") {
      for (const QString& prn : sats) alignedAC[prn] = acClocks(prn, offset, true);
    }

    // Differences of each selected stream to its reference
    for (const QString& staID : staList) {
      if (!selected(staID.toStdString()) || staID == _reference) {
        continue;
      }
      const QString ownAC = acOf(staID);
      QMap<QString, double> dClk;
      QMap<QString, ColumnVector> dRac;
      for (const QString& prn : sats) {
        if (!crd[staID].contains(prn)) continue;
        const t_satCrd& sc = crd[staID][prn];
        QVector<double> clks, xs[3];
        if (_reference == "CONSENSUS") {
          QMapIterator<QString, double> itA(alignedAC[prn]);
          while (itA.hasNext()) {
            itA.next();
            if (itA.key() != ownAC) clks << itA.value();
          }
          // orbits: per AC the median of its streams with the same reference point
          QMap<QString, QVector<double> > perAC[3];
          for (const QString& oo : staList) {
            if (!votes(oo) || acOf(oo) == ownAC || !crd[oo].contains(prn)) continue;
            const t_satCrd& rc = crd[oo][prn];
            if (rc.refPoint != sc.refPoint || sc.refPoint == '?') continue;
            for (int ii = 0; ii < 3; ii++) perAC[ii][acOf(oo)] << rc.xx[ii];
          }
          for (const QString& ac : perAC[0].keys()) {
            for (int ii = 0; ii < 3; ii++) xs[ii] << median(perAC[ii][ac]);
          }
          if (clks.size() < 2) clks.clear();
          if (xs[0].size() < 2) for (int ii = 0; ii < 3; ii++) xs[ii].clear();
        }
        else {
          if (!crd.value(_reference).contains(prn)) continue;
          const t_satCrd& rc = crd[_reference][prn];
          if (rc.clkOK) clks << rc.clk;
          if (rc.refPoint == sc.refPoint && sc.refPoint != '?') {
            for (int ii = 0; ii < 3; ii++) xs[ii] << rc.xx[ii];
          }
        }
        if (sc.clkOK && offset.contains(staID) && !clks.isEmpty()) {
          dClk[prn] = sc.clk - offset[staID] - median(clks);
        }
        else if (!sc.clkOK) {
          _cmpStat[staID][prn].clkNotReferred += 1;
        }
        if (!xs[0].isEmpty()) {
          ColumnVector dx(3);
          for (int ii = 0; ii < 3; ii++) dx[ii] = sc.xx[ii] - median(xs[ii]);
          ColumnVector er = sc.xx / sc.xx.NormFrobenius();
          ColumnVector ec = crossproduct(sc.xx, sc.vv);
          ec /= ec.NormFrobenius();
          ColumnVector ea = crossproduct(ec, er);
          ColumnVector rac(3);
          rac[0] = DotProduct(dx, er);
          rac[1] = DotProduct(dx, ea);
          rac[2] = DotProduct(dx, ec);
          dRac[prn] = rac;
        }
      }
      if (dClk.isEmpty() && dRac.isEmpty()) {
        continue;
      }
      // remaining common clock offset of this stream and epoch
      double off = dClk.isEmpty() ? 0.0 : median(dClk.values().toVector());
      _cmpEpochs[staID] += 1;
      QSet<QString> prns;
      for (const QString& prn : dClk.keys()) prns.insert(prn);
      for (const QString& prn : dRac.keys()) prns.insert(prn);
      for (const QString& prn : prns) {
        t_cmpStat& stat = _cmpStat[staID][prn];
        if (dRac.contains(prn)) {
          const ColumnVector& rac = dRac[prn];
          stat.numOrb += 1;
          stat.sR   += rac[0] * rac[0];
          stat.sA   += rac[1] * rac[1];
          stat.sC   += rac[2] * rac[2];
        }
        if (dClk.contains(prn)) {
          double dc = dClk[prn] - off;
          stat.numClk += 1;
          stat.sClk   += dc * dc;
          stat.mClk   += dc;
          if (dRac.contains(prn)) {
            double sis = dRac[prn][0] - dc;
            stat.numSis += 1;
            stat.sSis   += sis * sis;
            stat.mSis   += sis;
          }
        }
      }
    }
  }
}

// Write the report of the current interval and reset the statistics
////////////////////////////////////////////////////////////////////////////
void t_ssrQc::writeReport(bool final) {
  if (_epoStat.isEmpty()) {
    return;
  }
  bncTime endTime = final ? _lastTime : _intervalStart + _interval;

  // Healthy and current broadcast ephemerides per system
  QMap<char, QSet<QString> > ephSats;
  for (const QString& prnInt : _ephUser.prnList()) {
    t_eph* eph = _ephUser.ephLast(prnInt);
    if (eph && eph->checkState() == t_eph::ok && !outDatedBcep(eph, endTime)) {
      ephSats[eph->prn().system()].insert(QString::fromStdString(eph->prn().toString()));
    }
  }

  // Streams and systems
  QMap<QString, QSet<char> > streams;
  QMapIterator<QString, t_epoStat> itEpo(_epoStat);
  while (itEpo.hasNext()) {
    itEpo.next();
    QStringList hlp = itEpo.key().split(' ');
    streams[hlp[0]].insert(hlp[1][0].toLatin1());
  }

  QMapIterator<QString, QSet<char> > itSta(streams);
  while (itSta.hasNext()) {
    itSta.next();
    const QString& staID = itSta.key();
    QList<char> systems = itSta.value().values();
    std::sort(systems.begin(), systems.end());
    ofstream& out = streamOut(staID, _intervalStart);
    if (final) {
      endTime = _staLast.value(staID, _lastTime);
    }
    out << endl << "=== SSR QC " << staID.toStdString() << ' '
        << _intervalStart.datestr() << ' ' << _intervalStart.timestr(0)
        << " - " << endTime.datestr() << ' ' << endTime.timestr(0)
        << (final ? " (end of data)" : "") << " ===" << endl
        << "  Sys Type        Epochs  Interval nom/med/max [s]  Gaps  Out-of-order"
          << "  Latency med/max [s]  Sats" << endl;
    for (char sys : systems) {
      for (int type = ORB; type <= PBIAS; type++) {
        QString kk = QString("%1 %2 %3").arg(staID).arg(sys).arg(type);
        if (!_epoStat.contains(kk)) {
          continue;
        }
        const t_epoStat& stat = _epoStat[kk];
        double medInt = median(stat.intervals);
        double maxInt = stat.intervals.isEmpty() ? 0.0 :
                        *std::max_element(stat.intervals.begin(), stat.intervals.end());
        double ref = max(double(stat.nominal), medInt);
        unsigned gaps = 0;
        for (double dt : stat.intervals) {
          if (ref > 0.0 && dt > 2.0 * ref) {
            ++gaps;
          }
        }
        double maxLat = stat.latencies.isEmpty() ? 0.0 :
                        *std::max_element(stat.latencies.begin(), stat.latencies.end());
        out << "  " << sys << "   " << left << setw(10) << TYPE_NAME[type] << right
              << setw(8) << stat.numEpo
              << setw(10) << stat.nominal << " /" << setw(6) << setprecision(1) << fixed << medInt
              << " /" << setw(7) << maxInt
              << setw(6) << gaps << setw(14) << stat.numOutOfOrder
              << setw(10) << median(stat.latencies) << " /" << setw(7) << maxLat
              << setw(7) << stat.sats.size() << endl;
      }
    }

    // Per system: completeness, ephemeris, biases, plausibility
    for (char sys : systems) {
      QMap<QString, unsigned> missOrb, missClk, missCB, missPB, iodUnknown, unhealthy, noEph,
                              jumps, missRef;
      double maxRadial = 0.0, maxAlongCross = 0.0;
      QString maxRadialPrn, maxAlongCrossPrn;
      QVector<double> biasIF;
      const QSet<QString>& ephs = ephSats[sys];
      for (int type = ORB; type <= PBIAS; type++) {
        QString kk = QString("%1 %2 %3").arg(staID).arg(sys).arg(type);
        if (!_epoStat.contains(kk)) {
          continue;
        }
        QMap<QString, unsigned>& miss = (type == ORB) ? missOrb : (type == CLK) ? missClk :
                                        (type == CBIAS) ? missCB : missPB;
        for (const QString& prn : ephs) {
          if (!_epoStat[kk].sats.contains(prn)) {
            miss[prn] = 0;
          }
        }
      }
      QMapIterator<QString, t_satStat> itSat(_satStat);
      while (itSat.hasNext()) {
        itSat.next();
        QStringList hlp = itSat.key().split(' ');
        if (hlp[0] != staID || hlp[1][0].toLatin1() != sys) {
          continue;
        }
        const QString& prn = hlp[1];
        const t_satStat& stat = itSat.value();
        if (stat.iodUnknown)   iodUnknown[prn] = stat.iodUnknown;
        if (stat.ephUnhealthy) unhealthy[prn]  = stat.ephUnhealthy;
        if (stat.noEph)        noEph[prn]      = stat.noEph;
        if (stat.clkJumps)     jumps[prn]      = stat.clkJumps;
        if (stat.maxRadial > maxRadial) {
          maxRadial    = stat.maxRadial;
          maxRadialPrn = prn;
        }
        if (stat.maxAlongCross > maxAlongCross) {
          maxAlongCross    = stat.maxAlongCross;
          maxAlongCrossPrn = prn;
        }
        if (!stat.missingRefSig.isEmpty()) {
          QStringList sigs = stat.missingRefSig.values();
          sigs.sort();
          missRef[prn + "[" + sigs.join("/") + "]"] = 0;
        }
        if (stat.refBiasOK) {
          biasIF.push_back(stat.biasIF);
        }
      }
      // only findings are listed
      QStringList lines;
      for (int type = ORB; type <= PBIAS; type++) {
        QString kk = QString("%1 %2 %3").arg(staID).arg(sys).arg(type);
        if (_epoStat.contains(kk)) {
          const t_epoStat& stat = _epoStat[kk];
          double medInt = median(stat.intervals);
          if (stat.nominal > 0 && stat.intervals.size() > 1 &&
              (medInt > 1.5 * stat.nominal || medInt < stat.nominal / 1.5)) {
            lines << QString("%1: declared update interval %2 s, actual %3 s")
                     .arg(TYPE_NAME[type]).arg(stat.nominal).arg(medInt, 0, 'f', 1);
          }
        }
      }
      auto addSats = [&lines](const QString& text, const QMap<QString, unsigned>& sats) {
        if (!sats.isEmpty()) {
          lines << text + ":" + satList(sats);
        }
      };
      addSats("orbit corrections missing (healthy broadcast ephemeris available)", missOrb);
      addSats("clock corrections missing (healthy broadcast ephemeris available)", missClk);
      addSats("code biases missing (healthy broadcast ephemeris available)", missCB);
      addSats("phase biases missing (healthy broadcast ephemeris available)", missPB);
      addSats("corrections referring to an IOD of no received broadcast ephemeris", iodUnknown);
      addSats("corrections for unhealthy broadcast ephemeris", unhealthy);
      addSats("corrections for satellites without broadcast ephemeris", noEph);
      QString staSys = QString("%1 %2").arg(staID).arg(sys);
      if (_codeSignals.contains(staSys)) {
        QStringList sigs = _codeSignals[staSys].values();
        sigs.sort();
        lines << "code bias signals: " + sigs.join(" ");
        addSats("code biases of the reference signals missing", missRef);
        if (biasIF.size() > 1) {
          double mean = 0.0;
          for (double bb : biasIF) mean += bb / biasIF.size();
          double var = 0.0;
          for (double bb : biasIF) var += (bb - mean) * (bb - mean) / biasIF.size();
          double maxAbs = max(fabs(*std::min_element(biasIF.begin(), biasIF.end())),
                              fabs(*std::max_element(biasIF.begin(), biasIF.end())));
          lines << QString("ionosphere-free combination B_IF of the reference signal biases: "
                           "std. dev. %1 m, max |B_IF| %2 m (%3 satellites)")
                   .arg(sqrt(var), 0, 'f', 3).arg(maxAbs, 0, 'f', 3).arg(biasIF.size());
        }
      }
      if (_phaseSignals.contains(staSys)) {
        QStringList sigs = _phaseSignals[staSys].values();
        sigs.sort();
        lines << "phase bias signals: " + sigs.join(" ");
      }
      addSats(QString("clock jumps > %1 m").arg(MAX_CLK_JUMP), jumps);
      if (!maxRadialPrn.isEmpty()) {
        lines << QString("largest orbit corrections: radial %1 m (%2), along/cross %3 m (%4)")
                 .arg(maxRadial, 0, 'f', 2).arg(maxRadialPrn).arg(maxAlongCross, 0, 'f', 2).arg(maxAlongCrossPrn);
      }
      // stream comparison
      if (_cmpStat.contains(staID)) {
        unsigned nSat = 0, nClk = 0, nOrb = 0, nSis = 0;
        double sClk = 0.0, sR = 0.0, sA = 0.0, sC = 0.0, sSis = 0.0;
        double vClk = 0.0, vSis = 0.0;  // sums of squares about the satellite means
        QMap<QString, double> rmsSis, rmsClk, rmsOrb;
        QMap<QString, unsigned> notReferred;
        QMapIterator<QString, t_cmpStat> itC(_cmpStat[staID]);
        while (itC.hasNext()) {
          itC.next();
          if (itC.key()[0].toLatin1() != sys) continue;
          const t_cmpStat& st = itC.value();
          if (st.numClk == 0 && st.numOrb == 0) {
            if (st.clkNotReferred) notReferred[itC.key()] = 0;
            continue;
          }
          ++nSat;
          if (st.clkNotReferred) notReferred[itC.key()] = 0;
          nClk += st.numClk; sClk += st.sClk;
          nOrb += st.numOrb; sR += st.sR; sA += st.sA; sC += st.sC;
          nSis += st.numSis; sSis += st.sSis;
          if (st.numClk) vClk += st.sClk - st.mClk * st.mClk / st.numClk;
          if (st.numSis) vSis += st.sSis - st.mSis * st.mSis / st.numSis;
          if (st.numSis) rmsSis[itC.key()] = sqrt(st.sSis / st.numSis);
          else if (st.numClk) rmsClk[itC.key()] = sqrt(st.sClk / st.numClk);
          if (st.numOrb) rmsOrb[itC.key()] = sqrt((st.sR + st.sA + st.sC) / st.numOrb);
        }
        if (nSat) {
          QString refText = (_reference == "CONSENSUS") ? QString("consensus of the other ACs")
                                                        : QString("stream %1").arg(_reference);
          QStringList parts;
          if (nClk) parts << QString("clock %1 (STD %2)").arg(sqrt(sClk / nClk), 0, 'f', 3)
                                                          .arg(sqrt(max(0.0, vClk) / nClk), 0, 'f', 3);
          if (nOrb) parts << QString("radial %1, along %2, cross %3")
                             .arg(sqrt(sR / nOrb), 0, 'f', 3).arg(sqrt(sA / nOrb), 0, 'f', 3)
                             .arg(sqrt(sC / nOrb), 0, 'f', 3);
          if (nSis) parts << QString("SIS (radial - clock) %1 (STD %2)").arg(sqrt(sSis / nSis), 0, 'f', 3)
                                                                         .arg(sqrt(max(0.0, vSis) / nSis), 0, 'f', 3);
          lines << QString("comparison with %1 (%2 satellites), RMS (STD about the satellite means) [m]: %3")
                   .arg(refText).arg(nSat).arg(parts.join(", "));
          if (!nOrb) {
            lines << "orbits not compared: no other stream with the same reference point (APC/CoM)";
          }
          addSats("clocks not compared (code biases of the reference signals missing)", notReferred);
          auto deviating = [&lines](const QString& what, const QMap<QString, double>& rms) {
            if (rms.size() < 2) return;
            double med = median(rms.values().toVector());
            QStringList dev;
            QMapIterator<QString, double> itR(rms);
            while (itR.hasNext()) {
              itR.next();
              if (itR.value() > max(DEV_MIN, DEV_FACTOR * med)) {
                dev << QString("%1(%2)").arg(itR.key()).arg(itR.value(), 0, 'f', 2);
              }
            }
            if (!dev.isEmpty()) {
              lines << QString("satellites deviating from the reference, RMS %1 [m]: %2").arg(what).arg(dev.join(" "));
            }
          };
          deviating("SIS", rmsSis);
          deviating("clock", rmsClk);
          deviating("orbit (3D)", rmsOrb);
        }
      }
      for (const QString& line : lines) {
        out << "  " << sys << " " << line.toStdString() << endl;
      }
    }
    out.flush();
  }

  emit newMessage(QString("SSR QC: report %1 %2 of %3 stream(s) written to %4")
                  .arg(QString::fromStdString(_intervalStart.datestr()))
                  .arg(QString::fromStdString(_intervalStart.timestr(0)))
                  .arg(streams.size()).arg(_dirName).toLatin1(), false);

  _cmpStat.clear();
  _cmpEpochs.clear();

  // reset, keeping the state needed across intervals
  for (QMap<QString, t_epoStat>::iterator it = _epoStat.begin(); it != _epoStat.end(); ++it) {
    t_epoStat& stat = it.value();
    bncTime last = stat.last;
    int nominal = stat.nominal;
    stat = t_epoStat();
    stat.last = last;
    stat.nominal = nominal;
  }
  for (QMap<QString, t_satStat>::iterator it = _satStat.begin(); it != _satStat.end(); ++it) {
    t_satStat& stat = it.value();
    t_satStat fresh;
    fresh.lastClkTime   = stat.lastClkTime;
    fresh.lastIod       = stat.lastIod;
    fresh.lastClk       = stat.lastClk;
    fresh.missingRefSig = stat.missingRefSig;
    fresh.refBiasOK     = stat.refBiasOK;
    fresh.biasIF        = stat.biasIF;
    stat = fresh;
  }
}