add parabola

.clang-format

@@ -1,3 +0,0 @@
----
-BasedOnStyle: Google
-IndentWidth: 2

.gitignore

@@ -4,9 +4,6 @@ result
 # cmake
 build
 
-# clangd
-.cache
-
 # direnv files
 .direnv
 .envrc

CMakeLists.txt

@@ -1,10 +1,8 @@
 cmake_minimum_required(VERSION 4.0)
 
-set(CMAKE_C_COMPILER clang)
+project(autoopt)
-set(CMAKE_CXX_COMPILER clang++)
-set(CMAKE_CXX_STANDARD 23)
 
-project(slopefit)
+set(CMAKE_CXX_STANDARD 23)
 
 add_compile_options(-Wall -Werror -Wpedantic)
 
@@ -17,32 +15,22 @@ execute_process(
   OUTPUT_STRIP_TRAILING_WHITESPACE OUTPUT_VARIABLE nanobind_ROOT)
 find_package(nanobind CONFIG REQUIRED)
 
-file(GLOB_RECURSE SOURCES lib/*.cpp)
+file(GLOB_RECURSE SOURCES src/*.cpp)
 
 add_library(autoopt STATIC ${SOURCES})
 target_include_directories(autoopt PUBLIC include)
 target_link_libraries(autoopt PUBLIC Eigen3::Eigen)
-target_link_libraries(autoopt PRIVATE fftw3)
 
-nanobind_add_module(slopefit NOMINSIZE src/module.cpp)
+install(DIRECTORY include/autoopt DESTINATION include)
-nanobind_add_stub(
-    slopefit_stub
-    MODULE slopefit
-    OUTPUT slopefit.pyi
-    PYTHON_PATH $<TARGET_FILE_DIR:slopefit>
-    DEPENDS slopefit
-)
 
+nanobind_add_module(slopefit NOMINSIZE pysrc/module.cpp)
 target_link_libraries(slopefit PRIVATE autoopt)
-install(TARGETS slopefit DESTINATION .)
+install(TARGETS slopefit DESTINATION lib/python3.13/site-packages)
-install(FILES $<TARGET_FILE_DIR:slopefit>/slopefit.pyi DESTINATION .)
 
 # add tests if gtest is found
-enable_testing()
-add_subdirectory(tests)
-
 find_library(GTestPackage gtest QUIET)
 if(GTestPackage)
-
+  enable_testing()
+  add_subdirectory(tests)
 endif()
 

README.md

@@ -1 +0,0 @@
-# slopefit

build.nix

@@ -1,44 +0,0 @@
-{
-  cmake,
-  ninja,
-  eigen,
-  fftw,
-  clang-tools,
-  python3Packages,
-  doCheck ? false,
-  gtest,
-}:
-python3Packages.buildPythonPackage {
-  pname = "slopefit";
-  version = "0.1.0";
-  pyproject = true;
-
-  dependencies = [
-    python3Packages.numpy
-    python3Packages.matplotlib
-    python3Packages.pyside6
-    python3Packages.pandas
-    python3Packages.pyqtgraph
-  ];
-
-  buildInputs = [
-    eigen
-    fftw
-  ]
-  ++ (if doCheck then [ gtest ] else [ ]);
-
-  nativeBuildInputs = [
-    cmake
-    clang-tools
-    ninja
-    python3Packages.scikit-build-core
-    python3Packages.nanobind
-  ];
-
-  buildPhase = ''
-    cd ..
-    pypaBuildPhase
-  '';
-
-  src = ./.;
-}

examples/test.py

@@ -0,0 +1,60 @@
+import sys
+import os.path as path
+
+build_dir = path.join(path.dirname(__file__), '..', 'build')
+print('Adding build directory to sys.path:', build_dir)
+sys.path.append(build_dir)
+
+import slopefit as sf
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+fl = 1200.0 # mm
+theta_mrad = 3.0 # mrad
+theta_rad = theta_mrad * 1e-3
+theta = theta_rad / np.pi * 180.0 # deg
+
+parab = sf.ParabolaParams(fl, theta)
+
+filename = 'microMAX_vsru_mispu_absum.dat'
+dirname = path.dirname(__file__)
+
+data = sf.import_dat_file(path.join(dirname, filename))
+
+# flip both x and y
+data[:, 0] = -data[:, 0]
+data[:, 1] = -data[:, 1]
+
+n_skip = 20
+data = data[n_skip:-n_skip, :]
+
+ref = parab(data[:, 0])
+
+# convert ref from slope to rad
+# ref = np.arctan(ref)
+# convert to arcsec
+# ref = ref / np.pi * 180.0 * 3600.0
+
+
+plt.plot(data[:, 0], data[:, 1], 'x')
+plt.plot(data[:, 0], ref, '-')
+plt.xlabel('x (mm)')
+plt.ylabel('slope (arcsec)')
+plt.show()
+
+d_parab = sf.ParabolaParams(100.0, 0.00001)
+
+fitted = sf.fit_parabola(data, parab, d_parab)
+
+print('Fitted parameters:')
+print('fl =', fitted.focal_length)
+print('theta =', fitted.theta * np.pi * 1e3 / 180.0, 'mrad')
+
+fitted_ref = fitted(data[:, 0])
+
+plt.plot(data[:, 0], data[:, 1], 'x')
+plt.plot(data[:, 0], fitted_ref, '-')
+plt.xlabel('x (mm)')
+plt.ylabel('slope (arcsec)')
+plt.show()

flake.lock

@@ -20,11 +20,11 @@
     },
     "nixpkgs": {
       "locked": {
-        "lastModified": 1789546076,
+        "lastModified": 1768564909,
-        "narHash": "sha256-zVxLZiSnmaaPLwnhj7pwmqe3axBg/C6nG5JZsJMh2g4=",
+        "narHash": "sha256-Kell/SpJYVkHWMvnhqJz/8DqQg2b6PguxVWOuadbHCc=",
         "owner": "nixos",
         "repo": "nixpkgs",
-        "rev": "b1b875982b17dabde9b4a37f3e229e74913e6db3",
+        "rev": "e4bae1bd10c9c57b2cf517953ab70060a828ee6f",
         "type": "github"
       },
       "original": {

flake.nix

@@ -1,5 +1,5 @@
 {
-  description = "slopefit";
+  description = "auto-opt-cpp";
 
   inputs = {
     nixpkgs.url = "github:nixos/nixpkgs?ref=nixos-unstable";
@@ -16,24 +16,66 @@
       system:
       let
         pkgs = nixpkgs.legacyPackages.${system};
-        pkg = import ./build.nix;
+        pkg =
+          {
+            stdenv,
+            cmake,
+            ninja,
+            lib,
+            eigen,
+            gtest,
+            python313,
+            python313Packages,
+            withTests ? false,
+          }:
+          let
+            fs = lib.fileset;
+          in
+          stdenv.mkDerivation {
+            pname = "autoopt-cpp";
+            version = "0.1.0";
+
+            src = fs.toSource {
+              root = ./.;
+              fileset = fs.unions [
+                ./include
+                ./CMakeLists.txt
+                ./tests
+                ./src
+                ./pysrc
+              ];
+            };
+
+            buildInputs = lib.optionals withTests [ gtest ];
+
+            nativeBuildInputs = [
+              cmake
+              ninja
+              eigen
+              (python313.withPackages (
+                ps: with ps; [
+                  nanobind
+                  matplotlib
+                  numpy
+                ]
+              ))
+            ];
+          };
       in
       rec {
         packages = rec {
-          slopefit = pkgs.callPackage pkg { };
+          autoopt-cpp = pkgs.callPackage pkg { };
-          slopefit-with-tests = pkgs.callPackage pkg { doCheck = true; };
+          autoopt-cpp-with-tests = pkgs.callPackage pkg { withTests = true; };
-          default = slopefit;
+          default = autoopt-cpp;
         };
         devShells = {
           default = pkgs.mkShell {
             inputsFrom = [
-              packages.slopefit-with-tests
+              packages.autoopt-cpp-with-tests
             ];
 
-            buildInputs = [ pkgs.python3Packages.pip ];
-
             shellHook = ''
-              echo "Development shell for slopefit"
+              echo "Development shell for autoopt-cpp"
             '';
           };
         };

include/autoopt/btls.hpp

@@ -1,5 +1,7 @@
 #pragma once
 
+#include <iostream>
+
 #include "autoopt/optimization_problem.hpp"
 
 namespace autoopt {

include/autoopt/ellipse.hpp

@@ -1,25 +1,16 @@
 #pragma once
 
-#include "autoopt/parameters.hpp"
 #include "autoopt/quadric.hpp"
 
 namespace autoopt {
 
-struct ellipse_new : curve_base {
+template <typename T>
-  // T focal_length;
+struct ellipse {
-  ellipse_new() : curve_base() {
+  T left_arm;
-    initial = Eigen::VectorXd::Zero(4);
+  T right_arm;
-    delta = Eigen::VectorXd::Zero(4);
+  T entrance_angle;
-    param_names =
-        std::vector<const char*>{"left arm", "right arm", "entrance angle", "rotation"};
-  }
 
-  template <typename T>
+  constexpr quadric<T> to_quadric() const {
-  quadric<T> grapher_from_params(const Eigen::VectorX<T>& params) const {
-    T left_arm = params(0);
-    T right_arm = params(1);
-    T entrance_angle = params(2);
-    T alpha = params(3);
     T a = (left_arm + right_arm) / T{2};
     T a2 = a * a;
     
@@ -52,7 +43,7 @@ struct ellipse_new : curve_base {
     T E = -B * c_x - T{2} * C * c_y;
     T F = A * c_x * c_x + B * c_x * c_y + C * c_y * c_y - a2 * b2;
     
-    return quadric<T>(A, B, C, D, E, F).rotated_by(alpha);
+    return quadric<T>(A, B, C, D, E, F);
   }
 };
 

include/autoopt/importer.hpp

@@ -1,23 +0,0 @@
-#pragma once
-
-#include <vector>
-#include <Eigen/Eigen>
-#include <generator>
-#include <filesystem>
-#include <fstream>
-
-namespace autoopt {
-struct raw_data {
-  Eigen::VectorXd _xs;
-  std::vector<Eigen::VectorXd> _yss;
-};
-
-std::generator<std::pair<double, double>> parse_dat_file(
-    const std::filesystem::path& filename);
-
-std::pair<std::vector<double>, std::vector<double>> parse_row(auto& iter,
-                                                              auto sentinel);
-
-raw_data read_data(const std::filesystem::path& filename);
-
-}

include/autoopt/interface.hpp

@@ -0,0 +1,16 @@
+#pragma once
+
+#include <eigen3/Eigen/Eigen>
+
+namespace autoopt {
+
+Eigen::VectorX<double> fit_ellipse(
+    const std::vector<std::pair<double, double>>& data,
+    const Eigen::VectorX<double>& inital_params,
+    const Eigen::VectorX<double>& delta);
+
+Eigen::VectorX<double> fit_parabola(
+    const std::vector<std::pair<double, double>>& data,
+    const Eigen::VectorX<double>& inital_params,
+    const Eigen::VectorX<double>& delta);
+}  // namespace autoopt

include/autoopt/optimization_problem.hpp

@@ -1,5 +1,7 @@
 #pragma once
 
+#include <cstddef>
+
 #include "autoopt/derivative.hpp"
 
 namespace autoopt {
@@ -74,7 +76,7 @@ struct log_barrier_optimization_problem
     auto base_grad = _base_problem.gradient(params);
     Eigen::VectorX<T> barrier_grad = autoopt::gradient<T>(
         [this]<typename U>(const Eigen::VectorX<U>& p) {
-          return this->barrier_term<U>(p);
+          return barrier_term<U>(p);
         },
         params);
     Eigen::VectorX<T> total_grad(params.size());
@@ -86,7 +88,7 @@ struct log_barrier_optimization_problem
     auto base_hess = _base_problem.hessian(params);
     Eigen::MatrixX<T> barrier_hess = autoopt::hessian<T>(
         [this]<typename U>(const Eigen::VectorX<U>& p) {
-          return this->barrier_term<U>(p);
+          return barrier_term<U>(p);
         },
         params);
     Eigen::MatrixX<T> total_hess(params.size(), params.size());

include/autoopt/optimizer.hpp

@@ -1,64 +0,0 @@
-#pragma once
-
-#include <filesystem>
-
-#include <functional>
-#include <variant>
-#include <vector>
-#include "autoopt/parabola_new.hpp"
-#include "autoopt/ellipse_new.hpp"
-#include "autoopt/importer.hpp"
-
-namespace autoopt {
-
-enum CurveType {
-  Parabola,
-  Ellipse
-};
-
-using curves = std::variant<parabola_new, ellipse_new>;
-
-struct optimizer {
-
-public:
-  void import_file(std::filesystem::path file_path);
-  void set_lower_bound(double lb);
-  void set_upper_bound(double ub);
-
-  void set_curve_type(int curve_type);
-  void set_initial(size_t i, double value);
-  void set_delta(size_t i, double value);
-  void flip_curve();
-  void calc_fit();
-
-  void set_plot_input_cb(std::function<void(Eigen::VectorXd xs, std::vector<Eigen::VectorXd> yss)> cb);
-  void set_plot_slope_cb(std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys, Eigen::VectorXd ys_curve)> cb);
-  void set_init_bounds_cb(std::function<void(double lb, double ub, double step)> cb);
-  void set_update_table_cb(std::function<void(std::vector<std::tuple<std::string, double, double>> entries)> cb);
-  void set_plot_error_cb(std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys)> cb);
-
-private:
-  struct data {
-    raw_data _raw;
-    double lb, ub;
-  };
-  std::optional<data> _data;
-  std::optional<test_data> _test_data;
-  curves _curve;
-  std::optional<std::pair<Eigen::VectorXd, double>> _fit_result;
-
-  std::function<void(Eigen::VectorXd xs, std::vector<Eigen::VectorXd> yss)> _plot_input_cb;
-  std::function<void(double lb, double ub, double step)> _init_bounds_cb;
-  std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys, Eigen::VectorXd ys_curve)> _plot_slope_cb;
-  std::function<void(std::vector<std::tuple<std::string, double, double>> entries)> _param_table_cb;
-  std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys)> _plot_error_cb;
-
-  void update_test_data();
-  void init_bounds_cb();
-  void plot_input_cb();
-  void plot_slope_cb();
-  void update_table_cb();
-  void plot_error_cb();
-};
-
-}  // namespace autoopt

include/autoopt/parabola.hpp

@@ -0,0 +1,27 @@
+#pragma once
+
+#include "autoopt/quadric.hpp"
+
+namespace autoopt
+{
+  
+template <typename T>
+struct parabola {
+  // T focal_length;
+  T exit_arm;
+  T entrance_angle;
+
+  quadric<T> to_quadric() const {
+    // T a = T{1} / (T{4} * focal_length);
+    // T x = T{2} * focal_length / std::tan(entrance_angle);
+    // T y = a * x * x;
+    T x = exit_arm * std::sin(T{2.0} * entrance_angle);
+    T f = T{0.5} * (T{1.0} - std::cos(T{2.0} * entrance_angle)) * exit_arm;
+    T a = T{1.0} / (T{4.0} * f);
+    T y = a * x * x;
+    quadric<T> q = quadric(a, T{0}, T{0}, T{0}, T{-1.0}, T{0});
+    return q.translated_by(-x, -y).rotated_by(entrance_angle - T{M_PI_2});
+  }
+};
+
+} // namespace autoopt

include/autoopt/parabola_new.hpp

@@ -1,33 +0,0 @@
-#pragma once
-
-#include "autoopt/parameters.hpp"
-#include "autoopt/quadric.hpp"
-
-namespace autoopt {
-
-struct parabola_new : curve_base {
-  // T focal_length;
-  parabola_new() : curve_base() {
-    initial = Eigen::VectorXd::Zero(3);
-    delta = Eigen::VectorXd::Zero(3);
-    param_names =
-        std::vector<const char*>{"focus", "entrance angle", "rotation"};
-  }
-
-  template <typename T>
-  quadric<T> grapher_from_params(const Eigen::VectorX<T>& params) const {
-    T exit_arm = params(0);
-    T entrance_angle = params(1);
-    T alpha = params(2);
-
-    T x = exit_arm * std::sin(T{2.0} * entrance_angle);
-    T f = T{0.5} * (T{1.0} - std::cos(T{2.0} * entrance_angle)) * exit_arm;
-    T a = T{1.0} / (T{4.0} * f);
-    T y = a * x * x;
-    quadric<T> q = quadric(a, T{0}, T{0}, T{0}, T{-1.0}, T{0});
-    return q.translated_by(-x, -y).rotated_by(entrance_angle - T{M_PI_2} +
-                                              alpha);
-  }
-};
-
-}  // namespace autoopt

include/autoopt/parameters.hpp

@@ -1,76 +0,0 @@
-#pragma once
-
-#include <Eigen/Eigen>
-#include <autoopt/btls.hpp>
-#include <autoopt/optimization_problem.hpp>
-#include <cmath>
-#include <cstddef>
-#include <optional>
-#include <utility>
-#include <vector>
-
-#include "autoopt/util.hpp"
-
-namespace autoopt {
-
-struct test_data {
-  Eigen::VectorXd _xs;
-  Eigen::VectorXd _ys;
-  double _center;
-
-  size_t size() const { return _xs.size(); }
-
-  std::pair<double, double> operator[](size_t i) const {
-    return {_xs(i) - _center, _ys(i)};
-  };
-};
-
-struct curve_base {
-  Eigen::VectorXd initial;
-  Eigen::VectorXd delta;
-  std::vector<const char*> param_names;
-  bool flip = false;
-
-  std::pair<Eigen::VectorXd, double> fit(this auto&& self,
-                                         const test_data& measurements) {
-    auto opt_func = [&]<typename T>(const Eigen::VectorX<T>& params) {
-      auto grapher = self.template grapher_from_params<T>(params);
-      T error{0.0};
-      for (size_t i = 0; i < measurements.size(); i++) {
-        auto [x, y] = measurements[i];
-        T distance = T{self.flip ? -1.0 : 1.0} *
-                         std::atan(grapher.slope_at(T{self.flip ? -x : x})) -
-                     T{arcsec2rad(y)};
-        error = error + distance * distance;
-      }
-      return error / T(measurements.size());
-    };
-
-    auto_diff_optimization_problem problem(opt_func, self.initial);
-
-    log_barrier_optimization_problem lb_problem(problem, self.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();
-    return {fitted_params, std::sqrt(problem.objective(fitted_params))};
-  }
-
-  Eigen::VectorXd plot(this auto&& self, const test_data& measurements,
-                       std::optional<Eigen::VectorXd> params = {}) {
-    auto grapher = self.template grapher_from_params<double>(
-        params.value_or(self.initial));
-    Eigen::VectorXd ys(measurements.size());
-    for (size_t i = 0; i < measurements.size(); i++) {
-      double x = measurements[i].first;
-      ys(i) = rad2arcsec((self.flip ? -1.0 : 1.0) *
-                         atan(grapher.slope_at(self.flip ? -x : x)));
-    }
-    return ys;
-  }
-};
-
-}  // namespace autoopt

include/autoopt/reactive.hpp

@@ -1,79 +0,0 @@
-
-#include <array>
-#include <cstdint>
-#include <functional>
-#include <tuple>
-
-namespace autoopt {
-
-template <typename T>
-struct source {
-  T _value;
-  std::uint64_t _version = 0;
-
-  source() = default;
-  explicit source(T value) : _value(std::move(value)) {}
-
-  virtual T& get() { return _value; };
-
-  bool set(T value) {
-    if (value == _value) return false;
-    _value = std::move(value);
-    ++_version;
-    return true;
-  }
-};
-
-template <typename O, typename... Is>
-struct transform : source<O> {
-  std::function<O(Is&...)> _func;
-  std::tuple<source<Is>&...> _sources;
-  std::array<std::uint64_t, sizeof...(Is)> _versions{};
-
-  transform(std::function<O(Is&...)> func, source<Is>&... sources)
-      : _func(func), _sources(sources...) {
-    this->_value = std::apply([&](source<Is>&... sources) {
-      return _func(sources.get()...);
-    }, _sources);
-    this->_versions = current_versions();
-  }
-
-  virtual O& get() override {
-    auto cur_versions = current_versions();
-    if (_versions == cur_versions) return this->_value;
-
-    O next = std::apply([&](source<Is>&... sources) {
-      return _func(sources.get()...);
-    }, _sources);
-
-    if (next != this->_value) {
-      this->_value = next;
-      ++this->_version;
-    }
-    return this->_value;
-  }
-
- private: 
-  std::array<std::uint64_t, sizeof...(Is)> current_versions() {
-    return std::apply([&](source<Is>&... sources){
-      
-      return std::array<std::uint64_t, sizeof...(Is)>{sources._version...};
-    }, _sources);
-  }
-};
-
-template <typename T>
-struct sink {
-  source<T>& _upstream;
-  std::uint64_t _version;
-  std::function<void(T&)> _cb;
-
-  void maybe_update() {
-    if (_version != _upstream._version) {
-      _cb(_upstream.get());
-      _version = _upstream._version;
-    }
-  }
-};
-
-}  // namespace autoopt

lib/importer.cpp

@@ -1,73 +0,0 @@
-#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/optimizer.cpp

@@ -1,224 +0,0 @@
-#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

pyproject.toml

@@ -1,36 +0,0 @@
-[build-system]
-requires = ["scikit-build-core>=0.11", "nanobind"]
-build-backend = "scikit_build_core.build"
-
-
-[tool.scikit-build]
-# Protect the configuration against future changes in scikit-build-core
-minimum-version = "0.4"
-
-# Setuptools-style build caching in a local directory
-build-dir = "build/{wheel_tag}"
-
-[project]
-name = "slopefit"
-version = "0.1.0"
-authors = [
-  { name="Valentin Leon Stöcker", email="valentin@stoeckermail.de" },
-]
-description = "fit conic sections"
-readme = "README.md"
-requires-python = ">=3.9"
-dependencies = [
-  "numpy",
-  "matplotlib"
-  "pyside6"
-]
-classifiers = [
-    "Programming Language :: Python :: 3",
-    "Operating System :: OS Independent",
-]
-license = "MIT"
-
-[project.urls]
-Homepage = "https://gitea.valentinstoecker.de/vls/slopfit"
-Issues = "https://gitea.valentinstoecker.de/vls/slopfit/issues"
-

pysrc/module.cpp

@@ -0,0 +1,266 @@
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/string.h>
+#include <nanobind/eigen/dense.h>
+
+#include <autoopt/btls.hpp>
+#include <autoopt/ellipse.hpp>
+#include <autoopt/parabola.hpp>
+#include <autoopt/interface.hpp>
+#include <autoopt/optimization_problem.hpp>
+#include <autoopt/util.hpp>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+namespace nb = nanobind;
+
+struct input_data {
+  size_t n_points;
+  double* data;
+
+  input_data(size_t n) : n_points(n) {
+    if (n_points == 0) {
+      data = nullptr;
+    } else {
+      data = new double[n_points * 2];
+    }
+  }
+
+  input_data(const input_data&) = delete;
+  input_data& operator=(const input_data&) = delete;
+
+  input_data(input_data&& other) noexcept
+      : n_points(other.n_points), data(other.data) {
+    other.data = nullptr;
+    other.n_points = 0;
+  }
+
+  input_data& operator=(input_data&& other) noexcept {
+    if (this != &other) {
+      delete[] data;
+      n_points = other.n_points;
+      data = other.data;
+      other.data = nullptr;
+      other.n_points = 0;
+    }
+    return *this;
+  }
+
+  ~input_data() { delete[] data; }
+
+  void set_point(size_t i, double x, double y) {
+    data[i * 2] = x;
+    data[i * 2 + 1] = y;
+  }
+};
+
+input_data read_data(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::unordered_map<double, size_t> index_map;
+  std::vector<double> x_values;
+  std::vector<std::vector<double>> data_points;
+
+  std::string line;
+
+  size_t index = 0;
+
+  double x_avg = 0.0;
+
+  while (std::getline(file, line)) {
+    if (line.empty() || line[0] == '#') {
+      std::cout << "Skipped: " << line << std::endl;
+      continue;
+    }
+    std::istringstream iss(line);
+    double x, y;
+    iss >> x;
+    for (size_t i = 1; i < 5; ++i) {
+      iss >> y;  // skip unused columns
+    }
+
+    y = autoopt::arcsec2rad(-y);
+
+    if (index_map.find(x) == index_map.end()) {
+      index_map[x] = index++;
+      x_values.push_back(x);
+      data_points.emplace_back();
+      data_points.back().push_back(y);
+      x_avg += x;
+      continue;
+    }
+    data_points[index_map[x]].push_back(y);
+  }
+  x_avg /= x_values.size();
+
+  input_data result(index);
+  for (size_t i = 0; i < x_values.size(); ++i) {
+    double x = x_values[i];
+    double y_avg = 0.0;
+    for (double y : data_points[i]) {
+      y_avg += y;
+    }
+    y_avg /= data_points[i].size();
+    result.set_point(i, x - x_avg, y_avg);
+  }
+  return result;
+}
+
+struct ellipse_params {
+  double left_arm;
+  double right_arm;
+  double theta;
+};
+
+struct parabola_params {
+  double focal_length;
+  double theta;
+};
+
+NB_MODULE(slopefit, m) {
+  m.def("import_dat_file", [](std::string filename) {
+    input_data data = read_data(filename);
+    std::cout << "First data point: (" << data.data[0] << ", " << data.data[1]
+              << ")" << std::endl;
+    nb::ndarray<double, nb::numpy> array(data.data, {data.n_points, 2});
+    return nb::cast(array);
+  });
+
+  nb::class_<ellipse_params>(m, "EllipseParams")
+      .def(nb::init<double, double, double>())
+      .def_rw("left_arm", &ellipse_params::left_arm)
+      .def_rw("right_arm", &ellipse_params::right_arm)
+      .def_rw("theta", &ellipse_params::theta)
+      .def("__repr__", [](const ellipse_params& params) {
+        std::ostringstream oss;
+        oss << "EllipseParams(left_arm=" << params.left_arm
+            << ", right_arm=" << params.right_arm << ", theta=" << params.theta
+            << ")";
+        return oss.str();
+      }).def("__call__", [](const ellipse_params& params, nb::ndarray<double, nb::ndim<1>> xs) {
+        auto e = autoopt::ellipse(params.left_arm, params.right_arm,
+                                  autoopt::deg2rad(params.theta));
+        autoopt::quadric<double> q = e.to_quadric();
+        Eigen::VectorXd ys(xs.shape(0));
+        for (size_t i = 0; i < xs.shape(0); ++i) {
+          ys(i) = q.slope_at(xs(i));
+        }
+        return ys;
+      })
+      ;
+    nb::class_<parabola_params>(m, "ParabolaParams")
+      .def(nb::init<double, double>())
+      .def_rw("focal_length", &parabola_params::focal_length)
+      .def_rw("theta", &parabola_params::theta)
+      .def("__repr__", [](const parabola_params& params) {
+        std::ostringstream oss;
+        oss << "ParabolaParams(focal_length=" << params.focal_length
+            << ", theta=" << params.theta << ")";
+        return oss.str();
+      })
+      .def("__call__", [](const parabola_params& params, nb::ndarray<double, nb::ndim<1>> xs) {
+        auto p = autoopt::parabola(params.focal_length, autoopt::deg2rad(params.theta));
+        autoopt::quadric<double> q = p.to_quadric();
+        Eigen::VectorXd ys(xs.shape(0));
+        for (size_t i = 0; i < xs.shape(0); ++i) {
+          ys(i) = q.slope_at(xs(i));
+        }
+        return ys;
+      });
+
+  m.def(
+      "fit_ellipse",
+      [](nb::ndarray<double, nb::ndim<2>> data, ellipse_params initial_params,
+         ellipse_params delta) -> ellipse_params {
+        std::cout << "Fitting ellipse to data with " << data.shape(0)
+                  << " points." << std::endl;
+
+        if (data.shape(1) != 2) {
+          throw std::runtime_error("Data array must have shape (n_points, 2)");
+        }
+
+        std::vector<std::pair<double, double>> data_vec;
+        for (size_t i = 0; i < data.shape(0); ++i) {
+          data_vec.emplace_back(data(i, 0), data(i, 1));
+        }
+
+        std::cout << "Initial parameters: "
+                  << "left_arm=" << initial_params.left_arm
+                  << ", right_arm=" << initial_params.right_arm
+                  << ", theta=" << initial_params.theta << std::endl;
+
+        double midpoint_y = data_vec[data_vec.size() / 2].second;
+
+        Eigen::VectorX<double> init_params(4);
+        init_params(0) = initial_params.left_arm;
+        init_params(1) = initial_params.right_arm;
+        init_params(2) = autoopt::deg2rad(initial_params.theta);
+        init_params(3) = midpoint_y;
+
+        Eigen::VectorX<double> deltas(4);
+        deltas(0) = delta.left_arm;
+        deltas(1) = delta.right_arm;
+        deltas(2) = autoopt::deg2rad(delta.theta);
+        deltas(3) = autoopt::deg2rad(0.1);
+
+        std::cout << "calculating fit..." << std::endl;
+        Eigen::VectorX<double> fitted_params =
+            autoopt::fit_ellipse(data_vec, init_params, deltas);
+        ellipse_params result;
+        result.left_arm = fitted_params(0);
+        result.right_arm = fitted_params(1);
+        result.theta = autoopt::rad2deg(fitted_params(2));
+        return result;
+      });
+
+  m.def(
+      "fit_parabola",
+      [](nb::ndarray<double, nb::ndim<2>> data, parabola_params initial_params,
+         parabola_params delta) -> parabola_params {
+        std::cout << "Fitting parabola to data with " << data.shape(0)
+                  << " points." << std::endl;
+
+        if (data.shape(1) != 2) {
+          throw std::runtime_error("Data array must have shape (n_points, 2)");
+        }
+
+        std::vector<std::pair<double, double>> data_vec;
+        for (size_t i = 0; i < data.shape(0); ++i) {
+          data_vec.emplace_back(data(i, 0), data(i, 1));
+        }
+
+        std::cout << "Initial parameters: "
+                  << "focal_length=" << initial_params.focal_length
+                  << ", theta=" << initial_params.theta  << std::endl;
+
+        double midpoint_y = data_vec[data_vec.size() / 2].second;
+
+        Eigen::VectorX<double> init_params(3);
+        init_params(0) = initial_params.focal_length;
+        init_params(1) = autoopt::deg2rad(initial_params.theta);
+        init_params(2) = midpoint_y;
+
+        Eigen::VectorX<double> deltas(3);
+        deltas(0) = delta.focal_length;
+        deltas(1) = autoopt::deg2rad(delta.theta);
+        deltas(2) = autoopt::deg2rad(0.1);
+
+        std::cout << "calculating fit..." << std::endl;
+        Eigen::VectorX<double> fitted_params =
+            autoopt::fit_parabola(data_vec, init_params, deltas);
+        parabola_params result;
+        result.focal_length = fitted_params(0);
+        result.theta = autoopt::rad2deg(fitted_params(1));
+        return result;
+      });
+}

src/interface.cpp

@@ -0,0 +1,89 @@
+#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;
+}

src/module.cpp

@@ -1,53 +0,0 @@
-#include <nanobind/eigen/dense.h>
-#include <nanobind/nanobind.h>
-#include <nanobind/ndarray.h>
-#include <nanobind/stl/filesystem.h>
-#include <nanobind/stl/optional.h>
-#include <nanobind/stl/pair.h>
-#include <nanobind/stl/tuple.h>
-#include <nanobind/stl/string.h>
-#include <nanobind/stl/vector.h>
-#include <nanobind/stl/function.h>
-
-#include <Eigen/Eigen>
-#include <autoopt/optimizer.hpp>
-#include <autoopt/util.hpp>
-#include <format>
-#include <optional>
-
-namespace nb = nanobind;
-
-namespace std {
-
-template <typename T>
-struct formatter<optional<T>> {
-  constexpr auto parse(format_parse_context& ctx) { return ctx.begin(); }
-
-  auto format(const optional<T>& opt, format_context& ctx) const {
-    if (opt) {
-      return format_to(ctx.out(), "{}", opt.value());
-    } else {
-      return format_to(ctx.out(), "None");
-    }
-  }
-};
-
-}  // namespace std
-
-NB_MODULE(slopefit, m) {
-  nb::class_<autoopt::optimizer>(m, "Optimizer")
-      .def(nb::init<>())
-      .def("import_file", &autoopt::optimizer::import_file)
-      .def("set_lower_bound", &autoopt::optimizer::set_lower_bound)
-      .def("set_upper_bound", &autoopt::optimizer::set_upper_bound)
-      .def("set_curve_type", &autoopt::optimizer::set_curve_type)
-      .def("set_initial", &autoopt::optimizer::set_initial)
-      .def("set_delta", &autoopt::optimizer::set_delta)
-      .def("flip_curve", &autoopt::optimizer::flip_curve)
-      .def("calc_fit", &autoopt::optimizer::calc_fit)
-      .def("set_plot_input_cb", &autoopt::optimizer::set_plot_input_cb)
-      .def("set_init_bounds_cb", &autoopt::optimizer::set_init_bounds_cb)
-      .def("set_plot_slope_cb", &autoopt::optimizer::set_plot_slope_cb)
-      .def("set_update_table_cb", &autoopt::optimizer::set_update_table_cb)
-      .def("set_plot_error_cb", &autoopt::optimizer::set_plot_error_cb);
-}

tests/ellipse.cpp

@@ -0,0 +1,79 @@
+#include <gtest/gtest.h>
+
+#include <autoopt/ellipse.hpp>
+#include <autoopt/optimization_problem.hpp>
+#include <autoopt/util.hpp>
+#include <autoopt/btls.hpp>
+#include <iomanip>
+#include <iostream>
+
+using namespace autoopt;
+
+TEST(Ellipse, Slope) {
+  ellipse<double> e{100, 1000, deg2rad(1.0)};  // entrance angle 1 degree
+  quadric<double> q = e.to_quadric();
+
+  EXPECT_NEAR(q.slope_at(-10), -0.0010305116165301856, 1e-9);
+  EXPECT_NEAR(q.slope_at(0), 0.0, 1e-9);
+  EXPECT_NEAR(q.slope_at(10), 0.00090001261192696272, 1e-9);
+}
+
+TEST(Ellipse, ParamGradient) {
+  std::vector<std::pair<double, double>> data_points = {
+      {-10.0, -0.00103}, {0.0, 0.0}, {10, 0.0009}};
+
+  Eigen::VectorX<double> params(4);
+  params << 100, 1000, deg2rad(1.0), 0.0;
+
+  auto loss_func = [&data_points]<typename T>(const Eigen::VectorX<T>& p) {
+    ellipse<T> e{T{p(0)}, T{p(1)}, T{p(2)}};
+    quadric<T> q = e.to_quadric().rotated_by(T{p(3)});
+    T loss = T{0};
+    for (const auto& [x, y_true] : data_points) {
+      T y_pred = q.slope_at(T{x});
+      T error = y_pred - T{y_true};
+      loss = loss + error * error;
+    }
+    return loss / T(data_points.size());
+  };
+
+  auto_diff_optimization_problem<double, decltype(loss_func)> problem(loss_func, params);
+
+  auto grad = problem.gradient(params);
+
+  std::cout << grad << std::endl;
+
+  EXPECT_NEAR(grad(0), -3.54845114759293e-12, 1e-15);  // d/d(left_arm)
+  EXPECT_NEAR(grad(1), -3.0016523630530093e-14, 1e-15);  // d/d(right_arm)
+  EXPECT_NEAR(grad(2), 2.0569619167404501e-08, 1e-15);   // d/d(entrance_angle)
+  EXPECT_NEAR(grad(3), -3.3267028547673413e-07, 1e-15);  // d/d(rotation_angle)
+
+  auto hess = problem.hessian(params);
+
+  std::cout << hess << std::endl;
+
+  Eigen::VectorX<double> params_delta(4);
+  params_delta << 10.0, 10.0, deg2rad(0.1), deg2rad(0.1);
+
+
+  // log barrier
+  log_barrier_optimization_problem<double> log_barrier_problem(
+      problem,
+      params_delta,
+      1e-5);
+  
+  while (log_barrier_problem._barrier_strength > 1e-20) {
+    btls(log_barrier_problem);
+    log_barrier_problem._barrier_strength *= 1e-2;
+  }
+
+  std::cout << "Optimum params:" << std::endl;
+  std::cout << "left_arm:       " << log_barrier_problem.x()(0) << std::endl;
+  std::cout << "right_arm:      " << log_barrier_problem.x()(1) << std::endl;
+  std::cout << "entrance_angle: " << rad2deg(log_barrier_problem.x()(2)) << std::endl;
+  std::cout << "rotation_angle: " << rad2deg(log_barrier_problem.x()(3)) << std::endl;
+  std::cout << "Optimum objective:" << std::endl;
+  std::cout << problem.objective(problem.x()) << std::endl;
+  std::cout << "Optimum grad:" << std::endl;
+  std::cout << problem.gradient(problem.x()) << std::endl;
+}

tests/quadric.cpp

@@ -1,6 +1,7 @@
 #include <gtest/gtest.h>
 
 #include "autoopt/quadric.hpp"
+#include <iostream>
 
 using namespace autoopt;
 

tests/reacitve.cpp

@@ -1,21 +0,0 @@
-#include <gtest/gtest.h>
-
-#include <autoopt/reactive.hpp>
-
-TEST(Reactive, Simple) {
-  using namespace autoopt;
-
-  source<int> s1(0);
-  source<int> s2(0);
-  source<int> s3(0);
-  transform<int, int, int> t1([](int x, int y) { return x + y; }, s1,s2);
-  transform<int, int, int> t2([](int x, int y) { return x + y; }, s2,s3);
-
-  ASSERT_EQ(t2.get(), 0);
-  ASSERT_TRUE(s1.set(1));
-  ASSERT_EQ(t2.get(), 1);
-  s2.set(2);
-  ASSERT_EQ(t2.get(), 3);
-  s3.set(3);
-  ASSERT_EQ(t2.get(), 6);
-}

tests/test_fit.cpp

@@ -0,0 +1,82 @@
+#include <gtest/gtest.h>
+
+#include <autoopt/ellipse.hpp>
+#include <autoopt/optimization_problem.hpp>
+#include <autoopt/util.hpp>
+#include <autoopt/btls.hpp>
+#include <fstream>
+#include <unordered_map>
+#include <autoopt/interface.hpp>
+#include <iomanip>
+
+
+using namespace autoopt;
+
+std::vector<std::pair<double, double>> read_data(const std::string& filename) {
+  std::fstream file(filename, std::ios::in);
+
+  std::unordered_map<double, size_t> index_map;
+  std::vector<double> x_values;
+  std::vector<std::vector<double>> data_points;
+
+  std::string line;
+
+  size_t index = 0;
+
+  double x_avg = 0.0;
+
+  for(;std::getline(file, line);) {
+    if (line.empty() || 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
+    }
+
+    y = arcsec2rad(-y);
+
+    if (index_map.find(x) == index_map.end()) {
+      index_map[x] = index++;
+      x_values.push_back(x);
+      data_points.emplace_back();
+      data_points.back().push_back(y);
+      x_avg += x;
+      continue;
+    }
+    data_points[index_map[x]].push_back(y);
+  }
+  x_avg /= x_values.size();
+  std::vector<std::pair<double, double>> result;
+  for (size_t i = 40; i < x_values.size() - 40; ++i) {
+    double x = x_values[i];
+    double y_avg = 0.0;
+    for (double y : data_points[i]) {
+      y_avg += y;
+    }
+    y_avg /= data_points[i].size();
+    result.emplace_back(x - x_avg, y_avg);
+  }
+  return result;
+}
+
+TEST(Fit, EllipseFit) {
+  auto data_points = read_data("tests/input/ellipse.dat");
+
+  std::cout << "Read " << data_points.size() << " data points." << std::endl;
+
+  auto mid_point = data_points[data_points.size() / 2];
+
+  std::cout << "Mid point: (" << mid_point.first << ", " << mid_point.second << ")" << std::endl;
+
+  Eigen::VectorX<double> initial_params(4);
+  initial_params << 6900.0, 500.0, deg2rad(2.0), mid_point.second;
+  
+  Eigen::VectorX<double> delta(4);
+  delta << 10.0, 10.0, deg2rad(0.1), deg2rad(0.1);
+
+  auto res = fit_ellipse(data_points, initial_params, delta);
+  std::cout << res.transpose() << std::endl;
+}

tests/test_rework.cpp

@@ -1,49 +0,0 @@
-#include <gtest/gtest.h>
-
-#include <autoopt/parabola_new.hpp>
-#include <print>
-#include "autoopt/util.hpp"
-
-
-using namespace autoopt;
-
-
-TEST(New, Parabola) {
-  parabola_new parab;
-  parab.flip = false;
-  parab.initial(0) = 50.0;
-  parab.initial(1) = autoopt::deg2rad(1.0);
-  parab.initial(2) = autoopt::deg2rad(0.05);
-
-  parab.delta(0) = 10.0;
-  parab.delta(1) = autoopt::deg2rad(0.1);
-  parab.delta(2) = autoopt::deg2rad(0.1);
-
-  test_data td;
-  td._xs = Eigen::VectorXd(21);
-  td._ys = Eigen::VectorXd(21);
-  td._center = 10.0;
-
-  for (size_t i = 0; i < 21; i++) {
-    td._xs(i) = static_cast<double>(i);
-  }
-
-  auto res = parab.plot(td);
-
-  for (auto y : res) {
-    std::println("y: {}", y);
-  }
-
-  td._ys = res;
-
-  parab.initial(0) = 54.0;
-  parab.initial(1) = autoopt::deg2rad(1.05);
-
-  auto [fit, rms] = parab.fit(td);
-
-  for (auto y : fit) {
-    std::println("fit: {}", y);
-  }
-
-  std::println("RMS: {}", rms);
-}

todo.md

@@ -1,9 +0,0 @@
-# main class
-
-methods:
-- import_file
-- adjust_bounds
-- select_curve
-- set_param
-- flip
-- calc_fit

ui/__init__.py

@@ -1 +0,0 @@
-

ui/app.py

@@ -1,221 +0,0 @@
-import PySide6.QtWidgets as qw
-import PySide6.QtCore as qc
-from PySide6.QtCore import Slot, Qt
-import PySide6.QtGui as qg
-import pyqtgraph as pg
-
-import sys
-
-sys.path.append('../build')
-
-import slopefit as sf
-
-class ParamTableModel(qc.QAbstractTableModel):
-  HEADERS = ['property', 'initial value', 'delta']
-  init_changed = qc.Signal(int, float)
-  delta_changed = qc.Signal(int, float)
-
-  def __init__(self, parent = None):
-    super().__init__(parent)
-    self._rows = []
-  
-  def rowCount(self, parent=qc.QModelIndex()):
-    return len(self._rows)
-  
-  def columnCount(self, parent=qc.QModelIndex()):
-    return len(self.HEADERS)
-  
-  def data(self, index, role=Qt.ItemDataRole.DisplayRole):
-    if not index.isValid():
-      return None
-    row, col = index.row(), index.column()
-    value = self._rows[row][col]
-
-    if role in (Qt.ItemDataRole.DisplayRole, Qt.ItemDataRole.EditRole):
-      return "" if value is None else str(value)
-    if role == Qt.ItemDataRole.TextAlignmentRole and col > 0:
-      return int(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter)
-    return None
-  
-  def setData(self, index, value, role=Qt.ItemDataRole.EditRole):
-    if not index.isValid() or role != Qt.ItemDataRole.EditRole:
-      return False
-    row, col = index.row(), index.column()
-    if col == 0:                       # name column: plain text
-      return False
-    else:                              # x / dx columns: float
-      try:
-        self._rows[row][col] = float(value) if value != "" else None
-      except ValueError:
-        return False
-    self.dataChanged.emit(index, index, [Qt.ItemDataRole.DisplayRole, Qt.ItemDataRole.EditRole])
-    if col == 1:
-      self.init_changed.emit(row, float(value))
-    elif col == 2:
-      self.delta_changed.emit(row, float(value))
-    return True
-  
-  def flags(self, index):
-    if not index.isValid() or index.column() == 0:
-      return Qt.ItemFlag.NoItemFlags
-    return Qt.ItemFlag.ItemIsEnabled | Qt.ItemFlag.ItemIsSelectable | Qt.ItemFlag.ItemIsEditable
-
-  def headerData(self, section, orientation, role=Qt.ItemDataRole.DisplayRole):
-    if role == Qt.ItemDataRole.DisplayRole and orientation == Qt.Orientation.Horizontal:
-      return self.HEADERS[section]
-    return None
-
-  def set_params(self, params):
-    self.beginResetModel()
-    self._rows = [[name, v, d] for (name, v, d) in params]
-    self.endResetModel()
-
-
-
-class MainWindow(qw.QMainWindow):
-  def __init__(self):
-    super().__init__()
-    self.setWindowTitle('slopefit')
-
-    
-    self.input_plot = pg.PlotWidget()
-
-    self.lower_bounds = qw.QDoubleSpinBox()
-    self.lower_bounds.setValue(0.0)
-    self.lower_bounds.valueChanged.connect(self.set_lb)
-
-    self.upper_bounds = qw.QDoubleSpinBox()
-    self.upper_bounds.setValue(0.0)
-    self.upper_bounds.valueChanged.connect(self.set_ub)
-
-
-    bounds_layout = qw.QHBoxLayout()
-    bounds_layout.addWidget(self.lower_bounds)
-    bounds_layout.addWidget(self.upper_bounds)
-
-    bounds_widget = qw.QWidget()
-    bounds_widget.setLayout(bounds_layout)
-
-    self.slope_plot = pg.PlotWidget()
-
-    self.curve_box = qw.QComboBox()
-    self.curve_box.addItems(['parabola', 'ellipse'])
-    self.curve_box.currentIndexChanged.connect(self.set_curve_i)
-
-    self.table = qw.QTableView()
-    self.param_model = ParamTableModel()
-    self.param_model.init_changed.connect(self.set_init)
-    self.param_model.delta_changed.connect(self.set_delta)
-    self.table.setModel(self.param_model)
-    self.table.horizontalHeader().setSectionResizeMode(qw.QHeaderView.Stretch)
-
-    self.flip_curve = qw.QCheckBox()
-    self.flip_curve.checkStateChanged.connect(self.set_check)
-
-    self.fit_button = qw.QPushButton()
-    self.fit_button.setText("calculate fit")
-    self.fit_button.clicked.connect(self.do_fit)
-
-    self.error_plot = pg.PlotWidget()
-    plot_group_1 = qw.QHBoxLayout()
-    plot_group_1.addWidget(self.input_plot)
-    plot_group_1.addWidget(self.slope_plot)
-
-    layout = qw.QVBoxLayout()
-    layout.addWidget(qw.QWidget(layout=plot_group_1))
-    layout.addWidget(bounds_widget)
-    layout.addWidget(self.flip_curve)
-    layout.addWidget(self.curve_box)
-    layout.addWidget(self.table)
-    layout.addWidget(self.fit_button)
-    layout.addWidget(self.error_plot)
-
-    central_widget = qw.QWidget()
-    central_widget.setLayout(layout)
-    self.setCentralWidget(central_widget)
-
-    self.init_menu()
-    
-    self.optimizer = sf.Optimizer()
-    self.optimizer.set_plot_input_cb(self.plot_input_cb)
-    self.optimizer.set_init_bounds_cb(self.init_bounds_cb)
-    self.optimizer.set_plot_slope_cb(self.plot_slope_cb)
-    self.optimizer.set_update_table_cb(self.update_table_cb)
-    self.optimizer.set_plot_error_cb(self.plot_error_cb)
-
-  @Slot(float)
-  def set_lb(self, lb):
-    self.optimizer.set_lower_bound(lb)
-    
-  @Slot(float)
-  def set_ub(self, ub):
-    self.optimizer.set_upper_bound(ub)
-  
-  @Slot(int)
-  def set_curve_i(self, i):
-    self.optimizer.set_curve_type(i)
-
-  @Slot(int, float)
-  def set_init(self, i, v):
-    self.optimizer.set_initial(i, v)
-  
-  @Slot(int, float)
-  def set_delta(self, i, v):
-    self.optimizer.set_delta(i, v)
-
-  @Slot(bool)
-  def set_check(self, flip):
-    self.optimizer.flip_curve()
-
-  @Slot()
-  def do_fit(self):
-    self.optimizer.calc_fit()
-
-  def init_menu(self):
-    import_action = qg.QAction('&import', self)
-    import_action.triggered.connect(self.import_file)
-
-    menu = self.menuBar()
-    file_menu = menu.addMenu('&file')
-    file_menu.addAction(import_action)
-
-  def import_file(self):
-    print("import file...")
-    fn, _ = qw.QFileDialog.getOpenFileName(self, caption='import file...', filter='data files (*.dat)')
-    print(f"open file {fn}")
-    self.optimizer.import_file(fn)
-
-  def plot_input_cb(self, xs, yss):
-    self.input_plot.clear()
-    for ys in yss:
-      self.input_plot.plot(xs, ys)
-
-  def init_bounds_cb(self, lb, ub, step):
-    self.lower_bounds.setRange(lb, ub)
-    self.lower_bounds.setSingleStep(step)
-    self.lower_bounds.setValue(lb)
-
-    self.upper_bounds.setRange(lb, ub)
-    self.upper_bounds.setSingleStep(step)
-    self.upper_bounds.setValue(ub)
-  
-  def plot_slope_cb(self, xs, ys_data, ys_curve):
-    self.slope_plot.clear()
-    self.slope_plot.plot(xs, ys_data)
-    self.slope_plot.plot(xs, ys_curve)
-  
-  def update_table_cb(self, entries):
-    self.param_model.set_params(entries)
-    print(entries)
-
-  def plot_error_cb(self, xs, ys):
-    self.error_plot.clear()
-    self.error_plot.plot(xs, ys)
-
-
-app = qw.QApplication(sys.argv)
-
-window = MainWindow()
-window.show()
-
-app.exec()

ui/model.py

@@ -1,4 +0,0 @@
-import PySide6.QtWidgets
-
-
-