## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval = FALSE-------------------------------------------------------------
# install.packages("TesiproV")

## ----eval = FALSE-------------------------------------------------------------
# remotes::install_gitlab("ls-massivbau/tesiprov", host = "gitlab.tu-dortmund.de")

## -----------------------------------------------------------------------------
library("TesiproV")

## -----------------------------------------------------------------------------
var_fck <- PROB_BASEVAR(
  Id = 1,
  Name = "f_ck",
  Description = "Concrete compressive strength",
  DistributionType = "lnorm",
  Mean = 47.35,
  Sd = 5.86
)

## -----------------------------------------------------------------------------
var.f_ck <- PROB_BASEVAR(Id = 1, Name = "f_ck", Description = "Char. Betondruckfestigkeit", Package = "TesiproV", DistributionType = "lt", DistributionParameters = c(3.85, 0.09, 3, 10))

## -----------------------------------------------------------------------------
pvar_d <- PARAM_BASEVAR(
  Id = 1,
  Name = "d",
  Description = "Effective depth",
  DistributionType = "norm",
  Sd = 10,
  ParamType = "Mean",
  ParamValues = c(seq(200, 800, 50), seq(1000, 3000, 100)) + 10
)

## -----------------------------------------------------------------------------
var.gamma_c <- PROB_DETVAR(Id = 1, Name = "gamma_c", Description = "PSF Concrete", Value = 1.5)

## -----------------------------------------------------------------------------
var.V_Ed <- PARAM_DETVAR(Id = 1, Name = "V_Ed", Description = "applied shear force", ParamValues = c(1, 2, 3, 4, seq(5, 10, 0.5)))

## -----------------------------------------------------------------------------
var1 <- PROB_BASEVAR(Name = "E", Description = "Effect", DistributionType = "norm", Mean = 10, Sd = 1)
var2 <- PROB_BASEVAR(Name = "R", Description = "Resistance", DistributionType = "norm", Mean = 15, Sd = 1.5)


lsf <- SYS_LSF(vars = list(var1, var2), name = "LSF XY")
lsf$func <- function(E, R) {
  R - E
}
# you can run this check to see if all transformations of the variables worked, all needed data is available and the set of vars fits to the limit state function
lsf$check()

## -----------------------------------------------------------------------------
machine <- PROB_MACHINE(name = "MVFOSM", fCall = "MVFOSM", options = list("isExpression" = FALSE, "h" = 0.1))

## -----------------------------------------------------------------------------
machine <- PROB_MACHINE(name = "FORM Rack.-Fies.", fCall = "FORM", options = list("n_optim" = 20, "loctol" = 0.001, "optim_type" = "rackfies"))

machine <- PROB_MACHINE(name = "FORM Lagrange.", fCall = "FORM", options = list("optim_type" = "auglag"))

## -----------------------------------------------------------------------------
machine <- PROB_MACHINE(name = "SORM", fCall = "SORM")

## -----------------------------------------------------------------------------
machine <- PROB_MACHINE(name = "MC CoV 0.05", fCall = "MC_CRUDE", options = list("n_max" = 1e6, "cov_user" = 0.05))

## ----eval = FALSE-------------------------------------------------------------
# help("MC_IS")

## -----------------------------------------------------------------------------
machine <- PROB_MACHINE(name = "MC IS", fCall = "MC_IS", options = list("cov_user" = 0.05, "n_max" = 300000))

## -----------------------------------------------------------------------------
machine <- PROB_MACHINE(name = "MC Subset Sampling", fCall = "MC_SubSam", options = list("variance" = "uniform"))

## ----results="hide"-----------------------------------------------------------
var1 <- PROB_BASEVAR(Name = "E", Description = "Effect", DistributionType = "norm", Mean = 10, Sd = 1)
var2 <- PROB_BASEVAR(Name = "R", Description = "Resistance", DistributionType = "norm", Mean = 15, Sd = 1.5)


lsf <- SYS_LSF(vars = list(var1, var2), name = "LSF XY")
lsf$func <- function(E, R) {
  R - E
}
# you can run this check to see if all transformations of the variables worked, all needed data is available and the set of vars fits to the limit state function
lsf$check()


machine <- PROB_MACHINE(name = "FORM Rack.-Fieß.", fCall = "FORM", options = list("n_optim" = 20, "loctol" = 0.001, "optim_type" = "rackfies"))

