TOF - iterative fitting

[chiaracarminati]

Implement schemes for iterative fitting

[chiaracarminati]

From Joachim's suggestion:

keeping some parameters fixed is not explicitly supported by lmfit,
but can be easily implemented on top of it: you just need to move
the free parameters from your full parameter vector to a dedicated
fit parameter vector, and move fit results back.

example:

 // Transfer parameters from C->P to Par
 //   (lmmin modifies Par; fit_evaluate transfers values back from Par to C->P)
 int np = fc->nP;
 vector<double> Par;
 for (int ip = 0; ip < np; ++ip)
     if (C->ParAttr[ip] != 'x')
         Par.push_back(C->P[ip]);
 int npfree = Par.size();
 vector<double> ParCovar(npfree * npfree);

 // Levenberg-Marquardt routine from library lmfit:

 data.timeout = time(nullptr) + std::thread::hardware_concurrency() * maxtime;

 lm_status_struct status;
 lmmin2(
     npfree, &(Par[0]), NULL, &(ParCovar[0]), nd, NULL, &data, fit_evaluate, &control,
     &status);
 if (status.outcome == 11)
     throw S("exception in fit function");

 // Transfer parameters back
 for (int i = 0; i < np * np; ++i)
     C->PCovar[i] = 0;
 int ipfit = 0;
 for (int ip = 0; ip < np; ++ip) {
     if (C->ParAttr[ip] != 'x') {
         C->P[ip] = Par[ipfit];
         int jpfit = 0;
         for (int jp = 0; jp < np; ++jp)
             if (C->ParAttr[jp] != 'x')
                 C->PCovar[ip * np + jp] = ParCovar[ipfit * npfree + (jpfit++)];
         ipfit++;
     }
 }

[chiaracarminati]

I'd recommend to use one fixed function fit_evaluate, which reacts flexibly
to the parameters provided through the *data pointer. The example below is
again from Frida, where fd points to a data set and fc to a curve (= function
definition plus parameters). This of course you will organize differently.
Look for the "PARAMETER TRANSLATION" paragraph for the key idea: you just
revert the encoding P -> par, depending on which parameters are free.

Hope this helps - Joachim

//! Data transfer between frida and lm_minimize callback functions.

typedef struct {
const COld* fd;
COlc* fc;
int k;
int j;
double timeout;
} FitDatTyp;

//! Callback routine for use in lm_minimize: function evaluation.

static void
fit_evaluate(const double* par, int m_dat, const void* data, double* fvec, int* userbreak)
// arguments (required by lmmin):
// par r/w function parameters, to be optimized
// m_dat r number of data points
// data r pointer to arbitrary data, here defined as FitDatTyp
// fvec w vector of residues, to be minimized
// userbreak w can be use to terminate the minimization
{
FitDatTyp* mydata {(FitDatTyp*)data};
COlc* fc = mydata->fc;
CCurve* c = fc->VC(mydata->j);
try {
if (mydata->timeout && time(nullptr) > mydata->timeout)
throw S("reached fit time out");
const COld* fd = mydata->fd;
COlc::TWgt wt = fc->weighing;

     /****  PARAMETER TRANSLATION  ****/
     int ip = 0;
     for (int i = 0; i < fc->nP; ++i)
         if (fc->VC(mydata->j)->ParAttr[i] != 'x')
             c->P[i] = par[ip++];

     const CSpec* s = fd->VS(mydata->j);
     bool want_error = wt == COlc::_VAR || wt == COlc::_VARC;

     RObjVec cu = fc->eval_curve(s->x, mydata->k, mydata->j, want_error);

     if (compute_residues(fvec, 0, mydata->j, s, cu, wt))
         *userbreak = -1;

 } catch (string& s) {
     cout << "fit_evaluate j=" << mydata->j;
     for (int i = 0; i < fc->nP; ++i)
         cout << " p" << i << "=" << c->P[i];
     cout << ": " << s << "\n";
     *userbreak = -1;
 } catch (...) {
     cout << "BUG: fit_evaluate catched invalid exception\n";
     *userbreak = -1;
 }

}