#pragma once
#define STRICT_R_HEADERS
#ifndef __rxode2_H__
#define __rxode2_H__

#define rxLlikSaveSize 9

#include <stdbool.h>
#include "rxode2parse.h"
#include <R.h>
#include <Rinternals.h>
#include <Rversion.h>
#include <Rmath.h>
#include <R_ext/Rdynload.h>

#include "rxode2parseSbuf.h"

#include <float.h>
#include <stdio.h>
#include <stdarg.h>


#define isDose(evid) ((evid) == 3 || (evid) >= 100)
#define isObs(evid) ((evid) == 0 || (evid) == 2 || ((evid) >= 9 && (evid) <= 99))

#include "rxode2_control.h"
#include <stdint.h>    // for uint64_t rather than unsigned long long

#ifndef __RXODE2PTR_H__  // directly refer to abi need to be excluded
// Effective neq for the SOLVE LOOP and per-event SOLVE BUFFER stride.
// Returns ind->neqOverride when set in [0, op->neq], else op->neq.
// Allocations remain sized for op->neq; the override is always <= op->neq
// (caller's contract).  Used by getSolve()/getAdvan() so that a per-
// individual pred-mode solve writes and reads at the same compact stride
// without mutating shared op->neq from a parallel worker thread.
// NOTE: getAdvan() + neqOverride is unsupported when op->numLin > 0
// (op->linOffset is computed from the full neq layout).  A linCmt() model
// mixed with ODEs does reach nlmixr2est's FOCEi flow with numLin > 0 on the
// predNoLhs model, so this is NOT covered by the model itself; nlmixr2est
// translates those models with linToOde() before estimating (its issue #286),
// which leaves numLin == 0 here.
static inline int rxEffNeq(const rx_solving_options_ind *ind,
                           const rx_solving_options *op) {
  int o = ind->neqOverride;
  return (o >= 0 && o <= op->neq) ? o : op->neq;
}
#define getSolve(idx) ind->solve + (rxEffNeq(ind, op))*(idx)
#define getAdvan(idx) ind->solve + op->linOffset + (rxEffNeq(ind, op))*(idx)
#endif

#ifndef max2
#define max2( a , b )  ( (a) > (b) ? (a) : (b) )
#endif
#define isSameTime(xout, xp) (fabs((xout)-(xp))  <= 2.0*DBL_EPSILON*max2(fabs(xout),fabs(xp)))
// use ~dop853 definition of same time
#define isSameTimeDop(xout, xp) (0.1 * fabs((xout)-(xp)) <= fabs(xout) * 2.3E-16)

#ifdef _isrxode2_

#else

#if defined(__cplusplus)
#include "rxode2_RcppExports.h"

#endif

#endif // _isrxode2_

#define rxNEG_LOG_SQRT_TWO_PI

