Index: trunk/npstat/stat/OrthoPolyGOFTest1D.icc
===================================================================
--- trunk/npstat/stat/OrthoPolyGOFTest1D.icc	(revision 748)
+++ trunk/npstat/stat/OrthoPolyGOFTest1D.icc	(revision 749)
@@ -1,175 +1,176 @@
 #include <cassert>
 #include <algorithm>
 #include <stdexcept>
 #include <cmath>
 
 #include "npstat/nm/truncatedInverseSqrt.hh"
 #include "npstat/nm/sumOfSquares.hh"
 
 namespace npstat {
     namespace Private {
         template<class Quadrature, bool>
         struct PolyQuadIntegral
         {
             inline static double integral(const Quadrature& q,
                                           const npstat::AbsClassicalOrthoPoly1D& poly,
                                           const unsigned i, const unsigned k,
                                           const unsigned j, const unsigned m)
             {
                 return poly.directQuadrature(q, i, k, j, m);
             }
         };
 
         template<class Quadrature>
         struct PolyQuadIntegral<Quadrature, true>
         {
             inline static double integral(const Quadrature& q,
                                           const npstat::AbsClassicalOrthoPoly1D& poly,
                                           const unsigned i, const unsigned k,
                                           const unsigned j, const unsigned m)
             {
                 return poly.jointAverage(q, i, k, j, m);
             }
         };
     }
 
     template <class Quadrature>
     OrthoPolyGOFTest1D::OrthoPolyGOFTest1D(
         const AbsDistribution1D& distro,
         const Quadrature& q,
         const AbsClassicalOrthoPoly1D& poly,
         const unsigned i_maxdeg,
         const unsigned char* mask, unsigned lenMask)
         : AbsUnbinnedGOFTest1D(distro),
           poly_(poly.clone()),
           buf_((i_maxdeg+1U)*(i_maxdeg+2U)/2U - 1U),
           statBuf_(2U*i_maxdeg),
           sqrBuf_(2U*i_maxdeg),
           mask_(2U*i_maxdeg, 0),
           maxdeg_(i_maxdeg),
           nUnmasked_(2U*i_maxdeg)
     {
         calculateNormalizingMatrix(q, &invsqr_);
 
         if (lenMask)
         {
             assert(mask);
             if (lenMask != 2U*maxdeg_) throw std::invalid_argument(
                 "In npstat::OrthoPolyGOFTest1D constructor: "
                 "incompatible mask length");
             nUnmasked_ = 0;
             for (unsigned i=0; i<lenMask; ++i)
             {
                 mask_[i] = mask[i];
                 if (!mask[i])
                     ++nUnmasked_;
             }
             if (!nUnmasked_) throw std::invalid_argument(
                 "In npstat::OrthoPolyGOFTest1D constructor: "
                 "all eigenvectors are masked");
         }
     }
 
     template <class Quadrature>
     void OrthoPolyGOFTest1D::calculateNormalizingMatrix(
         const Quadrature& q, Matrix<lapack_double>* result) const
     {
         assert(result);
 
         const unsigned maxdegP1 = maxdeg_ + 1U;
         const unsigned nRows = maxdegP1*(maxdegP1+1U)/2U - 1U;
         Matrix<lapack_double> covmat(nRows, nRows);
 
         unsigned iRow = 0;
         for (unsigned i=0; i<=maxdeg_; ++i)
         {
             for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
             {
                 assert(iRow < nRows);
                 const double k_del_ik = k == i ? 1.0 : 0.0;
                 unsigned iCol = 0;
                 for (unsigned j=0; j<=maxdeg_; ++j)
                 {
                     for (unsigned m=j ? j : 1U; m<=maxdeg_; ++m, ++iCol)
                     {
                         if (iCol >= iRow)
                         {
                             assert(iCol < nRows);
                             const double k_del_jm = j == m ? 1.0 : 0.0;
                             double cov =
                                 Private::PolyQuadIntegral<Quadrature,Quadrature::interval_quadrature>::integral(q, *poly_, i, k, j, m) -
                                 k_del_ik*k_del_jm;
                             if (i == j && k == m && cov < 0.0)
                                 cov = 0.0;
                             covmat[iRow][iCol] = cov;
                         }
                     }
                 }
                 assert(iCol == nRows);
             }
         }
         assert(iRow == nRows);
 
         iRow = 0;
         for (unsigned i=0; i<=maxdeg_; ++i)
         {
             for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
             {
                 unsigned iCol = 0;
                 for (unsigned j=0; j<=maxdeg_; ++j)
                 {
                     for (unsigned m=j ? j : 1U; m<=maxdeg_; ++m, ++iCol)
                         if (iCol < iRow)
                             covmat[iRow][iCol] = covmat[iCol][iRow];
                 }       
             }
         }
 
         truncatedInverseSqrt(covmat, 2U*maxdeg_, 0.0, result);
         assert(result->nRows() == 2U*maxdeg_);
+        assert(result->nColumns() == nRows);
     }
 
     template<typename Numeric>
     unsigned OrthoPolyGOFTest1D::normalizedDeviations(
         const Numeric* data, const unsigned long lenData,
         double* deviations, const unsigned lenDeviations) const
     {
         assert(lenData);
         assert(data);
         assert(lenDeviations >= nUnmasked_);
         assert(deviations);
 
         const Matrix<double>& unnorm = poly_-> sampleProductAverages(
             data, lenData, maxdeg_);
         double* buf = const_cast<double*>(&buf_[0]);
 
         unsigned iRow = 0;
         for (unsigned i=0; i<=maxdeg_; ++i)
             for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
                 buf[iRow] = unnorm[i][k] - (i == k ? 1.0 : 0.0);
 
         assert(iRow == buf_.size());
         const unsigned twomaxdeg = 2U*maxdeg_;
         assert(sqrBuf_.size() == twomaxdeg);
         double* sqrBuf = const_cast<double*>(&sqrBuf_[0]);
         invsqr_.timesVector(buf, iRow, sqrBuf, twomaxdeg);
 
         const double sqrlen = std::sqrt(static_cast<double>(lenData));
         unsigned idev=0;
         for (unsigned i=0; i<twomaxdeg; ++i)
             if (!mask_[i])
                 deviations[idev++] = sqrBuf[i]*sqrlen;
         assert(idev == nUnmasked_);
         return nUnmasked_;
     }
 
     template<typename Numeric>
     double OrthoPolyGOFTest1D::testStat(
         const Numeric* data, const unsigned long lenData) const
     {
         double* statbuf = const_cast<double*>(&statBuf_[0]);
         const unsigned ndevs = normalizedDeviations(
             data, lenData, statbuf, statBuf_.size());
         return sumOfSquares(statbuf, ndevs);
     }
 }
Index: trunk/npstat/stat/OrthoPolyGOFTest1D.cc
===================================================================
--- trunk/npstat/stat/OrthoPolyGOFTest1D.cc	(revision 748)
+++ trunk/npstat/stat/OrthoPolyGOFTest1D.cc	(revision 749)
@@ -1,80 +1,97 @@
 #include <sstream>
+#include <cmath>
+#include <cfloat>
 
 #include "npstat/stat/OrthoPolyGOFTest1D.hh"
 #include "npstat/stat/Distributions1D.hh"
 
 namespace npstat {
     double OrthoPolyGOFTest1D::analyticPValue(
         const double stat, unsigned long /* sz */) const
     {
         Gamma1D chisq(0.0, 2.0, nUnmasked_/2.0);
         return chisq.exceedance(stat);
     }
 
     std::string OrthoPolyGOFTest1D::shortName() const
     {
         const unsigned twodeg = 2U*maxdeg_;
 
         std::ostringstream os;
         os << "OP_" << twodeg;
         if (twodeg != nUnmasked_)
         {
             // Figure out the number of the last masked eigenvector
             unsigned lastMasked = twodeg;
             for (unsigned i=0; i<twodeg; ++i)
                 if (mask_[i])
                     lastMasked = i;
             assert(lastMasked < twodeg);
 
             // Add mask info
             os << '_';
             for (unsigned i=0; i<=lastMasked; ++i)
             {
                 if (mask_[i])
                     os << '1';
                 else
                     os << '0';
             }
         }
         return os.str();
     }
 
     double OrthoPolyGOFTest1D::inverseExceedance(
         const double pvalue, const unsigned long sz,
         const double smin, const double smax) const
     {
         if (pvalue <= 0.0 || pvalue >= 1.0) throw std::invalid_argument(
             "In npstat::OrthoPolyGOFTest1D::inverseExceedance: "
             "p-value argument must be inside (0, 1) interval");
 
         if (pvalue > 1.0e-10)
         {
             Gamma1D chisq(0.0, 2.0, nUnmasked_/2.0);
             return chisq.quantile(1.0 - pvalue);
         }
         else
             return AbsUnbinnedGOFTest1D::inverseExceedance(
                 pvalue, sz, smin, smax);
     }
 
     double OrthoPolyGOFTest1D::eigenUncertainty(const unsigned i) const
     {
         if (i >= invsqr_.nRows()) return 0.0;
         const long double normsq = sumOfSquares(invsqr_[i], invsqr_.nColumns());
         assert(normsq > 0.0L);
         return 1.0L/sqrtl(normsq);
     }
 
     std::vector<double> OrthoPolyGOFTest1D::covmatEigenvector(const unsigned i) const
     {
         if (i >= invsqr_.nRows()) throw std::invalid_argument(
             "In npstat::OrthoPolyGOFTest1D::covmatEigenvector: "
             "eigenvector number is out of range");
         const unsigned ncol = invsqr_.nColumns();
         std::vector<double> vec(invsqr_[i], invsqr_[i]+ncol);
         const double sigma = eigenUncertainty(i);
         for (unsigned i=0; i<ncol; ++i)
             vec[i] *= sigma;
         return vec;
     }
+
+    std::vector<OrthoPolyGOFTest1D::EigenInfo> OrthoPolyGOFTest1D::eigenComponents(
+        const unsigned i, const double threshold) const
+    {
+        const std::vector<double>& vec = covmatEigenvector(i);
+        std::vector<EigenInfo> result;
+
+        unsigned iRow = 0;
+        for (unsigned i=0; i<=maxdeg_; ++i)
+            for (unsigned k=i ? i : 1U; k<=maxdeg_; ++k, ++iRow)
+                if (std::abs(vec[iRow]) > threshold)
+                    result.push_back(EigenInfo(vec[iRow], UUPair(i, k)));
+        assert(iRow == vec.size());
+        return result;
+    }
 }
Index: trunk/npstat/stat/OrthoPolyGOFTest1D.hh
===================================================================
--- trunk/npstat/stat/OrthoPolyGOFTest1D.hh	(revision 748)
+++ trunk/npstat/stat/OrthoPolyGOFTest1D.hh	(revision 749)
@@ -1,110 +1,124 @@
 #ifndef NPSTAT_ORTHOPOLYGOFTEST1D_HH_
 #define NPSTAT_ORTHOPOLYGOFTEST1D_HH_
 
 /*!
 // \file OrthoPolyGOFTest1D.hh
 //
 // \brief Orthogonal polynomial goodness-of-fit test
 //
 // Author: I. Volobouev
 //
 // November 2020
 */
 
 #include <memory>
 
 #include "npstat/nm/Matrix.hh"
 #include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
+#include "npstat/nm/matrixIndexPairs.hh"
 
 #include "npstat/stat/AbsUnbinnedGOFTest1D.hh"
 
 namespace npstat {
     class OrthoPolyGOFTest1D : public AbsUnbinnedGOFTest1D
     {
     public:
+        typedef std::pair<double,UUPair> EigenInfo;
+
         // The mask, if provided, should contain 2*maxdeg elements.
         // Puts some number in the mask in order to discard the
         // corresponding eigenvector and put 0 there to keep it.
         template <class Quadrature>
         OrthoPolyGOFTest1D(const AbsDistribution1D& distro,
                            const Quadrature& q,
                            const AbsClassicalOrthoPoly1D& poly,
                            unsigned maxdeg,
                            const unsigned char* mask=0, unsigned lenMask=0);
 
         inline virtual ~OrthoPolyGOFTest1D() {}
 
         virtual std::string shortName() const;
 
         inline unsigned maxDegree() const {return maxdeg_;}
         inline unsigned numDeviations() const {return nUnmasked_;}
         inline bool isEigenvectorUsed(const unsigned i) const
             {return i >= 2U*maxdeg_ ? false : !mask_[i];}
 
         /**
-        // Get the uncertainty (i.e., the square root of the covariance matrix
-        // eigenvalue) associated with a particular eigenvector. This function
-        // ignores the mask. The return value is 0 if the eigenvector number is
-        // out of range.
+        // Get the uncertainty (i.e., the square root of the covariance
+        // matrix eigenvalue) associated with a particular eigenvector.
+        // This function ignores the mask. The return value is 0 if the
+        // eigenvector number is out of range.
         */
         double eigenUncertainty(unsigned i) const;
 
         /**
         // Get the covariance matrix eigenvector with the given number. If
         // the vector number is out of range, an exception will be thrown.
         // Interpretation of these eigenvectors is for experts only. Make
-        // sure you know what you are doing.
+        // sure that you know what you are doing.
         */
         std::vector<double> covmatEigenvector(unsigned i) const;
 
         /**
+        // Collect the information about the given eigenvector.
+        // The first element in the "EigenInfo" pair is the
+        // eigenvector value and the second element is the
+        // pair of polynomials corresponding to that value.
+        // Components are included in the result only if the
+        // magnitude of the eigenvector value is above the threshold.
+        */
+        std::vector<EigenInfo> eigenComponents(
+            unsigned i, double threshold) const;
+
+        /**
         // Get the value of the test statistic with all
         // terms weighted equally
         */
         inline virtual double testStatistic(
             const double* data, const unsigned long sz, bool) const
             {return testStat(data, sz);}
 
         inline virtual double testStatistic(
             const float* data, const unsigned long sz, bool) const
             {return testStat(data, sz);}
 
         inline virtual bool hasAnalyticPValue() const {return true;}
 
         virtual double analyticPValue(double stat, unsigned long sz) const;
 
         virtual double inverseExceedance(double pvalue, unsigned long sz,
                                          double smin, double smax) const;
 
         /**
         // Get individual normalized deviations from which
         // the test statistic is constructed. Only unmasked
         // deviations are returned. The return value of this
         // function is the number of such deviations.
         */
         template<typename Numeric>
         unsigned normalizedDeviations(
             const Numeric* data, unsigned long lenData,
             double* deviations, unsigned lenDeviations) const;
 
     private:
         template <class Quadrature>
         void calculateNormalizingMatrix(const Quadrature& q,
                                         Matrix<lapack_double>* m) const;
         template<typename Numeric>
         double testStat(const Numeric* data, unsigned long lenData) const;
 
         std::shared_ptr<const AbsClassicalOrthoPoly1D> poly_;
         Matrix<lapack_double> invsqr_;
         std::vector<double> buf_;
         std::vector<double> statBuf_;
         std::vector<double> sqrBuf_;
         std::vector<unsigned char> mask_;
         unsigned maxdeg_;
         unsigned nUnmasked_;
     };
 }
 
 #include "npstat/stat/OrthoPolyGOFTest1D.icc"
 
 #endif // NPSTAT_ORTHOPOLYGOFTEST1D_HH_