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

Question

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

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)))

该函数没有正确理解条件，只为所有eci 值返回ec，而不管它应该只在eci < abs(ec) 时返回ec 的条件

以下是我尝试在 R 中重新创建的示例。

Answer 1

这是一种更加矢量化的方法 - 它一次计算所有 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

Answer 2

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

另外，如前所述，eci 不能是长度为 1000 的向量。它必须从 1E-6 开始，并以 1E-6 递增每个循环（根据该图）。您的退出标准很复杂，因为我们不知道0in 是什么，所以无法正确重构。

最后，在apply 中，您可以调用用户定义的函数，但您必须向它传递一些东西。

话虽如此，这就是你所追求的吗？

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 条件。根据ec 和eci 的值，该函数返回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