From d22c11c61f60affb54756f2dbb3267743f1d24a0 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 10 Feb 2025 14:37:30 -0300
Subject: Made some optimizations to the code.

---
 fits.cpp | 30 +++++++++++++++++-------------
 1 file changed, 17 insertions(+), 13 deletions(-)

(limited to 'fits.cpp')

diff --git a/fits.cpp b/fits.cpp
index 62d75c9..3843dcf 100644
--- a/fits.cpp
+++ b/fits.cpp
@@ -42,9 +42,10 @@ Real value(const Vector& coeffs, Real x) {
   return v;
 }
 
-Real cost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) {
+Real cost(const std::vector<std::tuple<Real, Real>>& data, const Vector& coeffs) {
   Real c = 0;
 
+#pragma omp parallel for reduction(+:c)
   for (const std::tuple<Real, Real>& xy : data) {
     Real x = std::get<0>(xy);
     Real y = std::get<1>(xy);
@@ -54,14 +55,16 @@ Real cost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) {
   return c;
 }
 
-Vector dCost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) {
+Vector dCost(const std::vector<std::tuple<Real, Real>>& data, const Vector& coeffs) {
   Vector dC = Vector::Zero(coeffs.size());
 
   for (const std::tuple<Real, Real>& xy : data) {
+    Real x = std::get<0>(xy);
+    Real y = std::get<1>(xy);
+    Real Δc = value(coeffs, x) - y;
+
+#pragma omp parallel for
     for (unsigned j = 0; j < coeffs.size(); j++) {
-      Real x = std::get<0>(xy);
-      Real y = std::get<1>(xy);
-      Real Δc = value(coeffs, x) - y;
       dC[j] += Δc * basisFunctions(coeffs.size(), j, x);
     }
   }
@@ -69,7 +72,7 @@ Vector dCost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs
   return dC;
 }
 
-Vector gradientDescent(const std::list<std::tuple<Real, Real>>& data, const Vector& a₀, unsigned maxSteps, Real ε = 1e-12) {
+Vector gradientDescent(const std::vector<std::tuple<Real, Real>>& data, const Vector& a₀, unsigned maxSteps, Real ε = 1e-12) {
   Vector xₜ = a₀;
   Real Hₜ = cost(data, a₀);
   Real α = 1.0;
@@ -100,7 +103,7 @@ Vector gradientDescent(const std::list<std::tuple<Real, Real>>& data, const Vect
   return xₜ;
 }
 
-Vector dCostRand(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs, Real batchP, Rng& r) {
+Vector dCostRand(const std::vector<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++) {
@@ -116,7 +119,7 @@ Vector dCostRand(const std::list<std::tuple<Real, Real>>& data, const Vector& co
   return dC;
 }
 
-Vector stochasticGradientDescent(std::list<std::tuple<Real, Real>>& data, const Vector& a₀, Rng& r, unsigned maxSteps, Real ε = 1e-12) {
+Vector stochasticGradientDescent(std::vector<std::tuple<Real, Real>>& data, const Vector& a₀, Rng& r, unsigned maxSteps, Real ε = 1e-12) {
   Vector xₜ = a₀;
   Real α = 1e-3;
   Real m;
@@ -134,8 +137,9 @@ Vector stochasticGradientDescent(std::list<std::tuple<Real, Real>>& data, const
   return xₜ;
 }
 
-std::list<std::tuple<Real, Real>> generateData(Real(*f)(Real), unsigned M, Real σ, Rng& r) {
-  std::list<std::tuple<Real, Real>> data;
+std::vector<std::tuple<Real, Real>> generateData(Real(*f)(Real), unsigned M, Real σ, Rng& r) {
+  std::vector<std::tuple<Real, Real>> data;
+  data.reserve(M);
 
   for (unsigned i = 0; i < M; i++) {
     Real x = ((Real)i) / (M - 1.0);
@@ -149,7 +153,7 @@ int main(int argc, char* argv[]) {
   unsigned N = 10;
   unsigned M = 10;
   Real σ = 0.2;
-  unsigned maxSteps = 1e8;
+  long unsigned maxSteps = 1e12;
 
   int opt;
 
@@ -165,7 +169,7 @@ int main(int argc, char* argv[]) {
       σ = atof(optarg);
       break;
     case 'S':
-      maxSteps = (unsigned)atof(optarg);
+      maxSteps = (long unsigned)atof(optarg);
       break;
     default:
       exit(1);
@@ -174,7 +178,7 @@ int main(int argc, char* argv[]) {
 
   Rng r;
 
-  std::list<std::tuple<Real, Real>> data = generateData([](Real x) {return std::cos(2 * M_PI * x);}, M, σ, r);
+  std::vector<std::tuple<Real, Real>> data = generateData([](Real x) {return std::cos(2 * M_PI * x);}, M, σ, r);
 
   std::cout << std::setprecision(15);
 
-- 
cgit v1.2.3-70-g09d2