fixed sortTimes code

sortTimes now takes dataframe and a column vector as an input

DESCRIPTION

@@ -23,7 +23,7 @@ Maintainer: Satyaprakash Nayak <[email protected]>
 URL: https://github.com/sn248/sundialr
 BugReports: https://github.com/sn248/sundialr/issues
 Description: Provides a way to call the functions in 'SUNDIALS' C ODE solving library (<https://computation.llnl.gov/projects/sundials>). Currently the serial version of ODE solver, 'CVODE' and sensitivity calculator 'CVODES' from the 'SUNDIALS' library are implemented. The package requires ODE to be written as an 'R' or 'Rcpp' function and does not require the 'SUNDIALS' library to be installed on the local machine.
-License: GPL (>=2)
+License: BSD_3_clause + file LICENSE
 LazyData: TRUE
 Imports:
     Rcpp (>= 0.12.5)

LICENSE

@@ -0,0 +1,3 @@
+YEAR: 2020
+COPYRIGHT HOLDER: Satyaprakash Nayak
+ORGANIZATION:

R/RcppExports.R

@@ -31,6 +31,23 @@ cvodes <- function(time_vector, IC, input_function, Parameters, reltolerance = 0
     .Call('_sundialr_cvodes', PACKAGE = 'sundialr', time_vector, IC, input_function, Parameters, reltolerance, abstolerance, SensType, ErrCon)
 }
 
+#'cvsolve
+#'
+#'solver to solve stiff ODEs with discontinuties
+#'@param time_vector time vector
+#'@param IC Initial Conditions
+#'@param input_function Right Hand Side function of ODEs
+#'@param Parameters Parameters input to ODEs
+#'@param events Discontinuities in the solution (a DataFrame, default value is NULL)
+#'@param Jacobian User-supplied Jacobian(a matrix, default value is NULL)
+#'@param reltolerance Relative Tolerance (a scalar, default value  = 1e-04)
+#'@param abstolerance Absolute Tolerance (a scalar or vector with length equal to ydot, default = 1e-04)
+#'@param constraints By default all the solution values are constrained to be non-negative
+#'@example /inst/examples/cv_Roberts_dns.r
+cvsolve <- function(time_vector, IC, input_function, Parameters, events = NULL, Jacobian = NULL, reltolerance = 0.0001, abstolerance = 0.0001) {
+    .Call('_sundialr_cvsolve', PACKAGE = 'sundialr', time_vector, IC, input_function, Parameters, events, Jacobian, reltolerance, abstolerance)
+}
+
 #'ida
 #'
 #' IDA solver to solve stiff DAEs

inst/include/rhs_func.h

@@ -0,0 +1,14 @@
+// File: rhs_func.h
+
+#ifndef RHS_FUNC_H
+#define RHS_FUNC_H
+
+// struct to use if R or Rcpp function is input as RHS function
+struct rhs_func{
+  Rcpp::Function rhs_eqn;
+  Rcpp::NumericVector params;
+};
+
+int rhs_function(realtype t, N_Vector y, N_Vector ydot, void* user_data);
+
+#endif /* rhs_func */

src/Makevars

@@ -5,7 +5,9 @@ LIBS=-L./ -L../inst/
 
 .PHONY: all ../inst/libsundials_all.a
 
-SOURCES= check_retval.cpp cvode.cpp cvodes.cpp ida.cpp RcppExports.cpp
+SOURCES= check_retval.cpp cvode.cpp cvodes.cpp \
+ida.cpp cvsolve.cpp rhs_func.cpp RcppExports.cpp
+
 OBJECTS= $(SOURCES:.cpp=.o)
 
 #SOURCES_CVODE= ./sundials/cvode/cvode.c ./sundials/cvode/cvode_bandpre.c \

src/RcppExports.cpp

@@ -40,6 +40,24 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cvsolve
+NumericMatrix cvsolve(NumericVector time_vector, NumericVector IC, SEXP input_function, NumericVector Parameters, NumericMatrix events, NumericMatrix Jacobian, double reltolerance, NumericVector abstolerance);
+RcppExport SEXP _sundialr_cvsolve(SEXP time_vectorSEXP, SEXP ICSEXP, SEXP input_functionSEXP, SEXP ParametersSEXP, SEXP eventsSEXP, SEXP JacobianSEXP, SEXP reltoleranceSEXP, SEXP abstoleranceSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericVector >::type time_vector(time_vectorSEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type IC(ICSEXP);
+    Rcpp::traits::input_parameter< SEXP >::type input_function(input_functionSEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type Parameters(ParametersSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type events(eventsSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type Jacobian(JacobianSEXP);
+    Rcpp::traits::input_parameter< double >::type reltolerance(reltoleranceSEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type abstolerance(abstoleranceSEXP);
+    rcpp_result_gen = Rcpp::wrap(cvsolve(time_vector, IC, input_function, Parameters, events, Jacobian, reltolerance, abstolerance));
+    return rcpp_result_gen;
+END_RCPP
+}
 // ida
 NumericMatrix ida(NumericVector time_vector, NumericVector IC, NumericVector IRes, SEXP input_function, NumericVector Parameters, double reltolerance, NumericVector abstolerance);
 RcppExport SEXP _sundialr_ida(SEXP time_vectorSEXP, SEXP ICSEXP, SEXP IResSEXP, SEXP input_functionSEXP, SEXP ParametersSEXP, SEXP reltoleranceSEXP, SEXP abstoleranceSEXP) {
@@ -61,6 +79,7 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_sundialr_cvode", (DL_FUNC) &_sundialr_cvode, 6},
     {"_sundialr_cvodes", (DL_FUNC) &_sundialr_cvodes, 8},
+    {"_sundialr_cvsolve", (DL_FUNC) &_sundialr_cvsolve, 8},
     {"_sundialr_ida", (DL_FUNC) &_sundialr_ida, 7},
     {NULL, NULL, 0}
 };

src/cvode.cpp

@@ -21,64 +21,12 @@
 #include <sunlinsol/sunlinsol_dense.h>
 
 #include <check_retval.h>
+#include <rhs_func.h>
 
 #define Ith(v,i)    NV_Ith_S(v,i-1)         /* i-th vector component i=1..NEQ */
 
 using namespace Rcpp;
 
-// struct to use if R or Rcpp function is input as RHS function
-struct rhs_func{
-  Function rhs_eqn;
-  NumericVector params;
-};
-
-// function called by CVodeInit if user inputs R function
-int rhs_function(realtype t, N_Vector y, N_Vector ydot, void* user_data){
-
-  // convert y to NumericVector y1
-  int y_len = NV_LENGTH_S(y);
-
-  NumericVector y1(y_len);    // filled with zeros
-  // realtype *y_ptr = N_VGetArrayPointer(y);
-  for (int i = 0; i < y_len; i++){
-    y1[i] = Ith(y,i+1); // y_ptr[i];
-  }
-
-  NumericVector ydot1(y_len);    // filled with zeros
-
-  if(!user_data){
-    stop("Something went wrong, stopping!");
-  }
-
-  // cast void pointer to pointer to struct and assign rhs to a Function
-  struct rhs_func *my_rhs_fun = (struct rhs_func*)user_data;
-
-  if(my_rhs_fun){
-    Function rhs_fun = (*my_rhs_fun).rhs_eqn;           // function
-    NumericVector p_values = (*my_rhs_fun).params;      // rate parameters
-
-    if (rhs_fun && (TYPEOF(rhs_fun) == CLOSXP)){
-      // use the function to calculate value of RHS ----
-      ydot1 = rhs_fun(t, y1, p_values);
-    }
-    else{
-      stop("Something went wrong, stopping!");
-    }
-  }
-  else {
-    stop("Something went wrong, stopping!");
-  }
-
-  // convert NumericVector ydot1 to N_Vector ydot
-  realtype *ydot_ptr = N_VGetArrayPointer(ydot);
-  for (int i = 0; i<  y_len; i++){
-    ydot_ptr[i] = ydot1[i];
-  }
-
-  // everything went smoothly
-  return(0);
-}
-//---RHS function definition ends ----------------------------------------------
 
 //------------------------------------------------------------------------------
 //'cvode
@@ -119,27 +67,26 @@ NumericMatrix cvode(NumericVector time_vector, NumericVector IC, SEXP input_func
     }
   }
   else if (abstol_len == y_len){
-
     for (int i = 0; i<abstol_len; i++){
       abstol_ptr[i] = abstolerance[i];
     }
   }
   else if(abstol_len != 1 || abstol_len != y_len){
-
     stop("Absolute tolerance must be a scalar or a vector of same length as IC \n");
   }
-  //----------------------------------------------------------------------------
-  // // Set the initial conditions-------------------------------------------------
+
+  // Set the initial conditions-------------------------------------------------
   N_Vector y0 = N_VNew_Serial(y_len);
   realtype *y0_ptr = N_VGetArrayPointer(y0);
   for (int i = 0; i<y_len; i++){
     y0_ptr[i] = IC[i]; // NV_Ith_S(y0, i)
   }
-  // //----------------------------------------------------------------------------
+
+  // Call CVodeCreate to create the solver memory and specify the Backward Differentiation Formula
   void *cvode_mem;
   cvode_mem = NULL;
 
-  // // Call CVodeCreate to create the solver memory and specify the Backward Differentiation Formula
+
   cvode_mem = CVodeCreate(CV_BDF);
   if (check_retval((void *) cvode_mem, "CVodeCreate", 0)) { stop("Stopping cvode!"); }
 
@@ -164,7 +111,6 @@ NumericMatrix cvode(NumericVector time_vector, NumericVector IC, SEXP input_func
     flag = CVodeSVtolerances(cvode_mem, reltol, abstol);
     if (check_retval(&flag, "CVodeSVtolerances", 1)) { stop("Stopping cvode, something went wrong in setting solver tolerances!"); }
 
-    //--required in the new version ----------------------------------------------
     // Create dense SUNMatrix for use in linear solves
     sunindextype y_len_M = y_len;
     SUNMatrix SM = SUNDenseMatrix(y_len_M, y_len_M);
@@ -204,7 +150,7 @@ NumericMatrix cvode(NumericVector time_vector, NumericVector IC, SEXP input_func
 
       flag = CVode(cvode_mem, tout, y0, &time, CV_NORMAL);
 
-      if (check_retval(&flag, "CVode", 1)) { stop("Stopping CVODE, something went wrong in solving the system of ODEs!"); break; } // Something went wrong in solving it!
+      if (check_retval(&flag, "CVode", 1)) { stop("Stopping CVODE, something went wrong in solving the system of ODEs!"); } // Something went wrong in solving it!
       if (flag == CV_SUCCESS) {
 
         // store results in soln matrix
@@ -221,10 +167,8 @@ NumericMatrix cvode(NumericVector time_vector, NumericVector IC, SEXP input_func
 
     // free integrator memory
     CVodeFree(&cvode_mem);
-
     // free the linear solver memory
     SUNLinSolFree(LS);
-
     // Free the matrix memory
     SUNMatDestroy(SM);
 
@@ -239,4 +183,4 @@ NumericMatrix cvode(NumericVector time_vector, NumericVector IC, SEXP input_func
   }
 
 }
-// //--- CVODE definition ends ----------------------------------------------------
+//--- cvode definition ends ----------------------------------------------------

src/cvodes.cpp

@@ -145,7 +145,6 @@ NumericMatrix cvodes(NumericVector time_vector, NumericVector IC, SEXP input_fun
   int abstol_len = abstolerance.length();
 
   int flag;
-  // realtype reltol = reltolerance;
 
   realtype T0 = RCONST(time_vector[0]);     //RCONST(0.0);  // Initial Time
 
@@ -206,7 +205,6 @@ NumericMatrix cvodes(NumericVector time_vector, NumericVector IC, SEXP input_fun
   cvode_mem = CVodeCreate(CV_BDF);
   if (check_retval((void *) cvode_mem, "CVodeCreate", 0)) { stop("Stopping cvodes, cannot allocate memory for CVODES!"); }
 
-
   //-- assign user input to the struct based on SEXP type of input_function
   if (!input_function){
     stop("Something is wrong with input function, stopping!");
@@ -272,7 +270,7 @@ NumericMatrix cvodes(NumericVector time_vector, NumericVector IC, SEXP input_fun
 
     /* Call CVodeSetSensParams to specify problem parameter information for
      sensitivity calculations */
-    struct rhs_func_sens *ptr = &my_rhs_function;        // struct UserData storing my data
+    // struct rhs_func_sens *ptr = &my_rhs_function;        // struct UserData storing my data
     // p = (my_rhs_function.params).begin();
     // Rcout << (my_rhs_function.params).begin() << "\n";
     // Rcout << &(ptr->params[0]);