IMP logo
IMP Reference Guide  2.6.1
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-2016 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 <IMP/algebra/eigen3/Eigen/Dense>
15 
16 using namespace IMP_Eigen;
17 
18 IMPSAXS_BEGIN_NAMESPACE
19 
20 //! non-negative least square fitting for profile weight solving problem
21 inline IMP_Eigen::VectorXf NNLS(const IMP_Eigen::MatrixXf& A,
22  const IMP_Eigen::VectorXf& b) {
23 
24  // TODO: make JacobiSVD a class object to avoid memory re-allocations
25  IMP_Eigen::JacobiSVD<IMP_Eigen::MatrixXf> svd(A, ComputeThinU | ComputeThinV);
26  IMP_Eigen::VectorXf x = svd.solve(b);
27 
28  // compute a small negative tolerance
29  double tol = 0;
30  int n = A.cols();
31  int m = A.rows();
32 
33  // count initial negatives
34  int negs = 0;
35  for (int i = 0; i < n; i++)
36  if (x[i] < 0.0) negs++;
37  if (negs <= 0) return x;
38 
39  int sip = int(negs / 100);
40  if (sip < 1) sip = 1;
41 
42  IMP_Eigen::VectorXf zeroed = IMP_Eigen::VectorXf::Zero(n);
43  IMP_Eigen::MatrixXf C = A;
44 
45  // iteratively zero some x values
46  for (int count = 0; count < n; count++) { // loop till no negatives found
47  // count negatives and choose how many to treat
48  negs = 0;
49  for (int i = 0; i < n; i++)
50  if (zeroed[i] < 1.0 && x[i] < tol) negs++;
51  if (negs <= 0) break;
52 
53  int gulp = std::max(negs / 20, sip);
54 
55  // zero the most negative solution values
56  for (int k = 1; k <= gulp; k++) {
57  int p = -1;
58  double worst = 0.0;
59  for (int j = 0; j < n; j++)
60  if (zeroed[j] < 1.0 && x[j] < worst) {
61  p = j;
62  worst = x[p];
63  }
64  if (p < 0) break;
65  for (int i = 0; i < m; i++) C(i, p) = 0.0;
66  zeroed[p] = 9;
67  }
68 
69  // re-solve
70  IMP_Eigen::JacobiSVD<IMP_Eigen::MatrixXf> svd(C,
71  ComputeThinU | ComputeThinV);
72  x = svd.solve(b);
73  }
74 
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;
77 
78  return x;
79 }
80 
81 IMPSAXS_END_NAMESPACE
82 
83 #endif /* IMPSAXS_NNLS_H */
IMP_Eigen::VectorXf NNLS(const IMP_Eigen::MatrixXf &A, const IMP_Eigen::VectorXf &b)
non-negative least square fitting for profile weight solving problem
Definition: nnls.h:21