如何使用 gsl_multifit_fdf 翻译 scipy.curve_fit

Question

我有一个带有特定包络的波形。我需要从那个信封中获取一些参数，所以我将曲线拟合到极值。极值的测量是熊猫数据框

peakDf

我先是用下面的python代码做的

import numpy as np
import pandas as pd
from scipy.optimize import curve_fit

def envelope(timestamps, start, startFlt, endFlt, end, startAmpl, endAmpl):
    startSlope = startAmpl/(startFlt-start)
    fltSlope = (endAmpl - startAmpl) / (endFlt - startFlt)
    endSlope = (-endAmpl/(end-endFlt))
    result = np.zeros(len(timestamps))
    for index, timestamp in enumerate(timestamps):
        if timestamp > start:
            if timestamp < startFlt:
                result[index] = (startSlope * (timestamp - start))
            elif timestamp < endFlt:
                result[index] = (fltSlope * (timestamp - startFlt) + startAmpl)
            elif timestamp < end:
                result[index] = (endSlope * (timestamp - endFlt) + endAmpl)
    return result

params, covs = curve_fit(envelope, peakDf.timestamps, np.abs(peakDf.value))

这给了我以下情节。您可以看到曲线的预期形状。开始、startFlt、endFlt 和结束时间用垂直线标记。

现在我需要使用 gsl 1.15 将其转换为 c++03。我尽力了，但由于某种原因，它在

fitEnvelope(...)

行

gslStatus = gsl_multifit_fdfsolver_set(solver, &envelopeFit, p);

中失败了 我进一步查看了代码，调用

gsl_multifit_fdfsolver_lmsder->set(...)

时似乎失败了。

我将它作为共享库提供给另一个程序，我不能与 gdb 一起使用，所以我不知道我得到了什么样的错误。

我已经在这方面工作了一段时间，但我不知道自己做错了什么。非常感谢关于我的翻译的任何提示：

struct ValueTimePair {
    std::vector<double> values;
    std::vector<double> timestamps;
};

std::vector<double> envelope(std::vector<double> timestamps, double start, double startFlt, double endFlt, double end, double startAmpl, double endAmpl) {
    std::vector<double> result(timestamps.size(), 0);
    double startSlope = startAmpl/(startFlt-start);
    double fltSlope = (endAmpl - startAmpl) / (endFlt - startFlt);
    double endSlope = (-endAmpl/(end-endFlt));
    for (size_t i=0; i<result.size(); i++) {
        double timestamp = timestamps[i];
        if (timestamp > start) {
            if (timestamp < startFlt) {
                result[i] = (startSlope * (timestamp - start));
            } else if (timestamp < endFlt) {
                result[i] = (fltSlope * (timestamp - startFlt) + startAmpl);
            } else if (timestamp < end) {
                result[i] = (endSlope * (timestamp - endFlt) + endAmpl);
            }
        }
    }
    return result;
}


int envelope_f(const gsl_vector* x, void* data, gsl_vector* f) {
    double start = gsl_vector_get(x, 0);
    double startFlt = gsl_vector_get(x, 1);
    double endFlt = gsl_vector_get(x, 2);
    double end = gsl_vector_get(x, 3);
    double startAmpl = gsl_vector_get(x, 4);
    double endAmpl = gsl_vector_get(x, 5);
    ValueTimePair* measurements = static_cast<ValueTimePair*>(data);

    std::vector<double> result = envelope(measurements->timestamps, start, startFlt, endFlt, end, startAmpl, endAmpl);
    
    for (size_t i = 0; i < result.size(); ++i) {
        gsl_vector_set(f, i, measurements->values[i]-result[i]);
    }
    return GSL_SUCCESS;
}

int envelope_df(const gsl_vector* x, void* data, gsl_matrix* J) {
    // TODO: work in progress. Check if Jakobian matrix is generated correctly
    double start = gsl_vector_get(x, 0);
    double startFlt = gsl_vector_get(x, 1);
    double endFlt = gsl_vector_get(x, 2);
    double end = gsl_vector_get(x, 3);
    double startAmpl = gsl_vector_get(x, 4);
    double endAmpl = gsl_vector_get(x, 5);
    ValueTimePair* measurements = static_cast<ValueTimePair*>(data);

    std::vector<double> envelopeValues = envelope(measurements->timestamps, start, startFlt, endFlt, end, startAmpl, endAmpl);
    
    for (size_t i = 0; i < envelopeValues.size(); ++i) {
        // calculate partial derivatives for each parameter at this point
        double timestamp = measurements->timestamps[i];
        double startSlope = timestamp >= start ? (timestamp < startFlt ? (timestamp - start) : 0.0) : 0.0;
        double fltSlope = timestamp >= startFlt && timestamp < endFlt ? (timestamp - startFlt) : 0.0;
        double endSlope = timestamp >= endFlt && timestamp < end ? (timestamp - end) : 0.0;
        
        double startAmplDeriv = startSlope;
        double endAmplDeriv = endSlope;
        double startDeriv = startSlope * startAmpl / (startFlt - start);
        double startFltDeriv = fltSlope * (startAmpl - endAmpl) / (endFlt - startFlt);
        double endFltDeriv = endSlope * (-endAmpl) / (end - endFlt);
        double endDeriv = endSlope * (-endAmpl) / (end - endFlt);
        // assign the partial derivatives to the Jacobian matrix
        gsl_matrix_set(J, i, 0, startDeriv);
        gsl_matrix_set(J, i, 1, startFltDeriv);
        gsl_matrix_set(J, i, 2, endFltDeriv);
        gsl_matrix_set(J, i, 3, endDeriv);
        gsl_matrix_set(J, i, 4, startAmplDeriv);
        gsl_matrix_set(J, i, 5, endAmplDeriv);
    }
    return GSL_SUCCESS;
}

void fitEnvelope(ValueTimePair& measurements, const EnvelopeParameters* guess, EnvelopeParameters* result, EnvelopeParameters* errors, size_t maxIter) {
    size_t n = measurements.timestamps.size();
    const int parameterN = 6;
    int gslStatus = 0;
    gsl_vector* p = gsl_vector_alloc(parameterN);
    gsl_vector_set(p, 0, guess->start);
    gsl_vector_set(p, 1, guess->startFlt);
    gsl_vector_set(p, 2, guess->endFlt);
    gsl_vector_set(p, 3, guess->end);
    gsl_vector_set(p, 4, guess->amplitudeStartFlt);
    gsl_vector_set(p, 5, guess->amplitudeEndFlt);
    if (!p) {
        throw std::runtime_error("Error allocating vector p");
    }
    gsl_multifit_function_fdf envelopeFit;
    envelopeFit.f = &envelope_f;
    envelopeFit.df = &envelope_df;
    envelopeFit.n = n;
    envelopeFit.p = parameterN;
    envelopeFit.params = &measurements;
    if (!envelopeFit.f || !envelopeFit.df || !envelopeFit.params) {
        LOG_ERROR_FORMAT_IF(logger_WaveformAnalyser,"Error allocating vector envelopeFit");
        throw std::runtime_error("WaveformAnalyser::fitEnvelope: Error initializing envelopeFit");
    }
    
    // Use lmsder because we calculated the jacobian matrix by hand. Use lmder otherwise
    const gsl_multifit_fdfsolver_type * T = gsl_multifit_fdfsolver_lmsder;
    gsl_multifit_fdfsolver* solver = gsl_multifit_fdfsolver_alloc(T, n, parameterN);
    if(!solver) {
        throw std::runtime_error("fitEnvelope: Could not initialize curve fit solver");
    }

    //FIXME:######## NB: HERE IS THE ERROR The following line exits the program. No error logs are returned
    gslStatus = gsl_multifit_fdfsolver_set(solver, &envelopeFit, p);
    
    if(gslStatus) {
        std::ostringstream oss;
        oss << "fitEnvelope::gsl_multifit_fdfsolver_set: GSL shows error: " << gsl_strerror(gslStatus);
        throw std::exception(oss.str());
    }

    size_t iter = 0;
    do {
        iter++;
        gslStatus = gsl_multifit_fdfsolver_iterate(solver);
        if (gslStatus) break;
        gslStatus = gsl_multifit_test_delta(solver->dx, solver->x, 1e-4, 1e-4);
    } while (gslStatus == GSL_CONTINUE && iter < maxIter);
    if (gslStatus == GSL_CONTINUE) {
        // Log hit max iterations
    } else if(gslStatus) {
        std::ostringstream oss;
        oss << "fitEnvelope::fit iteration: GSL shows error: " << gsl_strerror(gslStatus);
        throw std::exception(oss.str());
    }
    
    gsl_matrix* pcov = gsl_matrix_alloc(parameterN, parameterN);
    gsl_multifit_covar(solver->J, 0.0, pcov);

    gsl_multifit_fdfsolver_free(solver);

    errors->start = sqrt(gsl_matrix_get(pcov, 0, 0));
    errors->startFlt = sqrt(gsl_matrix_get(pcov, 1, 1));
    errors->endFlt = sqrt(gsl_matrix_get(pcov, 2, 2));
    errors->end = sqrt(gsl_matrix_get(pcov, 3, 3));
    errors->amplitudeStartFlt = sqrt(gsl_matrix_get(pcov, 4, 4));
    errors->amplitudeEndFlt = sqrt(gsl_matrix_get(pcov, 5, 5));
    gsl_matrix_free(pcov);

    result->start = gsl_vector_get(p, 0);
    result->startFlt = gsl_vector_get(p, 1);
    result->endFlt = gsl_vector_get(p, 2);
    result->end = gsl_vector_get(p, 3);
    result->amplitudeStartFlt = gsl_vector_get(p, 4);
    result->amplitudeEndFlt = gsl_vector_get(p, 5);
    gsl_vector_free(p);
}