Bounds Error with Multiple dosing

I am running into the error:

BoundsError: attempt to access 2-element Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}} at index [3]

When fitting one subject that has multiple doses at steady state.
This is what their dataframe looks like:


When I change both SS to 1, the fitting does not give the error, but the fit is not correct since the compartments are initialized at the second dose. I have been able to fit this subject in the past, but it is not working now.
Here is the model I am using:

model_2stage_indiv = @model begin
    @param begin
        i_dur      ∈ RealDomain(lower=0.01, upper=24)
        i_ka_slow  ∈ RealDomain(lower=0.01)
        i_vc_F     ∈ RealDomain(lower=0.1)
        i_cl_F     ∈ RealDomain(lower=0.1)
        #RUV
        i_σ²_add   ∈ RealDomain(lower=0.0001)
        i_σ²_prop  ∈ RealDomain(lower=0.0001) 
    end
    @pre begin  
        dur      = i_dur
        ka_slow  = i_ka_slow
        vc_F     = i_vc_F
        cl_F     = i_cl_F
        σ²_add   = i_σ²_add
        σ²_prop  = i_σ²_prop
    end
    @dosecontrol begin
        duration = (;Depot = i_dur)
    end
    @dynamics begin 
        Depot'   = -ka_slow*Depot
        Central' =  ka_slow*Depot - Central*(cl_F/vc_F)       
    end
    @derived begin
        Cp  = @. Central/vc_F
        CONC ~ @. Normal(Cp, sqrt(((Cp^2)*i_σ²_prop) + i_σ²_add))
    end
end

Hi!

Based on your error message, it seems you do not use the latest release of Pumas. Could you try rerunning your analysis with Pumas 2.6 if you have access to it?

Based on the data you showed I tried to run the following script:

using Pumas, DataFrames, PumasPlots

model_2stage_indiv = @model begin
    @param begin
        i_dur ∈ RealDomain(lower = 0.01, upper = 24)
        i_ka_slow ∈ RealDomain(lower = 0.01)
        i_vc_F ∈ RealDomain(lower = 0.1)
        i_cl_F ∈ RealDomain(lower = 0.1)
        #RUV
        i_σ²_add ∈ RealDomain(lower = 0.0001)
        i_σ²_prop ∈ RealDomain(lower = 0.0001)
    end
    @pre begin
        dur = i_dur
        ka_slow = i_ka_slow
        vc_F = i_vc_F
        cl_F = i_cl_F
        σ²_add = i_σ²_add
        σ²_prop = i_σ²_prop
    end
    @dosecontrol begin
        duration = (; Depot = i_dur)
    end
    @dynamics begin
        Depot' = -ka_slow * Depot
        Central' = ka_slow * Depot - Central * (cl_F / vc_F)
    end
    @derived begin
        Cp = @. Central / vc_F
        CONC ~ @. Normal(Cp, sqrt(((Cp^2) * i_σ²_prop) + i_σ²_add))
    end
end

# Population (single subject)
df = DataFrame(;
    ID = fill("sub", 22),
    TIME = [
        0,
        0.0001,
        0.5,
        0.8,
        1,
        1.5,
        2,
        3,
        4,
        5.5,
        5.5001,
        6,
        6.3,
        6.5,
        7,
        7.5,
        8.5,
        9.5,
        11.5,
        13.5,
        17.5,
        21.5,
    ],
    CONC = [
        missing,
        2.207877,
        2.652662,
        5.760089,
        56.12353,
        104.4228,
        134.7264,
        127.8454,
        95.49685,
        missing,
        71.41494,
        57.27378,
        106.7809,
        178.8866,
        172.984,
        164.1932,
        131.9889,
        104.5845,
        65.16986,
        30.63678,
        11.5158,
        4.608377,
    ],
    EVID = [1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    AMT = [
        18000,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        18000,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
    ],
    CMT = [
        1,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        1,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
    ],
    RATE = [
        -2,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        -2,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
    ],
    SS = [
        1,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        2,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
    ],
    II = [
        24,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        24,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
    ],
)
pop = read_pumas(
    df;
    id = :ID,
    time = :TIME,
    observations = [:CONC],
    evid = :EVID,
    amt = :AMT,
    cmt = :CMT,
    rate = :RATE,
    ss = :SS,
    ii = :II,
)

# Fit the model
ft = fit(model_2stage_indiv, pop, init_params(model_2stage_indiv), NaivePooled())
sim_plot(ft)

On Pumas 2.5.1, it throws a BoundsError, as you reported above. On Pumas 2.6.1 and our development version, however, I was able to fit the model successfully: