Skip to content

Commit

Permalink
add Cpp function to CMD sum calculation for better speedz
Browse files Browse the repository at this point in the history
  • Loading branch information
kdaust committed Nov 15, 2024
1 parent 9cd04ef commit 610bde0
Showing 1 changed file with 39 additions and 12 deletions.
51 changes: 39 additions & 12 deletions data_processing/Extreme_Vars.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -48,23 +48,50 @@ bio_stats <- dcast(temp, id ~ variable + Quantile)

## CMD Variables

```{r}
calc_cmd <- function(vals){
res <- numeric(length(vals))
res[1] <- vals[1]
for(i in 2:length(vals)){
prev <- res[i-1]
```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;
using namespace std;
// [[Rcpp::export]]
NumericVector cmd_sum(NumericVector vals, double cutoff){
NumericVector res = NumericVector(vals.size());
res[0] = vals[0];
double prev, tmp;
for(int i = 1; i < vals.size(); i++){
prev = res[i-1];
if(vals[i] < 0){
if(prev > 500) prev <- 500
if(prev > cutoff){
prev = cutoff;
}
}
tmp = prev + vals[i];
if(tmp < 0){
tmp = 0;
}
tmp <- prev + vals[i]
if(tmp < 0) tmp <- 0
res[i] <- tmp
res[i] = tmp;
}
return(res)
return(res);
}
```

```{r}
# calc_cmd <- function(vals){
# res <- numeric(length(vals))
# res[1] <- vals[1]
# for(i in 2:length(vals)){
# prev <- res[i-1]
# if(vals[i] < 0){
# if(prev > 500) prev <- 500
# }
# tmp <- prev + vals[i]
# if(tmp < 0) tmp <- 0
# res[i] <- tmp
# }
# return(res)
# }
```

```{r}
res <- downscale(pnts, which_refmap = "refmap_climr",
obs_years = 1961:1990,
Expand All @@ -83,7 +110,7 @@ cmd_res <- res2[Var == "CMD",.(id,PERIOD,Var,Month,value)]
setorder(cmd_res, id, PERIOD, Month)
cmd_res[precip_recharge, precip_sum := precip_sum, on = c("id","PERIOD","Month")]
cmd_res[!is.na(precip_sum),value := value - precip_sum]
cmd_res[,cmd_sum := calc_cmd(value), by = .(id)]
cmd_res[,cmd_sum := cmd_sum(value, cutoff = 500), by = .(id)]
cmd_res[,Date := as.Date(paste0(PERIOD,"-",Month,"-01"))]
Expand Down

0 comments on commit 610bde0

Please sign in to comment.