/*!
 * \file lmfit_test.cpp
 * \author Konrad Werys
 * \date 2019/08/08
 */

#include "CmakeConfigForTomato.h"
#ifdef USE_LMFIT

#include "gtest/gtest.h"
#include "OxTestData.h"

#include "lmmin.h"

void evaluate_nonlin1(
        const double *p, int n, const void *data, double *f, int *info )
{
    f[0] = p[0]*p[0] + p[1]*p[1] - 1; /* unit circle    x^2+y^2=1 */
    f[1] = p[1] - p[0]*p[0];          /* standard parabola  y=x^2 */
}

TEST(lmfit, example_copied_from_lmfit_code) {

    bool verbose = false;

    int n = 2;   /* dimension of the problem */
    double p[2]; /* parameter vector p=(x,y) */

    /* auxiliary parameters */
    lm_control_struct control = lm_control_double;
    lm_status_struct  status;
    if (verbose) {
        control.verbosity = 3;
    }

    p[0] = 2;
    p[1] = 2;

    /* the minimization */
    if (verbose) {
        printf("Minimization:\n");
    }

    lmmin(n, p, n, NULL, NULL, evaluate_nonlin1, &control, &status);

    if (verbose) {
        /* print results */
        printf("\n");
        printf("lmmin status after %d function evaluations:\n  %s\n",
               status.nfev, lm_infmsg[status.outcome]);

        printf("\n");
        printf("Solution:\n");
        printf("  x = %19.11f\n", p[0]);
        printf("  y = %19.11f\n", p[1]);
        printf("  d = %19.11f => ", status.fnorm);

        /* convergence of lmfit is not enough to ensure validity of the solution */
        if (status.fnorm >= control.ftol)
            printf("not a valid solution, try other starting values\n");
        else
            printf("valid, though not the only solution: "
                   "try other starting values\n");
    }
    EXPECT_NEAR(p[0], 0.786, 1e-2);
    EXPECT_NEAR(p[1], 0.618, 1e-2);
}



void evaluate_T1(const double *parameters, int nSamples, const void *data_void, double *residuals, int *info )
{

    double A = parameters[0];
    double B = parameters[1];
    double T1star = parameters[2];

    // casting the pointer data_void
    double *data = (double*)data_void;

    for (unsigned int i = 0; i < nSamples; i++) {
        // decoding the data from the data_void pointer
        double invTime = data[i];
        double measured = data[i+nSamples];
        double calculated = 0;

        calculated = A - B * exp( -invTime/ T1star );

        residuals[i] = calculated - measured;
    }
}

TEST(lmfit, simpleT1_test) {

    bool verbose = false;

    typedef double TYPE;

    char filePath [] = "testData/myocardium.yaml";
    Ox::TestData<TYPE> testData(filePath);
    int nSamples = testData.getNSamples();

    // you can pass a single data pointer to lmfit models, here I use a single array
    double data[14];
    KWUtil::copyArrayToArray(nSamples, data, testData.getInvTimesPtr());
    KWUtil::copyArrayToArray(nSamples, &data[7], testData.getSignalPtr());

    int n = 3;   /* dimension of the problem */
    double parameters[3]; /* parameter vector p=(x,y) */

    /* auxiliary parameters */
    lm_control_struct control = lm_control_double;
    lm_status_struct  status;
    if (verbose) {
        control.verbosity = 3;
    }

    // starting point
    parameters[0] = 100;
    parameters[1] = 200;
    parameters[2] = 1000;

    /* the minimization */
    if (verbose){
        printf("Minimization:\n");
    }

    lmmin(n, parameters, nSamples, NULL, data, evaluate_T1, &control, &status);

    if (verbose) {
        /* print results */
        printf("\n");
        printf("lmmin status after %d function evaluations:\n  %s\n",
               status.nfev, lm_infmsg[status.outcome]);

        printf("\n");
        printf("Solution:\n");
        printf("  x = %19.11f\n", parameters[0]);
        printf("  y = %19.11f\n", parameters[1]);
        printf("  y = %19.11f\n", parameters[2]);
        printf("  d = %19.11f => ", status.fnorm);

        /* convergence of lmfit is not enough to ensure validity of the solution */
        if (status.fnorm >= control.ftol)
            printf("not a valid solution, try other starting values\n");
        else
            printf("valid, though not the only solution: "
                   "try other starting values\n");

        std::cout << "control.patience: " << control.patience << std::endl;
    }

    EXPECT_NEAR(parameters[0], testData.getResultsMolli()[0], 1e-2);
    EXPECT_NEAR(parameters[1], testData.getResultsMolli()[1], 1e-2);
    EXPECT_NEAR(parameters[2], testData.getResultsMolli()[2], 1e-2);
}

#endif // USE_LMFIT