LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Legendre2D.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 #include <algorithm>
26 #include <string>
27 #include <sstream>
28 #include <stdexcept>
29 
33 
34 namespace lsst {
35 namespace meas {
36 namespace algorithms {
37 namespace shapelet {
38 
39  static std::string makeRangeExceptionMessage(
40  const Position& p, const Bounds& b)
41  {
42  std::ostringstream s;
43  s << "Range error: Position "<<p<<" is not in Bounds "<<b;
44  return s.str();
45  }
46 
48  std::runtime_error(makeRangeExceptionMessage(p,b)), _p(p), _b(b)
49  {}
50 
52  int xOrder, int yOrder, const DVector& fVect)
53  {
54  if (_xOrder != xOrder || _yOrder != yOrder) {
55  _xOrder = xOrder; _yOrder = yOrder;
56  _coeffs.reset(new DMatrix(_xOrder+1,_yOrder+1));
57  _coeffs->setZero();
58  }
59  int k=0;
60  for(int m=0; m <= std::max(_xOrder,_yOrder); ++m) {
61  for(int i=std::min(m,_xOrder);m-i<=std::min(m,_yOrder);--i) {
62  (*_coeffs)(i,m-i) = fVect(k++);
63  }
64  }
65  Assert(k==(int)fVect.size());
66  }
67 
68  Legendre2D::Legendre2D(std::istream& fin) : Function2D()
69  {
70  // Order of parameters is same as for Polynomial2D. Difference
71  // is that instead of x,x^2,x^3,etc., we use P1(x),P2(x),P3(x),etc.
72  fin >> _xOrder >> _yOrder >> _bounds;
73  if (!fin) throw std::runtime_error("reading order, bounds");
74  _coeffs.reset(new DMatrix(_xOrder+1,_yOrder+1));
75  _coeffs->setZero();
76  int maxOrder = std::max(_xOrder,_yOrder);
77  for(int m=0; m<=maxOrder; ++m) {
78  for(int i=std::min(m,_xOrder); m-i<=std::min(m,_yOrder); --i) {
79  fin >> (*_coeffs)(i,m-i);
80  }
81  }
82  if (!fin) throw std::runtime_error("reading (legendre) function");
83  }
84 
85  void Legendre2D::write(std::ostream& fout) const
86  {
87  int oldprec = fout.precision(6);
88  std::ios::fmtflags oldf =
89  fout.setf(std::ios::scientific,std::ios::floatfield);
90  int maxOrder = std::max(_xOrder,_yOrder);
91  if (maxOrder == 0) {
92  fout << "C " << (*_coeffs)(0,0) << std::endl;
93  } else {
94  fout << "L " << _xOrder << ' ' << _yOrder << ' ' << _bounds << ' ';
95  for(int m=0; m<=maxOrder; ++m) {
96  for(int i=std::min(m,_xOrder);m-i<=std::min(m,_yOrder);--i) {
97  fout << (*_coeffs)(i,m-i) << ' ';
98  }
99  }
100  fout << std::endl;
101  }
102  if (!fout) throw std::runtime_error("writing (legendre) function");
103  fout.precision(oldprec);
104  fout.flags(oldf);
105  }
106 
107  void Legendre2D::addLinear(double a, double b, double c)
108  {
109  double xAve = (getXMin() + getXMax())/2.;
110  double yAve = (getYMin() + getYMax())/2.;
111  double xRange = getXMax() - getXMin();
112  double yRange = getYMax() - getYMin();
113 
114  (*_coeffs)(0,0) += a + b*xAve + c*yAve;
115  (*_coeffs)(1,0) += b*xRange/2.;
116  (*_coeffs)(0,1) += c*yRange/2.;
117  }
118 
120  double , double , double , double , double , double )
121  {
122  // Not implemented yet.
123  Assert(false);
124  }
125 
127  {
128  const Legendre2D* lrhs = dynamic_cast<const Legendre2D*>(&rhs);
129  Assert(lrhs);
130  Assert(getBounds() == lrhs->getBounds());
131  if (_xOrder == lrhs->_xOrder && _yOrder == lrhs->_yOrder) {
132  *_coeffs += *lrhs->_coeffs;
133  } else {
134  int newXOrder = std::max(_xOrder,lrhs->_xOrder);
135  int newYOrder = std::max(_yOrder,lrhs->_yOrder);
136  std::auto_ptr<DMatrix > newc(
137  new DMatrix(newXOrder+1,newYOrder+1));
138  newc->setZero();
139  newc->TMV_subMatrix(0,_xOrder+1,0,_yOrder+1) = *_coeffs;
140  newc->TMV_subMatrix(0,lrhs->_xOrder+1,0,lrhs->_yOrder+1) += *lrhs->_coeffs;
141  _coeffs = newc;
142  _xOrder = newXOrder;
143  _yOrder = newYOrder;
144  }
145  }
146 
147  // dP_2n(x)/dx = Sum_k=0..n-1 (4k+3) P_2k+1(x)
148  // dP_2n+1(x)/dx = Sum_k=0..n (4k+1) P_2k(x)
149  std::auto_ptr<Function2D > Legendre2D::dFdX() const
150  {
151  if (_xOrder == 0) {
152  return std::auto_ptr<Function2D >(new Constant2D());
153  }
154  if (_xOrder == 1 && _yOrder == 0) {
155  return std::auto_ptr<Function2D >(
156  new Constant2D((*_coeffs)(1,0)*2./(getXMax()-getXMin())));
157  }
158 
159  int newXOrder = _xOrder-1;
160  int newYOrder = _xOrder > _yOrder ? _yOrder : _yOrder-1;
161 
162  std::auto_ptr<Legendre2D > temp(
163  new Legendre2D(newXOrder,newYOrder,_bounds));
164  // initialized to 0's
165 
166  int maxOrder = std::max(_xOrder,_yOrder);
167  for(int i=_xOrder;i>=1;--i) {
168  for(int j=std::min(maxOrder-i,_yOrder);j>=0;--j) {
169  if (i%2 == 0) {
170  for(int k=0;k<=i/2-1;++k)
171  (*temp->_coeffs)(2*k+1,j) += (4.*k+3.)*(*_coeffs)(i,j);
172  } else {
173  for(int k=0;k<=(i-1)/2;++k)
174  (*temp->_coeffs)(2*k,j) += (4.*k+1.)*(*_coeffs)(i,j);
175  }
176  }
177  }
178  maxOrder = std::max(newXOrder,newYOrder);
179  for(int i=newXOrder;i>=0;--i) {
180  for(int j=std::min(maxOrder-i,newYOrder);j>=0;--j) {
181  (*temp->_coeffs)(i,j) *= 2./(getXMax()-getXMin());
182  }
183  }
184  return std::auto_ptr<Function2D >(temp);
185  }
186 
187  // dP_2n(x)/dx = Sum_k=0..n-1 (4k+3) P_2k+1(x)
188  // dP_2n+1(x)/dx = Sum_k=0..n-1 (4k+1) P_2k(x)
189  std::auto_ptr<Function2D > Legendre2D::dFdY() const
190  {
191  if (_yOrder == 0) {
192  return std::auto_ptr<Function2D >(new Constant2D());
193  }
194  if (_yOrder == 1 && _xOrder == 0) {
195  return std::auto_ptr<Function2D >(new Constant2D(
196  (*_coeffs)(0,1)*2./(getYMax()-getYMin())));
197  }
198 
199  int newXOrder = _yOrder > _xOrder ? _xOrder : _xOrder-1;
200  int newYOrder = _yOrder-1;
201 
202  std::auto_ptr<Legendre2D > temp(
203  new Legendre2D(newXOrder,newYOrder,_bounds));
204  // initialized to 0's
205 
206  int maxOrder = std::max(_xOrder,_yOrder);
207  for(int j=_yOrder;j>=1;--j) {
208  for(int i=std::min(maxOrder-j,_xOrder);i>=0;--i) {
209  if (j%2 == 0) {
210  for(int k=0;k<=(j-2)/2;++k)
211  (*temp->_coeffs)(i,2*k+1) += (4.*k+3.)*(*_coeffs)(i,j);
212  } else {
213  for(int k=0;k<=(j-1)/2;++k)
214  (*temp->_coeffs)(i,2*k) += (4.*k+1.)*(*_coeffs)(i,j);
215  }
216  }
217  }
218  maxOrder = std::max(newXOrder,newYOrder);
219  for(int j=newYOrder;j>=0;--j) {
220  for(int i=std::min(maxOrder-j,newXOrder);i>=0;--i) {
221  (*temp->_coeffs)(i,j) *= 2./(getYMax()-getYMin());
222  }
223  }
224  return std::auto_ptr<Function2D >(temp);
225  }
226 
228  int order, double xy, double min, double max) const
229  {
230  DVector temp(order+1);
231  double xyNew = (2.*xy-min-max)/(max-min);
232  temp[0] = 1.; if(order>0) temp[1] = xyNew;
233  for(int i=2;i<=order;++i)
234  temp[i] = ((2.*i-1.)*xyNew*temp[i-1] - (i-1.)*temp[i-2])/i;
235  return temp;
236  }
237 
238  DVector Legendre2D::definePX(int order, double x) const
239  {
240  if (x < getXMin() || x > getXMax()) {
241  dbg<<"Error: x = "<<x<<
242  ", min,max = "<<getXMin()<<','<<getXMax()<<std::endl;
244  }
245  return definePXY(order,x,getXMin(),getXMax());
246  }
247 
248  DVector Legendre2D::definePY(int order, double y) const
249  {
250  if (y < getYMin() || y > getYMax()) {
251  dbg<<"Error: y = "<<y<<
252  ", min,max = "<<getYMin()<<','<<getYMax()<<std::endl;
254  }
255  return definePXY(order,y,getYMin(),getYMax());
256  }
257 
258 
259 
260 }}}}
int y
virtual void write(std::ostream &fout) const
Definition: Legendre2D.cc:85
virtual std::auto_ptr< Function2D > dFdY() const
Definition: Legendre2D.cc:189
virtual DVector definePY(int order, double y) const
Definition: Legendre2D.cc:248
DVector definePXY(int order, double xy, double min, double max) const
Definition: Legendre2D.cc:227
virtual void addLinear(double a, double b, double c)
Definition: Legendre2D.cc:107
int x
Eigen::VectorXd _b
tuple m
Definition: lsstimport.py:48
virtual void linearPreTransform(double a, double b, double c, double d, double e, double f)
Definition: Legendre2D.cc:119
virtual DVector definePX(int order, double x) const
Definition: Legendre2D.cc:238
afw::table::Key< double > b
RangeException(const Position &p, const Bounds &b)
Definition: Legendre2D.cc:47
#define dbg
Definition: dbg.h:69
virtual std::auto_ptr< Function2D > dFdX() const
Definition: Legendre2D.cc:149
virtual void operator+=(const Function2D &rhs)
Definition: Legendre2D.cc:126
virtual void setFunction(int xorder, int yorder, const DVector &fvect)
Definition: Legendre2D.cc:51
#define Assert(x)
Definition: dbg.h:73