LSST Applications g1cfbe01172+01aa18f939,g20cdd03214+31e6b93548,g28da252d5a+ea8665a95b,g2bbee38e9b+9ec6cc348d,g2bc492864f+9ec6cc348d,g347aa1857d+9ec6cc348d,g3a166c0a6a+9ec6cc348d,g4322eb9e3a+65eff1e020,g461a3dce89+b86e4b8053,g50ff169b8f+f991eae79d,g52b1c1532d+b86e4b8053,g607f77f49a+31e6b93548,g78056777b3+8ae2798781,g858d7b2824+31e6b93548,g8cd86fa7b1+4851e61ca4,g9ddcbc5298+f24b38b85a,ga1e77700b3+3309dba821,gae0086650b+b86e4b8053,gb0e22166c9+6076c0b52b,gbb886bcc26+dccb771098,gbd462c55f0+dc07f8e65d,gc0c51c7ec2+31e6b93548,gc120e1dc64+a417ce3171,gc28159a63d+9ec6cc348d,gc2a6998b7e+f95f64aeae,gcdd4ae20e8+507450c4cd,gcf0d15dbbd+507450c4cd,gd1535ee943+bcf88ba65f,gd598c5cd71+66126f91fb,gdaeeff99f8+006e14e809,gdbce86181e+39d5515b1a,ge3d4d395c2+b12d4d6a95,ge79ae78c31+9ec6cc348d,gf048a9a2f4+d9c36e6b63,gfbcc870c63+ea41c4420b,w.2024.27
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
39#include "lsst/afw/table/io/python.h" // for declarePersistableFacade
40
41
42namespace lsst { namespace meas { namespace modelfit {
43
48public:
49
51 int getDimension() const { return _mu.size(); }
52
55
57
58 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
94private:
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
112public:
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
124private:
125 int _dim;
126};
127
129public:
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
164
166 std::shared_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
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
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
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
339 ndarray::Array<Scalar const,2,1> const & x,
340 UpdateRestriction const & restriction,
341 Scalar tau1=0.0, Scalar tau2=0.5
342 );
343
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
366protected:
367
368 std::string getPythonModule() const override { return "lsst.meas.modelfit"; }
369
371
372 void write(OutputArchiveHandle & handle) const override;
373
374private:
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
afw::table::Key< double > sigma
table::Key< table::Array< int > > components
double z
Definition Match.cc:44
std::ostream * os
Definition Schema.cc:557
Basic LSST definitions.
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.
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
friend std::ostream & operator<<(std::ostream &os, MixtureComponent const &self)
Definition Mixture.h:89
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
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
MixtureComponent & operator=(MixtureComponent const &other)
int getDimension() const
Return the number of dimensions.
Definition Mixture.h:51
MixtureComponent(int dim)
Default-construct a mixture component with weight=1, mu=0, sigma=identity.
Component const & operator[](std::size_t i) const
Definition Mixture.h:153
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...
std::size_t size() const
Return the number of components.
Definition Mixture.h:157
friend std::ostream & operator<<(std::ostream &os, Mixture const &self)
Definition Mixture.h:359
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.
std::shared_ptr< Mixture > project(int dim) const
Project the distribution onto the given dimensions (marginalize over all others)
virtual bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
Definition Mixture.h:364
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
Component & operator[](std::size_t i)
Definition Mixture.h:152
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
virtual std::shared_ptr< Mixture > clone() const
Polymorphic deep copy.
iterator begin()
Iterator and indexed access to components.
Definition Mixture.h:146
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
std::shared_ptr< Mixture > project(int dim1, int dim2) const
Project the distribution onto the given dimensions (marginalize over all others)
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.
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
T size(T... args)
double w
Definition CoaddPsf.cc:70