Commit 8ac19439 authored by luroth's avatar luroth
Browse files

spline fitting optimized

parent dac1cd47
......@@ -32,6 +32,7 @@ library(scam)
coefs
}
fit_scam_spline <- function(x, y, k = NA, bs = "mpi", label = NULL, optimizer = "bfgs") {
if (!is.null(label)) { print(label) }
if (is.na(k)) k <- round(length(x) * 3 / 4)
......@@ -39,32 +40,51 @@ fit_scam_spline <- function(x, y, k = NA, bs = "mpi", label = NULL, optimizer =
spline <- NULL
try(
system.time(spline <- scam(y ~ s(as.numeric(x), k = k, bs = bs), optimizer = optimizer))
spline <- R.utils::withTimeout(scam(y ~ s(as.numeric(x), k = k, bs = bs), optimizer = optimizer), timeout=10)
)
if (is.null(spline)) {
print("decreasing k")
spline <- scam(y ~ s(as.numeric(x), k = k - 1, bs = bs))
try(
spline <- R.utils::withTimeout(scam(y ~ s(as.numeric(x), k = k - 1, bs = bs)), timeout=10)
)
if (is.null(spline)) {
print("2. time decreasing k")
try(
spline <- R.utils::withTimeout(scam(y ~ s(as.numeric(x), k = k - 2, bs = bs)), timeout=10)
)
}
}
return(spline)
}
fit_scam_spline_weights <- function(x, y, w, k = NA, bs = "mpi", label = NULL, optimizer = "bfgs") {
w <- 1 / (w^2)
if (!is.null(label)) { print(label) }
if (is.na(k)) k <- round(length(x) * 3 / 4)
if (k > 20) k <- 20
spline <- NULL
try(
system.time(spline <- scam(y ~ s(as.numeric(x), k = k, bs = bs), optimizer = optimizer, weights=w))
spline <- R.utils::withTimeout(scam(y ~ s(as.numeric(x), k = k, bs = bs), optimizer = optimizer, weights = w), timeout=1)
)
if (is.null(spline)) {
print("decreasing k")
spline <- scam(y ~ s(as.numeric(x), k = k - 1, bs = bs), weights=w)
try(
spline <- R.utils::withTimeout(scam(y ~ s(as.numeric(x), k = k - 1, bs = bs), weights = w), timeout=2)
)
if (is.null(spline)) {
print("2. time decreasing k")
try(
spline <- R.utils::withTimeout(scam(y ~ s(as.numeric(x), k = k - 2, bs = bs)), timeout=4)
)
}
}
return(spline)
}
predict_scam_spline <- function(spline, x_, deriv = NULL, se = FALSE)
{
......@@ -123,7 +143,7 @@ find_start_stop_growth_phase <- function(df, text = NA, threshold_start = 1 / 4,
# Extract maximum growth phase
max_growth <- max(df$predict_deriv)
if (max_growth == 0) {
print("Spline has zero growth!")
#print("Spline has zero growth!")
return(tibble(predict_start_growth = NA,
predict_start_growth_se = NA,
predict_stop_growth = NA,
......@@ -233,6 +253,7 @@ find_start_stop_growth_phase_posterior <- function(df, text = NA, threshold_star
dplyr::select(ends_with("sepost"))
predict <- bind_cols(predict, df_predict_posterior_se)
print(text)
return(predict)
}
\ No newline at end of file
# Working directory with temperature data
path_home <- 'C:/Users/luroth/PycharmProjects/htfp_wheat_canopy_height_processing'
path_simulation <- 'E:/Simulation/Runs_Posterior'
setwd(path_home)
# Libraries to use
library(readr)
library(tidyr)
library(purrr)
library(ggplot2)
library(lubridate)
library(gridExtra)
library(plyr)
library(stringr)
library(xtable)
library(dplyr)
source("R/Model/Graphs.R")
my_label_parsed <- function (variable, value) {
llply(as.character(value), function(x) parse(text = str_replace(x, " ", "~~")))
}
start_run <- 1
max_runs <- 10
df_correlation_plot_predict_true <- data.frame()
df_correlation_genotype_predict_true <- data.frame()
read_files_runs <- function(add_to_df, pattern, col_types =NULL) {
files <- Sys.glob(paste0(path_simulation, "/", run, "/", pattern))
for(file in files) {
df <- read_csv(file, col_types = col_types)
df$run <- run
add_to_df <- bind_rows(add_to_df, df)
}
return(add_to_df)
}
for(run in start_run:max_runs) {
files <- Sys.glob(paste0(path_simulation, "/", run, "/*_plot_true_versus_predict.csv"))
if(length(files) >= 1) {
print(run)
df_correlation_plot_predict_true <- read_files_runs(
df_correlation_plot_predict_true, "*_plot_true_versus_predict.csv",
col_types = cols(
.default = col_double(),
plot.UID = col_character(),
parameter = col_character(),
year_site.UID = col_character(),
model = col_character(),
site.label = col_character(),
label = col_character(),
set = col_character()
))
}
}
df_correlation_plot_predict_true <- df_correlation_plot_predict_true %>% filter(set == "set3d")
df_correlation_plot_predict_true_ <- df_correlation_plot_predict_true %>% filter(model %in% c("spline"))
df_correlation_plot_predict_true_ <- df_correlation_plot_predict_true_ %>%
mutate(
model = case_when(
model == "bplm_simple" ~ "Break-point T[mean]",
model == "percentile" ~ "Percentiles",
TRUE ~ "P-spline/QMER"
),
parameter = case_when(
parameter == "bplm_Asym" ~"G",
parameter == "bplm_c0" ~"c[0]",
parameter == "bplm_cOpt" ~"c[Opt]",
parameter == "final_height" ~"PH[max]",
parameter == "start_growth" ~"t[PH[start]]",
parameter == "stop_growth" ~"t[PH[stop]]",
)
)
ggplot(data=df_correlation_plot_predict_true_, aes(x=se, y=sepost)) +
geom_point() +
facet_wrap(parameter ~ year_site.UID, scales="free")
df_correlation_plot_predict_true_ %>% group_by(parameter, year_site.UID) %>%
summarise(corr = cor(se, sepost, use="pairwise.complete.obs"))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment