Nothing
inst/optimal-control.R
In mosqcontrol: Mosquito Control Resource Optimization
# install.packages("pracma")
# install.packages("matlib")
# install.packages("NlcOptim")
library("pracma")
library("matlib")
# sample data
# mosquito count from data
y_in <- c(
15, 40, 45, 88, 99, 145, 111, 132, 177, 97, 94, 145, 123, 111,
125, 115, 155, 160, 143, 132, 126, 125, 105, 98, 87, 54, 55, 8
)
# day of year in which measurments making up yin occur, Jan 1st = day 1
t_in_user <- c(
93, 100, 107, 114, 121, 128, 135, 142, 149, 163, 170, 177,
184, 191, 198, 205, 212, 219, 226, 233, 240, 247, 254, 261,
267, 274, 281, 288
)
if (length(y_in) != length(t_in_user)) {
stop("input data mismatch")
}
# User inputs (aside from trap data inputs)
# mu = natural mosquito death rate
mu <- 1 / 14 # program should set mu = 1 / tau_mosq, where tau_mosq is the mosquito lifetime input by user in units of days
# tau_mosq should be bounded between 1 and 30
# m = number of mosquito lifetimes for population decay between seasons..
m <- 3 # this will be input by the user, and must be an integer greater than 2, upperbound 100, default value 3
# N_lam = max fourier mode order to calculate
N_lam <- 25
# user input between 1 and 200, default value 25 (I referred to this as jmax over Skype - JD)
# setting up time and trap count vectors to fit the model to
# first round user time inputs to integer values
t_in <- round(t_in_user)
# defining N0 = number of input data points (not a user input)
N0 <- length(y_in)
# defining vector of time differences to be used later (not a user input)
delta_t_in <- numeric(N0 - 1)
for (i in 1:(N0 - 1)) {
delta_t_in[i] <- t_in[i + 1] - t_in[i]
if (delta_t_in[i] < 0) {
stop("non-increasing measurement times")
}
}
# N_pts = number of points in extended data set, N_pts = N0 + 3 always (not a user input)
N_pts <- N0 + 3
# defining extended trap count and time data vectors
t_dat <- numeric(N_pts)
if (t_in[1] - m / mu > 0) {
y_dat <- c(0, 0, y_in, 0)
t_dat[1] <- 0
t_dat[2] <- (m - 1) / mu
for (i in 3:(N_pts - 1)) {
t_dat[i] <- t_in[i - 2] - t_in[1] + m / mu
}
t_dat[N_pts] <- t_dat[N_pts - 1] + (m) / mu
t_dat <- matrix(t_dat, ncol = 1)
} else {
y_dat <- c(0, y_in, 0, 0)
t_dat[1] <- 0
for (i in 2:(N_pts - 2)) {
t_dat[i] <- t_in[i - 1]
}
t_dat[N_pts - 1] <- t_dat[N_pts - 2] + (2 * (m - 1)) / mu - t_dat[2]
t_dat[N_pts] <- t_dat[N_pts - 1] + 1 / mu
}
# extended vector of time differences
delta_t_dat <- numeric(N_pts)
for (i in 1:(N_pts - 1)) {
delta_t_dat[i] <- t_dat[i + 1] - t_dat[i]
}
delta_t_dat[N_pts] <- 0
delta_t_dat <- matrix(delta_t_dat, ncol = 1)
# tau = length of 'season' in days
tau <- t_dat[N_pts]
# Defining matrices 'J' and 'M'
# J matrix and inverse
J <- matrix(0, nrow = N_pts, ncol = N_pts)
for (i in 1:(N_pts - 1)) {
J[i, i] <- 1
J[i + 1, i] <- -exp(-mu * delta_t_dat[i])
}
J[N_pts, N_pts] <- 1
J[1, N_pts] <- -exp(-mu * delta_t_dat[N_pts])
# M matrix and inverse
M <- matrix(0, nrow = N_pts, ncol = N_pts)
for (i in 2:(N_pts - 1)) {
M[i, i] <- (1 - exp(-mu * delta_t_dat[i - 1])) / (mu * delta_t_dat[i - 1]) - exp(-mu * delta_t_dat[i - 1])
M[i + 1, i] <- 1 - (1 - exp(-mu * delta_t_dat[i])) / (mu * delta_t_dat[i])
}
M[N_pts, N_pts] <- (1 - exp((-1) * mu * delta_t_dat[N_pts - 1])) / (mu * delta_t_dat[(N_pts - 1)]) - exp(-mu * delta_t_dat[(N_pts - 1)])
M[1, N_pts] <- 0
M[1, 1] <- 0
M[2, 1] <- 1 - (1 - exp(-mu * delta_t_dat[1])) / (mu * delta_t_dat[(1)])
###### weighting code code
# Weight matrix: diagonal entries to be filled with sqrt(number of counts on day) except for the artifical 0 measurments on end of data
W <- matrix(0, nrow = N_pts, ncol = N_pts)
if (t_in[(1)] - m / mu > 0) {
W[1, 1] <- 1
W[2, 2] <- 1
W[N_pts, N_pts] <- 1
for (i in 3:N_pts - 1) {
W[i, i] <- sqrt(1) # fill sqrt(number of counts) here
}
} else {
W[1, 1] <- 1
W[N_pts - 1, N_pts - 1] <- 1
W[N_pts, N_pts] <- 1
for (i in 2:N_pts - 2) {
W[i, i] <- sqrt(1) # fill sqrt(number of counts) here
}
}
# Obtaining emergence rate vector lambda_dat by contrained optimization
Aeq <- matrix(numeric(N_pts), nrow = 1)
Aeq[(1)] <- 1
Aeq[(N_pts)] <- -1
beq <- 0
bineq <- numeric(N_pts)
lb <- numeric(N_pts)
###################
objective_fun <- function(x) {
t(W %*% ((1 / mu) * (mldivide(J, (M %*% x), pinv = TRUE)) - y_dat)) %*% (W %*% ((1 / mu) * (mldivide(J, (M %*% x), pinv = TRUE)) - y_dat))
}
lambda_dat <- fmincon(numeric(N_pts), objective_fun, gr = NULL, method = "SQP", A = (-(1 / mu) * (inv(J) %*% (M))), b = bineq, Aeq = Aeq, beq = beq, lb = lb, ub = NULL)
#########
# Calculating emergence rate fourier modes lam_fourier
lam_fourier <- numeric(N_lam + 1)
for (j in 1:(N_pts - 1)) {
lam_fourier[1] <- lam_fourier[1] + (1 / tau) * (lambda_dat$par[j] + lambda_dat$par[j + 1]) * delta_t_dat[j] / 2
} # zero mode
for (k in 2:(N_lam + 1)) # non-zero postive modes (negative modes given by complex conjugates of positve modes)
{
for (j in 1:(N_pts - 1)) {
lam_fourier[(k)] <- lam_fourier[(k)] + ((lambda_dat$par[(j + 1)] - lambda_dat$par[(j)]) / delta_t_dat[(j)]) * (tau / (2 * pi * (k - 1))) * ((exp(-2 * pi * complex(real = 0, imaginary = 1) * (k - 1) * t_dat[(j + 1)] / tau) - exp(-2 * pi * complex(real = 0, imaginary = 1) * (k - 1) * t_dat[(j)] / tau)) / (2 * pi * (k - 1))) + complex(real = 0, imaginary = 1) * (lambda_dat$par[(j + 1)] * exp(-2 * pi * complex(real = 0, imaginary = 1) * (k - 1) * t_dat[(j + 1)] / tau) - lambda_dat$par[(j)] * exp(-2 * pi * complex(real = 0, imaginary = 1) * (k - 1) * t_dat[(j)] / tau)) / (2 * pi * (k - 1))
}
}
# defing functions for optimization code
# goto line 1670 or so to skip past all of these definitons
# Permutaion function: For a given list of numbers, this function outputs a matrix, where each row is a unique permutation of the list
uperm <- function(d) {
dat <- factor(d)
N <- length(dat)
n <- tabulate(dat)
ng <- length(n)
if (ng == 1) {
return(d)
}
a <- N - c(0, cumsum(n))[-(ng + 1)]
foo <- lapply(1:ng, function(i) matrix(combn(a[i], n[i]), nrow = n[i]))
out <- matrix(NA, nrow = N, ncol = prod(sapply(foo, ncol)))
xxx <- c(0, cumsum(sapply(foo, nrow)))
xxx <- cbind(xxx[-length(xxx)] + 1, xxx[-1])
miss <- matrix(1:N, ncol = 1)
for (i in seq_len(length(foo) - 1)) {
l1 <- foo[[i]]
nn <- ncol(miss)
miss <- matrix(rep(miss, ncol(l1)), nrow = nrow(miss))
k <- (rep(0:(ncol(miss) - 1), each = nrow(l1))) * nrow(miss) +
l1[, rep(1:ncol(l1), each = nn)]
out[xxx[i, 1]:xxx[i, 2], ] <- matrix(miss[k], ncol = ncol(miss))
miss <- matrix(miss[-k], ncol = ncol(miss))
}
k <- length(foo)
out[xxx[k, 1]:xxx[k, 2], ] <- miss
out <- out[rank(as.numeric(dat), ties = "first"), ]
foo <- cbind(as.vector(out), as.vector(col(out)))
out[foo] <- d
t(out)
}
# this function computes the "P" product and exponentials for a given k
P_vals <- function(z, rho, Npulse, tau, k) {
P_func <- function(x) {
rho * 1 / (log(1 / (1 - rho)) - 2 * pi * complex(real = 0, imaginary = 1) * x)
}
z_vec <- matrix(0, nrow = Npulse, ncol = 1)
for (i in 1:Npulse) {
z_vec[i, 1] <- z[i]
}
# determine number of integer vectors to permute (will be smaller for
# only two and three pulses). number of combos = N_terms
if (k == 0) {
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 7
} else {
N_terms <- 8
}
# create a matrix to hold the integer vectors to permute, and
# initialize to zero. Each column will give an integer vector
int_vecs <- matrix(0, nrow = Npulse, ncol = N_terms)
# the first 5 integer vectors have at most two non-zero elements.
# This loop replaces the first two zeros in the first 5 columns of int_vecs
# with the appropriate integers.
for (i in 1:5) {
int_vecs[1:2, i] <- c((i - 1), -(i - 1))
}
# if there are are three or more pulses, we need to include
# integer vectors with three non-zero terms
if (Npulse >= 3) {
int_vecs[1:3, 6] <- c(2, -1, -1)
int_vecs[1:3, 7] <- c(-2, 1, 1)
}
# if there are four or more pulses, then we need to include an
# integer vector with four non-zero terms
if (Npulse >= 4) {
int_vecs[1:4, 8] <- c(1, 1, -1, -1)
}
# calculates the "P" functions for each column of integer combos in
# int_vec. "P" function values are independet of integer
# permutations. pre_factors is a list which holds the "P"-function
# values for each integer combo
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
# Now calculate the exponential terms which will multiply the "P"functions in the pre_factors list
# post_factors = list of the exponential terms to multiply the
# corresponding "P" functions in the pre_factors list
post_factors <- numeric(N_terms) # initialize to zero
for (i in 1:N_terms) { # for each integer combo given by the columns of int_vecs...
# calculates the unique permutations of the ith integer combo
# int_vecs(:, i) = ith column of int_vecs matrix
# perm_list = Nperms x Npulse matrix, where Nperms is the number
# of unique permutations of the ith integer combo
# Each row of perm_list gives a unique permutaion of the ith integer combo
perm_list <- uperm(int_vecs[, i])
if (i == 1) {
Nperms <- 1
} else {
Nperms <- length(perm_list[, 1]) # number of unique permutations of the ith integer combo
}
for (m in 1:Nperms) { # for each uniqe permutation of the ith integer combo, add the corresponding exponential to the running total of the corresponding post_factor value
if (i == 1) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * (perm_list %*% z) / tau)
} else {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * (perm_list[m, ] %*% z) / tau)
}
} # perm_list[m, ] = mth row of perm_list = mth unique permutaton of ith integer combo
}
val <- 0
for (i in 1:N_terms) { # add together the (pre_factor * post_factor) for each integer combo
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 1) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 8
}
int_vecs <- matrix(0, nrow = Npulse, ncol = N_terms)
for (i in 1:5) {
int_vecs[1:2, i] <- sign(k) * c(kk + (i - 1), -(i - 1))
}
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(1, 1, -1)
int_vecs[1:3, 7] <- sign(k) * c(3, -1, -1)
int_vecs[1:3, 8] <- sign(k) * c(2, -2, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
}
else if (abs(k) == 2) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 6
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk + 2, -2)
int_vecs[1:2, 5] <- sign(k) * c(kk + 3, -3)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(2, -1, 1)
}
if (Npulse >= 4) {
int_vecs[1:4, 7] <- sign(k) * c(1, 1, 1, -1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 2) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 3) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 8
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk + 2, -2)
int_vecs[1:2, 5] <- sign(k) * c(kk + 3, -3)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(1, 1, 1)
int_vecs[1:3, 7] <- sign(k) * c(3, 1, -1)
int_vecs[1:3, 8] <- sign(k) * c(2, 2, -1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 3 && i == 6) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 4) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 7
} else {
N_terms <- 8
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(4, -1, 1)
int_vecs[1:3, 7] <- sign(k) * c(2, 1, 1)
}
if (Npulse >= 4) {
int_vecs[1:4, 8] <- sign(k) * c(1, 1, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 4 && i == 8) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else if (Npulse == 2 && i == 3) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 5) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(2, 2, 1)
int_vecs[1:3, 7] <- sign(k) * c(3, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 6) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 6
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk - 3, 3)
int_vecs[1:2, 5] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 6] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 7] <- sign(k) * c(4, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 4) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 7) {
kk <- abs(k)
N_terms <- 6
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk - 3, 3)
int_vecs[1:2, 5] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 6] <- sign(k) * c(kk + 2, -2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 8) {
kk <- abs(k)
N_terms <- 6
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk - 3, +3)
int_vecs[1:2, 6] <- sign(k) * c(kk - 4, +4)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 6) {
Nperm <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 9) {
kk <- abs(k)
N_terms <- 5
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 5] <- sign(k) * c(kk - 3, +3)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[i] * post_factors[i]
}
return(val)
} else if (abs(k) == 10) {
kk <- abs(k)
N_terms <- 4
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 11) {
kk <- abs(k)
N_terms <- 4
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
out <- val
} else {
kk <- abs(k)
N_terms <- 3
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[i] * post_factors[i]
}
return(val)
}
}
# this function computes the "Q" product and exponentials with the fourier modes for a given k
Q_vals <- function(z, rho, Npulse, tau, j, jj) {
Q_func <- function(x) {
(rho / (1 - rho)) * 1 / (log(1 / (1 - rho)) + 2 * pi * complex(real = 0, imaginary = 1) * x)
}
k <- -(j + jj)
# determine number of integer vectors to permute (will be smaller for
# only two and three pulses). number of combos = N_terms
if (k == 0) {
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 7
} else {
N_terms <- 8
}
# create a matrix to hold the integer vectors to permute, and
# initialize to zero. Each column will give an integer vector
int_vecs <- matrix(0, nrow = Npulse, ncol = N_terms)
# the first 5 integer vectors have at most two non-zero elements.
# This loop replaces the first two zeros in the first 5 columns of int_vecs
# with the appropriate integers.
for (i in 1:5) {
int_vecs[1:2, i] <- c((i - 1), -(i - 1))
}
# if there are are three or more pulses, we need to include
# integer vectors with three non-zero terms
if (Npulse >= 3) {
int_vecs[1:3, 6] <- c(2, -1, -1)
int_vecs[1:3, 7] <- c(-2, 1, 1)
}
# if there are four or more pulses, then we need to include an
# integer vector with four non-zero terms
if (Npulse >= 4) {
int_vecs[1:4, 8] <- c(1, 1, -1, -1)
}
# calculates the "P" functions for each column of integer combos in
# int_vec. "P" function values are independet of integer
# permutations. pre_factors is a list which holds the "P"-function
# values for each integer combo
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
# Now calculate the exponential terms which will multiply the "P"functions in the pre_factors list
# post_factors = list of the exponential terms to multiply the
# corresponding "P" functions in the pre_factors list
post_factors <- numeric(N_terms) # initialize to zero
for (i in 1:N_terms) { # for each integer combo given by the columns of int_vecs...
# calculates the unique permutations of the ith integer combo
# int_vecs(:, i) = ith column of int_vecs matrix
# perm_list = Nperms x Npulse matrix, where Nperms is the number
# of unique permutations of the ith integer combo
# Each row of perm_list gives a unique permutaion of the ith integer combo
perm_list <- uperm(int_vecs[, i])
if (i == 1) {
Nperms <- 1
} else {
Nperms <- length(perm_list[, 1]) # number of unique permutations of the ith integer combo
}
for (m in 1:Nperms) { # for each uniqe permutation of the ith integer combo, add the corresponding exponential to the running total of the corresponding post_factor value
if (i == 1) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * (perm_list %*% z) / tau)
} else {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * (perm_list[m, ] %*% z) / tau)
}
} # perm_list[m, ] = mth row of perm_list = mth unique permutaton of ith integer combo
}
val <- 0
for (i in 1:N_terms) { # add together the (pre_factor * post_factor) for each integer combo
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 1) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 8
}
int_vecs <- matrix(0, nrow = Npulse, ncol = N_terms)
for (i in 1:5) {
int_vecs[1:2, i] <- sign(k) * c(kk + (i - 1), -(i - 1))
}
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(1, 1, -1)
int_vecs[1:3, 7] <- sign(k) * c(3, -1, -1)
int_vecs[1:3, 8] <- sign(k) * c(2, -2, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
}
else if (abs(k) == 2) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 6
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk + 2, -2)
int_vecs[1:2, 5] <- sign(k) * c(kk + 3, -3)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(2, -1, 1)
}
if (Npulse >= 4) {
int_vecs[1:4, 7] <- sign(k) * c(1, 1, 1, -1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 2) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 3) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 8
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk + 2, -2)
int_vecs[1:2, 5] <- sign(k) * c(kk + 3, -3)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(1, 1, 1)
int_vecs[1:3, 7] <- sign(k) * c(3, 1, -1)
int_vecs[1:3, 8] <- sign(k) * c(2, 2, -1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 3 && i == 6) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 4) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 7
} else {
N_terms <- 8
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(4, -1, 1)
int_vecs[1:3, 7] <- sign(k) * c(2, 1, 1)
}
if (Npulse >= 4) {
int_vecs[1:4, 8] <- sign(k) * c(1, 1, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 4 && i == 8) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else if (Npulse == 2 && i == 3) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 5) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(2, 2, 1)
int_vecs[1:3, 7] <- sign(k) * c(3, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 6) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 6
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk - 3, 3)
int_vecs[1:2, 5] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 6] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 7] <- sign(k) * c(4, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 4) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 7) {
kk <- abs(k)
N_terms <- 6
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk - 3, 3)
int_vecs[1:2, 5] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 6] <- sign(k) * c(kk + 2, -2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 8) {
kk <- abs(k)
N_terms <- 6
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk - 3, +3)
int_vecs[1:2, 6] <- sign(k) * c(kk - 4, +4)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 6) {
Nperm <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 9) {
kk <- abs(k)
N_terms <- 5
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 5] <- sign(k) * c(kk - 3, +3)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[i] * post_factors[i]
}
return(val)
} else if (abs(k) == 10) {
kk <- abs(k)
N_terms <- 4
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 11) {
kk <- abs(k)
N_terms <- 4
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
out <- val
} else {
kk <- abs(k)
N_terms <- 3
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[i] * post_factors[i]
}
return(val)
}
}
# this function computes the fourier sum for 1 pulse
sum_1_pulse <- function(z, rho, mu, tau, kmax, modes) {
jmax <- length(modes) # modes is the list of fourier modes, jmax = number of list elements
out_val <- 0 # running total for output value
P_func <- function(x) {
rho * 1 / (log(1 / (1 - rho)) - 2 * pi * complex(real = 0, imaginary = 1) * x)
}
Q_func <- function(x) {
(rho / (1 - rho)) * 1 / (log(1 / (1 - rho)) + 2 * pi * complex(real = 0, imaginary = 1) * x)
}
for (k in -kmax:kmax) {
pre_factor <- P_func(k)
post_factor <- exp(-2 * pi * complex(real = 0, imaginary = 1) * k * z / tau)
P_val <- pre_factor * post_factor
Q_val <- 0
for (j in (-jmax + 1):(jmax - 1)) {
pre_factor <- Q_func(-k - j)
post_factor <- exp(-2 * pi * complex(real = 0, imaginary = 1) * (-k - j) * z / tau)
if (j == 0) {
Q_val <- Q_val + pre_factor * post_factor * modes[1] / mu
} else if (j < 0) {
Q_val <- Q_val + pre_factor * post_factor * Conj(modes[-(j - 1)]) / mu
} else if (j > 0) {
Q_val <- Q_val + pre_factor * post_factor * modes[j + 1] / mu
}
}
out_val <- out_val + (P_val) * (Q_val) / (1 - log(1 - rho) / (mu * tau) - 2 * pi * complex(real = 0, imaginary = 1) * k / (mu * tau))
}
return(Re(out_val))
}
# this function computes the Fourier sum for N pulses, where the "P" and "Q" products and exponentials are computeded by the "P_vals" and "Q_vals" functions
sum_N_pulse <- function(z, rho, Npulse, mu, tau, kmax, modes) {
jmax <- length(modes)
out_val <- 0 # running total for output
for (k in -kmax:kmax) {
P_val <- P_vals(z, rho, Npulse, tau, k) # compute "P" product and exponentials at this k value
Q_val <- 0 # running total for "Q" product and exponentials
# add together "Q" products and exponentials at (k + j) value for all j emergence fourier modes
for (j in (-jmax + 1):(jmax - 1)) {
# if (abs(k+j) < kmax){
if (j == 0) {
Q_val <- Q_val + Q_vals(z, rho, Npulse, tau, k, j) * modes[1] / mu # zero fourier mode
} else if (j < 0) {
Q_val <- Q_val + Q_vals(z, rho, Npulse, tau, k, j) * Conj(modes[-(j - 1)]) / mu # negative fourier modes
} else if (j > 0) {
Q_val <- Q_val + Q_vals(z, rho, Npulse, tau, k, j) * modes[j + 1] / mu # positive fourier modes
}
# }
}
# add contribution from this k value to running output total
out_val <- out_val + (P_val) * (Q_val) / (1 - Npulse * log(1 - rho) / (mu * tau) - 2 * pi * complex(real = 0, imaginary = 1) * k / (mu * tau))
}
out <- Re(out_val)
return(out)
}
## optimization parameters set by user
kmax <- 20 # max number of dynamics fourier modes to use in calculating fourier sum (different than N_lam = max emergence fourier mode set by user for curve fitting portion of the code. I've often refered to N_lam as jmax when we've talked. Hopefully not too confusing). kmax should be an integer between 2 and 200, defualt at 20
global_opt <- 0 # set to 0 if user chooses local optimum, 1 if user chooses golbal GN_DIRECT_L_RAND method, 2 if user chooses global GN_ISRES mthod
Npulse <- 4 # number of pulses, set by user, integer between 1 and 10
rho <- .30 # percent knockdown (user set between .01 and .30, e.g. 1% to 30% knockdown)
days_between <- 3 # minimum number of days allowed between pulses set by user (integer bewtween 0 and 30 days)
# code to execute optimization algorithms
# install.packages("nloptr")
library("nloptr") # used for local cobyla optimum function and global optimum functions
if (Npulse == 1) { # 1 pulse
fun1 <- function(x) sum_1_pulse(x, rho, mu, tau, kmax, lam_fourier)
guess <- tau / 2
if (global_opt == 0) { # local optimum
opt_out <- cobyla(guess, fun1, lower = 0, upper = tau)
times <- opt_out$par
ave_pop_fourier <- opt_out$value
} else if (global_opt == 1) { # global optimum directL
opt_out <- directL(fun1, lower = 0, upper = tau, randomized = TRUE)
times <- opt_out$par
ave_pop_fourier <- opt_out$value
} else if (global_opt == 2) { # global optimum isres
opt_out <- nloptr(x0 = guess, eval_f = fun1, lb = 0, ub = tau, eval_g_ineq = NULL, eval_g_eq = NULL, opts = list(algorithm = "NLOPT_GN_ISRES", maxeval = 10000))
times <- opt_out$solution
ave_pop_fourier <- opt_out$objective
}
} else { # 2 or more pulses
funN <- function(x) sum_N_pulse(x, rho, Npulse, mu, tau, kmax, lam_fourier)
# matrix and vector for inequality constraints to enforce minimum days between
A_mat <- matrix(0, nrow = Npulse - 1, ncol = Npulse)
b_vec <- (-1) * days_between * matrix(1, Npulse - 1, 1)
for (i in 1:(Npulse - 1)) {
A_mat[i, i] <- 1
A_mat[i, i + 1] <- -1
}
l_bound <- numeric(Npulse) # column vector of zeros for pulse lower bound
u_bound <- tau * matrix(1, Npulse, 1) # column vector of taus for pulse upper bound
guess <- numeric(Npulse) # vector for initial guesses for fmincon
for (i in 1:Npulse) {
if (t_in[1] - m / mu > 0) {
guess[i] <- i * (t_dat[N_pts - 1] - t_dat[3]) / (Npulse + 1) + t_dat[3]
} else {
guess[i] <- i * (t_dat[N_pts - 2] - t_dat[2]) / (Npulse + 1) + t_dat[2]
}
}
if (global_opt == 0) { # local optimum
opt_out <- fmincon(guess, funN, gr = NULL, method = "SQP", A = A_mat, b = b_vec, Aeq = NULL, beq = NULL, lb = l_bound, ub = u_bound)
times <- opt_out$par
ave_pop_fourier <- opt_out$value
} else if (global_opt == 1) { # global optimum directL
opt_out <- directL(funN, lower = l_bound, upper = u_bound, randomized = TRUE)
times <- opt_out$par
ave_pop_fourier <- opt_out$value
} else if (global_opt == 2) { # global optimum isres
ineq <- function(x) {
(-1) * A_mat %*% x - b_vec
}
opt_out <- nloptr(x0 = guess, eval_f = funN, lb = l_bound, ub = u_bound, eval_g_ineq = ineq, eval_g_eq = NULL, opts = list(algorithm = "NLOPT_GN_ISRES", maxeval = 10000))
times <- opt_out$solution
ave_pop_fourier <- opt_out$objective
}
}
# code to shift times and generate plot data
# shift pulse times back to original times input by user
if (t_in[1] - (m / mu) > 0) {
pulse_times_output <- times - t_dat[3] + t_in[1]
} else {
pulse_times_output <- times - t_dat[2] + t_in[1]
}
# code to plot results:
t_steps <- 6 * tau * 10 + 1 # ten time steps per day, six total seasons (need to integrate over many seasons to reach periodic population curves)
t_vec <- linspace(0, 6 * tau, n = t_steps)
# list between 0 and 6*tau, t_steps long
y_fourier_controlled <- numeric(t_steps) # list (row vector) for controlled population values at each time step
y_fourier_uncontrolled <- numeric(t_steps) # list (row vector) for uncontrolled population values at each time step
# simple integrator to calculate population values
for (i in 2:t_steps) {
for (j in 1:N_pts - 1) {
if (mod(t_vec[i - 1], tau) >= t_dat[j] && mod(t_vec[i - 1], tau) < t_dat[j + 1]) {
y_fourier_controlled[i] <- y_fourier_controlled[i - 1] + (-mu * y_fourier_controlled[i - 1] + lambda_dat$par[j + 1] * (mod(t_vec[i - 1], tau) - t_dat[j]) / delta_t_dat[j] + lambda_dat$par[j] * (t_dat[j + 1] - mod(t_vec[i - 1], tau)) / delta_t_dat[j]) * (t_vec[i] - t_vec[i - 1])
y_fourier_uncontrolled[i] <- y_fourier_uncontrolled[i - 1] + (-mu * y_fourier_uncontrolled[i - 1] + lambda_dat$par[j + 1] * (mod(t_vec[i - 1], tau) - t_dat[j]) / delta_t_dat[j] + lambda_dat$par[j] * (t_dat[j + 1] - mod(t_vec[i - 1], tau)) / delta_t_dat[j]) * (t_vec[i] - t_vec[i - 1])
}
}
# impulses for the controlled population
for (j in 1:Npulse) {
if (mod(t_vec[i], tau) >= times[j] && mod(t_vec[i - 1], tau) < times[j]) {
y_fourier_controlled[i] <- (1 - rho) * y_fourier_controlled[i]
}
}
}
pop_cont <- numeric(10 * tau + 1) # list for controlled population vlaues, one season long
pop_un_cont <- numeric(10 * tau + 1) # list for uncontrolled population vlaues, one season long
t_vec_plot <- numeric(10 * tau + 1) # list for time values, one season long
for (i in 1:(10 * tau + 1)) {
pop_cont[i] <- y_fourier_controlled[i + 10 * tau * 5] # take controlled population values from the integration for the final season
pop_un_cont[i] <- y_fourier_uncontrolled[i + 10 * tau * 5] # take uncontrolled population values from the integration for the final season
# shift plot times for population curves back to original times input by user
if (t_in[1] - m / mu > 0) {
t_vec_plot[i] <- t_vec[i] - t_dat[(3)] + t_in[(1)]
} else {
t_vec_plot[i] <- t_vec[i] - t_dat[(2)] + t_in[(1)]
}
}
# shift data times back to orginal times input by user
if (t_in[(1)] - (m / mu) > 0) {
t_dat_plot <- t_dat - t_dat[(3)] + t_in[(1)]
} else {
t_dat_plot <- t_dat - t_dat[2] + t_in[1]
}
# code to generate outputs to user
# install.packages("sfsmisc")
library("sfsmisc") # for the integrate.xy function
ave_pop_un_cont <- integrate.xy(t_vec_plot, pop_un_cont) / (tau)
ave_pop_cont <- integrate.xy(t_vec_plot, pop_cont) / (tau)
percent_reduction <- (ave_pop_un_cont - ave_pop_cont) / ave_pop_un_cont
accuracy_measure <- (ave_pop_fourier - ave_pop_cont) / ave_pop_cont
# outputs to display to user
print(pulse_times_output) # times of optimal pulses in units of day of year
print(ave_pop_un_cont) # average population over interval (0, tau) under uncontrolledfourier dynamics
print(ave_pop_cont) # average population over interval (0, tau) under controlled fourier dynamics
print(percent_reduction) # fractional reduction in average controlled population relative to uncontrolled
print(accuracy_measure) # reliability measuere - if the magnitude of this quantity is larger than .05 or .1 or so, then the fourier sum may be of unrelaiable accuracy, and the user should try increasing kmax
# plots to display to user
plot(t_dat_plot, y_dat,
ylim = c(0, 1.5 * max(y_dat)), xlim = c(t_dat_plot[1], t_dat_plot[1] + tau), col = "blue",
main = "Fitted Population Model",
ylab = "Mosquito population count",
xlab = "Day of year"
)
lines(t_vec_plot, pop_un_cont, col = "red")
lines(t_vec_plot, pop_cont, col = "orange")
legend("topleft", c("Data", "Uncontrolled Fourier Approximation", "Controlled Fourier Approximation"), fill = c("blue", "red", "orange"))
Try the mosqcontrol package in your browser
Any scripts or data that you put into this service are public.
mosqcontrol documentation built on Aug. 7, 2020, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# install.packages("pracma")
# install.packages("matlib")
# install.packages("NlcOptim")
library("pracma")
library("matlib")
# sample data
# mosquito count from data
y_in <- c(
15, 40, 45, 88, 99, 145, 111, 132, 177, 97, 94, 145, 123, 111,
125, 115, 155, 160, 143, 132, 126, 125, 105, 98, 87, 54, 55, 8
)
# day of year in which measurments making up yin occur, Jan 1st = day 1
t_in_user <- c(
93, 100, 107, 114, 121, 128, 135, 142, 149, 163, 170, 177,
184, 191, 198, 205, 212, 219, 226, 233, 240, 247, 254, 261,
267, 274, 281, 288
)
if (length(y_in) != length(t_in_user)) {
stop("input data mismatch")
}
# User inputs (aside from trap data inputs)
# mu = natural mosquito death rate
mu <- 1 / 14 # program should set mu = 1 / tau_mosq, where tau_mosq is the mosquito lifetime input by user in units of days
# tau_mosq should be bounded between 1 and 30
# m = number of mosquito lifetimes for population decay between seasons..
m <- 3 # this will be input by the user, and must be an integer greater than 2, upperbound 100, default value 3
# N_lam = max fourier mode order to calculate
N_lam <- 25
# user input between 1 and 200, default value 25 (I referred to this as jmax over Skype - JD)
# setting up time and trap count vectors to fit the model to
# first round user time inputs to integer values
t_in <- round(t_in_user)
# defining N0 = number of input data points (not a user input)
N0 <- length(y_in)
# defining vector of time differences to be used later (not a user input)
delta_t_in <- numeric(N0 - 1)
for (i in 1:(N0 - 1)) {
delta_t_in[i] <- t_in[i + 1] - t_in[i]
if (delta_t_in[i] < 0) {
stop("non-increasing measurement times")
}
}
# N_pts = number of points in extended data set, N_pts = N0 + 3 always (not a user input)
N_pts <- N0 + 3
# defining extended trap count and time data vectors
t_dat <- numeric(N_pts)
if (t_in[1] - m / mu > 0) {
y_dat <- c(0, 0, y_in, 0)
t_dat[1] <- 0
t_dat[2] <- (m - 1) / mu
for (i in 3:(N_pts - 1)) {
t_dat[i] <- t_in[i - 2] - t_in[1] + m / mu
}
t_dat[N_pts] <- t_dat[N_pts - 1] + (m) / mu
t_dat <- matrix(t_dat, ncol = 1)
} else {
y_dat <- c(0, y_in, 0, 0)
t_dat[1] <- 0
for (i in 2:(N_pts - 2)) {
t_dat[i] <- t_in[i - 1]
}
t_dat[N_pts - 1] <- t_dat[N_pts - 2] + (2 * (m - 1)) / mu - t_dat[2]
t_dat[N_pts] <- t_dat[N_pts - 1] + 1 / mu
}
# extended vector of time differences
delta_t_dat <- numeric(N_pts)
for (i in 1:(N_pts - 1)) {
delta_t_dat[i] <- t_dat[i + 1] - t_dat[i]
}
delta_t_dat[N_pts] <- 0
delta_t_dat <- matrix(delta_t_dat, ncol = 1)
# tau = length of 'season' in days
tau <- t_dat[N_pts]
# Defining matrices 'J' and 'M'
# J matrix and inverse
J <- matrix(0, nrow = N_pts, ncol = N_pts)
for (i in 1:(N_pts - 1)) {
J[i, i] <- 1
J[i + 1, i] <- -exp(-mu * delta_t_dat[i])
}
J[N_pts, N_pts] <- 1
J[1, N_pts] <- -exp(-mu * delta_t_dat[N_pts])
# M matrix and inverse
M <- matrix(0, nrow = N_pts, ncol = N_pts)
for (i in 2:(N_pts - 1)) {
M[i, i] <- (1 - exp(-mu * delta_t_dat[i - 1])) / (mu * delta_t_dat[i - 1]) - exp(-mu * delta_t_dat[i - 1])
M[i + 1, i] <- 1 - (1 - exp(-mu * delta_t_dat[i])) / (mu * delta_t_dat[i])
}
M[N_pts, N_pts] <- (1 - exp((-1) * mu * delta_t_dat[N_pts - 1])) / (mu * delta_t_dat[(N_pts - 1)]) - exp(-mu * delta_t_dat[(N_pts - 1)])
M[1, N_pts] <- 0
M[1, 1] <- 0
M[2, 1] <- 1 - (1 - exp(-mu * delta_t_dat[1])) / (mu * delta_t_dat[(1)])
###### weighting code code
# Weight matrix: diagonal entries to be filled with sqrt(number of counts on day) except for the artifical 0 measurments on end of data
W <- matrix(0, nrow = N_pts, ncol = N_pts)
if (t_in[(1)] - m / mu > 0) {
W[1, 1] <- 1
W[2, 2] <- 1
W[N_pts, N_pts] <- 1
for (i in 3:N_pts - 1) {
W[i, i] <- sqrt(1) # fill sqrt(number of counts) here
}
} else {
W[1, 1] <- 1
W[N_pts - 1, N_pts - 1] <- 1
W[N_pts, N_pts] <- 1
for (i in 2:N_pts - 2) {
W[i, i] <- sqrt(1) # fill sqrt(number of counts) here
}
}
# Obtaining emergence rate vector lambda_dat by contrained optimization
Aeq <- matrix(numeric(N_pts), nrow = 1)
Aeq[(1)] <- 1
Aeq[(N_pts)] <- -1
beq <- 0
bineq <- numeric(N_pts)
lb <- numeric(N_pts)
###################
objective_fun <- function(x) {
t(W %*% ((1 / mu) * (mldivide(J, (M %*% x), pinv = TRUE)) - y_dat)) %*% (W %*% ((1 / mu) * (mldivide(J, (M %*% x), pinv = TRUE)) - y_dat))
}
lambda_dat <- fmincon(numeric(N_pts), objective_fun, gr = NULL, method = "SQP", A = (-(1 / mu) * (inv(J) %*% (M))), b = bineq, Aeq = Aeq, beq = beq, lb = lb, ub = NULL)
#########
# Calculating emergence rate fourier modes lam_fourier
lam_fourier <- numeric(N_lam + 1)
for (j in 1:(N_pts - 1)) {
lam_fourier[1] <- lam_fourier[1] + (1 / tau) * (lambda_dat$par[j] + lambda_dat$par[j + 1]) * delta_t_dat[j] / 2
} # zero mode
for (k in 2:(N_lam + 1)) # non-zero postive modes (negative modes given by complex conjugates of positve modes)
{
for (j in 1:(N_pts - 1)) {
lam_fourier[(k)] <- lam_fourier[(k)] + ((lambda_dat$par[(j + 1)] - lambda_dat$par[(j)]) / delta_t_dat[(j)]) * (tau / (2 * pi * (k - 1))) * ((exp(-2 * pi * complex(real = 0, imaginary = 1) * (k - 1) * t_dat[(j + 1)] / tau) - exp(-2 * pi * complex(real = 0, imaginary = 1) * (k - 1) * t_dat[(j)] / tau)) / (2 * pi * (k - 1))) + complex(real = 0, imaginary = 1) * (lambda_dat$par[(j + 1)] * exp(-2 * pi * complex(real = 0, imaginary = 1) * (k - 1) * t_dat[(j + 1)] / tau) - lambda_dat$par[(j)] * exp(-2 * pi * complex(real = 0, imaginary = 1) * (k - 1) * t_dat[(j)] / tau)) / (2 * pi * (k - 1))
}
}
# defing functions for optimization code
# goto line 1670 or so to skip past all of these definitons
# Permutaion function: For a given list of numbers, this function outputs a matrix, where each row is a unique permutation of the list
uperm <- function(d) {
dat <- factor(d)
N <- length(dat)
n <- tabulate(dat)
ng <- length(n)
if (ng == 1) {
return(d)
}
a <- N - c(0, cumsum(n))[-(ng + 1)]
foo <- lapply(1:ng, function(i) matrix(combn(a[i], n[i]), nrow = n[i]))
out <- matrix(NA, nrow = N, ncol = prod(sapply(foo, ncol)))
xxx <- c(0, cumsum(sapply(foo, nrow)))
xxx <- cbind(xxx[-length(xxx)] + 1, xxx[-1])
miss <- matrix(1:N, ncol = 1)
for (i in seq_len(length(foo) - 1)) {
l1 <- foo[[i]]
nn <- ncol(miss)
miss <- matrix(rep(miss, ncol(l1)), nrow = nrow(miss))
k <- (rep(0:(ncol(miss) - 1), each = nrow(l1))) * nrow(miss) +
l1[, rep(1:ncol(l1), each = nn)]
out[xxx[i, 1]:xxx[i, 2], ] <- matrix(miss[k], ncol = ncol(miss))
miss <- matrix(miss[-k], ncol = ncol(miss))
}
k <- length(foo)
out[xxx[k, 1]:xxx[k, 2], ] <- miss
out <- out[rank(as.numeric(dat), ties = "first"), ]
foo <- cbind(as.vector(out), as.vector(col(out)))
out[foo] <- d
t(out)
}
# this function computes the "P" product and exponentials for a given k
P_vals <- function(z, rho, Npulse, tau, k) {
P_func <- function(x) {
rho * 1 / (log(1 / (1 - rho)) - 2 * pi * complex(real = 0, imaginary = 1) * x)
}
z_vec <- matrix(0, nrow = Npulse, ncol = 1)
for (i in 1:Npulse) {
z_vec[i, 1] <- z[i]
}
# determine number of integer vectors to permute (will be smaller for
# only two and three pulses). number of combos = N_terms
if (k == 0) {
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 7
} else {
N_terms <- 8
}
# create a matrix to hold the integer vectors to permute, and
# initialize to zero. Each column will give an integer vector
int_vecs <- matrix(0, nrow = Npulse, ncol = N_terms)
# the first 5 integer vectors have at most two non-zero elements.
# This loop replaces the first two zeros in the first 5 columns of int_vecs
# with the appropriate integers.
for (i in 1:5) {
int_vecs[1:2, i] <- c((i - 1), -(i - 1))
}
# if there are are three or more pulses, we need to include
# integer vectors with three non-zero terms
if (Npulse >= 3) {
int_vecs[1:3, 6] <- c(2, -1, -1)
int_vecs[1:3, 7] <- c(-2, 1, 1)
}
# if there are four or more pulses, then we need to include an
# integer vector with four non-zero terms
if (Npulse >= 4) {
int_vecs[1:4, 8] <- c(1, 1, -1, -1)
}
# calculates the "P" functions for each column of integer combos in
# int_vec. "P" function values are independet of integer
# permutations. pre_factors is a list which holds the "P"-function
# values for each integer combo
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
# Now calculate the exponential terms which will multiply the "P"functions in the pre_factors list
# post_factors = list of the exponential terms to multiply the
# corresponding "P" functions in the pre_factors list
post_factors <- numeric(N_terms) # initialize to zero
for (i in 1:N_terms) { # for each integer combo given by the columns of int_vecs...
# calculates the unique permutations of the ith integer combo
# int_vecs(:, i) = ith column of int_vecs matrix
# perm_list = Nperms x Npulse matrix, where Nperms is the number
# of unique permutations of the ith integer combo
# Each row of perm_list gives a unique permutaion of the ith integer combo
perm_list <- uperm(int_vecs[, i])
if (i == 1) {
Nperms <- 1
} else {
Nperms <- length(perm_list[, 1]) # number of unique permutations of the ith integer combo
}
for (m in 1:Nperms) { # for each uniqe permutation of the ith integer combo, add the corresponding exponential to the running total of the corresponding post_factor value
if (i == 1) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * (perm_list %*% z) / tau)
} else {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * (perm_list[m, ] %*% z) / tau)
}
} # perm_list[m, ] = mth row of perm_list = mth unique permutaton of ith integer combo
}
val <- 0
for (i in 1:N_terms) { # add together the (pre_factor * post_factor) for each integer combo
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 1) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 8
}
int_vecs <- matrix(0, nrow = Npulse, ncol = N_terms)
for (i in 1:5) {
int_vecs[1:2, i] <- sign(k) * c(kk + (i - 1), -(i - 1))
}
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(1, 1, -1)
int_vecs[1:3, 7] <- sign(k) * c(3, -1, -1)
int_vecs[1:3, 8] <- sign(k) * c(2, -2, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
}
else if (abs(k) == 2) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 6
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk + 2, -2)
int_vecs[1:2, 5] <- sign(k) * c(kk + 3, -3)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(2, -1, 1)
}
if (Npulse >= 4) {
int_vecs[1:4, 7] <- sign(k) * c(1, 1, 1, -1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 2) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 3) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 8
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk + 2, -2)
int_vecs[1:2, 5] <- sign(k) * c(kk + 3, -3)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(1, 1, 1)
int_vecs[1:3, 7] <- sign(k) * c(3, 1, -1)
int_vecs[1:3, 8] <- sign(k) * c(2, 2, -1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 3 && i == 6) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 4) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 7
} else {
N_terms <- 8
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(4, -1, 1)
int_vecs[1:3, 7] <- sign(k) * c(2, 1, 1)
}
if (Npulse >= 4) {
int_vecs[1:4, 8] <- sign(k) * c(1, 1, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 4 && i == 8) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else if (Npulse == 2 && i == 3) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 5) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(2, 2, 1)
int_vecs[1:3, 7] <- sign(k) * c(3, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 6) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 6
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk - 3, 3)
int_vecs[1:2, 5] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 6] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 7] <- sign(k) * c(4, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 4) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 7) {
kk <- abs(k)
N_terms <- 6
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk - 3, 3)
int_vecs[1:2, 5] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 6] <- sign(k) * c(kk + 2, -2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 8) {
kk <- abs(k)
N_terms <- 6
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk - 3, +3)
int_vecs[1:2, 6] <- sign(k) * c(kk - 4, +4)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 6) {
Nperm <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 9) {
kk <- abs(k)
N_terms <- 5
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 5] <- sign(k) * c(kk - 3, +3)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[i] * post_factors[i]
}
return(val)
} else if (abs(k) == 10) {
kk <- abs(k)
N_terms <- 4
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 11) {
kk <- abs(k)
N_terms <- 4
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
out <- val
} else {
kk <- abs(k)
N_terms <- 3
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * P_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[i] * post_factors[i]
}
return(val)
}
}
# this function computes the "Q" product and exponentials with the fourier modes for a given k
Q_vals <- function(z, rho, Npulse, tau, j, jj) {
Q_func <- function(x) {
(rho / (1 - rho)) * 1 / (log(1 / (1 - rho)) + 2 * pi * complex(real = 0, imaginary = 1) * x)
}
k <- -(j + jj)
# determine number of integer vectors to permute (will be smaller for
# only two and three pulses). number of combos = N_terms
if (k == 0) {
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 7
} else {
N_terms <- 8
}
# create a matrix to hold the integer vectors to permute, and
# initialize to zero. Each column will give an integer vector
int_vecs <- matrix(0, nrow = Npulse, ncol = N_terms)
# the first 5 integer vectors have at most two non-zero elements.
# This loop replaces the first two zeros in the first 5 columns of int_vecs
# with the appropriate integers.
for (i in 1:5) {
int_vecs[1:2, i] <- c((i - 1), -(i - 1))
}
# if there are are three or more pulses, we need to include
# integer vectors with three non-zero terms
if (Npulse >= 3) {
int_vecs[1:3, 6] <- c(2, -1, -1)
int_vecs[1:3, 7] <- c(-2, 1, 1)
}
# if there are four or more pulses, then we need to include an
# integer vector with four non-zero terms
if (Npulse >= 4) {
int_vecs[1:4, 8] <- c(1, 1, -1, -1)
}
# calculates the "P" functions for each column of integer combos in
# int_vec. "P" function values are independet of integer
# permutations. pre_factors is a list which holds the "P"-function
# values for each integer combo
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
# Now calculate the exponential terms which will multiply the "P"functions in the pre_factors list
# post_factors = list of the exponential terms to multiply the
# corresponding "P" functions in the pre_factors list
post_factors <- numeric(N_terms) # initialize to zero
for (i in 1:N_terms) { # for each integer combo given by the columns of int_vecs...
# calculates the unique permutations of the ith integer combo
# int_vecs(:, i) = ith column of int_vecs matrix
# perm_list = Nperms x Npulse matrix, where Nperms is the number
# of unique permutations of the ith integer combo
# Each row of perm_list gives a unique permutaion of the ith integer combo
perm_list <- uperm(int_vecs[, i])
if (i == 1) {
Nperms <- 1
} else {
Nperms <- length(perm_list[, 1]) # number of unique permutations of the ith integer combo
}
for (m in 1:Nperms) { # for each uniqe permutation of the ith integer combo, add the corresponding exponential to the running total of the corresponding post_factor value
if (i == 1) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * (perm_list %*% z) / tau)
} else {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * (perm_list[m, ] %*% z) / tau)
}
} # perm_list[m, ] = mth row of perm_list = mth unique permutaton of ith integer combo
}
val <- 0
for (i in 1:N_terms) { # add together the (pre_factor * post_factor) for each integer combo
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 1) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 8
}
int_vecs <- matrix(0, nrow = Npulse, ncol = N_terms)
for (i in 1:5) {
int_vecs[1:2, i] <- sign(k) * c(kk + (i - 1), -(i - 1))
}
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(1, 1, -1)
int_vecs[1:3, 7] <- sign(k) * c(3, -1, -1)
int_vecs[1:3, 8] <- sign(k) * c(2, -2, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
}
else if (abs(k) == 2) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 6
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk + 2, -2)
int_vecs[1:2, 5] <- sign(k) * c(kk + 3, -3)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(2, -1, 1)
}
if (Npulse >= 4) {
int_vecs[1:4, 7] <- sign(k) * c(1, 1, 1, -1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 2) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 3) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 8
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk + 2, -2)
int_vecs[1:2, 5] <- sign(k) * c(kk + 3, -3)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(1, 1, 1)
int_vecs[1:3, 7] <- sign(k) * c(3, 1, -1)
int_vecs[1:3, 8] <- sign(k) * c(2, 2, -1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 3 && i == 6) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 4) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else if (Npulse == 3) {
N_terms <- 7
} else {
N_terms <- 8
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(4, -1, 1)
int_vecs[1:3, 7] <- sign(k) * c(2, 1, 1)
}
if (Npulse >= 4) {
int_vecs[1:4, 8] <- sign(k) * c(1, 1, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 4 && i == 8) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else if (Npulse == 2 && i == 3) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 5) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 5
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 6] <- sign(k) * c(2, 2, 1)
int_vecs[1:3, 7] <- sign(k) * c(3, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 6) {
kk <- abs(k)
if (Npulse == 2) {
N_terms <- 6
} else {
N_terms <- 7
}
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk - 3, 3)
int_vecs[1:2, 5] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 6] <- sign(k) * c(kk + 2, -2)
if (Npulse >= 3) {
int_vecs[1:3, 7] <- sign(k) * c(4, 1, 1)
}
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 4) {
Nperms <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 7) {
kk <- abs(k)
N_terms <- 6
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk - 3, 3)
int_vecs[1:2, 5] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 6] <- sign(k) * c(kk + 2, -2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 8) {
kk <- abs(k)
N_terms <- 6
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 4] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 5] <- sign(k) * c(kk - 3, +3)
int_vecs[1:2, 6] <- sign(k) * c(kk - 4, +4)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
if (Npulse == 2 && i == 6) {
Nperm <- 1
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list %*% z / tau)
} else {
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 9) {
kk <- abs(k)
N_terms <- 5
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
int_vecs[1:2, 5] <- sign(k) * c(kk - 3, +3)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[i] * post_factors[i]
}
return(val)
} else if (abs(k) == 10) {
kk <- abs(k)
N_terms <- 4
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
return(val)
} else if (abs(k) == 11) {
kk <- abs(k)
N_terms <- 4
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
int_vecs[1:2, 4] <- sign(k) * c(kk - 2, 2)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[1, i] * post_factors[i]
}
out <- val
} else {
kk <- abs(k)
N_terms <- 3
int_vecs <- matrix(0, Npulse, N_terms)
int_vecs[1:2, 1] <- sign(k) * c(kk + 0, 0)
int_vecs[1:2, 2] <- sign(k) * c(kk - 1, 1)
int_vecs[1:2, 3] <- sign(k) * c(kk + 1, -1)
pre_factors <- matrix(1, ncol = N_terms, nrow = 1)
for (i in 1:N_terms) {
for (m in 1:Npulse) {
pre_factors[1, i] <- pre_factors[1, i] * Q_func(int_vecs[m, i])
}
}
post_factors <- numeric(N_terms)
for (i in 1:N_terms) {
perm_list <- uperm(int_vecs[, i])
Nperms <- dim(perm_list)[1]
for (m in 1:Nperms) {
post_factors[i] <- post_factors[i] + exp(-2 * pi * complex(real = 0, imaginary = 1) * perm_list[m, ] %*% z / tau)
}
}
val <- 0
for (i in 1:N_terms) {
val <- val + pre_factors[i] * post_factors[i]
}
return(val)
}
}
# this function computes the fourier sum for 1 pulse
sum_1_pulse <- function(z, rho, mu, tau, kmax, modes) {
jmax <- length(modes) # modes is the list of fourier modes, jmax = number of list elements
out_val <- 0 # running total for output value
P_func <- function(x) {
rho * 1 / (log(1 / (1 - rho)) - 2 * pi * complex(real = 0, imaginary = 1) * x)
}
Q_func <- function(x) {
(rho / (1 - rho)) * 1 / (log(1 / (1 - rho)) + 2 * pi * complex(real = 0, imaginary = 1) * x)
}
for (k in -kmax:kmax) {
pre_factor <- P_func(k)
post_factor <- exp(-2 * pi * complex(real = 0, imaginary = 1) * k * z / tau)
P_val <- pre_factor * post_factor
Q_val <- 0
for (j in (-jmax + 1):(jmax - 1)) {
pre_factor <- Q_func(-k - j)
post_factor <- exp(-2 * pi * complex(real = 0, imaginary = 1) * (-k - j) * z / tau)
if (j == 0) {
Q_val <- Q_val + pre_factor * post_factor * modes[1] / mu
} else if (j < 0) {
Q_val <- Q_val + pre_factor * post_factor * Conj(modes[-(j - 1)]) / mu
} else if (j > 0) {
Q_val <- Q_val + pre_factor * post_factor * modes[j + 1] / mu
}
}
out_val <- out_val + (P_val) * (Q_val) / (1 - log(1 - rho) / (mu * tau) - 2 * pi * complex(real = 0, imaginary = 1) * k / (mu * tau))
}
return(Re(out_val))
}
# this function computes the Fourier sum for N pulses, where the "P" and "Q" products and exponentials are computeded by the "P_vals" and "Q_vals" functions
sum_N_pulse <- function(z, rho, Npulse, mu, tau, kmax, modes) {
jmax <- length(modes)
out_val <- 0 # running total for output
for (k in -kmax:kmax) {
P_val <- P_vals(z, rho, Npulse, tau, k) # compute "P" product and exponentials at this k value
Q_val <- 0 # running total for "Q" product and exponentials
# add together "Q" products and exponentials at (k + j) value for all j emergence fourier modes
for (j in (-jmax + 1):(jmax - 1)) {
# if (abs(k+j) < kmax){
if (j == 0) {
Q_val <- Q_val + Q_vals(z, rho, Npulse, tau, k, j) * modes[1] / mu # zero fourier mode
} else if (j < 0) {
Q_val <- Q_val + Q_vals(z, rho, Npulse, tau, k, j) * Conj(modes[-(j - 1)]) / mu # negative fourier modes
} else if (j > 0) {
Q_val <- Q_val + Q_vals(z, rho, Npulse, tau, k, j) * modes[j + 1] / mu # positive fourier modes
}
# }
}
# add contribution from this k value to running output total
out_val <- out_val + (P_val) * (Q_val) / (1 - Npulse * log(1 - rho) / (mu * tau) - 2 * pi * complex(real = 0, imaginary = 1) * k / (mu * tau))
}
out <- Re(out_val)
return(out)
}
## optimization parameters set by user
kmax <- 20 # max number of dynamics fourier modes to use in calculating fourier sum (different than N_lam = max emergence fourier mode set by user for curve fitting portion of the code. I've often refered to N_lam as jmax when we've talked. Hopefully not too confusing). kmax should be an integer between 2 and 200, defualt at 20
global_opt <- 0 # set to 0 if user chooses local optimum, 1 if user chooses golbal GN_DIRECT_L_RAND method, 2 if user chooses global GN_ISRES mthod
Npulse <- 4 # number of pulses, set by user, integer between 1 and 10
rho <- .30 # percent knockdown (user set between .01 and .30, e.g. 1% to 30% knockdown)
days_between <- 3 # minimum number of days allowed between pulses set by user (integer bewtween 0 and 30 days)
# code to execute optimization algorithms
# install.packages("nloptr")
library("nloptr") # used for local cobyla optimum function and global optimum functions
if (Npulse == 1) { # 1 pulse
fun1 <- function(x) sum_1_pulse(x, rho, mu, tau, kmax, lam_fourier)
guess <- tau / 2
if (global_opt == 0) { # local optimum
opt_out <- cobyla(guess, fun1, lower = 0, upper = tau)
times <- opt_out$par
ave_pop_fourier <- opt_out$value
} else if (global_opt == 1) { # global optimum directL
opt_out <- directL(fun1, lower = 0, upper = tau, randomized = TRUE)
times <- opt_out$par
ave_pop_fourier <- opt_out$value
} else if (global_opt == 2) { # global optimum isres
opt_out <- nloptr(x0 = guess, eval_f = fun1, lb = 0, ub = tau, eval_g_ineq = NULL, eval_g_eq = NULL, opts = list(algorithm = "NLOPT_GN_ISRES", maxeval = 10000))
times <- opt_out$solution
ave_pop_fourier <- opt_out$objective
}
} else { # 2 or more pulses
funN <- function(x) sum_N_pulse(x, rho, Npulse, mu, tau, kmax, lam_fourier)
# matrix and vector for inequality constraints to enforce minimum days between
A_mat <- matrix(0, nrow = Npulse - 1, ncol = Npulse)
b_vec <- (-1) * days_between * matrix(1, Npulse - 1, 1)
for (i in 1:(Npulse - 1)) {
A_mat[i, i] <- 1
A_mat[i, i + 1] <- -1
}
l_bound <- numeric(Npulse) # column vector of zeros for pulse lower bound
u_bound <- tau * matrix(1, Npulse, 1) # column vector of taus for pulse upper bound
guess <- numeric(Npulse) # vector for initial guesses for fmincon
for (i in 1:Npulse) {
if (t_in[1] - m / mu > 0) {
guess[i] <- i * (t_dat[N_pts - 1] - t_dat[3]) / (Npulse + 1) + t_dat[3]
} else {
guess[i] <- i * (t_dat[N_pts - 2] - t_dat[2]) / (Npulse + 1) + t_dat[2]
}
}
if (global_opt == 0) { # local optimum
opt_out <- fmincon(guess, funN, gr = NULL, method = "SQP", A = A_mat, b = b_vec, Aeq = NULL, beq = NULL, lb = l_bound, ub = u_bound)
times <- opt_out$par
ave_pop_fourier <- opt_out$value
} else if (global_opt == 1) { # global optimum directL
opt_out <- directL(funN, lower = l_bound, upper = u_bound, randomized = TRUE)
times <- opt_out$par
ave_pop_fourier <- opt_out$value
} else if (global_opt == 2) { # global optimum isres
ineq <- function(x) {
(-1) * A_mat %*% x - b_vec
}
opt_out <- nloptr(x0 = guess, eval_f = funN, lb = l_bound, ub = u_bound, eval_g_ineq = ineq, eval_g_eq = NULL, opts = list(algorithm = "NLOPT_GN_ISRES", maxeval = 10000))
times <- opt_out$solution
ave_pop_fourier <- opt_out$objective
}
}
# code to shift times and generate plot data
# shift pulse times back to original times input by user
if (t_in[1] - (m / mu) > 0) {
pulse_times_output <- times - t_dat[3] + t_in[1]
} else {
pulse_times_output <- times - t_dat[2] + t_in[1]
}
# code to plot results:
t_steps <- 6 * tau * 10 + 1 # ten time steps per day, six total seasons (need to integrate over many seasons to reach periodic population curves)
t_vec <- linspace(0, 6 * tau, n = t_steps)
# list between 0 and 6*tau, t_steps long
y_fourier_controlled <- numeric(t_steps) # list (row vector) for controlled population values at each time step
y_fourier_uncontrolled <- numeric(t_steps) # list (row vector) for uncontrolled population values at each time step
# simple integrator to calculate population values
for (i in 2:t_steps) {
for (j in 1:N_pts - 1) {
if (mod(t_vec[i - 1], tau) >= t_dat[j] && mod(t_vec[i - 1], tau) < t_dat[j + 1]) {
y_fourier_controlled[i] <- y_fourier_controlled[i - 1] + (-mu * y_fourier_controlled[i - 1] + lambda_dat$par[j + 1] * (mod(t_vec[i - 1], tau) - t_dat[j]) / delta_t_dat[j] + lambda_dat$par[j] * (t_dat[j + 1] - mod(t_vec[i - 1], tau)) / delta_t_dat[j]) * (t_vec[i] - t_vec[i - 1])
y_fourier_uncontrolled[i] <- y_fourier_uncontrolled[i - 1] + (-mu * y_fourier_uncontrolled[i - 1] + lambda_dat$par[j + 1] * (mod(t_vec[i - 1], tau) - t_dat[j]) / delta_t_dat[j] + lambda_dat$par[j] * (t_dat[j + 1] - mod(t_vec[i - 1], tau)) / delta_t_dat[j]) * (t_vec[i] - t_vec[i - 1])
}
}
# impulses for the controlled population
for (j in 1:Npulse) {
if (mod(t_vec[i], tau) >= times[j] && mod(t_vec[i - 1], tau) < times[j]) {
y_fourier_controlled[i] <- (1 - rho) * y_fourier_controlled[i]
}
}
}
pop_cont <- numeric(10 * tau + 1) # list for controlled population vlaues, one season long
pop_un_cont <- numeric(10 * tau + 1) # list for uncontrolled population vlaues, one season long
t_vec_plot <- numeric(10 * tau + 1) # list for time values, one season long
for (i in 1:(10 * tau + 1)) {
pop_cont[i] <- y_fourier_controlled[i + 10 * tau * 5] # take controlled population values from the integration for the final season
pop_un_cont[i] <- y_fourier_uncontrolled[i + 10 * tau * 5] # take uncontrolled population values from the integration for the final season
# shift plot times for population curves back to original times input by user
if (t_in[1] - m / mu > 0) {
t_vec_plot[i] <- t_vec[i] - t_dat[(3)] + t_in[(1)]
} else {
t_vec_plot[i] <- t_vec[i] - t_dat[(2)] + t_in[(1)]
}
}
# shift data times back to orginal times input by user
if (t_in[(1)] - (m / mu) > 0) {
t_dat_plot <- t_dat - t_dat[(3)] + t_in[(1)]
} else {
t_dat_plot <- t_dat - t_dat[2] + t_in[1]
}
# code to generate outputs to user
# install.packages("sfsmisc")
library("sfsmisc") # for the integrate.xy function
ave_pop_un_cont <- integrate.xy(t_vec_plot, pop_un_cont) / (tau)
ave_pop_cont <- integrate.xy(t_vec_plot, pop_cont) / (tau)
percent_reduction <- (ave_pop_un_cont - ave_pop_cont) / ave_pop_un_cont
accuracy_measure <- (ave_pop_fourier - ave_pop_cont) / ave_pop_cont
# outputs to display to user
print(pulse_times_output) # times of optimal pulses in units of day of year
print(ave_pop_un_cont) # average population over interval (0, tau) under uncontrolledfourier dynamics
print(ave_pop_cont) # average population over interval (0, tau) under controlled fourier dynamics
print(percent_reduction) # fractional reduction in average controlled population relative to uncontrolled
print(accuracy_measure) # reliability measuere - if the magnitude of this quantity is larger than .05 or .1 or so, then the fourier sum may be of unrelaiable accuracy, and the user should try increasing kmax
# plots to display to user
plot(t_dat_plot, y_dat,
ylim = c(0, 1.5 * max(y_dat)), xlim = c(t_dat_plot[1], t_dat_plot[1] + tau), col = "blue",
main = "Fitted Population Model",
ylab = "Mosquito population count",
xlab = "Day of year"
)
lines(t_vec_plot, pop_un_cont, col = "red")
lines(t_vec_plot, pop_cont, col = "orange")
legend("topleft", c("Data", "Uncontrolled Fourier Approximation", "Controlled Fourier Approximation"), fill = c("blue", "red", "orange"))
Try the mosqcontrol package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.