text: // Automatically generated by odin 1.1.15 - do not edit
  #include <R.h>
  #include <Rmath.h>
  #include <Rinternals.h>
  #include <stdbool.h>
  #include <R_ext/Rdynload.h>
  #ifndef CINTERPOLTE_CINTERPOLATE_H_
  #define CINTERPOLTE_CINTERPOLATE_H_
  
  // Allow use from C++
  #ifdef __cplusplus
  extern "C" {
  #endif
  
  // There are only three functions in the interface; allocation,
  // evaluation and freeing.
  
  // Allocate an interpolation object.
  //
  //   type: The mode of interpolation. Must be one of "constant",
  //       "linear" or "spline" (an R error is thrown if a different
  //       value is given).
  //
  //   n: The number of `x` points to interpolate over
  //
  //   ny: the number of `y` points per `x` point.  This is 1 in the
  //       case of zimple interpolation as used by Rs `interpolate()`
  //
  //   x: an array of `x` values of length `n`
  //
  //   y: an array of `ny` sets of `y` values.  This is in R's matrix
  //       order (i.e., the first `n` values are the first series to
  //       interpolate over).
  //
  //   fail_on_extrapolate: if true, when an extrapolation occurs throw
  //       an error; if false return NA_REAL
  //
  //   auto_free: automatically clean up the interpolation object on
  //       return to R. This uses `R_alloc` for allocations rather than
  //       `Calloc` so freeing will always happen (even on error
  //       elsewhere in the code). However, this prevents returning back
  //       a pointer to R that will last longer than the call into C
  //       code.
  //
  // The return value is an opaque pointer that can be passed through to
  // `cinterpolate_eval` and `cinterpolate_free`
  void *cinterpolate_alloc(const char *type, size_t n, size_t ny,
                           double *x, double *y, bool fail_on_extrapolate,
                           bool auto_free);
  
  // Evaluate the interpolated function at a new `x` point.
  //
  //   x: A new, single, `x` point to interpolate `y` values to
  //
  //   obj: The interpolation object, as returned by `cinterpolate_alloc`
  //
  //   y: An array of length `ny` to store the interpolated values
  //
  // The return value is 0 if the interpolation is successful (with x
  // lying within the range of values that the interpolation function
  // supports), -1 otherwise
  int cinterpolate_eval(double x, void *obj, double *y);
  
  // Clean up all allocated memory
  //
  //   obj: The interpolation object, as returned by `cinterpolate_alloc`
  void cinterpolate_free(void *obj);
  
  #ifdef __cplusplus
  }
  #endif
  
  #endif
  typedef struct agesirdGenerator_internal {
    double Alpha;
    double *Beta;
    double *Betas;
    double *changeTimes;
    double *D0;
    int dim_Beta;
    int dim_Beta_1;
    int dim_Beta_2;
    int dim_Betas;
    int dim_Betas_1;
    int dim_Betas_12;
    int dim_Betas_2;
    int dim_Betas_3;
    int dim_changeTimes;
    int dim_D;
    int dim_D0;
    int dim_I;
    int dim_I0;
    int dim_infectiousWeight;
    int dim_infectiousWeight_1;
    int dim_infectiousWeight_2;
    int dim_N;
    int dim_R;
    int dim_R0;
    int dim_S;
    int dim_S0;
    double Gamma;
    double *I0;
    double *infectiousWeight;
    double *initial_D;
    double *initial_I;
    double *initial_R;
    double *initial_S;
    void *interpolate_Beta;
    double *N;
    int offset_variable_D;
    int offset_variable_R;
    double *R0;
    double *S0;
  } agesirdGenerator_internal;
  typedef struct seirdGenerator_internal {
    double Alpha;
    double *Betas;
    double *changeTimes;
    double D0;
    int dim_Betas;
    int dim_changeTimes;
    double E0;
    double Gamma;
    double I0;
    double initial_D;
    double initial_E;
    double initial_I;
    double initial_R;
    double initial_S;
    void *interpolate_Beta;
    double Lambda;
    double R0;
    double S0;
  } seirdGenerator_internal;
  typedef struct sirdGenerator_internal {
    double Alpha;
    double *Betas;
    double *changeTimes;
    double D0;
    int dim_Betas;
    int dim_changeTimes;
    double Gamma;
    double I0;
    double initial_D;
    double initial_I;
    double initial_R;
    double initial_S;
    void *interpolate_Beta;
    double R0;
    double S0;
  } sirdGenerator_internal;
  agesirdGenerator_internal* agesirdGenerator_get_internal(SEXP internal_p, int closed_error);
  static void agesirdGenerator_finalise(SEXP internal_p);
  SEXP agesirdGenerator_create(SEXP user);
  void agesirdGenerator_initmod_desolve(void(* odeparms) (int *, double *));
  SEXP agesirdGenerator_contents(SEXP internal_p);
  SEXP agesirdGenerator_set_user(SEXP internal_p, SEXP user);
  SEXP agesirdGenerator_set_initial(SEXP internal_p, SEXP t_ptr, SEXP state_ptr, SEXP agesirdGenerator_use_dde_ptr);
  SEXP agesirdGenerator_metadata(SEXP internal_p);
  SEXP agesirdGenerator_initial_conditions(SEXP internal_p, SEXP t_ptr);
  void agesirdGenerator_rhs(agesirdGenerator_internal* internal, double t, double * state, double * dstatedt, double * output);
  void agesirdGenerator_rhs_dde(size_t neq, double t, double * state, double * dstatedt, void * internal);
  void agesirdGenerator_rhs_desolve(int * neq, double * t, double * state, double * dstatedt, double * output, int * np);
  void agesirdGenerator_output_dde(size_t n_eq, double t, double * state, size_t n_output, double * output, void * internal_p);
  SEXP agesirdGenerator_rhs_r(SEXP internal_p, SEXP t, SEXP state);
  seirdGenerator_internal* seirdGenerator_get_internal(SEXP internal_p, int closed_error);
  static void seirdGenerator_finalise(SEXP internal_p);
  SEXP seirdGenerator_create(SEXP user);
  void seirdGenerator_initmod_desolve(void(* odeparms) (int *, double *));
  SEXP seirdGenerator_contents(SEXP internal_p);
  SEXP seirdGenerator_set_user(SEXP internal_p, SEXP user);
  SEXP seirdGenerator_set_initial(SEXP internal_p, SEXP t_ptr, SEXP state_ptr, SEXP seirdGenerator_use_dde_ptr);
  SEXP seirdGenerator_metadata(SEXP internal_p);
  SEXP seirdGenerator_initial_conditions(SEXP internal_p, SEXP t_ptr);
  void seirdGenerator_rhs(seirdGenerator_internal* internal, double t, double * state, double * dstatedt, double * output);
  void seirdGenerator_rhs_dde(size_t neq, double t, double * state, double * dstatedt, void * internal);
  void seirdGenerator_rhs_desolve(int * neq, double * t, double * state, double * dstatedt, double * output, int * np);
  void seirdGenerator_output_dde(size_t n_eq, double t, double * state, size_t n_output, double * output, void * internal_p);
  SEXP seirdGenerator_rhs_r(SEXP internal_p, SEXP t, SEXP state);
  sirdGenerator_internal* sirdGenerator_get_internal(SEXP internal_p, int closed_error);
  static void sirdGenerator_finalise(SEXP internal_p);
  SEXP sirdGenerator_create(SEXP user);
  void sirdGenerator_initmod_desolve(void(* odeparms) (int *, double *));
  SEXP sirdGenerator_contents(SEXP internal_p);
  SEXP sirdGenerator_set_user(SEXP internal_p, SEXP user);
  SEXP sirdGenerator_set_initial(SEXP internal_p, SEXP t_ptr, SEXP state_ptr, SEXP sirdGenerator_use_dde_ptr);
  SEXP sirdGenerator_metadata(SEXP internal_p);
  SEXP sirdGenerator_initial_conditions(SEXP internal_p, SEXP t_ptr);
  void sirdGenerator_rhs(sirdGenerator_internal* internal, double t, double * state, double * dstatedt, double * output);
  void sirdGenerator_rhs_dde(size_t neq, double t, double * state, double * dstatedt, void * internal);
  void sirdGenerator_rhs_desolve(int * neq, double * t, double * state, double * dstatedt, double * output, int * np);
  void sirdGenerator_output_dde(size_t n_eq, double t, double * state, size_t n_output, double * output, void * internal_p);
  SEXP sirdGenerator_rhs_r(SEXP internal_p, SEXP t, SEXP state);
  double user_get_scalar_double(SEXP user, const char *name,
                                double default_value, double min, double max);
  int user_get_scalar_int(SEXP user, const char *name,
                          int default_value, double min, double max);
  void user_check_values_double(double * value, size_t len,
                                    double min, double max, const char *name);
  void user_check_values_int(int * value, size_t len,
                                 double min, double max, const char *name);
  void user_check_values(SEXP value, double min, double max,
                             const char *name);
  SEXP user_list_element(SEXP list, const char *name);
  void odin_set_dim(SEXP target, int rank, ...);
  void* user_get_array_dim(SEXP user, bool is_integer, void * previous,
                           const char *name, int rank,
                           double min, double max, int *dest_dim);
  void* user_get_array(SEXP user, bool is_integer, void * previous,
                       const char *name, double min, double max,
                       int rank, ...);
  SEXP user_get_array_check(SEXP el, bool is_integer, const char *name,
                            double min, double max);
  SEXP user_get_array_check_rank(SEXP user, const char *name, int rank,
                                 bool required);
  void interpolate_check_y(size_t nx, size_t ny, size_t i, const char *name_arg, const char *name_target);
  double scalar_real(SEXP x, const char * name);
  double odin_sum1(double *x, size_t from, size_t to);
  double odin_sum2(double* x, int from_i, int to_i, int from_j, int to_j, int dim_x_1);
  double odin_sum3(double* x, int from_i, int to_i, int from_j, int to_j, int from_k, int to_k, int dim_x_1, int dim_x_12);
  agesirdGenerator_internal* agesirdGenerator_get_internal(SEXP internal_p, int closed_error) {
    agesirdGenerator_internal *internal = NULL;
    if (TYPEOF(internal_p) != EXTPTRSXP) {
      Rf_error("Expected an external pointer");
    }
    internal = (agesirdGenerator_internal*) R_ExternalPtrAddr(internal_p);
    if (!internal && closed_error) {
      Rf_error("Pointer has been invalidated");
    }
    return internal;
  }
  void agesirdGenerator_finalise(SEXP internal_p) {
    agesirdGenerator_internal *internal = agesirdGenerator_get_internal(internal_p, 0);
    if (internal_p) {
      cinterpolate_free(internal->interpolate_Beta);
      internal->interpolate_Beta = NULL;
      Free(internal->Beta);
      Free(internal->Betas);
      Free(internal->changeTimes);
      Free(internal->D0);
      Free(internal->I0);
      Free(internal->infectiousWeight);
      Free(internal->initial_D);
      Free(internal->initial_I);
      Free(internal->initial_R);
      Free(internal->initial_S);
      Free(internal->N);
      Free(internal->R0);
      Free(internal->S0);
      Free(internal);
      R_ClearExternalPtr(internal_p);
    }
  }
  SEXP agesirdGenerator_create(SEXP user) {
    agesirdGenerator_internal *internal = (agesirdGenerator_internal*) Calloc(1, agesirdGenerator_internal);
    internal->Beta = NULL;
    internal->Betas = NULL;
    internal->changeTimes = NULL;
    internal->D0 = NULL;
    internal->I0 = NULL;
    internal->infectiousWeight = NULL;
    internal->initial_D = NULL;
    internal->initial_I = NULL;
    internal->initial_R = NULL;
    internal->initial_S = NULL;
    internal->interpolate_Beta = NULL;
    internal->N = NULL;
    internal->R0 = NULL;
    internal->S0 = NULL;
    internal->Alpha = NA_REAL;
    internal->Betas = NULL;
    internal->changeTimes = NULL;
    internal->D0 = NULL;
    internal->Gamma = NA_REAL;
    internal->I0 = NULL;
    internal->R0 = NULL;
    internal->S0 = NULL;
    SEXP ptr = PROTECT(R_MakeExternalPtr(internal, R_NilValue, R_NilValue));
    R_RegisterCFinalizer(ptr, agesirdGenerator_finalise);
    UNPROTECT(1);
    return ptr;
  }
  static agesirdGenerator_internal *agesirdGenerator_internal_ds;
  void agesirdGenerator_initmod_desolve(void(* odeparms) (int *, double *)) {
    static DL_FUNC get_desolve_gparms = NULL;
    if (get_desolve_gparms == NULL) {
      get_desolve_gparms =
        R_GetCCallable("deSolve", "get_deSolve_gparms");
    }
    agesirdGenerator_internal_ds = agesirdGenerator_get_internal(get_desolve_gparms(), 1);
  }
  SEXP agesirdGenerator_contents(SEXP internal_p) {
    agesirdGenerator_internal *internal = agesirdGenerator_get_internal(internal_p, 1);
    SEXP contents = PROTECT(allocVector(VECSXP, 39));
    SET_VECTOR_ELT(contents, 0, ScalarReal(internal->Alpha));
    SEXP Beta = PROTECT(allocVector(REALSXP, internal->dim_Beta));
    memcpy(REAL(Beta), internal->Beta, internal->dim_Beta * sizeof(double));
    odin_set_dim(Beta, 2, internal->dim_Beta_1, internal->dim_Beta_2);
    SET_VECTOR_ELT(contents, 1, Beta);
    SEXP Betas = PROTECT(allocVector(REALSXP, internal->dim_Betas));
    memcpy(REAL(Betas), internal->Betas, internal->dim_Betas * sizeof(double));
    odin_set_dim(Betas, 3, internal->dim_Betas_1, internal->dim_Betas_2, internal->dim_Betas_3);
    SET_VECTOR_ELT(contents, 2, Betas);
    SEXP changeTimes = PROTECT(allocVector(REALSXP, internal->dim_changeTimes));
    memcpy(REAL(changeTimes), internal->changeTimes, internal->dim_changeTimes * sizeof(double));
    SET_VECTOR_ELT(contents, 3, changeTimes);
    SEXP D0 = PROTECT(allocVector(REALSXP, internal->dim_D0));
    memcpy(REAL(D0), internal->D0, internal->dim_D0 * sizeof(double));
    SET_VECTOR_ELT(contents, 4, D0);
    SET_VECTOR_ELT(contents, 5, ScalarInteger(internal->dim_Beta));
    SET_VECTOR_ELT(contents, 6, ScalarInteger(internal->dim_Beta_1));
    SET_VECTOR_ELT(contents, 7, ScalarInteger(internal->dim_Beta_2));
    SET_VECTOR_ELT(contents, 8, ScalarInteger(internal->dim_Betas));
    SET_VECTOR_ELT(contents, 9, ScalarInteger(internal->dim_Betas_1));
    SET_VECTOR_ELT(contents, 10, ScalarInteger(internal->dim_Betas_12));
    SET_VECTOR_ELT(contents, 11, ScalarInteger(internal->dim_Betas_2));
    SET_VECTOR_ELT(contents, 12, ScalarInteger(internal->dim_Betas_3));
    SET_VECTOR_ELT(contents, 13, ScalarInteger(internal->dim_changeTimes));
    SET_VECTOR_ELT(contents, 14, ScalarInteger(internal->dim_D));
    SET_VECTOR_ELT(contents, 15, ScalarInteger(internal->dim_D0));
    SET_VECTOR_ELT(contents, 16, ScalarInteger(internal->dim_I));
    SET_VECTOR_ELT(contents, 17, ScalarInteger(internal->dim_I0));
    SET_VECTOR_ELT(contents, 18, ScalarInteger(internal->dim_infectiousWeight));
    SET_VECTOR_ELT(contents, 19, ScalarInteger(internal->dim_infectiousWeight_1));
    SET_VECTOR_ELT(contents, 20, ScalarInteger(internal->dim_infectiousWeight_2));
    SET_VECTOR_ELT(contents, 21, ScalarInteger(internal->dim_N));
    SET_VECTOR_ELT(contents, 22, ScalarInteger(internal->dim_R));
    SET_VECTOR_ELT(contents, 23, ScalarInteger(internal->dim_R0));
    SET_VECTOR_ELT(contents, 24, ScalarInteger(internal->dim_S));
    SET_VECTOR_ELT(contents, 25, ScalarInteger(internal->dim_S0));
    SET_VECTOR_ELT(contents, 26, ScalarReal(internal->Gamma));
    SEXP I0 = PROTECT(allocVector(REALSXP, internal->dim_I0));
    memcpy(REAL(I0), internal->I0, internal->dim_I0 * sizeof(double));
    SET_VECTOR_ELT(contents, 27, I0);
    SEXP infectiousWeight = PROTECT(allocVector(REALSXP, internal->dim_infectiousWeight));
    memcpy(REAL(infectiousWeight), internal->infectiousWeight, internal->dim_infectiousWeight * sizeof(double));
    odin_set_dim(infectiousWeight, 2, internal->dim_infectiousWeight_1, internal->dim_infectiousWeight_2);
    SET_VECTOR_ELT(contents, 28, infectiousWeight);
    SEXP initial_D = PROTECT(allocVector(REALSXP, internal->dim_D));
    memcpy(REAL(initial_D), internal->initial_D, internal->dim_D * sizeof(double));
    SET_VECTOR_ELT(contents, 29, initial_D);
    SEXP initial_I = PROTECT(allocVector(REALSXP, internal->dim_I));
    memcpy(REAL(initial_I), internal->initial_I, internal->dim_I * sizeof(double));
    SET_VECTOR_ELT(contents, 30, initial_I);
    SEXP initial_R = PROTECT(allocVector(REALSXP, internal->dim_R));
    memcpy(REAL(initial_R), internal->initial_R, internal->dim_R * sizeof(double));
    SET_VECTOR_ELT(contents, 31, initial_R);
    SEXP initial_S = PROTECT(allocVector(REALSXP, internal->dim_S));
    memcpy(REAL(initial_S), internal->initial_S, internal->dim_S * sizeof(double));
    SET_VECTOR_ELT(contents, 32, initial_S);
    SEXP N = PROTECT(allocVector(REALSXP, internal->dim_N));
    memcpy(REAL(N), internal->N, internal->dim_N * sizeof(double));
    SET_VECTOR_ELT(contents, 34, N);
    SET_VECTOR_ELT(contents, 35, ScalarInteger(internal->offset_variable_D));
    SET_VECTOR_ELT(contents, 36, ScalarInteger(internal->offset_variable_R));
    SEXP R0 = PROTECT(allocVector(REALSXP, internal->dim_R0));
    memcpy(REAL(R0), internal->R0, internal->dim_R0 * sizeof(double));
    SET_VECTOR_ELT(contents, 37, R0);
    SEXP S0 = PROTECT(allocVector(REALSXP, internal->dim_S0));
    memcpy(REAL(S0), internal->S0, internal->dim_S0 * sizeof(double));
    SET_VECTOR_ELT(contents, 38, S0);
    SEXP nms = PROTECT(allocVector(STRSXP, 39));
    SET_STRING_ELT(nms, 0, mkChar("Alpha"));
    SET_STRING_ELT(nms, 1, mkChar("Beta"));
    SET_STRING_ELT(nms, 2, mkChar("Betas"));
    SET_STRING_ELT(nms, 3, mkChar("changeTimes"));
    SET_STRING_ELT(nms, 4, mkChar("D0"));
    SET_STRING_ELT(nms, 5, mkChar("dim_Beta"));
    SET_STRING_ELT(nms, 6, mkChar("dim_Beta_1"));
    SET_STRING_ELT(nms, 7, mkChar("dim_Beta_2"));
    SET_STRING_ELT(nms, 8, mkChar("dim_Betas"));
    SET_STRING_ELT(nms, 9, mkChar("dim_Betas_1"));
    SET_STRING_ELT(nms, 10, mkChar("dim_Betas_12"));
    SET_STRING_ELT(nms, 11, mkChar("dim_Betas_2"));
    SET_STRING_ELT(nms, 12, mkChar("dim_Betas_3"));
    SET_STRING_ELT(nms, 13, mkChar("dim_changeTimes"));
    SET_STRING_ELT(nms, 14, mkChar("dim_D"));
    SET_STRING_ELT(nms, 15, mkChar("dim_D0"));
    SET_STRING_ELT(nms, 16, mkChar("dim_I"));
    SET_STRING_ELT(nms, 17, mkChar("dim_I0"));
    SET_STRING_ELT(nms, 18, mkChar("dim_infectiousWeight"));
    SET_STRING_ELT(nms, 19, mkChar("dim_infectiousWeight_1"));
    SET_STRING_ELT(nms, 20, mkChar("dim_infectiousWeight_2"));
    SET_STRING_ELT(nms, 21, mkChar("dim_N"));
    SET_STRING_ELT(nms, 22, mkChar("dim_R"));
    SET_STRING_ELT(nms, 23, mkChar("dim_R0"));
    SET_STRING_ELT(nms, 24, mkChar("dim_S"));
    SET_STRING_ELT(nms, 25, mkChar("dim_S0"));
    SET_STRING_ELT(nms, 26, mkChar("Gamma"));
    SET_STRING_ELT(nms, 27, mkChar("I0"));
    SET_STRING_ELT(nms, 28, mkChar("infectiousWeight"));
    SET_STRING_ELT(nms, 29, mkChar("initial_D"));
    SET_STRING_ELT(nms, 30, mkChar("initial_I"));
    SET_STRING_ELT(nms, 31, mkChar("initial_R"));
    SET_STRING_ELT(nms, 32, mkChar("initial_S"));
    SET_STRING_ELT(nms, 33, mkChar("interpolate_Beta"));
    SET_STRING_ELT(nms, 34, mkChar("N"));
    SET_STRING_ELT(nms, 35, mkChar("offset_variable_D"));
    SET_STRING_ELT(nms, 36, mkChar("offset_variable_R"));
    SET_STRING_ELT(nms, 37, mkChar("R0"));
    SET_STRING_ELT(nms, 38, mkChar("S0"));
    setAttrib(contents, R_NamesSymbol, nms);
    UNPROTECT(15);
    return contents;
  }
  SEXP agesirdGenerator_set_user(SEXP internal_p, SEXP user) {
    agesirdGenerator_internal *internal = agesirdGenerator_get_internal(internal_p, 1);
    internal->Alpha = user_get_scalar_double(user, "Alpha", internal->Alpha, NA_REAL, NA_REAL);
    int dim_Betas[4];
    internal->Betas = (double*) user_get_array_dim(user, false, internal->Betas, "Betas", 3, NA_REAL, NA_REAL, dim_Betas);
    internal->dim_Betas = dim_Betas[0];
    internal->dim_Betas_1 = dim_Betas[1];
    internal->dim_Betas_2 = dim_Betas[2];
    internal->dim_Betas_3 = dim_Betas[3];
    internal->changeTimes = (double*) user_get_array_dim(user, false, internal->changeTimes, "changeTimes", 1, NA_REAL, NA_REAL, &internal->dim_changeTimes);
    internal->Gamma = user_get_scalar_double(user, "Gamma", internal->Gamma, NA_REAL, NA_REAL);
    internal->S0 = (double*) user_get_array_dim(user, false, internal->S0, "S0", 1, NA_REAL, NA_REAL, &internal->dim_S0);
    internal->dim_Beta_1 = internal->dim_S0;
    internal->dim_Beta_2 = internal->dim_S0;
    internal->dim_Betas_12 = internal->dim_Betas_1 * internal->dim_Betas_2;
    internal->dim_D = internal->dim_S0;
    internal->dim_D0 = internal->dim_S0;
    internal->dim_I = internal->dim_S0;
    internal->dim_I0 = internal->dim_S0;
    internal->dim_infectiousWeight_1 = internal->dim_S0;
    internal->dim_infectiousWeight_2 = internal->dim_S0;
    internal->dim_N = internal->dim_S0;
    internal->dim_R = internal->dim_S0;
    internal->dim_R0 = internal->dim_S0;
    internal->dim_S = internal->dim_S0;
    Free(internal->initial_D);
    internal->initial_D = (double*) Calloc(internal->dim_D, double);
    Free(internal->initial_I);
    internal->initial_I = (double*) Calloc(internal->dim_I, double);
    Free(internal->initial_R);
    internal->initial_R = (double*) Calloc(internal->dim_R, double);
    Free(internal->initial_S);
    internal->initial_S = (double*) Calloc(internal->dim_S, double);
    Free(internal->N);
    internal->N = (double*) Calloc(internal->dim_N, double);
    internal->D0 = (double*) user_get_array(user, false, internal->D0, "D0", NA_REAL, NA_REAL, 1, internal->dim_D0);
    internal->dim_Beta = internal->dim_Beta_1 * internal->dim_Beta_2;
    internal->dim_infectiousWeight = internal->dim_infectiousWeight_1 * internal->dim_infectiousWeight_2;
    internal->I0 = (double*) user_get_array(user, false, internal->I0, "I0", NA_REAL, NA_REAL, 1, internal->dim_I0);
    for (int i = 1; i <= internal->dim_S; ++i) {
      internal->initial_S[i - 1] = internal->S0[i - 1];
    }
    internal->offset_variable_D = internal->dim_I + internal->dim_R + internal->dim_S;
    internal->offset_variable_R = internal->dim_I + internal->dim_S;
    internal->R0 = (double*) user_get_array(user, false, internal->R0, "R0", NA_REAL, NA_REAL, 1, internal->dim_R0);
    Free(internal->Beta);
    internal->Beta = (double*) Calloc(internal->dim_Beta, double);
    Free(internal->infectiousWeight);
    internal->infectiousWeight = (double*) Calloc(internal->dim_infectiousWeight, double);
    for (int i = 1; i <= internal->dim_D; ++i) {
      internal->initial_D[i - 1] = internal->D0[i - 1];
    }
    for (int i = 1; i <= internal->dim_I; ++i) {
      internal->initial_I[i - 1] = internal->I0[i - 1];
    }
    for (int i = 1; i <= internal->dim_R; ++i) {
      internal->initial_R[i - 1] = internal->R0[i - 1];
    }
    interpolate_check_y(internal->dim_changeTimes, internal->dim_Betas_1, 1, "Betas", "Beta");
    interpolate_check_y(internal->dim_Beta_1, internal->dim_Betas_2, 2, "Betas", "Beta");
    interpolate_check_y(internal->dim_Beta_2, internal->dim_Betas_3, 3, "Betas", "Beta");
    cinterpolate_free(internal->interpolate_Beta);
    internal->interpolate_Beta = cinterpolate_alloc("constant", internal->dim_changeTimes, internal->dim_Beta, internal->changeTimes, internal->Betas, true, false);
    return R_NilValue;
  }
  SEXP agesirdGenerator_set_initial(SEXP internal_p, SEXP t_ptr, SEXP state_ptr, SEXP agesirdGenerator_use_dde_ptr) {
    return R_NilValue;
  }
  SEXP agesirdGenerator_metadata(SEXP internal_p) {
    agesirdGenerator_internal *internal = agesirdGenerator_get_internal(internal_p, 1);
    SEXP ret = PROTECT(allocVector(VECSXP, 4));
    SEXP nms = PROTECT(allocVector(STRSXP, 4));
    SET_STRING_ELT(nms, 0, mkChar("variable_order"));
    SET_STRING_ELT(nms, 1, mkChar("output_order"));
    SET_STRING_ELT(nms, 2, mkChar("n_out"));
    SET_STRING_ELT(nms, 3, mkChar("interpolate_t"));
    setAttrib(ret, R_NamesSymbol, nms);
    SEXP variable_length = PROTECT(allocVector(VECSXP, 4));
    SEXP variable_names = PROTECT(allocVector(STRSXP, 4));
    setAttrib(variable_length, R_NamesSymbol, variable_names);
    SET_VECTOR_ELT(variable_length, 0, ScalarInteger(internal->dim_S));
    SET_VECTOR_ELT(variable_length, 1, ScalarInteger(internal->dim_I));
    SET_VECTOR_ELT(variable_length, 2, ScalarInteger(internal->dim_R));
    SET_VECTOR_ELT(variable_length, 3, ScalarInteger(internal->dim_D));
    SET_STRING_ELT(variable_names, 0, mkChar("S"));
    SET_STRING_ELT(variable_names, 1, mkChar("I"));
    SET_STRING_ELT(variable_names, 2, mkChar("R"));
    SET_STRING_ELT(variable_names, 3, mkChar("D"));
    SET_VECTOR_ELT(ret, 0, variable_length);
    UNPROTECT(2);
    SEXP output_length = PROTECT(allocVector(VECSXP, 1));
    SEXP output_names = PROTECT(allocVector(STRSXP, 1));
    setAttrib(output_length, R_NamesSymbol, output_names);
    SET_VECTOR_ELT(output_length, 0, allocVector(INTSXP, 2));
    int * dim_Beta = INTEGER(VECTOR_ELT(output_length, 0));
    dim_Beta[0] = internal->dim_Beta_1;
    dim_Beta[1] = internal->dim_Beta_2;
    SET_STRING_ELT(output_names, 0, mkChar("Beta"));
    SET_VECTOR_ELT(ret, 1, output_length);
    UNPROTECT(2);
    SET_VECTOR_ELT(ret, 2, ScalarInteger(internal->dim_Beta));
    SEXP interpolate_t = PROTECT(allocVector(VECSXP, 3));
    SEXP interpolate_t_nms = PROTECT(allocVector(STRSXP, 3));
    setAttrib(interpolate_t, R_NamesSymbol, interpolate_t_nms);
    SET_VECTOR_ELT(interpolate_t, 0, ScalarReal(internal->changeTimes[0]));
    SET_VECTOR_ELT(interpolate_t, 1, ScalarReal(R_PosInf));
    SET_STRING_ELT(interpolate_t_nms, 0, mkChar("min"));
    SET_STRING_ELT(interpolate_t_nms, 1, mkChar("max"));
    SET_VECTOR_ELT(ret, 3, interpolate_t);
    UNPROTECT(2);
    UNPROTECT(2);
    return ret;
  }
  SEXP agesirdGenerator_initial_conditions(SEXP internal_p, SEXP t_ptr) {
    agesirdGenerator_internal *internal = agesirdGenerator_get_internal(internal_p, 1);
    SEXP r_state = PROTECT(allocVector(REALSXP, internal->dim_D + internal->dim_I + internal->dim_R + internal->dim_S));
    double * state = REAL(r_state);
    memcpy(state + 0, internal->initial_S, internal->dim_S * sizeof(double));
    memcpy(state + internal->dim_S, internal->initial_I, internal->dim_I * sizeof(double));
    memcpy(state + internal->offset_variable_R, internal->initial_R, internal->dim_R * sizeof(double));
    memcpy(state + internal->offset_variable_D, internal->initial_D, internal->dim_D * sizeof(double));
    UNPROTECT(1);
    return r_state;
  }
  void agesirdGenerator_rhs(agesirdGenerator_internal* internal, double t, double * state, double * dstatedt, double * output) {
    double * S = state + 0;
    double * I = state + internal->dim_S;
    double * R = state + internal->offset_variable_R;
    double * D = state + internal->offset_variable_D;
    for (int i = 1; i <= internal->dim_D; ++i) {
      dstatedt[internal->offset_variable_D + i - 1] = I[i - 1] * internal->Alpha;
    }
    for (int i = 1; i <= internal->dim_R; ++i) {
      dstatedt[internal->offset_variable_R + i - 1] = I[i - 1] * internal->Gamma;
    }
    for (int i = 1; i <= internal->dim_N; ++i) {
      internal->N[i - 1] = S[i - 1] + I[i - 1] + R[i - 1] + D[i - 1];
    }
    cinterpolate_eval(t, internal->interpolate_Beta, internal->Beta);
    for (int i = 1; i <= internal->dim_infectiousWeight_1; ++i) {
      for (int j = 1; j <= internal->dim_infectiousWeight_2; ++j) {
        internal->infectiousWeight[i - 1 + internal->dim_infectiousWeight_1 * (j - 1)] = internal->Beta[internal->dim_Beta_1 * (j - 1) + i - 1] * I[i - 1] / (double) internal->N[i - 1];
      }
    }
    for (int i = 1; i <= internal->dim_I; ++i) {
      dstatedt[internal->dim_S + i - 1] = S[i - 1] * odin_sum2(internal->infectiousWeight, 0, internal->dim_infectiousWeight_1, i - 1, i, internal->dim_infectiousWeight_1) - I[i - 1] * internal->Gamma - I[i - 1] * internal->Alpha;
    }
    for (int i = 1; i <= internal->dim_S; ++i) {
      dstatedt[0 + i - 1] = -(S[i - 1]) * odin_sum2(internal->infectiousWeight, 0, internal->dim_infectiousWeight_1, i - 1, i, internal->dim_infectiousWeight_1);
    }
    if (output) {
      memcpy(output + 0, internal->Beta, internal->dim_Beta * sizeof(double));
    }
  }
  void agesirdGenerator_rhs_dde(size_t neq, double t, double * state, double * dstatedt, void * internal) {
    agesirdGenerator_rhs((agesirdGenerator_internal*)internal, t, state, dstatedt, NULL);
  }
  void agesirdGenerator_rhs_desolve(int * neq, double * t, double * state, double * dstatedt, double * output, int * np) {
    agesirdGenerator_rhs(agesirdGenerator_internal_ds, *t, state, dstatedt, output);
  }
  void agesirdGenerator_output_dde(size_t n_eq, double t, double * state, size_t n_output, double * output, void * internal_p) {
    agesirdGenerator_internal *internal = (agesirdGenerator_internal*) internal_p;
    cinterpolate_eval(t, internal->interpolate_Beta, internal->Beta);
    memcpy(output + 0, internal->Beta, internal->dim_Beta * sizeof(double));
  }
  SEXP agesirdGenerator_rhs_r(SEXP internal_p, SEXP t, SEXP state) {
    SEXP dstatedt = PROTECT(allocVector(REALSXP, LENGTH(state)));
    agesirdGenerator_internal *internal = agesirdGenerator_get_internal(internal_p, 1);
    SEXP output_ptr = PROTECT(allocVector(REALSXP, internal->dim_Beta));
    setAttrib(dstatedt, install("output"), output_ptr);
    UNPROTECT(1);
    double *output = REAL(output_ptr);
    agesirdGenerator_rhs(internal, scalar_real(t, "t"), REAL(state), REAL(dstatedt), output);
    UNPROTECT(1);
    return dstatedt;
  }
  seirdGenerator_internal* seirdGenerator_get_internal(SEXP internal_p, int closed_error) {
    seirdGenerator_internal *internal = NULL;
    if (TYPEOF(internal_p) != EXTPTRSXP) {
      Rf_error("Expected an external pointer");
    }
    internal = (seirdGenerator_internal*) R_ExternalPtrAddr(internal_p);
    if (!internal && closed_error) {
      Rf_error("Pointer has been invalidated");
    }
    return internal;
  }
  void seirdGenerator_finalise(SEXP internal_p) {
    seirdGenerator_internal *internal = seirdGenerator_get_internal(internal_p, 0);
    if (internal_p) {
      cinterpolate_free(internal->interpolate_Beta);
      internal->interpolate_Beta = NULL;
      Free(internal->Betas);
      Free(internal->changeTimes);
      Free(internal);
      R_ClearExternalPtr(internal_p);
    }
  }
  SEXP seirdGenerator_create(SEXP user) {
    seirdGenerator_internal *internal = (seirdGenerator_internal*) Calloc(1, seirdGenerator_internal);
    internal->Betas = NULL;
    internal->changeTimes = NULL;
    internal->Alpha = NA_REAL;
    internal->Betas = NULL;
    internal->changeTimes = NULL;
    internal->Gamma = NA_REAL;
    internal->I0 = NA_REAL;
    internal->Lambda = NA_REAL;
    internal->S0 = NA_REAL;
    internal->D0 = 0;
    internal->E0 = 0;
    internal->R0 = 0;
    SEXP ptr = PROTECT(R_MakeExternalPtr(internal, R_NilValue, R_NilValue));
    R_RegisterCFinalizer(ptr, seirdGenerator_finalise);
    UNPROTECT(1);
    return ptr;
  }
  static seirdGenerator_internal *seirdGenerator_internal_ds;
  void seirdGenerator_initmod_desolve(void(* odeparms) (int *, double *)) {
    static DL_FUNC get_desolve_gparms = NULL;
    if (get_desolve_gparms == NULL) {
      get_desolve_gparms =
        R_GetCCallable("deSolve", "get_deSolve_gparms");
    }
    seirdGenerator_internal_ds = seirdGenerator_get_internal(get_desolve_gparms(), 1);
  }
  SEXP seirdGenerator_contents(SEXP internal_p) {
    seirdGenerator_internal *internal = seirdGenerator_get_internal(internal_p, 1);
    SEXP contents = PROTECT(allocVector(VECSXP, 18));
    SET_VECTOR_ELT(contents, 0, ScalarReal(internal->Alpha));
    SEXP Betas = PROTECT(allocVector(REALSXP, internal->dim_Betas));
    memcpy(REAL(Betas), internal->Betas, internal->dim_Betas * sizeof(double));
    SET_VECTOR_ELT(contents, 1, Betas);
    SEXP changeTimes = PROTECT(allocVector(REALSXP, internal->dim_changeTimes));
    memcpy(REAL(changeTimes), internal->changeTimes, internal->dim_changeTimes * sizeof(double));
    SET_VECTOR_ELT(contents, 2, changeTimes);
    SET_VECTOR_ELT(contents, 3, ScalarReal(internal->D0));
    SET_VECTOR_ELT(contents, 4, ScalarInteger(internal->dim_Betas));
    SET_VECTOR_ELT(contents, 5, ScalarInteger(internal->dim_changeTimes));
    SET_VECTOR_ELT(contents, 6, ScalarReal(internal->E0));
    SET_VECTOR_ELT(contents, 7, ScalarReal(internal->Gamma));
    SET_VECTOR_ELT(contents, 8, ScalarReal(internal->I0));
    SET_VECTOR_ELT(contents, 9, ScalarReal(internal->initial_D));
    SET_VECTOR_ELT(contents, 10, ScalarReal(internal->initial_E));
    SET_VECTOR_ELT(contents, 11, ScalarReal(internal->initial_I));
    SET_VECTOR_ELT(contents, 12, ScalarReal(internal->initial_R));
    SET_VECTOR_ELT(contents, 13, ScalarReal(internal->initial_S));
    SET_VECTOR_ELT(contents, 15, ScalarReal(internal->Lambda));
    SET_VECTOR_ELT(contents, 16, ScalarReal(internal->R0));
    SET_VECTOR_ELT(contents, 17, ScalarReal(internal->S0));
    SEXP nms = PROTECT(allocVector(STRSXP, 18));
    SET_STRING_ELT(nms, 0, mkChar("Alpha"));
    SET_STRING_ELT(nms, 1, mkChar("Betas"));
    SET_STRING_ELT(nms, 2, mkChar("changeTimes"));
    SET_STRING_ELT(nms, 3, mkChar("D0"));
    SET_STRING_ELT(nms, 4, mkChar("dim_Betas"));
    SET_STRING_ELT(nms, 5, mkChar("dim_changeTimes"));
    SET_STRING_ELT(nms, 6, mkChar("E0"));
    SET_STRING_ELT(nms, 7, mkChar("Gamma"));
    SET_STRING_ELT(nms, 8, mkChar("I0"));
    SET_STRING_ELT(nms, 9, mkChar("initial_D"));
    SET_STRING_ELT(nms, 10, mkChar("initial_E"));
    SET_STRING_ELT(nms, 11, mkChar("initial_I"));
    SET_STRING_ELT(nms, 12, mkChar("initial_R"));
    SET_STRING_ELT(nms, 13, mkChar("initial_S"));
    SET_STRING_ELT(nms, 14, mkChar("interpolate_Beta"));
    SET_STRING_ELT(nms, 15, mkChar("Lambda"));
    SET_STRING_ELT(nms, 16, mkChar("R0"));
    SET_STRING_ELT(nms, 17, mkChar("S0"));
    setAttrib(contents, R_NamesSymbol, nms);
    UNPROTECT(4);
    return contents;
  }
  SEXP seirdGenerator_set_user(SEXP internal_p, SEXP user) {
    seirdGenerator_internal *internal = seirdGenerator_get_internal(internal_p, 1);
    internal->Alpha = user_get_scalar_double(user, "Alpha", internal->Alpha, NA_REAL, NA_REAL);
    internal->Betas = (double*) user_get_array_dim(user, false, internal->Betas, "Betas", 1, NA_REAL, NA_REAL, &internal->dim_Betas);
    internal->D0 = user_get_scalar_double(user, "D0", internal->D0, NA_REAL, NA_REAL);
    internal->E0 = user_get_scalar_double(user, "E0", internal->E0, NA_REAL, NA_REAL);
    internal->Gamma = user_get_scalar_double(user, "Gamma", internal->Gamma, NA_REAL, NA_REAL);
    internal->I0 = user_get_scalar_double(user, "I0", internal->I0, NA_REAL, NA_REAL);
    internal->Lambda = user_get_scalar_double(user, "Lambda", internal->Lambda, NA_REAL, NA_REAL);
    internal->R0 = user_get_scalar_double(user, "R0", internal->R0, NA_REAL, NA_REAL);
    internal->S0 = user_get_scalar_double(user, "S0", internal->S0, NA_REAL, NA_REAL);
    internal->initial_D = internal->D0;
    internal->initial_E = internal->E0;
    internal->initial_I = internal->I0;
    internal->initial_R = internal->R0;
    internal->initial_S = internal->S0;
    internal->dim_changeTimes = internal->dim_Betas;
    internal->changeTimes = (double*) user_get_array(user, false, internal->changeTimes, "changeTimes", NA_REAL, NA_REAL, 1, internal->dim_changeTimes);
    interpolate_check_y(internal->dim_changeTimes, internal->dim_Betas, 0, "Betas", "Beta");
    cinterpolate_free(internal->interpolate_Beta);
    internal->interpolate_Beta = cinterpolate_alloc("constant", internal->dim_changeTimes, 1, internal->changeTimes, internal->Betas, true, false);
    return R_NilValue;
  }
  SEXP seirdGenerator_set_initial(SEXP internal_p, SEXP t_ptr, SEXP state_ptr, SEXP seirdGenerator_use_dde_ptr) {
    return R_NilValue;
  }
  SEXP seirdGenerator_metadata(SEXP internal_p) {
    seirdGenerator_internal *internal = seirdGenerator_get_internal(internal_p, 1);
    SEXP ret = PROTECT(allocVector(VECSXP, 4));
    SEXP nms = PROTECT(allocVector(STRSXP, 4));
    SET_STRING_ELT(nms, 0, mkChar("variable_order"));
    SET_STRING_ELT(nms, 1, mkChar("output_order"));
    SET_STRING_ELT(nms, 2, mkChar("n_out"));
    SET_STRING_ELT(nms, 3, mkChar("interpolate_t"));
    setAttrib(ret, R_NamesSymbol, nms);
    SEXP variable_length = PROTECT(allocVector(VECSXP, 5));
    SEXP variable_names = PROTECT(allocVector(STRSXP, 5));
    setAttrib(variable_length, R_NamesSymbol, variable_names);
    SET_VECTOR_ELT(variable_length, 0, R_NilValue);
    SET_VECTOR_ELT(variable_length, 1, R_NilValue);
    SET_VECTOR_ELT(variable_length, 2, R_NilValue);
    SET_VECTOR_ELT(variable_length, 3, R_NilValue);
    SET_VECTOR_ELT(variable_length, 4, R_NilValue);
    SET_STRING_ELT(variable_names, 0, mkChar("S"));
    SET_STRING_ELT(variable_names, 1, mkChar("E"));
    SET_STRING_ELT(variable_names, 2, mkChar("I"));
    SET_STRING_ELT(variable_names, 3, mkChar("R"));
    SET_STRING_ELT(variable_names, 4, mkChar("D"));
    SET_VECTOR_ELT(ret, 0, variable_length);
    UNPROTECT(2);
    SEXP output_length = PROTECT(allocVector(VECSXP, 1));
    SEXP output_names = PROTECT(allocVector(STRSXP, 1));
    setAttrib(output_length, R_NamesSymbol, output_names);
    SET_VECTOR_ELT(output_length, 0, R_NilValue);
    SET_STRING_ELT(output_names, 0, mkChar("Beta"));
    SET_VECTOR_ELT(ret, 1, output_length);
    UNPROTECT(2);
    SET_VECTOR_ELT(ret, 2, ScalarInteger(1));
    SEXP interpolate_t = PROTECT(allocVector(VECSXP, 3));
    SEXP interpolate_t_nms = PROTECT(allocVector(STRSXP, 3));
    setAttrib(interpolate_t, R_NamesSymbol, interpolate_t_nms);
    SET_VECTOR_ELT(interpolate_t, 0, ScalarReal(internal->changeTimes[0]));
    SET_VECTOR_ELT(interpolate_t, 1, ScalarReal(R_PosInf));
    SET_STRING_ELT(interpolate_t_nms, 0, mkChar("min"));
    SET_STRING_ELT(interpolate_t_nms, 1, mkChar("max"));
    SET_VECTOR_ELT(ret, 3, interpolate_t);
    UNPROTECT(2);
    UNPROTECT(2);
    return ret;
  }
  SEXP seirdGenerator_initial_conditions(SEXP internal_p, SEXP t_ptr) {
    seirdGenerator_internal *internal = seirdGenerator_get_internal(internal_p, 1);
    SEXP r_state = PROTECT(allocVector(REALSXP, 5));
    double * state = REAL(r_state);
    state[0] = internal->initial_S;
    state[1] = internal->initial_E;
    state[2] = internal->initial_I;
    state[3] = internal->initial_R;
    state[4] = internal->initial_D;
    UNPROTECT(1);
    return r_state;
  }
  void seirdGenerator_rhs(seirdGenerator_internal* internal, double t, double * state, double * dstatedt, double * output) {
    double S = state[0];
    double E = state[1];
    double I = state[2];
    double R = state[3];
    double D = state[4];
    dstatedt[4] = I * internal->Alpha;
    dstatedt[2] = E * internal->Lambda - I * internal->Gamma - I * internal->Alpha;
    dstatedt[3] = I * internal->Gamma;
    double N = S + E + I + R + D;
    double Beta = 0.0;
    cinterpolate_eval(t, internal->interpolate_Beta, &Beta);
    dstatedt[1] = Beta * S * I / (double) N - E * internal->Lambda;
    dstatedt[0] = -(Beta) * S * I / (double) N;
    if (output) {
      output[0] = Beta;
    }
  }
  void seirdGenerator_rhs_dde(size_t neq, double t, double * state, double * dstatedt, void * internal) {
    seirdGenerator_rhs((seirdGenerator_internal*)internal, t, state, dstatedt, NULL);
  }
  void seirdGenerator_rhs_desolve(int * neq, double * t, double * state, double * dstatedt, double * output, int * np) {
    seirdGenerator_rhs(seirdGenerator_internal_ds, *t, state, dstatedt, output);
  }
  void seirdGenerator_output_dde(size_t n_eq, double t, double * state, size_t n_output, double * output, void * internal_p) {
    seirdGenerator_internal *internal = (seirdGenerator_internal*) internal_p;
    double Beta = 0.0;
    cinterpolate_eval(t, internal->interpolate_Beta, &Beta);
    output[0] = Beta;
  }
  SEXP seirdGenerator_rhs_r(SEXP internal_p, SEXP t, SEXP state) {
    SEXP dstatedt = PROTECT(allocVector(REALSXP, LENGTH(state)));
    seirdGenerator_internal *internal = seirdGenerator_get_internal(internal_p, 1);
    SEXP output_ptr = PROTECT(allocVector(REALSXP, 1));
    setAttrib(dstatedt, install("output"), output_ptr);
    UNPROTECT(1);
    double *output = REAL(output_ptr);
    seirdGenerator_rhs(internal, scalar_real(t, "t"), REAL(state), REAL(dstatedt), output);
    UNPROTECT(1);
    return dstatedt;
  }
  sirdGenerator_internal* sirdGenerator_get_internal(SEXP internal_p, int closed_error) {
    sirdGenerator_internal *internal = NULL;
    if (TYPEOF(internal_p) != EXTPTRSXP) {
      Rf_error("Expected an external pointer");
    }
    internal = (sirdGenerator_internal*) R_ExternalPtrAddr(internal_p);
    if (!internal && closed_error) {
      Rf_error("Pointer has been invalidated");
    }
    return internal;
  }
  void sirdGenerator_finalise(SEXP internal_p) {
    sirdGenerator_internal *internal = sirdGenerator_get_internal(internal_p, 0);
    if (internal_p) {
      cinterpolate_free(internal->interpolate_Beta);
      internal->interpolate_Beta = NULL;
      Free(internal->Betas);
      Free(internal->changeTimes);
      Free(internal);
      R_ClearExternalPtr(internal_p);
    }
  }
  SEXP sirdGenerator_create(SEXP user) {
    sirdGenerator_internal *internal = (sirdGenerator_internal*) Calloc(1, sirdGenerator_internal);
    internal->Betas = NULL;
    internal->changeTimes = NULL;
    internal->Alpha = NA_REAL;
    internal->Betas = NULL;
    internal->changeTimes = NULL;
    internal->Gamma = NA_REAL;
    internal->I0 = NA_REAL;
    internal->S0 = NA_REAL;
    internal->D0 = 0;
    internal->R0 = 0;
    SEXP ptr = PROTECT(R_MakeExternalPtr(internal, R_NilValue, R_NilValue));
    R_RegisterCFinalizer(ptr, sirdGenerator_finalise);
    UNPROTECT(1);
    return ptr;
  }
  static sirdGenerator_internal *sirdGenerator_internal_ds;
  void sirdGenerator_initmod_desolve(void(* odeparms) (int *, double *)) {
    static DL_FUNC get_desolve_gparms = NULL;
    if (get_desolve_gparms == NULL) {
      get_desolve_gparms =
        R_GetCCallable("deSolve", "get_deSolve_gparms");
    }
    sirdGenerator_internal_ds = sirdGenerator_get_internal(get_desolve_gparms(), 1);
  }
  SEXP sirdGenerator_contents(SEXP internal_p) {
    sirdGenerator_internal *internal = sirdGenerator_get_internal(internal_p, 1);
    SEXP contents = PROTECT(allocVector(VECSXP, 15));
    SET_VECTOR_ELT(contents, 0, ScalarReal(internal->Alpha));
    SEXP Betas = PROTECT(allocVector(REALSXP, internal->dim_Betas));
    memcpy(REAL(Betas), internal->Betas, internal->dim_Betas * sizeof(double));
    SET_VECTOR_ELT(contents, 1, Betas);
    SEXP changeTimes = PROTECT(allocVector(REALSXP, internal->dim_changeTimes));
    memcpy(REAL(changeTimes), internal->changeTimes, internal->dim_changeTimes * sizeof(double));
    SET_VECTOR_ELT(contents, 2, changeTimes);
    SET_VECTOR_ELT(contents, 3, ScalarReal(internal->D0));
    SET_VECTOR_ELT(contents, 4, ScalarInteger(internal->dim_Betas));
    SET_VECTOR_ELT(contents, 5, ScalarInteger(internal->dim_changeTimes));
    SET_VECTOR_ELT(contents, 6, ScalarReal(internal->Gamma));
    SET_VECTOR_ELT(contents, 7, ScalarReal(internal->I0));
    SET_VECTOR_ELT(contents, 8, ScalarReal(internal->initial_D));
    SET_VECTOR_ELT(contents, 9, ScalarReal(internal->initial_I));
    SET_VECTOR_ELT(contents, 10, ScalarReal(internal->initial_R));
    SET_VECTOR_ELT(contents, 11, ScalarReal(internal->initial_S));
    SET_VECTOR_ELT(contents, 13, ScalarReal(internal->R0));
    SET_VECTOR_ELT(contents, 14, ScalarReal(internal->S0));
    SEXP nms = PROTECT(allocVector(STRSXP, 15));
    SET_STRING_ELT(nms, 0, mkChar("Alpha"));
    SET_STRING_ELT(nms, 1, mkChar("Betas"));
    SET_STRING_ELT(nms, 2, mkChar("changeTimes"));
    SET_STRING_ELT(nms, 3, mkChar("D0"));
    SET_STRING_ELT(nms, 4, mkChar("dim_Betas"));
    SET_STRING_ELT(nms, 5, mkChar("dim_changeTimes"));
    SET_STRING_ELT(nms, 6, mkChar("Gamma"));
    SET_STRING_ELT(nms, 7, mkChar("I0"));
    SET_STRING_ELT(nms, 8, mkChar("initial_D"));
    SET_STRING_ELT(nms, 9, mkChar("initial_I"));
    SET_STRING_ELT(nms, 10, mkChar("initial_R"));
    SET_STRING_ELT(nms, 11, mkChar("initial_S"));
    SET_STRING_ELT(nms, 12, mkChar("interpolate_Beta"));
    SET_STRING_ELT(nms, 13, mkChar("R0"));
    SET_STRING_ELT(nms, 14, mkChar("S0"));
    setAttrib(contents, R_NamesSymbol, nms);
    UNPROTECT(4);
    return contents;
  }
  SEXP sirdGenerator_set_user(SEXP internal_p, SEXP user) {
    sirdGenerator_internal *internal = sirdGenerator_get_internal(internal_p, 1);
    internal->Alpha = user_get_scalar_double(user, "Alpha", internal->Alpha, NA_REAL, NA_REAL);
    internal->Betas = (double*) user_get_array_dim(user, false, internal->Betas, "Betas", 1, NA_REAL, NA_REAL, &internal->dim_Betas);
    internal->D0 = user_get_scalar_double(user, "D0", internal->D0, NA_REAL, NA_REAL);
    internal->Gamma = user_get_scalar_double(user, "Gamma", internal->Gamma, NA_REAL, NA_REAL);
    internal->I0 = user_get_scalar_double(user, "I0", internal->I0, NA_REAL, NA_REAL);
    internal->R0 = user_get_scalar_double(user, "R0", internal->R0, NA_REAL, NA_REAL);
    internal->S0 = user_get_scalar_double(user, "S0", internal->S0, NA_REAL, NA_REAL);
    internal->initial_D = internal->D0;
    internal->initial_I = internal->I0;
    internal->initial_R = internal->R0;
    internal->initial_S = internal->S0;
    internal->dim_changeTimes = internal->dim_Betas;
    internal->changeTimes = (double*) user_get_array(user, false, internal->changeTimes, "changeTimes", NA_REAL, NA_REAL, 1, internal->dim_changeTimes);
    interpolate_check_y(internal->dim_changeTimes, internal->dim_Betas, 0, "Betas", "Beta");
    cinterpolate_free(internal->interpolate_Beta);
    internal->interpolate_Beta = cinterpolate_alloc("constant", internal->dim_changeTimes, 1, internal->changeTimes, internal->Betas, true, false);
    return R_NilValue;
  }
  SEXP sirdGenerator_set_initial(SEXP internal_p, SEXP t_ptr, SEXP state_ptr, SEXP sirdGenerator_use_dde_ptr) {
    return R_NilValue;
  }
  SEXP sirdGenerator_metadata(SEXP internal_p) {
    sirdGenerator_internal *internal = sirdGenerator_get_internal(internal_p, 1);
    SEXP ret = PROTECT(allocVector(VECSXP, 4));
    SEXP nms = PROTECT(allocVector(STRSXP, 4));
    SET_STRING_ELT(nms, 0, mkChar("variable_order"));
    SET_STRING_ELT(nms, 1, mkChar("output_order"));
    SET_STRING_ELT(nms, 2, mkChar("n_out"));
    SET_STRING_ELT(nms, 3, mkChar("interpolate_t"));
    setAttrib(ret, R_NamesSymbol, nms);
    SEXP variable_length = PROTECT(allocVector(VECSXP, 4));
    SEXP variable_names = PROTECT(allocVector(STRSXP, 4));
    setAttrib(variable_length, R_NamesSymbol, variable_names);
    SET_VECTOR_ELT(variable_length, 0, R_NilValue);
    SET_VECTOR_ELT(variable_length, 1, R_NilValue);
    SET_VECTOR_ELT(variable_length, 2, R_NilValue);
    SET_VECTOR_ELT(variable_length, 3, R_NilValue);
    SET_STRING_ELT(variable_names, 0, mkChar("S"));
    SET_STRING_ELT(variable_names, 1, mkChar("I"));
    SET_STRING_ELT(variable_names, 2, mkChar("R"));
    SET_STRING_ELT(variable_names, 3, mkChar("D"));
    SET_VECTOR_ELT(ret, 0, variable_length);
    UNPROTECT(2);
    SEXP output_length = PROTECT(allocVector(VECSXP, 1));
    SEXP output_names = PROTECT(allocVector(STRSXP, 1));
    setAttrib(output_length, R_NamesSymbol, output_names);
    SET_VECTOR_ELT(output_length, 0, R_NilValue);
    SET_STRING_ELT(output_names, 0, mkChar("Beta"));
    SET_VECTOR_ELT(ret, 1, output_length);
    UNPROTECT(2);
    SET_VECTOR_ELT(ret, 2, ScalarInteger(1));
    SEXP interpolate_t = PROTECT(allocVector(VECSXP, 3));
    SEXP interpolate_t_nms = PROTECT(allocVector(STRSXP, 3));
    setAttrib(interpolate_t, R_NamesSymbol, interpolate_t_nms);
    SET_VECTOR_ELT(interpolate_t, 0, ScalarReal(internal->changeTimes[0]));
    SET_VECTOR_ELT(interpolate_t, 1, ScalarReal(R_PosInf));
    SET_STRING_ELT(interpolate_t_nms, 0, mkChar("min"));
    SET_STRING_ELT(interpolate_t_nms, 1, mkChar("max"));
    SET_VECTOR_ELT(ret, 3, interpolate_t);
    UNPROTECT(2);
    UNPROTECT(2);
    return ret;
  }
  SEXP sirdGenerator_initial_conditions(SEXP internal_p, SEXP t_ptr) {
    sirdGenerator_internal *internal = sirdGenerator_get_internal(internal_p, 1);
    SEXP r_state = PROTECT(allocVector(REALSXP, 4));
    double * state = REAL(r_state);
    state[0] = internal->initial_S;
    state[1] = internal->initial_I;
    state[2] = internal->initial_R;
    state[3] = internal->initial_D;
    UNPROTECT(1);
    return r_state;
  }
  void sirdGenerator_rhs(sirdGenerator_internal* internal, double t, double * state, double * dstatedt, double * output) {
    double S = state[0];
    double I = state[1];
    double R = state[2];
    double D = state[3];
    dstatedt[3] = I * internal->Alpha;
    dstatedt[2] = I * internal->Gamma;
    double N = S + I + R + D;
    double Beta = 0.0;
    cinterpolate_eval(t, internal->interpolate_Beta, &Beta);
    dstatedt[1] = Beta * S * I / (double) N - I * internal->Gamma - I * internal->Alpha;
    dstatedt[0] = -(Beta) * S * I / (double) N;
    if (output) {
      output[0] = Beta;
    }
  }
  void sirdGenerator_rhs_dde(size_t neq, double t, double * state, double * dstatedt, void * internal) {
    sirdGenerator_rhs((sirdGenerator_internal*)internal, t, state, dstatedt, NULL);
  }
  void sirdGenerator_rhs_desolve(int * neq, double * t, double * state, double * dstatedt, double * output, int * np) {
    sirdGenerator_rhs(sirdGenerator_internal_ds, *t, state, dstatedt, output);
  }
  void sirdGenerator_output_dde(size_t n_eq, double t, double * state, size_t n_output, double * output, void * internal_p) {
    sirdGenerator_internal *internal = (sirdGenerator_internal*) internal_p;
    double Beta = 0.0;
    cinterpolate_eval(t, internal->interpolate_Beta, &Beta);
    output[0] = Beta;
  }
  SEXP sirdGenerator_rhs_r(SEXP internal_p, SEXP t, SEXP state) {
    SEXP dstatedt = PROTECT(allocVector(REALSXP, LENGTH(state)));
    sirdGenerator_internal *internal = sirdGenerator_get_internal(internal_p, 1);
    SEXP output_ptr = PROTECT(allocVector(REALSXP, 1));
    setAttrib(dstatedt, install("output"), output_ptr);
    UNPROTECT(1);
    double *output = REAL(output_ptr);
    sirdGenerator_rhs(internal, scalar_real(t, "t"), REAL(state), REAL(dstatedt), output);
    UNPROTECT(1);
    return dstatedt;
  }
  double user_get_scalar_double(SEXP user, const char *name,
                                double default_value, double min, double max) {
    double ret = default_value;
    SEXP el = user_list_element(user, name);
    if (el != R_NilValue) {
      if (length(el) != 1) {
        Rf_error("Expected a scalar numeric for '%s'", name);
      }
      if (TYPEOF(el) == REALSXP) {
        ret = REAL(el)[0];
      } else if (TYPEOF(el) == INTSXP) {
        ret = INTEGER(el)[0];
      } else {
        Rf_error("Expected a numeric value for %s", name);
      }
    }
    if (ISNA(ret)) {
      Rf_error("Expected a value for '%s'", name);
    }
    user_check_values_double(&ret, 1, min, max, name);
    return ret;
  }
  int user_get_scalar_int(SEXP user, const char *name,
                          int default_value, double min, double max) {
    int ret = default_value;
    SEXP el = user_list_element(user, name);
    if (el != R_NilValue) {
      if (length(el) != 1) {
        Rf_error("Expected scalar integer for %d", name);
      }
      if (TYPEOF(el) == REALSXP) {
        double tmp = REAL(el)[0];
        if (fabs(tmp - round(tmp)) > 2e-8) {
          Rf_error("Expected '%s' to be integer-like", name);
        }
      }
      ret = INTEGER(coerceVector(el, INTSXP))[0];
    }
    if (ret == NA_INTEGER) {
      Rf_error("Expected a value for '%s'", name);
    }
    user_check_values_int(&ret, 1, min, max, name);
    return ret;
  }
  void user_check_values_double(double * value, size_t len,
                                    double min, double max, const char *name) {
    for (size_t i = 0; i < len; ++i) {
      if (ISNA(value[i])) {
        Rf_error("'%s' must not contain any NA values", name);
      }
    }
    if (min != NA_REAL) {
      for (size_t i = 0; i < len; ++i) {
        if (value[i] < min) {
          Rf_error("Expected '%s' to be at least %g", name, min);
        }
      }
    }
    if (max != NA_REAL) {
      for (size_t i = 0; i < len; ++i) {
        if (value[i] > max) {
          Rf_error("Expected '%s' to be at most %g", name, max);
        }
      }
    }
  }
  void user_check_values_int(int * value, size_t len,
                                 double min, double max, const char *name) {
    for (size_t i = 0; i < len; ++i) {
      if (ISNA(value[i])) {
        Rf_error("'%s' must not contain any NA values", name);
      }
    }
    if (min != NA_REAL) {
      for (size_t i = 0; i < len; ++i) {
        if (value[i] < min) {
          Rf_error("Expected '%s' to be at least %g", name, min);
        }
      }
    }
    if (max != NA_REAL) {
      for (size_t i = 0; i < len; ++i) {
        if (value[i] > max) {
          Rf_error("Expected '%s' to be at most %g", name, max);
        }
      }
    }
  }
  void user_check_values(SEXP value, double min, double max,
                             const char *name) {
    size_t len = (size_t)length(value);
    if (TYPEOF(value) == INTSXP) {
      user_check_values_int(INTEGER(value), len, min, max, name);
    } else {
      user_check_values_double(REAL(value), len, min, max, name);
    }
  }
  SEXP user_list_element(SEXP list, const char *name) {
    SEXP ret = R_NilValue, names = getAttrib(list, R_NamesSymbol);
    for (int i = 0; i < length(list); ++i) {
      if (strcmp(CHAR(STRING_ELT(names, i)), name) == 0) {
        ret = VECTOR_ELT(list, i);
        break;
      }
    }
    return ret;
  }
  void odin_set_dim(SEXP target, int rank, ...) {
    SEXP r_dim = PROTECT(allocVector(INTSXP, rank));
    int *dim = INTEGER(r_dim);
  
    va_list ap;
    va_start(ap, rank);
    for (size_t i = 0; i < (size_t)rank; ++i) {
      dim[i] = va_arg(ap, int);
    }
    va_end(ap);
  
    setAttrib(target, R_DimSymbol, r_dim);
    UNPROTECT(1);
  }
  void* user_get_array_dim(SEXP user, bool is_integer, void * previous,
                           const char *name, int rank,
                           double min, double max, int *dest_dim) {
    SEXP el = user_get_array_check_rank(user, name, rank, previous == NULL);
    if (el == R_NilValue) {
      return previous;
    }
  
    dest_dim[0] = LENGTH(el);
    if (rank > 1) {
      SEXP r_dim = PROTECT(coerceVector(getAttrib(el, R_DimSymbol), INTSXP));
      int *dim = INTEGER(r_dim);
  
      for (size_t i = 0; i < (size_t) rank; ++i) {
        dest_dim[i + 1] = dim[i];
      }
  
      UNPROTECT(1);
    }
  
    el = PROTECT(user_get_array_check(el, is_integer, name, min, max));
  
    int len = LENGTH(el);
    void *dest = NULL;
    if (is_integer) {
      dest = Calloc(len, int);
      memcpy(dest, INTEGER(el), len * sizeof(int));
    } else {
      dest = Calloc(len, double);
      memcpy(dest, REAL(el), len * sizeof(double));
    }
    Free(previous);
  
    UNPROTECT(1);
  
    return dest;
  }
  void* user_get_array(SEXP user, bool is_integer, void * previous,
                       const char *name, double min, double max,
                       int rank, ...) {
    SEXP el = user_get_array_check_rank(user, name, rank, previous == NULL);
    if (el == R_NilValue) {
      return previous;
    }
  
    SEXP r_dim;
    int *dim;
  
    size_t len = LENGTH(el);
    if (rank == 1) {
      r_dim = PROTECT(ScalarInteger(len));
    } else {
      r_dim = PROTECT(coerceVector(getAttrib(el, R_DimSymbol), INTSXP));
    }
    dim = INTEGER(r_dim);
  
    va_list ap;
    va_start(ap, rank);
    for (size_t i = 0; i < (size_t) rank; ++i) {
      int dim_expected = va_arg(ap, int);
      if (dim[i] != dim_expected) {
        va_end(ap); // avoid a leak
        if (rank == 1) {
          Rf_error("Expected length %d value for %s", dim_expected, name);
        } else {
          Rf_error("Incorrect size of dimension %d of %s (expected %d)",
                   i + 1, name, dim_expected);
        }
      }
    }
    va_end(ap);
    UNPROTECT(1);
  
    el = PROTECT(user_get_array_check(el, is_integer, name, min, max));
  
    void *dest = NULL;
    if (is_integer) {
      dest = Calloc(len, int);
      memcpy(dest, INTEGER(el), len * sizeof(int));
    } else {
      dest = Calloc(len, double);
      memcpy(dest, REAL(el), len * sizeof(double));
    }
    Free(previous);
  
    UNPROTECT(1);
  
    return dest;
  }
  SEXP user_get_array_check(SEXP el, bool is_integer, const char *name,
                            double min, double max) {
    size_t len = (size_t) length(el);
    if (is_integer) {
      if (TYPEOF(el) == INTSXP) {
        user_check_values_int(INTEGER(el), len, min, max, name);
      } else if (TYPEOF(el) == REALSXP) {
        el = PROTECT(coerceVector(el, INTSXP));
        user_check_values_int(INTEGER(el), len, min, max, name);
        UNPROTECT(1);
      } else {
        Rf_error("Expected a integer value for %s", name);
      }
    } else {
      if (TYPEOF(el) == INTSXP) {
        el = PROTECT(coerceVector(el, REALSXP));
        user_check_values_double(REAL(el), len, min, max, name);
        UNPROTECT(1);
      } else if (TYPEOF(el) == REALSXP) {
        user_check_values_double(REAL(el), len, min, max, name);
      } else {
        Rf_error("Expected a numeric value for %s", name);
      }
    }
    return el;
  }
  SEXP user_get_array_check_rank(SEXP user, const char *name, int rank,
                                 bool required) {
    SEXP el = user_list_element(user, name);
    if (el == R_NilValue) {
      if (required) {
        Rf_error("Expected a value for '%s'", name);
      }
    } else {
      if (rank == 1) {
        if (isArray(el)) {
          Rf_error("Expected a numeric vector for '%s'", name);
        }
      } else {
        SEXP r_dim = getAttrib(el, R_DimSymbol);
        if (r_dim == R_NilValue || LENGTH(r_dim) != rank) {
          if (rank == 2) {
            Rf_error("Expected a numeric matrix for '%s'", name);
          } else {
            Rf_error("Expected a numeric array of rank %d for '%s'", rank, name);
          }
        }
      }
    }
    return el;
  }
  void interpolate_check_y(size_t nx, size_t ny, size_t i, const char *name_arg, const char *name_target) {
    if (nx != ny) {
      if (i == 0) {
        // vector case
        Rf_error("Expected %s to have length %d (for %s)",
                 name_arg, nx, name_target);
      } else {
        // array case
        Rf_error("Expected dimension %d of %s to have size %d (for %s)",
                 i, name_arg, nx, name_target);
      }
    }
  }
  double scalar_real(SEXP x, const char * name) {
    if (Rf_length(x) != 1) {
      Rf_error("Expected a scalar for %s", name);
    }
    double ret = 0.0;
    if (TYPEOF(x) == INTSXP) {
      ret = INTEGER(x)[0];
    } else if (TYPEOF(x) == REALSXP) {
      ret = REAL(x)[0];
    } else {
      Rf_error("Expected a numeric value for %s", name);
    }
    return ret;
  }
  double odin_sum1(double *x, size_t from, size_t to) {
    double tot = 0.0;
    for (size_t i = from; i < to; ++i) {
      tot += x[i];
    }
    return tot;
  }
  double odin_sum2(double* x, int from_i, int to_i, int from_j, int to_j, int dim_x_1) {
    double tot = 0.0;
    for (int j = from_j; j < to_j; ++j) {
      int jj = j * dim_x_1;
      for (int i = from_i; i < to_i; ++i) {
        tot += x[i + jj];
      }
    }
    return tot;
  }
  double odin_sum3(double* x, int from_i, int to_i, int from_j, int to_j, int from_k, int to_k, int dim_x_1, int dim_x_12) {
    double tot = 0.0;
    for (int k = from_k; k < to_k; ++k) {
      int kk = k * dim_x_12;
      for (int j = from_j; j < to_j; ++j) {
        int jj = j * dim_x_1 + kk;
        for (int i = from_i; i < to_i; ++i) {
          tot += x[i + jj];
        }
      }
    }
    return tot;
  }
  // This construction is to help odin
  #ifndef CINTERPOLTE_CINTERPOLATE_H_
  #endif
  
  
  void * cinterpolate_alloc(const char *type, size_t n, size_t ny,
                            double *x, double *y, bool fail_on_extrapolate,
                            bool auto_clean) {
    typedef void* interpolate_alloc_t(const char *, size_t, size_t,
                                      double*, double*, bool, bool);
    static interpolate_alloc_t *fun;
    if (fun == NULL) {
      fun = (interpolate_alloc_t*)
        R_GetCCallable("cinterpolate", "interpolate_alloc");
    }
    return fun(type, n, ny, x, y, fail_on_extrapolate, auto_clean);
  }
  
  
  int cinterpolate_eval(double x, void *obj, double *y) {
    typedef int interpolate_eval_t(double, void*, double*);
    static interpolate_eval_t *fun;
    if (fun == NULL) {
      fun = (interpolate_eval_t*)
        R_GetCCallable("cinterpolate", "interpolate_eval");
    }
    return fun(x, obj, y);
  }
  
  
  void cinterpolate_free(void *obj) {
    typedef int interpolate_free_t(void*);
    static interpolate_free_t *fun;
    if (fun == NULL) {
      fun = (interpolate_free_t*)
        R_GetCCallable("cinterpolate", "interpolate_free");
    }
    fun(obj);
  }
  