10 #ifndef IMPSAXS_NNLS_H
11 #define IMPSAXS_NNLS_H
13 #include <IMP/saxs/saxs_config.h>
14 #include <Eigen/Dense>
16 IMPSAXS_BEGIN_NAMESPACE
19 inline Eigen::VectorXf
NNLS(
const Eigen::MatrixXf& A,
20 const Eigen::VectorXf& b) {
23 Eigen::JacobiSVD<Eigen::MatrixXf> svd(A, Eigen::ComputeThinU
24 | Eigen::ComputeThinV);
25 Eigen::VectorXf x = svd.solve(b);
34 for (
int i = 0; i < n; i++)
35 if (x[i] < 0.0) negs++;
36 if (negs <= 0)
return x;
38 int sip = int(negs / 100);
41 Eigen::VectorXf zeroed = Eigen::VectorXf::Zero(n);
42 Eigen::MatrixXf C = A;
45 for (
int count = 0; count < n; count++) {
48 for (
int i = 0; i < n; i++)
49 if (zeroed[i] < 1.0 && x[i] < tol) negs++;
52 int gulp = std::max(negs / 20, sip);
55 for (
int k = 1; k <= gulp; k++) {
58 for (
int j = 0; j < n; j++)
59 if (zeroed[j] < 1.0 && x[j] < worst) {
64 for (
int i = 0; i < m; i++) C(i, p) = 0.0;
69 Eigen::JacobiSVD<Eigen::MatrixXf> svd(C, Eigen::ComputeThinU
70 | Eigen::ComputeThinV);
74 for (
int j = 0; j < n; j++)
75 if (x[j] < 0.0 && std::abs(x[j]) <= std::abs(tol)) x[j] = 0.0;
Eigen::VectorXf NNLS(const Eigen::MatrixXf &A, const Eigen::VectorXf &b)
non-negative least square fitting for profile weight solving problem