Skip to content

Commit

Permalink
Eliminate shared ptrs to Eigen matrices & vectors
Browse files Browse the repository at this point in the history 
Eliminate all use of shared pointers to Eigen matrices
and vectors, especially as arguments and return values.
The impetus was that pybind11 did not wrap the return
value from the following properly (it build but could
not turn the result into a Python object, and the error
message made it look like the result was an Eigen::MatrixXd
instead of a shared pointer to the same):
```
std::shared_ptr<Eigen::MatrixXd> makeRegularizationMatrix
```
  • Loading branch information
@r-owen @TallJimbo
r-owen authored and TallJimbo committed Mar 6, 2017
1 parent 8227541 commit 06a92e4
Show file tree
Hide file tree
Showing 13 changed files with 249 additions and 298 deletions.
28 changes: 14 additions & 14 deletions examples/pixelAccess.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,16 +36,16 @@ Eigen::MatrixXd test(afwImage::Image<ImageT> varianceEstimate,

if (cswitch == 3) {
/* a list of images - in diffim each one of these is associated with a basis function */
std::vector<std::shared_ptr<Eigen::VectorXd> > imageList(nParameters);
typename std::vector<std::shared_ptr<Eigen::VectorXd> >::iterator eiter = imageList.begin();
std::vector<Eigen::VectorXd> imageList(nParameters);
typename std::vector<Eigen::VectorXd>::iterator eiter = imageList.begin();
afwImage::Image<ImageT> cimage(varianceEstimate.getDimensions());
for (int i = 1; eiter != imageList.end(); ++eiter, ++i) {
cimage = i; /* give it a value */
Eigen::MatrixXd cMat = diffim::imageToEigenMatrix(cimage).block(startRow, startCol,
endRow-startRow,
endCol-startCol);
cMat.resize(cMat.rows()*cMat.cols(), 1);
std::shared_ptr<Eigen::VectorXd> vMat (new Eigen::VectorXd(cMat.col(0)));
Eigen::VectorXd vMat = cMat.col(0);
*eiter = vMat;
}

Expand All @@ -58,11 +58,11 @@ Eigen::MatrixXd test(afwImage::Image<ImageT> varianceEstimate,
Eigen::VectorXd eigeniVarianceV = eigeniVarianceM.col(0);

Eigen::MatrixXd cMat(eigeniVarianceV.size(), nParameters);
typename std::vector<std::shared_ptr<Eigen::VectorXd> >::iterator eiterj = imageList.begin();
typename std::vector<std::shared_ptr<Eigen::VectorXd> >::iterator eiterE = imageList.end();
typename std::vector<Eigen::VectorXd>::iterator eiterj = imageList.begin();
typename std::vector<Eigen::VectorXd>::iterator eiterE = imageList.end();
for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
cMat.block(0, kidxj, eigeniVarianceV.size(), 1) =
Eigen::MatrixXd(**eiterj).block(0, 0, eigeniVarianceV.size(), 1);
Eigen::MatrixXd(*eiterj).block(0, 0, eigeniVarianceV.size(), 1);
}

// Caculate the variance-weighted pixel values
Expand All @@ -75,8 +75,8 @@ Eigen::MatrixXd test(afwImage::Image<ImageT> varianceEstimate,
}
else if (cswitch == 2) {
/* a list of images - in diffim each one of these is associated with a basis function */
std::vector<std::shared_ptr<Eigen::VectorXd> > imageList(nParameters);
typename std::vector<std::shared_ptr<Eigen::VectorXd> >::iterator eiter = imageList.begin();
std::vector<Eigen::VectorXd> imageList(nParameters);
typename std::vector<Eigen::VectorXd>::iterator eiter = imageList.begin();
afwImage::Image<ImageT> cimage(varianceEstimate.getDimensions());
for (int i = 1; eiter != imageList.end(); ++eiter, ++i) {
cimage = i; /* give it a value */
Expand All @@ -85,7 +85,7 @@ Eigen::MatrixXd test(afwImage::Image<ImageT> varianceEstimate,
endRow-startRow,
endCol-startCol);
cMat.resize(cMat.rows()*cMat.cols(), 1);
std::shared_ptr<Eigen::VectorXd> vMat (new Eigen::VectorXd(cMat.col(0)));
Eigen::VectorXd vMat = cMat.col(0);
*eiter = vMat;
}

Expand All @@ -97,14 +97,14 @@ Eigen::MatrixXd test(afwImage::Image<ImageT> varianceEstimate,
eigeniVarianceM.resize(eigeniVarianceM.rows()*eigeniVarianceM.cols(), 1);
Eigen::VectorXd eigeniVarianceV = eigeniVarianceM.col(0);

typename std::vector<std::shared_ptr<Eigen::VectorXd> >::iterator eiteri = imageList.begin();
typename std::vector<std::shared_ptr<Eigen::VectorXd> >::iterator eiterE = imageList.end();
typename std::vector<Eigen::VectorXd>::iterator eiteri = imageList.begin();
typename std::vector<Eigen::VectorXd>::iterator eiterE = imageList.end();
for (unsigned int kidxi = 0; eiteri != eiterE; eiteri++, kidxi++) {
Eigen::VectorXd eiteriDotiVariance = ((*eiteri)->array() * eigeniVarianceV.array()).matrix();
Eigen::VectorXd eiteriDotiVariance = (eiteri->array() * eigeniVarianceV.array()).matrix();

typename std::vector<std::shared_ptr<Eigen::VectorXd> >::iterator eiterj = eiteri;
typename std::vector<Eigen::VectorXd>::iterator eiterj = eiteri;
for (unsigned int kidxj = kidxi; eiterj != eiterE; eiterj++, kidxj++) {
mMat(kidxi, kidxj) = (eiteriDotiVariance.array() * (**eiterj).array()).sum();
mMat(kidxi, kidxj) = (eiteriDotiVariance.array() * eiterj->array()).sum();
mMat(kidxj, kidxi) = mMat(kidxi, kidxj);
}
}
Expand Down
21 changes: 11 additions & 10 deletions include/lsst/ip/diffim/BasisLists.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,11 @@
* @ingroup ip_diffim
*/

#ifndef LSST_IP_DIFFIM_BASISSETS_H
#define LSST_IP_DIFFIM_BASISSETS_H
#ifndef LSST_IP_DIFFIM_BASISSLISTS_H
#define LSST_IP_DIFFIM_BASISSLISTS_H

#include <memory>
#include <vector>

#include "Eigen/Core"

Expand Down Expand Up @@ -48,7 +49,7 @@ namespace diffim {
* @note Calls either makeForwardDifferenceMatrix or
* makeCentralDifferenceMatrix based on the policy file.
*/
std::shared_ptr<Eigen::MatrixXd> makeRegularizationMatrix(
Eigen::MatrixXd makeRegularizationMatrix(
lsst::pex::policy::Policy policy
);

Expand All @@ -63,10 +64,10 @@ namespace diffim {
*
* @ingroup ip_diffim
*/
std::shared_ptr<Eigen::MatrixXd> makeForwardDifferenceMatrix(
Eigen::MatrixXd makeForwardDifferenceMatrix(
int width,
int height,
std::vector<int> const& orders,
std::vector<int> const & orders,
float borderPenalty,
bool fitForBackground
);
Expand All @@ -82,7 +83,7 @@ namespace diffim {
*
* @ingroup ip_diffim
*/
std::shared_ptr<Eigen::MatrixXd> makeCentralDifferenceMatrix(
Eigen::MatrixXd makeCentralDifferenceMatrix(
int width,
int height,
int stencil,
Expand Down Expand Up @@ -123,10 +124,10 @@ namespace diffim {
* @ingroup ip_diffim
*/
lsst::afw::math::KernelList makeAlardLuptonBasisList(
int halfWidth, ///< size is 2*N + 1
int nGauss, ///< number of gaussians
std::vector<double> const& sigGauss, ///< width of the gaussians
std::vector<int> const& degGauss ///< local spatial variation of gaussians
int halfWidth,
int nGauss,
std::vector<double> const& sigGauss,
std::vector<int> const& degGauss
);

}}} // end of namespace lsst::ip::diffim
Expand Down
6 changes: 3 additions & 3 deletions include/lsst/ip/diffim/BuildSingleKernelVisitor.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ namespace detail {
BuildSingleKernelVisitor(
lsst::afw::math::KernelList const& basisList,
lsst::pex::policy::Policy const& policy,
std::shared_ptr<Eigen::MatrixXd> hMat
Eigen::MatrixXd const& hMat
);
virtual ~BuildSingleKernelVisitor() {};

Expand All @@ -61,7 +61,7 @@ namespace detail {
private:
lsst::afw::math::KernelList const _basisList; ///< Basis set
lsst::pex::policy::Policy _policy; ///< Policy controlling behavior
std::shared_ptr<Eigen::MatrixXd> _hMat; ///< Regularization matrix
Eigen::MatrixXd const _hMat; ///< Regularization matrix
ImageStatistics<PixelT> _imstats; ///< To calculate statistics of difference image
bool _skipBuilt; ///< Skip over built candidates during processCandidate()
int _nRejected; ///< Number of candidates rejected during processCandidate()
Expand Down Expand Up @@ -89,7 +89,7 @@ namespace detail {
makeBuildSingleKernelVisitor(
lsst::afw::math::KernelList const& basisList,
lsst::pex::policy::Policy const& policy,
std::shared_ptr<Eigen::MatrixXd> hMat
Eigen::MatrixXd const & hMat
) {

return typename BuildSingleKernelVisitor<PixelT>::Ptr(
Expand Down
4 changes: 2 additions & 2 deletions include/lsst/ip/diffim/KernelCandidate.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,7 @@ namespace diffim {
);
void build(
afw::math::KernelList const& basisList,
std::shared_ptr<Eigen::MatrixXd> hMat
Eigen::MatrixXd const& hMat
);

private:
Expand All @@ -199,7 +199,7 @@ namespace diffim {
std::shared_ptr<StaticKernelSolution<PixelT> > _kernelSolutionPca; ///< Most recent solution

void _buildKernelSolution(afw::math::KernelList const& basisList,
std::shared_ptr<Eigen::MatrixXd> hMat);
Eigen::MatrixXd const& hMat);
};


Expand Down
42 changes: 21 additions & 21 deletions include/lsst/ip/diffim/KernelSolution.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,32 +46,32 @@ namespace diffim {
SVD = 1
};

explicit KernelSolution(std::shared_ptr<Eigen::MatrixXd> mMat,
std::shared_ptr<Eigen::VectorXd> bVec,
explicit KernelSolution(Eigen::MatrixXd mMat,
Eigen::VectorXd bVec,
bool fitForBackground);
explicit KernelSolution(bool fitForBackground);
explicit KernelSolution();

virtual ~KernelSolution() {};
virtual void solve();
virtual void solve(Eigen::MatrixXd mMat,
Eigen::VectorXd bVec);
virtual void solve(Eigen::MatrixXd const& mMat,
Eigen::VectorXd const& bVec);
KernelSolvedBy getSolvedBy() {return _solvedBy;}
virtual double getConditionNumber(ConditionNumberType conditionType);
virtual double getConditionNumber(Eigen::MatrixXd mMat, ConditionNumberType conditionType);
virtual double getConditionNumber(Eigen::MatrixXd const& mMat, ConditionNumberType conditionType);

inline std::shared_ptr<Eigen::MatrixXd> getM() {return _mMat;}
inline std::shared_ptr<Eigen::VectorXd> getB() {return _bVec;}
void printM() {std::cout << *_mMat << std::endl;}
void printB() {std::cout << *_bVec << std::endl;}
void printA() {std::cout << *_aVec << std::endl;}
inline Eigen::MatrixXd const& getM() {return _mMat;}
inline Eigen::VectorXd const& getB() {return _bVec;}
void printM() {std::cout << _mMat << std::endl;}
void printB() {std::cout << _bVec << std::endl;}
void printA() {std::cout << _aVec << std::endl;}
inline int getId() const { return _id; }

protected:
int _id; ///< Unique ID for object
std::shared_ptr<Eigen::MatrixXd> _mMat; ///< Derived least squares M matrix
std::shared_ptr<Eigen::VectorXd> _bVec; ///< Derived least squares B vector
std::shared_ptr<Eigen::VectorXd> _aVec; ///< Derived least squares solution matrix
Eigen::MatrixXd _mMat; ///< Derived least squares M matrix
Eigen::VectorXd _bVec; ///< Derived least squares B vector
Eigen::VectorXd _aVec; ///< Derived least squares solution matrix
KernelSolvedBy _solvedBy; ///< Type of algorithm used to make solution
bool _fitForBackground; ///< Background terms included in fit
static int _SolutionId; ///< Unique identifier for solution
Expand Down Expand Up @@ -101,9 +101,9 @@ namespace diffim {
virtual std::pair<std::shared_ptr<lsst::afw::math::Kernel>, double> getSolutionPair();

protected:
std::shared_ptr<Eigen::MatrixXd> _cMat; ///< K_i x R
std::shared_ptr<Eigen::VectorXd> _iVec; ///< Vectorized I
std::shared_ptr<Eigen::VectorXd> _ivVec; ///< Inverse variance
Eigen::MatrixXd _cMat; ///< K_i x R
Eigen::VectorXd _iVec; ///< Vectorized I
Eigen::VectorXd _ivVec; ///< Inverse variance

lsst::afw::math::Kernel::Ptr _kernel; ///< Derived single-object convolution kernel
double _background; ///< Derived differential background estimate
Expand Down Expand Up @@ -150,7 +150,7 @@ namespace diffim {

RegularizedKernelSolution(lsst::afw::math::KernelList const& basisList,
bool fitForBackground,
std::shared_ptr<Eigen::MatrixXd> hMat,
Eigen::MatrixXd const& hMat,
lsst::pex::policy::Policy policy
);
virtual ~RegularizedKernelSolution() {};
Expand All @@ -159,10 +159,10 @@ namespace diffim {
double estimateRisk(double maxCond);

/* Include additive term (_lambda * _hMat) in M matrix? */
std::shared_ptr<Eigen::MatrixXd> getM(bool includeHmat = true);
Eigen::MatrixXd getM(bool includeHmat = true);

private:
std::shared_ptr<Eigen::MatrixXd> _hMat; ///< Regularization weights
Eigen::MatrixXd const _hMat; ///< Regularization weights
double _lambda; ///< Overall regularization strength
lsst::pex::policy::Policy _policy;

Expand All @@ -184,8 +184,8 @@ namespace diffim {
virtual ~SpatialKernelSolution() {};

void addConstraint(float xCenter, float yCenter,
std::shared_ptr<Eigen::MatrixXd> qMat,
std::shared_ptr<Eigen::VectorXd> wVec);
Eigen::MatrixXd const& qMat,
Eigen::VectorXd const& wVec);

void solve();
lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel>::Ptr makeKernelImage(lsst::afw::geom::Point2D const& pos);
Expand Down
1 change: 1 addition & 0 deletions include/lsst/ip/diffim/KernelSumVisitor.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#define LSST_IP_DIFFIM_KERNELSUMVISITOR_H

#include <memory>

#include "lsst/afw/math.h"
#include "lsst/afw/image.h"
#include "lsst/pex/policy/Policy.h"
Expand Down
26 changes: 14 additions & 12 deletions python/lsst/ip/diffim/detail/buildSingleKernelVisitor.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,26 +54,27 @@ template <typename PixelT>
void declareBuildSingleKernelVisitor(py::module& mod, std::string const& suffix) {
py::class_<BuildSingleKernelVisitor<PixelT>, std::shared_ptr<BuildSingleKernelVisitor<PixelT>>,
afw::math::CandidateVisitor>
cls(mod, ("BuildSingleKernelVisitor" + suffix).c_str());
cls(mod, ("BuildSingleKernelVisitor" + suffix).c_str());

cls.def(py::init<afw::math::KernelList, pex::policy::Policy>(), "basisList"_a, "policy"_a);
cls.def(py::init<afw::math::KernelList, pex::policy::Policy, std::shared_ptr<Eigen::MatrixXd>>(),
"basisList"_a, "policy"_a, "hMat"_a);
cls.def(py::init<afw::math::KernelList, pex::policy::Policy, Eigen::MatrixXd const&>(), "basisList"_a,
"policy"_a, "hMat"_a);

cls.def("setSkipBuilt", &BuildSingleKernelVisitor<PixelT>::setSkipBuilt, "skip"_a);
cls.def("getNRejected", &BuildSingleKernelVisitor<PixelT>::getNRejected);
cls.def("getNProcessed", &BuildSingleKernelVisitor<PixelT>::getNProcessed);
cls.def("reset", &BuildSingleKernelVisitor<PixelT>::reset);
cls.def("processCandidate", &BuildSingleKernelVisitor<PixelT>::processCandidate, "candidate"_a);

mod.def("makeBuildSingleKernelVisitor", (std::shared_ptr<BuildSingleKernelVisitor<PixelT>>(*)(
afw::math::KernelList const&, pex::policy::Policy const&)) &
makeBuildSingleKernelVisitor<PixelT>,
mod.def("makeBuildSingleKernelVisitor",
(std::shared_ptr<BuildSingleKernelVisitor<PixelT>>(*)(afw::math::KernelList const&,
pex::policy::Policy const&)) &
makeBuildSingleKernelVisitor<PixelT>,
"basisList"_a, "policy"_a);
mod.def("makeBuildSingleKernelVisitor",
(std::shared_ptr<BuildSingleKernelVisitor<PixelT>>(*)(
afw::math::KernelList const&, pex::policy::Policy const&, std::shared_ptr<Eigen::MatrixXd>)) &
makeBuildSingleKernelVisitor<PixelT>,
afw::math::KernelList const&, pex::policy::Policy const&, Eigen::MatrixXd const&)) &
makeBuildSingleKernelVisitor<PixelT>,
"basisList"_a, "policy"_a, "hMat"_a);
}

Expand All @@ -92,7 +93,8 @@ PYBIND11_PLUGIN(_buildSingleKernelVisitor) {

return mod.ptr();
}
}
}
}
} // lsst::ip::diffim::detail

} // detail
} // diffim
} // ip
} // lsst

0 comments on commit 06a92e4

Please sign in to comment.