ps <- SYS_PROB(
  sys_input = list(lsf),
  probMachines = list(machine),
  debug.level = 0
)

ps$runMachines()

## -----------------------------------------------------------------------------
ps$beta_single

## -----------------------------------------------------------------------------
# Systemberechnungn (System aus zwei gleichen LSF´s)
ps_sys <- SYS_PROB(
  sys_input = list(lsf, lsf),
  probMachines = list(machine),
  sys_type = "serial"
)

ps_sys$runMachines()
ps_sys$calculateSystemProbability()

## -----------------------------------------------------------------------------
ps_sys$beta_sys

## -----------------------------------------------------------------------------
# Force in kN
pvar_F_mod <- PARAM_BASEVAR(
  Name = "Force", DistributionType = "gumbel", Sd = 2, Package = "evd",
  ParamType = "Mean", ParamValues = c(17, 18, 19, 20) 
)

# yield strength in kN, for Mean the mean of the unshifted distribution must be specified.
fy <- PROB_BASEVAR(
  Name = "fy", DistributionType = "slnorm", Package = "brms", Mean = 265000, Sd = 25000, x0 = 160000
)

# plastic moment of resistance in m³ - considered deterministic in this example
W_pl <- 2.5e-4 
# length of cantilever beam in m - considered deterministic in this example
l    <- 2     

lsf_param <- SYS_LSF(vars = list(pvar_F_mod, fy), name = "Spaethe (1992) - Beispiel 2")
lsf_param$func <- function(Force,fy) {
  return(W_pl*fy-l*Force)
}

machine_form <- PROB_MACHINE(name = "FORM", fCall = "FORM")

ps_param <- SYS_PARAM(
  sys_input = list(lsf_param),
  probMachines = list(machine_form)
)

ps_param$runMachines()
print(ps_param$beta_params)

## -----------------------------------------------------------------------------
para_values <- c(17, 18, 19, 20) 
beta_values <- ps_param$beta_params
plot(para_values, beta_values, type = "l", xlab = "Force F [kN]", ylab = "Reliability index β [-]")

## ----eval=FALSE---------------------------------------------------------------
# MC_IS(..., seed = 123, use_threads = 1)
# MC_IS(..., seed = 123, use_threads = 8)

## ----eval=FALSE---------------------------------------------------------------
# MC_IS(..., backend = "parallel", use_threads = 8)

## ----eval=FALSE---------------------------------------------------------------
# library(future)
# plan(multicore, workers = 8)
# 
# MC_IS(..., backend = "future")

## ----eval=FALSE---------------------------------------------------------------
# library(future)
# plan(multisession, workers = 8)
# 
# MC_IS(..., backend = "future")If no future plan is set, execution will be sequential.

## ----eval=FALSE---------------------------------------------------------------
# 250,000 -- 500,000

## ----eval=FALSE---------------------------------------------------------------
# library(future)
# plan(multicore, workers = 16)
# 
# res <- MC_IS(
#   lsf = my_lsf,
#   lDistr = my_distr,
#   n_batch = 300000,
#   backend = "future",
#   seed = 12345
# )

## ----eval = FALSE-------------------------------------------------------------
# ps <- SYS_PARAM(
#   sys_input = list(lsf),
#   probMachines = list(machine)
# )
# 
# ps$runMachines()
# 
# local_beta <- ps$beta
# local_runtime <- ps$res_single[[1]][[1]]$runtime
# 
# # define path with setwd("...")
# ps$printResults("result_file") # Prints a report file in .txt format with all informations about the calculation
# ps$saveProject(4, "project_file") # stores the project in a file

## ----eval = FALSE-------------------------------------------------------------
# readRDS("example_proj_res.rds")
# load("example_proj.RData")

## ----results="hide"-----------------------------------------------------------
Var1 <- PROB_BASEVAR(Id = 1, Name = "X1", DistributionType = "norm", Mean = 0.25, Sd = 1) # kN/m²
Var2 <- PROB_BASEVAR(Id = 2, Name = "X2", DistributionType = "norm", Mean = 0.25, Sd = 1) # m

lsf <- SYS_LSF(vars = list(Var1, Var2), name = "UQLab 2d_hat")
lsf$func <- function(X1, X2) {
  return(20 - (X1 - X2)^2 - 8 * (X1 + X2 - 4)^3)
}

