Add unit tests, Speed up dexp1, dexp2, sinc, sinchw

[katrinabrock]

Switching out ifelse for if(...){...}else{...}. Added test to ensure outputs are identical.

Functionality

This works only if the input is of length 1 (so a single number, not a vector). This assumption seems to be met. I determined this by running the below snippet from the vignettes with if(length(x) > 1) print(x) added inside these functions. Looking for input on whether this is a thorough enough test or alternatively if there may be vector inputs to dexp1 and dexp2 when they get called from within functions other than ctmm.fit.

    M.IID <- ctmm.fit(Cilla,m.iid)
    M.OU <- ctmm.fit(Cilla,m.ou)
    M.OUF <- ctmm.fit(Cilla,m.ouf

Performance

I was thinking to take this a first piece to translate to CPP, but based on benchmarks, this version is faster than what I was able to do in cpp. Compared the current version (orig_dexp2), version in this PR (if_else_dexp2), and a cpp version (cpp_dexp2) with this script:

orig_dexp2 <- function(x,Exp=exp(-x)) { ifelse(Exp<0.7071068,1-Exp^2,2*Exp*sinh(x)) }
if_else_dexp2 <- function(x,Exp=exp(-x)) { if(Exp<0.7071068) return(1-Exp^2) else return(2*Exp*sinh(x))}

Rcpp::cppFunction('
double cpp_dexp2(double x, NumericVector Exp_in = NumericVector::create()){
  double Exp = Exp_in.size() == 0 ? exp(-x) : Exp_in[0];
  if(Exp<0.7071068){
      return 1 - pow(Exp, 2.0);
  } else {
      return 2.0 * Exp * sinh(x);
  }
}')

run_over_inputs <- function(fn, precompute_Exp = FALSE){
    input <- c(100, 10, 1, .1, .01, 0, -.01, -.1, -10, -100)
    if(precompute_Exp){
        return(sapply(input, function(x) return(fn(x, exp(-x)))))
    } else {
        return(sapply(input, fn))
    }
}

result <- microbenchmark::microbenchmark(
    run_over_inputs(orig_dexp2),
    run_over_inputs(if_else_dexp2),
    run_over_inputs(cpp_dexp2),
    check = 'identical'
)

print(result)

result <- microbenchmark::microbenchmark(
    run_over_inputs(orig_dexp2, TRUE),
    run_over_inputs(if_else_dexp2, TRUE),
    run_over_inputs(cpp_dexp2, TRUE),
    check = 'identical'
)

print(result)

Sample Results

Unit: microseconds
                           expr    min      lq     mean  median     uq     max
    run_over_inputs(orig_dexp2) 43.680 44.7145 47.45272 45.9695 47.759 124.034
 run_over_inputs(if_else_dexp2) 25.152 26.0370 26.68957 26.3665 26.867  30.277
     run_over_inputs(cpp_dexp2) 45.460 47.5845 50.07772 49.2625 50.967  76.968
 neval cld
   100 a  
   100  b 
   100   c
Unit: microseconds
                                 expr    min      lq     mean  median      uq
    run_over_inputs(orig_dexp2, TRUE) 46.782 47.9780 49.19511 48.6045 49.3140
 run_over_inputs(if_else_dexp2, TRUE) 28.847 29.4125 30.82980 29.8760 30.5215
     run_over_inputs(cpp_dexp2, TRUE) 50.126 51.7805 53.68817 53.2860 54.9645
    max neval cld
 68.066   100 a  
 88.777   100  b 
 66.803   100   c

Simple if else contained in this PR seems to be the winner.

I don't expect this to have a huge impact on the overall performance of the package. Using it more as a proof of concept in adding tests and benchmarking.

[chfleming]

I'm finishing up the code I told you about where langevin() will now only be called for unique time-lags instead of being lazily recycled, so this function's call should become much less frequent on GPS data. I also made the time-lag sorting happen at a higher level before calling the Kalman filter over and over within the optimizer - this was the annoying part.

I will check if those functions are called on vectors in a couple of days.

[katrinabrock]

I added a speedup for sinc/sinch as well, replacing Vectorize({if...else}) with ifelse. This is a 5x speedup of these functions based on my tests (below), and sinch is a strong contributor to overall slowness based on my profiling. If these functions are also taking single values, not vectors, we could speed up further by replacing with if...else.. without vectorize or ifelse. I didn't go this far because I think these functions are used with longer vectors at some point, so retained all existing functionality.

test result

Unit: microseconds
                                 expr     min       lq     mean   median
   run_over_inputs(orig_sinch, input) 391.899 718.2425 941.4729 951.4340
 run_over_inputs(ifelse_sinch, input)  56.143  91.6640 142.1830 142.8185
       uq      max neval cld
 1059.483 7372.384   100  a 
  181.906  437.511   100   b
Unit: microseconds
                                                   expr     min       lq
   run_over_inputs(orig_sinch, input, precomputed_sinh) 671.441 817.0655
 run_over_inputs(ifelse_sinch, input, precomputed_sinh)  99.245 118.8870
     mean  median      uq      max neval cld
 853.2558 853.534 886.392 1270.305   100  a 
 138.8795 139.271 154.753  193.450   100   b

test script

```{r} orig_sinch <- Vectorize( function(x, SINH=sinh(x)) { if(x==0) { return(1) } else { return(SINH/x) } } ) ifelse_sinch <- function(x, SINH=sinh(x)) { return(ifelse(x==0, 1, SINH/x)) } input <- c(100, 10, 1, .1, .01, 0, -.01, -.1, -10, -100) precomputed_sinh <- sinh(input) run_over_inputs <- function(fn, input, precomputed = NA){ if(all(is.na(precomputed))){ return(list( indiv = sapply(input, fn), vect = fn(input) )) } else { return(list( indiv = sapply(1:length(input), function(idx) return(fn(input[idx], precomputed[idx]))), vect = fn(input, precomputed) )) } } result <- microbenchmark::microbenchmark( run_over_inputs(orig_sinch, input), run_over_inputs(ifelse_sinch, input), check = 'identical' ) testthat::expect_identical(run_over_inputs(orig_sinch, input), run_over_inputs(ifelse_sinch, input)) print(result) result <- microbenchmark::microbenchmark( run_over_inputs(orig_sinch, input, precomputed_sinh), run_over_inputs(ifelse_sinch, input, precomputed_sinh), check = 'identical' ) print(result) ```

[chfleming]

Thanks @katrinabrock