Commit 4d0a1132 authored by luroth's avatar luroth
Browse files

spline fitting optimized

parent ab2467ff
# Working directory with temperature data # Working directory with temperature data
path_home <- 'E:/Scripts/htfp_data_processing' path_home <- 'E:/Scripts/htfp_data_processing'
path_simulation <- 'E:/Simulation/Runs' path_simulation <- 'E:/Simulation/Runs_Posterior'
setwd(path_home) setwd(path_home)
...@@ -16,6 +16,7 @@ library(zoo) ...@@ -16,6 +16,7 @@ library(zoo)
library(plyr) library(plyr)
library(stringr) library(stringr)
library(SpATS) library(SpATS)
library(fs)
library(foreach) library(foreach)
library(dplyr) library(dplyr)
...@@ -27,9 +28,10 @@ source("R/Model/Dose_response.R") ...@@ -27,9 +28,10 @@ source("R/Model/Dose_response.R")
source("R/Model/FitSpATS.R") source("R/Model/FitSpATS.R")
source("R/Model/Graphs.R") source("R/Model/Graphs.R")
# Missing from euler: 1, 3, 8
start_run <- 1 start_run <- 1
max_runs <- 1 max_runs <- 500
number_of_cpus <- 5 number_of_cpus <- 29
sigma_error <- 10 sigma_error <- 10
...@@ -103,15 +105,15 @@ df_designs <- df_designs %>% mutate(plot.discrete_x = plot.row * if_else(plot.re ...@@ -103,15 +105,15 @@ df_designs <- df_designs %>% mutate(plot.discrete_x = plot.row * if_else(plot.re
df_genotypes <- read_csv('Simulation/genotypes.csv') df_genotypes <- read_csv('Simulation/genotypes.csv')
df_temp <- read_csv('Simulation/covariate_temp.csv') df_temp <- read_csv('Simulation/covariate_temp.csv')
run <- 1 # run <- 1
for(run in start_run:max_runs) { # for(run in start_run:max_runs) {
# cl <- parallel::makeCluster(number_of_cpus) cl <- parallel::makeCluster(number_of_cpus)
# doParallel::registerDoParallel(cl) doParallel::registerDoParallel(cl)
# foreach(run = start_run:max_runs, .verbose = TRUE, foreach(run = start_run:max_runs, .verbose = TRUE,
# .packages = c("readr", "tidyr", "purrr", "ggplot2", "lubridate", "gridExtra", "plyr", "stringr", "SpATS", "dplyr", "scam", "MASS") .packages = c("fs", "readr", "tidyr", "purrr", "ggplot2", "lubridate", "gridExtra", "plyr", "stringr", "SpATS", "dplyr", "scam", "MASS")
# ) %dopar% { ) %dopar% {
# Validation data # Validation data
...@@ -156,11 +158,12 @@ for(run in start_run:max_runs) { ...@@ -156,11 +158,12 @@ for(run in start_run:max_runs) {
if(length(files)==0) { if(length(files)==0) {
#try({ #try({
fs::file_touch(paste0(path_simulation, "/", run, "/", set, "_MILESTONE_start_plot.txt"))
df_values_for_fit_orig_ <- df_values_for_fit_orig_ %>% df_values_for_fit_orig_ <- df_values_for_fit_orig_ %>%
filter(format(timestamp, "%m%d") %in% measurement_dates_set) filter(format(timestamp, "%m%d") %in% measurement_dates_set)
# Initial correction to test van Eeujwick 2018: # Initial correction to test van Eeujwick 2018:
df_BLUEs_for_fit_orig_ <- df_values_for_fit_orig_ %>% df_BLUEs_for_fit_orig_ <- df_values_for_fit_orig_ %>%
mutate(se=1, year_site.UID_ = year_site.UID) %>% mutate(se=1, year_site.UID_ = year_site.UID) %>%
group_by(year_site.UID_, timestamp) %>% group_by(year_site.UID_, timestamp) %>%
...@@ -170,6 +173,7 @@ for(run in start_run:max_runs) { ...@@ -170,6 +173,7 @@ for(run in start_run:max_runs) {
df_BLUEs_for_fit_orig_ <- df_BLUEs_for_fit_orig_ %>% df_BLUEs_for_fit_orig_ <- df_BLUEs_for_fit_orig_ %>%
mutate(year_site.UID = paste0(year_site.UID_, "_corrected"), mutate(year_site.UID = paste0(year_site.UID_, "_corrected"),
value = BLUE, value_se = BLUE_SE, plot.UID = paste0("BLUE_", genotype.id, year_site.UID)) value = BLUE, value_se = BLUE_SE, plot.UID = paste0("BLUE_", genotype.id, year_site.UID))
file_touch(paste0(path_simulation, "/", run, "/", set, "_MILESTONE_reverse_SpATS.txt"))
df_values_for_fit_orig_$value_se <- 1 df_values_for_fit_orig_$value_se <- 1
df_values_for_fit_orig <- bind_rows(df_values_for_fit_orig_, df_BLUEs_for_fit_orig_) df_values_for_fit_orig <- bind_rows(df_values_for_fit_orig_, df_BLUEs_for_fit_orig_)
...@@ -202,6 +206,8 @@ for(run in start_run:max_runs) { ...@@ -202,6 +206,8 @@ for(run in start_run:max_runs) {
~fit_scam_spline_weights(.$timestamp, .$value, .$value_se))) ~fit_scam_spline_weights(.$timestamp, .$value, .$value_se)))
df_spline_model <- df_spline_model %>% dplyr::select(-data) df_spline_model <- df_spline_model %>% dplyr::select(-data)
file_touch(paste0(path_simulation, "/", run, "/", set, "_MILESTONE_fit_spline.txt"))
# Predict value with spline model # Predict value with spline model
time_interval <- 60*60*12 time_interval <- 60*60*12
...@@ -222,6 +228,8 @@ for(run in start_run:max_runs) { ...@@ -222,6 +228,8 @@ for(run in start_run:max_runs) {
predict_se = predict_scam_spline(.$spline_model[[1]], list(x=.$prediction_timepoint[[1]]), se = T), predict_se = predict_scam_spline(.$spline_model[[1]], list(x=.$prediction_timepoint[[1]]), se = T),
predict_deriv = predict_scam_spline(.$spline_model[[1]], list(x=.$prediction_timepoint[[1]]), deriv = T), predict_deriv = predict_scam_spline(.$spline_model[[1]], list(x=.$prediction_timepoint[[1]]), deriv = T),
timestamp = .$prediction_timepoint[[1]])) timestamp = .$prediction_timepoint[[1]]))
file_touch(paste0(path_simulation, "/", run, "/", set, "_MILESTONE_predict_spline.txt"))
if (run == 1 & set == "set3d") { if (run == 1 & set == "set3d") {
df_spline_predicts_export <- df_spline_predicts %>% df_spline_predicts_export <- df_spline_predicts %>%
...@@ -234,6 +242,8 @@ for(run in start_run:max_runs) { ...@@ -234,6 +242,8 @@ for(run in start_run:max_runs) {
print("Find start/stop") print("Find start/stop")
df_growth_phase_predicts <- df_spline_predicts %>% ungroup() %>% group_by(plot.UID) %>% df_growth_phase_predicts <- df_spline_predicts %>% ungroup() %>% group_by(plot.UID) %>%
mutate(growth_phase_params = map(spline_predicts, find_start_stop_growth_phase_posterior)) mutate(growth_phase_params = map(spline_predicts, find_start_stop_growth_phase_posterior))
file_touch(paste0(path_simulation, "/", run, "/", set, "_MILESTONE_extract_traits.txt"))
# Extract predicted values # Extract predicted values
df_growth_phase_predicts <- df_growth_phase_predicts %>% df_growth_phase_predicts <- df_growth_phase_predicts %>%
......
Supports Markdown
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