Age_OSLC14() ... Compilation error on line 72

[RLumSK]

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]

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]

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:

BayLum/R/Age_OSLC14.R

Lines 363 to 373 in 6db4dc7

	#--- Index preparation
	ind_OSL=which(SampleNature[1,]==1)
	CS_OSL=cumsum(SampleNature[1,])
	ind_C14=which(SampleNature[2,]==1)
	CS_C14=cumsum(SampleNature[2,])
	ind_change=c(1)
	for(i in 2:(Nb_sample-1)){
	if(SampleNature[1,i]!=SampleNature[1,i+1]){
	ind_change=c(ind_change,i)
	}
	}

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]

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

• OSL <--> C-14 <--> OSL has two index changes (recorded as c(1, 2))
• C-14 <--> OSL <--> OSL has only one index change (recorded as c(1) -> the function crashes because subsequent index get 0

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.