Merge pull request #42 from lochbika/dev

Performance optimizations (redux)

R/buildup_index.r

@@ -22,14 +22,17 @@
 
 buildup_index <- function(dmc, dc) {
   # Eq. 27a
-  bui1 <- ifelse(dmc == 0 & dc == 0, 0, 0.8 * dc * dmc / (dmc + 0.4 * dc))
+  bui1 <- 0.8 * dc * dmc / (dmc + 0.4 * dc)
+  bui1[dmc == 0 & dc == 0] <- 0
   # Eq. 27b - next 3 lines
-  p <- ifelse(dmc == 0, 0, (dmc - bui1) / dmc)
+  p <- rep(0, length(dmc))
+  p[dmc != 0] <- (dmc[dmc != 0] - bui1[dmc != 0]) / dmc[dmc != 0]
   cc <- 0.92 + ((0.0114 * dmc)^1.7)
   bui0 <- dmc - cc * p
   # Constraints
-  bui0 <- ifelse(bui0 < 0, 0, bui0)
-  bui1 <- ifelse(bui1 < dmc, bui0, bui1)
+  bui0[bui0 < 0] <- 0
+  bui1[bui1 < dmc] <- bui0[bui1 < dmc]
+  bui1[bui1 >= dmc] <- bui1[bui1 >= dmc]
   return(bui1)
 }
 

R/drought_code.r

@@ -36,30 +36,31 @@ drought_code <- function(
   fl02 <- c(6.4, 5, 2.4, 0.4, -1.6, -1.6, -1.6, -1.6, -1.6, 0.9, 3.8, 5.8)
   # Near the equator, we just use 1.4 for all months.
   # Constrain temperature
-  temp <- ifelse(temp < (-2.8), -2.8, temp)
+  temp[temp < (-2.8)] <- -2.8
 
   # Eq. 22 - Potential Evapotranspiration
   pe <- (0.36 * (temp + 2.8) + fl01[mon]) / 2
   # Daylength factor adjustment by latitude for Potential Evapotranspiration
   if (lat.adjust) {
-    pe <- ifelse(lat <= -20, (0.36 * (temp + 2.8) + fl02[mon]) / 2, pe)
-    pe <- ifelse(lat > -20 & lat <= 20, (0.36 * (temp + 2.8) + 1.4) / 2, pe)
+    pe[lat <= -20] <- (0.36 * (temp[lat <= -20] + 2.8) + fl02[mon]) / 2
+    pe[lat > -20 & lat <= 20] <- (0.36 * (temp[lat > -20 & lat <= 20] + 2.8) + 1.4) / 2
   }
   # Cap potential evapotranspiration at 0 for negative winter DC values
-  pe <- ifelse(pe < 0, 0, pe)
+  pe[pe < 0] <- 0
   ra <- prec
   # Eq. 18 - Effective Rainfall
   rw <- 0.83 * ra - 1.27
   # Eq. 19
   smi <- 800 * exp(-1 * dc_yda / 400)
   # Alteration to Eq. 21
   dr0 <- dc_yda - 400 * log(1 + 3.937 * rw / smi)
-  dr0 <- ifelse(dr0 < 0, 0, dr0)
+  dr0[dr0 < 0] <- 0
   # if precip is less than 2.8 then use yesterday's DC
-  dr <- ifelse(prec <= 2.8, dc_yda, dr0)
+  dr <- dr0
+  dr[prec <= 2.8] <- dc_yda[prec <= 2.8]
   # Alteration to Eq. 23
   dc1 <- dr + pe
-  dc1 <- ifelse(dc1 < 0, 0, dc1)
+  dc1[dc1 < 0] <- 0
   return(dc1)
 }
 

R/duff_moisture_code.r

@@ -47,62 +47,48 @@ duff_moisture_code <- function(
   ell04 <- c(11.5, 10.5, 9.2, 7.9, 6.8, 6.2, 6.5, 7.4, 8.7, 10, 11.2, 11.8)
   # For latitude near the equator, we simple use a factor of 9 for all months
   # constrain low end of temperature
-  temp <- ifelse(temp < (-1.1), -1.1, temp)
+  temp[temp < (-1.1)] <- -1.1
   # Eq. 16 - The log drying rate
   rk <- 1.894 * (temp + 1.1) * (100 - rh) * ell01[mon] * 1e-04
   # Adjust the day length and thus the drying r, based on latitude and month
   if (lat.adjust) {
-    rk <- ifelse(
-      lat <= 30 & lat > 10,
-      1.894 * (temp + 1.1) * (100 - rh) * ell02[mon] * 1e-04,
-      rk
-    )
-    rk <- ifelse(
-      lat <= -10 & lat > -30,
-      1.894 * (temp + 1.1) * (100 - rh) * ell03[mon] * 1e-04,
-      rk
-    )
-    rk <- ifelse(
-      lat <= -30 & lat >= -90,
-      1.894 * (temp + 1.1) * (100 - rh) * ell04[mon] * 1e-04,
-      rk
-    )
-    rk <- ifelse(
-      lat <= 10 & lat > -10,
-      1.894 * (temp + 1.1) * (100 - rh) * 9 * 1e-04,
-      rk
-    )
+    lat.sel <- lat <= 30 & lat > 10
+    if(any(lat.sel))rk[lat.sel] <- 1.894 * (temp[lat.sel] + 1.1) * (100 - rh[lat.sel]) * ell02[mon] * 1e-04
+    lat.sel <- lat <= -10 & lat > -30
+    if(any(lat.sel))rk[lat.sel] <- 1.894 * (temp[lat.sel] + 1.1) * (100 - rh[lat.sel]) * ell03[mon] * 1e-04
+    lat.sel <- lat <= -30 & lat >= -90
+    if(any(lat.sel))rk[lat.sel] <- 1.894 * (temp[lat.sel] + 1.1) * (100 - rh[lat.sel]) * ell04[mon] * 1e-04
+    lat.sel <- lat <= 10 & lat > -10
+    if(any(lat.sel))rk[lat.sel] <- 1.894 * (temp[lat.sel] + 1.1) * (100 - rh[lat.sel]) * 9 * 1e-04
   }
   # Constrain P
-  pr <- ifelse(
-    prec <= 1.5,
-    dmc_yda,
-    {
-      ra <- prec
-      # Eq. 11 - Net rain amount
-      rw <- 0.92 * ra - 1.27
-      # Alteration to Eq. 12 to calculate more accurately
-      wmi <- 20 + 280 / exp(0.023 * dmc_yda)
-      # Eqs. 13a, 13b, 13c
-      b <- ifelse(
-        dmc_yda <= 33,
-        100 / (0.5 + 0.3 * dmc_yda),
-        ifelse(
-          dmc_yda <= 65,
-          14 - 1.3 * log(dmc_yda),
-          6.2 * log(dmc_yda) - 17.2
-        )
-      )
-      # Eq. 14 - Moisture content after rain
-      wmr <- wmi + 1000 * rw / (48.77 + b * rw)
-      # Alteration to Eq. 15 to calculate more accurately
-      43.43 * (5.6348 - log(wmr - 20))
-    }
-  )
-  pr <- ifelse(pr < 0, 0, pr)
+  # if leq 1.5
+  pr <- prec
+  prec.sel <- prec <= 1.5
+  if(any(prec.sel))pr[prec.sel] <- dmc_yda[prec.sel]
+  # else
+  if(any(!prec.sel)){
+    b <- dmc_yda[!prec.sel]
+    dmc_yda.sel <- dmc_yda[!prec.sel]
+    b[dmc_yda.sel <= 33] <- 100 / (0.5 + 0.3 * dmc_yda.sel[dmc_yda.sel <= 33])
+    b[dmc_yda.sel > 33 & dmc_yda.sel <= 65] <- 14 - 1.3 * log(dmc_yda.sel[dmc_yda.sel > 33 & dmc_yda.sel <= 65])
+    b[dmc_yda.sel > 65] <- 6.2 * log(dmc_yda.sel[dmc_yda.sel > 65]) - 17.2
+
+    ra <- prec[!prec.sel]
+    # Eq. 11 - Net rain amount
+    rw <- 0.92 * ra - 1.27
+    # Alteration to Eq. 12 to calculate more accurately
+    wmi <- 20 + 280 / exp(0.023 * dmc_yda.sel)
+    # Eqs. 13a, 13b, 13c
+    # Eq. 14 - Moisture content after rain
+    wmr <- wmi + 1000 * rw / (48.77 + b * rw)
+    # Alteration to Eq. 15 to calculate more accurately
+    pr[!prec.sel] <- 43.43 * (5.6348 - log(wmr - 20))
+  }
+  pr[pr < 0] <- 0
   # Calculate final P (DMC)
   dmc1 <- pr + rk
-  dmc1 <- ifelse(dmc1 < 0, 0, dmc1)
+  dmc1[dmc1 < 0] <- 0
   return(dmc1)
 }
 

