crp2a / BayLum

Chronological Bayesian models integrated for Optically Stimulated Luminescence (OSL) Dating
https://crp2a.github.io/BayLum/
GNU General Public License v3.0
9 stars 1 forks source link

Age_OSLC14() ... Compilation error on line 72 #14

Closed RLumSK closed 2 years ago

RLumSK commented 5 years ago

Observed behaviour

Error in rjags::jags.model(file = con, data = dataList, n.chains = n.chains,  : 
 RUNTIME ERROR:
Compilation error on line 72.
Unknown variable k
Either supply values for this variable with the data
or define it on the left hand side of a relation. 

Running mini example

The following example sets is a copy & paste example from the manual. However, here the example exchanges thesamplenature. Now the first sample is the C-14 sample, not the OSL sample (as it is in the original sample).

library(BayLum)
## Load data
# OSL data
data(DATA1,envir = environment())
data(DATA2,envir = environment())
Data <- combine_DataFiles(DATA2,DATA1)

# 14C data
C14Cal <- DATA_C14$C14[1,1]
SigmaC14Cal <- DATA_C14$C14[1,2]
Names <- DATA_C14$Names[1]

# Prior Age
prior=rep(c(1,60),3)
samplenature = matrix(
  data = c(0, 1, 1, 
           1, 0, 0),
  ncol = 3,
  nrow = 2,
  byrow = TRUE
)
SC <- matrix(data=c(1,1,1,0,1,1,0,0,1,0,0,0),ncol=3,nrow=4,byrow=TRUE)

## Age computation of samples
Age <-
  Age_OSLC14(
    DATA = Data,
    Data_C14Cal = C14Cal,
    Data_SigmaC14Cal = SigmaC14Cal,
    SampleNames = c(Names, "GDB5", "GDB3"),
    Nb_sample = 3,
    SampleNature = samplenature,
    PriorAge = prior,
    StratiConstraints = SC,
    Iter = 50,
    n.chains = 2
  )

Solution?

Lines 364 onwards in Age_OSLC14() contains probably the root of the problem. It calculates the variable ind_change, which is preset to 1 and expresses the index change in the sample nature matrix. If like in the example above, C-14 comes first, later the function jumps into the wrong if expression and sets the object q to 0. If that happens, a loop in the BUGS model fails, and the variable k gets never set. Unfortunately, it appears that this is a bug that is not straightforward to fix. Luckily it seems to disappear when the number of samples is increased to at least four. Nevertheless, given the odd code, the function may cause other unwanted effects, and a fix is strongly recommended.

guiguerin commented 4 years ago

In the example above there is an inconsistency between:

'samplenature = matrix( data = c(0, 1, 1, 1, 0, 0)'...

which means C14, OSL, OSL

and Age <- Age_OSLC14( DATA = Data, Data_C14Cal = C14Cal, Data_SigmaC14Cal = SigmaC14Cal, SampleNames = c("GDB5", Names, "GDB3"),

where the name order corresponds to OSL, C14, OSL.

In other words can you reproduce the error by replacing the last line by SampleNames = c(Names, "GDB5", "GDB3")? Could this explain the issue? If not, then perhaps the function needs at least 2 C14 (and/or 2 OSL) samples to run...

RLumSK commented 4 years ago

Unfortunately, this has no effect. Though, I tested it and corrected the example above to avoid this misunderstanding. As written in the initial analysis, the problem starts all here: https://github.com/crp2a/BayLum/blob/6db4dc7e785261656cd962e28e0030c88d280152/R/Age_OSLC14.R#L363-L373

I will have another look into the code.

Update

The function requires at least two changes. In other words, the example will also fail in the case when the C-14 sample is at the end. However, if the sample is in the middle, we have two changes and it will still work.

RLumSK commented 4 years ago

I had a thorough look at it. The function expects at least two index changes for the sample nature. Example:

Oddly enough, OSL <--> OSL <--> C-14 has only one index change but is recorded as c(1,3).

After giving it some thought, I found it easier and safer to opt for a proper error message instead of a fix, because the index counting extends into the BUGS model. In other words, more code would need attention. I found that too risky because I do not overlook all consequences and it is just an edge case. The new error message should help the user to understand the underlying problem.

If ok for you, I would close the issue.