Skip to content

Commit f50ad2f

Browse files
akenmorrisakenmorris
authored
Merge pull request #2445 from SCIInstitute/early_stop_with_mahalanobis
Python API for Morphological Deviation Score and improve precision matrix computation
2 parents 5c971ff + 24a3533 commit f50ad2f

File tree

4 files changed

+124
-111
lines changed

4 files changed

+124
-111
lines changed

‎Libs/Optimize/Function/EarlyStop/EarlyStopping.cpp‎

Lines changed: 3 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,6 @@ void EarlyStopping::reset() {
4242
ema_initialized_ = false;
4343
ema_diff_ = Eigen::VectorXd();
4444
last_checked_iter_ = -1;
45-
// stop_flag_ = false;
4645
stop_flag_.store(false);
4746
}
4847

@@ -81,16 +80,12 @@ void EarlyStopping::update(int iteration, const ParticleSystem* p) {
8180

8281
}
8382

84-
// bool EarlyStopping::ShouldStop() const {
85-
// return stop_flag_;
86-
// }
87-
8883
//---------------------------------------------------------------------------
8984
bool EarlyStopping::ShouldStop() const { return stop_flag_.load(); }
9085

9186
//---------------------------------------------------------------------------
9287
Eigen::VectorXd EarlyStopping::ComputeScore(const Eigen::MatrixXd& X) {
93-
return score_func_.GetMahalanobisDistance(X);
88+
return score_func_.GetMorphoDevScore(X);
9489
}
9590

9691
//---------------------------------------------------------------------------
@@ -202,8 +197,7 @@ bool EarlyStopping::SetControlShapes(const ParticleSystem* p) {
202197
Eigen::RowVectorXd shape_vector(VDimension * num_points);
203198

204199
for (int k = 0; k < num_points; ++k) {
205-
// PointType pt = p->GetTransformedPosition(dom, k);
206-
PointType pt = p->GetPositions(dom)->Get(k);
200+
PointType pt = p->GetTransformedPosition(k, dom);
207201

208202
shape_vector(VDimension * k + 0) = pt[0];
209203
shape_vector(VDimension * k + 1) = pt[1];
@@ -252,8 +246,7 @@ Eigen::MatrixXd EarlyStopping::GetTestShapes(const ParticleSystem* p) {
252246
Eigen::RowVectorXd shape_vector(VDimension * num_points);
253247

254248
for (int k = 0; k < num_points; ++k) {
255-
// PointType pt = p->GetTransformedPosition(dom, k);
256-
PointType pt = p->GetPositions(dom)->Get(k);
249+
PointType pt = p->GetTransformedPosition(k, dom);
257250
shape_vector(3 * k + 0) = pt[0];
258251
shape_vector(3 * k + 1) = pt[1];
259252
shape_vector(3 * k + 2) = pt[2];

‎Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.cpp‎

Lines changed: 50 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,5 @@
11
#include "MorphologicalDeviationScore.h"
22
#include <Logging.h>
3-
43
namespace shapeworks {
54

65
MorphologicalDeviationScore::MorphologicalDeviationScore() = default;
@@ -32,44 +31,31 @@ bool MorphologicalDeviationScore::FitPPCA(const Eigen::MatrixXd& X) {
3231
X_c, Eigen::ComputeThinU | Eigen::ComputeThinV);
3332
Eigen::VectorXd s = svd.singularValues(); // (r,)
3433
Eigen::MatrixXd V = svd.matrixV(); // (d x r)
35-
Eigen::VectorXd eigvals =
36-
(s.array().square()) / std::max(1, n - 1); // (r,)
37-
38-
double total_var = eigvals.sum();
39-
if (total_var <= 0.0) {
40-
throw std::runtime_error(
41-
"Total variance is non-positive. Cannot fit PCA on control shapes");
42-
}
43-
44-
Eigen::VectorXd cum_var = eigvals;
45-
for (int i = 1; i < cum_var.size(); ++i) {
46-
cum_var(i) += cum_var(i - 1);
47-
}
48-
cum_var /= total_var;
49-
50-
int q = 0;
51-
for (int i = 0; i < cum_var.size(); ++i) {
52-
if (cum_var(i) > 0.95) {
53-
q = i + 1;
54-
break;
55-
}
56-
}
57-
if (q == 0) q = cum_var.size();
58-
59-
eigvals_ = eigvals.head(q); // (q,)
60-
components_ = V.leftCols(q); // (d x q)
61-
62-
if (q < d) {
63-
double rem_sum = eigvals.tail(eigvals.size() - q).sum();
64-
noise_variance_ = rem_sum / (d - q);
65-
} else {
66-
noise_variance_ = 0.0;
67-
}
34+
Eigen::VectorXd eigvals = (s.array().square()) / (n - 1);
35+
36+
// save_matrix(X, std::string(DEBUG_FILES_PTH) + "/X.npy");
37+
// save_matrix(X_c, std::string(DEBUG_FILES_PTH) + "/X_c.npy");
38+
// save_matrix(V, std::string(DEBUG_FILES_PTH) + "/V.npy");
39+
// save_vector(eigvals, std::string(DEBUG_FILES_PTH) + "/eigvals.npy");
40+
41+
Eigen::VectorXd cumsum = eigvals;
42+
for (int i = 1; i < eigvals.size(); ++i) cumsum[i] += cumsum[i - 1];
43+
Eigen::VectorXd cumvar = cumsum / eigvals.sum();
44+
n_components_ = 0;
45+
while (n_components_ < cumvar.size() && cumvar[n_components_] <= retained_variance_ratio_)
46+
++n_components_;
47+
int q = n_components_;
48+
principal_components_variance_ = eigvals.head(q); // (q,)
49+
components_ = V.leftCols(q); // (d x q)
50+
// save_matrix(components_, std::string(DEBUG_FILES_PTH) + "/components.npy");
51+
// save_vector(principal_components_variance_,
52+
// std::string(DEBUG_FILES_PTH) + "/eigvals.npy");
53+
noise_variance_ = eigvals.tail(eigvals.size() - q).sum() / double(d - q);
6854
return true;
6955

7056
} catch (std::exception& e) {
7157
SW_ERROR(
72-
"Exception in SVD computation for early stopping score function {}",
58+
"SVD computation failed for MorphologicalDeviationScore: {}",
7359
e.what());
7460
return false;
7561
}
@@ -84,7 +70,7 @@ Eigen::MatrixXd MorphologicalDeviationScore::ComputeCovarianceMatrix() {
8470
return noise_variance_ * Eigen::MatrixXd::Identity(d, d);
8571
}
8672

87-
Eigen::MatrixXd diag_lambda = eigvals_.asDiagonal(); // (q x q)
73+
Eigen::MatrixXd diag_lambda = principal_components_variance_.asDiagonal(); // (q x q)
8874
Eigen::MatrixXd cov =
8975
components_ *
9076
(diag_lambda - noise_variance_ * Eigen::MatrixXd::Identity(q, q)) *
@@ -99,56 +85,52 @@ Eigen::MatrixXd MorphologicalDeviationScore::ComputePrecisionMatrix() {
9985
const int d = components_.rows();
10086
const int q = components_.cols();
10187

102-
if (q == 0) {
103-
if (noise_variance_ <= 0.0)
104-
throw std::runtime_error(
105-
"Noise variance is zero; precision undefined.");
106-
return (1.0 / noise_variance_) * Eigen::MatrixXd::Identity(d, d);
107-
}
88+
Eigen::MatrixXd A_inv = (1 / noise_variance_) * Eigen::MatrixXd::Identity(d, d);
10889

109-
Eigen::MatrixXd P = components_ * components_.transpose(); // (d x d)
110-
Eigen::MatrixXd term_principal = components_ *
111-
eigvals_.cwiseInverse().asDiagonal() *
112-
components_.transpose(); // (d x d)
90+
Eigen::MatrixXd U = components_;
91+
Eigen::MatrixXd C_inv = principal_components_variance_.cwiseInverse().asDiagonal();
92+
Eigen::MatrixXd V = components_.transpose();
93+
// Woodbury matrix identity
94+
// (A + U C V)^(-1)
95+
// = A^(-1)
96+
// - A^(-1) * U * (C^(-1) + V * A^(-1) * U)^(-1) * V * A^(-1)
11397

114-
if (noise_variance_ <= 0.0) {
115-
if (q != d)
116-
throw std::runtime_error(
117-
"Noise variance is zero and q < d; covariance is singular.");
118-
return term_principal;
119-
}
98+
Eigen::MatrixXd precision = A_inv - A_inv * U * (C_inv + V * A_inv * U).inverse() * V * A_inv;
99+
// save_matrix(precision, std::string(DEBUG_FILES_PTH)+"/precision.npy");
120100

121-
Eigen::MatrixXd precision =
122-
term_principal + (1.0 / noise_variance_) *
123-
(Eigen::MatrixXd::Identity(d, d) - P); // (d x d)
124101
return precision;
125-
126102
} catch (std::exception& e) {
127103
SW_ERROR(
128-
"Exception in computation of precision matrix in early stopping score "
129-
"function {}",
104+
"Failed to compute precision matrix in the early stopping score function: {}",
130105
e.what());
131106
return Eigen::MatrixXd();
132107
}
133108
}
134109

135110
//---------------------------------------------------------------------------
136-
Eigen::VectorXd MorphologicalDeviationScore::GetMahalanobisDistance(
111+
Eigen::VectorXd MorphologicalDeviationScore::GetMorphoDevScore(
137112
const Eigen::MatrixXd& X) {
138113
try {
114+
139115
if (!is_fitted_) {
140116
throw std::runtime_error(
141-
"PCA model has not been fitted yet on control shapes.");
117+
"PPCA model is not fitted on control shapes.");
142118
}
143119
const int n = X.rows();
144-
Eigen::MatrixXd X_c = X.rowwise() - mean_;
145-
Eigen::VectorXd dist(n);
146-
for (int i = 0; i < n; ++i) {
147-
Eigen::RowVectorXd xi = X_c.row(i);
148-
dist(i) = std::sqrt(
149-
std::max(0.0, (xi * precision_matrix_ * xi.transpose())(0)));
150-
}
151-
return dist;
120+
// Eigen::MatrixXd X_c = X.rowwise() - mean_;
121+
// Eigen::VectorXd dist(n);
122+
// for (int i = 0; i < n; ++i) {
123+
// Eigen::RowVectorXd xi = X_c.row(i);
124+
// dist(i) = std::sqrt(
125+
// std::max(0.0, (xi * precision_matrix_ * xi.transpose())(0)));
126+
// }
127+
128+
Eigen::MatrixXd X_bar = X.rowwise() - mean_;
129+
Eigen::MatrixXd Y = X_bar * precision_matrix_; // (n_samples × n_features)
130+
Eigen::VectorXd sq_mahal = (X_bar.array() * Y.array()).rowwise().sum(); // quadratic forms
131+
Eigen::VectorXd mahalanobis = sq_mahal.array().sqrt(); // final distances
132+
return mahalanobis;
133+
152134
} catch (std::exception& e) {
153135
SW_ERROR(
154136
"Exception in computing Mahalanobis distance for early stopping score "
@@ -158,15 +140,4 @@ Eigen::VectorXd MorphologicalDeviationScore::GetMahalanobisDistance(
158140
}
159141
}
160142

161-
//---------------------------------------------------------------------------
162-
// double MorphologicalDeviationScore::GetDeviationScore(const Eigen::MatrixXd&
163-
// X_test) {
164-
// if (!is_fitted_) {
165-
// throw std::runtime_error("PPCA model has not been fitted. Call
166-
// SetControlShapes() first.");
167-
// }
168-
// Eigen::VectorXd dists = GetMahalanobisDistance(X_test);
169-
// return dists(0);
170-
// }
171-
172143
} // namespace shapeworks

‎Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.h‎

Lines changed: 36 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,30 +1,45 @@
11
#pragma once
22

33
#include <Eigen/Dense>
4+
// #include "cnpy.h"
45

56
namespace shapeworks {
67
class MorphologicalDeviationScore {
7-
public:
8-
MorphologicalDeviationScore();
9-
// Fit PPCA model on control shapes
10-
bool SetControlShapes(const Eigen::MatrixXd& X); // (n_samples x n_features)
11-
/// Get Mahalanobis-based deviation score for test samples (non-fixed shapes/domains)
12-
Eigen::VectorXd GetMahalanobisDistance(const Eigen::MatrixXd& X); // (n,)
8+
public:
9+
MorphologicalDeviationScore();
10+
// Fit PPCA model on control shapes
11+
bool SetControlShapes(const Eigen::MatrixXd& X); // (n_samples x n_features)
12+
/// Get Mahalanobis-based deviation score for test samples (non-fixed
13+
/// shapes/domains)
14+
Eigen::VectorXd GetMorphoDevScore(const Eigen::MatrixXd& X); // (n,)
1315

14-
private:
15-
/// Flag to ensure control shapes are set and PCA model is in place
16-
bool is_fitted_ = false;
16+
private:
17+
/// Flag to ensure control shapes are set and PCA model is in place
18+
bool is_fitted_ = false;
1719

18-
// Fitted model parameters
19-
Eigen::RowVectorXd mean_; // (1 x d)
20-
Eigen::MatrixXd components_; // (d x q)
21-
Eigen::VectorXd eigvals_; // (q,)
22-
double noise_variance_;
23-
// Derived matrices
24-
Eigen::MatrixXd precision_matrix_; // (d x d)
25-
// Helper functions
26-
bool FitPPCA(const Eigen::MatrixXd& X);
27-
Eigen::MatrixXd ComputeCovarianceMatrix();
28-
Eigen::MatrixXd ComputePrecisionMatrix();
20+
// Fitted model parameters
21+
Eigen::RowVectorXd mean_; // (1 x d)
22+
Eigen::MatrixXd components_; // (d x q)
23+
Eigen::VectorXd principal_components_variance_; // (q,)
24+
int n_components_ = 0;
25+
double noise_variance_ = 0.0;
26+
double retained_variance_ratio_ = 0.95;
27+
// Derived matrices
28+
Eigen::MatrixXd precision_matrix_; // (d x d)
29+
// Helper functions
30+
bool FitPPCA(const Eigen::MatrixXd& X);
31+
Eigen::MatrixXd ComputeCovarianceMatrix();
32+
Eigen::MatrixXd ComputePrecisionMatrix();
33+
// inline static void save_vector(const Eigen::VectorXd& v,
34+
// const std::string& fname) {
35+
// cnpy::npy_save(fname, v.data(), {(size_t)v.size()}, "w");
36+
// };
37+
// inline static void save_matrix(const Eigen::MatrixXd m,
38+
// const std::string& fname) {
39+
// Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
40+
// X_rm = m;
41+
// cnpy::npy_save(fname, X_rm.data(), {(size_t)m.rows(), (size_t)m.cols()},
42+
// "w");
43+
// };
2944
};
30-
} // namespace shapeworks
45+
} // namespace shapeworks

‎Libs/Python/ShapeworksPython.cpp‎

Lines changed: 35 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@ using namespace pybind11::literals;
4545
#include "Variant.h"
4646
#include "VectorImage.h"
4747
#include "pybind_utils.h"
48-
48+
#include <Optimize/Function/EarlyStop/MorphologicalDeviationScore.h>
4949
namespace fs = boost::filesystem;
5050

5151
using namespace shapeworks;
@@ -1701,4 +1701,38 @@ PYBIND11_MODULE(shapeworks_py, m) {
17011701
"Return the variant string content")
17021702

17031703
; // Variant
1704+
1705+
py::class_<MorphologicalDeviationScore>(m, "MorphologicalDeviationScore")
1706+
.def(py::init<>(), "Create an instance for MorphologicalDeviationScore")
1707+
.def("SetControlShapes", &MorphologicalDeviationScore::SetControlShapes,
1708+
py::arg("X"),
1709+
R"pbdoc(
1710+
Fit PPCA model on control shapes.
1711+
1712+
Parameters
1713+
----------
1714+
X : numpy.ndarray
1715+
Matrix of control shapes (n_samples x n_features)
1716+
1717+
Returns
1718+
-------
1719+
bool
1720+
True if fitting was successful, False otherwise.
1721+
)pbdoc")
1722+
1723+
.def("GetMorphoDevScore",
1724+
&MorphologicalDeviationScore::GetMorphoDevScore, py::arg("X"),
1725+
R"pbdoc(
1726+
Compute Mahalanobis-based deviation score for test samples.
1727+
1728+
Parameters
1729+
----------
1730+
X : numpy.ndarray
1731+
Matrix of test samples (n_samples x n_features)
1732+
1733+
Returns
1734+
-------
1735+
numpy.ndarray
1736+
Vector of Mahalanobis distances for each sample.
1737+
)pbdoc");
17041738
} // PYBIND11_MODULE(shapeworks_py)

0 commit comments

Comments
 (0)