R/fine_fuel_moisture_code.r

@@ -32,54 +32,47 @@ fine_fuel_moisture_code <- function(ffmc_yda, temp, rh, ws, prec) {
   wmo <- FFMC_COEFFICIENT * (101 - ffmc_yda) / (59.5 + ffmc_yda)
   # Eq. 2 Rain reduction to allow for loss in
   #  overhead canopy
-  ra <- ifelse(prec > 0.5, prec - 0.5, prec)
+  ra <- prec
+  ra[prec > 0.5] <- prec[prec > 0.5] - 0.5
   # Eqs. 3a & 3b
-  wmo <- ifelse(
-    prec > 0.5,
-    ifelse(
-      wmo > 150,
-      (wmo + 0.0015 * (wmo - 150) * (wmo - 150) * sqrt(ra)
-        + 42.5 * ra * exp(-100 / (251 - wmo)) * (1 - exp(-6.93 / ra))),
-      wmo + 42.5 * ra * exp(-100 / (251 - wmo)) * (1 - exp(-6.93 / ra))
-    ),
-    wmo
-  )
+  wmo.sel <- prec > 0.5 & wmo > 150 
+  wmo[wmo.sel] <- (wmo[wmo.sel] + 0.0015 * (wmo[wmo.sel] - 150) * (wmo[wmo.sel] - 150) * sqrt(ra[wmo.sel]) + 
+                     42.5 * ra[wmo.sel] * exp(-100 / (251 - wmo[wmo.sel])) * (1 - exp(-6.93 / ra[wmo.sel])))
+  wmo.sel <- prec > 0.5 & wmo <= 150
+  wmo[wmo.sel] <- wmo[wmo.sel] + 42.5 * ra[wmo.sel] * exp(-100 / (251 - wmo[wmo.sel])) * (1 - exp(-6.93 / ra[wmo.sel]))
   # The real moisture content of pine litter ranges up to about 250 percent,
   # so we cap it at 250
-  wmo <- ifelse(wmo > 250, 250, wmo)
+  wmo[wmo > 250] <- 250
   # Eq. 4 Equilibrium moisture content from drying
   ed <- (0.942 * (rh^0.679) + (11 * exp((rh - 100) / 10))
     + 0.18 * (21.1 - temp) * (1 - 1 / exp(rh * 0.115)))
   # Eq. 5 Equilibrium moisture content from wetting
   ew <- (0.618 * (rh^0.753) + (10 * exp((rh - 100) / 10))
     + 0.18 * (21.1 - temp) * (1 - 1 / exp(rh * 0.115)))
   # Eq. 6a (ko) Log drying rate at the normal temperature of 21.1 C
-  z <- ifelse(
-    wmo < ed & wmo < ew,
-    (0.424 * (1 - (((100 - rh) / 100)^1.7))
-      + 0.0694 * sqrt(ws) * (1 - ((100 - rh) / 100)^8)),
-    0
-  )
+  z <- rep(0.0, length(wmo))
+  z.sel <- wmo < ed & wmo < ew
+  z[z.sel] <- (0.424 * (1 - (((100 - rh[z.sel]) / 100)^1.7))
+               + 0.0694 * sqrt(ws[z.sel]) * (1 - ((100 - rh[z.sel]) / 100)^8))
   # Eq. 6b Affect of temperature on  drying rate
   x <- z * 0.581 * exp(0.0365 * temp)
   # Eq. 8
-  wm <- ifelse(wmo < ed & wmo < ew, ew - (ew - wmo) / (10^x), wmo)
+  wm <- wmo
+  wm.sel <- wmo < ed & wmo < ew
+  wm[wm.sel] <- ew[wm.sel] - (ew[wm.sel] - wmo[wm.sel]) / (10^x[wm.sel])
   # Eq. 7a (ko) Log wetting rate at the normal temperature of 21.1 C
-  z <- ifelse(
-    wmo > ed,
-    (0.424 * (1 - (rh / 100)^1.7)
-      + 0.0694 * sqrt(ws) * (1 - (rh / 100)^8)),
-    z
-  )
+  z.sel <- wmo > ed
+  z[z.sel] <- (0.424 * (1 - (rh[z.sel] / 100)^1.7) + 0.0694 * sqrt(ws[z.sel]) * (1 - (rh[z.sel] / 100)^8))
   # Eq. 7b Affect of temperature on  wetting rate
   x <- z * 0.581 * exp(0.0365 * temp)
   # Eq. 9
-  wm <- ifelse(wmo > ed, ed + (wmo - ed) / (10^x), wm)
+  wm[wmo > ed] <- ed[wmo > ed] + (wmo[wmo > ed] - ed[wmo > ed]) / (10^x[wmo > ed])
+  wm[is.na(wmo > ed)] <- NA
   # Eq. 10 Final ffmc calculation
   ffmc1 <- (59.5 * (250 - wm)) / (FFMC_COEFFICIENT + wm)
   # Constraints
-  ffmc1 <- ifelse(ffmc1 > 101, 101, ffmc1)
-  ffmc1 <- ifelse(ffmc1 < 0, 0, ffmc1)
+  ffmc1[ffmc1 > 101] <- 101
+  ffmc1[ffmc1 < 0] <- 0
   return(ffmc1)
 }
 

