2

I need to fit data to a curve that is chosen by the user as the sum of pre-defined functions. The idea is to create something similar to what peak-o-mat does to spectral data.

Example: Fit some data to a function that is the sum of one linear function (two parameters) and one cauchy function (two parameters).

If I know in compilation time, that my curve to fit is the sum of a linear and a cauchy, I could create a specific residual for that. But the problem is that I only know it on runtime.

What I would like to do is:

1) Define residuals based on a list of functions I expose to the user:

struct Residual {
  Residual(double x, double y) : m_x(x), m_y(y) {}

 protected:
  // Observations for a sample.
  const double m_x;
  const double m_y;
};

struct LorentzianResidual : Residual {
  LorentzianResidual(double x, double y) : Residual(x, y) {}

  template <typename T>
  bool operator()(const T* const m, const T* const c, T* residual) const {
    residual[0] = T(m_y) - (1 / M_PI) * (0.5 * c[0]) /
                               ((T(m_x) - m[0]) * (T(m_x) - m[0]) +
                                (0.5 * c[0]) * (0.5 * c[0]));
    return true;
  }
};

struct LinearResidual : Residual {
  LinearResidual(double x, double y) : Residual(x, y) {}

  template <typename T>
  bool operator()(const T* const m, const T* const c, T* residual) const {
    residual[0] = T(m_y) - (m[0] * T(m_x) + c[0]);
    return true;
  }
};

2) Solve a ceres::Problem adding residual blocks, based on the combination of functions chosen by the user.

I was thinking two alternatives:

a) Create a curve fitting class, with a member that has all the chosen functions and an array of parameters for each of the functions. Then I would create a Residual inside this class, that has x and y as parameters, but would loop through this functions and parameters, returning the sum of their outputs. The problem is that I would not be able to add the residual block using: 

ceres::CostFunction* cost_function1 =
      new ceres::AutoDiffCostFunction<Residual, 1, 1, 1>(
          new Residual(xdata[i], ydata[i]));

Because of the template arguments <Residual, 1, 1, 1> that I would know only on runtime (and there is a maximum of 9).

So I was looking for a different alternative, in which I could add each residuals separately. I tried one simple example with a linear and cauchy (lorentzian) residuals, but it is not working.

    std::vector<double> linear_coeffs{0.0, 0.0};
    std::vector<double> lorentz_coeffs{0.0, 0.0};

    ceres::Problem problem;
    for (size_t i = 0; i < xdata.size(); ++i) {

      ceres::CostFunction* cost_function1 =
          new ceres::AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
              new ExponentialResidual(xdata[i], ydata[i]));
      problem.AddResidualBlock(cost_function1, nullptr, &linear_coeffs[0],
                               &linear_coeffs[1]);
      ceres::CostFunction* cost_function2 =
          new ceres::AutoDiffCostFunction<LinearResidual, 1, 1, 1>(
              new LinearResidual(xdata[i], ydata[i]));
      problem.AddResidualBlock(cost_function2, nullptr, &lorentz_coeffs[0],
                               &lorentz_coeffs[1]);
    }

Error in evaluating the ResidualBlock.

There are two possible reasons. Either the CostFunction did not evaluate and fill all    
residual and jacobians that were requested or there was a non-finite value (nan/infinite)
generated during the or jacobian computation. 

Residual Block size: 2 parameter blocks x 1 residuals

I checked documentation, but I didn't find anything on curve fitting with combination of different functions.

Does anyone have an idea of how to make it work?

Taiguara Tupinambás
  • 175
  • 12
  • 1
    If you can fit 1.0 or 0.0 times the function, for example "fitted_parameter * 1.0 * log(x)" or "fitted_parameter * 0.0 * sin(x)", then you could use all of the different functions as an array of pointers with an index that points both to a function and to an array of doubles that contains iser-selectable ones or zeroes. – James Phillips Oct 11 '19 at 21:19
  • It would be an interesting workaround. But I guess it doesn't apply to my problem, because the user might select multiple functions of the same type. I managed to make it work with the `ceres::DynamicAutoDiffCostFunction`, but in order to generate the residual, I am using switch cases and it's not looking good. If I find an elegant solution, I will post here. – Taiguara Tupinambás Oct 14 '19 at 17:12

0 Answers0