diff options
Diffstat (limited to 'fits.cpp')
| -rw-r--r-- | fits.cpp | 114 |
1 files changed, 49 insertions, 65 deletions
@@ -1,6 +1,5 @@ #include <getopt.h> #include <iomanip> -#include <list> #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" @@ -11,26 +10,19 @@ using Rng = randutils::random_generator<pcg32>; using Real = double; using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>; +using Data = std::vector<std::tuple<Real, Real>>; -/* -inline Real basisFunctions(unsigned N, unsigned i, Real x) { - if (i == 0) { - return 1; - } else { - return pow(x, i); - } -} -*/ - +#ifndef FITS_ABS_BASIS inline Real basisFunctions(unsigned N, unsigned i, Real x) { return std::legendre(i, 2.0 * (x - 0.5)); } +#endif -/* +#ifdef FITS_ABS_BASIS inline Real basisFunctions(unsigned N, unsigned i, Real x) { return abs(x - i / (N - 1.0)); } -*/ +#endif Real value(const Vector& coeffs, Real x) { Real v = 0; @@ -42,25 +34,25 @@ Real value(const Vector& coeffs, Real x) { return v; } -Real cost(const std::vector<std::tuple<Real, Real>>& data, const Vector& coeffs) { +Real cost(Data::const_iterator batchBegin, Data::const_iterator batchEnd, 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); + for (Data::const_iterator it = batchBegin; it != batchEnd; it++) { + Real x = std::get<0>(*it); + Real y = std::get<1>(*it); c += 0.5 * pow(value(coeffs, x) - y, 2); } return c; } -Vector dCost(const std::vector<std::tuple<Real, Real>>& data, const Vector& coeffs) { +Vector dCost(Data::const_iterator batchBegin, Data::const_iterator batchEnd, 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); + for (Data::const_iterator it = batchBegin; it != batchEnd; it++) { + Real x = std::get<0>(*it); + Real y = std::get<1>(*it); Real Δc = value(coeffs, x) - y; #pragma omp parallel for @@ -72,73 +64,56 @@ Vector dCost(const std::vector<std::tuple<Real, Real>>& data, const Vector& coef return dC; } -Vector gradientDescent(const std::vector<std::tuple<Real, Real>>& data, const Vector& a₀, unsigned maxSteps, Real ε = 1e-12) { +Vector stochasticGradientDescent(const Data& data, const Vector& a₀, unsigned nBatches, long unsigned maxSteps, Real ε = 1e-12) { Vector xₜ = a₀; - Real Hₜ = cost(data, a₀); + Real Hₜ; Real α = 1.0; Real m; Vector ∇H; - unsigned steps = 0; + long unsigned steps = 0; + + unsigned batchSize = data.size() / nBatches; + Data::const_iterator batchBegin = data.begin(); + Data::const_iterator batchEnd = data.begin() + batchSize; + Data::const_iterator effectiveEnd = data.begin() + batchSize * nBatches; while ( - ∇H = dCost(data, xₜ), m = ∇H.squaredNorm(), + Hₜ = cost(batchBegin, batchEnd, xₜ), + ∇H = dCost(batchBegin, batchEnd, xₜ), + m = ∇H.squaredNorm(), m > ε && steps < maxSteps ) { Vector xₜ₊₁; Real Hₜ₊₁; while ( - xₜ₊₁ = xₜ - α * ∇H, Hₜ₊₁ = cost(data, xₜ₊₁), + xₜ₊₁ = xₜ - α * ∇H, Hₜ₊₁ = cost(batchBegin, batchEnd, xₜ₊₁), Hₜ₊₁ > Hₜ - 0.25 * α * m ) { α /= 2; } xₜ = xₜ₊₁; - Hₜ = Hₜ₊₁; α *= 1.25; steps++; - } - - return xₜ; -} - -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++) { - 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); - } + if (batchEnd == data.end()) { + batchBegin = data.begin(); + batchEnd = data.begin() + batchSize; + } else if (batchEnd == effectiveEnd) { + batchBegin = effectiveEnd; + batchEnd = data.end(); + } else { + batchBegin += batchSize; + batchEnd += batchSize; } } - return dC; -} - -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; - Vector ∇H; - unsigned steps = 0; - - while ( - ∇H = dCostRand(data, xₜ, 0.1, r), m = ∇H.squaredNorm(), - m > ε && steps < maxSteps - ) { - xₜ = xₜ - α * ∇H; - steps++; - } - return xₜ; } -std::vector<std::tuple<Real, Real>> generateData(Real(*f)(Real), unsigned M, Real σ, Rng& r) { - std::vector<std::tuple<Real, Real>> data; +Data generateData(Real(*f)(Real), unsigned M, Real σ, Rng& r) { + Data data; data.reserve(M); for (unsigned i = 0; i < M; i++) { @@ -152,12 +127,14 @@ std::vector<std::tuple<Real, Real>> generateData(Real(*f)(Real), unsigned M, Rea int main(int argc, char* argv[]) { unsigned N = 10; unsigned M = 10; + unsigned nBatches = 5; Real σ = 0.2; + Real iniVar = 0.0; long unsigned maxSteps = 1e12; int opt; - while ((opt = getopt(argc, argv, "N:M:s:S:")) != -1) { + while ((opt = getopt(argc, argv, "N:M:s:S:B:i:")) != -1) { switch (opt) { case 'N': N = (unsigned)atof(optarg); @@ -171,6 +148,12 @@ int main(int argc, char* argv[]) { case 'S': maxSteps = (long unsigned)atof(optarg); break; + case 'B': + nBatches = (unsigned)atof(optarg); + break; + case 'i': + iniVar = atof(optarg); + break; default: exit(1); } @@ -178,7 +161,7 @@ int main(int argc, char* argv[]) { Rng r; - std::vector<std::tuple<Real, Real>> data = generateData([](Real x) {return std::cos(2 * M_PI * x);}, M, σ, r); + Data data = generateData([](Real x) {return std::sin(2 * M_PI * x);}, M, σ, r); std::cout << std::setprecision(15); @@ -189,9 +172,10 @@ int main(int argc, char* argv[]) { Vector a₀ = Vector::Zero(N); for (Real& aa : a₀) { - aa = r.variate<Real, std::normal_distribution>(0, 0); + aa = r.variate<Real, std::normal_distribution>(0, iniVar); } - Vector a = gradientDescent(data, a₀, maxSteps); + + Vector a = stochasticGradientDescent(data, a₀, nBatches, maxSteps); for (unsigned i = 0; i < N; i++) { std::cout << a[i] << " "; |