japilo / colorednoise

@@ -85,11 +85,14 @@
Loading
85 85
//' @param covMatrix A valid covariance matrix. The number of rows/columns must match the length of the mu, sigma, and phi vectors.
86 86
//' @return A matrix with as many rows as timesteps and as many columns as mu/sigma/phi values.
87 87
//' @examples
88 -
//' cov <- matrix(c(0.037, 0.044, -0.048, 0.044, 0.247, -0.008, -0.047, -0.008, 0.074), nrow = 3)
88 +
//' cov <- matrix(c(1, 0.53, 0.73, 0.53, 1, 0.44, 0.73, 0.44, 1), nrow = 3)
89 89
//' test <- colored_multi_rnorm(100, c(0, 3, 5), c(1, 0.5, 1), c(0.5, -0.3, 0), cov)
90 90
//' var(test)
91 -
//' library(dplyr)
92 -
//' test %>% as.data.frame() %>% summarize_all(.funs = c("mean", "sd", "autocorrelation"))
91 +
//' library(data.table)
92 +
//' as.data.table(test)[, .(V1_mean = mean(V1), V2_mean = mean(V2), V3_mean = mean(V3),
93 +
//' V1_sd = sd(V1), V2_sd = sd(V2), V3_sd = sd(V3),
94 +
//' V1_autocorrelation = autocorrelation(V1), V2_autocorrelation = autocorrelation(V2),
95 +
//' V3_autocorrelation = autocorrelation(V3))]
93 96
//' @export
94 97
// [[Rcpp::export]]
95 98
Rcpp::NumericMatrix colored_multi_rnorm(int timesteps, Rcpp::NumericVector mean, Rcpp::NumericVector sd, Rcpp::NumericVector phi, Rcpp::NumericMatrix covMatrix) {

@@ -57,8 +57,8 @@
Loading
57 57
        }
58 58
    }
59 59
    # Unnests the list and adds estimates of survival and fertility
60 -
    sims <- labeled_sims %>% flatten() %>% map(~mutate(., est_surv = survivors/population,
61 -
        est_fecund = newborns/survivors))
60 +
    sims <- labeled_sims %>% flatten() %>% map(setDT) %>%
61 +
      map(~.[, `:=`(est_surv = survivors/population, est_fecund = newborns/survivors)])
62 62
    return(sims)
63 63
}
64 64
@@ -117,10 +117,10 @@
Loading
117 117
    colNames = NULL, matrixStructure = NULL, repeatElements = NULL,
