Changed code for making presentation figures.

1 files changed, 32 insertions, 2 deletions

diff --git a/fits.cpp b/fits.cpp
index 9696e16..09f9be5 100644
--- a/fits.cpp
+++ b/fits.cpp
@@ -10,6 +10,7 @@ using Rng = randutils::random_generator<pcg32>;
 
 using Real = double;
 using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
+using Matrix = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
 using Data = std::vector<std::tuple<Real, Real>>;
 
 #ifndef FITS_ABS_BASIS
@@ -64,6 +65,22 @@ Vector dCost(Data::const_iterator batchBegin, Data::const_iterator batchEnd, con
   return dC;
 }
 
+Vector underSolve(const Data& data, unsigned N) {
+  unsigned M = data.size();
+  Matrix A(M, N);
+  Vector b(M);
+  for (unsigned i = 0; i < M; i++) {
+    Real x = std::get<0>(data[i]);
+    Real y = std::get<1>(data[i]);
+    b[i] = y;
+    for (unsigned j = 0; j < N; j++) {
+      A(i, j) = basisFunctions(N, j, x);
+    }
+  }
+
+  return A.colPivHouseholderQr().solve(b);
+}
+
 Vector stochasticGradientDescent(const Data& data, const Vector& a₀, unsigned nBatches, long unsigned maxSteps, Real ε = 1e-12) {
   Vector xₜ = a₀;
   Real Hₜ;
@@ -125,6 +142,7 @@ Data generateData(Real(*f)(Real), unsigned M, Real σ, Rng& r) {
 }
 
 int main(int argc, char* argv[]) {
+  unsigned Nend = 1;
   unsigned M = 10;
   unsigned nBatches = 5;
   Real σ = 0.2;
@@ -133,8 +151,11 @@ int main(int argc, char* argv[]) {
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "M:s:S:B:i:")) != -1) {
+  while ((opt = getopt(argc, argv, "M:s:S:B:i:N:")) != -1) {
     switch (opt) {
+    case 'N':
+      Nend = (unsigned)atof(optarg);
+      break;
     case 'M':
       M = (unsigned)atof(optarg);
       break;
@@ -166,7 +187,16 @@ int main(int argc, char* argv[]) {
   }
   std::cout << std::endl;
 
-  for (unsigned N = 1; N <= 2 * M; N++) {
+  for (unsigned N = 1; N <= M; N++) {
+    Vector a = underSolve(data, N);
+
+    for (Real ai : a) {
+      std::cout << ai << " ";
+    }
+    std::cout << std::endl;
+  }
+
+  for (unsigned N = Nend; N > M; N--) {
     Vector a₀ = Vector::Zero(N);
     for (Real& aa : a₀) {
       aa = r.variate<Real, std::normal_distribution>(0, iniVar);