#if defined(__cplusplus)
extern "C" {
#endif


#ifdef _isrxode2_
  void parseFree(int last);
#endif

typedef void (*t_dydt)(int *neq, double t, double *A, double *DADT);
typedef void (*t_calc_jac)(int *neq, double t, double *A, double *JAC, unsigned int __NROWPD__);
typedef void (*t_calc_lhs)(int cSub, double t, double *A, double *lhs);
typedef void (*t_update_inis)(int cSub, double *);
typedef void (*t_dydt_lsoda_dum)(int *neq, double *t, double *A, double *DADT);
typedef void (*t_jdum_lsoda)(int *neq, double *t, double *A,int *ml, int *mu, double *JAC, int *nrowpd);
typedef int (*t_dydt_liblsoda)(double t, double *y, double *ydot, void *data);
typedef void (*t_ode_current)(void);

typedef void (*t_set_solve)(rx_solve *);

typedef rx_solve *(*t_get_solve)(void);

typedef void *(*t_assignFuns)(void);

// Adjoint objective-gradient backward sweep (src/adjoint.cpp).  Its address is
// exported to downstream packages through the rxode2 function-pointer table
// (see _rxode2_rxode2Ptr in src/init.c and rxode2ptr.h); this direct declaration
// is used only when building rxode2 itself (guarded off once rxode2ptr.h has
// redeclared the name as a function pointer).
#ifndef __RXODE2PTR_H__
void rxode2AdjointSweep(double *tg, double *J, double *dP, double *cover,
                        int *obsK, int ns, int np, int nt, int nobs, double *out,
                        int nCj, int *cjK, int *cjCmt, double *cjAlpha,
                        int nDual, int *dualK, double *dualW, double *dualC);
void rxode2AdjointTrajSweep(double *tg, double *J, double *dP, int ns, int np,
                            int nt, int *outK, int nOut, int *stateIdx,
                            int nStates, double *result, int nCj, int *cjK,
                            int *cjCmt, double *cjAlpha, int nDual, int *dualK,
                            double *dualW, double *dualC);
rx_solve *getRxSolve_(void);
void rxSetSolveAtolRtol(double atol, double rtol);
void rxGetSolveAtolRtol(double *atol, double *rtol);
int getIndCmt(rx_solving_options* op, rx_solving_options_ind* ind, int kk);
void setIndCmt(rx_solving_options* op, rx_solving_options_ind* ind, int kk, int cmt);
// Event ("jump") sensitivity shape (src/par_solve.cpp); see the block comment
// there for how a downstream package swaps shapes around a batch of solves.
int rxode2EventSensShapeSize(void);
int rxode2EventSensShapeSave(void *buf, int bufSize);
int rxode2EventSensShapeRestore(const void *buf, int bufSize);
void rxode2EventSensLoadFull(SEXP trans, int active, int nState, int nParam,
                             int nParam2, int nParam3, int useCalcJac);
void rxode2EventSensGetDims(int *active, int *nState, int *nParam, int *nParam2,
                            int *nParam3, int *useCalcJac);
void rxode2EventSensSetDims(int active, int nState, int nParam, int nParam2,
                            int nParam3, int useCalcJac);
void rxode2EventSensSetActive(int active);
void rxode2EventSensDeactivate(void);
#endif
rx_solve *getRxSolve2_(void);
rx_solve *getRxSolve(SEXP ptr);

#ifndef __RXODE2PTR_H__
void par_solve(rx_solve *rx);
#endif

rx_solving_options *getRxOp(rx_solve *rx);

SEXP rxode2_par_df(void);

void rxOptionsIniEnsure(int mx, int cores);

void rxUpdateFuns(SEXP trans);

#define _eps sqrt(DBL_EPSILON)

static inline double erfinv(double x)  __attribute__((unused));
static inline double erfinv(double x) {
  return Rf_qnorm5((1 + x)/2.0, 0, 1, 1, 0)*M_SQRT1_2;
}
#define rxDistributionNorm     1
#define rxDistributionPois     2
#define rxDistributionBinom    3
#define rxDistributionBeta     4
#define rxDistributionT        5
#define rxDistributionChisq    6
#define rxDistributionDexp     7
#define rxDistributionF        8
#define rxDistributionGeom     9
#define rxDistributionHyper   10
#define rxDistributionUnif    11
#define rxDistributionWeibull 12
#define rxDistributionCauchy  13
#define rxDistributionGamma   14
#define rxDistributionOrdinal 15
#define rxDistributionN2ll    16
#define rxDistributionDnorm   17

static inline void _splitYj(int *yj, int *dist,  int *trans) {
  *dist  = *yj/10;
  *trans = *yj - *dist*10;
  *dist  = *dist + 1;
}
// Inverse
static inline double _powerDi(double x, double lambda, int yj0, double low, double high)  __attribute__((unused));
static inline double _powerDi(double x, double lambda, int yj0, double low, double high){
  if (!R_finite(x)) return NA_REAL;
  double x0=x, ret, l2, yjd;
  int yj, dist;
  _splitYj(&yj0, &dist,  &yj);
  switch(yj){
  case 7: // inverse-Yeo Johnson followed by pnorm
    yjd = _powerDi(x, lambda, 1, low, high);
    return (high-low)*Rf_pnorm5(yjd, 0, 1, 1, 0)+low;
  case 6: // probitInverse
    return (high-low)*Rf_pnorm5(x, 0, 1, 1, 0)+low;
  case 5: // inverse-Yeo-Johnson followed by expit
    yjd = _powerDi(x, lambda, 1, low, high);
    return (high-low)/(1+exp(-yjd))+low;
  case 4:
    return (high-low)/(1+exp(-x))+low; // expit
  case 3:
    return exp(x);
  case 2:
    return x;
  case 0:
    if (lambda == 1.0) return (x+1.0);
    if (lambda == 0) return exp(x);
    // (x^lambda-1)/lambda=y
    // (lambda*y+1)^(1/lambda)
    x0 = x*lambda+1.0;
    if (x0 <= _eps) return _eps;
    ret = pow(x0, 1.0/lambda);
    if (ISNA(ret)) {
      // Warning?
      return _eps;
    }
    return ret;
  case 1:
    if (lambda == 1.0) return x;
    if (x >= 0){
      // log(x+1)= y; exp(y)-1=x
      if (lambda == 0) return expm1(x);
      // ((x+1)^lambda-1)/lambda=y
      // (y*lambda+1)^(1/y)-1=y
      return pow(x*lambda+1.0, 1.0/lambda)-1.0;
    } else {
      // (-(1-x)^(2-lambda)-1)/(2-lambda)
      if (lambda ==  2.0) return -expm1(-x);
      // (-(1-x)^(2-lambda)-1)/(2-lambda) = y
      l2 = (2.0 - lambda);
      return 1.0 - pow(1.0 - l2*x, 1.0/l2);
    }
  }
  return NA_REAL;
}

static inline double _powerD(double x, double lambda, int yj0, double low, double high)  __attribute__((unused));
static inline double _powerD(double x, double lambda, int yj0, double low, double high) {
  if (!R_finite(x)) return NA_REAL;
  double x0=x, l2, p;
  int yj, dist;
  _splitYj(&yj0, &dist,  &yj);
  switch (yj) {
  case 7:
    p = (x-low)/(high-low);
    if (p >= 1) return R_NaN;
    if (p <= 0) return R_NaN;
    /* REprintf("%f %f %f\n", x, p, -log(1/p-1)); */
    p = Rf_qnorm5(p, 0, 1, 1, 0);
    if (lambda == 1.0) return p;
    if (p >= 0){
      if (lambda == 0) return log1p(p);
      return (pow(p + 1.0, lambda) - 1.0)/lambda;
    } else {
      if (lambda == 2.0) return -log1p(-p);
      l2 = 2.0 - lambda;
      return (1.0 - pow(1.0 - p, l2))/l2;
    }
  case 6: // probitNorm
    p = (x-low)/(high-low);
    if (p >= 1) return R_NaN;
    if (p <= 0) return R_NaN;
    /* REprintf("%f %f %f\n", x, p, -log(1/p-1)); */
    return Rf_qnorm5(p, 0, 1, 1, 0);
  case 5: // logit followed by yeo-johnson
    p = (x-low)/(high-low);
    if (p >= 1) return R_NaN;
    if (p <= 0) return R_NaN;
    p = -log(1/p-1);
    if (lambda == 1.0) return p;
    if (p >= 0){
      if (lambda == 0) return log1p(p);
      return (pow(p + 1.0, lambda) - 1.0)/lambda;
    } else {
      if (lambda == 2.0) return -log1p(-p);
      l2 = 2.0 - lambda;
      return (1.0 - pow(1.0 - p, l2))/l2;
    }
  case 4: // logitNorm
    p = (x-low)/(high-low);
    if (p >= 1) return R_NaN;
    if (p <= 0) return R_NaN;
    /* REprintf("%f %f %f\n", x, p, -log(1/p-1)); */
    return -log(1/p-1);
  case 3: // logNorm
    if (x <= _eps) x0= _eps;
    return log(x0);
  case 2: // norm
    return x;
  case 0: // boxCoxNorm
    if (lambda == 1.0) return x-1.0;
    if (x <= _eps) x0= _eps;
    if (lambda ==  0.0) return log(x0);
    return (pow(x0, lambda) - 1.0)/lambda;
  case 1: // yeoJohnsonNorm
    if (lambda == 1.0) return x;
    if (x >= 0){
      if (lambda == 0) return log1p(x);
      return (pow(x + 1.0, lambda) - 1.0)/lambda;
    } else {
      if (lambda == 2.0) return -log1p(-x);
      l2 = 2.0 - lambda;
      return (1.0 - pow(1.0 - x, l2))/l2;
    }
  }
  return NA_REAL;
}

static inline double _powerDD(double x, double lambda, int yj0, double low, double high)  __attribute__((unused));
static inline double _powerDD(double x, double lambda, int yj0, double low, double high){
  if (!R_finite(x)) return NA_REAL;
  double x0 = x, xl, hl,eri;
  int yj, dist;
  _splitYj(&yj0, &dist,  &yj);
  switch(yj){
  case 7:
    // Subs(Derivative(yeoJohnson(_xi_1), _xi_1), (_xi_1), (logit(x)))*Derivative(logit(x), x)
    return _powerDD(_powerD(x, lambda, 6, low, high), lambda, 1, low, high)*_powerDD(x, lambda, 6, low, high);
  case 6: // derivative
    // 2.82842712474619*M_SQRT_PI/2*exp((erfinv(-1+2*(-low+x)/(high-low)))^2)/(high-low)
    hl = (high-low);
    eri = erfinv(-1+2*(-low+x)/hl);
    return 2.506628274631000241612*exp(eri*eri)/hl;
  case 5: // logit followed by yeo-johnson  yeoJohnson(logit(x))
    // Subs(Derivative(yeoJohnson(_xi_1), _xi_1), (_xi_1), (logit(x)))*Derivative(logit(x), x)
    return _powerDD(_powerD(x, lambda, 4, low, high), lambda, 1, low, high)*_powerDD(x, lambda, 4, low, high);
  case 4: // logitNorm
    xl = (x-low);
    hl = (high - low);
    return hl/(xl*(hl-xl));
  case 3:
    if (x <= _eps) x0 = _eps;
    return 1/x0;
  case 2:
    return 1.0;
  case 0:
    if (lambda == 1.0) return 1.0;
    if (x <= _eps) x0 = _eps;
    if (lambda == 0.0) return 1/x0;
    // pow(x,lambda)/lambda - 1/lambda
    return pow(x0, lambda-1);
  case 1: // Yeo Johnson derivative
    if (lambda ==  1.0) return 1.0;
    if (x >= 0){
      if (lambda == 0.0) return 1.0/(x + 1.0);
      return pow(x + 1.0, lambda-1.0);
    } else {
      if (lambda == 2.0) return 1/(1.0 - x);
      return pow(1.0 - x, 1.0-lambda);
    }
  }
  return NA_REAL;
}

static inline double _powerDDD(double x, double lambda, int yj0, double low, double high) __attribute__((unused));
static inline double _powerDDD(double x, double lambda, int yj0, double low, double high){
  if (!R_finite(x)) return NA_REAL;
  double x0 = x, hl, xl,  t1, dL, eri;
  int yj, dist;
  _splitYj(&yj0, &dist,  &yj);
  switch(yj){
  case 7:
    // yj''(probit(x))*probit'(x)^2 + yj'(probit(x))*probit''(x)
    dL = _powerDD(x, lambda, 6, low, high);
    t1 = _powerD(x, lambda, 6, low, high);
    return dL*dL*_powerDDD(t1, lambda, 1, low, high) +
      _powerDD(t1, lambda, 1, low, high)*_powerDDD(x, lambda, 6, low, high);
  case 6: // derivative
    //10.026513098524*exp(erfinv(-1+2*(-low+x)/(high-low))^2)*M_SQRT_PI/2*exp((erfinv(-1+2*(-low+x)/(high-low)))^2)*erfinv(-1+2*(-low+x)/(high-low))/(high-low)^2
    hl = (high-low);
    eri = erfinv(-1+2*(-low+x)/hl);
    return 8.885765876316728650863*exp(2*eri*eri)*eri/(hl*hl);
  case 5:
    // yj''(logit(x))*logit'(x)^2 + yj'(logit(x))*logit''(x)
    dL = _powerDD(x, lambda, 4, low, high);
    t1 = _powerD(x, lambda, 4, low, high);
    return dL*dL*_powerDDD(t1, lambda, 1, low, high) +
      _powerDD(t1, lambda, 1, low, high)*_powerDDD(x, lambda, 4, low, high);
  case 4: // logit
    // d2/dx2 logit((x-low)/(high-low)) = 1/(high-x)^2 - 1/(x-low)^2
    hl = (high - low);
    xl = (-low + x);
    t1 = hl - xl;
    return 1.0/(t1*t1) - 1.0/(xl*xl);
  case 3:
    if (x <= _eps) x0 = _eps;
    return -1/(x0*x0);
  case 2:
    return 0;
  case 0:
    if (lambda == 1.0) return 0;
    if (x <= _eps) x0 = _eps;
    if (lambda == 0.0) return -1/(x0*x0);
    // pow(x,lambda)/lambda - 1/lambda
    return (lambda-1)*pow(x0, lambda-2);
  case 1:
    if (lambda == 1.0) return 0;
    if (x >= 0){
      if (lambda ==  0.0) return -1/((x + 1.0)*(x + 1.0));
      return (lambda-1.0)*pow(x + 1.0, lambda-2.0);
    } else {
      if (lambda == 2.0) return 1/((1.0 - x)*(1.0 - x));
      return -(1.0-lambda)*pow(1.0 - x, -lambda);
    }
  }
  return NA_REAL;
}

// Box-Cox VALUE lambda-derivatives at base b (>0), exponent lam.  y=(b^lam-1)/lam;
// near lam=0, y ~ ln b + lam*(ln b)^2/2 + lam^2*(ln b)^3/6 gives the removable-singularity limits.
static inline double _bcDL(double b, double lam) __attribute__((unused));
static inline double _bcDL(double b, double lam){ // dy/dlam
  if (b <= _eps) b = _eps;
  double lb = log(b);
  if (lam == 0.0) return 0.5*lb*lb;
  double bl = pow(b, lam);
  return bl*lb/lam - (bl - 1.0)/(lam*lam);
}
static inline double _bcDLx(double b, double lam) __attribute__((unused));
static inline double _bcDLx(double b, double lam){ // d2y/(dlam db)
  if (b <= _eps) b = _eps;
  return pow(b, lam - 1.0)*log(b);
}
static inline double _bcDL2(double b, double lam) __attribute__((unused));
static inline double _bcDL2(double b, double lam){ // d2y/dlam2
  if (b <= _eps) b = _eps;
  double lb = log(b);
  if (lam == 0.0) return lb*lb*lb/3.0;
  double bl = pow(b, lam);
  return bl*lb*lb/lam - 2.0*bl*lb/(lam*lam) + 2.0*(bl - 1.0)/(lam*lam*lam);
}
// Yeo-Johnson VALUE lambda-derivatives (x>=0: base x+1, exp lam; x<0: base 1-x, exp 2-lam,
// with the extra chain-rule signs from d(2-lam)/dlam=-1 and d(1-x)/dx=-1).
static inline double _yjDL(double x, double lam) __attribute__((unused));
static inline double _yjDL(double x, double lam){
  if (x >= 0) return _bcDL(x + 1.0, lam);
  return _bcDL(1.0 - x, 2.0 - lam);
}
static inline double _yjDLx(double x, double lam) __attribute__((unused));
static inline double _yjDLx(double x, double lam){
  if (x >= 0) return _bcDLx(x + 1.0, lam);
  return -_bcDLx(1.0 - x, 2.0 - lam);
}
static inline double _yjDL2(double x, double lam) __attribute__((unused));
static inline double _yjDL2(double x, double lam){
  if (x >= 0) return _bcDL2(x + 1.0, lam);
  return -_bcDL2(1.0 - x, 2.0 - lam);
}

// d(tbs)/dlambda -- the VALUE lambda-derivative of the both-sides transform (rxTBSdL).
static inline double _powerDLambda(double x, double lambda, int yj0, double low, double high) __attribute__((unused));
static inline double _powerDLambda(double x, double lambda, int yj0, double low, double high){
  if (!R_finite(x)) return NA_REAL;
  double p; int yj, dist;
  _splitYj(&yj0, &dist, &yj);
  switch(yj){
  case 0: return _bcDL(x, lambda);                     // boxCox
  case 1: return _yjDL(x, lambda);                     // yeoJohnson
  case 2: case 3: case 4: case 6: return 0.0;          // norm/logNorm/logit/probit: no lambda
  case 5: // logit then yeoJohnson
    p = (x-low)/(high-low); if (p >= 1 || p <= 0) return R_NaN; p = -log(1/p-1);
    return _yjDL(p, lambda);
  case 7: // probit then yeoJohnson
    p = (x-low)/(high-low); if (p >= 1 || p <= 0) return R_NaN; p = Rf_qnorm5(p, 0, 1, 1, 0);
    return _yjDL(p, lambda);
  }
  return NA_REAL;
}
// d2(tbs)/dlambda2 (rxTBSdL2).
static inline double _powerDLambda2(double x, double lambda, int yj0, double low, double high) __attribute__((unused));
static inline double _powerDLambda2(double x, double lambda, int yj0, double low, double high){
  if (!R_finite(x)) return NA_REAL;
  double p; int yj, dist;
  _splitYj(&yj0, &dist, &yj);
  switch(yj){
  case 0: return _bcDL2(x, lambda);
  case 1: return _yjDL2(x, lambda);
  case 2: case 3: case 4: case 6: return 0.0;
  case 5:
    p = (x-low)/(high-low); if (p >= 1 || p <= 0) return R_NaN; p = -log(1/p-1);
    return _yjDL2(p, lambda);
  case 7:
    p = (x-low)/(high-low); if (p >= 1 || p <= 0) return R_NaN; p = Rf_qnorm5(p, 0, 1, 1, 0);
    return _yjDL2(p, lambda);
  }
  return NA_REAL;
}
// d2(tbs)/(dlambda dx) (rxTBSdLx).  For the logit/probit-composed cases the inner transform p
// is lambda-free, so d/dx chains through dp/dx = _powerDD(x, lambda, <innerYj>, low, high).
static inline double _powerDLambdaX(double x, double lambda, int yj0, double low, double high) __attribute__((unused));
static inline double _powerDLambdaX(double x, double lambda, int yj0, double low, double high){
  if (!R_finite(x)) return NA_REAL;
  double p; int yj, dist;
  _splitYj(&yj0, &dist, &yj);
  switch(yj){
  case 0: return _bcDLx(x, lambda);
  case 1: return _yjDLx(x, lambda);
  case 2: case 3: case 4: case 6: return 0.0;
  case 5:
    p = (x-low)/(high-low); if (p >= 1 || p <= 0) return R_NaN; p = -log(1/p-1);
    return _yjDLx(p, lambda)*_powerDD(x, lambda, 4, low, high);
  case 7:
    p = (x-low)/(high-low); if (p >= 1 || p <= 0) return R_NaN; p = Rf_qnorm5(p, 0, 1, 1, 0);
    return _yjDLx(p, lambda)*_powerDD(x, lambda, 6, low, high);
  }
  return NA_REAL;
}

static inline double _powerL(double x, double lambda, int yj0, double low, double high) __attribute__((unused));
static inline double _powerL(double x, double lambda, int yj0, double low, double high){
  if (!R_finite(x)) return NA_REAL;
  double x0 = x, hl, xl, hl2, eri;
  int yj, dist;
  _splitYj(&yj0, &dist,  &yj);
  switch(yj){
  case 7:
    return log(_powerDD(_powerD(x, lambda, 6, low, high), lambda, 1, low, high))+log(_powerDD(x, lambda, 6, low, high));
  case 6:
    hl = (high-low);
    eri = erfinv(-1+2*(-low+x)/hl);
    return 0.918938533204672669541 +eri*eri-log(hl);
  case 5:
    // Subs(Derivative(yeoJohnson(_xi_1), _xi_1), (_xi_1), (logit(x)))*Derivative(logit(x), x)
    return log(_powerDD(_powerD(x, lambda, 4, low, high), lambda, 1, low, high))+log(_powerDD(x, lambda, 4, low, high));
  case 4: // logit d/dx(logit(x))
    xl = (x-low);
    if (xl <= _eps) xl = _eps;
    hl = (high - low);
    hl2 = hl-xl;
    if (hl2 <= _eps) hl2 = _eps;
    return log(hl)-log(xl)-log(hl2);
    /* return 0; */
  case 3:
    if (x <= _eps) x0 = _eps;
    return -log(x0);
  case 2:
    return 0;
  case 0:
    if (lambda == 1.0) return 0;
    if (x <= _eps) x0 = _eps;
    return (lambda - 1.0)*log(x0);
  case 1:
    if (x >= 0) return (lambda - 1.0)*log1p(x);
    return (1.0-lambda)*log1p(-x);
  }
  return NA_REAL;
  // d = 0.0 for cox box
  // d = 1.0 fo  Yeo- Johnson
  // logLik approximation
  // y^(lambda)/lambda - 1/lambda
  // dh/dy = y^(lambda-1)
  // log(dh/dy) = (lambda-1)*log(y) + log(lambda)
  //
  // (x + 1.0)^(lambda)/lambda - 1/lambda
  // dh/dy = (x+1.0)^(lambda-1)
  // log(dh/dy) = (lambda-1)*log(x+1.0)

  // For negative values yj becomes
  // (-x+1)^(2-lambda)/(2-lambda) - 1/(2-lambda)
  // dh/dy = (-x+1)^(1-lambda)
  // log(dh/dy) = (1-lambda)*log(-x+1)
}

// extra liklihood
static inline double _powerDL(double x, double lambda, int yj0, double low, double hi) __attribute__((unused));
static inline double _powerDL(double x, double lambda, int yj0, double low, double hi) {
  // d(logLik/dlambda)
  if (!R_finite(x)) return NA_REAL;
  double x0 = x;
  int yj, dist;
  _splitYj(&yj0, &dist,  &yj);
  switch (yj){
  case 7:
    return _powerDL(_powerD(x, lambda, 6, low, hi), lambda, 1, low, hi);
  case 6:
    return 0; // does not depend on lambda
  case 5:
    return _powerDL(_powerD(x, lambda, 4, low, hi), lambda, 1, low, hi);
  case 4:
    // For logit norm, no dependence on lambda
    return 0;
  case 3:
    // log-Jacobian of logNorm is -log(x): no dependence on lambda
    return 0;
  case 2:
    // For normal transform no dependence of lambda
    return 0;
  case 0:
    // d/dlambda[(lambda-1)*log(x)] = log(x) for every lambda
    if (x <= _eps) x0 = _eps;
    return log(x0);
  case 1:
    // d/dlambda[(lambda-1)*log1p(x)] and d/dlambda[(1-lambda)*log1p(-x)]
    if (x >= 0) return log1p(x);
    return -log1p(-x);
  }
  return NA_REAL;
  // d = 0.0 for cox box
  // d = 1.0 fo  Yeo- Johnson
  // logLik approximation
  // y^(lambda)/lambda - 1/lambda
  // dh/dy = y^(lambda-1)
  // log(dh/dy) = (lambda-1)*log(y) + log(lambda)
  //
  // (x + 1.0)^(lambda)/lambda - 1/lambda
  // dh/dy = (x+1.0)^(lambda-1)
  // log(dh/dy) = (lambda-1)*log(x+1.0)

  // For negative values yj becomes
  // (-x+1)^(2-lambda)/(2-lambda) - 1/(2-lambda)
  // dh/dy = (-x+1)^(1-lambda)
  // log(dh/dy) = (1-lambda)*log(-x+1)
}

static inline double abs1(double x){
  if (x == 0.0) return 1.0;
  return fabs(x);
}
static inline double dabs1(double x){
  return (x>0)-(x<0);
}

static inline double dabs(double x) {
  return (x>0)-(x<0);
}

static inline double dabs2(double x) {
  return 0.0;
}

#if defined(__cplusplus)
  extern "C" rx_solve rx_global;
  extern "C" rx_solving_options op_global;
  extern "C" rx_solving_options_ind *inds_global;
  extern "C" rx_solving_options_ind *inds_thread;
  // Event ("jump") sensitivities: model dF/dLag/dydt functions + runtime dims
  // (defined in par_solve.cpp); used by handle_evid to inject dosing-parameter
  // jumps.  dydt provides the Jacobian column by central difference.
  extern "C" t_dF dF;
  extern "C" t_dLag dLagEs;
  extern "C" t_dRate dRateEs;
  extern "C" t_dDur dDurEs;
  extern "C" t_dF d2FEs;
  extern "C" t_dLag d2LagEs;
  extern "C" t_dRate d2RateEs;
  extern "C" t_dDur d2DurEs;
  extern "C" t_dF d3FEs;
  extern "C" t_dF dFQEs;
  extern "C" t_dLag dLagJacEs;
  extern "C" t_dLag dLagQEs;
  extern "C" t_dDur dDurQEs;
  extern "C" t_DUR durEsFn;
  extern "C" t_dydt dydtEs;
  // matExp()/indLin() models have no functional dydt() (the primal system is
  // solved by matrix-exponential propagation, not RHS evaluation), so the
  // dtau/lag jump row's usual central-difference-of-dydt Jacobian column is
  // always zero for them. `_rxEsUseCalcJac` (set alongside the dims, from the
  // model's `mv$indLin`) tells handle_evid to read the Jacobian column from
  // `calc_jac` instead for such models (populated via explicit df/dy lines
  // rxSensMatExp() emits from its already-known Jacobian).
  extern "C" t_calc_jac calc_jac;
  extern "C" int _rxEsUseCalcJac;
  extern "C" int _rxEsActive;
  extern "C" int _rxEsNState;
  extern "C" int _rxEsNParam;
  extern "C" int _rxEsNParam2;
  extern "C" int _rxEsNParam3;
  extern "C" int _esSSDurOffCmt;
  extern "C" int _esSSRateOffCmt;
#else
  extern rx_solve rx_global;
  extern rx_solving_options op_global;
  extern rx_solving_options_ind *inds_global;
  extern rx_solving_options_ind *inds_thread;
  extern t_dF dF;
  extern t_dLag dLagEs;
  extern t_dRate dRateEs;
  extern t_dDur dDurEs;
  extern t_dF d2FEs;
  extern t_dLag d2LagEs;
  extern t_dRate d2RateEs;
  extern t_dDur d2DurEs;
  extern t_dF d3FEs;
  extern t_dF dFQEs;
  extern t_dLag dLagJacEs;
  extern t_dLag dLagQEs;
  extern t_dDur dDurQEs;
  extern t_DUR durEsFn;
  extern t_dydt dydtEs;
  extern t_calc_jac calc_jac;
  extern int _rxEsUseCalcJac;
  extern int _rxEsActive;
  extern int _rxEsNState;
  extern int _rxEsNParam;
  extern int _rxEsNParam2;
  extern int _rxEsNParam3;
  extern int _esSSDurOffCmt;
  extern int _esSSRateOffCmt;
#endif



#endif
#if defined(__cplusplus)
}
#endif
