MKL fitting example

//
//  C-style wrapper for MinPack callback function
//
class FitFunctionWrapper
{
public:
    static void FitFunction_MKL(MKL_INT *m, MKL_INT *n, double *x, double *f, void *priv)
    {
        if (priv)
        {
            CFitFunction* pFFObj = reinterpret_cast<CFitFunction*>(priv);
            return pFFObj->MKL_Evaluate_function(m, n, x, f);
        }
    }
};


bool CFitFunction::MKL_Fit()
{
    /* n - number of function variables
       m - dimension of function value */
    MKL_INT         n = m_parameters.size(), m = m_data_size;
    /* precisions for stop-criteria (see manual for more details) */
    std::vector<double>eps(6,1E-8);

    /* iter1 - maximum number of iterations
       iter2 - maximum number of iterations of calculation of trial-step */
    MKL_INT         iter1 = m_max_iterations, iter2 = 100;
    /* initial step bound */
    double  rs = 1.0;
    /* reverse communication interface parameter */
    MKL_INT         RCI_Request;
    /* controls of rci cycle */
    MKL_INT         successful;
    /* number of iterations */
    MKL_INT         iter;
    /* number of stop-criterion */
    MKL_INT         st_cr;
    /* initial and final residauls */
    double r1=1E12, r2=1E12;
    /* TR solver handle */
    _TRNSPBC_HANDLE_t handle;
    MKL_INT i;

    /* initial parameter values */
    std::vector<double> pars = GetParameterValuesVector();

    // intermediate function values
    std::vector<double> fvec(m_data_size,0);
    // jacobi matrix
    std::vector<double> fjac(m_data_size,0);

    /* lower and upper bounds */
    std::vector<double> lowerBounds = GetLowerBoundsVector();
    std::vector<double> upperBounds = GetUpperBoundsVector();


    /* initialize solver (allocate mamory, set initial values)
        handle  in/out: TR solver handle
        n       in:     number of function variables
        m       in:     dimension of function value
        x       in:     solution vector. contains values x for f(x)
        eps     in:     precisions for stop-criteria
        iter1   in:     maximum number of iterations
        iter2   in:     maximum number of iterations of calculation of trial-step
        rs      in:     initial step bound */
    if (dtrnlspbc_init (&handle, &n, &m, &pars[0], &lowerBounds[0], &upperBounds[0], &eps[0], &iter1, &iter2, &rs) !=
        TR_SUCCESS)
    {
        /* if function does not complete successful then print error message */
        DM::Debug("| error in dtrnlspbc_init\n");
      /* Release internal MKL memory that might be used for computations         */
      /* NOTE: It is important to call the routine below to avoid memory leaks   */
      /* unless you disable MKL Memory Manager                                   */
      MKL_FreeBuffers();
        /* and exit */
        return 0;
    }
    /* set initial rci cycle variables */
    RCI_Request = 0;
    successful = 0;
    /* rci cycle */
    while (successful == 0)
    {
        /* call tr solver
            handle      in/out: tr solver handle
            fvec        in:     vector
            fjac        in:     jacobi matrix
            RCI_request in/out: return number which denote next step for performing */
        if (dtrnlspbc_solve (&handle, &fvec[0], &fjac[0], &RCI_Request) != TR_SUCCESS)
        {
            /* if function does not complete successful then print error message */
            DM::Debug("| error in dtrnlspbc_solve\n");
            /* Release internal MKL memory that might be used for computations         */
            /* NOTE: It is important to call the routine below to avoid memory leaks   */
            /* unless you disable MKL Memory Manager                                   */
            MKL_FreeBuffers();
            /* and exit */
            return 0;
        }
        /* according with rci_request value we do next step */
        if (RCI_Request == -1 ||
            RCI_Request == -2 ||
            RCI_Request == -3 ||
            RCI_Request == -4 ||
            RCI_Request == -5 ||
            RCI_Request == -6)
            /* exit rci cycle */
            successful = 1;
        if (RCI_Request == 1)
        {
            /* recalculate function value
                m       in:     dimension of function value
                n       in:     number of function variables
                x       in:     solution vector
                fvec    out:    function value f(x) */
            MKL_Evaluate_function(&m, &n, &pars[0], &fvec[0]);
        }
        if (RCI_Request == 2)
        {
            /* compute jacobi matrix
                MKL_Evaluate_function   in:     external objective function
                n               in:     number of function variables
                m               in:     dimension of function value
                fjac            out:    jacobi matrix
                x               in:     solution vector
                jac_eps         in:     jacobi calculation precision */
            if (djacobix (&FitFunctionWrapper::FitFunction_MKL, &n, &m, &fjac[0], &pars[0], &eps[0], (void *)this) != TR_SUCCESS)
            {
                /* if function does not complete successful then print error message */
                DM::Debug("| error in djacobi\n");
                /* Release internal MKL memory that might be used for computations         */
                /* NOTE: It is important to call the routine below to avoid memory leaks   */
                /* unless you disable MKL Memory Manager                                   */
                MKL_FreeBuffers();
                /* and exit */
                return 0;
            }
        }
    }
    /* get solution statuses
        handle            in:   TR solver handle
        iter              out:  number of iterations
        st_cr             out:  number of stop criterion
        r1                out:  initial residuals
        r2                out:  final residuals */
    if (dtrnlspbc_get (&handle, &iter, &st_cr, &r1, &r2) != TR_SUCCESS)
    {
        /* if function does not complete successful then print error message */
        DM::Debug("| error in dtrnlspbc_get\n");
        /* Release internal MKL memory that might be used for computations         */
        /* NOTE: It is important to call the routine below to avoid memory leaks   */
        /* unless you disable MKL Memory Manager                                   */
        MKL_FreeBuffers();
        /* and exit */
        return 0;
    }
    /* free handle memory */
    if (dtrnlspbc_delete (&handle) != TR_SUCCESS)
    {
        /* if function does not complete successful then print error message */
        DM::Debug("| error in dtrnlspbc_delete\n");
        /* Release internal MKL memory that might be used for computations         */
        /* NOTE: It is important to call the routine below to avoid memory leaks   */
        /* unless you disable MKL Memory Manager                                   */
        MKL_FreeBuffers();
        /* and exit */
        return 0;
    }
    SetParameterValues(pars);

    //bool converged = !m_abort_fit
    fvec.clear();
    fjac.clear();
    eps.clear();
    pars.clear();
    lowerBounds.clear();
    upperBounds.clear();

    /* Release internal MKL memory that might be used for computations         */
    /* NOTE: It is important to call the routine below to avoid memory leaks   */
    /* unless you disable MKL Memory Manager                                   */
    MKL_FreeBuffers();

    // If we've gotten this far, then things worked.  Don't check the residual value as a means of convergence,
    //    because we should always allow poor fits (it's pretty hard for us to get really excellent fits on
    //    noisy data with not-completely-accurate models)
    return 1;
}

// MKL computes the Jacobian with this.  MKL uses OpenMP to compute this for multiple
//     elements of the Jacobian simultaneously.  DO NOT use the ComputeFunctionValues function,
//     as this would be changing the member temporary parameter values from multiple threads simultaneously.
//     This function should be much more thread-friendly.
void CFitFunction::MKL_Evaluate_function(MKL_INT *m, MKL_INT *n, double *x, double *f)
{
    for (MKL_INT i=0; i<*m; i++)
    {
        f[i]=(m_data_ptr[i]-FunctionValue(0,i,x))/m_errors_ptr[i];
    }
}

double CFitFunction::FunctionValue(ulong iValue, ulong i, double *pars)
{
    double N = pars[0];
    double mu = pars[1];
    double sigma = pars[2];
    double x;
    bool inROI;

    if (iValue != 0) return 0;
    inROI=CalculateCoordinate(i,x);
    if (!inROI && !GetCalculateOutsideROI()) return 0;

    double tmp = (x - mu)/sigma;
    return N * exp( - tmp*tmp);
}