#define wnlin_N         40

#define wnlin_P         3

/* initial guess should be chosen so that Jacobian has full rank,

 * or some solvers will fail */

static double wnlin_x0[wnlin_P] = { 1.0, 0.9, 0.0 };

static double wnlin_epsrel = 1.0e-7;

static int wnlin_internal_weight = 1;

/* data */

static double wnlin_Y[wnlin_N] = {

  6.08035e+00, 5.47552e+00, 5.94654e+00, 5.04920e+00, 4.78568e+00,

  3.51748e+00, 2.84671e+00, 3.24634e+00, 3.23395e+00, 3.30385e+00,

  2.83439e+00, 2.31891e+00, 2.33858e+00, 2.40559e+00, 2.41856e+00,

  1.99966e+00, 1.88127e+00, 1.91477e+00, 1.70415e+00, 1.60316e+00,

  1.77937e+00, 1.55302e+00, 1.50903e+00, 1.36364e+00, 1.36873e+00,

  1.41954e+00, 1.37778e+00, 1.23573e+00, 1.28524e+00, 1.46327e+00,

  1.22315e+00, 1.19330e+00, 1.18717e+00, 8.83172e-01, 1.23424e+00,

  1.14683e+00, 1.11091e+00, 1.20396e+00, 1.28722e+00, 1.05801e+00

};

/* weights */

static double wnlin_W[wnlin_N] = {

  2.77778e+00, 3.27690e+00, 3.85426e+00, 4.51906e+00, 5.28083e+00,

  6.14919e+00, 7.13370e+00, 8.24349e+00, 9.48703e+00, 1.08717e+01,

  1.24036e+01, 1.40869e+01, 1.59238e+01, 1.79142e+01, 2.00553e+01,

  2.23415e+01, 2.47646e+01, 2.73137e+01, 2.99753e+01, 3.27337e+01,

  3.55714e+01, 3.84696e+01, 4.14085e+01, 4.43678e+01, 4.73278e+01,

  5.02690e+01, 5.31731e+01, 5.60234e+01, 5.88046e+01, 6.15036e+01,

  6.41092e+01, 6.66121e+01, 6.90054e+01, 7.12839e+01, 7.34442e+01,

  7.54848e+01, 7.74053e+01, 7.92069e+01, 8.08918e+01, 8.24632e+01

};

static void

wnlin_checksol(const double x[], const double sumsq,

               const double epsrel, const char *sname,

               const char *pname)

{

  size_t i;

  const double sumsq_exact = 29.7481259665713758;

  const double wnlin_x[wnlin_P] = { 5.17378551196259195,

                                    0.111041758006851149,

                                    1.05282724070446099 };

  gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",

               sname, pname);

  for (i = 0; i < wnlin_P; ++i)

    {

      gsl_test_rel(x[i], wnlin_x[i], epsrel, "%s/%s i=%zu",

                   sname, pname, i);

    }

}

static int

wnlin_f (const gsl_vector *x, void *params, gsl_vector *f)

{

  int *iptr = (int *) params;

  int doweight = iptr ? *iptr : 0;

  double A = gsl_vector_get (x, 0);

  double lambda = gsl_vector_get (x, 1);

  double b = gsl_vector_get (x, 2);

  size_t i;

  /* model Yi = A * exp(-lambda * i) + b */

  for (i = 0; i < wnlin_N; i++)

    {

      double ti = i;

      double yi = wnlin_Y[i];

      double swi = sqrt(wnlin_W[i]);

      double Mi = A * exp (-lambda * ti) + b;

      if (doweight)

        gsl_vector_set (f, i, swi * (Mi - yi));

      else

        gsl_vector_set (f, i, Mi - yi);

    }

  return GSL_SUCCESS;

}

static int

wnlin_df (const gsl_vector *x, void *params, gsl_matrix *df)

{

  int *iptr = (int *) params;

  int doweight = iptr ? *iptr : 0;

  double A = gsl_vector_get (x, 0);

  double lambda = gsl_vector_get (x, 1);

  size_t i;

  for (i = 0; i < wnlin_N; i++)

    {

      gsl_vector_view v = gsl_matrix_row(df, i);

      double ti = i;

      double swi = sqrt(wnlin_W[i]);

      double e = exp(-lambda * ti);

      gsl_vector_set(&v.vector, 0, e);

      gsl_vector_set(&v.vector, 1, -ti * A * e);

      gsl_vector_set(&v.vector, 2, 1.0);

      if (doweight)

        gsl_vector_scale(&v.vector, swi);

    }

  return GSL_SUCCESS;

}

static int

wnlin_fvv (const gsl_vector * x, const gsl_vector * v,

           void *params, gsl_vector * fvv)

{

  int *iptr = (int *) params;

  int doweight = iptr ? *iptr : 0;

  double A = gsl_vector_get (x, 0);

  double lambda = gsl_vector_get (x, 1);

  double v1 = gsl_vector_get(v, 0);

  double v2 = gsl_vector_get(v, 1);

  size_t i;

  for (i = 0; i < wnlin_N; i++)

    {

      double ti = i;

      double swi = sqrt(wnlin_W[i]);

      double fvvi;

      fvvi = exp(-ti*lambda)*ti*v2 * (-2*v1 + ti*v2*A);

      if (doweight)

        gsl_vector_set(fvv, i, swi * fvvi);

      else

        gsl_vector_set(fvv, i, fvvi);

    }

  return GSL_SUCCESS;

}

static gsl_multifit_nlinear_fdf wnlin_func1 =

{

  wnlin_f,

  wnlin_df,

  wnlin_fvv,

  wnlin_N,

  wnlin_P,

  (void *) &wnlin_internal_weight,

  0,

  0,

  0

};

static gsl_multifit_nlinear_fdf wnlin_func2 =

{

  wnlin_f,

  wnlin_df,

  wnlin_fvv,

  wnlin_N,

  wnlin_P,

  NULL,

  0,

  0,

  0

};

static test_fdf_problem wnlin_problem1 =

{

  "wnlin_internal_weights",

  wnlin_x0,

  NULL,

  NULL,

  &wnlin_epsrel,

  &wnlin_checksol,

  &wnlin_func1

};

static test_fdf_problem wnlin_problem2 =

{

  "wnlin_external_weights",

  wnlin_x0,

  wnlin_W,

  NULL,

  &wnlin_epsrel,

  &wnlin_checksol,

  &wnlin_func2

};