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
Function2D.h
Go to the documentation of this file.
1 #ifndef MeasAlgoShapeletFunction2D_H
2 #define MeasAlgoShapeletFunction2D_H
3 
4 #include <iostream>
5 #include <functional>
6 #include <memory>
10 
11 namespace lsst {
12 namespace meas {
13 namespace algorithms {
14 namespace shapelet {
15 
16  struct RangeException : public std::runtime_error
17  {
18  public :
19  RangeException(const Position& p, const Bounds& b);
20  ~RangeException() throw() {}
21 
22  const Bounds& getBounds() const { return _b; }
23  const Position& getPosition() const { return _p; }
24 
25  private :
28 
29  };
30 
31  class Function2D
32  {
33  public:
34 
35  Function2D() : _xOrder(0), _yOrder(0), _coeffs(new DMatrix(1,1))
36  { (*_coeffs)(0,0) = 0.0; }
37 
38  Function2D(int xo, int yo) :
39  _xOrder(xo), _yOrder(yo), _coeffs(new DMatrix(xo+1,yo+1))
40  { _coeffs->setZero(); }
41 
42  Function2D(int xo, int yo, const DMatrix& c) :
43  _xOrder(xo),_yOrder(yo),_coeffs(new DMatrix(c)) {}
44 
45  Function2D(const Function2D& rhs) :
46  _xOrder(rhs._xOrder), _yOrder(rhs._yOrder),
47  _coeffs(new DMatrix(*rhs._coeffs)) {}
48 
49  virtual ~Function2D() {}
50 
51  virtual void write(std::ostream& fout) const =0;
52 
53  static std::auto_ptr<Function2D> read(std::istream& fin);
54 
55  virtual std::auto_ptr<Function2D> copy() const =0;
56 
57  // Adds a + bx + cy to function
58  virtual void addLinear(double a, double b, double c) = 0;
59 
60  // Converts function to f(a+bx+cy,d+ex+fy)
61  virtual void linearPreTransform(
62  double a, double b, double c, double d, double e, double f) = 0;
63 
64  // returns new function x derivative of f
65  virtual std::auto_ptr<Function2D> dFdX() const = 0;
66 
67  // returns new function y derivative of f
68  virtual std::auto_ptr<Function2D> dFdY() const = 0;
69 
70  virtual std::auto_ptr<Function2D> conj() const;
71 
72  virtual void operator*=(double scale)
73  { *_coeffs *= scale; }
74 
75  virtual double operator()(double x,double y) const;
76 
77  double operator()(const Position& p) const
78  { return operator()(p.getX(),p.getY()); }
79 
80  virtual void setTo(double value)
81  {
82  if (_xOrder || _yOrder) {
83  _xOrder = 0; _yOrder = 0;
84  _coeffs.reset(new DMatrix(1,1));
85  }
86  (*_coeffs)(0,0) = value;
87  }
88 
89  bool isNonZero() const
90  { return (*_coeffs)(0,0) != 0.0 || _xOrder!=0 || _yOrder!=0; }
91 
92  int getXOrder() const { return _xOrder; }
93 
94  int getYOrder() const { return _yOrder; }
95 
96  const DMatrix& getCoeffs() const { return *_coeffs; }
97 
98  // Sets function to fit of f(pos_i) = v_i using i only if
99  // shouldUse[i] = true
100  virtual void simpleFit(
101  int order, const std::vector<Position>& pos,
102  const std::vector<double>& v, const std::vector<bool>& shouldUse,
103  const std::vector<double>* sigList=0,
104  double* chisqOut = 0, int* dofOut=0, DMatrix* cov=0);
105 
106  // Sets function to fit of f(pos_i) = v_i using i if fit is within
107  // nsig sigma. *shouldUse is returned as list of good i's.
108  virtual void outlierFit(
109  int order,double nsig,
110  const std::vector<Position>& pos, const std::vector<double>& v,
111  std::vector<bool>* shouldUse,
112  const std::vector<double>* sigList=0,
113  double* chisqout = 0, int* dofout = 0, DMatrix* cov=0);
114 
115  // Sets function to fit of f(pos_i) = v_i reducing the order as far
116  // as possible keeping quality of fit the same as for maxOrder
117  // equivProb = rejection percentile. eg. 0.9 means a low order fit is
118  // rejected if it is statistically rejected at the 90th percentile
119  // Higher values of equivProb result in lower order fits.
120  virtual void orderFit(
121  int maxOrder, double equivProb,
122  const std::vector<Position>& pos, const std::vector<double>& v,
123  const std::vector<bool>& shouldUse, const std::vector<double>* sigList=0,
124  double *chisqout = 0, int *dofout = 0, DMatrix* cov=0);
125 
126  // Sets function to h(x,y) = a + b*f(x,y) + c*g(x,y)
127  virtual void linearTransform(
128  double a, double b, double c,
129  const Function2D& f, const Function2D& g);
130 
131  // Sets function to g(x,y) = a + b*f(x,y)
132  virtual void linearTransform(double a, double b, const Function2D& f)
133  { linearTransform(a,b,0.,f,f); }
134 
135  virtual void operator+=(const Function2D& rhs) = 0;
136 
137  virtual void setFunction(
138  int _xOrder, int _yOrder, const DVector& fVect) = 0;
139 
140  protected:
141 
142  virtual DVector definePX(int order, double x) const = 0;
143  virtual DVector definePY(int order, double y) const = 0;
144 
146  std::auto_ptr<DMatrix> _coeffs;
147 
148  private:
149 
150  void doSimpleFit(
151  int xOrder, int yOrder,
152  const std::vector<Position>& pos, const std::vector<double>& v,
153  const std::vector<bool>& shouldUse, DVector *f,
154  const std::vector<double>* sigList=0, int *dof=0,
155  DVector *diff=0, DMatrix* cov=0);
156  };
157 
158  inline std::ostream& operator<<(std::ostream& fout, const Function2D& f)
159  { f.write(fout); return fout; }
160 
161  inline std::istream& operator>>(
162  std::istream& fin, std::auto_ptr<Function2D>& f)
163  { f = Function2D::read(fin); return fin; }
164 
165  class Constant2D : public Function2D
166  {
167 
168  public:
169 
171 
172  Constant2D(const Constant2D& rhs) : Function2D(rhs) {}
173 
174  Constant2D(double value) : Function2D(0,0)
175  { Function2D::setTo(value); }
176 
177  Constant2D(std::istream& fin);
178 
179  virtual ~Constant2D() {}
180 
181  virtual double operator()(double ,double ) const
182  { return (*this->_coeffs)(0,0); }
183 
184  virtual void write(std::ostream& fout) const;
185 
186  virtual std::auto_ptr<Function2D> dFdX() const
187  { return std::auto_ptr<Function2D>(new Constant2D()); }
188 
189  virtual std::auto_ptr<Function2D> dFdY() const
190  { return std::auto_ptr<Function2D>(new Constant2D()); }
191 
192  std::auto_ptr<Function2D> copy() const
193  { return std::auto_ptr<Function2D>(new Constant2D(*this)); }
194 
195  virtual void addLinear(double a, double b, double c)
196  {
197  Assert(b == 0.);
198  Assert(c == 0.);
199  (*this->_coeffs)(0,0) += a;
200  }
201 
202  virtual void linearPreTransform(
203  double , double , double , double , double , double )
204  { Assert(false); }
205 
206  virtual void operator+=(const Function2D& rhs);
207 
208  virtual void setFunction(
209  int /*xOrder*/, int /*yOrder*/, const DVector& fVect)
210  {
211  Assert(_xOrder == 0);
212  Assert(_yOrder == 0);
213  Assert(fVect.size()==1);
214  (*this->_coeffs)(0,0) = fVect(0);
215  }
216 
217  private:
218 
219  virtual DVector definePX(int order, double ) const
220  {
221  Assert(order == 0);
222  DVector temp(1);
223  temp(0) = 1.0;
224  return temp;
225  }
226  virtual DVector definePY(int order, double y) const
227  { return definePX(order,y); }
228 
229  using Function2D::_coeffs;
230  using Function2D::_xOrder;
231  using Function2D::_yOrder;
232  };
233 
234  class Polynomial2D : public Function2D
235  {
236 
237  public:
238 
240 
242  Function2D(rhs),_scale(rhs._scale) {}
243 
244  Polynomial2D(const DMatrix& a, double scale=1.) :
245  Function2D(a.TMV_colsize()-1,a.TMV_rowsize()-1,a), _scale(scale) {}
246 
247  Polynomial2D(std::istream& fin);
248 
249  Polynomial2D(int xo,int yo,double scale=1.) :
250  Function2D(xo,yo), _scale(scale) {}
251 
252  virtual ~Polynomial2D() {}
253 
254  virtual void write(std::ostream& fout) const;
255 
256  virtual std::auto_ptr<Function2D> dFdX() const;
257 
258  virtual std::auto_ptr<Function2D> dFdY() const;
259 
260  std::auto_ptr<Function2D> copy() const
261  { return std::auto_ptr<Function2D>(new Polynomial2D(*this)); }
262 
263  virtual void addLinear(double a, double b, double c);
264 
265  virtual void linearPreTransform(
266  double a, double b, double c, double d, double e, double f);
267 
268  virtual void operator+=(const Function2D& rhs);
269 
270  void makeProductOf(const Polynomial2D& f, const Polynomial2D& g);
271 
272  virtual void setFunction(
273  int xOrder, int yOrder, const DVector& fVect);
274 
275  private:
276 
277  double _scale;
278 
279  virtual DVector definePX(int order, double x) const
280  {
281  DVector temp(order+1);
282  temp(0) = 1.0;
283  for(int i=1;i<=order;++i) temp(i) = temp(i-1)*x/_scale;
284  return temp;
285  }
286  virtual DVector definePY(int order, double y) const
287  { return definePX(order,y); }
288 
289  using Function2D::_coeffs;
290  using Function2D::_xOrder;
291  using Function2D::_yOrder;
292  };
293 
294 }}}}
295 
296 #endif
int y
virtual void setTo(double value)
Definition: Function2D.h:80
virtual std::auto_ptr< Function2D > dFdX() const
Definition: Function2D.h:186
virtual std::auto_ptr< Function2D > dFdY() const =0
virtual double operator()(double, double) const
Definition: Function2D.h:181
virtual DVector definePY(int order, double y) const =0
Polynomial2D(const DMatrix &a, double scale=1.)
Definition: Function2D.h:244
#define TMV_rowsize()
Definition: MyMatrix.h:174
virtual DVector definePX(int order, double x) const =0
std::auto_ptr< Function2D > copy() const
Definition: Function2D.h:260
std::auto_ptr< Function2D > copy() const
Definition: Function2D.h:192
virtual std::auto_ptr< Function2D > conj() const
Definition: Function2D.cc:315
virtual DVector definePX(int order, double) const
Definition: Function2D.h:219
virtual void linearTransform(double a, double b, const Function2D &f)
Definition: Function2D.h:132
virtual void addLinear(double a, double b, double c)=0
virtual void setFunction(int _xOrder, int _yOrder, const DVector &fVect)=0
virtual void operator+=(const Function2D &rhs)=0
std::ostream & operator<<(std::ostream &os, const Position &pos)
Definition: Bounds.h:78
virtual void outlierFit(int order, double nsig, const std::vector< Position > &pos, const std::vector< double > &v, std::vector< bool > *shouldUse, const std::vector< double > *sigList=0, double *chisqout=0, int *dofout=0, DMatrix *cov=0)
Definition: Function2D.cc:534
#define TMV_colsize()
Definition: MyMatrix.h:175
virtual void linearPreTransform(double, double, double, double, double, double)
Definition: Function2D.h:202
virtual void write(std::ostream &fout) const
Definition: Function2D.cc:44
virtual void linearPreTransform(double a, double b, double c, double d, double e, double f)
Definition: Function2D.cc:156
Function2D(int xo, int yo, const DMatrix &c)
Definition: Function2D.h:42
virtual std::auto_ptr< Function2D > dFdY() const
Definition: Function2D.cc:288
virtual std::auto_ptr< Function2D > dFdX() const =0
virtual void operator+=(const Function2D &rhs)
Definition: Function2D.cc:214
virtual void write(std::ostream &fout) const
Definition: Function2D.cc:122
virtual DVector definePY(int order, double y) const
Definition: Function2D.h:286
virtual std::auto_ptr< Function2D > dFdX() const
Definition: Function2D.cc:260
virtual void operator*=(double scale)
Definition: Function2D.h:72
double operator()(const Position &p) const
Definition: Function2D.h:77
virtual void addLinear(double a, double b, double c)
Definition: Function2D.cc:149
int x
std::istream & operator>>(std::istream &os, Position &pos)
Definition: Bounds.h:81
void makeProductOf(const Polynomial2D &f, const Polynomial2D &g)
Definition: Function2D.cc:236
virtual void linearTransform(double a, double b, double c, const Function2D &f, const Function2D &g)
Definition: Function2D.cc:355
virtual double operator()(double x, double y) const
Definition: Function2D.cc:322
afw::table::Key< double > b
virtual std::auto_ptr< Function2D > dFdY() const
Definition: Function2D.h:189
virtual void linearPreTransform(double a, double b, double c, double d, double e, double f)=0
virtual void simpleFit(int order, const std::vector< Position > &pos, const std::vector< double > &v, const std::vector< bool > &shouldUse, const std::vector< double > *sigList=0, double *chisqOut=0, int *dofOut=0, DMatrix *cov=0)
Definition: Function2D.cc:513
virtual void setFunction(int xOrder, int yOrder, const DVector &fVect)
Definition: Function2D.cc:63
void doSimpleFit(int xOrder, int yOrder, const std::vector< Position > &pos, const std::vector< double > &v, const std::vector< bool > &shouldUse, DVector *f, const std::vector< double > *sigList=0, int *dof=0, DVector *diff=0, DMatrix *cov=0)
Definition: Function2D.cc:404
RangeException(const Position &p, const Bounds &b)
Definition: Legendre2D.cc:47
virtual void setFunction(int, int, const DVector &fVect)
Definition: Function2D.h:208
virtual void addLinear(double a, double b, double c)
Definition: Function2D.h:195
virtual void write(std::ostream &fout) const =0
virtual std::auto_ptr< Function2D > copy() const =0
virtual void operator+=(const Function2D &rhs)
Definition: Function2D.cc:55
#define Assert(x)
Definition: dbg.h:73
virtual DVector definePX(int order, double x) const
Definition: Function2D.h:279
Polynomial2D(int xo, int yo, double scale=1.)
Definition: Function2D.h:249
static std::auto_ptr< Function2D > read(std::istream &fin)
Definition: Function2D.cc:330
virtual void orderFit(int maxOrder, double equivProb, const std::vector< Position > &pos, const std::vector< double > &v, const std::vector< bool > &shouldUse, const std::vector< double > *sigList=0, double *chisqout=0, int *dofout=0, DMatrix *cov=0)
Definition: Function2D.cc:602
virtual DVector definePY(int order, double y) const
Definition: Function2D.h:226