IMP logo
IMP Reference Guide  2.16.0
The Integrative Modeling Platform
nnls.h
Go to the documentation of this file.
1 /**
2  * \file IMP/saxs/nnls.h
3  * \brief
4  *
5  * \authors Dina Schneidman
6  * Copyright 2007-2021 IMP Inventors. All rights reserved.
7  *
8  */
9 
10 #ifndef IMPSAXS_NNLS_H
11 #define IMPSAXS_NNLS_H
12 
13 #include <IMP/saxs/saxs_config.h>
14 #include <Eigen/Dense>
15 
16 IMPSAXS_BEGIN_NAMESPACE
17 
18 //! non-negative least square fitting for profile weight solving problem
19 inline Eigen::VectorXf NNLS(const Eigen::MatrixXf& A,
20  const Eigen::VectorXf& b) {
21 
22  // TODO: make JacobiSVD a class object to avoid memory re-allocations
23  Eigen::JacobiSVD<Eigen::MatrixXf> svd(A, Eigen::ComputeThinU
24  | Eigen::ComputeThinV);
25  Eigen::VectorXf x = svd.solve(b);
26 
27  // compute a small negative tolerance
28  double tol = 0;
29  int n = A.cols();
30  int m = A.rows();
31 
32  // count initial negatives
33  int negs = 0;
34  for (int i = 0; i < n; i++)
35  if (x[i] < 0.0) negs++;
36  if (negs <= 0) return x;
37 
38  int sip = int(negs / 100);
39  if (sip < 1) sip = 1;
40 
41  Eigen::VectorXf zeroed = Eigen::VectorXf::Zero(n);
42  Eigen::MatrixXf C = A;
43 
44  // iteratively zero some x values
45  for (int count = 0; count < n; count++) { // loop till no negatives found
46  // count negatives and choose how many to treat
47  negs = 0;
48  for (int i = 0; i < n; i++)
49  if (zeroed[i] < 1.0 && x[i] < tol) negs++;
50  if (negs <= 0) break;
51 
52  int gulp = std::max(negs / 20, sip);
53 
54  // zero the most negative solution values
55  for (int k = 1; k <= gulp; k++) {
56  int p = -1;
57  double worst = 0.0;
58  for (int j = 0; j < n; j++)
59  if (zeroed[j] < 1.0 && x[j] < worst) {
60  p = j;
61  worst = x[p];
62  }
63  if (p < 0) break;
64  for (int i = 0; i < m; i++) C(i, p) = 0.0;
65  zeroed[p] = 9;
66  }
67 
68  // re-solve
69  Eigen::JacobiSVD<Eigen::MatrixXf> svd(C, Eigen::ComputeThinU
70  | Eigen::ComputeThinV);
71  x = svd.solve(b);
72  }
73 
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;
76 
77  return x;
78 }
79 
80 IMPSAXS_END_NAMESPACE
81 
82 #endif /* IMPSAXS_NNLS_H */
Eigen::VectorXf NNLS(const Eigen::MatrixXf &A, const Eigen::VectorXf &b)
non-negative least square fitting for profile weight solving problem
Definition: nnls.h:19