R/fwi.r

@@ -403,13 +403,13 @@ fwi <- function(
   # Length of weather run
   n0 <- length(temp) / n
   # Initialize variables
-  ffmc <- dmc <- dc <- isi <- bui <- fwi <- dsr <- NULL
+  ffmc <- dmc <- dc <- isi <- bui <- fwi <- dsr <- rep(NA_real_, length(temp))
   # For each day in the run
   for (i in 1:n0) {
     # k is the data for all stations by day
     k <- ((i - 1) * n + 1):(i * n)
     # constrain relative humidity
-    rh[k] <- ifelse(rh[k] >= 100, 99.9999, rh[k])
+    if(any(rh[k] >= 100))rh[k][rh[k] >= 100] <- 99.9999
     ###########################################################################
     # Fine Fuel Moisture Code (FFMC)
     ###########################################################################
@@ -452,13 +452,13 @@ fwi <- function(
     dsr1 <- 0.0272 * (fwi1^1.77)
 
     # Concatenate values
-    ffmc <- c(ffmc, ffmc1)
-    dmc <- c(dmc, dmc1)
-    dc <- c(dc, dc1)
-    isi <- c(isi, isi1)
-    bui <- c(bui, bui1)
-    fwi <- c(fwi, fwi1)
-    dsr <- c(dsr, dsr1)
+    ffmc[k] <- ffmc1
+    dmc[k] <- dmc1
+    dc[k] <- dc1
+    isi[k] <- isi1
+    bui[k] <- bui1
+    fwi[k] <- fwi1
+    dsr[k] <- dsr1
     ffmc_yda <- ffmc1
     dmc_yda <- dmc1
     dc_yda <- dc1

R/initial_spread_index.r

@@ -31,11 +31,15 @@ initial_spread_index <- function(
   # Eq. 24 - Wind Effect
   # the ifelse, also takes care of the ISI modification for the fbp functions
   # This modification is Equation 53a in FCFDG (1992)
-  fW <- ifelse(
-    ws >= 40 & fbpMod == TRUE,
-    12 * (1 - exp(-0.0818 * (ws - 28))),
-    exp(0.05039 * ws)
-  )
+  if(fbpMod){
+    fW <- rep(NA_real_, length(ws))
+    fW.sel <- ws >= 40 & !is.na(ws)
+    fW[fW.sel] <- 12 * (1 - exp(-0.0818 * (ws[fW.sel] - 28)))
+    fW[!fW.sel] <- exp(0.05039 * ws[!fW.sel])
+    fW[is.na(fW.sel)] <- NA
+  }else{
+    fW <- exp(0.05039 * ws)
+  }
   # Eq. 25 - Fine Fuel Moisture
   fF <- 91.9 * exp(-0.1386 * fm) * (1 + (fm^5.31) / 49300000)
   # Eq. 26 - Spread Index Equation

tests/testthat/test_fwi.r

@@ -76,11 +76,8 @@ test_that("fwi_10", {
 test_that("fwi_11", {
   fct_test_fwi("fwi_11",
            function(test_fwi) {
-             expect_warning({
-               expect_warning(
+             expect_warning(
                  {actual <- fwi(test_fwi, batch = FALSE)},
-                 "*NaNs produced*")
-               },
                "Same initial data were used for multiple weather stations")
            return(actual)
            })