form <- PROB_MACHINE(name = "FORM Rack.-Fieß.", fCall = "FORM")

ps <- SYS_PROB(
  sys_input = list(lsf),
  probMachines = list(form)
)
ps$runMachines()

## -----------------------------------------------------------------------------
ps$beta_single

## ----results ='hide'----------------------------------------------------------
library("evd")
X1 <- PROB_BASEVAR(Id = 1, Name = "Z", DistributionType = "norm", Mean = 100, Cov = 0.04) # kN/m²
X2 <- PROB_BASEVAR(Id = 2, Name = "Fy", DistributionType = "lnorm", Mean = 40, Cov = 0.1) # m
X3 <- PROB_BASEVAR(Id = 2, Name = "M", DistributionType = "gumbel", Package = "evd", Mean = 2000, Cov = 0.1) # m

lsf <- SYS_LSF(vars = list(X1, X2, X3), name = "Nowak & Collins Exp5.11_1")
lsf$func <- function(M, Fy, Z) {
  return(Z * Fy - M)
}

form <- PROB_MACHINE(name = "FORM Rack.-Fieß.", fCall = "FORM")

ps <- SYS_PROB(
  sys_input = list(lsf),
  probMachines = list(form)
)
ps$runMachines()

## -----------------------------------------------------------------------------
ps$beta_single

## ----results="hide"-----------------------------------------------------------
Var1 <- PROB_BASEVAR(Id = 1, Name = "X1", DistributionType = "norm", Mean = 0.25, Sd = 1) # kN/m²
Var2 <- PROB_BASEVAR(Id = 2, Name = "X2", DistributionType = "norm", Mean = 0.25, Sd = 1) # m

lsf1 <- SYS_LSF(vars = list(Var1, Var2), name = "UQLab 2d_hat")
lsf1$func <- function(X1, X2) {
  return(20 - (X1 - X2)^2 - 8 * (X1 + X2 - 4)^3)
}


X1 <- PROB_BASEVAR(Id = 4, Name = "Z", DistributionType = "norm", Mean = 100, Cov = 0.04) # kN/m²
X2 <- PROB_BASEVAR(Id = 5, Name = "Fy", DistributionType = "lnorm", Mean = 40, Cov = 0.1) # m
X3 <- PROB_BASEVAR(Id = 6, Name = "M", DistributionType = "gumbel", Package = "evd", Mean = 2000, Cov = 0.1) # m

lsf2 <- SYS_LSF(vars = list(X1, X2, X3), name = "Nowak & Collins Exp5.11_1")
lsf2$func <- function(M, Fy, Z) {
  return(Z * Fy - M)
}

form <- PROB_MACHINE(name = "FORM Rack.-Fieß.", fCall = "FORM")

ps <- SYS_PROB(
  sys_input = list(lsf1, lsf2),
  sys_type = "serial",
  probMachines = list(form)
)
ps$runMachines()

## -----------------------------------------------------------------------------
ps$beta_single

ps$calculateSystemProbability()
ps$beta_sys

## ----eval = FALSE-------------------------------------------------------------
# params <- list("cov_user" = 0.05, "n_max" = 50000)
# ps$calculateSystemProbability("MC", params)
# 
# params <- list("cov_user" = 0.05, "n_max" = 50000)
# ps$calculateSystemProbability("MCIS", params)
# 
# ps$calculateSystemProbability("MCSUS", params)

## ----results="hide", warning=FALSE--------------------------------------------
var.d <- PARAM_BASEVAR(Id = 1, Name = "d", Description = "Stat. Nutzhöhe", DistributionType = "norm", Sd = 10, ParamType = "Mean", ParamValues = c(seq(200, 800, 50), seq(1000, 3000, 100)) + 10)
pre.d <- c(seq(200, 800, 50), seq(1000, 3000, 100))

var.f_ck <- PROB_BASEVAR(Id = 2, Name = "f_ck", Description = "Char.Betondruckfestigkeit", Package = "TesiproV", DistributionType = "lt", DistributionParameters = c(3.85, 0.09, 3, 10))
pre.f_ck <- 35

var.f_ywk <- PROB_BASEVAR(Id = 3, Name = "f_ywk", Description = "Zugfestigkeit Bügelbewehrung", DistributionType = "norm", Mean = 560, Cov = 0.02)
pre.f_ywk <- 500

var.f_yk <- PROB_BASEVAR(Id = 4, Name = "f_yk", Description = "Zugfestigkeit Längsbewehrung", DistributionType = "norm", Mean = 560, Sd = 30)
pre.f_yk <- 500

var.rho_l <- PROB_BASEVAR(Id = 5, Name = "rho_l", Description = "Biegebewehrungsgrad", DistributionType = "norm", Mean = 0.01, Cov = 0.02)
pre.rho_l <- var.rho_l$Mean

var.c_D <- PROB_BASEVAR(Id = 6, Name = "c_D", Description = "Stützendurchmesser", DistributionType = "norm", Mean = 500 + 0.003 * 500, Sd = 4 + 0.006 * 500)
pre.c_D <- 500

var.theta_c <- PROB_BASEVAR(Id = 7, Name = "theta_c", Description = "Modellunsicherheit c nach DIBT", DistributionType = "lnorm", Mean = 1.0644, Cov = 0.1679)

var.gamma_c <- PROB_DETVAR(Id = 8, Name = "gamma_c", Description = "TSB Beton", Value = 1.5)
var.gamma_s <- PROB_DETVAR(Id = 9, Name = "gamma_s", Description = "TSB Stahl", Value = 1.15)
var.alpha_cc <- PROB_DETVAR(Id = 10, Name = "alpha_cc", Description = "Langzeitfaktor AlphaCC", Value = 0.85)

V_Ed <- c(
  0.6412447, 0.8760515, 1.1424185, 1.4399554, 1.7554079, 2.0188587, 2.3003489, 2.5997038, 2.9167825, 3.2514671, 3.6036567, 3.9732630,
  4.3602076, 6.0800504, 7.0424326, 8.0726216, 9.1702994, 10.3351829, 11.5670173, 12.8655715, 14.2306347, 15.6620132, 17.1595283, 18.7230145,
  20.3523175, 22.0472935, 23.8078076, 25.6337332, 27.5249511, 29.4813487, 31.5028193, 33.5892621, 35.7405811, 37.9566849
)
var.V_Ed <- PARAM_DETVAR(Id = 11, Name = "V_Ed", Description = "Einwirkende Querkraft", ParamValues = V_Ed)

lsf1 <- SYS_LSF(
  vars = list(
    var.d, var.f_ck, var.f_yk, var.f_ywk, var.rho_l, var.c_D,
    var.theta_c, var.V_Ed, var.gamma_c, var.gamma_s, var.alpha_cc
  ),
  name = "Durchstanzwiderstand ohne Bewehrung"
)
lsf1$func <- function(d, f_ck, f_yk, f_ywk, rho_l, c_D, theta_c, V_Ed, gamma_c, gamma_s, alpha_cc) {
  f_cd <- f_ck * alpha_cc / gamma_c
  f_yd <- f_yk / gamma_s
  f_ywd <- f_ywk / gamma_s

  u_0 <- c_D * pi
  u_1 <- (c_D / 2 + 2 * d) * 2 * pi
  u_0.5 <- (c_D / 2 + 0.5 * d) * 2 * pi

  k <- min(1 + sqrt(200 / d), 2)

  if ((u_0 / d) < 4) {
    C_Rd_c <- 0.18 * (0.1 * u_0 / d + 0.6)
  } else {
    C_Rd_c <- 0.18
  }
  rho_l <- min(rho_l, 0.5 * f_cd / f_yd, 0.02)
  C_Rd_c1 <- C_Rd_c * k * (100 * rho_l * f_ck)^(1 / 3)

  v_Rd_c <- C_Rd_c1
  V_Rd_c <- v_Rd_c * d * u_1 * 1e-6 * theta_c
  return(V_Rd_c - V_Ed)
}

form <- PROB_MACHINE(name = "FORM RF.", fCall = "FORM", options = list("n_optim" = 20, "loctol" = 0.0001, "optim_type" = "rackfies"))

ps <- SYS_PARAM(
  sys_input = list(lsf1),
  probMachines = list(form)
)
ps$runMachines()

## ----results="hide", warning=FALSE--------------------------------------------
d <- pre.d
beta <- unlist(unname(ps$beta_params))
plot(x = d, y = beta, type = "l")