Skip to content

Commit

Permalink
Merge c7c3165 into ce19d54
Browse files Browse the repository at this point in the history
  • Loading branch information
hdugan committed Aug 14, 2020
2 parents ce19d54 + c7c3165 commit 174c0d4
Show file tree
Hide file tree
Showing 61 changed files with 658 additions and 122 deletions.
3 changes: 3 additions & 0 deletions .Rbuildignore
Expand Up @@ -4,3 +4,6 @@
appveyor.yml
^README\.Rmd$
^README-.*\.png$
^doc$
^Meta$
glm4.nml
11 changes: 8 additions & 3 deletions .travis.yml
Expand Up @@ -18,8 +18,13 @@ r:
- release
- devel

r_build_args: --no-build-vignettes --no-manual
r_check_args: --no-build-vignettes --no-manual --as-cran
r_build_args: --no-manual
r_check_args: --no-manual --as-cran

after_success:
- Rscript -e 'covr::coveralls()'

branches:
only:
- /.*/

13 changes: 7 additions & 6 deletions DESCRIPTION
Expand Up @@ -37,15 +37,16 @@ Imports:
ncdf4,
readr,
reshape2,
tidyr,
tools
tidyr (>= 1.0.0),
tools,
rlang
Suggests:
testthat,
vdiffr,
knitr,
rmarkdown
VignetteBuilder: knitr
BuildVignettes: true
LazyLoad: yes
VignetteBuilder:
knitr
LazyData: yes
RoxygenNote: 7.0.2
RoxygenNote: 7.1.1
Encoding: UTF-8
8 changes: 8 additions & 0 deletions NAMESPACE
Expand Up @@ -55,6 +55,7 @@ importFrom(GLM3r,nml_template_path)
importFrom(GLM3r,run_glm)
importFrom(akima,interp)
importFrom(akima,interp2xyz)
importFrom(graphics,legend)
importFrom(lazyeval,lazy_dots)
importFrom(lazyeval,lazy_eval)
importFrom(methods,getPackageName)
Expand All @@ -70,8 +71,15 @@ importFrom(rLakeAnalyzer,layer.temperature)
importFrom(readr,parse_number)
importFrom(readr,read_csv)
importFrom(reshape2,melt)
importFrom(rlang,.data)
importFrom(stats,approx)
importFrom(stats,cor.test)
importFrom(stats,dnorm)
importFrom(stats,median)
importFrom(stats,qqnorm)
importFrom(stats,sd)
importFrom(tidyr,gather)
importFrom(tidyr,pivot_longer)
importFrom(tools,file_ext)
importFrom(utils,capture.output)
importFrom(utils,lsf.str)
Expand Down
50 changes: 36 additions & 14 deletions R/calib_helpers.R
@@ -1,27 +1,29 @@
#'@param var character variable string
#'@noRd
#'@import adagio
#'@importFrom utils write.csv
calib_GLM <- function(var, ub, lb, init.val, obs, method, glmcmd,
metric, target.fit, target.iter, nml_file, path, scaling, verbose){
path <<- path

metric, target.fit, target.iter, nml_file, path, scaling,
verbose, pars){
if (method == 'CMA-ES'){
glmOPT <- pureCMAES(par = init.val, fun = glmFUN, lower = rep(0,length(init.val)),
upper = rep(10,length(init.val)),
sigma = 0.5,
stopfitness = target.fit,
stopeval = target.iter,
glmcmd = glmcmd, nml_file = nml_file, var = var,
scaling = scaling, metric = metric, verbose = verbose)
scaling = scaling, metric = metric, verbose = verbose,
ub = ub, lb = lb, pars = pars, path = path, obs = obs)
} else if (method == 'Nelder-Mead'){
glmOPT <- neldermeadb(fn = glmFUN, init.val, lower = rep(0,length(init.val)),
glmOPT <- neldermeadb(fn = glmFUN, init.val, ub = ub, lb = lb, lower = rep(0,length(init.val)),
upper = rep(10,length(init.val)),
adapt = TRUE,
tol = 1e-10,
maxfeval = target.iter, glmcmd = glmcmd, nml_file = nml_file, var = var,
scaling = scaling, metric = metric, verbose = verbose )
}
maxfeval = target.iter, glmcmd = glmcmd, nml_file = nml_file, var = var, ub = ub, lb = lb, path = path, obs = obs)
}

glmFUN(p = glmOPT$xmin, nml_file = nml_file, glmcmd = glmcmd, var = var, scaling, metric, verbose, ub = ub, lb = lb, path = path, pars = pars, obs = obs)

glmFUN(p = glmOPT$xmin,nml_file = nml_file,glmcmd = glmcmd, var = var,scaling,metric,verbose)
# glmFUN(glmOPT$xmin, glmcmd, scaling, metric, verbose)
calib <- read.csv(paste0(path,'/calib_results_',metric,'_',var,'.csv'))
eval(parse(text = paste0('best_par <- calib[which.min(calib$',metric,'),]')))
Expand Down Expand Up @@ -88,7 +90,9 @@ calib_GLM <- function(var, ub, lb, init.val, obs, method, glmcmd,
#Run GLM
run_glmcmd(glmcmd, path, verbose)

g1 <- diag.plots(mod2obs(paste0(path,'/output/output.nc'), obs, reference = 'surface', var), obs)
g1 <- diag.plots(
mod2obs(paste0(path,'/output/output.nc'), obs, reference = 'surface', var),
obs, ggplot=F)

ggsave(filename = paste0(path,'/diagnostics_',method,'_',var,'.png'), plot = g1,
dpi = 300,width = 384,height = 216, units = 'mm')
Expand All @@ -97,7 +101,9 @@ calib_GLM <- function(var, ub, lb, init.val, obs, method, glmcmd,
}


glmFUN <- function(p, glmcmd, nml_file, var, scaling, metric, verbose){
#'@param var character
#'@noRd
glmFUN <- function(p, glmcmd, nml_file, var, scaling, metric, verbose, ub, lb, pars, path, obs){
#Catch non-numeric arguments
if(!is.numeric(p)){
p = values.optim
Expand All @@ -107,7 +113,6 @@ glmFUN <- function(p, glmcmd, nml_file, var, scaling, metric, verbose){
}

eg_nml <- read_nml(nml_file)

for(i in 1:length(pars[!duplicated(pars)])){
if (any(pars[!duplicated(pars)][i] == pars[duplicated(pars)])){
eg_nml <- set_nml(eg_nml, pars[!duplicated(pars)][i],
Expand All @@ -121,6 +126,7 @@ glmFUN <- function(p, glmcmd, nml_file, var, scaling, metric, verbose){
write_nml(eg_nml, file = write_path)

error <- try(run_glmcmd(glmcmd, path, verbose))

while (error != 0){
error <- try(run_glmcmd(glmcmd, path, verbose))
}
Expand Down Expand Up @@ -173,6 +179,8 @@ get_rmse <- function(mods, obs){
return(rmse)
}

#'@param var character
#'@noRd
#'@importFrom reshape2 melt
mod2obs <- function(mod_nc, obs, reference = 'surface', var){
deps = unique(obs[,2])
Expand Down Expand Up @@ -217,6 +225,7 @@ match.tstep <- function(df1, df2){
}

#'@import graphics
#'@importFrom stats dnorm median qqnorm sd
diag.plots <- function(mod, obs, ggplot = T){
stats = sum_stat(mod, obs, depth = T)
if(max(mod[,2]) >= 0){ #Makes depths negative
Expand Down Expand Up @@ -294,7 +303,7 @@ diag.plots <- function(mod, obs, ggplot = T){
gp=grid::gpar(col="black", fontsize=10)))

#Plots
p1 <-ggplot(mod, aes(x = res)) +
p1 <- ggplot(mod, aes(x = .data$res)) +
geom_histogram(fill = "lightblue4", colour = 'black', breaks = seq(min.res, max.res, bw)) +
stat_function(
fun = function(x, mean, sd, n, bw){
Expand Down Expand Up @@ -355,7 +364,7 @@ diag.plots <- function(mod, obs, ggplot = T){
theme_bw()
p5 <- p5 + annotation_custom(grob2) + annotation_custom(grob3)

p6 <- ggplot(mod, aes(sample = res))+
p6 <- ggplot(mod, aes(sample = .data$res))+
stat_qq()+
geom_abline(slope = 1, intercept = 0, size =1, linetype = 'dashed')+
xlab('Sample Quantiles')+
Expand Down Expand Up @@ -387,6 +396,7 @@ get_nse <- function(x, y){
}

# gotmtools.R
#' @importFrom stats cor.test
sum_stat <- function(mod, obs, depth =F,na.rm =T, depth.range =NULL){
if(depth == T){
if(!is.null(depth.range)){
Expand Down Expand Up @@ -440,3 +450,15 @@ checkHourFormat <- function(timest, hourst){
}
return(as.POSIXct(newTimest))
}

#' Add text to the corner of a plot
#'
#' Add text to a location within a plot. <https://github.com/aemon-j/gotmtools/blob/master/R/Corner_text.R>
#'
#' @param text text to be added to the plot as a language object
#' @param location location within the plot for the text to be placed
#' @return data
#' @importFrom graphics legend
Corner_text <- function(text, location="topleft"){
legend(location,legend=text, bty ="n", pch=NA)
}
24 changes: 13 additions & 11 deletions R/calibrate_sim.R
Expand Up @@ -40,6 +40,7 @@
#'output = file.path(sim_folder, 'output/output.nc')
#'
#'var = 'temp' # variable to apply the calibration procedure
#'\dontrun{
#'calibrate_sim(var = var, path = sim_folder, field_file = field_file,
#' nml_file = nml_file, calib_setup = calib_setup,
#' glmcmd = NULL,
Expand All @@ -49,6 +50,7 @@
#' metric = 'RMSE',plotting = FALSE,
#' target.fit = 1.5,
#' target.iter = 50, output = output)
#' }
#'@import adagio
#'@import GLM3r
#'@import ggplot2
Expand Down Expand Up @@ -91,11 +93,10 @@ calibrate_sim <- function(var = 'temp',
if (is.null(calib_setup)){
calib_setup <- get_calib_setup()
}
pars <<- as.character(calib_setup$pars)
ub <<- calib_setup$ub
lb <<- calib_setup$lb
variable <<- var

pars <- as.character(calib_setup$pars)
ub <- calib_setup$ub
lb <- calib_setup$lb
variable <- var

if (scaling){
init.val <- (calib_setup$x0 - lb) *10 /(ub-lb)
Expand All @@ -107,15 +108,16 @@ calibrate_sim <- function(var = 'temp',
write_nml(nml,nml_file)
}

# path <<- path
obs <<- read_field_obs(field_file)
obs <- read_field_obs(field_file)
calib_GLM(var, ub, lb, init.val, obs, method, glmcmd,
metric, target.fit, target.iter, nml_file, path, scaling, verbose)
metric, target.fit, target.iter, nml_file, path, scaling, verbose, pars)

# loads all iterations
results <- read.csv(paste0(path,'/calib_results_RMSE_temp.csv'))
results$DateTime <- as.POSIXct(results$DateTime)
g1 <- ggplot(results, aes(nrow(results):1, RMSE)) +

g1 <- ggplot(results,
aes(nrow(results):1, .data$RMSE)) +
geom_point() +
geom_smooth(se = FALSE, method = "gam", formula = y ~ s(x), color = 'lightblue4') +
theme_bw() + xlab('Iterations') +
Expand Down Expand Up @@ -183,7 +185,7 @@ calibrate_sim <- function(var = 'temp',
print(paste('calibration:',round(temp_rmse1,2),'deg C RMSE'))
print(paste('validation:',round(temp_rmse2,2),'deg C RMSE'))
print(paste('total time period:',round(temp_rmse3,2),'deg C RMSE'))
return(print(calibrated_results))

print(calibrated_results)
return(calibrated_results)
}

2 changes: 1 addition & 1 deletion R/compare_to_field.R
Expand Up @@ -17,7 +17,7 @@
#'@author
#'Jordan S. Read
#'@examples
#' nc_file <- system.file("extdata", "output.nc", package = "glmtools")
#' nc_file <- system.file("extdata", "output/output.nc", package = "glmtools")
#' field_file <- system.file("extdata", "LakeMendota_field_data_hours.csv", package = "glmtools")
#'
#' thermo_values <- compare_to_field(nc_file, field_file,
Expand Down
2 changes: 1 addition & 1 deletion R/get_calib_init_validation.R
Expand Up @@ -9,7 +9,7 @@
#'Robert Ladwig
#'
#'@examples
#'nc_file <- system.file("extdata", "output.nc", package = "glmtools")
#'nc_file <- system.file("extdata", "output/output.nc", package = "glmtools")
#'nml_file <- system.file("extdata", "glm3.nml", package = "glmtools")
#'initvalues <- get_calib_init_validation(nml = nml_file, output = nc_file)
#'@export
Expand Down
2 changes: 1 addition & 1 deletion R/get_evaporation.R
Expand Up @@ -9,7 +9,7 @@
#'Jordan S. Read
#'@examples
#'sim_folder <- run_example_sim(verbose = FALSE)
#'nc_file <- file.path(sim_folder, 'output.nc')
#'nc_file <- file.path(sim_folder, 'output/output.nc')
#'evap <- get_evaporation(nc_file)
#'@export
get_evaporation <- function(file){
Expand Down
2 changes: 1 addition & 1 deletion R/get_ice.R
Expand Up @@ -13,7 +13,7 @@
#'Luke A. Winslow, Jordan S. Read
#'@examples
#'sim_folder <- run_example_sim(verbose = FALSE)
#'nc_file <- file.path(sim_folder, 'output.nc')
#'nc_file <- file.path(sim_folder, 'output/output.nc')
#'ice <- get_ice(nc_file)
#'ice_and_snow <- get_ice(nc_file, snow.rm = FALSE)
#'plot(ice)
Expand Down
2 changes: 1 addition & 1 deletion R/get_raw.R
Expand Up @@ -29,7 +29,7 @@
#'Luke A. Winslow, Jordan S. Read
#'@examples
#'sim_folder <- run_example_sim(verbose = FALSE)
#'nc_file <- file.path(sim_folder, 'output.nc')
#'nc_file <- file.path(sim_folder, 'output/output.nc')
#'wind <- get_raw(nc_file, 'wind') #raw wind data
#'@importFrom ncdf4 ncvar_get
#'@seealso \code{\link{sim_vars}}
Expand Down
2 changes: 1 addition & 1 deletion R/get_surface_height.R
Expand Up @@ -13,7 +13,7 @@
#'Jordan S. Read, Luke A. Winslow
#'@examples
#'sim_folder <- run_example_sim(verbose = FALSE)
#'nc_file <- file.path(sim_folder, 'output.nc')
#'nc_file <- file.path(sim_folder, 'output/output.nc')
#'surface <- get_surface_height(file = nc_file)
#'surface_w_ice <- get_surface_height(file = nc_file, ice.rm = FALSE, snow.rm = FALSE)
#'@importFrom ncdf4 ncvar_get
Expand Down
2 changes: 1 addition & 1 deletion R/get_temp.R
Expand Up @@ -18,7 +18,7 @@
#'Jordan S. Read, Luke A. Winslow
#'@examples
#'sim_folder <- run_example_sim(verbose = FALSE)
#'nc_file <- file.path(sim_folder, 'output.nc')
#'nc_file <- file.path(sim_folder, 'output/output.nc')
#'temp_surf <- get_temp(nc_file, reference = 'surface', z_out = c(0,1,2))
#'temp_bot <- get_temp(nc_file, reference = 'bottom', z_out = c(0,1,2))
#'temp_bot <- get_temp(nc_file)
Expand Down
2 changes: 1 addition & 1 deletion R/get_var.R
Expand Up @@ -21,7 +21,7 @@
#'Jordan S. Read, Luke A. Winslow
#'@examples
#'sim_folder <- run_example_sim(verbose = FALSE)
#'nc_file <- file.path(sim_folder, 'output.nc')
#'nc_file <- file.path(sim_folder, 'output/output.nc')
#'# list variables in the netcdf output from GLM:
#'print(sim_vars(nc_file))
#'evaporation <- get_var(nc_file, var_name = "evap")
Expand Down
2 changes: 1 addition & 1 deletion R/get_wind.R
Expand Up @@ -11,7 +11,7 @@
#'Luke A. Winslow, Jordan S. Read
#'@examples
#'sim_folder <- run_example_sim(verbose = FALSE)
#'nc_file <- file.path(sim_folder, 'output.nc')
#'nc_file <- file.path(sim_folder, 'output/output.nc')
#'wind <- get_wind(file = nc_file)
#'
#'@export
Expand Down
2 changes: 1 addition & 1 deletion R/nc_helpers.R
Expand Up @@ -16,7 +16,7 @@ close_glm_nc <- function(glm_nc){
# Summary: Returns the converted time vector in R format
#' @importFrom ncdf4 ncvar_get
get_time <- function(glm_nc){
hours_since <- ncvar_get(glm_nc, "time")
hours_since <- as.numeric(ncvar_get(glm_nc, "time"))
time_info <- get_time_info(glm_nc)

time <- time_info$startDate + time_info$time_unit * hours_since * 60*60*24
Expand Down
8 changes: 4 additions & 4 deletions R/plot_compare_stage.R
Expand Up @@ -8,13 +8,12 @@
#'@examples
#'
#'sim_folder <- run_example_sim(verbose = FALSE)
#'nc_file <- file.path(sim_folder, 'output.nc')
#'nc_file <- file.path(sim_folder, 'output/output.nc')
#'
#'field_file <- file.path(sim_folder, 'LakeMendota_stage_USGS05428000.csv')
#'
#'plot_compare_stage(nc_file, field_file) ##makes a plot!
#'
#'
#'@export
plot_compare_stage = function(nc_file, field_file, fig_path = NULL, ...){

Expand All @@ -32,8 +31,9 @@ plot_compare_stage = function(nc_file, field_file, fig_path = NULL, ...){
# Bind modeled and observed data
stage_all = bind_rows(stage_mod,stage_obs)

h3 = ggplot(stage_all) + geom_path(aes(x = DateTime, y = Stage, color = Group)) +
geom_point(aes(x = DateTime, y = Stage, color = Group)) +
h3 = ggplot(stage_all) +
geom_path(aes(x = .data$DateTime, y = .data$Stage, color = .data$Group)) +
geom_point(aes(x = .data$DateTime, y = .data$Stage, color = .data$Group)) +
scale_color_manual(values = c('black','lightblue4')) +
theme_bw() + theme(legend.title = element_blank())

Expand Down

0 comments on commit 174c0d4

Please sign in to comment.