rework

lib/importer.cpp

@@ -0,0 +1,73 @@
+#include "autoopt/importer.hpp"
+
+namespace autoopt {
+std::generator<std::pair<double, double>> parse_dat_file(
+    const std::filesystem::path& filename) {
+  // check if file exists
+  if (!std::filesystem::exists(filename)) {
+    throw std::runtime_error("File does not exist: " + filename.string());
+  }
+
+  std::fstream file(filename, std::ios::in);
+
+  std::string line;
+  while (std::getline(file, line)) {
+    if (line.empty() || line == "\r" || line[0] == '#') {
+      continue;
+    }
+    std::istringstream iss(line);
+    double x, y;
+    iss >> x;
+    for (size_t i = 1; i < 5; ++i) {
+      iss >> y;  // skip unused columns
+    }
+
+    co_yield {x, -y};
+  }
+  co_return;
+}
+
+std::pair<std::vector<double>, std::vector<double>> parse_row(auto& iter,
+                                                              auto sentinel) {
+  std::vector<double> xs;
+  std::vector<double> ys;
+  xs.push_back((*iter).first);
+  ys.push_back((*iter).second);
+  ++iter;
+  double last_x = xs[0];
+  bool forward = last_x < (*iter).first;
+  while (iter != sentinel && ((forward && last_x < (*iter).first) ||
+                              (!forward && last_x > (*iter).first))) {
+    xs.push_back((*iter).first);
+    ys.push_back((*iter).second);
+    last_x = (*iter).first;
+    ++iter;
+  }
+  if (!forward) {
+    std::reverse(xs.begin(), xs.end());
+    std::reverse(ys.begin(), ys.end());
+  }
+  return std::make_pair(std::move(xs), std::move(ys));
+}
+
+raw_data read_data(const std::filesystem::path &filename){
+  auto lines = parse_dat_file(filename);
+
+  raw_data result;
+
+  auto it = lines.begin();
+  while (it != lines.end()) {
+    auto [xs, ys] = parse_row(it, lines.end());
+    if (result._xs.rows() == 0) {
+      result._xs = Eigen::VectorXd(xs.size());
+      std::copy(xs.begin(), xs.end(), result._xs.begin());
+    }
+    Eigen::VectorXd ys_vec = Eigen::VectorXd(ys.size());
+    std::copy(ys.begin(), ys.end(), ys_vec.begin());
+    result._yss.push_back(ys_vec);
+  }
+
+  return result;
+}
+
+}  // namespace autoopt

lib/interface.cpp

@@ -1,89 +0,0 @@
-#include <autoopt/btls.hpp>
-#include <autoopt/ellipse.hpp>
-#include <autoopt/parabola.hpp>
-#include <autoopt/interface.hpp>
-#include <autoopt/optimization_problem.hpp>
-#include <iomanip>
-#include <iostream>
-
-Eigen::VectorX<double> autoopt::fit_ellipse(
-    const std::vector<std::pair<double, double>>& data,
-    const Eigen::VectorX<double>& initial_params,
-    const Eigen::VectorX<double>& delta) {
-  auto opt_func = [&]<typename T>(const Eigen::VectorX<T>& params) {
-    ellipse<T> ellip(params(0), params(1), params(2));
-    quadric<T> quad = ellip.to_quadric().rotated_by(params(3));
-    T error = T(0);
-    for (const auto& [x, y] : data) {
-      T y_fit = quad.slope_at(T(x));
-      error = error + (y_fit - T(y)) * (y_fit - T(y));
-    }
-    return error / T(data.size());
-  };
-
-  auto_diff_optimization_problem problem(opt_func, initial_params);
-
-  log_barrier_optimization_problem lb_problem(problem, delta, 1e-4);
-
-  while (lb_problem._barrier_strength > 1e-20) {
-    btls(lb_problem);
-    lb_problem._barrier_strength *= 1e-2;
-  }
-
-  Eigen::VectorX<double> fitted_params = problem.x();
-
-  std::cout << "Fitted parameters: " << std::setprecision(10)
-            << fitted_params.transpose() << std::endl;
-
-  double obj_value = problem.objective(fitted_params);
-  std::cout << "Objective value: " << obj_value << std::endl;
-
-  // rms in radians
-  std::cout << "RMS error: " << std::sqrt(obj_value) << " radians" << std::endl;
-  // rms in arcsec
-  double rms_arcsec = std::sqrt(obj_value) * (3600.0 * 180.0 / M_PI);
-  std::cout << "RMS error: " << rms_arcsec << " arcsec" << std::endl;
-
-  return fitted_params;
-}
-
-Eigen::VectorX<double> autoopt::fit_parabola(
-    const std::vector<std::pair<double, double>>& data,
-    const Eigen::VectorX<double>& initial_params,
-    const Eigen::VectorX<double>& delta) {
-  auto opt_func = [&]<typename T>(const Eigen::VectorX<T>& params) {
-    parabola<T> parab(params(0), params(1));
-    quadric<T> quad = parab.to_quadric().rotated_by(params(2));
-    T error = T(0);
-    for (const auto& [x, y] : data) {
-      T y_fit = quad.slope_at(T(x));
-      error = error + (y_fit - T(y)) * (y_fit - T(y));
-    }
-    return error / T(data.size());
-  };
-
-  auto_diff_optimization_problem problem(opt_func, initial_params);
-
-  log_barrier_optimization_problem lb_problem(problem, delta, 1e-4);
-
-  while (lb_problem._barrier_strength > 1e-20) {
-    btls(lb_problem);
-    lb_problem._barrier_strength *= 1e-2;
-  }
-
-  Eigen::VectorX<double> fitted_params = problem.x();
-
-  std::cout << "Fitted parameters: " << std::setprecision(10)
-            << fitted_params.transpose() << std::endl;
-
-  double obj_value = problem.objective(fitted_params);
-  std::cout << "Objective value: " << obj_value << std::endl;
-
-  // rms in radians
-  std::cout << "RMS error: " << std::sqrt(obj_value) << " radians" << std::endl;
-  // rms in arcsec
-  double rms_arcsec = std::sqrt(obj_value) * (3600.0 * 180.0 / M_PI);
-  std::cout << "RMS error: " << rms_arcsec << " arcsec" << std::endl;
-
-  return fitted_params;
-}

