Nothing
#' @import compiler
#' @import JuliaConnectoR
#' @import ggplot2
#' @importFrom stats median optim runif
#' @importFrom utils stack write.csv
#'
NULL
#_______________________________________________________________________________
#' @title Check Julia setup
#'
#' @description Checks installation of Julia and install the needed packages.
#' May take little time to run.Only run the first time you use Jflimo.
#'
#' @return Boolean. True if correct setup, False else.
#'
#' @export
julia_setup <- function(){
requireNamespace("JuliaConnectoR")
if (juliaSetupOk()){
juliaEval('
import Pkg
Pkg.add(url = "https://git.metabarcoding.org/lecasofts/flimo/jflimo.git")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("Distributions")
Pkg.add("ForwardDiff")
Pkg.add("Optim")
Pkg.add("Random")
')
return(TRUE)
}
else {
warning("FAILURE :
Uncomplete Julia setup - See JuliaConnectoR's documentation")
return(FALSE)
}
}
#_______________________________________________________________________________
#' @title Load Julia
#'
#' @description Load needed Julia packages. Run to use Jflimo.
#'
#' @return Boolean. True if load is done correctly
#'
#' @export
julia_load <- function(){
tryCatch({
juliaEval('
using Jflimo
using CSV
using DataFrames
using Distributions
using ForwardDiff
using Optim
using Random
')
return(TRUE)
},
error=function(cond) {
warning("FAILURE :
You should run flimo.julia_setup() once
to complete installation and reload flimo.")
return(FALSE)
}
)
}
#_______________________________________________________________________________
#' @title Objective function minimized by flimo
#'
#' @description Computes the summary statistics between simulations w.r.t.
#' Theta and data. This function is to be minimized by flimoptim.
#'
#' @param Theta 1D array. parameters for the simulations.
#' @param quantiles 2D array containing values drawn in U(0,1).
#' Row number = number of simulations.
#' Column number = number of random variables to draw in one simulation.
#' @param data 1D array containing the observations.
#' @param dsumstats Function computing the distance between simulations and data
#' of form dsumstats(simulations, data) where
#' simulations : 2D array and data : 1D array.
#' ncol(simulations) = length(data) mandatory.
#' @param simulatorQ Function of type simulatorQ(Theta, quantiles)
#' where Theta is the parameter set for the simulations
#' and quantiles are drawn in U(0,1).
#' See README for details.
#'
#' @return Numeric value. Distance between summary statistics of data
#' and simulations w.r.t. Theta.
#'
#' @examples
#'
#'quantiles <- matrix(runif(50), nrow = 10)
#'
#'data <- rep(100, 5)
#'
#'dsumstats <- function(simulations, data){
#' mean_simu <- mean(rowMeans(simulations))
#' mean_data <- mean(data)
#' (mean_simu-mean_data)^2
#'}
#'
#'simulatorQ <- function(Theta, quantiles){
#' qpois(quantiles, lambda = Theta)
#'}
#'
#'flimobjective(100, quantiles, data, dsumstats, simulatorQ)
#'
#' @export
flimobjective <- function(Theta, quantiles, data, dsumstats, simulatorQ){
# simulations <- t(apply(quantiles, 1,
# simulatorQ, Theta = Theta))
first_sim <- simulatorQ(Theta, quantiles[1,])
simulations <- matrix(0, nrow = nrow(quantiles), ncol = length(first_sim))
simulations[1,] <- first_sim
if (nrow(quantiles)>=2){
for (i in 2:nrow(quantiles)){
simulations[i,] <- simulatorQ(Theta, quantiles[i,])
}
}
dsumstats(simulations, data)
}
#_______________________________________________________________________________
#' @title Main function to use flimo inference
#'
#' @description Computes several parameter inferences with R optimizer or
#' Julia optimizer in a full Julia mode.
#' In R mode (default) : L-BFGS-B optimization or other methods available for
#' the base::optim function.
#' In Julia mode : either IPNewton with or without Automatic Differentiation,
#' Nelder-Mead or Brent optimization.
#' Argument ndraw is mandatory.
#' You need either to provide data, dsumstats AND simulatorQ
#' OR obj.
#' @param ndraw Integer. Number of random variables to draw
#' for one simulation of the model.
#' @param data 1D array containing the observations.
#' @param dsumstats Summary statistics to measure distance
#' between simulations and data.
#' In R mode : R function of type dsumstats(simulations, data)
#' where simulations : 2D array and data : 1D array.
#' ncol(simulations) = length(data) mandatory.
#' In Julia mode : a string containing the script of the Julia function
#' dsumstats(simulations, data). The name "dsumstats" is mandatory.
#' @param simulatorQ Simulator of the stochastic process with fixed quantiles
#' (see README).
# Either an R function of type simulatorQ(Theta, quantiles) (in mode "R")
#' or a string (in mode "Julia") containing the script of the Julia function
#' simulatorQ(Theta, quantiles).
#' In Julia mode, the name "simulatorQ" is mandatory.
#' Theta is the parameter set for the simulations and
#' quantiles are drawn in U(0,1).
#' @param obj Objective function to minimize.
#' Default : is directly computed from dsumstats and simulatorQ.
#' Either an R function of type objective(Theta, quantiles) (in mode "R")
#' or a string (in mode "Julia") containing the script of the Julia function
#' julia_obj(Theta, quantiles).
#' Warning : could be tricky if mode = "Julia" to call data.
#' In Julia mode, the name "julia_obj" is mandatory.
#' @param nsim Integer. Number of simulations to run for each step
#' of the optimization algorithm.
#' Computation time grows linearly with this number. Default to 10.
#' @param ninfer Integer. Number of independent inferences to run. Default to 1.
#' @param lower 1D array. Lower bounds for parameters. Same length as upper.
#' With Nelder-Mead in Julia mode: only used for starting point.
#' @param upper 1D array. Upper bounds for parameters. Same length as lower.
#' With Nelder-Mead in Julia mode: only used for starting point.
#' @param Theta0 1D array. Initial values of the parameters.
#' Default : mean(lower, upper).
#' @param randomTheta0 Boolean.
#' If True, Theta0 is randomly drawn between lower and upper bounds.
#' @param mode String. "R" (default) or "Julia". See README.
#' @param AD Boolean.
#' Only in Julia mode, uses Automatic Differentiation with IPNewton method.
#' Default to true.
#' @param method String.
#' In Julia mode, allows to choose the optimization method :
#' "IPNewton", "Brent" or "NelderMead". Default : IPNewton.
#' In R mode, allows to choose any of the optimization methods used by
#' base::optim. Default is L-BFGS-B. Random methods do not work with flimo.
#' Bounded methods are L-BFGS-B and Brent.
#' @param obj_threshold Float. Threshold score. If final value of objective is
#' bigger, relaunch the inference if number_tries is not reached.
#' The purpose is to avoid local minima. Default to Inf (no threshold).
#' @param number_tries Integer. Number of tries (inferences) for the objective
#' value to reach a point lower than obj_threshold. Default to 1.
#' @param maxit Integer. Max number of iterations during optimization.
#' Default to 1000.
#' @param time_limit Float. Time limit in second for each inference.
#' Default to no limit. Not available for R mode and Brent method in Julia mode.
#' @param factr Float.
#' In R-mode : control parameter for L-BFGS-B method in stats::optim.
#' Default to 1e7.
#' @param pgtol Float.
#' In R-mode : control parameter for L-BFGS-B method in stats::optim.
#' Default to 0.
#' @param xtol Float.
#' In Julia mode with IPNewton method : xtol option in Optim.Options.
#' Default to 0.
#' @param ftol Float.
#'In Julia mode with IPNewton method : ftol option in Optim.Options.
#Default to 0.
#' @param gtol Float.
#' In Julia mode with IPNewton method : gtol option in Optim.Options.
#' Default to 1e-8.
#' @param reltol Float.
#' In Julia mode with Brent method : reltol of Optim.optimize.
#' Default is sqrt(.Machine$double.eps), about 1e-8.
#' @param abstol Float.
#' In Julia mode with Brent method : abstol of Optim.optimize.
#' Default is .Machine$double.eps, about 1e-16.
#' @param show_trace Boolean. If true, shows standard trace. Default to false.
#' @param store_trace Boolean.
#' If true, stores standard trace as an array of strings.
#' Default to false. Not available for R mode.
#' @param store_quantiles Boolean.
#' If true, stores every quantiles used for inference, to reproduce the results.
#'Default to false.
#' @param par_names vector of names for parameters.
#' Default is "par1", ..., "parn".
#' @param load_julia Boolean. If true, run julia_load.
#' It can take few seconds. Default to False.
#'
#' @return Object of class flimo_result (list)
#' (converted from Julia object in Julia mode) containing every information
#' about convergence results.
#'
#' @examples
#'data <- rep(100, 5)
#'
#'simulatorQ <- function(Theta, quantiles){
#' qpois(quantiles, lambda = Theta)
#'}
#'dsumstats <- function(simulations, data){
#' mean_simu <- mean(rowMeans(simulations))
#' mean_data <- mean(data)
#' (mean_simu-mean_data)^2
#'}
#'
#' flimoptim(5, data, dsumstats, simulatorQ,
#' lower = 50,
#' upper = 150)
#'
#' @export
flimoptim <- function(ndraw,
data = NULL,
dsumstats = NULL,
simulatorQ = NULL,
obj = NULL,
nsim = 10,
ninfer = 1,
lower = 0,
upper = 1,
Theta0 = (lower+upper)/2,
randomTheta0 = FALSE,
mode = c("R", "Julia"),
AD = TRUE,
method = "",
obj_threshold = Inf,
number_tries = 1,
maxit = 1e3,
time_limit = NaN,
factr = 1e7,
pgtol = 0,
xtol = 0,
ftol = 0,
gtol = 1e-8,
reltol = sqrt(.Machine$double.eps),
abstol = .Machine$double.eps,
show_trace = FALSE,
store_trace = FALSE,
store_quantiles = FALSE,
par_names = NULL,
load_julia = FALSE){
if (length(mode) > 1){
mode <- mode[1] #default is R
}
if (mode == "Julia"){
flimoptim_Julia(ndraw, data, dsumstats, simulatorQ,
julia_obj = obj,
nsim = nsim,
ninfer = ninfer,
lower = lower,
upper = upper,
Theta0 = Theta0,
randomTheta0 = randomTheta0,
AD = AD,
method = method,
obj_threshold = obj_threshold,
number_tries = number_tries,
maxit = maxit,
time_limit = time_limit,
xtol = xtol,
ftol = ftol,
gtol = gtol,
reltol = reltol,
abstol = abstol,
show_trace = show_trace,
store_trace = store_trace,
store_quantiles = store_quantiles,
par_names = par_names,
load_julia = load_julia)
}
else {
if (mode != "R"){
message("Default mode : full R optimization")
}
flimoptim_R(ndraw, data, dsumstats, simulatorQ,
obj = obj,
nsim = nsim,
ninfer = ninfer,
lower = lower,
upper = upper,
Theta0 = Theta0,
randomTheta0 = randomTheta0,
obj_threshold = obj_threshold,
method = method,
number_tries = number_tries,
maxit = maxit,
factr = factr,
pgtol = pgtol,
show_trace = show_trace,
store_quantiles = store_quantiles,
par_names = par_names)
}
}
#_______________________________________________________________________________
#' @title Sample function with fixed quantiles
#'
#' @description Replace the sample function in the context of fixed quantiles.
#' Warning : first argument has less features.
#'
#' @param x a vector of one or more elements from which to choose.
#' @param quantiles 1D array containing values drawn in U(0,1).
#' Length has to be equal to size.
#' @param size a non-negative integer giving the number of items to choose.
#' @param replace should sampling be with replacement ?
#' @param prob a vector of probability weights for obtaining
#' the elements of the vector being sampled.
#'
#' @return a vector of length size with elements drawn from x.
#'
#' @examples
#'
#' set.seed(1)
#' quantiles <- runif(40)
#' sampleQ(1:10, 10, quantiles[1:10])
#' sampleQ(1:10, 10, quantiles[11:20])
#' sampleQ(11:20, 10, quantiles[1:10])
#' sampleQ(1:10, 20, quantiles[21:40], replace = TRUE)
#'
#' @export
sampleQ <- function(x, size, quantiles,
replace = FALSE, prob = NULL){
res <- NULL
if (is.null(prob)) prob <- rep(1, length(x))
prob <- prob/sum(prob)
if (replace){
Sprob <- cumsum(prob)/sum(prob)
fsample <- function(i){x[which(quantiles[i]<Sprob)[1]]}
res <- rep(0, size)
for (i in 1:size){
res[i] <- fsample(i)
}
return(res)
}
else {
for (n in 1:size){
Sprob <- cumsum(prob)/sum(prob)
index <- which(quantiles[n]<Sprob)[1]
i <- x[index]
x <- x[-index]
prob <- prob[-index]
res <- c(res, i)
}
return(res)
}
}
#_______________________________________________________________________________
#' @title Internal flimoptim function in R mode
#'
#' @description Computes several parameter inferences with R optimizer.
#'
#' @param ndraw Integer. Number of random variables to draw
#' for one simulation of the model.
#' @param data 1D array containing the observations.
#' @param dsumstats Summary statistics to measure distance
#' between simulations and data.
#' R function of type dsumstats(simulations, data)
#' where simulations : 2D array and data : 1D array.
#' ncol(simulations) = length(data) mandatory.
#' @param simulatorQ Simulator of the stochastic process with fixed quantiles
#' (see README).
# R function of type simulatorQ(Theta, quantiles)
#' Theta is the parameter set for the simulations and
#' quantiles are drawn in U(0,1).
#' @param obj Objective function to minimize.
#' Default : is directly computed from dsumstats and simulatorQ.
#' R function of type objective(Theta, quantiles)
#' @param nsim Integer. Number of simulations to run for each step
#' of the optimization algorithm.
#' Computation time grows linearly with this number. Default to 10.
#' @param ninfer Integer. Number of independent inferences to run. Default to 1.
#' @param lower 1D array. Lower bounds for parameters. Same length as upper.
#' @param upper 1D array. Upper bounds for parameters. Same length as lower.
#' @param Theta0 1D array. Initial values of the parameters.
#' Default : mean(lower, upper).
#' @param randomTheta0 Boolean.
#' If True, Theta0 is randomly drawn between lower and upper bounds.
#' @param obj_threshold Float. Threshold score. If Final value of objective is
#' bigger, relaunch the inference if number_tries is not reached.
#' The purpose is to avoid local minima. Default to Inf (no threshold).
#' @param method String. Either "L-BFGS-B" (default) or any other method used by
#' the base function optim. Stochastic methods do not work with flimo.
#' If you want to provide bounds, you need to use L-BFGS-B or Brent.
#' @param number_tries Integer. Number of tries (inferences) for the objective
#' value to reach a point lower than obj_threshold. Default to 1.
#' @param maxit Integer. Max number of iterations during optimization.
#' Default to 1000.
#' @param factr Float. Control parameter for L-BFGS-B method in stats::optim.
#'Default to 1e7.
#' @param pgtol Float. Control parameter for L-BFGS-B method in stats::optim.
#'Default to 0.
#' @param show_trace Boolean. If true, shows standard trace. Default to false.
#' @param store_quantiles Boolean.
#'If true, stores every quantiles used for inference, to reproduce the results.
# Default to False.
#' @param par_names vector of names for parameters.
#' Default is "par1", ..., "parn".
#'
#' @return Object of class flimo_result (list) containing every information
#' about convergence results.
flimoptim_R <- function(ndraw,
data = NULL,
dsumstats = NULL,
simulatorQ = NULL,
obj = NULL,
nsim = 10,
ninfer = 1,
lower = 0,
upper = 1,
Theta0 = (lower+upper)/2,
randomTheta0 = FALSE,
obj_threshold = Inf,
method = "L-BFGS-B",
number_tries = 1,
maxit = 1e3,
factr = 1e7,
pgtol = 0,
show_trace = FALSE,
store_quantiles = FALSE,
par_names = NULL){
if (method == ""){method <- "L-BFGS-B"}
minimizer <- matrix(rep(NA, ninfer*length(lower)), nrow = ninfer)
if (is.null(par_names)) colnames(minimizer) <- paste0("par", 1:length(lower))
else colnames(minimizer) <- par_names
minimum <- rep(NA, ninfer)
f_calls <- rep(NA, ninfer)
g_calls <- rep(NA, ninfer)
initial_x <- matrix(rep(NA, ninfer*length(lower)), nrow = ninfer)
converged <- rep(NA, ninfer)
message <- rep(NA, ninfer)
if (store_quantiles){
all_quantiles <- array(rep(NA, ninfer*nsim*ndraw),
dim=c(ninfer, nsim, ndraw))
}
time_run <- rep(NA, ninfer)
for (infer in 1:ninfer){
obj_value <- obj_threshold + 1
if (is.finite(obj_threshold)) tries <- 0
else tries <- number_tries - 1
t <- 0 #time
while ((obj_value >= obj_threshold) && (tries < number_tries)){
quantiles <- matrix(runif(ndraw*nsim), nrow = nsim)
if (is.null(obj)){
intern_obj <- function(Theta){
return(flimobjective(Theta, quantiles, data, dsumstats, simulatorQ))
}
}
else {
intern_obj <- function(Theta){
obj(Theta, quantiles)
}
}
if (randomTheta0){
Theta0 <- runif(length(lower))*(upper-lower)+lower
}
start_time <- Sys.time()
if (!any(method == c("L-BFGS-B", "Brent"))){
opt <- stats::optim(par = Theta0, fn = intern_obj,
method = method,
control = list(trace = show_trace,
maxit = maxit,
factr = factr,
pgtol = pgtol))
}
else {
opt <- stats::optim(par = Theta0, fn = intern_obj,
method = method,
lower = lower, upper = upper,
control = list(trace = show_trace,
maxit = maxit,
factr = factr,
pgtol = pgtol))
}
end_time <- Sys.time()
t <- t + end_time - start_time
obj_value <- opt$value
tries <- tries + 1
}
time_run[infer] <- t
initial_x[infer,] <- Theta0
minimizer[infer,] <- opt$par
minimum[infer] <- opt$value
converged[infer] <- opt$convergence == 0
f_calls[infer] <- opt$counts[1]
g_calls[infer] <- opt$counts[2]
if (!is.null(opt$message)){message[infer] <- opt$message}
if (store_quantiles){
all_quantiles[infer,,] <- quantiles
}
}
optim_result <- NULL
optim_result$mode <- "R"
optim_result$method <- method
optim_result$AD <- FALSE
optim_result$minimizer <- minimizer
optim_result$minimum <- minimum
optim_result$converged <- converged
optim_result$initial_x <- initial_x
optim_result$f_calls <- f_calls
optim_result$g_calls <- g_calls
optim_result$message <- message
if (store_quantiles){
optim_result$quantiles <- all_quantiles
}
optim_result$time_run <- time_run
class(optim_result) <- "flimo_result"
optim_result
}
#_______________________________________________________________________________
#' @title Internal flimoptim function in Julia mode
#'
#' @description Computes several parameter inferences with Julia optimizer and
#' either IPNewton with or without Automatic Differentiation, Nelder-Mead
#' or Brent method.
#'
#' @param ndraw Integer. Number of random variables to draw
#' for one simulation of the model.
#' @param data 1D array containing the observations.
#' @param dsumstats Summary statistics to measure distance
#' between simulations and data.
#' String containing the script of the Julia function
#' dsumstats(simulations, data).
#' The name "dsumstats" is mandatory.
#' @param simulatorQ Simulator of the stochastic process
#' with fixed quantiles (see README).
# String containing the script of the Julia function
#' simulatorQ(Theta, quantiles).
#' The name "simulatorQ" is mandatory.
#' Theta is the parameter set for the simulations and
#' quantiles are drawn in U(0,1).
#' @param julia_obj Objective function to minimize.
#' Default : is directly computed from dsumstats and simulatorQ.
#' String containing the script of the Julia function
#' julia_obj(Theta, quantiles).
#' Warning : can be tricky to call data.
#' The name "julia_obj" is mandatory.
#' @param nsim Integer. Number of simulations to run for each step
#' of the optimization algorithm.
#' Computation time grows linearly with this number. Default to 10.
#' @param ninfer Integer. Number of independent inferences to run. Default to 1.
#' @param lower 1D array. Lower bounds for parameters. Same length as upper.
#' @param upper 1D array. Upper bounds for parameters. Same length as lower.
#' @param Theta0 1D array. Initial values of the parameters.
#' Default : mean(lower, upper).
#' @param randomTheta0 Boolean.
#' If True, Theta0 is randomly drawn between lower and upper bounds.
#' @param AD Boolean.
#' Only in Julia mod, uses Automatic Differentiation with IPNewton method.
#' Default to true.
#' @param method String.
#' Allows to choose the optimization method :
#' "Brent", "IPNewton" or "NelderMead". Default : IPNewton.
#' @param obj_threshold Float. Threshold score. If Final value of objective is
#' bigger, relaunch the inference if number_tries is not reached.
#' The purpose is to avoid local minima. Default to Inf (no threshold).
#' @param number_tries Integer. Number of tries (inferences) for the objective
#' value to reach a point lower than obj_threshold. Default to 1.
#' @param maxit Integer. Max number of iterations during optimization.
#' Default to 1000.
#' @param time_limit Float. Time limit in second for each inference.
#' Default to no limit. Not available for Brent method.
#' @param xtol Float. With IPNewton method : xtol option in Optim.Options.
#' Default to 0.
#' @param ftol Float. With IPNewton method : ftol option in Optim.Options.
#' Default to 0.
#' @param gtol Float. With IPNewton method : gtol option in Optim.Options.
#' Default to 1e-8.
#' @param reltol Float. With Brent method : reltol of Optim.optimize.
#' Default is sqrt(.Machine$double.eps), about 1e-8.
#' @param abstol Float. With Brent method : abstol of Optim.optimize.
#' Default is .Machine$double.eps, about 1e-16.
#' @param show_trace Boolean. If true, shows standard trace. Default to false.
#' @param store_trace Boolean.
#' If true, stores standard trace as an array of strings.
#' Default to false. Not available for R mod.
#' @param store_quantiles Boolean.
#' If true, stores every quantiles used for inference, to reproduce the results.
#' Default to false.
#' @param par_names vector of names for parameters.
#' Default is "par1", ..., "parn".
#' @param load_julia Boolean. If true, run julia_load. It can take few seconds.
#' Default to False.
#'
#' @return Object of class flimo_result (list) converted from Julia object
#' containing every information about convergence results.
#'
flimoptim_Julia <- function(ndraw,
data = NULL,
dsumstats = NULL,
simulatorQ = NULL,
julia_obj = NULL,
nsim = 10,
ninfer = 1,
lower = 0,
upper = 1,
Theta0 = (lower+upper)/2,
randomTheta0 = FALSE,
AD = TRUE,
method = "",
obj_threshold = Inf,
number_tries = 1,
maxit = 1e3,
time_limit = NULL,
xtol = 0,
ftol = 0,
gtol = 1e-8,
reltol = sqrt(.Machine$double.eps),
abstol = .Machine$double.eps,
show_trace = FALSE,
store_trace = FALSE,
store_quantiles = FALSE,
par_names = NULL,
load_julia = FALSE){
if (load_julia){
julia_load()
}
if (!is.null(data)){
file_data <- tempfile(pattern = "data_flimoptim", fileext = ".csv")
write.csv(as.matrix(data), file_data,
row.names = FALSE)
cmd_julia <- ""
juliaEval(paste0('julia_data = CSV.read("',file_data,'", DataFrame)[:,1]\n',
'julia_data = convert(Array{Float64,1}, julia_data)'))
unlink(file_data)
}
else{
cmd_julia <- 'julia_data = nothing'
}
cmd_julia <- paste0(cmd_julia,
paste0('julia_ndraw = ',ndraw,
'\njulia_nsim = ',nsim,
'\njulia_ninfer = ',ninfer,
'\njulia_method = "',method,'"',
'\njulia_obj_threshold = ', obj_threshold,
'\njulia_number_tries = ', number_tries,
'\njulia_maxit = ',maxit,
'\njulia_Theta0 = ',"[",paste(Theta0, collapse = ","),"]",
'\njulia_lower = ',"[",paste(lower, collapse = ","),"]",
'\njulia_upper = ',"[",paste(upper, collapse = ","),"]",
'\njulia_xtol = ',xtol,
'\njulia_ftol = ',ftol,
'\njulia_gtol = ',gtol,
'\njulia_reltol = ',reltol,
'\njulia_abstol = ',abstol),
sep = '\n')
if (is.null(time_limit)){cmd_julia <- paste0(cmd_julia, 'julia_time_limit = NaN',
sep = '\n')}
else {cmd_julia <- paste0(cmd_julia, paste0('julia_time_limit = ',time_limit),
sep = '\n')}
if (AD){cmd_julia <- paste0(cmd_julia, 'julia_AD = true', sep = '\n')}
else {cmd_julia <- paste0(cmd_julia, 'julia_AD = false', sep = '\n')}
if (randomTheta0){cmd_julia <- paste0(cmd_julia, 'julia_randomTheta0 = true',
sep = '\n')}
else {cmd_julia <- paste0(cmd_julia, 'julia_randomTheta0 = false', sep = '\n')}
if (show_trace){cmd_julia <- paste0(cmd_julia, 'julia_show_trace = true',
sep = '\n')}
else {cmd_julia <- paste0(cmd_julia, 'julia_show_trace = false', sep = '\n')}
if (store_trace){cmd_julia <- paste0(cmd_julia, 'julia_store_trace = true',
sep = '\n')}
else {cmd_julia <- paste0(cmd_julia, 'julia_store_trace = false', sep = '\n')}
if (store_quantiles){cmd_julia <- paste0(cmd_julia, 'julia_store_quantiles = true',
sep = '\n')}
else {cmd_julia <- paste0(cmd_julia, 'julia_store_quantiles = false', sep = '\n')}
if (is.null(julia_obj)){
cmd_julia <- paste0(cmd_julia,
paste0('julia_obj = nothing','\n', dsumstats, '\n', simulatorQ),
sep = '\n')
}
else {
cmd_julia <- paste0(cmd_julia,
paste0('julia_obj = nothing','\n', dsumstats, '\n', simulatorQ),
sep = '\n')
}
cmd_julia <- paste0(cmd_julia,
"
julia_xtol = convert(Float64, julia_xtol)
julia_ftol = convert(Float64, julia_ftol)
julia_gtol = convert(Float64, julia_gtol)
julia_reltol = convert(Float64, julia_reltol)
julia_abstol = convert(Float64, julia_abstol)
julia_Theta0 = convert(Array{Float64,1}, julia_Theta0)
julia_lower = convert(Array{Float64,1}, julia_lower)
julia_upper = convert(Array{Float64,1}, julia_upper)
opt = Jflimo.flimoptim(julia_ndraw,
data = julia_data,
dsumstats = dsumstats,
simulatorQ = simulatorQ,
obj = julia_obj,
nsim = julia_nsim,
ninfer = julia_ninfer,
lower = julia_lower,
upper = julia_upper,
Theta0 = julia_Theta0,
randomTheta0 = julia_randomTheta0,
AD = julia_AD,
method = julia_method,
obj_threshold = julia_obj_threshold,
number_tries = julia_number_tries,
maxit = julia_maxit,
xtol = julia_xtol,
ftol = julia_ftol,
gtol = julia_gtol,
reltol = julia_reltol,
abstol = julia_abstol,
time_limit = julia_time_limit,
show_trace = julia_show_trace,
store_trace = julia_store_trace,
store_quantiles = julia_store_quantiles)", sep = '\n')
optim_result <- juliaGet(juliaEval(cmd_julia))
class(optim_result) <- "flimo_result"
optim_result$mode <- "Julia"
optim_result$AD <- AD
if (is.null(par_names))
colnames(optim_result$minimizer) <- paste0("par", 1:length(lower))
else colnames(optim_result$minimizer) <- par_names
names(optim_result)[names(optim_result) == "iteration_converged"] <-
"iteration_limit_reached"
optim_result
}
#_______________________________________________________________________________
#' @title Print flimo results
#'
#' @description Prints most important information about inference results.
#'
#' @param x Object of class flimo_result from any mod/method algorithm
#' of the flimo package.
#' @param ... optional args for generic method
#'
#' @return String containing most important information about argument
#' of class flimo_result.
#'
#' @export
#'
print.flimo_result <- function(x, ...){
value <- 'optimization Result\n'
value <- paste0(value, "Mode : ", x$mod,"\n")
if (length(as.character(x$method)) == 1){
value <- paste0(value, "method : ", as.character(x$method),"\n")
}
else{
value <- paste0(value, "method : ", attr(x$method,"JLTYPE"),"\n")
}
value <- paste0(value, "Number of inferences : ", length(x$minimum),"\n")
value <- paste0(value, "Convergence : ", length(which(x$converged)),"/",
length(x$minimum),"\n")
value <- paste0(value,"Mean of minimizer :\n")
value <- paste0(value, paste(colMeans(x$minimizer), collapse = " "),"\n")
value <- paste0(value, "Best minimizer :\n")
value <- paste0(value,
paste(x$minimizer[which(x$minimum == min(x$minimum))[1],],
collapse = " "))
value <- paste0(value, "\nReached minima : from ", min(x$minimum), " to ",
max(x$minimum))
value <- paste0(value, "\nMedian time by inference ", median(x$time_run))
cat(value)
value
}
#' @title Summary of flimo results
#'
#' @description Most important information about inference results.
#'
#' @param object Object of class flimo_result from any mode/method algorithm
#' of the Jflimo package.
#' @param ... optional args for generic method summary
#'
#' @return List containing most important information about argument
#' of class flimo_result.
#' @export
#'
summary.flimo_result <- function(object, ...){
value <- list()
value$Mode <- object$mode
if (length(as.character(object$method)) == 1){
value$method <- as.character(object$method)
}
else{
value$method <- attr(object$method,"JLTYPE")
}
value$number_inferences <- length(object$minimum)
if (!is.null(attr(object$method,"JLTYPE"))){
value$number_converged <- length(which(object$converged))
}
else {
value$number_converged <- length(which(object$converged))
}
value$minimizer <- summary(object$minimizer)
value$minimum <- summary(object$minimum)
value$median_time_inference <- median(object$time_run)
value
}
#_______________________________________________________________________________
#' @title Check if simulator with fixed quantiles is well implemented
#'
#' @description Run simulations to catch random variations.
#' Warning : does not check it formally.
#' Warning : does not check if quantiles are used several times.
#'
#' @param simulatorQ Function of type simulatorQ(Theta, quantiles)
#' where Theta is the parameter set for the simulations and
#' quantiles are drawn in U(0,1).
#' @param ndraw Integer. Number of random variables to draw
#' for one simulation of the model.
#' @param Theta_lower 1D numeric array. Lower bounds of Theta parameters.
#' @param Theta_upper 1D numeric array. Upper bounds of Theta parameters.
#' @param ntheta Integer. Number of Theta parameters to test.
#' @param nruns Integer. For each Theta, number of simulations to run.
#'
#' @return Boolean. True if no random effect was detected, False else.
#'
#' @examples
#' simulatorQ <- function(Theta, quantiles){
#' qpois(quantiles, lambda = Theta)
#'}
#'check_simulator(simulatorQ, 5,
#'Theta_lower = 50, Theta_upper = 150)
#'
#' @export
#'
check_simulator <- function(simulatorQ, ndraw, Theta_lower = 0, Theta_upper = 1,
ntheta = 5, nruns = 3){
if (length(Theta_lower) != length(Theta_upper)){
warning("Warning : no correct dimension")
stop()
}
if (nruns <= 1){
message("check_simulator is useless with nruns <= 1. nruns set to 2.")
nruns <- 2
}
for (i in 1:ntheta){
Theta <- runif(length(Theta_lower))*(Theta_upper-Theta_lower)+Theta_lower
quantiles_1sim <- runif(ndraw)
quantiles <- matrix(rep(quantiles_1sim, nruns), nrow = nruns, byrow = TRUE)
simulations <- tryCatch({
matrix(simulatorQ(Theta, quantiles), nrow = nruns)
}, error = function(e){
message("This simulator is not conceived to work with quantiles")
return(FALSE)
})
tryCatch({if (!all(apply(simulations, 2, function(x)
length(unique(x)) == 1))){
return(FALSE)
}},error = function(e){
return(FALSE)
})
}
TRUE
}
#_______________________________________________________________________________
#' @title Plot main flimo results
#'
#' @description Shows the plots for most important inference results.
#' Default only shows normalized boxplots for each inferred parameter.
#'
#' @param x Object of class flimo_result.
#' @param y unused generic argument.
#' @param ... optional args for generic method
#' @param bins Integer. Number of bins if hist is True.
#' @param hist Boolean. If True, plots the histogram of each inferred parameter.
#' Default to false.
#' @param par_minimum Boolean.
#' If True, plots each inferred parameter by reached minimum. Default to false.
#' @param pairwise_par Boolean.
#' If True, plots each pairs of inferred parameters.Default to false.
#' @param boxplot Boolean.
#' If True, plots the boxplots of each inferred parameter scaled by their mean.
#' Default to true.
#' @param par_names Vector of names for parameters.
#' Default is "par1", ..., "parn".
#'
#' @return Nothing. Prints the asked ggplot objects.
#'
#' @export
plot.flimo_result <- function(x, y, ...,
hist = FALSE,
bins = 1+as.integer(nrow(x$minimizer)^(1/3)),
par_minimum = FALSE,
pairwise_par = FALSE,
boxplot = TRUE,
par_names = NULL){
ind <- NULL #to remove R CHECK NOTE
values <- NULL #to remove R CHECK NOTE
npar <- ncol(x$minimizer)
for (par in 1:npar){
if (hist){ #Histogram for each inferred parameter
print(ggplot2::ggplot()+
geom_histogram(aes(x$minimizer[,par]), bins = bins)+
labs(x = paste0("parameter ", par))+
ggtitle(paste0("Histogram of inferred parameter ", par)))
}
if (par_minimum){#Plot inferred par = f(minimum)
print(ggplot2::ggplot()+
geom_point(aes(x$minimizer[,par], x$minimum))+
scale_y_log10()+
labs(x = paste0("parameter ", par), y = "reached minimum")+
ggtitle(paste0("inferred parameter ", par,
" by reached minimum")))
}
if (pairwise_par){ #Plot par1 = f(par2)
if (npar > 1 & par < npar){
for (par2 in (par+1):npar){
print(ggplot2::ggplot()+
geom_point(aes(x$minimizer[,par], x$minimizer[,par2]))+
labs(x = paste0("Parameter ", par),
y = paste0("Parameter ", par2))+
ggtitle(paste0("Parameter ", par2, " by parameter ", par)))
}
}
}
}
if (boxplot){ #Normalized boxplots
aux <- as.data.frame(x$minimizer)
if (is.null(par_names)) colnames(aux) <- paste0("par", 1:ncol(aux))
else colnames(aux) <- par_names
aux <- stack(aux)
for (par in unique(aux$ind)){
aux[aux$ind == par,]$values <-
aux[aux$ind == par,]$values/mean(aux[aux$ind == par,]$values)
}
print(ggplot2::ggplot(aux)+
geom_boxplot(aes(x = ind, y = values))+
labs(x = "Parameters", y = "Value/mean(Value)")+
ggtitle(paste0("Normalized boxplots\nfor each inferred parameter")))
}
}
#_______________________________________________________________________________
#' @title Plot the objective to be minimized using flimo
#'
#' @description Plot of the objective function with one parameter moving
#' (objective = f(theta_index)).
#' You need either to provide data, dsumstats AND simulatorQ
#' OR obj.
#'
#' @param ndraw Integer.Number of random variables to draw
#' for one simulation of the model.
#' @param data 1D array containing the observations.
#' @param dsumstats Function computing the distance
#' between simulations and data of form dsumstats(simulations, data)
#' where simulations : 2D array and data : 1D array.
#' ncol(simulations) = length(data) mandatory.
#' @param simulatorQ Function of type simulatorQ(Theta, quantiles)
#' where Theta is the parameter set for the simulations and
#' quantiles are drawn in U(0,1).
#' @param obj objective function of type objective(Theta, quantiles).
#' Default : directly computed with "dsumstats" and "simulatorQ".
#' @param quantiles 2D array containing values drawn in U(0,1).
#' Row number = number of simulations. Default: simulated within the function.
#' Column number = number of random variables to draw in one simulation.
#' @param index Integer. Index of the moving parameter.
#' @param other_param Other parameters of the model. If NULL : assume 1D-model.
#' If numeric : 2D-model, one curve.
#' If 1D-array and dim2 is True (default) :
#' 2D-model, one curve by value in other_param.
#' If 1D-array and dim2 is False or 2D-array : (n>2)D-model,
#' one curve by row in other_param.
#' If your model has n>2 dimensions,
#' you should define other_param as a matrix even if
#' you have only one parameter set to test
#' (with as.matrix(t(vect_param)) where vect_param is a 1D-array).
#' @param nsim Integer. Number of simulations to run for each step
#' of the optimization algorithm.
#' Computation time grows linearly with this number. Default to 10.
#' @param lower Numeric. Lower value of the plot.
#' @param upper Numeric. Upper value of the plot.
#' @param dim2 Boolean. True if model is 2-dimensional.
#' @param visualize_min Boolean. If True, show explicitly the minimum point.
#' @param plot_legend Boolean. If True (default), plots the legend.
#' @param npoints Integer. Number of points evaluated. Default = 300.
#' @param add_to_plot ggplot object. If not NULL,
#' will add all curves/points on previous plot instead of creating a new one.
#' Does not change title/labels/limits defined in previous plot.
#'
#' @return ggplot object representing the objective function to be minimized.
#'
#' @examples
#'data <- rep(100, 5)
#'
#'dsumstats <- function(simulations, data){
#' mean_simu <- mean(rowMeans(simulations))
#' mean_data <- mean(data)
#' (mean_simu-mean_data)^2
#'}
#'
#'simulatorQ <- function(Theta, quantiles){
#' qpois(quantiles, lambda = Theta)
#'}
#'
#' plot_objective(5, data, dsumstats, simulatorQ,
#' lower = 0, upper = 200)
#'
#' @export
plot_objective <- function(ndraw,
data = NULL,
dsumstats = NULL,
simulatorQ = NULL,
obj = NULL,
quantiles = NULL,
index = NULL,
other_param = NULL,
nsim = 10,
lower = 0,
upper = 1,
dim2 = TRUE,
visualize_min = TRUE,
plot_legend = TRUE,
npoints = 300,
add_to_plot = NULL){
if (is.null(quantiles)){q <- matrix(runif(ndraw*nsim), nrow = nsim)}
else {q <- quantiles}
if (is.null(obj)){intern_obj <- function(Theta) flimobjective(Theta, q, data,
dsumstats, simulatorQ)}
else {intern_obj <- function(Theta) obj(Theta, q)}
if (is.null(other_param)){
return(plot_objective1D(intern_obj, lower = lower, upper = upper,
visualize_min = visualize_min,
npoints = npoints,
add_to_plot = add_to_plot))
}
else if (is.null(dim(other_param))){
if (dim2){
other_param <- as.matrix(other_param) #convert to length*1 matrix
return(plot_objectivenD(intern_obj, index, other_param,
lower = lower, upper = upper,
visualize_min = visualize_min,
plot_legend = plot_legend,
npoints = npoints,
add_to_plot = add_to_plot))
}
else{
other_param <- as.matrix(t(other_param)) #convert to 1*length matrix
return(plot_objectivenD(intern_obj, index, other_param,
lower = lower,
upper = upper,
visualize_min = visualize_min,
plot_legend = plot_legend,
npoints = npoints,
add_to_plot = add_to_plot))
}
}
else {
return(plot_objectivenD(intern_obj, index, other_param,
lower = lower,
upper = upper,
visualize_min = visualize_min,
plot_legend = plot_legend,
npoints = npoints,
add_to_plot = add_to_plot))
}
}
plot_objective1D <- function(obj,
lower = 0,
upper = 1,
visualize_min = TRUE,
npoints = 300,
add_to_plot = NULL){
#plot_objective when model has 1 parameter
#internal function
x <- seq(lower, upper, length.out = npoints)
y <- sapply(x, FUN = obj)
if (!is.null(add_to_plot)){
p <- add_to_plot
}
else{
p<-ggplot2::ggplot()
}
p <- p+geom_line(aes(x,y))
if (is.null(add_to_plot)){
p <- p + xlim(lower, upper)+
ggtitle(paste("Objective function to minimize by parameter 1"))+
labs(x = paste0("theta_1"), y = "objective value")
}
if (visualize_min){
#approaching minimum with y values
ymin <- min(y, na.rm = TRUE)
xmin <- x[which.min(y)]
p <- p + geom_point(aes(xmin, ymin)) +
geom_segment(aes(x = xmin, xend = xmin, y = 0, yend = ymin), size = 0.25)
}
p
}
plot_objectivenD <- function(obj,
index,
other_param,
lower = 0,
upper = 1,
visualize_min = TRUE,
plot_legend = TRUE,
npoints = 300,
add_to_plot = NULL){
#plot_objective when model has n>=2 parameters
#internal function
#other_param has to be a 2D-array (conversion done in plot_objective function)
other_par <- NULL #remove R CHECK NOTE
xmin <- #remove R CHECK NOTE
ymin <- #remove R CHECK NOTE
if (index < 1 | index > nrow(other_param)+1){
warning("FAILURE : Index not valid")
stop()
}
x <- seq(lower, upper, length.out = npoints)
obj_values <- NULL
if (visualize_min){
minima <- NULL
}
aux <- function(x, other_par){
Theta <- rep(NA, (length(other_par)+1))
Theta[index] <- x
if (index > 1){
Theta[1:(index-1)] <- other_par[1:(index-1)]
}
if (index < length(other_par)+1){
Theta[(index+1):(length(other_par)+1)] <-
other_par[index:length(other_par)]
}
obj(Theta)
}
for (rpar in 1:nrow(other_param)){
par <- other_param[rpar,]
par_ref <- rep(NA, length(par)+1)
par_ref[index] <- "x"
if (index > 1){
par_ref[1:(index-1)] <- as.character(par[1:(index-1)])
}
if (index < length(par)+1){
par_ref[(index+1):(length(par)+1)] <- as.character(par[index:length(par)])
}
par_ref <- paste0(par_ref, collapse = ";")
obj_values <- rbind(obj_values,
data.frame(x = x,
y = sapply(x, FUN = aux, other_par = par),
other_par = as.factor(par_ref)))
if (visualize_min){
x <- obj_values[obj_values$other_par == par_ref,]$x
y <- obj_values[obj_values$other_par == par_ref,]$y
minima <- rbind(minima, data.frame(xmin = x[which.min(y)],
ymin = min(y, na.rm = TRUE),
other_par = as.factor(par_ref)))
}
}
if (!is.null(add_to_plot)){
p <- add_to_plot
}
else{
p <- ggplot2::ggplot(obj_values)
}
p <- p + geom_line(data = obj_values, aes(x,y, color = other_par))+
xlim(lower, upper)+
ggtitle(paste("Objective function to minimize by parameter", index))+
labs(x = paste0("theta_", index), y = "objective value",
color = "Param.")
if (is.null(add_to_plot)){
p <- p +
ggtitle(paste("Objective function to minimize by parameter",index))+
labs(x = paste0("theta_", index), y = "objective value")
}
if (visualize_min){
p <- p + geom_point(data = minima, aes(xmin, ymin, color = other_par))+
geom_segment(data = minima, aes(x = xmin, xend = xmin,
y = 0, yend = ymin, color = other_par),
size = 0.25)
}
if (!plot_legend){
p <- p + theme(legend.position = "none")
}
p
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.