ERROR : R Session Aborted. R encountered a fatal error. The session was terminated

[kei]

Hi there!

I have a problem when I do runMCMC with particlefilter.

I use this model and code, 

## my function

calc_volume = nimbleFunction(  
  run = function(distance = double(0), radius = double(0), mat = double(2),A = double(1), AA = double(1),AAA = double(2)) {
     
    rad_big = min(which(A - radius > 0))
    rad_sma = rad_big-1
    small_rad_line = AAA[rad_sma,]
    big_rad_line = AAA[rad_big,]
    dist_big1 = min(which(small_rad_line - distance > 0))
    dist_big2 = min(which(big_rad_line - distance > 0))
    ### distanceが0の時にエラーが出るかわからないからエラー処理は保留
    dist_small1 = dist_big1-1
    dist_small2 = dist_big2-1
    a1 = AAA[rad_sma,dist_small1] # x of left low
    a2 = AAA[rad_sma,dist_big1] # x of right low
    A1 = AAA[rad_big,dist_small2] # x of left high
    A2 = AAA[rad_big,dist_big2] # x of right high
    aa1 = mat[rad_sma, dist_small1] # y of left low
    aa2 = mat[rad_sma, dist_big1]  # y of right low
    AA1 = mat[rad_big, dist_small2] # y of left high
    AA2 = mat[rad_big, dist_big2]  # y of right high
    rate1_low = (distance - a1)/(a2 - a1)
    rate1_high = (distance - A1)/(A2 - A1)
    rate2 = (radius - A[rad_sma])/(A[rad_big] - A[rad_sma])
    x1 = (a2 - a1)*rate1_low + a1
    x2 = (A2 - A1)*rate1_high + A1
    y1 = (aa2 - aa1)*rate1_low + aa1
    y2 = (AA2 - AA1)*rate1_high + AA1
    x_mid = (x2 - x1)*rate2 + x1
    y_mid = (y2 - y1)*rate2 + y1
  returnType(double(0))
  return(y_mid)
})

### model
## latent variableをx一つにまとめる
Particle_code <- nimbleCode({
  ## produce init particles
  x[1,1] ~ dunif(-1500, 3000)
  x[1,2] ~ dunif(-3000, 2000) # y
  x[1,3] ~ dunif(0, 3500)     # r
 
  for(i in 2:t){
    ## forcast
    x[i,1] ~ T(dnorm(x[i-1, 1], sd = sqrt(std_xy)) ,x_min, x_max) 
    x[i,2] ~ T(dnorm(x[i-1, 2], sd = sqrt(std_xy)) ,y_min, y_max) # invGは分散に対して。
    x[i,3] ~ T(dnorm(x[i-1, 3], sd = sqrt(std_r)) ,0,100000)  # 下限が0.1の時は収束した(普通のMCMCで)
    d[i,1:5] <- sqrt((XX[1:5]-x[i,1])^2+(YY[1:5]-x[i,2])^2)     

    for (j in 1:n_mic) {  # ここでerror  'list' object cannot be coerced to type 'double'
      res_box[i,j] <- calc_volume(distance = d[i,j], radius = x[i,3], mat = mat[1:297,1:171], A = A[1:297], AA = AA[1:297], AAA = AAA[1:297,1:171])
    }
    vol[i, 1:5] <- 120 + 20*log(rho*res_box[i, 1:5])/log(10)  # log10(x) = log(x)/log(10)
    for (k in 1:n_mic) {
      Data[i,k] ~ ddexp(location = vol[i,k], scale = std2)
    }
  }
 
  # priors
  std_xy ~ dunif(0, 1000)
  std_r ~ dunif(0, 1000)
  std2 ~ dunif(0,100)
  rho ~ dunif(0, 100)
 
})


mat = as.matrix(read.csv("sound_volume_matrix_Log10_2.csv", header = FALSE))
A = unlist(read.csv("calculated_radius_matrix_Log10_2.csv", header=FALSE))
AA = unlist(read.csv("number_of_distance_Log10_2.csv", header=FALSE))
AAA = as.matrix(read.csv("radius_distance_matrix_Log10_2.csv", header=FALSE))

x1 = 0   # マイクロホン座標
y1 = 0
x2 = 941.1171
y2 = 583.7835
x3 = 863.6608
y3 = -1402.1910
x4 = 1325.5312
y4 = -1126.3762
x5 = 1801.3538
y5 = -843.9334
XX = c(x1, x2, x3, x4, x5)
YY = c(y1, y2, y3, y4, y5)
x_min = min(XX) - 3000
x_max = max(XX) + 3000
y_min = min(YY) - 3000
y_max = max(YY) + 3000
res_box = array(rep(0,5))
vol = array(rep(0,5))

df = read.csv("df8_4_q25Log.csv", header = FALSE)[,-1] 
df[1,] = 1
df2 <- data.frame( lapply(df, as.numeric) )  # dataをnumericに
time_point = nrow(df)  #

n_mic = 5
n_par = 500
Inits = list(std_xy = 1, std2 = 1, std_r = 1, rho = 1,
             x = matrix(rep(1, 3*time_point), ncol = 3))
Const = list(x_min = x_min, x_max = x_max, y_min = y_min, y_max = y_max, XX = XX, YY = YY, mat = mat, A = A, AA = AA, AAA = AAA, n_mic = n_mic, t = time_point)

Particle_model <- nimbleModel(Particle_code, constants = Const,
                      data <- list(Data = df2),inits = Inits,
                      check = TRUE)

particleconf <- configureMCMC(Particle_model,
                              node = NULL,  print = TRUE, enableWAIC = TRUE) #
particleconf$addMonitors(c("std_xy", "std2", "std_r", "rho"))
particleconf$addSampler(target = c("std_xy", "std2", "std_r", "rho"), type = "RW_PF_block",
                      control = list(adaptive = TRUE, propCov= diag(4), adaptScaleOnly = FALSE,adaptInterval = 50, 
                                     pfNparticles = n_par, pfType = "bootstrap",
                                     latents = c("x"),
                                     pfOptimizeNparticles = FALSE)) 

## compile my function
cCalc_volume = compileNimble(calc_volume)  # 別々でコンパイル
cParticle_model <- compileNimble(Particle_model)
Particleconf <- buildMCMC(particleconf)  # 
cParticleMCMC <- compileNimble(Particleconf, project = cParticle_model, resetFunctions = TRUE)   
res = runMCMC(cParticleMCMC, niter = 2000, nchains = 1, nburnin = 500, summary = TRUE, WAIC = TRUE)    ←←←← Thiscause R session error!!

When I try runMCMC(), there are no error code, just error pop up happens!

"R Session Aborted. R encountered a fatal error. The session was terminated"  will be shown and new R session starts. 

I have no idea why this problem happens because I can calculate nimbleMCMC and bootstrap filter with the same model I used in runMCMC.

But only PMCMC have this trouble.

How can I solve this problem?

By the way I can calculate PMCMC(doing runMCMC) with sample code in nimble user manual, so I think my computer has a potential to do PMCMC.

please give me any suggestions and comments! 

Thank you for reading!

Best, 

Kei 

[kei]

Hi all,

To solve this problem, I checked my defined function(calc_volume).
Then, aborted error didn't happen when I commented out after small_rad_line = AAA[rad_sma,], this is the third line of my function.

My function can work until the second line.

Thus, small_rad_line = AAA[rad_sma,] is the error source.

How can I solve this terminating error?

I have no idea to debug this error because there is no error message, and this section can work as normal R code, not nimble code.

BTW, this AAA is the matrix that have 297 rows and 171 cols.

Thanks,

kei

2024年5月18日土曜日 21:06:13 UTC+9 kei:

[Perry de Valpine]

Kei,

My guess is that rad_sma is getting an invalid value. R execution is much more forgiving than C++ execution for that. Could rad_big ever be set to 1, making rad_sma be 0, which would invalid? I'm also guessing that the problem could arise in PMCMC because it proposes parameter values and then given those values it runs the particle filter. Compared to running a bootstrap filter on its own or running an MCMC on all model parameters and states, the PMCMC could have more possibility of running the particle filter with very bad parameter values and hence getting into some corner case that causes a problem.

Here are some debugging strategy ideas:

For running uncompiled, you could add after line 2 of calc_volume:

if((rad_sma < 1) || (rad_sma > 297)) browser()

That will let you investigate what input values are giving invalid  rad_sma, if uncompiled execution gets to such a situation fast enough to find it.

If you need to be running in compiled execution to get to the problem, you could add at the same point:

print(distance, " ", radius)

This will show you all the distance and radius values with which calc_volume is being called. With luck it will show you the values involved in the crash just before the crash occurs.

If you get familiar with nimbleRcall, you can even make a nimbleRcall to an R function, passing the values of e.g. distance, radius, and rad_sma and then in R doing only:

if((rad_sma < 1) || (rad_sma > 297)) browser()

which would again allow you to stop and look at the distance and radius values associated with any invalid rad_sma values.

Another strategy (which I'd suggest in combination with the above) is to use set.seed(<any arbitrary seed>) before you run, and see if the iteration on which the crash occurs is repeatable. This can help you get close to the problem.

Since running MCMCs over and over again can be a pain, another strategy would be to compile the model object and then manually change values of the parameters and/or x (e.g. cmodel$x <- your_values) and call 

cmodel$calculate("res_box") or cmodel$calculate("res_box[3, 5]") for whatever specific indices you want to check,

and look at the result (cmodel$res_box). In this way you can try to figure out any parameter and/or x settings that cause a problem.

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/6854221f-5777-4563-9a69-e7990f4247ffn%40googlegroups.com.

[kei]

Dear Perry and All nimble users

Thank you for your suggestion!

It is really useful for debugging.

I tried debugging with your "print() method".

Actually I can't understand why and where this error happens.

I set print() function anywhere in my code, then I found some curious situations.

1. Error occurs different places even the same seed value.

I attach two screenshots to this mail, they are the screenshots of when the error occurred.

As you can see, error occurred at different points.

2. Error happens between two print()s even these are next to each other!!!

In my code below, error occurs 

print("vol ",vol, "  i ",i)
print("std2_1 ", std, "  i ",i)

between these 2 prints. I guess so because vol was printed but std was not printed just before Aborting error.

3. (I guess this can't be the source of error) "for function" acts in curious order.

In my code, 

for() calculates calc_volume() from i=1 to i=t at first, then it calculates calc_2() from i=1 to i=t (I guess so from the order of print() in screenshots).

I'm actually despaired now.

It looks impossible to debug, I thought.

Thank you for reading!

It would be help if any comments or suggestions!

Thanks,

kei

This is my code, 

calc_volume = nimbleFunction(  
  run = function(distance = double(0), radius = double(0), mat = double(2),A = double(1), AA = double(1),AAA = double(2),i = integer(0),j = integer(0)) {

    rad_big = min(which(A - radius > 0))
    rad_sma = rad_big-1  

    #if((rad_sma < 1) || (rad_sma > 297)) browser()
    print(distance, " ", radius, "  i ",i,"  j ",j)

    small_rad_line = AAA[rad_sma,]
    big_rad_line = AAA[rad_big,]
   

    dist_big1 = min(which(small_rad_line - distance > 0))

    dist_big2 = min(which(big_rad_line - distance > 0))
   
    ### distanceが0の時にエラーが出るかわからないからエラー処理は保留
    dist_small1 = dist_big1-1
    dist_small2 = dist_big2-1
   
    a1 = AAA[rad_sma,dist_small1] # x of left low
    a2 = AAA[rad_sma,dist_big1] # x of right low
    A1 = AAA[rad_big,dist_small2] # x of left high
    A2 = AAA[rad_big,dist_big2] # x of right high
    aa1 = mat[rad_sma, dist_small1] # y of left low
    aa2 = mat[rad_sma, dist_big1]  # y of right low
    AA1 = mat[rad_big, dist_small2] # y of left high
    AA2 = mat[rad_big, dist_big2]  # y of right high
   
    rate1_low = (distance - a1)/(a2 - a1)
    rate1_high = (distance - A1)/(A2 - A1)
    rate2 = (radius - A[rad_sma])/(A[rad_big] - A[rad_sma])
    x1 = (a2 - a1)*rate1_low + a1
    x2 = (A2 - A1)*rate1_high + A1
    y1 = (aa2 - aa1)*rate1_low + aa1
    y2 = (AA2 - AA1)*rate1_high + AA1
    x_mid = (x2 - x1)*rate2 + x1
    y_mid = (y2 - y1)*rate2 + y1

    print("y_mid ",y_mid, "  i ",i,"  j ",j)
    returnType(double(0))
    return(y_mid)
  })

calc_2 = nimbleFunction(  # AAAとかmatもdoubleにしないとエラー
  run = function(res_box = double(2), rho = double(0),i = integer(0), std = double(0)) {
    print("res_box   ",res_box[i, 1:5], "  i ",i)
   
    print("rho  ",rho, "  i ",i)
    vol <- 120 + 20*log(rho*res_box[i, 1:5])/log(10)  # log10(x) = log(x)/log(10)
   
    print("vol ",vol, "  i ",i)
    print("std2_1 ", std, "  i ",i)        ←←← Error occurred here!!
   
    returnType(double(1))
    return(vol)
  })

### My model
Particle_code_std2s <- nimbleCode({

  ## produce init particles
  x[1,1] ~ dunif(-1500, 3000)
  x[1,2] ~ dunif(-3000, 2000) # y
  x[1,3] ~ dunif(0, 3500)     # r
 
  for(i in 2:t){
    ## forcast    

    x[i,1] ~ T(dnorm(x[i-1, 1], sd = sqrt(std_xy)) ,x_min, x_max) # nimbleでは明示的にsd= or tau=にしないといけない

    x[i,2] ~ T(dnorm(x[i-1, 2], sd = sqrt(std_xy)) ,y_min, y_max) # invGは分散に対して。
    x[i,3] ~ T(dnorm(x[i-1, 3], sd = sqrt(std_r)) ,0,100000)  # 下限が0.1の時は収束した(普通のMCMCで)
   

    d[i,1:5] <- sqrt((XX[1:5]-x[i,1])^2+(YY[1:5]-x[i,2])^2)  # ちゃんとベクトルを扱うときは全部を指定する

   
    for (j in 1:n_mic) {  # ここでerror  'list' object cannot be coerced to type 'double'

      res_box[i,j] <- calc_volume(distance = d[i,j], radius = x[i,3], mat = mat[1:297,1:171], A = A[1:297], AA = AA[1:297], AAA = AAA[1:297,1:171], i = i, j = j)
   
    }

    vol[i, 1:5] <- calc_2(res_box = res_box[1:t,1:5], rho = rho, i = i, std = std2_1)

    Data[i,1] ~ ddexp(location = vol[i,1], scale = std2_1)
    Data[i,2] ~ ddexp(location = vol[i,2], scale = std2_1)
    Data[i,3] ~ ddexp(location = vol[i,3], scale = std2_1)
    Data[i,4] ~ ddexp(location = vol[i,4], scale = std2_1)
    Data[i,5] ~ ddexp(location = vol[i,5], scale = std2_1)
  }
 
  # priors
  std_xy ~ dunif(0, 1000000)
  std_r ~ dunif(0, 1000)
  std2_1 ~ dunif(0,100)

  rho ~ dunif(0, 100)
 
})
mat = as.matrix(read.csv("sound_volume_matrix_Log10_2.csv", header = FALSE))
A = unlist(read.csv("calculated_radius_matrix_Log10_2.csv", header=FALSE))
AA = unlist(read.csv("number_of_distance_Log10_2.csv", header=FALSE))
AAA = as.matrix(read.csv("radius_distance_matrix_Log10_2.csv", header=FALSE))

AAA[is.na(AAA)] = 10000000#NULL
#AAA

x1 = 0   # マイクロホン座標
y1 = 0
x2 = 941.1171
y2 = 583.7835
x3 = 863.6608
y3 = -1402.1910
x4 = 1325.5312
y4 = -1126.3762
x5 = 1801.3538
y5 = -843.9334
XX = c(x1, x2, x3, x4, x5)
YY = c(y1, y2, y3, y4, y5)
x_min = min(XX) - 3000
x_max = max(XX) + 3000
y_min = min(YY) - 3000
y_max = max(YY) + 3000
res_box = array(rep(0,5))
vol = array(rep(0,5))

df = read.csv("df8_4_q25Log.csv", header = FALSE)[,-1]  # 1列目はいらない列(何時のデータかが書いてある列)
df[1,] = 1# あえて1行目は使わない行が書いてある(尤度の比較とかが2行目からだから)

df2 <- data.frame( lapply(df, as.numeric) )  # dataをnumericに

Data = list(Data = df2)

time_point = nrow(df)  #

n_mic = 5

Inits = list(std_xy = 1, std2_1 = 1, std2_2 = 1, std2_3 = 1, std2_4 = 1, std2_5 = 1,

             std_r = 1, rho = 1, x = matrix(rep(1, 3*time_point), ncol = 3))
Const = list(x_min = x_min, x_max = x_max, y_min = y_min, y_max = y_max, XX = XX, YY = YY, mat = mat, A = A, AA = AA, AAA = AAA, n_mic = n_mic, t = time_point)

Particle_model <- nimbleModel(Particle_code_std2s, constants = Const,
                              data <- Data,inits = Inits,
                              check = TRUE)

cParticle_model <- compileNimble(Particle_model)

## watch the action
#res_box[i,j] <- calc_volume(distance = d[i,j], radius = x[i,3], mat = mat[1:297,1:171], A = A[1:297], AA = AA[1:297], AAA = AAA[1:297,1:171])
i = 19
#cParticle_model$vol[,i] = c()
cParticle_model$x[2,i] = 721.435
cParticle_model$rho
#cParticle_model$rad_big
cParticle_model$calculate("res_box[2, 3]")
cParticle_model$res_box[2, 3] = 35.1879
cParticle_model$calculate("vol[2, 3]")
cParticle_model$vol[2,]


n_par = 1  # the number of particles

particleconf <- configureMCMC(Particle_model,
                              node = NULL,  # 
                              print = TRUE, enableWAIC = TRUE) # enableWAIC = TRUEにしないと、runMCMCでWAIC=TRUEにできない
particleconf$addMonitors(c("std_xy", "std2_1","std_r", "rho")) #
## 相関が大きいからblock_samplerした方が収束が良くなる！alphaとbetaの時間自己相関も減少
particleconf$addSampler(target = c("std_xy", "std2_1","std_r", "rho"), type = "RW_PF_block",
                        control = list(adaptive = TRUE, propCov= diag(4), adaptScaleOnly = FALSE,adaptInterval = 50, # adaptScaleOnlyがFalseだとadoptiveがscaleだけじゃなくてpropcovも行う。

                                       pfNparticles = n_par, pfType = "bootstrap",
                                       latents = c("x"),

                                       pfOptimizeNparticles = FALSE))  # adoptIntervalは小さいと意味ない,でも大きすぎてもいけない


###  calc
Particleconf <- buildMCMC(particleconf)  # xyrをxにまとめることでbuildできるようになった！
#cParticleconf <- compileNimble(Particleconf)

cParticleMCMC <- compileNimble(Particleconf, project = cParticle_model, resetFunctions = TRUE)   

iter = 20
burn_in = 0

res = runMCMC(cParticleMCMC, niter = iter, nchains = 1, nburnin = burn_in, summary = TRUE, WAIC = TRUE, setSeed = 1)  ←←←← Aborting error occourred when I did this line.

2024年5月20日月曜日 23:25:18 UTC+9 Perry de Valpine:

[Perry de Valpine]

If you are willing to share your csv files (off-list if you prefer) so that we can run your code, we will try to take a look at it. 

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/1a409aa5-e32c-4caa-90b9-d8dacbca0bd5n%40googlegroups.com.

[kei]

Hi all nimble users!

I sent my R script and data only for Perry.

However, I'll send my data and script that are changed for other users if you're interested in this bug or are kind enough to help me fix bugs.

I keep trying to debug, however, it is still unsolved problem.

I'll be really happy if you give me any suggestions or comments!

Thanks, 

kei

2024年5月24日金曜日 2:15:35 UTC+9 Perry de Valpine:

[Perry de Valpine]

Hello Kei,

Thanks for sharing your code. It was reproducible and you did a nice job adopting the nimbleRcall as a way to catch a problem as soon as it occurred.

Once I was in a browser() with the caught problem in calc_volume2, I took a look at various variables in the model. What I see is that cParticle_model$std_r is negative, which is obviously invalid and results in an invalid value in x[i, 3]. The problem is that we had some old logic in the sampler_RW_PF_block that ran the particle filter even if the parameter proposal resulted in prior logProb of -Inf, i.e. even if the proposal for the parameters was invalid (outside the support of the prior). The newer logic that some time ago we put throughout nimble's samplers for both efficiency and more stability is to calculate the prior logProb first and then not proceed with dependency calculations (or, in this case, the PF) if the parameters being considered (e.g. the proposal parameters) are outside the support of their priors. The samplers in nimbleSMC were delayed in being updated to this more robust logic. As a result, a negative value could be proposed for std_r and then, instead of being rejected immediately, the PF was run even with the invalid proposed parameter. This resulted in invalid simulations in the PF and the resulting problems you encountered. I recognize that has been frustrating!

Please try installing nimbleSMC as follows to obtain a fixed version:

library(devtools)

install_github("nimble-dev/nimbleSMC", ref = "update_RW_PF_block", subdir = "packages/nimbleSMC")

Cheers,

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/0a051b08-cc08-49fd-bf1b-9f006993daeen%40googlegroups.com.

[kei]

Dear Perry and all nimble users

Thanks to your kind advise, I got to remove the terrible error!

Thank you very much!

Thanks,

kei

2024年5月30日木曜日 0:54:12 UTC+9 Perry de Valpine: