LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Public Types | Public Member Functions | Protected Attributes | Static Protected Attributes | List of all members
lsst::ip::diffim::SpatialKernelSolution Class Reference

#include <KernelSolution.h>

Inheritance diagram for lsst::ip::diffim::SpatialKernelSolution:
lsst::ip::diffim::KernelSolution

Public Types

typedef std::shared_ptr< SpatialKernelSolutionPtr
 
enum  KernelSolvedBy {
  NONE = 0 , CHOLESKY_LDLT = 1 , CHOLESKY_LLT = 2 , LU = 3 ,
  EIGENVECTOR = 4
}
 
enum  ConditionNumberType { EIGENVALUE = 0 , SVD = 1 }
 
typedef lsst::afw::math::Kernel::Pixel PixelT
 
typedef lsst::afw::image::Image< lsst::afw::math::Kernel::PixelImageT
 

Public Member Functions

 SpatialKernelSolution (lsst::afw::math::KernelList const &basisList, lsst::afw::math::Kernel::SpatialFunctionPtr spatialKernelFunction, lsst::afw::math::Kernel::SpatialFunctionPtr background, lsst::daf::base::PropertySet const &ps)
 
virtual ~SpatialKernelSolution ()
 
void addConstraint (float xCenter, float yCenter, Eigen::MatrixXd const &qMat, Eigen::VectorXd const &wVec)
 
void solve ()
 
std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage (lsst::geom::Point2D const &pos)
 
std::pair< std::shared_ptr< lsst::afw::math::LinearCombinationKernel >, lsst::afw::math::Kernel::SpatialFunctionPtrgetSolutionPair ()
 
virtual void solve (Eigen::MatrixXd const &mMat, Eigen::VectorXd const &bVec)
 
KernelSolvedBy getSolvedBy ()
 
virtual double getConditionNumber (ConditionNumberType conditionType)
 
virtual double getConditionNumber (Eigen::MatrixXd const &mMat, ConditionNumberType conditionType)
 
Eigen::MatrixXd const & getM ()
 
Eigen::VectorXd const & getB ()
 
void printM ()
 
void printB ()
 
void printA ()
 
int getId () const
 

Protected Attributes

int _id
 Unique ID for object. More...
 
Eigen::MatrixXd _mMat
 Derived least squares M matrix. More...
 
Eigen::VectorXd _bVec
 Derived least squares B vector. More...
 
Eigen::VectorXd _aVec
 Derived least squares solution matrix. More...
 
KernelSolvedBy _solvedBy
 Type of algorithm used to make solution. More...
 
bool _fitForBackground
 Background terms included in fit. More...
 

Static Protected Attributes

static int _SolutionId = 0
 Unique identifier for solution. More...
 

Detailed Description

Definition at line 174 of file KernelSolution.h.

Member Typedef Documentation

◆ ImageT

Definition at line 35 of file KernelSolution.h.

◆ PixelT

Definition at line 34 of file KernelSolution.h.

◆ Ptr

Definition at line 176 of file KernelSolution.h.

Member Enumeration Documentation

◆ ConditionNumberType

Enumerator
EIGENVALUE 
SVD 

Definition at line 45 of file KernelSolution.h.

◆ KernelSolvedBy

Enumerator
NONE 
CHOLESKY_LDLT 
CHOLESKY_LLT 
LU 
EIGENVECTOR 

Definition at line 37 of file KernelSolution.h.

Constructor & Destructor Documentation

◆ SpatialKernelSolution()

lsst::ip::diffim::SpatialKernelSolution::SpatialKernelSolution ( lsst::afw::math::KernelList const &  basisList,
lsst::afw::math::Kernel::SpatialFunctionPtr  spatialKernelFunction,
lsst::afw::math::Kernel::SpatialFunctionPtr  background,
lsst::daf::base::PropertySet const &  ps 
)

Definition at line 1291 of file KernelSolution.cc.

1296  :
1297  KernelSolution(),
1298  _spatialKernelFunction(spatialKernelFunction),
1299  _constantFirstTerm(false),
1300  _kernel(),
1301  _background(background),
1302  _kSum(0.0),
1303  _ps(ps.deepCopy()),
1304  _nbases(0),
1305  _nkt(0),
1306  _nbt(0),
1307  _nt(0) {
1308 
1309  bool isAlardLupton = _ps->getAsString("kernelBasisSet") == "alard-lupton";
1310  bool usePca = _ps->getAsBool("usePcaForSpatialKernel");
1311  if (isAlardLupton || usePca) {
1312  _constantFirstTerm = true;
1313  }
1314  this->_fitForBackground = _ps->getAsBool("fitForBackground");
1315 
1316  _nbases = basisList.size();
1317  _nkt = _spatialKernelFunction->getParameters().size();
1318  _nbt = _fitForBackground ? _background->getParameters().size() : 0;
1319  _nt = 0;
1320  if (_constantFirstTerm) {
1321  _nt = (_nbases - 1) * _nkt + 1 + _nbt;
1322  } else {
1323  _nt = _nbases * _nkt + _nbt;
1324  }
1325 
1326  Eigen::MatrixXd mMat(_nt, _nt);
1327  Eigen::VectorXd bVec(_nt);
1328  mMat.setZero();
1329  bVec.setZero();
1330 
1331  this->_mMat = mMat;
1332  this->_bVec = bVec;
1333 
1335  new afwMath::LinearCombinationKernel(basisList, *_spatialKernelFunction)
1336  );
1337 
1338  LOGL_DEBUG("TRACE3.ip.diffim.SpatialKernelSolution",
1339  "Initializing with size %d %d %d and constant first term = %s",
1340  _nkt, _nbt, _nt,
1341  _constantFirstTerm ? "true" : "false");
1342 
1343  }
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:704
bool _fitForBackground
Background terms included in fit.
Eigen::VectorXd _bVec
Derived least squares B vector.
Eigen::MatrixXd _mMat
Derived least squares M matrix.

◆ ~SpatialKernelSolution()

virtual lsst::ip::diffim::SpatialKernelSolution::~SpatialKernelSolution ( )
inlinevirtual

Definition at line 185 of file KernelSolution.h.

185 {};

Member Function Documentation

◆ addConstraint()

void lsst::ip::diffim::SpatialKernelSolution::addConstraint ( float  xCenter,
float  yCenter,
Eigen::MatrixXd const &  qMat,
Eigen::VectorXd const &  wVec 
)

Definition at line 1345 of file KernelSolution.cc.

1347  {
1348 
1349  LOGL_DEBUG("TRACE5.ip.diffim.SpatialKernelSolution.addConstraint",
1350  "Adding candidate at %f, %f", xCenter, yCenter);
1351 
1352  /* Calculate P matrices */
1353  /* Pure kernel terms */
1354  Eigen::VectorXd pK(_nkt);
1355  std::vector<double> paramsK = _spatialKernelFunction->getParameters();
1356  for (int idx = 0; idx < _nkt; idx++) { paramsK[idx] = 0.0; }
1357  for (int idx = 0; idx < _nkt; idx++) {
1358  paramsK[idx] = 1.0;
1359  _spatialKernelFunction->setParameters(paramsK);
1360  pK(idx) = (*_spatialKernelFunction)(xCenter, yCenter); /* Assume things don't vary over stamp */
1361  paramsK[idx] = 0.0;
1362  }
1363  Eigen::MatrixXd pKpKt = (pK * pK.transpose());
1364 
1365  Eigen::VectorXd pB;
1366  Eigen::MatrixXd pBpBt;
1367  Eigen::MatrixXd pKpBt;
1368  if (_fitForBackground) {
1369  pB = Eigen::VectorXd(_nbt);
1370 
1371  /* Pure background terms */
1372  std::vector<double> paramsB = _background->getParameters();
1373  for (int idx = 0; idx < _nbt; idx++) { paramsB[idx] = 0.0; }
1374  for (int idx = 0; idx < _nbt; idx++) {
1375  paramsB[idx] = 1.0;
1376  _background->setParameters(paramsB);
1377  pB(idx) = (*_background)(xCenter, yCenter); /* Assume things don't vary over stamp */
1378  paramsB[idx] = 0.0;
1379  }
1380  pBpBt = (pB * pB.transpose());
1381 
1382  /* Cross terms */
1383  pKpBt = (pK * pB.transpose());
1384  }
1385 
1386  if (DEBUG_MATRIX) {
1387  std::cout << "Spatial weights" << std::endl;
1388  std::cout << "pKpKt " << pKpKt << std::endl;
1389  if (_fitForBackground) {
1390  std::cout << "pBpBt " << pBpBt << std::endl;
1391  std::cout << "pKpBt " << pKpBt << std::endl;
1392  }
1393  }
1394 
1395  if (DEBUG_MATRIX) {
1396  std::cout << "Spatial matrix inputs" << std::endl;
1397  std::cout << "M " << qMat << std::endl;
1398  std::cout << "B " << wVec << std::endl;
1399  }
1400 
1401  /* first index to start the spatial blocks; default=0 for non-constant first term */
1402  int m0 = 0;
1403  /* how many rows/cols to adjust the matrices/vectors; default=0 for non-constant first term */
1404  int dm = 0;
1405  /* where to start the background terms; this is always true */
1406  int mb = _nt - _nbt;
1407 
1408  if (_constantFirstTerm) {
1409  m0 = 1; /* we need to manually fill in the first (non-spatial) terms below */
1410  dm = _nkt-1; /* need to shift terms due to lack of spatial variation in first term */
1411 
1412  _mMat(0, 0) += qMat(0,0);
1413  for(int m2 = 1; m2 < _nbases; m2++) {
1414  _mMat.block(0, m2*_nkt-dm, 1, _nkt) += qMat(0,m2) * pK.transpose();
1415  }
1416  _bVec(0) += wVec(0);
1417 
1418  if (_fitForBackground) {
1419  _mMat.block(0, mb, 1, _nbt) += qMat(0,_nbases) * pB.transpose();
1420  }
1421  }
1422 
1423  /* Fill in the spatial blocks */
1424  for(int m1 = m0; m1 < _nbases; m1++) {
1425  /* Diagonal kernel-kernel term; only use upper triangular part of pKpKt */
1426  _mMat.block(m1*_nkt-dm, m1*_nkt-dm, _nkt, _nkt) +=
1427  (pKpKt * qMat(m1,m1)).triangularView<Eigen::Upper>();
1428 
1429  /* Kernel-kernel terms */
1430  for(int m2 = m1+1; m2 < _nbases; m2++) {
1431  _mMat.block(m1*_nkt-dm, m2*_nkt-dm, _nkt, _nkt) += qMat(m1,m2) * pKpKt;
1432  }
1433 
1434  if (_fitForBackground) {
1435  /* Kernel cross terms with background */
1436  _mMat.block(m1*_nkt-dm, mb, _nkt, _nbt) += qMat(m1,_nbases) * pKpBt;
1437  }
1438 
1439  /* B vector */
1440  _bVec.segment(m1*_nkt-dm, _nkt) += wVec(m1) * pK;
1441  }
1442 
1443  if (_fitForBackground) {
1444  /* Background-background terms only */
1445  _mMat.block(mb, mb, _nbt, _nbt) +=
1446  (pBpBt * qMat(_nbases,_nbases)).triangularView<Eigen::Upper>();
1447  _bVec.segment(mb, _nbt) += wVec(_nbases) * pB;
1448  }
1449 
1450  if (DEBUG_MATRIX) {
1451  std::cout << "Spatial matrix outputs" << std::endl;
1452  std::cout << "mMat " << _mMat << std::endl;
1453  std::cout << "bVec " << _bVec << std::endl;
1454  }
1455 
1456  }
T endl(T... args)
#define DEBUG_MATRIX

◆ getB()

Eigen::VectorXd const& lsst::ip::diffim::KernelSolution::getB ( )
inlineinherited

Definition at line 65 of file KernelSolution.h.

65 {return _bVec;}

◆ getConditionNumber() [1/2]

double lsst::ip::diffim::KernelSolution::getConditionNumber ( ConditionNumberType  conditionType)
virtualinherited

Definition at line 94 of file KernelSolution.cc.

94  {
95  return getConditionNumber(_mMat, conditionType);
96  }
virtual double getConditionNumber(ConditionNumberType conditionType)

◆ getConditionNumber() [2/2]

double lsst::ip::diffim::KernelSolution::getConditionNumber ( Eigen::MatrixXd const &  mMat,
ConditionNumberType  conditionType 
)
virtualinherited

Definition at line 98 of file KernelSolution.cc.

99  {
100  switch (conditionType) {
101  case EIGENVALUE:
102  {
103  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
104  Eigen::VectorXd eValues = eVecValues.eigenvalues();
105  double eMax = eValues.maxCoeff();
106  double eMin = eValues.minCoeff();
107  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
108  "EIGENVALUE eMax / eMin = %.3e", eMax / eMin);
109  return (eMax / eMin);
110  break;
111  }
112  case SVD:
113  {
114  Eigen::VectorXd sValues = mMat.jacobiSvd().singularValues();
115  double sMax = sValues.maxCoeff();
116  double sMin = sValues.minCoeff();
117  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
118  "SVD eMax / eMin = %.3e", sMax / sMin);
119  return (sMax / sMin);
120  break;
121  }
122  default:
123  {
125  "Undefined ConditionNumberType : only EIGENVALUE, SVD allowed.");
126  break;
127  }
128  }
129  }
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Reports invalid arguments.
Definition: Runtime.h:66

◆ getId()

int lsst::ip::diffim::KernelSolution::getId ( ) const
inlineinherited

Definition at line 69 of file KernelSolution.h.

69 { return _id; }
int _id
Unique ID for object.

◆ getM()

Eigen::MatrixXd const& lsst::ip::diffim::KernelSolution::getM ( )
inlineinherited

Definition at line 64 of file KernelSolution.h.

64 {return _mMat;}

◆ getSolutionPair()

std::pair< std::shared_ptr< afwMath::LinearCombinationKernel >, afwMath::Kernel::SpatialFunctionPtr > lsst::ip::diffim::SpatialKernelSolution::getSolutionPair ( )

Definition at line 1488 of file KernelSolution.cc.

1488  {
1490  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
1491  }
1492 
1493  return std::make_pair(_kernel, _background);
1494  }
KernelSolvedBy _solvedBy
Type of algorithm used to make solution.
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
T make_pair(T... args)

◆ getSolvedBy()

KernelSolvedBy lsst::ip::diffim::KernelSolution::getSolvedBy ( )
inlineinherited

Definition at line 60 of file KernelSolution.h.

60 {return _solvedBy;}

◆ makeKernelImage()

std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > lsst::ip::diffim::SpatialKernelSolution::makeKernelImage ( lsst::geom::Point2D const &  pos)

Definition at line 1458 of file KernelSolution.cc.

1458  {
1460  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
1461  }
1464  );
1465  (void)_kernel->computeImage(*image, false, pos[0], pos[1]);
1466  return image;
1467  }
afw::table::Key< afw::table::Array< ImagePixelT > > image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
lsst::geom::Extent2I const getDimensions() const
Return the Kernel's dimensions (width, height)
Definition: Kernel.h:212
double computeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize, double x=0.0, double y=0.0) const
Compute an image (pixellized representation of the kernel) in place.
Definition: Kernel.cc:76
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.

◆ printA()

void lsst::ip::diffim::KernelSolution::printA ( )
inlineinherited

Definition at line 68 of file KernelSolution.h.

68 {std::cout << _aVec << std::endl;}
Eigen::VectorXd _aVec
Derived least squares solution matrix.

◆ printB()

void lsst::ip::diffim::KernelSolution::printB ( )
inlineinherited

Definition at line 67 of file KernelSolution.h.

67 {std::cout << _bVec << std::endl;}

◆ printM()

void lsst::ip::diffim::KernelSolution::printM ( )
inlineinherited

Definition at line 66 of file KernelSolution.h.

66 {std::cout << _mMat << std::endl;}

◆ solve() [1/2]

void lsst::ip::diffim::SpatialKernelSolution::solve ( )
virtual

Reimplemented from lsst::ip::diffim::KernelSolution.

Definition at line 1469 of file KernelSolution.cc.

1469  {
1470  /* Fill in the other half of mMat */
1471  for (int i = 0; i < _nt; i++) {
1472  for (int j = i+1; j < _nt; j++) {
1473  _mMat(j,i) = _mMat(i,j);
1474  }
1475  }
1476 
1477  try {
1479  } catch (pexExcept::Exception &e) {
1480  LSST_EXCEPT_ADD(e, "Unable to solve spatial kernel matrix");
1481  throw e;
1482  }
1483  /* Turn matrices into _kernel and _background */
1484  _setKernel();
1485  }
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition: Exception.h:54

◆ solve() [2/2]

void lsst::ip::diffim::KernelSolution::solve ( Eigen::MatrixXd const &  mMat,
Eigen::VectorXd const &  bVec 
)
virtualinherited

Definition at line 131 of file KernelSolution.cc.

132  {
133 
134  if (DEBUG_MATRIX) {
135  std::cout << "M " << std::endl;
136  std::cout << mMat << std::endl;
137  std::cout << "B " << std::endl;
138  std::cout << bVec << std::endl;
139  }
140 
141  Eigen::VectorXd aVec = Eigen::VectorXd::Zero(bVec.size());
142 
143  boost::timer t;
144  t.restart();
145 
146  LOGL_DEBUG("TRACE2.ip.diffim.KernelSolution.solve",
147  "Solving for kernel");
148  _solvedBy = LU;
149  Eigen::FullPivLU<Eigen::MatrixXd> lu(mMat);
150  if (lu.isInvertible()) {
151  aVec = lu.solve(bVec);
152  } else {
153  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
154  "Unable to determine kernel via LU");
155  /* LAST RESORT */
156  try {
157 
159  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
160  Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
161  Eigen::VectorXd eValues = eVecValues.eigenvalues();
162 
163  for (int i = 0; i != eValues.rows(); ++i) {
164  if (eValues(i) != 0.0) {
165  eValues(i) = 1.0/eValues(i);
166  }
167  }
168 
169  aVec = rMat * eValues.asDiagonal() * rMat.transpose() * bVec;
170  } catch (pexExcept::Exception& e) {
171 
172  _solvedBy = NONE;
173  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
174  "Unable to determine kernel via eigen-values");
175 
176  throw LSST_EXCEPT(pexExcept::Exception, "Unable to determine kernel solution");
177  }
178  }
179 
180  double time = t.elapsed();
181  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
182  "Compute time for matrix math : %.2f s", time);
183 
184  if (DEBUG_MATRIX) {
185  std::cout << "A " << std::endl;
186  std::cout << aVec << std::endl;
187  }
188 
189  _aVec = aVec;
190  }
T time(T... args)

Member Data Documentation

◆ _aVec

Eigen::VectorXd lsst::ip::diffim::KernelSolution::_aVec
protectedinherited

Derived least squares solution matrix.

Definition at line 75 of file KernelSolution.h.

◆ _bVec

Eigen::VectorXd lsst::ip::diffim::KernelSolution::_bVec
protectedinherited

Derived least squares B vector.

Definition at line 74 of file KernelSolution.h.

◆ _fitForBackground

bool lsst::ip::diffim::KernelSolution::_fitForBackground
protectedinherited

Background terms included in fit.

Definition at line 77 of file KernelSolution.h.

◆ _id

int lsst::ip::diffim::KernelSolution::_id
protectedinherited

Unique ID for object.

Definition at line 72 of file KernelSolution.h.

◆ _mMat

Eigen::MatrixXd lsst::ip::diffim::KernelSolution::_mMat
protectedinherited

Derived least squares M matrix.

Definition at line 73 of file KernelSolution.h.

◆ _SolutionId

int lsst::ip::diffim::KernelSolution::_SolutionId = 0
staticprotectedinherited

Unique identifier for solution.

Definition at line 78 of file KernelSolution.h.

◆ _solvedBy

KernelSolvedBy lsst::ip::diffim::KernelSolution::_solvedBy
protectedinherited

Type of algorithm used to make solution.

Definition at line 76 of file KernelSolution.h.


The documentation for this class was generated from the following files: