Сохраняйте значения в матрице сразу, а не печатайте их внутри цикла for после каждой итерации - PullRequest
0 голосов
/ 28 февраля 2019

* Забыл пакеты в предыдущей версии этого вопроса. У меня есть рабочий код, который включает в себя несколько входных параметров, 2 основные функции, которые используют эти параметры, и, наконец, цикл for, который я использую для создания моего окончательного кадра данных с именем temp1.Проблема в том, что я использую цикл for для печати значений после каждой итерации, что делает мое блестящее приложение очень медленным, особенно если у меня много итераций.Альтернатива, о которой я мог бы подумать, - заменить этот цикл for методом, который сохраняет то, что происходит внутри цикла for, в кадр данных или матрицу, а затем использует его для создания temp1.Я пробовал что-то вроде

# If function is not vectorised
    df %>% 
      mutate(Dose = map2_dbl(Kd1Par, Kd2Par, MyFunction))

, но это не работает должным образом.

library(deSolve)
library(minpack.lm)
library(reshape2)
library(pracma)


###READ THE PARAMETERS
reqMinInh <- 90 # (%) Min inhibition of Target
nd <- 4 # Number of doses
tau <-7
endTime <- (nd+1)*tau
BW <- 70
MW <- 150e3
nIter <- 3
Kd1Par <- logspace(-1.98,1.698,n = nIter)
Kd2Par <- logspace(-1.98,1.698,n = nIter)
myDose <- matrix(c(0), nrow= length(Kd1Par), ncol = length(Kd2Par))
Kon_m1 <- 1.3824 # (1/(nmol/L)/day)
Kon_m2 <- 1.3824 # (1/(nmol/L)/day)
Base1 <- 0.1
Base2 <- 0.1
HL1 <- 100
HL2 <- 100
Kint_m1  <- 0.693*60*24/HL1 # (1/day)
Kint_m2  <- 0.693*60*24/HL2 # (1/day)
Kdeg_m1  <- Kint_m1 # (1/day)
Kdeg_m2  <- Kint_m2 # (1/day)
Ksyn_m1  <- Base1*Kdeg_m1 # (nmol/L/day)
Ksyn_m2  <- Base2*Kdeg_m2 # (nmol/L/day)
Vp  <- 3 # (L) Ref: Vaishali et al. 2015
Vph  <- 3.1 # (L)  Ref: Tiwari et al. 2016
Vt  <- 0.192 # (L) Spleen, Ref: Davis et al. 1993
k_01  <- 1 # (1/day)  Ref: Leonid Gibiansky
CL  <- 0.24 # (L/day)  Ref: Leonid Gibiansky
K_el  <- CL/Vp # (1/day)
k_pph  <- 0.186 # (1/day) Ref: Tiwari et al. 2016
k_php  <- 0.184 # (1/day) Ref: Tiwari et al. 2016
Ktp  <-  0.26 # (1/day)
Kpt  <- 0.004992 # (1/day)
times <- seq(from = 0, to = endTime, by =0.1)
yInit <- c(Ap = 0.0, Dp = 0.0, Dt = 0.0, 
           M1 = Base1, M2 = Base2,
           DtM1 = 0.0, DtM2 = 0.0, DtM1M2 = 0.0, Dph = 0.0) 

##FUNCTIONS

#FUNCTION 1
derivs_pk1 <- function(t, y, parms) {
  with(as.list(c(y,parms)),{
    dAp_dt <- -k_01*Ap
    dDp_dt <- k_01*Ap/Vp -K_el*Dp +Vt/Vp*Ktp*Dt -Kpt*Dp +Vph/Vp*k_php*Dph -k_pph*Dp
    dDt_dt <- Vp/Vt*Kpt*Dp -Ktp*Dt -Kon_m1*Dt*M1 +Koff_m1*DtM1 -Kon_m2*Dt*M2 +Koff_m2*DtM2
    dM1_dt <- Ksyn_m1 -Kdeg_m1*M1 -Kon_m1*Dt*M1 +Koff_m1*DtM1 -Kon_m1*DtM2*M1 +Koff_m1*DtM1M2
    dM2_dt <- Ksyn_m2 -Kdeg_m2*M2 -Kon_m2*Dt*M2 +Koff_m2*DtM2 -Kon_m2*DtM1*M2 +Koff_m2*DtM1M2
    dDtM1_dt <- -Kint_m1*DtM1 -Koff_m1*DtM1 +Kon_m1*Dt*M1 -Kon_m2*DtM1*M2 +Koff_m2*DtM1M2
    dDtM2_dt <- -Kint_m2*DtM2 -Koff_m2*DtM2 +Kon_m2*Dt*M2 -Kon_m1*DtM2*M1 +Koff_m1*DtM1M2
    dDtM1M2_dt <- Kon_m2*DtM1*M2 -Koff_m2*DtM1M2 +Kon_m1*DtM2*M1 -Koff_m1*DtM1M2 -Kint_m1*DtM1M2 -Kint_m2*DtM1M2
    dDph_dt <- Vp/Vph*k_pph*Dp - k_php*Dph
    list(c(dAp_dt,dDp_dt,dDt_dt,dM1_dt,dM2_dt,dDtM1_dt,dDtM2_dt,dDtM1M2_dt,dDph_dt))
  })
}

#FUNCTION2
ssq <- function(parmsToOptm){

  Dose <- parmsToOptm[1]

  injectEvents <- data.frame(var = "Ap",
                             time = seq(0,tau*(nd-1),tau),
                             value = Dose*1e6*BW/MW, # (nmol)
                             method = "add")
  pars_pk1 <- c()

  qss_pk10<-ode(times = times, y = yInit, func =derivs_pk1, parms = pars_pk1,events = list(data = injectEvents))
  qss_pk1<- data.frame(qss_pk10)

  temp <- qss_pk1[qss_pk1$time>tau*(nd-2)&qss_pk1$time<tau*(nd-1),]
  inh1 <- (1-temp$M1/Base1)*100
  inh2 <- (1-temp$M2/Base2)*100
  if(min(inh1,inh2) %in% inh1) {
    currMinInh <- inh1
  } else {currMinInh <-inh2}

  ssqres = currMinInh - reqMinInh
  return(ssqres)

}

###FOR LOOP
for (i in 1:length(Kd1Par)){
  for (j in 1:length(Kd2Par)){

    Kd1 <- Kd1Par[i]
    Kd2 <- Kd2Par[j]
    print(c(Kd1 = Kd1Par[i], Kd2 = Kd2Par[j]))
    Koff_m1 <- Kon_m1*Kd1 # (1/day)
    Koff_m2 = Kon_m2*Kd2 # (1/day)

    # Initial guess
    parmsToOptm <- c(10)

    fitval<-nls.lm(par=parmsToOptm,fn=ssq,control = nls.lm.control(ftol = sqrt(.Machine$double.eps),
                                                                   ptol = sqrt(.Machine$double.eps), gtol = 0, diag = list(), epsfcn = parmsToOptm[1]/100,
                                                                   factor = 100, maxfev = integer(), maxiter = 50, nprint = 0))
    myDose[i, j] <- c(coef(fitval))

    print(c(Dose = myDose[i,j]))
  }
}
###Final DF
temp1  <- melt(myDose)
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...