I want to calculate 'log10' in nimbleCode

[kei]

Hi all,

 

I am a beginner in NIMBLE. 

I am writing a model with nimbleCode and there is a part of it that calculates log10. However,I know  I can't use the functions provided by R in nimbleCode, so I wrote the code to calculate log10 like this. 

log_10_nimble <- nimbleFunction(
  run = function(mat = double(2)) {
    returnType(double(2))
    ln10 = log10(mat)
    return(ln10)
  })

However, when I tried to define this function, 

[Note] Detected possible use of R functions in nimbleFunction run code.
         For this nimbleFunction to compile, 'log10' must be defined as a nimbleFunction, nimbleFunctionList, or nimbleList. 

this message appeared.Calculating log10 There is no point in re-calculating log10 for a function written for a model that calculates log10. How can I calculate log10 in my model?

Kei

[Perry de Valpine]

Hi Kei,

That's a good point that we don't have a good way to do log10 directly.

I guess you can simply use log10(x) = log(x) / log(10). And that could be coded straight into a model without needing a nimbleFunction call.

HTH

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/496f33bc-be21-440f-9484-e5153df00858n%40googlegroups.com.

[kei]

Dear Perry, all nimble users

In the previous mail,  there are some not English comments, and may cause corruption of text.

So, I email again. The content is the same as in the previous section.

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Thank you for your kind teaching. I was surprised to find such a simple solution!

In fact, I have other errors now.
If you don't mind, I would like your help to solve that error.

The area where the error is occurring is coloured red.

# doing if-then-else for not only constants
ifjudge2 <- nimbleFunction(
  run = function(Nw_copy = double(), indexes = integer(2),ind_k = integer(0),N_par = integer(0)) {
    for(i in 1:N_par){
      if (Nw_copy[i,] != 0){
        for(j in 1:Nw_copy[i,]) { 
          indexes[ind_k] <- i
          ind_k <- ind_k+1
        }
      }
    }
    #returnType(integer(2))
    return(nimList(indexes, ind_k))
  })

####################    Model           ########################################
particlessm <- nimbleCode({   

  # priors
  std1 ~ dunif(0, 3000)  
  std2 ~ dunif(0, 1000)
 
  for (i in 1:Npar){  # generate first particles
    x[i,] ~ dunif(min = -1000, max = 2500)  
    y[i,] ~ dunif(min = -2500, max = 1500)
    r[i,] ~ dunif(min = 1, max = 800)
    a[i,] ~ dunif(min = 0, max = 60)
  }
  for (par_iter in 1:t) {
    for(i in 1:Nsite) {
      x_[,i] <- x[,1]
      y_[,i] <- y[,1]
      r_[,i] <- r[,1]
      a_[,i] <- a[,1]
    }
    for (i in 1:Nsite) {  
      d[,i] <- pmax(sqrt((X[,i]-x_[,i])^2+(Y[,i]-y_[,i])^2)-r_[,i], 1) + 1e-10
    }
   
    for(i in 1:Npar){
      Obs[i,] <- obes[par_iter,1:Nsite]
    }
    for (i in 1:Npar) { 
      diff_obe[i] <- sum(abs(a_[i,] - (Obs[i,] + k*log(d[i,])/log(10))))  # log10(x) = log(x) / log(10)
    }
   
    weight_log[1:Npar] <-  -log(std2) - diff_obe/std2     # calculate log_likelihood of each particle. abs ver.

    weight_max <- max(weight_log)
    weight <- exp(weight_log - weight_max)
    weight_normed <- weight/sum(weight)

    loglikeli[par_iter] <- log(sum(weight)/Npar) + weight_max  # log likelihood of time point "t"

    ## residual resampling ---
    Nw <- Npar*weight_normed
   
    Nw_copy[1:Npar,] <- trunc(Nw)  # truncate smaller than 1 
    ind_k <- 1

    reslist <- ifjudge2(Nw_copy = Nw_copy, indexes = indexes, ind_k = ind_k, N_par = Npar) 

    indexes <- reslist[1]
    ind_k <- reslist[2]
  }
})

#end of the model################################################################

# Const
Nparticle = 1000
Consts <- list(t = nrow(df), # the number of time point "t"
               k = 10,
               Npar = Nparticle,
               Nsite = 5,  
               X = matrix(rep(c(0,941.1171,863.6608,1325.5312,1801.3538), times = Nparticle), nrow = Nparticle, byrow = TRUE),  
               Y = matrix(rep(c(0,583.7835,-1402.1910,-1126.3762,-843.9334), times = Nparticle), nrow = Nparticle, byrow = TRUE))

#Data
Data <- list(obs = df[,2:ncol(df)])
# Initial 
Inits <- with(Consts, list(std1 = 50,
                           std2 = 10,
                           x = matrix(0,nrow = Nparticle, ncol = 1),
                           y = matrix(0,nrow = Nparticle, ncol = 1),
                           r = matrix(0,nrow = Nparticle, ncol = 1),
                           a = matrix(0,nrow = Nparticle, ncol = 1),
                           x_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           y_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           r_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           a_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           d = matrix(0,nrow = Nparticle, ncol = Nsite),
                           diff_obe = rep(0, Nparticle),
                           Obs = matrix(0,nrow = Nparticle, ncol = Nsite),
                           weight_log = rep(0, Nparticle),
                           indexes = matrix(0,nrow = Nparticle, ncol = 1),
                           random_uni = matrix(0,nrow = Nparticle, ncol = 1),
                           Nw_copy = matrix(0,nrow = Nparticle, ncol = 1)
                           ))

particleModel <- nimbleModel(particlessm, constants = Consts,
                          data = Data, check = FALSE, inits = Inits)


These two errors have the same content, dimension error.

In the first red-letter place, this is the error.

Error in genVarInfo3() : Dimension of diff_obe is 1, which does not match its usage in 'weight_log[1:1000] <- -log(std2) - diff_obe/std2'.

If you comment out the above error and continue, you will get the same error in the next red section.

Error in genVarInfo3() : Dimension of Nw_copy is 2, which does not match its usage in 'reslist <- ifjudge2(Nw_copy = Nw_copy, indexes = indexes, ind_k = ind_k, N_par = 1000)'.

I think the error is probably caused by a similar reason, but I have no idea what either cause is and I am very desperate!

There are some more signs of a solution if we can find a way to know the size of the variables inside the model, but I can't find a way to do that...

Any help would be very much appreciated!

Many thanks,

kei

2023年9月20日水曜日 15:51:08 UTC+9 kei:

Dear Perry

Thank you for your kind teaching. I was surprised to find such a simple solution!

In fact, I have other errors now.
If you don't mind, I would like your help to solve that error.

The area where the error is occurring is coloured red.

# doing if-then-else for not only constants
ifjudge2 <- nimbleFunction(
  run = function(Nw_copy = double(), indexes = integer(2),ind_k = integer(0),N_par = Npar) {
    for(i in 1:N_par){
      if (Nw_copy[i,] != 0){
        for(j in 1:Nw_copy[i,]) { # num_copies[i] 個だけ i を選択する。
          indexes[ind_k] <- i
          ind_k <- ind_k+1
        }
      }
    }
    #returnType(integer(2))
    return(nimList(indexes, ind_k))
  })

####################    Model           ########################################
particlessm <- nimbleCode({   
  # priors
  std1 ~ dunif(0, 3000)  # これは事前分布(MCMCで推定するパラメータの分布)
  std2 ~ dunif(0, 1000)
 
  for (i in 1:Npar){  # generate first particles
    x[i,] ~ dunif(min = -1000, max = 2500)  # 明示的に～行目に入れると書かないとエラー
    y[i,] ~ dunif(min = -2500, max = 1500)
    r[i,] ~ dunif(min = 1, max = 800)
    a[i,] ~ dunif(min = 0, max = 60)
  }
  for (par_iter in 1:t) {
    for(i in 1:Nsite) {
      x_[,i] <- x[,1]
      y_[,i] <- y[,1]
      r_[,i] <- r[,1]
      a_[,i] <- a[,1]
    }
    for (i in 1:Nsite) {  # たぶんpmaxがverctor or scalarに対してしかできないから、ベクトルに変換してから行う
      d[,i] <- pmax(sqrt((X[,i]-x_[,i])^2+(Y[,i]-y_[,i])^2)-r_[,i], 1) + 1e-10
    }
   
    for(i in 1:Npar){
      Obs[i,] <- obes[par_iter,1:Nsite]
    }
    for (i in 1:Npar) { # ここと上は統合できそう
      diff_obe[i] <- sum(abs(a_[i,] - (Obs[i,] + k*log(d[i,])/log(10))))  # 各行の合計 log10(x) = log(x) / log(10)
    }
   
    weight_log[1:Npar] <- -log(std2) - diff_obe/std2     # calculate log_likelihood of each particle. abs ver.

    weight_max <- max(weight_log)
    weight <- exp(weight_log - weight_max)
    weight_normed <- weight/sum(weight)

    loglikeli[par_iter] <- log(sum(weight)/Npar) + weight_max  # log likelihood of time point "t"

    ## residual resampling ---
    Nw <- Npar*weight_normed
   
    Nw_copy[1:Npar,] <- trunc(Nw)  # truncate smaller than 1 # この行に入れろと明示的に言わないとエラー
    ind_k <- 1

    reslist <- ifjudge2(Nw_copy = Nw_copy, indexes = indexes, ind_k = ind_k, N_par = Npar) 

    indexes <- reslist[1]
    ind_k <- reslist[2]
  }
})

#end of the model################################################################

# Const
Nparticle = 1000
Consts <- list(t = nrow(df), # データの行数(タイムポイント数)
               k = 10,
               Npar = Nparticle,
               Nsite = 5,  # マイク本数
               X = matrix(rep(c(0,941.1171,863.6608,1325.5312,1801.3538), times = Nparticle), nrow = Nparticle, byrow = TRUE),  # マイクの緯度経度
               Y = matrix(rep(c(0,583.7835,-1402.1910,-1126.3762,-843.9334), times = Nparticle), nrow = Nparticle, byrow = TRUE))

#Data
Data <- list(obs = df[,2:ncol(df)])
# Initial 
Inits <- with(Consts, list(std1 = 50,
                           std2 = 10,
                           x = matrix(0,nrow = Nparticle, ncol = 1),
                           y = matrix(0,nrow = Nparticle, ncol = 1),
                           r = matrix(0,nrow = Nparticle, ncol = 1),
                           a = matrix(0,nrow = Nparticle, ncol = 1),
                           x_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           y_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           r_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           a_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           d = matrix(0,nrow = Nparticle, ncol = Nsite),
                           diff_obe = rep(0, Nparticle),
                           Obs = matrix(0,nrow = Nparticle, ncol = Nsite),
                           weight_log = rep(0, Nparticle),
                           indexes = matrix(0,nrow = Nparticle, ncol = 1),
                           random_uni = matrix(0,nrow = Nparticle, ncol = 1),
                           Nw_copy = matrix(0,nrow = Nparticle, ncol = 1)
                           #Nw_copy = matrix(0,nrow = Npar, ncol = 1)  # 正しい次元に修正
                           ))

particleModel <- nimbleModel(particlessm, constants = Consts,
                          data = Data, check = FALSE, inits = Inits)


These two errors have the same content, dimension error.

In the first red-letter place, this is the error.

Error in genVarInfo3() : Dimension of diff_obe is 1, which does not match its usage in 'weight_log[1:1000] <- -log(std2) - diff_obe/std2'.

If you comment out the above error and continue, you will get the same error in the next red section.

Error in genVarInfo3() : Dimension of Nw_copy is 2, which does not match its usage in 'reslist <- ifjudge2(Nw_copy = Nw_copy, indexes = indexes, ind_k = ind_k, N_par = 1000)'.

I think the error is probably caused by a similar reason, but I have no idea what either cause is and I am very desperate!

There are some more signs of a solution if we can find a way to know the size of the variables inside the model, but I can't find a way to do that...

Any help would be very much appreciated!

Many thanks,

kei

2023年9月20日水曜日 11:51:16 UTC+9 Perry de Valpine:

[heather...@gmail.com]

Nimble wants you to specify dimensions for all vectors. So when you type:

weight_log[1:Npar] <-  -log(std2) - diff_obe/std2

It thinks diff_one is just one value. Instead, try 

weight_log[1:Npar] <-  -log(std2) - diff_obe[1:Npar]/std2

For your other error, I suspect a similar fix is needed. 

Note that when you defined weight you use weight_log, which nimble wants you to write as weight_log[1:Npar]. So make sure weight is also written as weight[1:Npar] and weight_normed becomes weight_normed[1:Npar]. Once this trickles down to NW, you should have NW[1:Npar] and your function will be happy. 

    ## residual resampling ---
    Nw[1:Npar]<- Npar*weight_normed[1:Npar]
   

Heather Gaya 

Sent from my iPhone; please excuse any typos

On Sep 20, 2023, at 5:38 AM, kei <keigoe...@gmail.com> wrote:

Dear Perry, all nimble users

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/9d653694-d170-49da-aba0-27200c929cc5n%40googlegroups.com.

[kei]

Dear  Heather Gaya , all nimble users

Thank you for your advice that is  precise and easy to understand !!

Thanks to your advice, I could solve all dimension errors and understand how indexes must be written!

However, as soon as I think the error has been resolved, another error happens...

Error in if (Nw_copy[i] != 0) { : missing value where TRUE/FALSE needed

in my-define-function, ifjudge2().

I guess this error occur because Nw_copy vector has NA.

The reason why I guess so is that when I set if (any(is.na(Nw_copy))){...}  in my function, ifjudge2() like below, 

NimListDef <- nimbleList(x = integer(2), y = integer(0))
ifjudge2 <- nimbleFunction(
  run = function(Nw_copy = double(2), indexes = integer(2),ind_k = integer(0),N_par = integer(0), reslist = NimListDef()) {
    for(i in 1:N_par){
      if (any(is.na(Nw_copy))){     ### "is.na" is here 
        if (Nw_copy[i] != 0){   
          for(j in 1:Nw_copy[i]) { 

            indexes[ind_k] <- i
            ind_k <- ind_k+1
          }
        }  
      }
    }

    reslist$x <- indexes[1:N_par]
    reslist$y <- ind_k
    returnType(NimListDef())
    return(reslist)
  })

I get the same error; if (Nw_copy[i] != 0) { : missing value where TRUE/FALSE needed.

But, if I set if (any(is.null(Nw_copy))){...} , if (any(is.nan(Nw_copy))){...},and so on, the error don't happen.

I don't know how check the source of the error and why NA happens....

I'll re-show my full model below because somewhere are changed.

 

NimListDef <- nimbleList(x = integer(2), y = integer(0))
ifjudge2 <- nimbleFunction(
  run = function(Nw_copy = double(2), indexes = integer(2),ind_k = integer(0),N_par = integer(0), reslist = NimListDef()) {
    for(i in 1:N_par){
        if (Nw_copy[i] != 0){   

          for(j in 1:Nw_copy[i]) { 

            indexes[ind_k] <- i
            ind_k <- ind_k+1
          }
      }
    }

    reslist$x <- indexes[1:N_par]
    reslist$y <- ind_k
    returnType(NimListDef())
    return(reslist)
  })


####  Model  #######################

particlessm <- nimbleCode({  
  # priors
  std1 ~ dunif(0, 3000)  
  std2 ~ dunif(0, 1000)
 
  for (i in 1:Npar){  # generate first particles
    x[i,] ~ dunif(min = -1000, max = 2500)  
    y[i,] ~ dunif(min = -2500, max = 1500)
    r[i,] ~ dunif(min = 1, max = 800)
    a[i,] ~ dunif(min = 0, max = 60)
  }
  for (par_iter in 1:t) {
    for(i in 1:Nsite) {
      x_[,i] <- x[,1]
      y_[,i] <- y[,1]
      r_[,i] <- r[,1]
      a_[,i] <- a[,1]
    }
    for (i in 1:Nsite) {  
      d[,i] <- pmax(sqrt((X[,i]-x_[,i])^2+(Y[,i]-y_[,i])^2)-r_[,i], 1) + 1e-10
    }
    for(i in 1:Npar){
      Obs[i,] <- obes[par_iter,1:Nsite]
    }
    for (i in 1:Npar) { 
      diff_obe[i] <- sum(abs(a_[i,] - (Obs[i,] + k*log(d[i,])/log(10))))  # log10(x) = log(x) / log(10)
    }

    ##Std2 <- nimRep(std2,Npar)
    weight_log[1:Npar] <- -log(std2) - diff_obe[1:Npar]/std2     # calculate log_likelihood of each particle. abs ver.
    weight_max <- max(weight_log[1:Npar]) 
    weight[1:Npar] <- exp(weight_log[1:Npar] - weight_max)
    weight_normed[1:Npar] <- weight[1:Npar]/sum(weight[1:Npar])
    loglikeli[par_iter] <- log(sum(weight[1:Npar])/Npar) + weight_max  # log likelihood of time point "t"
    ## residual resampling ---

    Nw[1:Npar] <- Npar*weight_normed[1:Npar]

    Nw_copy[1:Npar,1] <- trunc(Nw[1:Npar])  # truncate smaller than 1 
    ind_k <- 1
    reslist <- ifjudge2(Nw_copy = Nw_copy[1:Npar,], indexes = indexes[1:Npar,], ind_k = ind_k, N_par = Npar, reslist = NimListDef) 

    indexes2[1:Npar] <- reslist$x  
    ind_k2 <- reslist$y
  }
})

# Const
Nparticle = 1000
Consts <- list(t = nrow(df), 

               k = 10,
               Npar = Nparticle,
               Nsite = 5,  
               X = matrix(rep(c(0,941.1171,863.6608,1325.5312,1801.3538), times = Nparticle), nrow = Nparticle, byrow = TRUE), 
               Y = matrix(rep(c(0,583.7835,-1402.1910,-1126.3762,-843.9334), times = Nparticle), nrow = Nparticle, byrow = TRUE)
               ) 

# data

Data <- list(obs = df[,2:ncol(df)])

# Initial value

Inits <- with(Consts, list(std1 = 50,
                           std2 = 10,
                           x = matrix(0,nrow = Nparticle, ncol = 1),
                           y = matrix(0,nrow = Nparticle, ncol = 1),
                           r = matrix(0,nrow = Nparticle, ncol = 1),
                           a = matrix(0,nrow = Nparticle, ncol = 1),
                           x_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           y_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           r_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           a_ = matrix(0,nrow = Nparticle, ncol = Nsite),
                           d = matrix(0,nrow = Nparticle, ncol = Nsite),
                           diff_obe = rep(0, Nparticle),
                           Obs = matrix(0,nrow = Nparticle, ncol = Nsite),
                           weight_log = rep(0, Nparticle),
                           indexes = matrix(0,nrow = Nparticle, ncol = 1),

                           #random_uni = matrix(0,nrow = Nparticle, ncol = 1),
                           random_uni = rep(0,Nparticle),
                           Nw_copy = matrix(0,nrow = Nparticle, ncol = 1),
                           indexes2 = rep(0, Nparticle)

                           ))

particleModel <- nimbleModel(particlessm, constants = Consts,
                          data = Data, check = FALSE, inits = Inits)

Thank you for your attention, and any help would be much appreciated!

 P.S 

Error also happens in last two lines of my model,

indexes2[1:Npar] <- reslist$x  

ind_k2 <- reslist$y

Error in model$indexes2[1:1000] <<- model$reslist[1]$x :

replacement has length zero

Maybe I also made a mistake for usage of nimbleList, but I don't know where is the mistake..

If it's an easy problem to figure out, we'd like to know.

kei

2023年9月20日水曜日 20:15:16 UTC+9 heather...@gmail.com:

[kei]

Hi all,

Sorry for asking so many questions...

I was able to resolve the error of "na", but the next error does not seem to be resolvable at all.

I would like to hear any advice you can give me...

I reshow error code and model code where considered to be the cause of error.

NimListDef <- nimbleList(x = integer(2), y = integer(0))
ifjudge2 <- nimbleFunction(  
  run = function(Nw_copy = double(2), indexes = integer(2),ind_k = integer(0),N_par = integer(0), reslist = NimListDef()) {
    for(i in 1:N_par){
        if (Nw_copy[i] != 0){   

          for(j in 1:Nw_copy[i]) { 

            indexes[ind_k] <- i
            ind_k <- ind_k+1
          }
      }
    }

    reslist$x <- indexes[1:N_par]
    reslist$y <- ind_k
    returnType(NimListDef())
    return(reslist)
  })

##### model ####

mymodel <- nimbleCode({ 

...

reslist <- ifjudge2(Nw_copy = Nw_copy[1:Npar,], indexes = indexes[1:Npar,], ind_k = ind_k, N_par = Npar, reslist = NimListDef) 

indexes2[1:Npar] <- reslist$x  
ind_k2 <- reslist$y

})

Error in model$indexes2[1:1000] <<- model$reslist[1]$x :

replacement has length zero

kei

2023年9月21日木曜日 12:59:37 UTC+9 kei:

[Perry de Valpine]

Dear Kei,

It is not allowed to use a nimbleList in model code.

HTH

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/855ff8bd-26ad-457f-9aec-952c50ddbcban%40googlegroups.com.