10 #ifndef IMPSAXS_NNLS_H
11 #define IMPSAXS_NNLS_H
13 #include <IMP/saxs/saxs_config.h>
14 #include <IMP/algebra/eigen3/Eigen/Dense>
16 using namespace IMP_Eigen;
18 IMPSAXS_BEGIN_NAMESPACE
21 inline IMP_Eigen::VectorXf
NNLS(
const IMP_Eigen::MatrixXf& A,
22 const IMP_Eigen::VectorXf& b) {
25 IMP_Eigen::JacobiSVD<IMP_Eigen::MatrixXf> svd(A, ComputeThinU | ComputeThinV);
26 IMP_Eigen::VectorXf x = svd.solve(b);
35 for (
int i = 0; i < n; i++)
36 if (x[i] < 0.0) negs++;
37 if (negs <= 0)
return x;
39 int sip = int(negs / 100);
42 IMP_Eigen::VectorXf zeroed = IMP_Eigen::VectorXf::Zero(n);
43 IMP_Eigen::MatrixXf C = A;
46 for (
int count = 0; count < n; count++) {
49 for (
int i = 0; i < n; i++)
50 if (zeroed[i] < 1.0 && x[i] < tol) negs++;
53 int gulp = std::max(negs / 20, sip);
56 for (
int k = 1; k <= gulp; k++) {
59 for (
int j = 0; j < n; j++)
60 if (zeroed[j] < 1.0 && x[j] < worst) {
65 for (
int i = 0; i < m; i++) C(i, p) = 0.0;
70 IMP_Eigen::JacobiSVD<IMP_Eigen::MatrixXf> svd(C,
71 ComputeThinU | ComputeThinV);
75 for (
int j = 0; j < n; j++)
76 if (x[j] < 0.0 && std::abs(x[j]) <= std::abs(tol)) x[j] = 0.0;
IMP_Eigen::VectorXf NNLS(const IMP_Eigen::MatrixXf &A, const IMP_Eigen::VectorXf &b)
non-negative least square fitting for profile weight solving problem