Index: trunk/tests/test_ScalableClassicalOrthoPoly1D.cc
===================================================================
--- trunk/tests/test_ScalableClassicalOrthoPoly1D.cc	(revision 884)
+++ trunk/tests/test_ScalableClassicalOrthoPoly1D.cc	(revision 885)
@@ -1,139 +1,150 @@
 #include <cmath>
 
 #include "UnitTest++.h"
 #include "test_utils.hh"
 
 #include "npstat/nm/ScalableClassicalOrthoPoly1D.hh"
 #include "npstat/nm/ClassicalOrthoPolys1D.hh"
 #include "npstat/nm/GaussLegendreQuadrature.hh"
+#include "npstat/nm/opsRootsFromJacobiMatrix.hh"
 
 #include "npstat/rng/MersenneTwister.hh"
 
 using namespace npstat;
 using namespace std;
 
 namespace {
     struct SomeFcn
     {
         inline double operator()(const double x) const
             {return std::pow(std::abs(x), 2.5);}
     };
 
     TEST(ScalableClassicalOrthoPoly1D_ref)
     {
         const double eps = 1.0e-14;
 
         LegendreOrthoPoly1D p0;
         ShiftedLegendreOrthoPoly1D ref;
         ScalableClassicalOrthoPoly1D spoly(p0, 0.5, 0.5);
 
         CHECK_EQUAL(0.5, spoly.location());
         CHECK_EQUAL(0.5, spoly.scale());
 
         CHECK_EQUAL(ref.xmin(), spoly.xmin());
         CHECK_EQUAL(ref.xmax(), spoly.xmax());
         CHECK_EQUAL(ref.maxDegree(), spoly.maxDegree());
 
         const unsigned maxdeg = 10;
         const unsigned ntries = 100;
         const unsigned npoints = 64;
 
+        for (unsigned ideg=1; ideg<=maxdeg; ++ideg)
+        {
+            const std::vector<double>& roots1 = opsRootsFromJacobiMatrix(ref, ideg);
+            const std::vector<double>& roots2 = opsRootsFromJacobiMatrix(spoly, ideg);
+            CHECK_EQUAL(ideg, roots1.size());
+            CHECK_EQUAL(ideg, roots2.size());
+            for (unsigned i=0; i<ideg; ++i)
+                CHECK_CLOSE(roots1[i], roots2[i], eps);
+        }
+
         std::vector<long double> valuesRef(maxdeg+1U);
         std::vector<long double> valuesTest(maxdeg+1U);
         std::vector<double> sample(npoints);
         std::vector<std::pair<double,double> > wsample(npoints);
         std::vector<std::pair<double,double> > wsample1(npoints);
         std::vector<double> coeffs(maxdeg+1U);
         std::vector<double> coeffs1(maxdeg+1U);
         std::vector<double> coeffsRef(maxdeg+1U);
         std::vector<double> coeffVars(maxdeg+1U);
         std::vector<double> coeffVars1(maxdeg+1U);
         std::vector<double> coeffVarsRef(maxdeg+1U);
 
         MersenneTwister rng;
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             const double x = rng();
             const unsigned ideg = rng()*(maxdeg + 1U);
             const unsigned ideg2 = rng()*(maxdeg + 1U);
 
             CHECK_CLOSE(ref.weight(x), spoly.weight(x), eps);
             CHECK_CLOSE(ref.poly(ideg, x), spoly.poly(ideg, x), eps);
 
             std::pair<long double,long double> p1(ref.twopoly(ideg, ideg2, x));
             std::pair<long double,long double> p2(spoly.twopoly(ideg, ideg2, x));
             CHECK_CLOSE(p1.first, p2.first, eps);
             CHECK_CLOSE(p1.second, p2.second, eps);
 
             ref.allpoly(x, &valuesRef[0], maxdeg);
             spoly.allpoly(x, &valuesTest[0], maxdeg);
             for (unsigned i=0; i<=maxdeg; ++i)
                 CHECK_CLOSE(valuesRef[i], valuesTest[i], eps);
 
             for (unsigned i=0; i<=ideg; ++i)
                 coeffs.at(i) = rng();
 
             CHECK_CLOSE(ref.series(&coeffs[0], ideg, x),
                         spoly.series(&coeffs[0], ideg, x), eps);
             CHECK_CLOSE(ref.integrateSeries(&coeffs[0], ideg),
                         spoly.integrateSeries(&coeffs[0], ideg), eps);
 
             for (unsigned i=0; i<npoints; ++i)
             {
                 sample[i] = rng();
                 wsample[i].first = sample[i];
                 wsample[i].second = rng();
                 wsample1[i].first = sample[i];
                 wsample1[i].second = 1.0;
             }
 
             ref.sampleCoeffs(&sample[0], npoints, &coeffsRef[0], maxdeg);
             spoly.sampleCoeffs(&sample[0], npoints, &coeffs[0], maxdeg);
             spoly.weightedSampleCoeffs(&wsample1[0], npoints, &coeffs1[0], maxdeg);
             ref.sampleCoeffVars(&sample[0], npoints, &coeffsRef[0], maxdeg, &coeffVarsRef[0]);
             spoly.sampleCoeffVars(&sample[0], npoints, &coeffs[0], maxdeg, &coeffVars[0]);
             spoly.weightedSampleCoeffVars(&wsample1[0], npoints, &coeffs1[0], maxdeg, &coeffVars1[0]);
 
             for (unsigned i=0; i<=maxdeg; ++i)
             {
                 CHECK_CLOSE(coeffsRef[i], coeffs[i], eps);
                 CHECK_CLOSE(coeffsRef[i], coeffs1[i], eps);
                 CHECK_CLOSE(coeffVarsRef[i], coeffVars[i], eps);
                 CHECK_CLOSE(coeffVarsRef[i], coeffVars1[i], eps);
             }
 
             const npstat::Matrix<double>& covref = ref.sampleCoeffCovariance(
                 &sample[0], npoints, &coeffsRef[0], maxdeg);
             const npstat::Matrix<double>& cov = spoly.sampleCoeffCovariance(
                 &sample[0], npoints, &coeffs[0], maxdeg);
             const npstat::Matrix<double>& cov3 = spoly.weightedSampleCoeffCovariance(
                 &wsample1[0], npoints, &coeffs1[0], maxdeg);
             CHECK((covref - cov).maxAbsValue() < eps);
             CHECK((covref - cov3).maxAbsValue() < eps);
 
             ref.weightedSampleCoeffs(&wsample[0], npoints, &coeffsRef[0], maxdeg);
             spoly.weightedSampleCoeffs(&wsample[0], npoints, &coeffs[0], maxdeg);
             ref.weightedSampleCoeffVars(&wsample[0], npoints, &coeffsRef[0], maxdeg, &coeffVarsRef[0]);
             spoly.weightedSampleCoeffVars(&wsample[0], npoints, &coeffs[0], maxdeg, &coeffVars[0]);
 
             for (unsigned i=0; i<=maxdeg; ++i)
             {
                 CHECK_CLOSE(coeffsRef[i], coeffs[i], eps);
                 CHECK_CLOSE(coeffVarsRef[i], coeffVars[i], eps);
             }
 
             const npstat::Matrix<double>& covref2 = ref.weightedSampleCoeffCovariance(
                 &wsample[0], npoints, &coeffsRef[0], maxdeg);
             const npstat::Matrix<double>& cov2 = spoly.weightedSampleCoeffCovariance(
                 &wsample[0], npoints, &coeffs[0], maxdeg);
             CHECK((covref2 - cov2).maxAbsValue() < eps);
         }
 
         SomeFcn fcn;
         GaussLegendreQuadrature glq(128);
         ref.calculateCoeffs(fcn, glq, &coeffsRef[0], maxdeg);
         spoly.calculateCoeffs(fcn, glq, &coeffs[0], maxdeg);
         for (unsigned i=0; i<=maxdeg; ++i)
             CHECK_CLOSE(coeffsRef[i], coeffs[i], eps);
     }
 }
Index: trunk/npstat/nm/ScalableClassicalOrthoPoly1D.hh
===================================================================
--- trunk/npstat/nm/ScalableClassicalOrthoPoly1D.hh	(revision 884)
+++ trunk/npstat/nm/ScalableClassicalOrthoPoly1D.hh	(revision 885)
@@ -1,316 +1,319 @@
 #ifndef NPSTAT_SCALABLECLASSICALORTHOPOLY1D_HH_
 #define NPSTAT_SCALABLECLASSICALORTHOPOLY1D_HH_
 
 /*!
 // \file ScalableClassicalOrthoPoly1D.hh
 //
 // \brief Class for shifting and scaling classical continuous orthonormal
 //        polynomials
 //
 // Author: I. Volobouev
 //
 // October 2020
 */
 
 #include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
 
 namespace npstat {
     class ScalableOrthoPoly1DDeg;
     class ScalableOrthoPoly1DWeight;
 
     class ScalableClassicalOrthoPoly1D
     {
     public:
         typedef ScalableOrthoPoly1DDeg degree_functor;
         typedef ScalableOrthoPoly1DWeight weight_functor;
 
         ScalableClassicalOrthoPoly1D(const AbsClassicalOrthoPoly1D& poly,
                                      long double location, long double scale);
 
         ScalableClassicalOrthoPoly1D(const AbsClassicalOrthoPoly1D& poly,
                                      double location, double scale);
 
         ScalableClassicalOrthoPoly1D(const ScalableClassicalOrthoPoly1D&);
         ScalableClassicalOrthoPoly1D& operator=(const ScalableClassicalOrthoPoly1D&);
 
         inline ScalableClassicalOrthoPoly1D* clone() const
             {return new ScalableClassicalOrthoPoly1D(*this);}
 
         inline ~ScalableClassicalOrthoPoly1D() {delete poly_;}
 
         inline long double location() const {return location_;}
         inline long double scale() const {return scale_;}
+
+        // Note that the Jacobi matrix principal minor is returned
+        // for the unscaled OPS
         inline Matrix<long double> jacobiMatrix(const unsigned n) const
             {return poly_->jacobiMatrix(n);}
 
         inline void setLocation(const long double v) {location_ = v;}
         inline void setScale(const long double v)
         {
             if (v <= 0.0) throw std::invalid_argument(
                 "In npstat::ScalableClassicalOrthoPoly1D::setScale: "
                 "scale parameter must be positive");
             scale_ = v;
         }
 
         inline long double weight(const long double x) const
             {return poly_->weight((x - location_)/scale_)/scale_;}
 
         inline double xmin() const {return scale_*poly_->xmin() + location_;}
 
         inline double xmax() const {return scale_*poly_->xmax() + location_;}
 
         inline unsigned maxDegree() const {return poly_->maxDegree();}
 
         inline long double poly(const unsigned deg, const long double x) const
             {return poly_->poly(deg, (x - location_)/scale_);}
 
         /** 
         // Values of two orthonormal polynomials.
         // Faster than calling "poly" two times.
         */
         inline std::pair<long double,long double> twopoly(
             const unsigned deg1, const unsigned deg2, const long double x) const
             {return poly_->twopoly(deg1, deg2, (x - location_)/scale_);}
 
         /**
         // Values of all orthonormal polynomials up to some degree.
         // Faster than calling "poly" multiple times. The size of
         // the "values" array should be at least maxdeg + 1.
         */
         inline void allpoly(const long double x,
                             long double* values, const unsigned maxdeg) const
             {poly_->allpoly((x - location_)/scale_, values, maxdeg);}
 
         inline double series(const double* coeffs, const unsigned maxdeg,
                              const double x) const
             {return poly_->series(coeffs, maxdeg, (x - location_)/scale_);}
 
         inline double integrateSeries(const double* coeffs, const unsigned maxdeg) const
             {return scale_*poly_->integrateSeries(coeffs, maxdeg);}
 
         /**
         // Build the coefficients of the orthogonal polynomial series
         // for the given function. The length of the array "coeffs"
         // should be at least maxdeg + 1. Note that the coefficients
         // are returned in the order of increasing degree (same order
         // is used by the "series" function).
         */
         template <class Functor, class Quadrature>
         void calculateCoeffs(const Functor& fcn, const Quadrature& quad,
                              double* coeffs, unsigned maxdeg) const;
 
         /**
         // Build the coefficients of the orthogonal polynomial series
         // for the given sample of points (empirical density function).
         // The length of the array "coeffs" should be at least maxdeg + 1.
         // Note that the coefficients are returned in the order of
         // increasing degree (same order is used by the "series" function).
         */
         template <class Numeric>
         inline void sampleCoeffs(const Numeric* coords, const unsigned long lenCoords,
                                  double* coeffs, const unsigned maxdeg) const
         {
             const Numeric mu = static_cast<Numeric>(location_);
             const Numeric s = static_cast<Numeric>(scale_);
             poly_->sampleCoeffs(coords, lenCoords, mu, s, coeffs, maxdeg);
         }
 
         template <class Numeric>
         inline void sampleCoeffVars(const Numeric* coords, const unsigned long lenCoords,
                                     const double* coeffs, const unsigned maxdeg,
                                     double* variances) const
         {
             const Numeric mu = static_cast<Numeric>(location_);
             const Numeric s = static_cast<Numeric>(scale_);
             poly_->sampleCoeffVars(coords, lenCoords, mu, s, coeffs, maxdeg, variances);
         }
 
         template <class Numeric>
         inline npstat::Matrix<double>
         sampleCoeffCovariance(const Numeric* coords, const unsigned long lenCoords,
                               const double* coeffs, const unsigned maxdeg) const
         {
             const Numeric mu = static_cast<Numeric>(location_);
             const Numeric s = static_cast<Numeric>(scale_);
             return poly_->sampleCoeffCovariance(coords, lenCoords, mu, s, coeffs, maxdeg);
         }
 
         /**
         // Build the coefficients of the orthogonal polynomial series for
         // the given sample of weighted points (empirical density function).
         // The first element of the pair will be treated as the coordinate
         // and the second element will be treated as weight. Weights must
         // not be negative. The length of the array "coeffs" should be
         // at least maxdeg + 1. Note that the coefficients are returned
         // in the order of increasing degree (same order is used by the
         // "series" function).
         */
         template <class Pair>
         inline void weightedSampleCoeffs(const Pair* points, const unsigned long numPoints,
                                          double* coeffs, const unsigned maxdeg) const
         {
             typedef typename Pair::first_type Numeric;
             const Numeric mu = static_cast<Numeric>(location_);
             const Numeric s = static_cast<Numeric>(scale_);
             return poly_->weightedSampleCoeffs(points, numPoints, mu, s, coeffs, maxdeg);
         }
 
         template <class Pair>
         inline void weightedSampleCoeffVars(const Pair* points, const unsigned long nPoints,
                                             const double* coeffs, const unsigned maxdeg,
                                             double* variances) const
         {
             typedef typename Pair::first_type Numeric;
             const Numeric mu = static_cast<Numeric>(location_);
             const Numeric s = static_cast<Numeric>(scale_);
             return poly_->weightedSampleCoeffVars(points,nPoints,mu,s,coeffs,maxdeg,variances);
         }
 
         template <class Pair>
         inline npstat::Matrix<double>
         weightedSampleCoeffCovariance(const Pair* points, const unsigned long nPoints,
                                       const double* coeffs, const unsigned maxdeg) const
         {
             typedef typename Pair::first_type Numeric;
             const Numeric mu = static_cast<Numeric>(location_);
             const Numeric s = static_cast<Numeric>(scale_);
             return poly_->weightedSampleCoeffCovariance(points, nPoints, mu, s, coeffs, maxdeg);
         }
 
         /**
         // Unweighted averages of the polynomial values for the given sample.
         // The length of array "averages" should be at least maxdeg + 1.
         */
         template <class Numeric>
         void sampleAverages(const Numeric* coords, unsigned long lenCoords,
                             double* averages, unsigned maxdeg) const;
 
         /**
         // Averages of the polynomial values weighted by an external
         // weight function. The length of array "averages" should be
         // at least maxdeg + 1.
         */
         template <class Functor>
         void extWeightAverages(const Functor& extWeight,
                                const AbsIntervalQuadrature1D& quad,
                                double* averages, unsigned maxdeg) const;
 
         /**
         // Pairwise products of the polynomial values weighted by
         // an external weight function. The returned matrix will be
         // symmetric and will have the dimensions (maxdeg + 1) x (maxdeg + 1).
         */
         template <class Functor>
         Matrix<double> extWeightProductAverages(const Functor& extWeight,
             const AbsIntervalQuadrature1D& quad, unsigned maxdeg) const;
 
     private:
         AbsClassicalOrthoPoly1D* poly_;
         long double location_;
         long double scale_;
 
 #ifdef SWIG
     public:
         inline void setLocation2(const double v)
             {setLocation(v);}
 
         inline void setScale2(const double v)
             {setScale(v);}
 
         inline double location2() const
             {return location();}
 
         inline double scale2() const
             {return scale();}
 
         inline double weight2(const double x) const
             {return weight(x);}
 
         inline double poly2(const unsigned deg, const double x) const
             {return poly(deg, x);}
 
         inline std::pair<double,double> twopoly2(const unsigned deg1,
                                                  const unsigned deg2,
                                                  const double x) const
         {
             const std::pair<long double,long double>& p = twopoly(deg1, deg2, x);
             return std::pair<double,double>(p.first, p.second);
         }
 
         inline std::vector<double> allpoly2(const unsigned maxdeg, const double x) const
         {
             std::vector<long double> p(maxdeg+1);
             allpoly(x, &p[0], p.size());
             return std::vector<double>(p.begin(), p.end());
         }
 #endif // SWIG
     };
 
     /**
     // A functor for the weight function of the given ortho poly system.
     // The poly system is not copied, only a reference is used. It is
     // a responsibility of the user to make sure that the lifetime of
     // the poly system object exceeds the lifetime of the functor.
     */
     class ScalableOrthoPoly1DWeight : public Functor1<long double, long double>
     {
     public:
         inline explicit ScalableOrthoPoly1DWeight(
             const ScalableClassicalOrthoPoly1D& fcn,
             const long double normfactor)
             : fcn_(fcn), norm_(normfactor) {}
 
         // Separate constructor here because swig gets confused
         // by long doubles with defaults
         inline explicit ScalableOrthoPoly1DWeight(
             const ScalableClassicalOrthoPoly1D& fcn)
             : fcn_(fcn), norm_(1.0L) {}
 
         inline virtual ~ScalableOrthoPoly1DWeight() {}
 
         inline virtual long double operator()(const long double& a) const
             {return norm_*fcn_.weight(a);}
 
     private:
         ScalableOrthoPoly1DWeight();
 
         const ScalableClassicalOrthoPoly1D& fcn_;
         long double norm_;
     };
 
     /**
     // A functor for a particular degree polynomial of the given ortho
     // poly system. The poly system is not copied, only a reference is used.
     // It is a responsibility of the user to make sure that the lifetime of
     // the poly system object exceeds the lifetime of the functor.
     */
     class ScalableOrthoPoly1DDeg : public Functor1<long double, long double>
     {
     public:
         inline ScalableOrthoPoly1DDeg(const ScalableClassicalOrthoPoly1D& fcn,
                                       const unsigned degree,
                                       const long double normfactor=1.0L)
             : fcn_(fcn), norm_(normfactor), deg_(degree)
         {
             if (deg_ > fcn_.maxDegree()) throw std::invalid_argument(
                 "In npstat::ScalableOrthoPoly1DDeg::constructor: "
                 "degree argument is out of range");
         }
 
         inline virtual ~ScalableOrthoPoly1DDeg() {}
 
         inline virtual long double operator()(const long double& a) const
             {return norm_*fcn_.poly(deg_, a);}
 
     private:
         ScalableOrthoPoly1DDeg();
 
         const ScalableClassicalOrthoPoly1D& fcn_;
         long double norm_;
         unsigned deg_;
     };
 }
 
 #include "npstat/nm/ScalableClassicalOrthoPoly1D.icc"
 
 #endif // NPSTAT_SCALABLECLASSICALORTHOPOLY1D_HH_