lib/optimizer.cpp

@@ -0,0 +1,224 @@
+#include "autoopt/optimizer.hpp"
+
+#include <Eigen/src/Core/Matrix.h>
+#include <fftw3.h>
+
+#include <algorithm>
+#include <cstddef>
+#include <iterator>
+#include <print>
+#include <stdexcept>
+#include <tuple>
+#include <variant>
+
+#include "autoopt/ellipse_new.hpp"
+#include "autoopt/parabola_new.hpp"
+
+namespace autoopt {
+
+void optimizer::import_file(std::filesystem::path file_path) {
+  _data = data{};
+  _data->_raw = read_data(file_path);
+  _data->lb = *_data->_raw._xs.begin();
+  _data->ub = *(_data->_raw._xs.end() - 1);
+  plot_input_cb();
+  init_bounds_cb();
+  update_test_data();
+  update_table_cb();
+}
+
+void optimizer::set_lower_bound(double lb) {
+  if (_data) {
+    _data->lb = lb;
+  }
+  update_test_data();
+}
+
+void optimizer::set_upper_bound(double ub) {
+  if (_data) {
+    _data->ub = ub;
+  }
+  update_test_data();
+}
+
+void optimizer::set_curve_type(int ct) {
+  if (ct == _curve.index()) return;
+  switch (ct) {
+    case CurveType::Parabola:
+      _curve = parabola_new{};
+      break;
+    case CurveType::Ellipse:
+      _curve = ellipse_new{};
+      break;
+  }
+  update_table_cb();
+}
+
+void optimizer::set_initial(size_t i, double value) {
+  _fit_result.reset();
+  std::visit(
+      [&](auto& curve) -> void {
+        if (i >= curve.initial.size()) {
+          throw std::runtime_error("index out of bounds");
+        }
+        curve.initial(i) = value;
+      },
+      _curve);
+  plot_slope_cb();
+}
+
+void optimizer::set_delta(size_t i, double value) {
+  _fit_result.reset();
+  std::visit(
+      [&](auto& curve) -> void {
+        if (i >= curve.delta.size()) {
+          throw std::runtime_error("index out of bounds");
+        }
+        curve.delta(i) = value;
+      },
+      _curve);
+  plot_slope_cb();
+}
+
+void optimizer::flip_curve() {
+  std::visit([](auto& curve) { curve.flip = !curve.flip; }, _curve);
+  plot_slope_cb();
+}
+
+void optimizer::calc_fit() {
+  _fit_result =
+      std::visit([&](auto& curve) { return curve.fit(*_test_data); }, _curve);
+  plot_slope_cb();
+  for (auto param : _fit_result->first) {
+    std::println("param: {}", param);
+  }
+  plot_error_cb();
+}
+
+void optimizer::set_plot_input_cb(
+    std::function<void(Eigen::VectorXd xs, std::vector<Eigen::VectorXd> yss)>
+        cb) {
+  _plot_input_cb = cb;
+}
+
+void optimizer::set_init_bounds_cb(
+    std::function<void(double lb, double ub, double step)> cb) {
+  _init_bounds_cb = cb;
+}
+
+void optimizer::set_plot_slope_cb(
+    std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys,
+                       Eigen::VectorXd ys_curve)>
+        cb) {
+  _plot_slope_cb = cb;
+}
+
+void optimizer::set_update_table_cb(
+    std::function<
+        void(std::vector<std::tuple<std::string, double, double>> entries)>
+        cb) {
+  _param_table_cb = cb;
+}
+
+void optimizer::set_plot_error_cb(
+    std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys)> cb) {
+  _plot_error_cb = cb;
+}
+
+void optimizer::update_test_data() {
+  _fit_result.reset();
+  if (!_data) {
+    _test_data.reset();
+    return;
+  }
+  auto lb_it = std::lower_bound(_data->_raw._xs.begin(), _data->_raw._xs.end(),
+                                _data->lb);
+  auto ub_it = std::upper_bound(_data->_raw._xs.begin(), _data->_raw._xs.end(),
+                                _data->ub);
+  ssize_t n_samples = std::distance(lb_it, ub_it);
+  if (n_samples < 1) {
+    _test_data.reset();
+    plot_input_cb();
+    return;
+  }
+
+  size_t start_index = std::distance(_data->_raw._xs.begin(), lb_it);
+
+  Eigen::VectorXd xs(n_samples);
+  Eigen::VectorXd ys = Eigen::VectorXd::Zero(n_samples);
+
+  std::copy(lb_it, ub_it, xs.begin());
+  for (size_t i = 0; i < n_samples; i++) {
+    for (auto r : _data->_raw._yss) {
+      ys(i) += r(start_index + i);
+    }
+    ys(i) /= _data->_raw._yss.size();
+  }
+  double center = (*lb_it + *std::prev(ub_it)) / 2;
+  _test_data.emplace(xs, ys, center);
+  plot_slope_cb();
+}
+
+void optimizer::init_bounds_cb() {
+  if (_init_bounds_cb) {
+    _init_bounds_cb(_data->_raw._xs(0), _data->_raw._xs.tail(1)(0),
+                    _data->_raw._xs(1) - _data->_raw._xs(0));
+  }
+}
+
+void optimizer::plot_input_cb() {
+  if (_plot_input_cb) {
+    _plot_input_cb(_data->_raw._xs, _data->_raw._yss);
+  }
+}
+
+void optimizer::plot_slope_cb() {
+  if (_plot_slope_cb) {
+    auto ys_curve =
+        std::visit([&](auto curve) { return curve.plot(*_test_data); }, _curve);
+    _plot_slope_cb(_test_data->_xs, _test_data->_ys, ys_curve);
+  }
+}
+
+void optimizer::update_table_cb() {
+  std::vector<std::tuple<std::string, double, double>> res;
+  std::visit(
+      [&]<typename T>(T& curve) {
+        for (size_t i = 0; i < curve.param_names.size(); i++) {
+          res.emplace_back(curve.param_names[i], curve.initial[i],
+                           curve.delta[i]);
+        }
+      },
+      _curve);
+  _param_table_cb(res);
+}
+
+void optimizer::plot_error_cb() {
+  if (_plot_error_cb) {
+    Eigen::VectorXd ys_curve = std::visit(
+        [&](auto curve) { return curve.plot(*_test_data, _fit_result->first); },
+        _curve);
+    Eigen::VectorXd error = _test_data->_ys - ys_curve;
+    _plot_error_cb(_test_data->_xs, error);
+
+    
+
+    fftw_complex* out = fftw_alloc_complex(error.size());
+    double* input = &error(0);
+    fftw_plan plan = fftw_plan_dft_r2c_1d(error.size(), input, out, FFTW_FORWARD);
+    fftw_execute(plan);
+    fftw_destroy_plan(plan);
+    Eigen::VectorXd out_abs(error.size());
+    for (size_t i = 0; i < error.size(); i++) {
+      out_abs(i) = out[i][0] * out[i][0] + out[i][1]*out[i][1];
+    }
+    fftw_free(out);
+    std::println("begin fftw");
+    for (auto v : out_abs) {
+      std::println("  {}", v);
+    }
+    std::println("end fftw");
+  }
+}
+
+}  // namespace autoopt