forum.alglib.net

void function_cx_1_func(const alglib::real_1d_array &c, const alglib::real_1d_array &x, double &func, void *ptr)
{
    //  The arguments follow the format of the alglib's lsfit subpackage
    //  This function requires the func to return a double.
    //  c is a 1d array of parameters
    //  x is temperature, which is the only variable in this case
    //  ptr is NULL in by the package's default.
    //  cout << "iterred" << endl;
    map <double,double> simulation_result;
    if (
        //(gmm != c[0]) ||
        (Tm != c[0]) ||
        (g31 != c[1]) ||
        (g32 != c[2]) ||
        (gg31 != c[3]) ||
        (gg32 != c[4]) ||
        (bn != c[5]) ||
        (bu != c[6]) ||
        (mn != c[7]) ||
        (mu != c[8])
        )
    {
        //gmm = c[0];
        Tm = c[0];
        g31 = c[1];
        g32 = c[2];
        gg31 = c[3];
        gg32 = c[4];
        bn = c[5];
        bu = c[6];
        mn = c[7];
        mu = c[8];
        Tjumpave = Tjump_ave();
    }
    double T = 5;
//    for (mat::iterator entry = Tjumpave.begin(); entry != Tjumpave.end(); ++entry) {
//        simulation_result[T] = *entry;
//        T += 5;
//    }
    for (int i = 0; i != 9; ++i)
    {
        simulation_result[T] = Tjumpave(0,i);
        T += 5;
    }
    func = simulation_result[x[0]];
}
vector<double> params_optim()
{
    //  Optimization using lsfit
    //  The params_optim function optimizes the values of the simulation parameters
    //  The parameters are: Tm, g31, g32, gg31, gg32, bn, bu, mn, mu (9)
    vector<double> par = {1,1,1,1,1,1,1,1,1};
    try
    {
        alglib::real_2d_array T = "[[5],[10],[15],[20],[25],[30],[35],[40],[45],[50],[55],[60],[65],[70],[75],[80],[85],[90],[95]]";  //  Temperatures in experiments (independent var)
        alglib::real_1d_array exp = "[346.301,346.312,346.432,346.394,346.471,346.605,346.797,346.948,347.121,347.384,347.626,348.08,348.561,349.333,350.404,351.761,352.975,354.17,354.809]";  //  Exp data
        alglib::real_1d_array params = "[80,550,400,0,0,347,354,0.013,0.012]"; //  in the order given above
        alglib::real_1d_array bndl = "[75,-INF,-INF,-INF,-INF,345,350,0.005,0.005]";  //  lower bounds of params
        alglib::real_1d_array bndu = "[83,+INF,+INF,+INF,+INF,349,360,0.030,0.030]";  //  upper bounds of params
        //  Define the stopping criteria
        double epsf = 0.0; //  stops when the parameter value of the next iter is less than epsf away from the previous's
        double epsx = 0.000001; //  stops when the scaled step vector (scaled by LSfitsetscale) is <= epsx
        alglib::real_1d_array s = "[1e+1,1e+2,1e+2,1e+2,1e+2,1e+2,1e+2,1e-2,1e-2]"; //  set the scales (orders of magnitude) of the parameter
        alglib::ae_int_t info; //  to store how the "more_information code number"
        alglib::ae_int_t maxits = 100; //  set the number of iteration to perform
        alglib::lsfitstate state; //  store the "state"
        alglib::lsfitreport rep;  //  store the report
        double diffstep = 5; //  step to take for each param: too high, miss points. too low, higher accuracy but more round-off error
       
        alglib::lsfitsetscale(state, s);
        alglib::ae_int_t n = 19;  //  number of points
        alglib::ae_int_t m = 1;  //  number of dimension (only temperature)
        alglib::ae_int_t k = 9;  //  number of parameters to fit (10)
        alglib::lsfitcreatef(T, exp, params, n,m,k, diffstep, state); //  create a fitting instance
        alglib::lsfitsetbc(state, bndl, bndu);
        alglib::lsfitsetcond(state, epsf, epsx, maxits); //  before fitting, set a stopping point
        alglib::lsfitfit(state, function_cx_1_func); //  fit
        alglib::lsfitresults(state, info, params, rep);
        par = { params[0],params[1],params[2],params[3],params[4],params[5],params[6],params[7],params[8] };
    }
    catch (const alglib::ap_error &e)
    {
        printf("error msg: %s\n", e.msg.c_str());
    }
    return par;
}