Modified C++ tool

1 files changed, 61 insertions, 58 deletions

diff --git a/fits.cpp b/fits.cpp
index 923712b..f2f3a1d 100644
--- a/fits.cpp
+++ b/fits.cpp
@@ -1,6 +1,6 @@
 #include <getopt.h>
 #include <iomanip>
-#include<list>
+#include <list>
 
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
@@ -12,80 +12,73 @@ using Rng = randutils::random_generator<pcg32>;
 using Real = double;
 using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
 
-class fitModel {
-public:
-  unsigned N;
+/*
+inline Real basisFunctions(unsigned i, Real x) const {
+  if (i == 0) {
+    return 1;
+  } else {
+    return pow(x, i);
+  }
+}
+*/
 
-  fitModel(unsigned N) : N(N) {}
+inline Real basisFunctions(unsigned N, unsigned i, Real x) {
+  return abs(x - i / (N - 1.0));
+}
 
-  /*
-  virtual Real basisFunctions(unsigned i, Real x) const {
-    if (i == 0) {
-      return 1;
-    } else {
-      return pow(x, i);
-    }
-  }
-  */
+Real value(const Vector& coeffs, Real x) {
+  Real v = 0;
 
-  Real basisFunctions(unsigned i, Real x) const {
-    return abs(x - i / (N - 1.0));
+  for (unsigned j = 0; j < coeffs.size(); j++) {
+    v += coeffs[j] * basisFunctions(coeffs.size(), j, x);
   }
 
-  Real value(const Vector& coeffs, Real x) const {
-    Real v = 0;
+  return v;
+}
 
-    for (unsigned j = 0; j < N; j++) {
-      v += coeffs[j] * basisFunctions(j, x);
-    }
+Real cost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) {
+  Real c = 0;
 
-    return v;
+  for (const std::tuple<Real, Real>& xy : data) {
+    Real x = std::get<0>(xy);
+    Real y = std::get<1>(xy);
+    c += 0.5 * pow(value(coeffs, x) - y, 2);
   }
 
-  Real cost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) const {
-    Real c = 0;
+  return c;
+}
 
-    for (std::tuple<Real, Real> xy : data) {
+Vector dCost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) {
+  Vector dC = Vector::Zero(coeffs.size());
+
+  for (unsigned j = 0; j < coeffs.size(); j++) {
+    for (const std::tuple<Real, Real>& xy : data) {
       Real x = std::get<0>(xy);
       Real y = std::get<1>(xy);
-      c += 0.5 * pow(value(coeffs, x) - y, 2);
+      dC[j] += (value(coeffs, x) - y) * basisFunctions(coeffs.size(), j, x);
     }
-
-    return c;
   }
 
-  Vector dCost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) const {
-    Vector dC = Vector::Zero(N);
-
-    for (unsigned j = 0; j < N; j++) {
-      for (std::tuple<Real, Real> xy : data) {
-        Real x = std::get<0>(xy);
-        Real y = std::get<1>(xy);
-        dC[j] += (value(coeffs, x) - y) * basisFunctions(j, x);
-      }
-    }
-
-    return dC;
-  }
-};
+  return dC;
+}
 
-Vector gradientDescent(const fitModel& M, const std::list<std::tuple<Real, Real>>& data, const Vector& a₀, unsigned maxSteps, Real ε = 1e-12) {
+Vector gradientDescent(const std::list<std::tuple<Real, Real>>& data, const Vector& a₀, unsigned maxSteps, Real ε = 1e-12) {
   Vector xₜ = a₀;
-  Real Hₜ = M.cost(data, a₀);
+  Real Hₜ = cost(data, a₀);
   Real α = 1.0;
   Real m;
   Vector ∇H;
   unsigned steps = 0;
 
   while (
-    ∇H = M.dCost(data, xₜ), m = ∇H.squaredNorm(),
+    ∇H = dCost(data, xₜ), m = ∇H.squaredNorm(),
     m > ε && steps < maxSteps
   ) {
     Vector xₜ₊₁;
     Real Hₜ₊₁;
 
     while (
-      xₜ₊₁ = xₜ - α * ∇H, Hₜ₊₁ = M.cost(data, xₜ₊₁),
+      xₜ₊₁ = xₜ - α * ∇H, Hₜ₊₁ = cost(data, xₜ₊₁),
       Hₜ₊₁ > Hₜ - 0.5 * α * m
     ) {
       α /= 2;
@@ -100,7 +93,23 @@ Vector gradientDescent(const fitModel& M, const std::list<std::tuple<Real, Real>
   return xₜ;
 }
 
-Vector stochasticGradientDescent(const fitModel& M, const std::list<std::tuple<Real, Real>>& data, const Vector& a₀, Rng& r, unsigned maxSteps, Real ε = 1e-12) {
+Vector dCostRand(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs, Real batchP, Rng& r) {
+  Vector dC = Vector::Zero(coeffs.size());
+
+  for (unsigned j = 0; j < coeffs.size(); j++) {
+    for (const std::tuple<Real, Real>& xy : data) {
+      if (batchP > r.uniform(0.0, 1.0)) {
+        Real x = std::get<0>(xy);
+        Real y = std::get<1>(xy);
+        dC[j] += (value(coeffs, x) - y) * basisFunctions(coeffs.size(), j, x);
+      }
+    }
+  }
+
+  return dC;
+}
+
+Vector stochasticGradientDescent(std::list<std::tuple<Real, Real>>& data, const Vector& a₀, Rng& r, unsigned maxSteps, Real ε = 1e-12) {
   Vector xₜ = a₀;
   Real α = 1e-3;
   Real m;
@@ -108,14 +117,10 @@ Vector stochasticGradientDescent(const fitModel& M, const std::list<std::tuple<R
   unsigned steps = 0;
 
   while (
-    ∇H = M.dCost(data, xₜ), m = ∇H.squaredNorm(),
+    ∇H = dCostRand(data, xₜ, 0.1, r), m = ∇H.squaredNorm(),
     m > ε && steps < maxSteps
   ) {
-    Vector γ(M.N);
-    for (Real& γi : γ) {
-      γi = r.variate<Real, std::normal_distribution>(0, 1e-4);
-    }
-    xₜ = xₜ - α * ∇H + γ;
+    xₜ = xₜ - α * ∇H;
     steps++;
   }
 
@@ -137,7 +142,7 @@ int main(int argc, char* argv[]) {
   unsigned N = 10;
   unsigned M = 10;
   Real σ = 0.2;
-  unsigned maxSteps = 1e5;
+  unsigned maxSteps = 1e8;
 
   int opt;
 
@@ -162,7 +167,7 @@ int main(int argc, char* argv[]) {
 
   Rng r;
 
-  std::list<std::tuple<Real, Real>> data = generateData([](Real x) {return sin(2 * M_PI * x);}, M, σ, r);
+  std::list<std::tuple<Real, Real>> data = generateData([](Real x) {return x;}, M, σ, r);
 
   std::cout << std::setprecision(15);
 
@@ -171,10 +176,8 @@ int main(int argc, char* argv[]) {
   }
   std::cout << std::endl;
 
-  fitModel model(N);
-
   Vector a₀ = Vector::Zero(N);
-  Vector a = gradientDescent(model, data, a₀, maxSteps);
+  Vector a = stochasticGradientDescent(data, a₀, r, maxSteps);
 
   for (unsigned i = 0; i < N; i++) {
     std::cout << a[i] << " ";