LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
Mixture.h
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008-2013 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 #ifndef LSST_MEAS_MODELFIT_Mixture_h_INCLUDED
25 #define LSST_MEAS_MODELFIT_Mixture_h_INCLUDED
26 
27 #include <limits>
28 
29 #include "Eigen/Cholesky"
30 #include "Eigen/StdVector"
31 #include "Eigen/Dense"
32 
33 #include "ndarray.h"
34 
35 #include "lsst/base.h"
36 #include "lsst/afw/math/Random.h"
39 #include "lsst/afw/table/io/python.h" // for declarePersistableFacade
40 
41 
42 namespace lsst { namespace meas { namespace modelfit {
43 
48 public:
49 
51  int getDimension() const { return _mu.size(); }
52 
55 
57  Vector getMu() const { return _mu; }
59  void setMu(Vector const & mu) { _mu = mu; }
61 
63 
71  Matrix getSigma() const { return _sigmaLLT.reconstructedMatrix(); }
72  void setSigma(Matrix const & sigma);
74 
76  MixtureComponent project(int dim) const;
77 
79  MixtureComponent project(int dim1, int dim2) const;
80 
82  explicit MixtureComponent(int dim);
83 
85  MixtureComponent(Scalar weight_, Vector const & mu, Matrix const & sigma);
86 
88 
90  self._stream(os);
91  return os;
92  }
93 
94 private:
95 
96  friend class Mixture;
97 
98  void _stream(std::ostream & os, int offset=0) const;
99 
100  Scalar _sqrtDet;
101  Vector _mu;
102  Eigen::LLT<Matrix> _sigmaLLT;
103 };
104 
112 public:
113 
114  int getDimension() const { return _dim; }
115 
116  virtual void restrictMu(Vector & mu) const {}
117 
118  virtual void restrictSigma(Matrix & sigma) const {}
119 
121 
122  explicit MixtureUpdateRestriction(int dim) : _dim(dim) {}
123 
124 private:
125  int _dim;
126 };
127 
129 public:
130 
134  typedef ComponentList::iterator iterator;
135  typedef ComponentList::const_iterator const_iterator;
136 
138 
146  iterator begin() { return _components.begin(); }
147  iterator end() { return _components.end(); }
148 
149  const_iterator begin() const { return _components.begin(); }
150  const_iterator end() const { return _components.end(); }
151 
152  Component & operator[](std::size_t i) { return _components[i]; }
153  Component const & operator[](std::size_t i) const { return _components[i]; }
155 
157  std::size_t size() const { return _components.size(); }
158 
160  virtual int getComponentCount() const { return size(); }
161 
163  PTR(Mixture) project(int dim) const;
164 
166  PTR(Mixture) project(int dim1, int dim2) const;
167 
169  int getDimension() const { return _dim; }
170 
172  void normalize();
173 
175  void shift(int dim, Scalar offset);
176 
184  std::size_t clip(Scalar threshold=0.0);
185 
187  Scalar getDegreesOfFreedom() const { return _df; }
188 
191 
197  template <typename Derived>
198  Scalar evaluate(Component const & component, Eigen::MatrixBase<Derived> const & x) const {
199  Scalar z = _computeZ(component, x);
200  return component.weight * _evaluate(z) / component._sqrtDet;
201  }
202 
208  template <typename Derived>
209  Scalar evaluate(Eigen::MatrixBase<Derived> const & x) const {
210  Scalar p = 0.0;
211  for (const_iterator i = begin(); i != end(); ++i) {
212  p += evaluate(*i, x);
213  }
214  return p;
215  }
216 
223  void evaluate(
224  ndarray::Array<Scalar const,2,1> const & x,
225  ndarray::Array<Scalar,1,0> const & p
226  ) const;
227 
235  ndarray::Array<Scalar const,2,1> const & x,
236  ndarray::Array<Scalar,2,1> const & p
237  ) const;
238 
247  ndarray::Array<Scalar const,1,1> const & x,
248  ndarray::Array<Scalar,1,1> const & gradient,
249  ndarray::Array<Scalar,2,1> const & hessian
250  ) const;
251 
260  Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & x,
261  Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & gradient,
262  Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & hessian
263  ) const;
264 
272  Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & x,
273  Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & gradient
274  ) const;
275 
282  void draw(afw::math::Random & rng, ndarray::Array<Scalar,2,1> const & x) const;
283 
306  void updateEM(
307  ndarray::Array<Scalar const,2,1> const & x,
308  ndarray::Array<Scalar const,1,0> const & w,
309  Scalar tau1=0.0, Scalar tau2=0.5
310  );
311 
322  void updateEM(
323  ndarray::Array<Scalar const,2,1> const & x,
324  ndarray::Array<Scalar const,1,0> const & w,
325  UpdateRestriction const & restriction,
326  Scalar tau1=0.0, Scalar tau2=0.5
327  );
328 
338  void updateEM(
339  ndarray::Array<Scalar const,2,1> const & x,
340  UpdateRestriction const & restriction,
341  Scalar tau1=0.0, Scalar tau2=0.5
342  );
343 
345  virtual PTR(Mixture) clone() const;
346 
358 
359  friend std::ostream & operator<<(std::ostream & os, Mixture const & self) {
360  self._stream(os);
361  return os;
362  }
363 
364  virtual bool isPersistable() const noexcept override { return true; }
365 
366 protected:
367 
368  std::string getPythonModule() const override { return "lsst.meas.modelfit"; }
369 
370  std::string getPersistenceName() const override;
371 
372  void write(OutputArchiveHandle & handle) const override;
373 
374 private:
375  template <typename A, typename B, typename C>
376  void _evaluateDerivativesImpl(A const & x,
377  B & gradient,
378  C * hessian,
379  bool computeHessian = true) const;
380 
381  template <typename Derived>
382  Scalar _computeZ(Component const & component, Eigen::MatrixBase<Derived> const & x) const {
383  _workspace = x - component._mu;
384  component._sigmaLLT.matrixL().solveInPlace(_workspace);
385  return _workspace.squaredNorm();
386  }
387 
388  // Helper function used in updateEM
389  void updateDampedSigma(int k, Matrix const & sigma, double tau1, double tau2);
390 
391  Scalar _evaluate(Scalar z) const;
392 
393  void _stream(std::ostream & os) const;
394 
395  bool _isGaussian;
396  int _dim;
397  Scalar _df;
398  Scalar _norm;
399  mutable Vector _workspace;
400  ComponentList _components;
401 };
402 
403 }}} // namespace lsst::meas::modelfit
404 
405 #endif // !LSST_MEAS_MODELFIT_Mixture_h_INCLUDED
double x
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
table::Key< table::Array< int > > components
double z
Definition: Match.cc:44
ItemVariant const * other
Definition: Schema.cc:56
std::ostream * os
Definition: Schema.cc:746
Basic LSST definitions.
#define PTR(...)
Definition: base.h:41
T begin(T... args)
A class that can be used to generate sequences of random numbers according to a number of different a...
Definition: Random.h:57
An object passed to Persistable::write to allow it to persist itself.
A CRTP facade class for subclasses of Persistable.
Definition: Persistable.h:176
A base class for objects that can be persisted via afw::table::io Archive classes.
Definition: Persistable.h:74
A weighted Student's T or Gaussian distribution used as a component in a Mixture.
Definition: Mixture.h:47
MixtureComponent(Scalar weight_, Vector const &mu, Matrix const &sigma)
Default-construct a mixture component with the given parameters.
MixtureComponent project(int dim) const
Project the distribution onto the given dimension (marginalize over all others)
MixtureComponent project(int dim1, int dim2) const
Project the distribution onto the given dimensions (marginalize over all others)
Vector getMu() const
Get/set the location parameter (mean/median/mode) of this component.
Definition: Mixture.h:58
Scalar weight
Weight of this distribution in the mixture.
Definition: Mixture.h:54
friend std::ostream & operator<<(std::ostream &os, MixtureComponent const &self)
Definition: Mixture.h:89
void setSigma(Matrix const &sigma)
void setMu(Vector const &mu)
Definition: Mixture.h:59
Matrix getSigma() const
Get/set the shape/size parameter.
Definition: Mixture.h:71
int getDimension() const
Return the number of dimensions.
Definition: Mixture.h:51
MixtureComponent & operator=(MixtureComponent const &other)
MixtureComponent(int dim)
Default-construct a mixture component with weight=1, mu=0, sigma=identity.
void updateEM(ndarray::Array< Scalar const, 2, 1 > const &x, ndarray::Array< Scalar const, 1, 0 > const &w, Scalar tau1=0.0, Scalar tau2=0.5)
Perform an Expectation-Maximization step, updating the component parameters to match the given weight...
Scalar getDegreesOfFreedom() const
Get the number of degrees of freedom in the component Student's T distributions (inf=Gaussian)
Definition: Mixture.h:187
void updateEM(ndarray::Array< Scalar const, 2, 1 > const &x, UpdateRestriction const &restriction, Scalar tau1=0.0, Scalar tau2=0.5)
Perform an Expectation-Maximization step, updating the component parameters to match the given unweig...
virtual boost::shared_ptr< Mixture > clone() const
Polymorphic deep copy.
std::size_t size() const
Return the number of components.
Definition: Mixture.h:157
const_iterator end() const
Definition: Mixture.h:150
Scalar evaluate(Eigen::MatrixBase< Derived > const &x) const
Evaluate the mixture distribution probability density function (PDF) at the given points.
Definition: Mixture.h:209
void evaluateDerivatives(Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > &x, Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > &gradient) const
Evaluate the derivative of the distribution at the given point.
virtual bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
Definition: Mixture.h:364
boost::shared_ptr< Mixture > project(int dim) const
Project the distribution onto the given dimensions (marginalize over all others)
Scalar evaluate(Component const &component, Eigen::MatrixBase< Derived > const &x) const
Evaluate the probability density at the given point for the given component distribution.
Definition: Mixture.h:198
const_iterator begin() const
Definition: Mixture.h:149
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
Definition: Mixture.h:368
void normalize()
Iterate over all components, rescaling their weights so they sum to one.
void evaluateComponents(ndarray::Array< Scalar const, 2, 1 > const &x, ndarray::Array< Scalar, 2, 1 > const &p) const
Evaluate the contributions of each component to the full probability at the given points.
void evaluateDerivatives(ndarray::Array< Scalar const, 1, 1 > const &x, ndarray::Array< Scalar, 1, 1 > const &gradient, ndarray::Array< Scalar, 2, 1 > const &hessian) const
Evaluate the derivative of the distribution at the given point.
ComponentList::iterator iterator
Definition: Mixture.h:134
void shift(int dim, Scalar offset)
Shift the mixture in the given dimension, adding the given offset to all mu vectors.
ComponentList::const_iterator const_iterator
Definition: Mixture.h:135
iterator begin()
Iterator and indexed access to components.
Definition: Mixture.h:146
Component const & operator[](std::size_t i) const
Definition: Mixture.h:153
friend std::ostream & operator<<(std::ostream &os, Mixture const &self)
Definition: Mixture.h:359
Component & operator[](std::size_t i)
Definition: Mixture.h:152
virtual int getComponentCount() const
Return the number of components.
Definition: Mixture.h:160
void evaluateDerivatives(Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > &x, Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > &gradient, Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > &hessian) const
Evaluate the derivative of the distribution at the given point.
Mixture(int dim, ComponentList &components, Scalar df=std::numeric_limits< Scalar >::infinity())
Construct a mixture model.
MixtureUpdateRestriction UpdateRestriction
Definition: Mixture.h:132
void draw(afw::math::Random &rng, ndarray::Array< Scalar, 2, 1 > const &x) const
Draw random variates from the distribution.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
void updateEM(ndarray::Array< Scalar const, 2, 1 > const &x, ndarray::Array< Scalar const, 1, 0 > const &w, UpdateRestriction const &restriction, Scalar tau1=0.0, Scalar tau2=0.5)
Perform an Expectation-Maximization step, updating the component parameters to match the given weight...
int getDimension() const
Return the number of dimensions.
Definition: Mixture.h:169
void setDegreesOfFreedom(Scalar df=std::numeric_limits< Scalar >::infinity())
Set the number of degrees of freedom in the component Student's T distributions (inf=Gaussian)
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
MixtureComponent Component
Definition: Mixture.h:131
void evaluate(ndarray::Array< Scalar const, 2, 1 > const &x, ndarray::Array< Scalar, 1, 0 > const &p) const
Evaluate the distribution probability density function (PDF) at the given points.
boost::shared_ptr< Mixture > project(int dim1, int dim2) const
Project the distribution onto the given dimensions (marginalize over all others)
std::size_t clip(Scalar threshold=0.0)
Iterate over all components, removing those with weight less than or equal to threshold.
std::vector< Component > ComponentList
Definition: Mixture.h:133
Helper class used to define restrictions to the form of the component parameters in Mixture::updateEM...
Definition: Mixture.h:111
virtual void restrictMu(Vector &mu) const
Definition: Mixture.h:116
virtual void restrictSigma(Matrix &sigma) const
Definition: Mixture.h:118
T end(T... args)
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > Vector
Definition: common.h:46
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > Matrix
Definition: common.h:45
double Scalar
Typedefs to be used for probability and parameter values.
Definition: common.h:44
A base class for image defects.
T size(T... args)
double w
Definition: CoaddPsf.cc:69