.Rhistory

library(flowkernel)
set.seed(123)
d <- 1; K <- 2; ntimes <- 200
ex1 <- list(
mu_function = function(times) {
mu <- array(NA, c(ntimes, K, d))
mu[, , 1] <- cbind(sin(2 * pi * times / 30), 2)
mu
},
Sigma_function = function(times) {
Sigma <- array(NA, c(ntimes, K, 1, 1))
Sigma[, , 1, 1] <- 0.25
Sigma
},
pi_function = function(times) {
pi1 <- seq(0.2, 0.8, length=length(times))
cbind(pi1, 1 - pi1)
},
num_points = rep(40, ntimes)
)
ex1$dat <- generate_smooth_gauss_mix(ex1$mu_function,
ex1$Sigma_function,
ex1$pi_function,
ex1$num_points)
d <- 3; K <- 4; ntimes <- 200
ex2 = list(
mu_function = function(times) {
mu <- array(NA, c(ntimes, K, d))
mu[, , 1] <- cbind(0.5*cos(2 * pi * times / 30), 0.3*sin(2 * pi * times / 30), sin(2 * pi * times / 30), -3)
mu[, , 2] = cbind (0.3*sin(2 * pi * times / 30), 2, -1, 0.6*cos(2 * pi * times / 30))
mu[, , 3] = cbind(2, 0.7*cos(2 * pi * times / 30), 0.4*sin(2 * pi * times / 30), 1)
mu
},
Sigma_function = function(times) {
Sigma <- array(NA, c(ntimes, K, d, d))
Sigma[, , 1, 1] <- 0.10
Sigma[, , 1, 2] <- 0
Sigma[, , 1, 3] <- 0
Sigma[, , 2, 1] <- 0
Sigma[, , 2, 2] <- 0.10
Sigma[, , 2, 3] <- 0
Sigma[, , 3, 1] <- 0
Sigma[, , 3, 2] <- 0
Sigma[, , 3, 3] <- 0.10
Sigma
},
pi_function = function(times) {
pi1 <- seq(0.2, 0.3, length=length(times))
cbind(pi1, pi1, 2*pi1/3, 1- (2*pi1 + 2*pi1/3))
},
num_points = rep(150, ntimes)
)
ex2$dat <- generate_smooth_gauss_mix(ex2$mu_function,
ex2$Sigma_function,
ex2$pi_function,
ex2$num_points)
plot_data_and_model(ex2$dat$y, ex2$dat$z, ex2$dat$mu, dim = 1)
plot_data_and_model(ex2$dat$y, ex2$dat$z, ex2$dat$mu, dim = 1, show_data = FALSE)
plot_data_and_model(ex2$dat$y, ex2$dat$z, ex2$dat$mu, dim = 2, show_data = FALSE)
plot_data_and_model(ex2$dat$y, ex2$dat$z, ex2$dat$mu, dim = 3, show_data = FALSE)
library(mclust)
library(plotly)
library(ggplot2)
library(tidyverse)
library(magrittr)
library(dplyr)
#change the directory below to where you have the data saved
grid_file_path <- "C:/Desktop Folders/USC/Academic/Research/paper-data/Justin-gridded-SeaFlow-2022-11-02/2022-11-01-diam-chl-pe-Qc-PSD-combined-popcycle-4.15.0-withoutBoundaryPoints.psd-grid.parquet"
# Read the parquet file
grid <- arrow::read_parquet(grid_file_path)
data_file_path <- "C:/Desktop Folders/USC/Academic/Research/paper-data/Justin-gridded-SeaFlow-2022-11-02/2022-11-01-diam-chl-pe-Qc-PSD-combined-popcycle-4.15.0-withoutBoundaryPoints.psd-hourly.parquet"
psd_hourly <- arrow::read_parquet(data_file_path)
# grid <- arrow::read_parquet("2022-11-01-diam-chl-pe-Qc-PSD-combined-popcycle-4.15.0-withoutBoundaryPoints.psd-grid.parquet")
# #MGL_grid = psd_hourly[grid$cruise == "MGL1704", ]
#
# psd_hourly <- arrow::read_parquet("2022-11-01-diam-chl-pe-Qc-PSD-combined-popcycle-4.15.0-withoutBoundaryPoints.psd-hourly.parquet")
View(psd_hourly)
psd_hourly <- psd_hourly %>%
mutate(diam = grid$diam[psd_hourly$diam_coord],
Qc = grid$Qc[psd_hourly$Qc_coord],
chl_small = grid$chl_small[psd_hourly$chl_small_coord],
pe = grid$pe[psd_hourly$pe_coord]) %>%
select(-c(diam_coord, Qc_coord, chl_small_coord, pe_coord))
MGL1704_data <- psd_hourly[psd_hourly$cruise == "MGL1704", ]
MGL1704_y_list <- MGL1704_data %>%
group_by(date) %>%
group_map(~ .x) %>%
map(~ .x %>% select(diam, chl_small, pe))
## end JB
MGL1704_y_list <- map(MGL1704_y_list, as.matrix)
MGL1704_biomass_list <- MGL1704_data %>%
group_by(date) %>%
group_map(~ .x) %>%
map(~ .x %>% select(Qc_sum_per_uL))
class(MGL1704_y_list[[1]])
# Convert each data frame to a numeric vector
MGL1704_biomass_list <- map(MGL1704_biomass_list, function(df) {
# Extract the column from the data frame and convert it to numeric
numeric_vector <- as.numeric(unlist(df))
return(numeric_vector)
})
fit_5 <- kernel_em(MGL1704_y_list, 5, 30, 30, 10, biomass = MGL1704_biomass_list)
plot_data_and_model(MGL1704_y_list, fit_5$zest, fit_5$mu, dim = 1, show_data = FALSE)
plot_pi(fit_5$pi)
plot_biomass(MGL1704_biomass_list, fit_5$resp)
plot_data_and_model(MGL1704_y_list, fit_5$zest, fit_5$mu, dim = 2, show_data = FALSE)
fit_5$mu[, 2, 1]
fit_5$mu[, 2, 2]
fit_5$mu[, 2, 3]
setwd("C:/Desktop Folders/USC/Academic/Research/flowkernel-project")
litr::render("create-flowkernel/index.Rmd")
fs::dir_copy("create-flowkernel/_book", "docs/create", overwrite = TRUE)
fs::dir_delete("create-flowkernel/_book")
fs::dir_copy("create-flowkernel/flowkernel", "flowkernel", overwrite = TRUE)
fs::dir_delete("create-flowkernel/flowkernel")
fs::dir_delete("create-flowkernel/_main_files/")
litr::render("create-flowkernel/index.Rmd")
litr::render("create-flowkernel/index.Rmd")
ls
dir()
dir("create-flowkernel")
litr::render("create-flowkernel/index.Rmd")
#' @param mu_function a function that maps a vector of times to a T-by-K-by-d
#' array of means
#' @param Sigma_function a function that maps a vector of times to a
#' T-K-by-d-by-d array of covariance matrices
#' @param pi_function a function that maps a vector of times to a T-by-K vector
#' of probabilities
#' @param num_points a T vector of integers giving the number of points n_t to
#' generate at each time point t.
#' @param start_date the starting date for the generated data. This is created for when running a test requires dates.
#' @export
generate_smooth_gauss_mix <- function(mu_function,
Sigma_function,
pi_function,
num_points,
start_date = "2017-05-28 21:00:00 UTC") {
times <- seq_along(num_points)
mu <- mu_function(times)
Sigma <- Sigma_function(times)
pi <- pi_function(times)
K <- ncol(pi) # number of components
d <- dim(mu)[3]
dimnames(mu) <- list(NULL, paste0("cluster", 1:K), NULL)
z <- list() # z[[t]][i] = class of point i at time t
y <- list() # y[[t]][i,] = d-vector of point i at time t
biomass <- list() # biomass[[t]] = biomass of particles in each bin at time t
for (t in times) {
z[[t]] <- apply(stats::rmultinom(num_points[t], 1, pi[t, ]) == 1, 2, which)
y[[t]] <- matrix(NA, num_points[t], d)
biomass[[t]] <- abs(stats::rnorm(num_points[t], mean = 0.01, sd = sqrt(0.001)))
for (k in 1:K) {
ii <- z[[t]] == k # index of points in component k at time t
if (sum(ii) == 0) next
if (d == 1)
y[[t]][ii, ] <- stats::rnorm(n = sum(ii),
mean = mu[t, k, ],
sd = Sigma[t, k, , ])
else
y[[t]][ii, ] <- mvtnorm::rmvnorm(n = sum(ii),
mean = mu[t, k, ],
sigma = Sigma[t, k, , ])
}
}
# Generate dates spaced 1 hour apart
start_time <- as.POSIXct(start_date, tz = "UTC")
dates <- seq(start_time, by = "hour", length.out = length(y))
list(y = y, z = z, mu = mu, Sigma = Sigma, pi = pi, biomass = biomass, dates = dates)
}
set.seed(123)
d <- 1; K <- 2; ntimes <- 200
ex1 <- list(
mu_function = function(times) {
mu <- array(NA, c(ntimes, K, d))
mu[, , 1] <- cbind(sin(2 * pi * times / 30), 2)
mu
},
Sigma_function = function(times) {
Sigma <- array(NA, c(ntimes, K, 1, 1))
Sigma[, , 1, 1] <- 0.25
Sigma
},
pi_function = function(times) {
pi1 <- seq(0.2, 0.8, length=length(times))
cbind(pi1, 1 - pi1)
},
num_points = rep(40, ntimes)
)
ex1$dat <- generate_smooth_gauss_mix(ex1$mu_function,
ex1$Sigma_function,
ex1$pi_function,
ex1$num_points)
d <- 3; K <- 4; ntimes <- 200
ex2 = list(
mu_function = function(times) {
mu <- array(NA, c(ntimes, K, d))
mu[, , 1] <- cbind(0.5*cos(2 * pi * times / 30), 0.3*sin(2 * pi * times / 30), sin(2 * pi * times / 30), -3)
mu[, , 2] = cbind (0.3*sin(2 * pi * times / 30), 2, -1, 0.6*cos(2 * pi * times / 30))
mu[, , 3] = cbind(2, 0.7*cos(2 * pi * times / 30), 0.4*sin(2 * pi * times / 30), 1)
mu
},
Sigma_function = function(times) {
Sigma <- array(NA, c(ntimes, K, d, d))
Sigma[, , 1, 1] <- 0.10
Sigma[, , 1, 2] <- 0
Sigma[, , 1, 3] <- 0
Sigma[, , 2, 1] <- 0
Sigma[, , 2, 2] <- 0.10
Sigma[, , 2, 3] <- 0
Sigma[, , 3, 1] <- 0
Sigma[, , 3, 2] <- 0
Sigma[, , 3, 3] <- 0.10
Sigma
},
pi_function = function(times) {
pi1 <- seq(0.2, 0.3, length=length(times))
cbind(pi1, pi1, 2*pi1/3, 1- (2*pi1 + 2*pi1/3))
},
num_points = rep(150, ntimes)
)
ex2$dat <- generate_smooth_gauss_mix(ex2$mu_function,
ex2$Sigma_function,
ex2$pi_function,
ex2$num_points)
#' @param hmu bandwidth for mu parameter
#' @param hSigma bandwidth for Sigma parameter
#' @param hpi bandwidth for pi parameter
#' @param dates A vector of POSIXct (or POSIXt) date-time objects of length T, where each element represents the timestamp corresponding to the t-th observation. If provided, these dates are converted to numeric values and rescaled (by dividing by 3600) for time-based smoothing operations. If omitted, the function defaults to using the sequence of observation indices.
#' @param num_iter number of iterations of EM to perform
#' @param biomass list of length T, where each element `biomass[[t]]` is a
#' numeric vector of length n_t containing the biomass (or count) of particles
#' in each bin
#' @param initial_fit a list of starting values for the parameters, responsibilities, and estimated cluster assignments
#' @export
kernel_em <- function (y, K, hmu, hSigma, hpi, dates = NULL, num_iter = 10,
biomass = default_biomass(y),
initial_fit = init_const(y, K, 50, 50)) {
num_times <- length(y)
d <- ncol(y[[1]])
mu <- initial_fit$mu
Sigma <- initial_fit$Sigma
pi <- initial_fit$pi
if (!is.null(dates)){
numeric_dates <- as.numeric(dates)
rescaled_dates <- (numeric_dates - min(numeric_dates)) / 3600
}
for (l in seq(num_iter)) {
<<E-step>>
#' @param fit A list returned by `kernel_em` containing the fitted parameters (mu, Sigma, pi, and responsibilities) from the training data.
#' @param test_dates A vector of POSIXct (or POSIXt) date-time objects of length T_test representing the time points at which predictions are required.
#' @param train_dates A vector of POSIXct (or POSIXt) date-time objects of length T_train representing the time points corresponding to the training data.
#' @param train_data A list of length T_train, where each element is a matrix containing the observations for the training set at the corresponding time point.
#' @param train_biomass A list of length T_train, where each element is a numeric vector of length n_t containing the biomass (or count) for the observations in the training set.
#' @param hmu Bandwidth parameter for smoothing the mu estimates.
#' @param hSigma Bandwidth parameter for smoothing the Sigma estimates.
#' @param hpi Bandwidth parameter for smoothing the pi (mixing proportions) estimates.
#' @param scale_dates Logical flag indicating whether to rescale the date-time objects. If TRUE, the dates are converted to numeric values and rescaled (by subtracting the minimum training date and dividing by 3600); if FALSE, the dates are used directly. Defaults to TRUE.
#' @export
kernel_em_predict <- function(fit,
test_dates,
train_dates,
train_data,
train_biomass,
hmu,
hSigma,
hpi,
scale_dates = TRUE) {
num_test <- length(test_dates)
num_train <- length(train_dates)
K <- ncol(fit$pi)
d <- dim(fit$mu)[3]
mu <- array(NA, c(num_test, K, d))
Sigma <- array(NA, c(num_test, K, d, d))
pi <- matrix(NA, num_test, K)
# Convert to numeric
numeric_train_dates <- as.numeric(train_dates)
numeric_test_dates  <- as.numeric(test_dates)
if (scale_dates) {
# If we are given real dates, then we do the original rescaling
min_date <- min(numeric_train_dates)
rescaled_train_dates <- (numeric_train_dates - min_date) / 3600
rescaled_test_dates  <- (numeric_test_dates  - min_date) / 3600
} else {
# If dates are already in hours or integer indices, treat them directly
rescaled_train_dates <- numeric_train_dates
rescaled_test_dates  <- numeric_test_dates
}
# 1) Predict pi (mixing proportions) ----------------------------------------
# Weighted responsibilities (weights = biomass)
resp_weighted <- purrr::map2(train_biomass, fit$resp, ~ .y * .x)
# Sum over each mixture component across all observations
resp_sum <- purrr::map(resp_weighted, ~ colSums(.x)) %>%
unlist() %>%
matrix(ncol = K, byrow = TRUE)
# Smooth those sums across time
resp_sum_smooth <- apply(
resp_sum, 2, function(x)
stats::ksmooth(rescaled_train_dates, x, kernel = "normal", bandwidth = hpi,
x.points = rescaled_test_dates)$y
)
# Normalize to get pi
pi <- resp_sum_smooth / rowSums(resp_sum_smooth)
# 2) M-step for mu ---------------------------------------------------------
# Weighted sum of y's for each mixture component
y_sum <- purrr::map2(resp_weighted, train_data, ~ crossprod(.x, .y)) %>%
unlist() %>%
array(c(K, d, num_train)) %>%
aperm(c(3,1,2))
# Smooth each dimension across time
y_sum_smoothed <- apply(
y_sum, 2:3, function(x)
stats::ksmooth(rescaled_train_dates, x, kernel = "normal", bandwidth = hmu,
x.points = rescaled_test_dates)$y
)
# Smooth the sum of responsibilities (again, but with bandwidth hmu)
resp_sum_smooth_mu <- apply(
resp_sum, 2, function(x)
stats::ksmooth(rescaled_train_dates, x, kernel = "normal", bandwidth = hmu,
x.points = rescaled_test_dates)$y
)
# Combine smoothed sums to get mu
for (j in seq(d)) {
mu[, , j] <- y_sum_smoothed[, , j] / resp_sum_smooth_mu
}
# 3) M-step for Sigma ------------------------------------------------------
mat_sum <- array(NA, c(num_train, K, d, d))
for (tt in seq(num_train)) {
# Prepare a matrix for each observation's difference from mu
yy <- matrix(NA, nrow(train_data[[tt]]), d)
for (k_idx in seq(K)) {
for (dd in seq(d)) {
yy[, dd] <- train_data[[tt]][, dd] - fit$mu[tt, k_idx, dd]
}
mat_sum[tt, k_idx, , ] <- crossprod(yy, yy * resp_weighted[[tt]][, k_idx])
}
}
# Smooth Sigma the same way
mat_sum_smoothed <- apply(
mat_sum, 2:4, function(x)
stats::ksmooth(rescaled_train_dates, x, kernel = "normal", bandwidth = hSigma,
x.points = rescaled_test_dates)$y
)
resp_sum_smooth_Sigma <- apply(
resp_sum, 2, function(x)
stats::ksmooth(rescaled_train_dates, x, kernel = "normal", bandwidth = hSigma,
x.points = rescaled_test_dates)$y
)
for (j in seq(d)) {
for (l in seq(d)) {
Sigma[, , j, l] <- mat_sum_smoothed[, , j, l] / resp_sum_smooth_Sigma
}
}
list(mu = mu, Sigma = Sigma, pi = pi)
}
fs::dir_copy("create-flowkernel/_book", "docs/create", overwrite = TRUE)
fs::dir_delete("create-flowkernel/_book")
fs::dir_copy("create-flowkernel/flowkernel", "flowkernel", overwrite = TRUE)
fs::dir_delete("create-flowkernel/flowkernel")
fs::dir_delete("create-flowkernel/_main_files/")
fs::dir_copy("create-flowkernel/_book", "docs/create", overwrite = TRUE)
fs::dir_delete("create-flowkernel/_book")
fs::dir_copy("create-flowkernel/flowkernel", "flowkernel", overwrite = TRUE)
fs::dir_delete("create-flowkernel/flowkernel")
fs::dir_delete("create-flowkernel/_main_files/")
fs::dir_copy("create-flowkernel/_book", "docs/create", overwrite = TRUE)
fs::dir_delete("create-flowkernel/_book")
fs::dir_copy("create-flowkernel/flowkernel", "flowkernel", overwrite = TRUE)
fs::dir_delete("create-flowkernel/flowkernel")
fs::dir_delete("create-flowkernel/_main_files/")
fs::dir_copy("create-flowkernel/_book", "docs/create", overwrite = TRUE)
fs::dir_delete("create-flowkernel/_book")
fs::dir_copy("create-flowkernel/flowkernel", "flowkernel", overwrite = TRUE)
fs::dir_delete("create-flowkernel/flowkernel")
fs::dir_delete("create-flowkernel/_main_files/")