/* nag::ad::e04us Passive Example Program.
 */

#include <dco.hpp>
#include <iostream>
#include <nagad.h>

// Function which calls NAG AD Library routines.
// Solve the sum-of-squares minimization problem
//    F(x) = .5 sum_{i=0}^{44} (y_i - f_i(x))^2
// for
//    f_i(x) = x_1 + (0.49 - x_1) * e^{-x_2 * (a_i - 8)}
//
// and a_i is given via ruser parameter.
template <typename T> void func(std::vector<T> &ruser, std::vector<T> &x);

int main()
{

  std::cout << " nag::ad::e04us Passive Example Program Results\n";

  // (Initial) Solution
  std::vector<double> x{0.4, 0.0};

  // Parameters a
  std::vector<double> ruser{
      8.0,  8.0,  10.0, 10.0, 10.0, 10.0, 12.0, 12.0, 12.0, 12.0, 14.0,
      14.0, 14.0, 16.0, 16.0, 16.0, 18.0, 18.0, 20.0, 20.0, 20.0, 22.0,
      22.0, 22.0, 24.0, 24.0, 24.0, 26.0, 26.0, 26.0, 28.0, 28.0, 30.0,
      30.0, 30.0, 32.0, 32.0, 34.0, 36.0, 36.0, 38.0, 38.0, 40.0, 42.0};

  func(ruser, x);

  // Print outputs
  std::cout.setf(std::ios::scientific, std::ios::floatfield);
  std::cout.precision(12);
  std::cout << "\n Solution:\n";
  std::cout << " x = [" << x[0] << ", " << x[1] << "]" << std::endl;
  return 0;
}

// Driver with the passive calls.
// Solve the sum-of-squares minimization problem
//    F(x) = .5 sum_{i=0}^{44} (y_i - f_i(x))^2
// for
//    f_i(x) = x_1 + (0.49 - x_1) * e^{-x_2 * (a_i - 8)}
//
// and a_i is given via ruser parameter.
template <typename T> void func(std::vector<T> &ruser, std::vector<T> &x)
{
  // Problem sizes
  Integer m = 44, n = 2, nclin = 1, ncnln = 1, lb = n + nclin + ncnln,
          liwork = 3 * n + nclin + 2 * ncnln, lda = nclin, sda = n,
          ldcj = ncnln, ldfj = m, ldr = n, lwork = 20 * n + m * (n + 3);

  lwork += 2 * n * n + 11 * nclin;
  lwork += n * nclin + (2 * n + 21) * ncnln;

  std::vector<T> c(ncnln), cjac(ncnln * n), f(m), fjac(m * n), clamda(lb),
      r(ldr * n), rwsav(475), work(lwork);
  std::vector<Integer> iwsav(610), iwork(liwork), istate(lb);
  std::vector<logical> lwsav(120);

  // Problem data
  std::vector<T> a{1.0, 1.0};
  std::vector<T> y{0.49, 0.49, 0.48, 0.47, 0.48, 0.47, 0.46, 0.46, 0.45,
                   0.43, 0.45, 0.43, 0.43, 0.44, 0.43, 0.43, 0.46, 0.45,
                   0.42, 0.42, 0.43, 0.41, 0.41, 0.40, 0.42, 0.40, 0.40,
                   0.41, 0.40, 0.41, 0.41, 0.40, 0.40, 0.40, 0.38, 0.41,
                   0.40, 0.40, 0.41, 0.38, 0.40, 0.40, 0.39, 0.39};
  std::vector<T> bl{0.4, -4.0, 1.0, 0.0};
  std::vector<T> bu{1.0E+25, 1.0E+25, 1.0E+25, 1.0E+25};

  // Initialize sav arrays
  Integer ifail = 0;
  char    cwsav;
  nag::ad::e04wb("E04USA", &cwsav, 1, lwsav.data(), lwsav.size(), iwsav.data(),
                 iwsav.size(), rwsav.data(), rwsav.size(), ifail);

  auto objfun = [&](nag::ad::handle_t &ad_handle,
                Integer &          mode,
                const Integer &    m,
                const Integer &    n,
                const Integer &    ldfj,
                const Integer &    needfi,
                const T            x[],
                T                  f[],
                T                  fjac[],
                const Integer &    nstate)
              {
                T x1{x[0]}, x2{x[1]};

              #define FJAC(i, j) fjac[(i) + (j)*m]

                if (mode == 0 && needfi > 0)
                {
                  // Update only one element f(i)
                  Integer i  = needfi - 1;
                  T       ai = ruser[i] - 8.0;
                  f[i]       = x1 + (0.49 - x1) * exp(-x2 * ai);
                }
                else
                {
                  for (int i = 0; i < m; ++i)
                  {
                    T ai = ruser[i] - 8.0;
                    if (mode == 0 || mode == 2)
                    {
                      // Update the whole of f(i)
                      f[i] = x1 + (0.49 - x1) * exp(-x2 * ai);
                    }
                    if (mode == 1 || mode == 2)
                    {
                      // Update df/dx1
                      FJAC(i, 0) = 1.0 - exp(-x2 * ai);
                      // Update df/dx2
                      FJAC(i, 1) = -(0.49 - x1) * ai * exp(-x2 * ai);
                    }
                  }
                }
              #undef FJAC
              };

  auto confun = [&](nag::ad::handle_t &ad_handle,
                  Integer &          mode,
                  const Integer &    ncnln,
                  const Integer &    n,
                  const Integer &    ldcj,
                  const Integer      needc[],
                  const T            x[],
                  T                  c[],
                  T                  cjac[],
                  const Integer &    nstate)
                {
                  // calling confun for the first time
                  if (nstate == 1)
                  {
                    for (int i = 0; i < ncnln * n; ++i)
                    {
                      cjac[i] = 0.0;
                    }
                  }
                  // evaluate c and cjac (depending on mode) at point x
                  if (needc[0] > 0)
                  {
                    // c required
                    if (mode == 0 || mode == 2)
                    {
                      c[0] = -0.09 - x[0] * x[1] + 0.49 * x[1];
                    }
                    // cjac required
                    if (mode == 1 || mode == 2)
                    {
                      cjac[0]     = -x[1];
                      cjac[ncnln] = -x[0] + 0.49;
                    }
                  }
                };

  // Create AD configuration data object
  nag::ad::handle_t ad_handle;
  ifail = 0;
  // Solve the problem
  Integer iter;
  T       objf;
  ifail = 0;
  nag::ad::e04us(ad_handle, m, n, nclin, ncnln, lda, ldcj, ldfj, ldr, a.data(),
                 bl.data(), bu.data(), y.data(), confun, objfun, iter,
                 istate.data(), c.data(), cjac.data(), f.data(), fjac.data(),
                 clamda.data(), objf, r.data(), x.data(), iwork.data(), liwork,
                 work.data(), lwork, lwsav.data(), iwsav.data(), rwsav.data(), ifail);
}