LSSTApplications  19.0.0+10,19.0.0+80,19.0.0-1-g20d9b18+31,19.0.0-1-g49a97f9+4,19.0.0-1-g8c57eb9+31,19.0.0-1-g9a028c0+13,19.0.0-1-ga72da6b+3,19.0.0-1-gb77924a+15,19.0.0-1-gbfe0924+66,19.0.0-1-gc3516c3,19.0.0-1-gefe1d0d+49,19.0.0-1-gf8cb8b4+3,19.0.0-14-g7511ce4+6,19.0.0-16-g3dc1a33c+6,19.0.0-17-g59f0e8a+4,19.0.0-17-g9c22e3c+9,19.0.0-18-g35bb99870+2,19.0.0-19-g2772d4a+9,19.0.0-2-g260436e+53,19.0.0-2-g31cdcee,19.0.0-2-g9675b69+3,19.0.0-2-gaa2795f,19.0.0-2-gbcc4de1,19.0.0-2-gd6f004e+6,19.0.0-2-gde8e5e3+5,19.0.0-2-gff6972b+15,19.0.0-21-ge306cd8,19.0.0-21-gf122e69+2,19.0.0-3-g2d43a51+2,19.0.0-3-gf3b1435+6,19.0.0-4-g56feb96,19.0.0-4-gac56cce+17,19.0.0-49-gce872c1+1,19.0.0-5-g66946eb+6,19.0.0-5-gd8897ba+6,19.0.0-51-gfc4a647e,19.0.0-7-g686a884+5,19.0.0-7-g886f805+5,19.0.0-8-g305ff64,w.2020.17
LSSTDataManagementBasePackage
ChebyshevBoundedField.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008-2014 LSST Corporation.
5  *
6  * This product includes software developed by the
7  * LSST Project (http://www.lsst.org/).
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the LSST License Statement and
20  * the GNU General Public License along with this program. If not,
21  * see <http://www.lsstcorp.org/LegalNotices/>.
22  */
23 
24 #include <memory>
25 
26 #include "ndarray/eigen.h"
35 
36 namespace lsst {
37 namespace afw {
38 
39 template std::shared_ptr<math::ChebyshevBoundedField> table::io::PersistableFacade<
40  math::ChebyshevBoundedField>::dynamicCast(std::shared_ptr<table::io::Persistable> const&);
41 
42 namespace math {
43 
45 
46 // ------------------ Constructors and helpers ---------------------------------------------------------------
47 
48 namespace {
49 
50 // Compute an affine transform that maps an arbitrary box to [-1,1]x[-1,1]
51 lsst::geom::AffineTransform makeChebyshevRangeTransform(lsst::geom::Box2D const bbox) {
54  lsst::geom::Extent2D(-(2.0 * bbox.getCenterX()) / bbox.getWidth(),
55  -(2.0 * bbox.getCenterY()) / bbox.getHeight()));
56 }
57 
58 } // namespace
59 
61  ndarray::Array<double const, 2, 2> const& coefficients)
62  : BoundedField(bbox),
63  _toChebyshevRange(makeChebyshevRangeTransform(lsst::geom::Box2D(bbox))),
64  _coefficients(coefficients) {}
65 
67  : BoundedField(bbox), _toChebyshevRange(makeChebyshevRangeTransform(lsst::geom::Box2D(bbox))) {}
68 
72 
73 // ------------------ fit() and helpers ---------------------------------------------------------------------
74 
75 namespace {
76 
78 typedef detail::TrapezoidalPacker Packer;
79 
80 // fill an array with 1-d Chebyshev functions of the 1st kind T(x), evaluated at the given point x
81 void evaluateBasis1d(ndarray::Array<double, 1, 1> const& t, double x) {
82  int const n = t.getSize<0>();
83  if (n > 0) {
84  t[0] = 1.0;
85  }
86  if (n > 1) {
87  t[1] = x;
88  }
89  for (int i = 2; i < n; ++i) {
90  t[i] = 2.0 * x * t[i - 1] - t[i - 2];
91  }
92 }
93 
94 // Create a matrix of 2-d Chebyshev functions evaluated at a set of positions, with
95 // Chebyshev order along columns and evaluation positions along rows. We pack the
96 // 2-d functions using the TrapezoidalPacker class, because we don't want any columns
97 // that correspond to coefficients that should be set to zero.
98 ndarray::Array<double, 2, 2> makeMatrix(ndarray::Array<double const, 1> const& x,
99  ndarray::Array<double const, 1> const& y,
100  lsst::geom::AffineTransform const& toChebyshevRange,
101  Packer const& packer, Control const& ctrl) {
102  int const nPoints = x.getSize<0>();
103  ndarray::Array<double, 1, 1> tx = ndarray::allocate(packer.nx);
104  ndarray::Array<double, 1, 1> ty = ndarray::allocate(packer.ny);
105  ndarray::Array<double, 2, 2> out = ndarray::allocate(nPoints, packer.size);
106  // Loop over x and y together, computing T_i(x) and T_j(y) arrays for each point,
107  // then packing them together.
108  for (int p = 0; p < nPoints; ++p) {
109  lsst::geom::Point2D sxy = toChebyshevRange(lsst::geom::Point2D(x[p], y[p]));
110  evaluateBasis1d(tx, sxy.getX());
111  evaluateBasis1d(ty, sxy.getY());
112  packer.pack(out[p], tx, ty); // this sets a row of out to the packed outer product of tx and ty
113  }
114  return out;
115 }
116 
117 // Create a matrix of 2-d Chebyshev functions evaluated on a grid of positions, with
118 // Chebyshev order along columns and evaluation positions along rows. We pack the
119 // 2-d functions using the TrapezoidalPacker class, because we don't want any columns
120 // that correspond to coefficients that should be set to zero.
121 ndarray::Array<double, 2, 2> makeMatrix(lsst::geom::Box2I const& bbox,
122  lsst::geom::AffineTransform const& toChebyshevRange,
123  Packer const& packer, Control const& ctrl) {
124  // Create a 2-d array that contains T_j(x) for each x value, with x values in rows and j in columns
125  ndarray::Array<double, 2, 2> tx = ndarray::allocate(bbox.getWidth(), packer.nx);
126  for (int x = bbox.getBeginX(), p = 0; p < bbox.getWidth(); ++p, ++x) {
127  evaluateBasis1d(tx[p], toChebyshevRange[lsst::geom::AffineTransform::XX] * x +
128  toChebyshevRange[lsst::geom::AffineTransform::X]);
129  }
130 
131  // Loop over y values, and at each point, compute T_i(y), then loop over x and multiply by the T_j(x)
132  // we already computed and stored above.
133  ndarray::Array<double, 2, 2> out = ndarray::allocate(bbox.getArea(), packer.size);
134  ndarray::Array<double, 2, 2>::Iterator outIter = out.begin();
135  ndarray::Array<double, 1, 1> ty = ndarray::allocate(ctrl.orderY + 1);
136  for (int y = bbox.getBeginY(), i = 0; i < bbox.getHeight(); ++i, ++y) {
137  evaluateBasis1d(ty, toChebyshevRange[lsst::geom::AffineTransform::YY] * y +
138  toChebyshevRange[lsst::geom::AffineTransform::Y]);
139  for (int j = 0; j < bbox.getWidth(); ++j, ++outIter) {
140  // this sets a row of out to the packed outer product of tx and ty
141  packer.pack(*outIter, tx[j], ty);
142  }
143  }
144  return out;
145 }
146 
147 } // namespace
148 
150  ndarray::Array<double const, 1> const& x,
151  ndarray::Array<double const, 1> const& y,
152  ndarray::Array<double const, 1> const& z,
153  Control const& ctrl) {
154  // Initialize the result object, so we can make use of the AffineTransform it builds
156  // This packer object knows how to map the 2-d Chebyshev functions onto a 1-d array,
157  // using only those that the control says should have nonzero coefficients.
158  Packer const packer(ctrl);
159  // Create a "design matrix" for the linear least squares problem (A in min||Ax-b||)
160  ndarray::Array<double, 2, 2> matrix = makeMatrix(x, y, result->_toChebyshevRange, packer, ctrl);
161  // Solve the linear least squares problem.
163  // Unpack the solution into a 2-d matrix, with zeros for values we didn't fit.
164  result->_coefficients = packer.unpack(lstsq.getSolution());
165  return result;
166 }
167 
169  ndarray::Array<double const, 1> const& x,
170  ndarray::Array<double const, 1> const& y,
171  ndarray::Array<double const, 1> const& z,
172  ndarray::Array<double const, 1> const& w,
173  Control const& ctrl) {
174  // Initialize the result object, so we can make use of the AffineTransform it builds
176  // This packer object knows how to map the 2-d Chebyshev functions onto a 1-d array,
177  // using only those that the control says should have nonzero coefficients.
178  Packer const packer(ctrl);
179  // Create a "design matrix" for the linear least squares problem ('A' in min||Ax-b||)
180  ndarray::Array<double, 2, 2> matrix = makeMatrix(x, y, result->_toChebyshevRange, packer, ctrl);
181  // We want to do weighted least squares, so we multiply both the data vector 'b' and the
182  // matrix 'A' by the weights.
183  ndarray::asEigenArray(matrix).colwise() *= ndarray::asEigenArray(w);
184  ndarray::Array<double, 1, 1> wz = ndarray::copy(z);
185  ndarray::asEigenArray(wz) *= ndarray::asEigenArray(w);
186  // Solve the linear least squares problem.
188  // Unpack the solution into a 2-d matrix, with zeros for values we didn't fit.
189  result->_coefficients = packer.unpack(lstsq.getSolution());
190  return result;
191 }
192 
193 template <typename T>
195  Control const& ctrl) {
196  // Initialize the result object, so we can make use of the AffineTransform it builds
199  // This packer object knows how to map the 2-d Chebyshev functions onto a 1-d array,
200  // using only those that the control says should have nonzero coefficients.
201  Packer const packer(ctrl);
202  ndarray::Array<double, 2, 2> matrix = makeMatrix(bbox, result->_toChebyshevRange, packer, ctrl);
203  // Flatten the data image into a 1-d vector.
204  ndarray::Array<double, 2, 2> imgCopy = ndarray::allocate(img.getArray().getShape());
205  imgCopy.deep() = img.getArray();
206  ndarray::Array<double const, 1, 1> z = ndarray::flatten<1>(imgCopy);
207  // Solve the linear least squares problem.
209  // Unpack the solution into a 2-d matrix, with zeros for values we didn't fit.
210  result->_coefficients = packer.unpack(lstsq.getSolution());
211  return result;
212 }
213 
214 // ------------------ modifier factories ---------------------------------------------------------------
215 
217  if (static_cast<std::size_t>(ctrl.orderX) >= _coefficients.getSize<1>()) {
219  (boost::format("New x order (%d) exceeds old x order (%d)") % ctrl.orderX %
220  (_coefficients.getSize<1>() - 1))
221  .str());
222  }
223  if (static_cast<std::size_t>(ctrl.orderY) >= _coefficients.getSize<0>()) {
225  (boost::format("New y order (%d) exceeds old y order (%d)") % ctrl.orderY %
226  (_coefficients.getSize<0>() - 1))
227  .str());
228  }
229  ndarray::Array<double, 2, 2> coefficients = ndarray::allocate(ctrl.orderY + 1, ctrl.orderX + 1);
230  coefficients.deep() = _coefficients[ndarray::view(0, ctrl.orderY + 1)(0, ctrl.orderX + 1)];
231  if (ctrl.triangular) {
232  Packer packer(ctrl);
233  ndarray::Array<double, 1, 1> packed = ndarray::allocate(packer.size);
234  packer.pack(packed, coefficients);
235  packer.unpack(coefficients, packed);
236  }
237  return std::make_shared<ChebyshevBoundedField>(getBBox(), coefficients);
238 }
239 
241  return std::make_shared<ChebyshevBoundedField>(bbox, _coefficients);
242 }
243 
244 // ------------------ evaluate() and helpers ---------------------------------------------------------------
245 
246 namespace {
247 
248 // To evaluate a 1-d Chebyshev function without needing to have workspace, we use the
249 // Clenshaw algorith, which is like going through the recurrence relation in reverse.
250 // The CoeffGetter argument g is something that behaves like an array, providing access
251 // to the coefficients.
252 template <typename CoeffGetter>
253 double evaluateFunction1d(CoeffGetter g, double x, int size) {
254  double b_kp2 = 0.0, b_kp1 = 0.0;
255  for (int k = (size - 1); k > 0; --k) {
256  double b_k = g[k] + 2 * x * b_kp1 - b_kp2;
257  b_kp2 = b_kp1;
258  b_kp1 = b_k;
259  }
260  return g[0] + x * b_kp1 - b_kp2;
261 }
262 
263 // This class imitates a 1-d array, by running evaluateFunction1d on a nested dimension;
264 // this lets us reuse the logic in evaluateFunction1d for both dimensions. Essentially,
265 // we run evaluateFunction1d on a column of coefficients to evaluate T_i(x), then pass
266 // the result of that to evaluateFunction1d with the results as the "coefficients" associated
267 // with the T_j(y) functions.
268 struct RecursionArrayImitator {
269  double operator[](int i) const {
270  return evaluateFunction1d(coefficients[i], x, coefficients.getSize<1>());
271  }
272 
273  RecursionArrayImitator(ndarray::Array<double const, 2, 2> const& coefficients_, double x_)
274  : coefficients(coefficients_), x(x_) {}
275 
276  ndarray::Array<double const, 2, 2> coefficients;
277  double x;
278 };
279 
280 } // namespace
281 
283  lsst::geom::Point2D p = _toChebyshevRange(position);
284  return evaluateFunction1d(RecursionArrayImitator(_coefficients, p.getX()), p.getY(),
285  _coefficients.getSize<0>());
286 }
287 
288 // The integral of T_n(x) over [-1,1]:
289 // https://en.wikipedia.org/wiki/Chebyshev_polynomials#Differentiation_and_integration
290 double integrateTn(int n) {
291  if (n % 2 == 1)
292  return 0;
293  else
294  return 2.0 / (1.0 - double(n * n));
295 }
296 
298  double result = 0;
299  double determinant = getBBox().getArea() / 4.0;
300  for (ndarray::Size j = 0; j < _coefficients.getSize<0>(); j++) {
301  for (ndarray::Size i = 0; i < _coefficients.getSize<1>(); i++) {
302  result += _coefficients[j][i] * integrateTn(i) * integrateTn(j);
303  }
304  }
305  return result * determinant;
306 }
307 
308 double ChebyshevBoundedField::mean() const { return integrate() / getBBox().getArea(); }
309 
310 // ------------------ persistence ---------------------------------------------------------------------------
311 
312 namespace {
313 
314 struct PersistenceHelper {
315  table::Schema schema;
316  table::Key<int> orderX;
318  table::Key<table::Array<double> > coefficients;
319 
320  PersistenceHelper(int nx, int ny)
321  : schema(),
322  orderX(schema.addField<int>("order_x", "maximum Chebyshev function order in x")),
323  bbox(table::Box2IKey::addFields(schema, "bbox", "bounding box", "pixel")),
324  coefficients(schema.addField<table::Array<double> >(
325  "coefficients", "Chebyshev function coefficients, ordered by y then x", nx * ny)) {}
326 
327  PersistenceHelper(table::Schema const& s)
328  : schema(s), orderX(s["order_x"]), bbox(s["bbox"]), coefficients(s["coefficients"]) {}
329 };
330 
331 class ChebyshevBoundedFieldFactory : public table::io::PersistableFactory {
332 public:
333  explicit ChebyshevBoundedFieldFactory(std::string const& name)
334  : afw::table::io::PersistableFactory(name) {}
335 
336  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
337  CatalogVector const& catalogs) const override {
338  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
339  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
340  table::BaseRecord const& record = catalogs.front().front();
341  PersistenceHelper const keys(record.getSchema());
342  lsst::geom::Box2I bbox(record.get(keys.bbox));
343  int nx = record.get(keys.orderX) + 1;
344  int ny = keys.coefficients.getSize() / nx;
345  LSST_ARCHIVE_ASSERT(nx * ny == keys.coefficients.getSize());
346  ndarray::Array<double, 2, 2> coefficients = ndarray::allocate(ny, nx);
347  ndarray::flatten<1>(coefficients) = record.get(keys.coefficients);
348  return std::make_shared<ChebyshevBoundedField>(bbox, coefficients);
349  }
350 };
351 
352 std::string getChebyshevBoundedFieldPersistenceName() { return "ChebyshevBoundedField"; }
353 
354 ChebyshevBoundedFieldFactory registration(getChebyshevBoundedFieldPersistenceName());
355 
356 } // namespace
357 
359  return getChebyshevBoundedFieldPersistenceName();
360 }
361 
362 std::string ChebyshevBoundedField::getPythonModule() const { return "lsst.afw.math"; }
363 
365  PersistenceHelper const keys(_coefficients.getSize<1>(), _coefficients.getSize<0>());
366  table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
367  std::shared_ptr<table::BaseRecord> record = catalog.addNew();
368  record->set(keys.orderX, _coefficients.getSize<1>() - 1);
369  record->set(keys.bbox, getBBox());
370  (*record)[keys.coefficients].deep() = ndarray::flatten<1>(_coefficients);
371  handle.saveCatalog(catalog);
372 }
373 
374 // ------------------ operators -----------------------------------------------------------------------------
375 
377  return std::make_shared<ChebyshevBoundedField>(getBBox(), ndarray::copy(getCoefficients() * scale));
378 }
379 
381  auto rhsCasted = dynamic_cast<ChebyshevBoundedField const*>(&rhs);
382  if (!rhsCasted) return false;
383 
384  return (getBBox() == rhsCasted->getBBox()) &&
385  (_coefficients.getShape() == rhsCasted->_coefficients.getShape()) &&
386  all(equal(_coefficients, rhsCasted->_coefficients));
387 }
388 
389 std::string ChebyshevBoundedField::toString() const {
391  os << "ChebyshevBoundedField (" << _coefficients.getShape() << " coefficients in y,x)";
392  return os.str();
393 }
394 
395 // ------------------ explicit instantiation ----------------------------------------------------------------
396 
397 #ifndef DOXYGEN
398 
399 #define INSTANTIATE(T) \
400  template std::shared_ptr<ChebyshevBoundedField> ChebyshevBoundedField::fit(image::Image<T> const& image, \
401  Control const& ctrl)
402 
403 INSTANTIATE(float);
404 INSTANTIATE(double);
405 
406 #endif
407 } // namespace math
408 } // namespace afw
409 } // namespace lsst
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
double getHeight() const noexcept
1-d interval accessors
Definition: Box.h:530
ChebyshevBoundedField(lsst::geom::Box2I const &bbox, ndarray::Array< double const, 2, 2 > const &coefficients)
Initialize the field from its bounding box an coefficients.
double getWidth() const noexcept
1-d interval accessors
Definition: Box.h:529
int getHeight() const noexcept
Definition: Box.h:188
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
An object passed to Persistable::write to allow it to persist itself.
An affine coordinate transformation consisting of a linear transformation and an offset.
BoxKey< lsst::geom::Box2I > Box2IKey
Definition: aggregates.h:201
double mean() const override
Compute the mean of this function over its bounding-box.
table::Box2IKey bbox
Reports attempts to exceed implementation-defined length limits for some classes. ...
Definition: Runtime.h:76
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:109
A control object used when fitting ChebyshevBoundedField to data (see ChebyshevBoundedField::fit) ...
int y
Definition: SpanSet.cc:49
Solver for linear least-squares problems.
Definition: LeastSquares.h:67
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
static LeastSquares fromDesignMatrix(ndarray::Array< T1, 2, C1 > const &design, ndarray::Array< T2, 1, C2 > const &data, Factorization factorization=NORMAL_EIGENSYSTEM)
Initialize from the design matrix and data vector given as ndarrays.
Definition: LeastSquares.h:100
double z
Definition: Match.cc:44
std::shared_ptr< ChebyshevBoundedField > truncate(Control const &ctrl) const
Return a new ChebyshevBoudedField with maximum orders set by the given control object.
double evaluate(lsst::geom::Point2D const &position) const override
Evaluate the field at the given point.
int getBeginY() const noexcept
Definition: Box.h:173
STL class.
FastFinder::Iterator Iterator
Definition: FastFinder.cc:178
bool operator==(BoundedField const &rhs) const override
BoundedFields (of the same sublcass) are equal if their bounding boxes and parameters are equal...
int computeSize() const
Return the number of nonzero coefficients in the Chebyshev function defined by this object...
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
BoundedField(BoundedField const &)=default
Use the normal equations with a symmetric Eigensystem decomposition.
Definition: LeastSquares.h:72
int orderY
"maximum Chebyshev function order in y" ;
double getCenterX() const noexcept
Definition: Box.h:552
A base class for image defects.
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition: ImageBase.h:482
int getArea() const
Definition: Box.h:189
int getBeginX() const noexcept
Definition: Box.h:172
lsst::geom::Box2I getBBox() const
Return the bounding box that defines the region where the field is valid.
Definition: BoundedField.h:112
T str(T... args)
int getWidth() const noexcept
Definition: Box.h:187
std::shared_ptr< ChebyshevBoundedField > relocate(lsst::geom::Box2I const &bbox) const
Return a new ChebyshevBoundedField with domain set to the given bounding box.
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
static LinearTransform makeScaling(double s) noexcept
double x
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
bool triangular
"if true, only include terms where the sum of the x and y order " "is less than or equal to max(order...
double w
Definition: CoaddPsf.cc:69
double getCenterY() const noexcept
Definition: Box.h:553
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
double integrate() const override
Compute the integral of this function over its bounding-box.
ndarray::Array< double const, 2, 2 > getCoefficients() const
Return the coefficient matrix.
A BoundedField based on 2-d Chebyshev polynomials of the first kind.
An abstract base class for 2-d functions defined on an integer bounding boxes.
Definition: BoundedField.h:55
ndarray::Array< double const, 1, 1 > getSolution()
Return the vector solution to the least squares problem.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
std::shared_ptr< BoundedField > operator*(double const scale) const override
Return a scaled BoundedField.
std::ostream * os
Definition: Schema.cc:746
An integer coordinate rectangle.
Definition: Box.h:55
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
py::object result
Definition: _schema.cc:429
table::Schema schema
A helper class ChebyshevBoundedField, for mapping trapezoidal matrices to 1-d arrays.
table::Key< int > orderX
#define INSTANTIATE(FROMSYS, TOSYS)
Definition: Detector.cc:484
ndarray::Array< double const, 2, 2 > coefficients
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:485
static std::shared_ptr< ChebyshevBoundedField > fit(lsst::geom::Box2I const &bbox, ndarray::Array< double const, 1 > const &x, ndarray::Array< double const, 1 > const &y, ndarray::Array< double const, 1 > const &z, Control const &ctrl)
Fit a Chebyshev approximation to non-gridded data with equal weights.
int orderX
"maximum Chebyshev function order in x" ;