118 118
    survivalOverflow = "scale") {
119 119
    stages <- length(initialPop)
120 -
    # Regularize all valid data inputs to the same format
120 +
    # Regularize all valid data input to the same format
121 121
    if (is.data.frame(data) == T) {
122 122
        if (is.null(colNames) == F) {
123 -
            data <- data[, colNames] %>% rename(!(!(!colNames)))
123 +
            data <- data %>% setDT() %>% setnames(old = colNames, new = names(colNames))
124 124
        }
125 125
        if (all(names(data) == c("mean", "sd", "autocorrelation")) ==
126 126
            F) {
@@ -135,7 +135,7 @@
Loading
135 135
        if (all(data$sd > 0) == F) {
136 136
            stop("Invalid values in SD column")
137 137
        }
138 -
        dat <- data %>% as_tibble()
138 +
        dat <- setDT(data)
139 139
    } else if (is.list(data) == T) {
140 140
        if (length(data) > 3) {
141 141
            stop("List data should only have 3 elements")
@@ -149,7 +149,7 @@
Loading
149 149
        if (length(unique(map_int(data, length)))>1) {
150 150
          stop("Matrices are not equal dimensions")
151 151
        }
152 -
        dat <- tibble(mean = as.vector(t(data[[1]])), sd = as.vector(t(data[[2]])),
152 +
        dat <- data.table(mean = as.vector(t(data[[1]])), sd = as.vector(t(data[[2]])),
153 153
            autocorrelation = as.vector(t(data[[3]])))
154 154
    } else {
155 155
        stop("Invalid data type. Must be a list of three matrices or a data frame with three columns.")
@@ -166,31 +166,35 @@
Loading
166 166
    if (is.null(repeatElements) == T) {
167 167
      repeatElements <- matrix(seq(1:stages^2), ncol = stages, byrow = T)
168 168
    }
169 -
    dat <- dat %>% mutate(dist = dists, zero = mean==0&sd==0,
170 -
                          ref = as.vector(t(repeatElements)))
169 +
    dat <- dat[, `:=`(dist = dists, zero = mean == 0 & sd == 0, ref = as.vector(t(repeatElements)))]
171 170
    repeats <- repeatElements == matrix(seq(1:stages^2), ncol = stages, byrow = T)
172 171
    # Create version of data that can be used to generate colored noise
173 -
    inputs <- dat %>% slice(which(t(repeats))) %>% filter(zero == F) %>%
174 -
      rowwise() %>% mutate(
175 -
        mean.trans = ifelse(mean == 0, 0, invoke(dist, list(mean))),
176 -
        sd.trans = ifelse(sd == 0, 0, variancefix(mean, sd, dist))
177 -
      )
172 +
    unique <- dat[which(t(repeats))][zero == F]
173 +
    unique$mean.trans <- map_dbl(c(1:nrow(unique)), function(x) {
174 +
      if(dat[x,]$mean == 0) {0} else {
175 +
        invoke(dat[x,]$dist, list(dat[x,]$mean))
176 +
        }
177 +
    })
178 +
    unique$sd.trans <- map_dbl(c(1:nrow(unique)), function(x) {
179 +
      if(dat[x,]$sd == 0) {0} else {
180 +
        stdev_transform(dat[x,]$mean, dat[x,]$sd, dat[x,]$dist)
181 +
      }
182 +
    })
178 183
    if (is.null(covMatrix) == T) {
179 -
      covMatrix <- cor2cov(inputs$sd, diag(nrow(inputs)))
184 +
      covMatrix <- cor2cov(unique$sd, diag(nrow(unique)))
180 185
    }
181 186
    # Create colored noise, discard if invalid matrix
182 187
    if(survivalOverflow == "redraw") {
183 188
      repeat {
184 -
        inputs$noise <- colored_multi_rnorm(timesteps, inputs$mean.trans, inputs$sd.trans,
185 -
                                            inputs$autocorrelation, covMatrix) %>% split(col(.))
186 -
        result <- left_join(dat, inputs, by = c("mean", "sd", "autocorrelation", "dist", "zero", "ref"))
189 +
        unique$noise <- colored_multi_rnorm(timesteps, unique$mean.trans, unique$sd.trans,
190 +
                                            unique$autocorrelation, covMatrix) %>% split(col(.))
191 +
        result <- dat[unique, on = .(mean, sd, autocorrelation, dist, zero, ref), allow.cartesian = TRUE]
187 192
        result$noise[dat$zero==T] <- rep(list(rep.int(0, timesteps)), sum(dat$zero==T))
188 193
        # checking for >1 probability
189 -
        result <- result %>% rowwise() %>% mutate(
190 -
          natural.noise = ifelse(zero == T, list(noise),
191 -
                                 ifelse(dist == "log", list(exp(noise)),
192 -
                                        list(plogis(noise))))
193 -
        )
194 +
        result$natural.noise <- map(c(1:nrow(result)), function(x) {
195 +
          if (result[x,]$zero) {result[x,]$noise} else if (result[x,]$dist=="log") {
196 +
            exp(result[x,]$noise[[1]])} else {plogis(result[x,]$noise[[1]])}
197 +
          })
194 198
        matrices <- map(1:timesteps, function(x) {
195 199
          matrix(map_dbl(result$natural.noise, x), byrow = T, ncol = stages)
196 200
        })
@@ -201,15 +205,15 @@
Loading
201 205
        }
202 206
      }
203 207
    } else if(survivalOverflow == "scale") {
204 -
      inputs$noise <- colored_multi_rnorm(timesteps, inputs$mean.trans, inputs$sd.trans,
205 -
                                          inputs$autocorrelation, covMatrix) %>% split(col(.))
206 -
      result <- left_join(dat, inputs, by = c("mean", "sd", "autocorrelation", "dist", "zero", "ref"))
208 +
      unique$noise <- colored_multi_rnorm(timesteps, unique$mean.trans, unique$sd.trans,
209 +
                                          unique$autocorrelation, covMatrix) %>% split(col(.))
210 +
      result <- dat[unique, on = .(mean, sd, autocorrelation, dist, zero, ref), allow.cartesian = TRUE]
207 211
      result$noise[dat$zero==T] <- rep(list(rep.int(0, timesteps)), sum(dat$zero==T))
208 212
      # checking for >1 probability
209 -
      result <- result %>% rowwise() %>% mutate(
210 -
        natural.noise = ifelse(zero == T, list(noise),
211 -
                               ifelse(dist == "log", list(exp(noise)),
212 -
                                      list(plogis(noise)))))
213 +
      result$natural.noise <- map(c(1:nrow(result)), function(x) {
214 +
        if (result[x,]$zero) {result[x,]$noise} else if (result[x,]$dist=="log") {
215 +
          exp(result[x,]$noise[[1]])} else {plogis(result[x,]$noise[[1]])}
216 +
      })
213 217
      matrices <- map(1:timesteps, function(x) {
214 218
        matrix(map_dbl(result$natural.noise, x), byrow = T, ncol = stages)
215 219
      })
@@ -226,7 +230,9 @@
Loading
226 230
      }
227 231
    } else {stop("survivalOverflow must be set to 'redraw' or 'scale'")}
228 232
    pop <- projection(initialPop, matrices)
229 -
    pop %>% map(as_tibble, .name_repair = ~ c(paste0("stage", 1:stages))) %>%
230 -
      bind_rows() %>% group_by(timestep = row_number()) %>% nest(data = -timestep) %>%
231 -
      mutate(total = map_dbl(data, sum)) %>% unnest(data)
233 +
    # CONVERT TO DATA.TABLE
234 +
    output <- pop %>% map(as.data.table) %>% rbindlist() %>%
235 +
      .[, `:=`(total = rowSums(.SD), timestep = seq_len(.N))] %>% setcolorder("timestep")
236 +
    names(output) <- c("timestep", map_chr(1:stages, ~paste0("stage", .)), "total")
237 +
    return(output)
232 238
}

@@ -2,12 +2,28 @@
Loading
2 2
#include "noise.h"
3 3
// [[Rcpp::depends(RcppArmadillo)]]
4 4
5 -
// Variance Fix
6 -
//
7 -
// This function changes the variance so that once it is back-transformed
8 -
// from the logit scale, the original variance is recovered.
5 +
//' Translate Standard Deviation from the Natural Scale to the Log or Logit Scale
6 +
//'
7 +
//' This function changes a given standard deviation so that when a vector of samples is drawn from the given distribution,
8 +
//' the original standard deviation will be recovered once it is back-transformed from the log or logit scale. In effect,
9 +
//' the function "translates" a standard deviation from the natural scale to the log or logit scale for the purposes of
10 +
//' random draws from a probability distribution.
11 +
//' @param mu The mean of the distribution on the natural scale.
12 +
//' @param sigma The standard devation of the distribution on the natural scale.
13 +
//' @param dist The distribution to which the standard deviation should be transformed.
14 +
//' @return The standard deviation translated to the log or logit scale.
15 +
//' @examples
16 +
//' mean <- 10
17 +
//' stdev <- 2
18 +
//' mean_trans <- log(mean)
19 +
//' stdev_trans <- stdev_transform(mean, stdev, "log")
20 +
//' draws <- rnorm(50, mean_trans, stdev_trans)
21 +
//' natural_scale <- exp(draws)
22 +
//' mean(draws)
23 +
//' sd(draws)
24 +
//' @export
9 25
// [[Rcpp::export]]
10 -
double variancefix(double mu, double sigma, std::string dist){
26 +
double stdev_transform(double mu, double sigma, std::string dist){
11 27
  if (sigma == 0) {
12 28
    return 0;
13 29
  }
@@ -68,11 +84,11 @@
Loading
68 84
Rcpp::DataFrame unstructured_pop(int start, int timesteps, double survPhi, double fecundPhi, double survMean, double survSd, double fecundMean, double fecundSd) {
69 85
  // These lines generate the temporally autocorrelated random numbers
70 86
  double survmu = R::qlogis(survMean, 0, 1, true, false);
71 -
  double survsigma = variancefix(survMean, survSd, "qlogis");
87 +
  double survsigma = stdev_transform(survMean, survSd, "qlogis");
72 88
  Rcpp::NumericVector survnoise = colored_noise(timesteps, survmu, survsigma, survPhi);
73 89
  Rcpp::NumericVector St = plogis(survnoise);
74 90
  double fecundmu = log(fecundMean);
75 -
  double fecundsigma = variancefix(fecundMean, fecundSd, "log");
91 +
  double fecundsigma = stdev_transform(fecundMean, fecundSd, "log");
76 92
  Rcpp::NumericVector fecundnoise = colored_noise(timesteps, fecundmu, fecundsigma, fecundPhi);
77 93
  Rcpp::NumericVector Ft = exp(fecundnoise);
78 94
  // This loop kills some individuals according to St probabilities, creates
Files Coverage
R 50.47%
src 77.06%
Project Totals (4 files) 63.89%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading