为具有矢量化输入的函数添加条件

Add a condition for a function that has vectorized input

我有以下带有变量的函数:

b <- 1  ; d_uhpc <- 4 
L_joint <- 8  ; A_bar <- 0.31
A_s <- (A_bar/L_joint)*12
L_unb <- 16 ; f_t <- 1.2
E_s <- 29000 ; E_uhpc <- 8000 

ec <- function(x){
  theta <- x[3]
  eci <- seq(10^-3,1,10^-3)
  while (TRUE) {
    fc <- eci * E_uhpc
    c <- (sqrt(A_s^2 * E_s^2 * eci^2 + fc * A_s * E_s * b * d_uhpc *eci + b^2 * f_t^2 * d_uhpc^2) +
               b * f_t * d_uhpc - A_s * E_s *eci)/(b * fc + 2 * b * f_t)
    ec <- (-2*theta*c)/L_unb
    if (eci > abs(ec)) return("Error") else return(ec)
  }
}

# sample rows
strain.analysis <- read.table(text="
L S     theta
1  60 6 0.3876484
2  70 6 0.3650651
3  80 6 0.3619457
4  90 6 0.3089947
5 100 6 0.3131211
6 110 6 0.3479024", header=TRUE)
strain.analysis1 <- cbind(strain.analysis, vars = t(apply(strain.analysis,1,ec)))

该函数没有正确理解条件,只是 returns ec 用于所有 eci 值,而不管条件它应该只 return ececi < abs(ec)

下面是我试图在 R 中重新创建的示例。

这是一种更加矢量化的方法 - 它会一次性计算所有 ec,然后确定哪个 ec 符合条件。这意味着它速度很快,尽管它比@Paul van Oppen 的解决方案占用更多内存。

eci <- seq(10^-6,1,10^-6)
fc <- eci * E_uhpc
const <- (sqrt(A_s^2 * E_s^2 * eci^2 + fc * A_s * E_s * b * d_uhpc *eci + b^2 * f_t^2 * d_uhpc^2) +
            b * f_t * d_uhpc - A_s * E_s *eci)/(b * fc + 2 * b * f_t)

ec <- function(x){
  theta = x[3]
  ec = (-2*theta*const)/L_unb
  comp = eci > abs(ec)
  if (any(comp)) { ## at least one of our eci > abs(ec)
    wm = which.max(comp) 
    if (comp[wm] == FALSE) ##which.max(c(FALSE, FALSE)) will still return something. We need NA_real_ in this case
      return (NA_real_)
    else
      return(ec[wm])
  }
  else
    return(ec[length(x)])
}

cbind(strain.analysis, V = apply(strain.analysis,1,ec))
#>     L S     theta           V
#> 1  60 6 0.3876484 -0.06845568
#> 2  70 6 0.3650651 -0.06447493
#> 3  80 6 0.3619457 -0.06392507
#> 4  90 6 0.3089947 -0.05459143
#> 5 100 6 0.3131211 -0.05531879
#> 6 110 6 0.3479024 -0.06144967

这就是我的想法:while 循环中的 return 语句不正确,因为它们退出函数 ec 而不仅仅是 while 循环。我认为你在 break 之后。当我 运行 你的代码,eci 作为向量固定时,它只对数据集中的每一行执行一个循环。

此外,如前所述,eci 不能是长度为 1000 的向量。它必须从 1E-6 开始,并在每个循环中递增 1E-6(如图所示)。您的退出标准很复杂,因为我们不知道 0in 是什么,所以无法正确重建。

最后,在 apply 中你可以调用一个 user-defined 函数但是你必须给它传递一些东西。

话虽如此,这就是你想要的吗?

b <- 1
d_uhpc <- 4 
L_joint <- 8
A_bar <- 0.31
A_s <- (A_bar/L_joint)*12
L_unb <- 16 
f_t <- 1.2
E_s <- 29000 
E_uhpc <- 8000 

eci <- 1E-6
ec <- function(x){
  #browser()
  theta <- x[3]
  while (TRUE) {
    fc <- eci * E_uhpc
    c <- (sqrt(A_s^2 * E_s^2 * eci^2 + fc * A_s * E_s * b * d_uhpc *eci + b^2 * f_t^2 * d_uhpc^2) +
            b * f_t * d_uhpc - A_s * E_s *eci)/(b * fc + 2 * b * f_t)
    ec <- (-2*theta*c)/L_unb
    if (eci > abs(ec)) break
    eci <- eci + 1E-6
  }
  
  if (c > d_uhpc | ((max(c(abs(ec), eci), na.rm = TRUE))/(min(c(abs(ec), eci), na.rm = TRUE))) - 1 > 0.05) return(NA_real_) else return(ec)
}


v <- apply(strain.analysis, 1, function(x) ec(x))
strain.analysis1 <- cbind(strain.analysis, v)

代码为strain,analysis的每一行循环多次,直到满足break条件。基于 ececi 的值,函数 returns 要么是 NA_real_.

的值

这是我的输出结果:

> strain.analysis1
    L S     theta           v
1  60 6 0.3876484 -0.06845568
2  70 6 0.3650651 -0.06447493
3  80 6 0.3619457 -0.06392507
4  90 6 0.3089947 -0.05459143
5 100 6 0.3131211 -0.05531879
6 110 6 0.3479024 -0.06144967