text:
println("===============================================")
println("Testing processing.jl")
println("===============================================")

function _resample_timeonly()
    d = timeline(t_start=0, t_end=1, step=1//4)
    d = resample(d,dt=0.1)
    (length(d.t)==11) && (d.t[2]==0.1)
end
@test _resample_timeonly()

function _resample_stressonly()
    d = RheoTimeData(σ=[1,2,4,7], t=[0,1,3,6])
    d = resample(d) # make timestep uniform without changing number of points
    (length(d.t)==4) && (d.t[2]==2) && (d.σ[2]≈3)

end
@test _resample_stressonly()

function _resample_strainonly()
    d = RheoTimeData(ϵ=[1,2,4,7], t=[0,1,3,6])
    d = resample(d,t=1:1//4:2)
    (length(d.t)==5) && (d.t[2]==5//4) && (d.ϵ[2]≈9//4)
end
@test _resample_strainonly()

function _resample_stressstrain()
    d = RheoTimeData(ϵ=[1,2,4,7], t=[0,1,3,6], σ=[7,4,2,1])
    d = resample(d,scale=2)
    (length(d.t)==7) && (d.t[4]≈1.875) && (d.σ[4]≈2.6328125) && (d.ϵ[4]≈2.875)
end
@test _resample_stressstrain()

# function _resample_strainonly()
#     t0 = collect(0.0:0.01:1.0)
#     ϵ0 = t0.^2
#
#     t1 = collect(0.0:0.1:1.0)
#     ϵ1 = t1.^2
#
#     data0 = RheoTimeData(t = t0, ϵ = ϵ0)
#
#     dataout = resample(data0, -10)
#
#     dataout.t==t1 && dataout.ϵ==ϵ1
# end
# @test _resample_strainonly()
#
# function _resample_stressonly()
#     t0 = collect(0.0:0.01:1.0)
#     σ0 = t0.^2
#
#     t1 = collect(0.0:0.1:1.0)
#     σ1 = t1.^2
#
#     data0 = RheoTimeData(t = t0, σ = σ0)
#
#     dataout = resample(data0, -10)
#
#     dataout.t==t1 && dataout.σ==σ1
# end
# @test _resample_stressonly()
#
# function _resample_stressandstrain()
#     t0 = collect(0.0:0.01:1.0)
#     ϵ0 = t0.^2
#     σ0 = t0.^3
#
#     t1 = collect(0.0:0.1:1.0)
#     ϵ1 = t1.^2
#     σ1 = t1.^3
#
#     data0 = RheoTimeData(t = t0, σ = σ0, ϵ = ϵ0)
#
#     dataout = resample(data0, -10)
#
#     dataout.t==t1 && dataout.ϵ==ϵ1 && dataout.σ==σ1
# end
# @test _resample_stressandstrain()
#
# function _resample_stressandstrain_multiplesections()
#     t0 = collect(-10.0:0.5:10.0)
#     ϵ0 = 2*t0
#     σ0 = 2*t0
#
#     # data for comparison composed
#     # of multiple sections
#     t1a = collect(-10.0:0.1:-5.1)
#     t1b = collect(-5.0:1.0:0.0)
#     t1c = collect(0.05:0.05:4.95)
#     t1d = collect(5.0:2.0:10.0)
#     t1e = [t0[end]]
#
#     t1 = vcat(t1a, t1b, t1c, t1d, t1e)
#     ϵ1 = 2*t1
#     σ1 = 2*t1
#
#     data0 = RheoTimeData(t = t0, σ = σ0, ϵ = ϵ0)
#
#     dataout = resample(data0, [5, -2, 10, -4]; time_boundaries = [-10.0, -5.0, 0.0, 5.0, 10.0])
#
#     dataout.t==t1 && dataout.ϵ==ϵ1 && dataout.σ==σ1
# end
# @test _resample_stressandstrain_multiplesections()

function _indexweight_oneregiondownsample_includelastel()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, -2)

    indices == [1, 3, 5, 7, 9, 11]
end
@test _indexweight_oneregiondownsample_includelastel()

function _indexweight_oneregiondownsample_nolastel()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, -2; includelast=false)

    indices == [1, 3, 5, 7, 9]
end
@test _indexweight_oneregiondownsample_nolastel()

function _indexweight_oneregionupsample_includelastel()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, 2)

    indices == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11]
end
@test _indexweight_oneregionupsample_includelastel()

function _indexweight_oneregionupsample_nolastel()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, 2; includelast=false)

    indices == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10]
end
@test _indexweight_oneregionupsample_nolastel()

function _indexweight_tworegiondownsample_includelastel()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, [-2, -3]; time_boundaries=[0.0, 5.0, 10.0])

    indices == [1, 3, 5, 6, 9, 11]
end
@test _indexweight_tworegiondownsample_includelastel()

function _indexweight_tworegiondownsample_nolastel()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, [-2, -3]; time_boundaries=[0.0, 5.0, 10.0], includelast=false)

    indices == [1, 3, 5, 6, 9]
end
@test _indexweight_tworegiondownsample_nolastel()

function _indexweight_tworegionupsample_includelastel()
    timedata = timeline(t_start=0.0, t_end=5.0, step=1.0)

    indices = indexweight(timedata, [2, 3]; time_boundaries=[0.0, 3.0, 4.0])

    indices == [1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5]
end
@test _indexweight_tworegionupsample_includelastel()

function _indexweight_tworegionupsample_nolastel()
    timedata = timeline(t_start=0.0, t_end=5.0, step=1.0)

    indices = indexweight(timedata, [2, 3]; time_boundaries=[0.0, 3.0, 4.0], includelast=false)

    indices == [1, 1, 2, 2, 3, 3, 4, 4, 4]
end
@test _indexweight_tworegionupsample_nolastel()

function _indexweight_mixed_includelastel_finishup()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, [-2, 3, -3, 4]; time_boundaries=[0.0, 3.0, 5.0, 9.0, 10.0])

    indices == [1, 3, 4, 4, 4, 5, 5, 5, 6, 9, 10, 10, 10, 10, 11, 11, 11, 11]
end
@test _indexweight_mixed_includelastel_finishup()

function _indexweight_mixed_includelastel_finishdown()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, [2, -3, 2, -2]; time_boundaries=[0.0, 3.0, 5.0, 7.0, 10.0])

    indices == [1, 1, 2, 2, 3, 3, 4, 6, 6, 7, 7, 8, 10, 11]
end
@test _indexweight_mixed_includelastel_finishdown()

function _indexweight_mixed_nolastel_finishup()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, [-2, 3, -3, 4]; time_boundaries=[0.0, 3.0, 5.0, 9.0, 10.0], includelast=false)

    indices == [1, 3, 4, 4, 4, 5, 5, 5, 6, 9, 10, 10, 10, 10]
end
@test _indexweight_mixed_nolastel_finishup()

function _indexweight_mixed_nolastel_finishdown()
    timedata = timeline(t_start=0.0, t_end=10.0, step=1.0)

    indices = indexweight(timedata, [2, -3, 2, -2]; time_boundaries=[0.0, 3.0, 5.0, 7.0, 10.0], includelast=false)

    indices == [1, 1, 2, 2, 3, 3, 4, 6, 6, 7, 7, 8, 10]
end
@test _indexweight_mixed_nolastel_finishdown()

function _cutting_strainonly()
    t0 = collect(0.0:1.0:5.0)
    ϵ0 = t0.^2

    t1 = collect(1.0:1.0:5.0)
    ϵ1 = t1.^2

    data0 = RheoTimeData(t = t0, ϵ = ϵ0)

    dataout = cutting(data0, 1.0, 5.0)

    dataout.t==t1 && dataout.ϵ==ϵ1
end
@test _cutting_strainonly()

function _cutting_stressonly()
    t0 = collect(0.0:1.0:5.0)
    σ0 = t0.^2

    t1 = collect(1.0:1.0:5.0)
    σ1 = t1.^2

    data0 = RheoTimeData(t = t0, σ = σ0)

    dataout = cutting(data0, 1.0, 5.0)

    dataout.t==t1 && dataout.σ==σ1
end
@test _cutting_stressonly()

function _cutting_stressandstrain()
    t0 = collect(0.0:1.0:5.0)
    ϵ0 = t0.^2
    σ0 = t0.^3

    t1 = collect(1.0:1.0:5.0)
    ϵ1 = t1.^2
    σ1 = t1.^3

    data0 = RheoTimeData(t = t0, σ = σ0, ϵ = ϵ0)

    dataout = cutting(data0, 1.0, 5.0)

    dataout.t==t1 && dataout.ϵ==ϵ1 && dataout.σ==σ1
end
@test _cutting_stressandstrain()

function _smooth_stressonly()
    f = 1.0 # Hz
    ω = 2*π*f # rad/s
    t0 = collect(0.0:0.01:4/f)
    σ0 = sin.(ω*t0)

    data0 = RheoTimeData(t = t0, σ = σ0)

    dataout = smooth(data0, 1/f; pad="circular")

    all(i -> isapprox(0.5*data0.σ[i], dataout.σ[i], atol=0.1), eachindex(data0.σ))
end
@test _smooth_stressonly()

function _smooth_strainonly()
    f = 1.0 # Hz
    ω = 2*π*f # rad/s
    t0 = collect(0.0:0.01:4/f)
    ϵ0 = sin.(ω*t0)

    data0 = RheoTimeData(t = t0, ϵ = ϵ0)

    dataout = smooth(data0, 1/f; pad="circular")

    all(i -> isapprox(0.5*data0.ϵ[i], dataout.ϵ[i], atol=0.1), eachindex(data0.ϵ))
end
@test _smooth_strainonly()

function _smooth_stressandstrain()
    f = 1.0 # Hz
    ω = 2*π*f # rad/s
    t0 = collect(0.0:0.01:4/f)
    ϵ0 = sin.(ω*t0)
    σ0 = sin.(ω*t0)

    data0 = RheoTimeData(t = t0, ϵ = ϵ0, σ = σ0)

    dataout = smooth(data0, 1/f; pad="circular")

    test1 = all(i -> isapprox(0.5*data0.ϵ[i], dataout.ϵ[i], atol=0.1), eachindex(data0.ϵ))
    test2 = all(i -> isapprox(0.5*data0.σ[i], dataout.σ[i], atol=0.1), eachindex(data0.σ))

    test1 && test2
end
@test _smooth_stressandstrain()

function _extract_timefromtimedata()
    f = 1.0 # Hz
    ω = 2*π*f # rad/s
    t0 = collect(0.0:0.01:4/f)
    ϵ0 = sin.(ω*t0)
    σ0 = sin.(ω*t0)

    data0 = RheoTimeData(t = t0, ϵ = ϵ0, σ = σ0)

    dataout = extract(data0, time_only)

    data0.t==dataout.t && dataout.σ==[] && dataout.ϵ==[]
end
@test _extract_timefromtimedata()

function _extract_timestressfromtimedata()
    f = 1.0 # Hz
    ω = 2*π*f # rad/s
    t0 = collect(0.0:0.01:4/f)
    ϵ0 = sin.(ω*t0)
    σ0 = sin.(ω*t0)

    data0 = RheoTimeData(t = t0, ϵ = ϵ0, σ = σ0)

    dataout = extract(data0, stress_only)

    data0.t==dataout.t && data0.σ==dataout.σ && dataout.ϵ==[]
end
@test _extract_timestressfromtimedata()

function _extract_timestrainfromtimedata()
    f = 1.0 # Hz
    ω = 2*π*f # rad/s
    t0 = collect(0.0:0.01:4/f)
    ϵ0 = sin.(ω*t0)
    σ0 = sin.(ω*t0)

    data0 = RheoTimeData(t = t0, ϵ = ϵ0, σ = σ0)

    dataout = extract(data0, strain_only)

    data0.t==dataout.t && dataout.σ==[] && data0.ϵ==dataout.ϵ
end
@test _extract_timestrainfromtimedata()

function _extract_freqfromfreqdata()
    ω0 = collect(0.0:0.01:10.0)
    Gp0 = 2*ω0
    Gpp0 = 3*ω0

    data0 = RheoFreqData(ω = ω0, Gp = Gp0, Gpp = Gpp0)

    dataout = extract(data0, freq_only)

    data0.ω==dataout.ω && dataout.Gp==[] && dataout.Gpp==[]
end
@test _extract_freqfromfreqdata()

function _fill_init_params_nothing()
    RHEOS.fill_init_params(Fract_Zener, nothing)==[0.5, 0.8, 0.5, 0.5, 0.5, 0.5]
end
@test _fill_init_params_nothing()

function _fill_init_params_fullyspecified()
    RHEOS.fill_init_params(Fract_Zener, (cₐ=1.0, a=0.51, cᵦ=1.5, β=0.4, cᵧ=2.0, γ=0.7))==[1.0, 0.51, 1.5, 0.4, 2.0, 0.7]
end
@test _fill_init_params_fullyspecified()

function _fill_init_params_semispecified1()
    RHEOS.fill_init_params(Fract_Zener, (a=0.51, cᵦ=1.5, β=0.4, cᵧ=2.0, γ=0.7))==[0.5, 0.51, 1.5, 0.4, 2.0, 0.7]
end
@test _fill_init_params_semispecified1()

function _fill_init_params_semispecified2()
    RHEOS.fill_init_params(Fract_Zener, (cᵦ=1.5, β=0.51, cᵧ=2.0, γ=0.7))==[0.5, 0.999, 1.5, 0.51, 2.0, 0.7]
end
@test _fill_init_params_semispecified2()

function _fill_init_params_semispecified3()
    RHEOS.fill_init_params(Fract_Zener, (cᵦ=1.5, β=0.4, cᵧ=2.0, γ=0.7))==[0.5, 0.8, 1.5, 0.4, 2.0, 0.7]
end
@test _fill_init_params_semispecified3()

function _fill_init_params_semispecified4()
    RHEOS.fill_init_params(Fract_Zener, (cₐ=1.0, a=0.4, cᵦ=1.5, cᵧ=2.0, γ=0.7))==[1.0, 0.4, 1.5, 0.2, 2.0, 0.7]
end
@test _fill_init_params_semispecified4()

function _fill_init_params_semispecified5()
    RHEOS.fill_init_params(Fract_Zener, (cₐ=1.0, cᵦ=1.5, cᵧ=2.0, γ=0.7))==[1.0, 0.8, 1.5, 0.5, 2.0, 0.7]
end
@test _fill_init_params_semispecified5()

function _fill_lower_bounds_nothing()
    isnothing(RHEOS.fill_lower_bounds(Fract_Zener, nothing))
end
@test _fill_lower_bounds_nothing()

function _fill_lower_bounds_fullyspecified()
    RHEOS.fill_lower_bounds(Fract_Zener, (cₐ=1.0, a=0.4, cᵦ=1.5, β=0.51, cᵧ=2.0, γ=0.7))==[1.0, 0.4, 1.5, 0.51, 2.0, 0.7]
end
@test _fill_lower_bounds_fullyspecified()

function _fill_lower_bounds_semispecified()
    RHEOS.fill_lower_bounds(Fract_Zener, (a=0.4, β=0.51, cᵧ=2.0, γ=0.7))==[0.0, 0.4, 0.0, 0.51, 2.0, 0.7]
end
@test _fill_lower_bounds_semispecified()

function _fill_upper_bounds_nothing()
    isnothing(RHEOS.fill_upper_bounds(Fract_Zener, nothing))
end
@test _fill_upper_bounds_nothing()

function _fill_upper_bounds_fullyspecified()
    RHEOS.fill_upper_bounds(Fract_Zener, (cₐ=1.0, a=0.4, cᵦ=1.5, β=0.51, cᵧ=2.0, γ=0.7))==[1.0, 0.4, 1.5, 0.51, 2.0, 0.7]
end
@test _fill_upper_bounds_fullyspecified()

function _fill_upper_bounds_semispecified()
    RHEOS.fill_upper_bounds(Fract_Zener, (a=0.4, β=0.51, cᵧ=2.0, γ=0.7))==[Inf, 0.4, Inf, 0.51, 2.0, 0.7]
end
@test _fill_upper_bounds_semispecified()

function _modelfit_const_ramp_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = ramp_loading, σ = exact_response)

    init_params = (α=1.0, β=1.0)
    modelout = modelfit(data0, model, strain_imposed, p0=init_params, lo=(α=0.9, β=0.9), hi=(α=1.1, β=1.1))

    found_params = modelout.params
    isapprox(collect(values(found_params)), collect(values(init_params)), atol = tol)
end
@test _modelfit_const_ramp_relax(tol)

function _modelfit_const_ramp_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    init_params = (α=1.0, β=1.0)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.9, β=0.9), hi=(α=1.1, β=1.1))

    found_params = modelout.params
    isapprox(collect(values(found_params)), collect(values(init_params)), atol = tol)
end
@test _modelfit_const_ramp_creep(tol)

function _modelfit_const_ramp_nobounds_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = ramp_loading, σ = exact_response)

    init_params = (α=1.0, β=1.0)
    modelout = modelfit(data0, model, strain_imposed, p0=init_params)

    found_params = modelout.params
    isapprox(collect(values(found_params)), collect(values(init_params)), atol = tol)
end
@test _modelfit_const_ramp_nobounds_relax(tol)

function _modelfit_const_ramp_nobounds_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    init_params = (α=1.0, β=1.0)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params)

    found_params = modelout.params
    isapprox(collect(values(found_params)), collect(values(init_params)), atol = tol)
end
@test _modelfit_const_ramp_nobounds_creep(tol)

function _modelfit_const_ramp_nobounds_singleparam_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t) end
    model = RheoModelClass(name = "testmodel", p = [:α], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = ramp_loading, σ = exact_response)

    init_params = (α=1.0,)
    modelout = modelfit(data0, model, strain_imposed, p0=init_params)

    found_params = modelout.params
    isapprox(collect(values(found_params)), collect(values(init_params)), atol = tol)
end
@test _modelfit_const_ramp_nobounds_singleparam_relax(tol)

function _modelfit_const_ramp_nobounds_singleparam_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t) end
    model = RheoModelClass(name = "testmodel", p = [:α], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    init_params = (α=1.0,)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params)

    found_params = modelout.params
    isapprox(collect(values(found_params)), collect(values(init_params)), atol = tol)
end
@test _modelfit_const_ramp_nobounds_singleparam_creep(tol)

function _modelfit_const_ramp_nobounds_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = ramp_loading, σ = exact_response)

    init_params = (α=1.0, β=1.0)
    modelout = modelfit(data0, model, strain_imposed, p0=init_params)

    found_params = modelout.params
    isapprox(collect(values(found_params)), collect(values(init_params)), atol = tol)
end
@test _modelfit_const_ramp_nobounds_relax(tol)

function _modelfit_const_ramp_nobounds_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    init_params = (α=1.0, β=1.0)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params)

    found_params = modelout.params
    isapprox(collect(values(found_params)), collect(values(init_params)), atol = tol)
end
@test _modelfit_const_ramp_nobounds_creep(tol)

function _modelfit_var_ramp_nonsing_relax(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = ramp_loading, σ = exact_response)

    init_params = (α=1.3, β=0.7)
    modelout = modelfit(data0, model, strain_imposed, p0=init_params, lo=(α=0.5, β=0.5), hi=(α=1.5, β=1.5))

    found_params = modelout.params
    actual_params = (α=1.0, β=1.0)

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _modelfit_var_ramp_nonsing_relax(tol)

function _modelfit_var_ramp_nonsing_creep(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    init_params = (α=1.3, β=0.7)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.5, β=0.5), hi=(α=1.5, β=1.5))

    found_params = modelout.params
    actual_params = (α=1.0, β=1.0)

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _modelfit_var_ramp_nonsing_creep(tol)

function _modelfit_const_ramp_sing_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=1.0, β=0.5)
    exact_response = actual_params.α*t.^(1.0 - actual_params.β) / (1.0 - actual_params.β)
    ramp_loading = t

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = ramp_loading, σ = exact_response)

    init_params = (α=1.3, β=0.7)
    modelout = modelfit(data0, model, strain_imposed, p0=init_params, lo=(α=0.5, β=0.2), hi=(α=1.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = 5*tol)
end
@test _modelfit_const_ramp_sing_relax(tol)

function _modelfit_const_ramp_sing_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=1.0, β=0.5)
    exact_response = actual_params.α*t.^(1.0 - actual_params.β) / (1.0 - actual_params.β)
    ramp_loading = t

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    init_params = (α=1.3, β=0.7)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.5, β=0.2), hi=(α=1.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = 5*tol)
end
@test _modelfit_const_ramp_sing_creep(tol)

function _modelfit_var_ramp_sing_relax(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    actual_params = (α=1.0, β=0.5)
    exact_response = actual_params.α*t.^(1.0 - actual_params.β) / (1.0 - actual_params.β)
    ramp_loading = t

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = ramp_loading, σ = exact_response)

    init_params = (α=1.3, β=0.7)
    modelout = modelfit(data0, model, strain_imposed, p0=init_params, lo=(α=0.5, β=0.2), hi=(α=1.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = 5*tol)
end
@test _modelfit_var_ramp_sing_relax(tol)

function _modelfit_const_ramp_creep_weighted_selftest(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    actual_params = (α=1.0, β=1.0)
    init_params = (α=1.2, β=0.8)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.7, β=0.7), hi=(α=1.5, β=1.5), weights=collect(Integer, 1:length(t)))

    found_params = modelout.params
    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _modelfit_const_ramp_creep_weighted_selftest(tol)

function _modelfit_const_ramp_sing_creep_weighted_selftest(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=1.0, β=0.5)
    exact_response = actual_params.α*t.^(1.0 - actual_params.β) / (1.0 - actual_params.β)
    ramp_loading = t

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    init_params = (α=1.3, β=0.7)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.5, β=0.2), hi=(α=1.5, β=1.5), weights=collect(Integer, 1:length(t)))

    found_params = modelout.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = 5*tol)
end
@test _modelfit_const_ramp_sing_creep_weighted_selftest(tol)

function _modelfit_const_ramp_creep_weighted_downsampled(tol)
    dt = 0.005
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0 .- exp.(-t)
    ramp_loading = t

    modulus = quote α*exp(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    indices = indexweight(data0, -100)

    init_params = (α=1.2, β=0.8)
    modelout_weighted = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.7, β=0.7), hi=(α=1.5, β=1.5), weights=indices)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.7, β=0.7), hi=(α=1.5, β=1.5))

    found_params_weighted = modelout_weighted.params
    found_params = modelout.params
    isapprox(collect(values(found_params_weighted)), collect(values(found_params)), atol = tol)
end
@test _modelfit_const_ramp_creep_weighted_downsampled(tol)

function _modelfit_const_ramp_sing_creep_weighted_downsampled(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=1.0, β=0.5)
    exact_response = actual_params.α*t.^(1.0 - actual_params.β) / (1.0 - actual_params.β)
    ramp_loading = t

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    indices = indexweight(data0, -100)

    init_params = (α=1.3, β=0.7)
    modelout_weighted = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.5, β=0.2), hi=(α=1.5, β=1.5), weights=indices)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.5, β=0.2), hi=(α=1.5, β=1.5))

    found_params_weighted = modelout_weighted.params
    found_params = modelout.params
    isapprox(collect(values(found_params_weighted)), collect(values(found_params)), atol = tol)
end
@test _modelfit_const_ramp_sing_creep_weighted_downsampled(tol)

function _modelfit_var_ramp_sing_creep(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    actual_params = (α=1.0, β=0.5)
    exact_response = actual_params.α*t.^(1.0 - actual_params.β) / (1.0 - actual_params.β)
    ramp_loading = t

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = ramp_loading)

    init_params = (α=1.3, β=0.7)
    modelout = modelfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.5, β=0.2), hi=(α=1.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = 5*tol)
end
@test _modelfit_var_ramp_sing_creep(tol)

function _modelpredict_nonsing_const_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1 .- exp.(-t)
    loading = t

    modulus = quote α*exp.(-t/β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=1.0)
    data0 = RheoTimeData(t = t, ϵ = loading)

    computed_response = modelpredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.σ[i], atol=tol), eachindex(exact_response))
end
@test _modelpredict_nonsing_const_relax(tol)

function _modelpredict_nonsing_const_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1 .- exp.(-t)
    loading = t

    modulus = quote α*exp.(-t/β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=1.0)
    data0 = RheoTimeData(t = t, σ = loading)

    computed_response = modelpredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.ϵ[i], atol=tol), eachindex(exact_response))
end
@test _modelpredict_nonsing_const_creep(tol)

function _modelpredict_nonsing_var_relax(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    exact_response = 1 .- exp.(-t)
    loading = t

    modulus = quote α*exp.(-t/β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=1.0)
    data0 = RheoTimeData(t = t, ϵ = loading)

    computed_response = modelpredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.σ[i], atol=tol), eachindex(exact_response))
end
@test _modelpredict_nonsing_var_relax(tol)

function _modelpredict_nonsing_var_creep(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    exact_response = 1 .- exp.(-t)
    loading = t

    modulus = quote α*exp.(-t/β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=1.0)
    data0 = RheoTimeData(t = t, σ = loading)

    computed_response = modelpredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.ϵ[i], atol=tol), eachindex(exact_response))
end
@test _modelpredict_nonsing_var_creep(tol)

function _modelpredict_sing_const_relax(tol)
    dt = 0.001
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)
    loading = t

    modulus = quote α*t.^(-β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=0.5)
    data0 = RheoTimeData(t = t, ϵ = loading)

    computed_response = modelpredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.σ[i], atol=5*tol), eachindex(exact_response))
end
@test _modelpredict_sing_const_relax(tol)

function _modelpredict_sing_const_creep(tol)
    dt = 0.001
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)
    loading = t

    modulus = quote α*t.^(-β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=0.5)
    data0 = RheoTimeData(t = t, σ = loading)

    computed_response = modelpredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.ϵ[i], atol=5*tol), eachindex(exact_response))
end
@test _modelpredict_sing_const_creep(tol)

function _modelpredict_sing_var_relax(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)
    loading = t

    modulus = quote α*t.^(-β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=0.5)
    data0 = RheoTimeData(t = t, ϵ = loading)

    computed_response = modelpredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.σ[i], atol=5*tol), eachindex(exact_response))
end
@test _modelpredict_sing_var_relax(tol)

function _modelpredict_sing_var_creep(tol)
    dt = 0.01
    t = [Vector{RheoFloat}(0.0:dt:(15.0-dt)); Vector{RheoFloat}(15.0:10*dt:20.0)]
    exact_response = t.^(1.0 - 0.5) / (1.0 - 0.5)
    loading = t

    modulus = quote α*t.^(-β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=0.5)
    data0 = RheoTimeData(t = t, σ = loading)

    computed_response = modelpredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.ϵ[i], atol=5*tol), eachindex(exact_response))
end
@test _modelpredict_sing_var_creep(tol)

function _modelstepfit_nonsing_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*exp.(-t/actual_params[2])

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = loading, σ = exact_response)

    init_params = (α=1.3, β=0.7)
    modelout = modelstepfit(data0, model, strain_imposed, p0=init_params, lo=(α=0.2, β=0.2), hi=(α=3.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_nonsing_relax(tol)

function _modelstepfit_nonsing_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*exp.(-t/actual_params[2])

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = loading)

    init_params = (α=1.3, β=0.7)
    modelout = modelstepfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.2, β=0.2), hi=(α=3.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_nonsing_creep(tol)

function _modelstepfit_nonsing_noinitparams_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*exp.(-t/actual_params[2])

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = loading, σ = exact_response)

    modelout = modelstepfit(data0, model, strain_imposed, lo=(α=0.2, β=0.2), hi=(α=3.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_nonsing_noinitparams_relax(tol)

function _modelstepfit_nonsing_noinitparams_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*exp.(-t/actual_params[2])

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = loading)

    modelout = modelstepfit(data0, model, stress_imposed, lo=(α=0.2, β=0.2), hi=(α=3.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_nonsing_noinitparams_creep(tol)

function _modelstepfit_nonsing_nobounds_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*exp.(-t/actual_params[2])

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = loading, σ = exact_response)

    init_params = (α=1.3, β=0.7)
    modelout = modelstepfit(data0, model, strain_imposed, p0=init_params)

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_nonsing_nobounds_relax(tol)

function _modelstepfit_nonsing_nobounds_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*exp.(-t/actual_params[2])

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = loading)

    init_params = (α=1.3, β=0.7)
    modelout = modelstepfit(data0, model, stress_imposed, p0=init_params)

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_nonsing_nobounds_creep(tol)

function _modelstepfit_sing_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*t.^(-actual_params[2])

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = loading, σ = exact_response)

    init_params = (α=1.3, β=0.7)
    modelout = modelstepfit(data0, model, strain_imposed, p0=init_params, lo=(α=0.2, β=0.2), hi=(α=3.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_sing_relax(tol)

function _modelstepfit_sing_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*t.^(-actual_params[2])

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = exact_response, σ = loading)

    init_params = (α=1.3, β=0.7)
    modelout = modelstepfit(data0, model, stress_imposed, p0=init_params, lo=(α=0.2, β=0.2), hi=(α=3.5, β=1.5))

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_sing_creep(tol)

function _modelstepfit_nonsing_nobounds_relax_weighted_selftest(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*exp.(-t/actual_params[2])

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = loading, σ = exact_response)
    indices = indexweight(data0, -1)

    init_params = (α=1.3, β=0.7)
    modelout = modelstepfit(data0, model, strain_imposed; p0=init_params, weights=indices)

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_nonsing_nobounds_relax_weighted_selftest(tol)

function _modelstepfit_sing_relax_weighted_selftest(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*t.^(-actual_params[2])

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = loading, σ = exact_response)
    indices = indexweight(data0, -1)

    init_params = (α=1.3, β=0.7)
    modelout = modelstepfit(data0, model, strain_imposed, p0=init_params, lo=(α=0.2, β=0.2), hi=(α=3.5, β=1.5), weights=indices)

    found_params = modelout.params

    isapprox(collect(values(actual_params)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_sing_relax_weighted_selftest(tol)

function _modelstepfit_nonsing_nobounds_relax_weighted_closeness(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*exp.(-t/actual_params[2])

    modulus = quote α*exp.(-t/β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = loading, σ = exact_response)
    indices = indexweight(data0, -100)

    init_params = (α=1.3, β=0.7)
    modelout_weighted = modelstepfit(data0, model, strain_imposed; p0=init_params, weights=indices)
    modelout = modelstepfit(data0, model, strain_imposed; p0=init_params)

    found_params_weighted = modelout_weighted.params
    found_params = modelout.params

    isapprox(collect(values(found_params_weighted)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_nonsing_nobounds_relax_weighted_closeness(tol)

function _modelstepfit_sing_relax_weighted_closeness(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    actual_params = (α=2.0, β=0.5)
    loading = ones(length(t))
    exact_response = loading[1]*actual_params[1]*t.^(-actual_params[2])

    modulus = quote α*t.^(-β) end
    model = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")

    data0 = RheoTimeData(t = t, ϵ = loading, σ = exact_response)
    indices = indexweight(data0, -100)

    init_params = (α=1.3, β=0.7)
    modelout_weighted = modelstepfit(data0, model, strain_imposed, p0=init_params, lo=(α=0.2, β=0.2), hi=(α=3.5, β=1.5), weights=indices)
    modelout = modelstepfit(data0, model, strain_imposed, p0=init_params, lo=(α=0.2, β=0.2), hi=(α=3.5, β=1.5))

    found_params_weighted = modelout_weighted.params
    found_params = modelout.params

    isapprox(collect(values(found_params_weighted)), collect(values(found_params)), atol=tol)
end
@test _modelstepfit_sing_relax_weighted_closeness(tol)

function _modelsteppredict_nonsing_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0*exp.(-t/1.0)
    loading = ones(RheoFloat, length(t))

    modulus = quote α*exp.(-t/β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=1.0)

    data0 = RheoTimeData(t = t, ϵ = loading)

    computed_response = modelsteppredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.σ[i], atol=tol), eachindex(exact_response))
end
@test _modelsteppredict_nonsing_relax(tol)

function _modelsteppredict_nonsing_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = 1.0*exp.(-t/1.0)
    loading = ones(RheoFloat, length(t))

    modulus = quote α*exp.(-t/β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=1.0)

    data0 = RheoTimeData(t = t, σ = loading)

    computed_response = modelsteppredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.ϵ[i], atol=tol), eachindex(exact_response))
end
@test _modelsteppredict_nonsing_creep(tol)

function _modelsteppredict_nonsing_shifted_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = [i>=5.0 ? 1.0*exp(-(i-5.0)) : 0.0 for i in t]
    loading = ones(RheoFloat, length(t))

    modulus = quote α*exp.(-t/β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=1.0)

    data0 = RheoTimeData(t = t, ϵ = loading)

    computed_response = modelsteppredict(data0, model; step_on=5.0)

    all(i -> isapprox(exact_response[i], computed_response.σ[i], atol=tol), eachindex(exact_response))
end
@test _modelsteppredict_nonsing_shifted_relax(tol)

function _modelsteppredict_nonsing_shifted_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = [i>=5.0 ? 1.0*exp(-(i-5.0)) : 0.0 for i in t]
    loading = ones(RheoFloat, length(t))

    modulus = quote α*exp.(-t/β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=1.0)

    data0 = RheoTimeData(t = t, σ = loading)

    computed_response = modelsteppredict(data0, model; step_on=5.0)

    all(i -> isapprox(exact_response[i], computed_response.ϵ[i], atol=tol), eachindex(exact_response))
end
@test _modelsteppredict_nonsing_shifted_creep(tol)

function _modelsteppredict_sing_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = t.^(-0.5)
    loading = ones(RheoFloat, length(t))

    modulus = quote α*t.^(-β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=0.5)
    data0 = RheoTimeData(t = t, ϵ = loading)

    computed_response = modelsteppredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.σ[i], atol=tol), 1:length(t))
end
@test _modelsteppredict_sing_relax(tol)

function _modelsteppredict_sing_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = t.^(-0.5)
    loading = ones(RheoFloat, length(t))

    modulus = quote α*t.^(-β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=0.5)
    data0 = RheoTimeData(t = t, σ = loading)

    computed_response = modelsteppredict(data0, model)

    all(i -> isapprox(exact_response[i], computed_response.ϵ[i], atol=tol), 1:length(t))
end
@test _modelsteppredict_sing_creep(tol)

function _modelsteppredict_sing_shifted_relax(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = [i>=5.0 ? (i - 5.0).^(-0.5) : 0.0 for i in t]
    loading = ones(RheoFloat, length(t))

    modulus = quote α*t.^(-β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], G = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=0.5)
    data0 = RheoTimeData(t = t, ϵ = loading)

    computed_response = modelsteppredict(data0, model, step_on=5.0)

    stepon_el = RHEOS.closestindex(t, 5.0)

    all(i -> isapprox(exact_response[i], computed_response.σ[i], atol=tol), 1:length(t))
end
@test _modelsteppredict_sing_shifted_relax(tol)

function _modelsteppredict_sing_shifted_creep(tol)
    dt = 0.01
    t = Vector{RheoFloat}(0.0:dt:20.0)
    exact_response = [i>=5.0 ? (i - 5.0).^(-0.5) : 0.0 for i in t]
    loading = ones(RheoFloat, length(t))

    modulus = quote α*t.^(-β) end
    modelclass = RheoModelClass(name = "testmodel", p = [:α, :β], J = modulus, info="none")
    model = RheoModel(modelclass, α=1.0, β=0.5)
    data0 = RheoTimeData(t = t, σ = loading)

    computed_response = modelsteppredict(data0, model, step_on=5.0)

    stepon_el = RHEOS.closestindex(t, 5.0)

    all(i -> isapprox(exact_response[i], computed_response.ϵ[i], atol=tol), 1:length(t))
end
@test _modelsteppredict_sing_shifted_creep(tol)

function _obj_dynamic(tol)
    params = [1.0, 0.1]
    ω = Vector{RheoFloat}(0.0:0.1:100.0)

    dataGp = 2*ω
    dataGpp = 3*ω

    moduliGp(ω, params) = params[1]*2*ω .+ params[2]
    moduliGpp(ω, params) = params[1]*3*ω .+ params[2]

    cost = RHEOS.obj_dynamic(params, nothing, ω, dataGp, dataGpp, moduliGp, moduliGpp)

    costGp = sum((dataGp - moduliGp(ω, params)).^2)
    costGpp = sum((dataGpp - moduliGpp(ω, params)).^2)

    isapprox(cost, (costGp + costGpp))
end
@test _obj_dynamic(tol)

function _obj_dynamic_mean(tol)
    params = [1.0, 0.1]
    ω = Vector{RheoFloat}(0.0:0.1:100.0)

    dataGp = 2*ω
    dataGpp = 3*ω

    meanGp = sum(dataGp)/length(dataGp)
    meanGpp = sum(dataGpp)/length(dataGpp)

    moduliGp(ω, params) = params[1]*2*ω .+ params[2]
    moduliGpp(ω, params) = params[1]*3*ω .+ params[2]

    cost = RHEOS.obj_dynamic_mean(params, nothing, ω, dataGp, dataGpp, moduliGp, moduliGpp, meanGp, meanGpp)

    costGp = sum(((dataGp - moduliGp(ω, params))./meanGp).^2)
    costGpp = sum(((dataGpp - moduliGpp(ω, params))./meanGpp).^2)

    isapprox(cost, (costGp + costGpp))
end
@test _obj_dynamic_mean(tol)

function _obj_dynamic_log(tol)
    params = [1.0, 0.2]
    ω = Vector{RheoFloat}(0.0:0.1:100.0)

    dataGp = 2*ω .+ 0.1
    dataGpp = 3*ω .+ 0.1

    moduliGp(ω, params) = params[1]*2*ω .+ params[2]
    moduliGpp(ω, params) = params[1]*3*ω .+ params[2]

    cost = RHEOS.obj_dynamic_log(params, nothing, ω, dataGp, dataGpp, moduliGp, moduliGpp)

    costGp = sum(log.(dataGp./moduliGp(ω, params)).^2)
    costGpp = sum(log.(dataGpp./moduliGpp(ω, params)).^2)

    isapprox(cost, (costGp + costGpp))
end
@test _obj_dynamic_log(tol)

function _obj_dynamic_local(tol)
    params = [1.0, 0.2]
    ω = Vector{RheoFloat}(0.0:0.1:100.0)

    dataGp = 2*ω .+ 0.1
    dataGpp = 3*ω .+ 0.1

    moduliGp(ω, params) = params[1]*2*ω .+ params[2]
    moduliGpp(ω, params) = params[1]*3*ω .+ params[2]

    cost = RHEOS.obj_dynamic_local(params, nothing, ω, dataGp, dataGpp, moduliGp, moduliGpp)

    costGp = sum((1.0 .- moduliGp(ω, params)./dataGp).^2)
    costGpp = sum((1.0 .- moduliGpp(ω, params)./dataGpp).^2)

    isapprox(cost, (costGp + costGpp))
end
@test _obj_dynamic_local(tol)

function _obj_dynamic_manual(tol)
    params = [1.0, 0.2]
    ω = Vector{RheoFloat}(0.0:0.1:100.0)

    dataGp = 2*ω .+ 0.1
    dataGpp = 3*ω .+ 0.1

    moduliGp(ω, params) = params[1]*2*ω .+ params[2]
    moduliGpp(ω, params) = params[1]*3*ω .+ params[2]

    weights = [10.125, 21.2532602]

    cost = RHEOS.obj_dynamic_manual(params, nothing, ω, dataGp, dataGpp, moduliGp, moduliGpp, weights)

    costGp = sum((dataGp - moduliGp(ω, params)).^2)
    costGpp = sum((dataGpp - moduliGpp(ω, params)).^2)

    isapprox(cost, (costGp*weights[1] + costGpp*weights[2]))
end
@test _obj_dynamic_manual(tol)

function _dynamicmodelfit(tol)
    ω = Vector{RheoFloat}(0.0:0.01:20.0)
    actual_model = RheoModel(SLS_Zener, η=5.0, kᵦ=2.5, kᵧ=7.5)

    dataGp = storagemod(actual_model,ω)
    dataGpp = lossmod(actual_model,ω)
    data0 = RheoFreqData(ω = ω, Gp = dataGp, Gpp = dataGpp)

    init_params = (η=2.0, kᵦ=10.5, kᵧ=3.0)
    modelout = dynamicmodelfit(data0, SLS_Zener, p0=init_params, lo=(η=0.1, kᵦ=0.2, kᵧ=0.3), hi=(η=12.0, kᵦ=Inf, kᵧ=27.0), weights="none", rel_tol=1e-6)

    found_params = modelout.params
    actual_params = actual_model.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _dynamicmodelfit(tol)

function _dynamicmodelfit_noinit(tol)
    ω = Vector{RheoFloat}(0.0:0.01:20.0)
    actual_model = RheoModel(SLS_Zener, η=5.0, kᵦ=2.5, kᵧ=7.5)

    dataGp = storagemod(actual_model,ω)
    dataGpp = lossmod(actual_model,ω)
    data0 = RheoFreqData(ω = ω, Gp = dataGp, Gpp = dataGpp)

    modelout = dynamicmodelfit(data0, SLS_Zener, lo=(η=0.1, kᵦ=0.2, kᵧ=0.3), hi=(η=12.0, kᵦ=Inf, kᵧ=27.0), weights="none", rel_tol=1e-6)

    found_params = modelout.params
    actual_params = actual_model.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _dynamicmodelfit_noinit(tol)

function _dynamicmodelfit_nobounds(tol)
    ω = Vector{RheoFloat}(0.0:0.01:20.0)
    actual_model = RheoModel(SLS_Zener, η=5.0, kᵦ=2.5, kᵧ=7.5)

    dataGp = storagemod(actual_model,ω)
    dataGpp = lossmod(actual_model,ω)
    data0 = RheoFreqData(ω = ω, Gp = dataGp, Gpp = dataGpp)

    init_params = (η=2.0, kᵦ=10.5, kᵧ=3.0)
    modelout = dynamicmodelfit(data0, SLS_Zener, p0=init_params, weights="none", rel_tol=1e-6)

    found_params = modelout.params
    actual_params = actual_model.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _dynamicmodelfit_nobounds(tol)

function _dynamicmodelfit_mean(tol)
    ω = Vector{RheoFloat}(0.0:0.01:20.0)
    actual_model = RheoModel(SLS_Zener, η=5.0, kᵦ=2.5, kᵧ=7.5)

    dataGp = storagemod(actual_model,ω)
    dataGpp = lossmod(actual_model,ω)
    data0 = RheoFreqData(ω = ω, Gp = dataGp, Gpp = dataGpp)

    init_params = (η=2.0, kᵦ=10.5, kᵧ=3.0)
    modelout = dynamicmodelfit(data0, SLS_Zener, p0=init_params, lo=(η=0.1, kᵦ=0.2, kᵧ=0.3), hi=(η=12.0, kᵦ=Inf, kᵧ=27.0), weights="mean", rel_tol=1e-6)

    found_params = modelout.params
    actual_params = actual_model.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _dynamicmodelfit_mean(tol)

function _dynamicmodelfit_log(tol)
    ω = Vector{RheoFloat}(0.01:0.01:20.0)
    actual_model = RheoModel(SLS_Zener, η=5.0, kᵦ=2.5, kᵧ=7.5)

    dataGp = storagemod(actual_model,ω)
    dataGpp = lossmod(actual_model,ω)
    data0 = RheoFreqData(ω = ω, Gp = dataGp, Gpp = dataGpp)

    init_params = (η=2.0, kᵦ=10.5, kᵧ=3.0)
    modelout = dynamicmodelfit(data0, SLS_Zener, p0=init_params, lo=(η=0.1, kᵦ=0.2, kᵧ=0.3), hi=(η=12.0, kᵦ=Inf, kᵧ=27.0), weights="log", rel_tol=1e-6)

    found_params = modelout.params
    actual_params = actual_model.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _dynamicmodelfit_log(tol)

function _dynamicmodelfit_local(tol)
    ω = Vector{RheoFloat}(0.01:0.01:20.0)
    actual_model = RheoModel(SLS_Zener, η=5.0, kᵦ=2.5, kᵧ=7.5)

    dataGp = storagemod(actual_model,ω)
    dataGpp = lossmod(actual_model,ω)
    data0 = RheoFreqData(ω = ω, Gp = dataGp, Gpp = dataGpp)

    init_params = (η=2.0, kᵦ=10.5, kᵧ=3.0)
    modelout = dynamicmodelfit(data0, SLS_Zener, p0=init_params, lo=(η=0.1, kᵦ=0.2, kᵧ=0.3), hi=(η=12.0, kᵦ=Inf, kᵧ=27.0), weights="local", rel_tol=1e-6)

    found_params = modelout.params
    actual_params = actual_model.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _dynamicmodelfit_local(tol)

function _dynamicmodelfit_manualweights(tol)
    ω = Vector{RheoFloat}(0.01:0.01:20.0)
    actual_model = RheoModel(SLS_Zener, η=5.0, kᵦ=2.5, kᵧ=7.5)

    dataGp = storagemod(actual_model,ω)
    dataGpp = lossmod(actual_model,ω)
    data0 = RheoFreqData(ω = ω, Gp = dataGp, Gpp = dataGpp)

    init_params = (η=2.0, kᵦ=10.5, kᵧ=3.0)
    modelout = dynamicmodelfit(data0, SLS_Zener, p0=init_params, lo=(η=0.1, kᵦ=0.2, kᵧ=0.3), hi=(η=12.0, kᵦ=Inf, kᵧ=27.0), weights=[0.5, 1.5], rel_tol=1e-6)

    found_params = modelout.params
    actual_params = actual_model.params

    isapprox(collect(values(found_params)), collect(values(actual_params)), atol = tol)
end
@test _dynamicmodelfit_manualweights(tol)

function _dynamicmodelpredict(tol)
    ω = Vector{RheoFloat}(0.0:0.01:20.0)
    actual_model = RheoModel(SLS_Zener, η=5.0, kᵦ=2.5, kᵧ=7.5)

    dataGp = storagemod(actual_model,ω)
    dataGpp = lossmod(actual_model,ω)
    data0 = RheoFreqData(ω = ω, Gp = dataGp, Gpp = dataGpp)

    datapredicted = dynamicmodelpredict(data0, actual_model)

    test1 = isapprox(dataGp, datapredicted.Gp)
    test2 = isapprox(dataGpp, datapredicted.Gpp)
    test3 = isapprox(ω, datapredicted.ω)

    test1 && test2 && test3
end
@test _dynamicmodelpredict(tol)


text:
module TestKroneckerProducts
using Test
using LinearAlgebra
using KroneckerProducts
using KroneckerProducts: kronecker, ⊗, KroneckerProduct, issquare, allsquare, isequiv, ispower, order, samesize
using Random
using LazyInverses

element_types = (Float64, ComplexF64)
for elty in element_types
    @testset "eltype = $elty" begin
        # IDEA: loop over kronecker matrices
        A = randn(elty, 3, 3)
        B = randn(elty, 2, 4)
        C = randn(elty, 4, 4)
        D = randn(elty, 5, 5)
        AA = A'A + I # to test positive definite matrices
        CC = C'C + I
        DD = D'D + I
        tuples = [(A, A), (A, A, copy(A)), (A, A, A), # square kronecker power
                (A, B), (A, C, D), (A, B, C, D), # square kronecker products
                (A, B, C), (B, B), (B, B, B), # rectangular kronecker products
                (AA, AA), (AA, CC), (AA, CC, DD)] # positive definite products
        @testset "smart constructor" begin
            @test kronecker(A) === A
            a = randn(elty)
            @test kronecker(a) == a
            @test Matrix(KroneckerProduct((A,))) ≈ A
            pow = 3
            KP = kronecker_power(A, pow)
            @test ispower(KP)
            @test isequiv(KP)
            @test order(KP) == pow
            @test order(A) == 1
            KT = KroneckerProduct((A,), 3, elty)
            @test eltype(KT.temporaries[1]) == elty

            KT2 = kronecker(KT, KT)
            @test KT2 isa KroneckerProduct
            @test length(KT2.factors) == 2length(KT.factors) # testing hierarchy got flattened

            KTB = kronecker(KT, B)
            @test length(KTB.factors) == length(KT.factors) + 1
            @test KTB.factors[end] === B

            BKT = kronecker(B, KT)
            @test length(BKT.factors) == length(KT.factors) + 1
            @test BKT.factors[1] === B
        end
        for (i, tuple) in enumerate(tuples)
            @testset "kronecker product $i" begin
                K = kronecker(tuple)
                isposdef_K = isposdef(K) # compute once since we need it multiple times
                M = kron(tuple...)
                @testset "basic properties" begin
                    @test K isa KroneckerProduct
                    @test size(K) == size(M)
                    @test order(K) == length(tuple)
                    @test IndexStyle(typeof(K)) == IndexCartesian()
                    @test Matrix(K) ≈ M
                    @test kronecker(K) === K

                    @test issquare(K) == (size(K, 1) == size(K, 2))
                    @test samesize(K) == all(==(size(tuple[1])), tuple)
                    @test allsquare(K) == all(issquare, tuple)
                    @test ispower(K) == all(==(tuple[1]), tuple)
                    @test isequiv(K) == all(A->A===tuple[1], tuple)
                    @test isposdef_K == isposdef(M)
                    @test issymmetric(K) == issymmetric(M)

                    # testing adjoint, transpose, conj
                    @test Matrix(adjoint(K)) ≈ adjoint(M)
                    @test Matrix(transpose(K)) ≈ transpose(M)
                    @test Matrix(conj(K)) ≈ conj(M)
                end

                @testset "determinants" begin
                    if issquare(K)
                        @test det(K) ≈ det(M)
                        @test logabsdet(K)[1] ≈ logabsdet(M)[1]
                        @test logabsdet(K)[2] ≈ logabsdet(M)[2]
                        @test tr(K) ≈ tr(M)
                    end
                    if isposdef_K
                        @test logdet(K) ≈ logdet(M)
                    end
                end

                @testset "factorization" begin
                    # factorizing generic kronecker product
                    F = factorize(K)
                    @test F isa FactorizedKroneckerProduct
                    qr_K = qr(K)
                    @test all(F -> F isa LinearAlgebra.QRCompactWY, qr_K.factors)

                    pivoted_qr_K = qr(K, ColumnNorm())
                    @test all(F -> F isa LinearAlgebra.QRCompactWY, qr_K.factors)
                    if isposdef_K
                        cholesky_K = cholesky(K)
                        @test issuccess(cholesky_K)
                        @test all(F -> F isa Cholesky, cholesky_K.factors)
                        pivoted_cholesky_K = cholesky(K, Val(true))
                        @test issuccess(pivoted_cholesky_K)
                        @test all(F -> F isa CholeskyPivoted, pivoted_cholesky_K.factors)
                    end
                    if issquare(K)
                        inv_M = inv(M)
                        @test Matrix(inv(K)) ≈ inv_M
                        @test Matrix(inverse(K)) ≈ inv_M

                        x = randn(elty, size(F, 2))
                        y = zero(x)
                        z = M \ x
                        @test ldiv!(y, F, x) ≈ z
                        @test ldiv!(F, x) ≈ z
                    end
                    pinv_M = pinv(M)
                    @test Matrix(pinv(K)) ≈ pinv_M
                    @test Matrix(pseudoinverse(K)) ≈ pinv_M
                end

                @testset "indexing" begin
                    for i in 1:size(K, 1), j in 1:size(K, 2)
                        @test K[i, j] == M[i, j]
                    end
                end

                @testset "mul!" begin # IDEA: test for more complex type combinations
                    x = randn(elty, size(K, 2))
                    y = randn(elty, size(K, 1))
                    mul!(y, K, x)
                    @test y ≈ M*x
                    α, β = randn(2)
                    r = α * M * x + β * y
                    mul!(y, K, x, α, β)
                    @test y ≈ r

                    # scalar lmul!
                    KC = copy(K)
                    @test Matrix(KC) ≈ M
                    @test Matrix(lmul!(α, KC)) ≈ α * M
                    KC = copy(K)
                    @test Matrix(rmul!(KC, α)) ≈ α * M

                    # ldiv! - need K to be factorized
                    # F = factorize(K)
                    # @test ldiv!(y, F, x) ≈ M \ x
                    # @test ldiv!(F, x) ≈ M \ x
                end

                @testset "algebra" begin
                    # right multiplication
                    x = randn(elty, size(K, 2))
                    @test K * x ≈ M * x

                    # left multiplication
                    x = randn(elty, size(K, 1))
                    @test x'K ≈ x'M

                    # matmat
                    k = 3
                    X = randn(elty, size(K, 2), k)
                    @test K*X ≈ M*X
                    X = randn(elty, size(K, 1), k)
                    @test X'*K ≈ X'*M
                    X = randn(elty, k, size(K, 1))
                    @test X*K ≈ X*M
                    X = randn(elty, k, size(K, 2))
                    @test K*X' ≈ M*X'

                    # right solve
                    if allsquare(K)
                        b = randn(elty, size(K, 1))
                        x = K \ b
                        @test K*x ≈ b
                    end
                    if isposdef_K
                        Chol = ⊗(cholesky, K)
                        @test isposdef(Chol)
                        b = randn(elty, size(Chol, 1))
                        @test K * (Chol \ b) ≈ b
                    end
                    # mixed product property
                    @test Matrix(K'K) ≈ M'M
                    a = randn(elty)
                    @test Matrix(a*K) ≈ a*M
                    @test Matrix(K*a) ≈ M*a
                end
            end # testset
        end # loop over kronecker_products
    end # testset
end # loop over elty

end # TestKroneckerProducts

ext: jl


text:
using GameZero
using Test

@testset "get positions" begin
    r = Rect(100, 200, 10, 20)
    @test r.pos == (100, 200)
    @test r.x == 100
    @test r.y == 200
    @test r.topleft == (100, 200)
    @test r.top == 200
    @test r.left == 100
    @test r.topright == (110, 200)
    @test r.right == 110
    @test r.bottomleft == (100, 220)
    @test r.bottom == 220
    @test r.bottomright == (110, 220)
    @test r.center == (105, 210)
    @test r.centerx == 105
    @test r.centery == 210
    @test r.centerleft == (100, 210)
    @test r.centerright == (110, 210)
    @test r.topcenter == (105, 200)
    @test r.bottomcenter == (105, 220)
end

@testset "set positions" begin

@testset    begin
    r = Rect(300, 400, 10, 20)
    @test r.pos == (300, 400)
    r.pos = (100,200)
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.x = 100
    r.y = 200
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.topleft = (100,200)
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.top = 200
    r.left = 100
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.topright = (110, 200)
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.bottomleft = (100, 220)
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.bottomright = (110, 220)
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.bottom = 220
    r.right = 110
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.center = (105, 210)
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.centerx = 105
    r.centery = 210
    @test r.pos == (100, 200)
end

@testset begin
    r = Rect(300, 400, 10, 20)
    r.centerleft = (100, 210)
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.centerright = (110, 210)
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.topcenter = (105, 200)
    @test r.pos == (100, 200)
end
@testset begin
    r = Rect(300, 400, 10, 20)
    r.bottomcenter = (105, 220)
    @test r.pos == (100, 200)
end
end

include("examples.jl")


text:
radius2(x, y) = x^2 + y^2
function loop_radius2(n)
    s = 0
    for i = 1:n
        s += radius2(1, i)
    end
    s
end

tmppath = ""
if isdefined(Main, :Revise)
    global tmppath
    tmppath, io = mktemp()
    print(io, """
    function jikwfunc(x, y=0; z="hello")
        a = x + y
        b = z^a
        return length(b)
    end
    """)
    close(io)
    includet(tmppath)
end

using JuliaInterpreter, Test

function stacklength(frame)
    n = 1
    frame = frame.callee
    while frame !== nothing
        n += 1
        frame = frame.callee
    end
    return n
end

@testset "Breakpoints" begin
    breakpoint(radius2)
    frame = JuliaInterpreter.enter_call(loop_radius2, 2)
    bp = JuliaInterpreter.finish_and_return!(frame)
    @test isa(bp, JuliaInterpreter.BreakpointRef)
    @test stacklength(frame) == 2
    @test leaf(frame).framecode.scope == @which radius2(0, 0)
    bp = JuliaInterpreter.finish_stack!(frame)
    @test isa(bp, JuliaInterpreter.BreakpointRef)
    @test stacklength(frame) == 2
    @test JuliaInterpreter.finish_stack!(frame) == loop_radius2(2)

    # Conditional breakpoints
    function runsimple()
        frame = JuliaInterpreter.enter_call(loop_radius2, 2)
        bp = JuliaInterpreter.finish_and_return!(frame)
        @test isa(bp, JuliaInterpreter.BreakpointRef)
        @test stacklength(frame) == 2
        @test leaf(frame).framecode.scope == @which radius2(0, 0)
        @test JuliaInterpreter.finish_stack!(frame) == loop_radius2(2)
    end
    remove()
    breakpoint(radius2, :(y > x))
    runsimple()
    remove()
    @breakpoint radius2(0,0) y>x
    runsimple()
    # Demonstrate the problem that we have with scope
    local_identity(x) = identity(x)
    remove()
    @breakpoint radius2(0,0) y>local_identity(x)
    @test_broken @interpret loop_radius2(2)

    # Conditional breakpoints on local variables
    remove()
    halfthresh = loop_radius2(5)
    @breakpoint loop_radius2(10) 5 s>$halfthresh
    frame, bp = @interpret loop_radius2(10)
    @test isa(bp, JuliaInterpreter.BreakpointRef)
    lframe = leaf(frame)
    s_extractor = eval(JuliaInterpreter.prepare_slotfunction(lframe.framecode, :s))
    @test s_extractor(lframe) == loop_radius2(6)
    JuliaInterpreter.finish_stack!(frame)
    @test s_extractor(lframe) == loop_radius2(7)
    disable(bp)
    @test JuliaInterpreter.finish_stack!(frame) == loop_radius2(10)

    # Next line with breakpoints
    function outer(x)
        inner(x)
    end
    function inner(x)
        return 2
    end
    breakpoint(inner)
    frame = JuliaInterpreter.enter_call(outer, 0)
    bp = JuliaInterpreter.next_line!(frame)
    @test isa(bp, JuliaInterpreter.BreakpointRef)
    @test JuliaInterpreter.finish_stack!(frame) == 2

    # Breakpoints by file/line
    if isdefined(Main, :Revise)
        remove()
        method = which(JuliaInterpreter.locals, Tuple{Frame})
        breakpoint(String(method.file), method.line+1)
        frame = JuliaInterpreter.enter_call(loop_radius2, 2)
        ret = @interpret JuliaInterpreter.locals(frame)
        @test isa(ret, Tuple{Frame,JuliaInterpreter.BreakpointRef})
        # Test kwarg method
        remove()
        bp = breakpoint(tmppath, 3)
        frame, bp2 = @interpret jikwfunc(2)
        @test bp2 == bp
        var = JuliaInterpreter.locals(leaf(frame))
        @test !any(v->v.name == :b, var)
        @test filter(v->v.name == :a, var)[1].value == 2
    else
        try
            breakpoint(pathof(JuliaInterpreter.CodeTracking), 5)
        catch err
            @test isa(err, ErrorException)
            @test occursin("Revise", err.msg)
        end
    end

    # Direct return
    @breakpoint gcd(1,1) a==5
    @test @interpret(gcd(10,20)) == 10
    # FIXME: even though they pass, these tests break Test!
    # frame, bp = @interpret gcd(5, 20)
    # @test stacklength(frame) == 1
    # @test isa(bp, JuliaInterpreter.BreakpointRef)
    remove()

    # break on error
    try
        @test_throws ArgumentError("unsupported state :missing") break_on(:missing)
        break_on(:error)

        inner(x) = error("oops")
        outer() = inner(1)
        frame = JuliaInterpreter.enter_call(outer)
        bp = JuliaInterpreter.finish_and_return!(frame)
        @test bp.err == ErrorException("oops")
        @test stacklength(frame) >= 2
        @test frame.framecode.scope.name == :outer
        cframe = frame.callee
        @test cframe.framecode.scope.name == :inner

        # Don't break on caught exceptions
        function f_exc_outer()
            try
                f_exc_inner()
            catch err
                return err
            end
        end
        function f_exc_inner()
            error()
        end
        frame = JuliaInterpreter.enter_call(f_exc_outer);
        v = JuliaInterpreter.finish_and_return!(frame)
        @test v isa ErrorException
        @test stacklength(frame) == 1
    finally
        break_off(:error)
    end

    # Breakpoint display
    io = IOBuffer()
    frame = JuliaInterpreter.enter_call(loop_radius2, 2)
    bp = JuliaInterpreter.BreakpointRef(frame.framecode, 1)
    show(io, bp)
    @test String(take!(io)) == "breakpoint(loop_radius2(n) in $(@__MODULE__) at $(@__FILE__):3, line 3)"
    bp = JuliaInterpreter.BreakpointRef(frame.framecode, 0)  # fictive breakpoint
    show(io, bp)
    @test String(take!(io)) == "breakpoint(loop_radius2(n) in $(@__MODULE__) at $(@__FILE__):3, %0)"
    bp = JuliaInterpreter.BreakpointRef(frame.framecode, 1, ArgumentError("whoops"))
    show(io, bp)
    @test String(take!(io)) == "breakpoint(loop_radius2(n) in $(@__MODULE__) at $(@__FILE__):3, line 3, ArgumentError(\"whoops\"))"

    # In source breakpointing
    f_outer_bp(x) = g_inner_bp(x)
    function g_inner_bp(x)
        sin(x)
        @bp
        @bp
        @bp
        x = 3
        return 2
    end
    fr, bp = @interpret f_outer_bp(3)
    @test leaf(fr).framecode.scope.name == :g_inner_bp
    @test bp.stmtidx == 3
end

if tmppath != ""
    rm(tmppath)
end


text:
@test mysignif(SpecificG(0.001, 0.534_433_241E3), TestDigits) ≈ mysignif(SpecificG_Ph(0.001, 3000), TestDigits)
@test mysignif(SpecificG(0.001u"MPa", 0.534_433_241E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ph(0.001u"MPa", 3000u"kJ/kg"), TestDigits)

@test mysignif(SpecificG(3.0, 0.575_373_370E3), TestDigits) ≈ mysignif(SpecificG_Ph(3.0, 3000), TestDigits)
@test mysignif(SpecificG(3.0u"MPa", 0.575_373_370E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ph(3.0u"MPa", 3000u"kJ/kg"), TestDigits)

@test mysignif(SpecificG(3.0, 0.101_077_577E4), TestDigits) ≈ mysignif(SpecificG_Ph(3.0, 4000), TestDigits)
@test mysignif(SpecificG(3.0u"MPa", 0.101_077_577E4u"K"), TestDigits) ≈ mysignif(SpecificG_Ph(3.0u"MPa", 4000u"kJ/kg"), TestDigits)

@test mysignif(SpecificF(0.001, 0.534_433_241E3), TestDigits) ≈ mysignif(SpecificF_Ph(0.001, 3000), TestDigits)
@test mysignif(SpecificF(0.001u"MPa", 0.534_433_241E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ph(0.001u"MPa", 3000u"kJ/kg"), TestDigits)

@test mysignif(SpecificF(3.0, 0.575_373_370E3), TestDigits) ≈ mysignif(SpecificF_Ph(3.0, 3000), TestDigits)
@test mysignif(SpecificF(3.0u"MPa", 0.575_373_370E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ph(3.0u"MPa", 3000u"kJ/kg"), TestDigits)

@test mysignif(SpecificF(3.0, 0.101_077_577E4), TestDigits) ≈ mysignif(SpecificF_Ph(3.0, 4000), TestDigits)
@test mysignif(SpecificF(3.0u"MPa", 0.101_077_577E4u"K"), TestDigits) ≈ mysignif(SpecificF_Ph(3.0u"MPa", 4000u"kJ/kg"), TestDigits)

@test SpecificV(0.001, 0.534_433_241E3) ≈ SpecificV_Ph(0.001, 3000)
@test SpecificV(0.001u"MPa", 0.534_433_241E3u"K") ≈ SpecificV_Ph(0.001u"MPa", 3000u"kJ/kg")

@test SpecificV(3.0, 0.575_373_370E3) ≈ SpecificV_Ph(3.0, 3000)
@test SpecificV(3.0u"MPa", 0.575_373_370E3u"K") ≈ SpecificV_Ph(3.0u"MPa", 3000u"kJ/kg")

@test SpecificV(3.0, 0.101_077_577E4) ≈ SpecificV_Ph(3.0, 4000)
@test SpecificV(3.0u"MPa", 0.101_077_577E4u"K") ≈ SpecificV_Ph(3.0u"MPa", 4000u"kJ/kg")

@test SpecificU(0.001, 0.534_433_241E3) ≈ SpecificU_Ph(0.001, 3000)
@test SpecificU(0.001u"MPa", 0.534_433_241E3u"K") ≈ SpecificU_Ph(0.001u"MPa", 3000u"kJ/kg")

@test SpecificU(3.0, 0.575_373_370E3) ≈ SpecificU_Ph(3.0, 3000)
@test SpecificU(3.0u"MPa", 0.575_373_370E3u"K") ≈ SpecificU_Ph(3.0u"MPa", 3000u"kJ/kg")

@test SpecificU(3.0, 0.101_077_577E4) ≈ SpecificU_Ph(3.0, 4000)
@test SpecificU(3.0u"MPa", 0.101_077_577E4u"K") ≈ SpecificU_Ph(3.0u"MPa", 4000u"kJ/kg")

@test SpecificS(0.001, 0.534_433_241E3) ≈ SpecificS_Ph(0.001, 3000)
@test SpecificS(0.001u"MPa", 0.534_433_241E3u"K") ≈ SpecificS_Ph(0.001u"MPa", 3000u"kJ/kg")

@test SpecificS(3.0, 0.575_373_370E3) ≈ SpecificS_Ph(3.0, 3000)
@test SpecificS(3.0u"MPa", 0.575_373_370E3u"K") ≈ SpecificS_Ph(3.0u"MPa", 3000u"kJ/kg")

@test SpecificS(3.0, 0.101_077_577E4) ≈ SpecificS_Ph(3.0, 4000)
@test SpecificS(3.0u"MPa", 0.101_077_577E4u"K") ≈ SpecificS_Ph(3.0u"MPa", 4000u"kJ/kg")

@test SpecificCP(0.001, 0.534_433_241E3) ≈ SpecificCP_Ph(0.001, 3000)
@test SpecificCP(0.001u"MPa", 0.534_433_241E3u"K") ≈ SpecificCP_Ph(0.001u"MPa", 3000u"kJ/kg")

@test SpecificCP(3.0, 0.575_373_370E3) ≈ SpecificCP_Ph(3.0, 3000)
@test SpecificCP(3.0u"MPa", 0.575_373_370E3u"K") ≈ SpecificCP_Ph(3.0u"MPa", 3000u"kJ/kg")

@test SpecificCP(3.0, 0.101_077_577E4) ≈ SpecificCP_Ph(3.0, 4000)
@test SpecificCP(3.0u"MPa", 0.101_077_577E4u"K") ≈ SpecificCP_Ph(3.0u"MPa", 4000u"kJ/kg")

@test SpecificCV(0.001, 0.534_433_241E3) ≈ SpecificCV_Ph(0.001, 3000)
@test SpecificCV(0.001u"MPa", 0.534_433_241E3u"K") ≈ SpecificCV_Ph(0.001u"MPa", 3000u"kJ/kg")

@test SpecificCV(3.0, 0.575_373_370E3) ≈ SpecificCV_Ph(3.0, 3000)
@test SpecificCV(3.0u"MPa", 0.575_373_370E3u"K") ≈ SpecificCV_Ph(3.0u"MPa", 3000u"kJ/kg")

@test SpecificCV(3.0, 0.101_077_577E4) ≈ SpecificCV_Ph(3.0, 4000)
@test SpecificCV(3.0u"MPa", 0.101_077_577E4u"K") ≈ SpecificCV_Ph(3.0u"MPa", 4000u"kJ/kg")

@test SpeedOfSound(0.001, 0.534_433_241E3) ≈ SpeedOfSound_Ph(0.001, 3000)
@test SpeedOfSound(0.001u"MPa", 0.534_433_241E3u"K") ≈ SpeedOfSound_Ph(0.001u"MPa", 3000u"kJ/kg")

@test SpeedOfSound(3.0, 0.575_373_370E3) ≈ SpeedOfSound_Ph(3.0, 3000)
@test SpeedOfSound(3.0u"MPa", 0.575_373_370E3u"K") ≈ SpeedOfSound_Ph(3.0u"MPa", 3000u"kJ/kg")

@test SpeedOfSound(3.0, 0.101_077_577E4) ≈ SpeedOfSound_Ph(3.0, 4000)
@test SpeedOfSound(3.0u"MPa", 0.101_077_577E4u"K") ≈ SpeedOfSound_Ph(3.0u"MPa", 4000u"kJ/kg")


text:
function ComplexComparator()
    return ComplexComparator(())
end

function compare(obj::ComplexComparator, arg0::JComplex, arg1::JComplex)
    return jcall(obj, "compare", jint, (JComplex, JComplex), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function reversed(obj::Comparator)
    return jcall(obj, "reversed", Comparator, ())
end

function then_comparing(obj::Comparator, arg0::Comparator)
    return jcall(obj, "thenComparing", Comparator, (Comparator,), arg0)
end

function then_comparing(obj::Comparator, arg0::Function)
    return jcall(obj, "thenComparing", Comparator, (Function,), arg0)
end

function then_comparing(obj::Comparator, arg0::Function, arg1::Comparator)
    return jcall(obj, "thenComparing", Comparator, (Function, Comparator), arg0, arg1)
end

function then_comparing_double(obj::Comparator, arg0::ToDoubleFunction)
    return jcall(obj, "thenComparingDouble", Comparator, (ToDoubleFunction,), arg0)
end

function then_comparing_int(obj::Comparator, arg0::ToIntFunction)
    return jcall(obj, "thenComparingInt", Comparator, (ToIntFunction,), arg0)
end

function then_comparing_long(obj::Comparator, arg0::ToLongFunction)
    return jcall(obj, "thenComparingLong", Comparator, (ToLongFunction,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
"""
    module CplTestingBL
 Defines an equation set (balance law) for testing coupling

 Defines kernels to evaluates RHS of

    ```
    ∂θ
    --  = - ∇ • ( u θ + κ ∇ θ ) = - ∇ • F
    ∂t
    ```

    Where

     - `θ` is the tracer (e.g. potential temperature)
     - `u` is the initial advection velocity (constant)
     - `κ` is the diffusivity tensor

 subject to specified prescribed gradient or flux bc's

 - [`CplMainBL`](@ref)
     balance law struct created by this module

"""
#module CplMainBL

export CplMainBL
export PenaltyNumFluxDiffusive
export ExteriorBoundary
export CoupledPrimaryBoundary, CoupledSecondaryBoundary

using ClimateMachine.BalanceLaws: Auxiliary, BalanceLaw, Gradient, GradientFlux, Prognostic

import ClimateMachine.BalanceLaws:
    boundary_conditions,
    boundary_state!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    init_state_prognostic!,
    nodal_init_state_auxiliary!,
    source!,
    vars_state,
    wavespeed

using ClimateMachine.Mesh.Geometry: LocalGeometry
using ClimateMachine.MPIStateArrays

using ClimateMachine.DGMethods.NumericalFluxes:
    CentralNumericalFluxGradient,
    CentralNumericalFluxSecondOrder,
    NumericalFluxFirstOrder,
    NumericalFluxSecondOrder,
    RusanovNumericalFlux

import ClimateMachine.DGMethods.NumericalFluxes:
    numerical_boundary_flux_second_order!,
    numerical_flux_second_order!,
    numerical_boundary_flux_first_order!,
    numerical_flux_first_order!

using ClimateMachine.VariableTemplates

using LinearAlgebra
using StaticArrays

"""
    CplMainBL <: BalanceLaw
    - used in both the atmos and ocean component models
"""
struct CplMainBL{BLP, BCS, PS} <: BalanceLaw
    bl_prop::BLP  # initial condition, ...
    boundaryconditions::BCS
    param_set::PS
end

l_type = CplMainBL

"""
Declare prognostic state variables
    - θ
    - F_accum: "Shadow" variable used to capture boundary fluxes that we want to accumulate over
    a single timestep and export to coupling as time integrals. We use a shadow variable
    because we want to integrate over whatever timestepper is being used.
    Eventually we should have the ability to potentially use a 2d field here.
    The shadow variable needs to be zeroed at the start of each
    coupling cycle for a component.
    - u: advective velocity - currently constant, but kept in prognostic variables to enable
    easier use of overintegration
"""
function vars_state(bl::l_type, ::Prognostic, FT)
    @vars begin
        θ::FT
        F_accum::FT
        u::SVector{3, FT}
    end
end

"""
Declare Aaxiliary state variables
  - for array index and real world coordinates and for
  θ value at reference time used to compute κ.
"""
function vars_state(bl::l_type, ::Auxiliary, FT)
    @vars begin
        npt::Int    # no. nodes
        elnum::Int  # no. elems

        xc::FT      # Cartesian x
        yc::FT      # Cartesian y
        zc::FT      # Cartesian z

        θⁱⁿⁱᵗ::FT   # unused in default setup
        θ_secondary::FT  # stores opposite face for primary (atmospheric import)
        F_prescribed::FT # stores prescribed flux for secondary (ocean import)
    end
end

"""
  Gradient computation stage input (and output) variable symbols
"""
function vars_state(bl::l_type, ::Gradient, FT)
    @vars begin
        ∇θ::FT
        ∇θⁱⁿⁱᵗ::FT # unused in default setup
    end
end

function vars_state(bl::l_type, ::GradientFlux, FT)
    @vars begin
        κ∇θ::SVector{3, FT}
    end
end

"""
  Initialize prognostic state variables
"""
function init_state_prognostic!(bl::l_type, Q::Vars, A::Vars, geom::LocalGeometry, FT)
    npt = getproperty(geom, :n)
    elnum = getproperty(geom, :e)
    x = A.xc
    y = A.yc
    z = A.zc

    Q.θ = bl.bl_prop.init_theta(npt, elnum, x, y, z)
    Q.F_accum = 0
    Q.u = bl.bl_prop.init_u(npt, elnum, x, y, z)
    nothing
end

"""
  Initialize auxiliary state variables
"""
function nodal_init_state_auxiliary!(bl::l_type, A::Vars, tmp::Vars, geom::LocalGeometry, _...)
    npt = getproperty(geom, :n)
    elnum = getproperty(geom, :e)
    x = geom.coord[1]
    y = geom.coord[2]
    z = geom.coord[3]

    A.npt, A.elnum, A.xc, A.yc, A.zc = bl.bl_prop.init_aux_geom(npt, elnum, x, y, z)

    A.θⁱⁿⁱᵗ = 0
    A.θ_secondary = 0
    A.F_prescribed = 0
    nothing
end

"""
Compute Kernels:
"""
#====

Atmos

----

Land

====#

"""
  Set source terms
  - for prognostic state external sources
  - for recording boundary flux terms into shadow variables for export to coupler
"""
function source!(bl::l_type, S::Vars, Q::Vars, G::Vars, A::Vars, _...)
    #S.θ=bl.bl_prop.source_theta(Q.θ,A.npt,A.elnum,A.xc,A.yc,A.zc,A.θ_secondary)
    # Record boundary condition fluxes as needed by adding to shadow
    # prognostic variable

    S.F_accum = (Q.θ - A.θ_secondary) * bl.bl_prop.coupling_lambda()
    nothing
end

"""
  Set values to have gradients computed.
"""
function compute_gradient_argument!(bl::l_type, G::Vars, Q::Vars, A::Vars, t)
    G.∇θ = Q.θ
    G.∇θⁱⁿⁱᵗ = A.θⁱⁿⁱᵗ
    nothing
end

"""
  Compute diffusivity tensor times computed gradient to give net gradient flux.
"""
function compute_gradient_flux!(bl::l_type, GF::Vars, G::Grad, Q::Vars, A::Vars, t)
    # "Non-linear" form (for time stepped)
    ### κ¹,κ²,κ³=bl.bl_prop.calc_kappa_diff(G.∇θ,A.npt,A.elnum,A.xc,A.yc,A.zc)
    # "Linear" form (for implicit)
    κ¹, κ², κ³ = bl.bl_prop.calc_kappa_diff(G.∇θⁱⁿⁱᵗ, A.npt, A.elnum, A.xc, A.yc, A.zc)
    # Maybe I should pass both G.∇θ and G.∇θⁱⁿⁱᵗ?

    # Messy diffusion tensor rotation - should find code for doing this succinctly!
    # we want κ¹, κ² to act in lat and lon directions respectively i.e. tangential
    # to shell, we want κ³ to act normal to shell i.e. z or r direction.
    ## don't use regular norm(), its not great on performance
    mynorm(x::Float64, y::Float64, z::Float64) = (x^2 + y^2 + z^2)^0.5
    ## excessively type stable forms!
    r̂ⁿᵒʳᵐ(x::Float64, y::Float64, z::Float64) = mynorm(x, y, z) ≈ 0 ? 1 : mynorm(x, y, z)^(-1)
    ϕ̂ⁿᵒʳᵐ(x::Float64, y::Float64, z::Float64) =
        mynorm(x, y, Float64(0)) ≈ 0 ? 1 : (mynorm(x, y, z) * mynorm(x, y, Float64(0)))^(-1)
    λ̂ⁿᵒʳᵐ(x::Float64, y::Float64, z::Float64) = mynorm(x, y, Float64(0)) ≈ 0 ? 1 : mynorm(x, y, Float64(0))^(-1)
    r̂(x::Float64, y::Float64, z::Float64) = r̂ⁿᵒʳᵐ(x, y, z) * @SVector([x, y, z])
    ϕ̂(x::Float64, y::Float64, z::Float64) = ϕ̂ⁿᵒʳᵐ(x, y, z) * @SVector [x * z, y * z, -(x^2 + y^2)]
    λ̂(x::Float64, y::Float64, z::Float64) = λ̂ⁿᵒʳᵐ(x, y, z) * @SVector [-y, x, 0]

    # Sphere surface bits
    x = A.xc
    y = A.yc
    z = A.zc
    Gλ = G.∇θ' * λ̂(x, y, z)
    Fλ = κ¹ * Gλ * λ̂(x, y, z)
    Gϕ = G.∇θ' * ϕ̂(x, y, z)
    Fϕ = κ² * Gϕ * ϕ̂(x, y, z)
    Gr = G.∇θ' * r̂(x, y, z)
    Fr = κ³ * Gr * r̂(x, y, z)
    GF.κ∇θ = Diagonal(@SVector([1, 1, 1])) * (Fλ + Fϕ + Fr)

    # Original cartesian form
    #GF.κ∇θ = Diagonal(@SVector([κ¹, κ², κ³])) * G.∇θ
    #GF.κ∇θ = Diagonal(@SVector([κ¹, κ¹, κ¹])) * G.∇θ
    nothing
end

"""
  Pass flux components for second order term into update kernel.
"""
function flux_second_order!(bl::l_type, F::Grad, Q::Vars, GF::Vars, H::Vars, A::Vars, t)
    F.θ += GF.κ∇θ
    nothing
end

# Boundary conditions
"""
  Define boundary condition flags/types to iterate over, for now keep it simple.
- 3 AbstractCouplerBoundary types:
    1. ExteriorBoundary
    2. CoupledPrimaryBoundary
    3. CoupledSecondaryBoundary
- PenaltyNumFluxDiffusive: additional NumericalFluxSecondOrder to account for an additional penalty term
"""
struct PenaltyNumFluxDiffusive <: NumericalFluxSecondOrder end

abstract type AbstractCouplerBoundary end

# ## 1. ExteriorBoundary
# flux is 0 across the boundary
struct ExteriorBoundary <: AbstractCouplerBoundary end

"""
  Zero normal gradient boundary condition.
"""

function boundary_state!(
    nF::Union{CentralNumericalFluxGradient},
    bc::ExteriorBoundary,
    bl::l_type,
    Q⁺::Vars,
    A⁺::Vars,
    n,
    Q⁻::Vars,
    A⁻::Vars,
    t,
    _...,
)
    Q⁺.θ = Q⁻.θ
    nothing
end
function numerical_boundary_flux_second_order!(
    numerical_flux::Union{NumericalFluxSecondOrder},
    bctype::ExteriorBoundary,
    balance_law::l_type,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, A, HD}

    fluxᵀn.θ += 0
end

# ## 2. CoupledPrimaryBoundary
# # compute flux based on opposite face
# # also need to accumulate net flux across boundary
struct CoupledPrimaryBoundary <: AbstractCouplerBoundary end

function boundary_state!(
    nF::Union{CentralNumericalFluxGradient},
    bc::CoupledPrimaryBoundary,
    bl::l_type,
    Q⁺::Vars,
    A⁺::Vars,
    n,
    Q⁻::Vars,
    A⁻::Vars,
    t,
    _...,
)

    Q⁺.θ = Q⁻.θ # Q⁺.θ=A⁺.θ_secondary

    nothing
end

function numerical_boundary_flux_second_order!(
    numerical_flux::Union{NumericalFluxSecondOrder},
    bctype::CoupledPrimaryBoundary,
    balance_law::l_type,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, A, HD}
    fluxᵀn.θ = (state_prognostic⁻.θ - state_auxiliary⁺.θ_secondary) * balance_law.bl_prop.coupling_lambda() # W/m^2

end

# θ: J / m^3
# ∇θ: J / m^4
# κ: (m^2/s)
# κ∇θ: J/(m^2 s) = W/m^2 = J/m^3 * m/s

#  - clean up and write primer
#  - imports and exports
#  - add boundary auxiliary interface
#  - add boundary flux accumulation interface
#  - vector quantity
#  - run on sphere to check normal terms are all correct
#  - single stack integration
#  - Held-Suarez with wind stress and heat

## 3. CoupledSecondaryBoundary
# use prescribed flux computed in primary
struct CoupledSecondaryBoundary <: AbstractCouplerBoundary end
function boundary_state!(
    nF::Union{CentralNumericalFluxGradient},
    bc::CoupledSecondaryBoundary,
    bl::l_type,
    Q⁺::Vars,
    A⁺::Vars,
    n,
    Q⁻::Vars,
    A⁻::Vars,
    t,
    _...,
)
    Q⁺.θ = Q⁻.θ
    nothing
end
function numerical_boundary_flux_second_order!(
    numerical_flux::Union{NumericalFluxSecondOrder},
    bctype::CoupledSecondaryBoundary,
    balance_law::l_type,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, A, HD}
    fluxᵀn.θ = -state_auxiliary⁺.F_prescribed
end

function wavespeed(bl::l_type, _...)
    # Used in Rusanov term.
    # Only active if there is a flux first order term?
    bl.bl_prop.get_wavespeed()
end

"""
  Penalty flux formulation of second order numerical flux. This formulation
  computes the CentralNumericalFluxSecondOrder term first (which is just the average
  of the + and - fluxes and an edge), and then adds a "penalty" flux that relaxes
  the edge state + and - toward each other.
"""
function numerical_flux_second_order!(
    ::PenaltyNumFluxDiffusive,
    bl::l_type,
    fluxᵀn::Vars{S},
    n::SVector,
    state⁻::Vars{S},
    diff⁻::Vars{D},
    hyperdiff⁻::Vars{HD},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    diff⁺::Vars{D},
    hyperdiff⁺::Vars{HD},
    aux⁺::Vars{A},
    t,
) where {S, HD, D, A}

    numerical_flux_second_order!(
        CentralNumericalFluxSecondOrder(),
        bl,
        fluxᵀn,
        n,
        state⁻,
        diff⁻,
        hyperdiff⁻,
        aux⁻,
        state⁺,
        diff⁺,
        hyperdiff⁺,
        aux⁺,
        t,
    )

    Fᵀn = parent(fluxᵀn)
    FT = eltype(Fᵀn)
    tau = bl.bl_prop.get_penalty_tau()
    Fᵀn .+= tau * (parent(state⁻) - parent(state⁺))
end

"""
  First order flux for advection equation

  NumericalFluxFirstOrder is an abstract type that currently generalizes
  RusanovNumericalFlux, CentralNumericalFluxFirstOrder, RoeNumericalFlux, HLLCNumericalFlux.
"""
function flux_first_order!(bl::l_type, F::Grad, Q::Vars, A::Vars, t::Real, directions)
    F.θ += Q.u * Q.θ
    nothing
end
function boundary_state!(
    nf::Union{NumericalFluxFirstOrder},
    bc::Union{CoupledPrimaryBoundary, CoupledSecondaryBoundary, ExteriorBoundary},
    atmos::l_type,
    Q⁺::Vars,
    A⁺::Vars,
    n,
    Q⁻::Vars,
    A⁻::Vars,
    t,
    state_int⁻,
    aux_int⁻,
)
    Q⁺.θ = Q⁻.θ
    nothing
end

function boundary_conditions(bl::l_type, _...)
    bl.boundaryconditions
end

# #this one is for the general adv only test
# function boundary_state!(
#     nF,
#     bc,
#     bl::l_type,
#     Q⁺,
#     state_gradient_flux⁺,
#     state_hyperdiffusive⁺,
#     state_auxiliary⁺,
#     normal_vector,
#     Q⁻,
#     state_gradient_flux⁻,
#     state_hyperdiffusive⁻,
#     state_auxiliary⁻,
#     t,
#     state1⁻,
#     diff1⁻,
#     aux1⁻,
# )
#     Q⁺.θ = Q⁻.θ # Q⁺.θ=A⁺.θ_secondary

#     nothing
# end

"""
  Set a default set of properties and their default values
  - init_aux_geom   :: function to initialize geometric terms stored in aux.
  - init_theta      :: function to set initial θ values.
  - source_theta    :: function to add a source term to θ.
  - calc_kappa_diff :: function to set diffusion coeffiecient(s).
  - get_wavespeed   :: function to return a wavespeed for Rusanov computations (there aren't any in this model)
  - get_penalty_tau :: function to set timescale on which to bring state+ and state- together
  - theta_shadow_boundary_flux :: function to set boundary flux into shadow variable for passing to coupler
"""
function prop_defaults()
    bl_prop = NamedTuple()

    function init_aux_geom(npt, elnum, x, y, z)
        return npt, elnum, x, y, z
    end
    # init_aux_geom(_...)=(return 0., 0., 0., 0., 0.)
    bl_prop = (bl_prop..., init_aux_geom = init_aux_geom)

    init_theta(_...) = (return 0.0)
    bl_prop = (bl_prop..., init_theta = init_theta)

    source_theta(_...) = (return 0.0)
    bl_prop = (bl_prop..., source_theta = source_theta)

    calc_kappa_diff(_...) = (return 0.0, 0.0, 0.0)
    bl_prop = (bl_prop..., calc_kappa_diff = calc_kappa_diff)

    get_wavespeed(_...) = (return 0.0)
    bl_prop = (bl_prop..., get_wavespeed = get_wavespeed)

    get_penalty_tau(_...) = (return 1.0)
    bl_prop = (bl_prop..., get_penalty_tau = get_penalty_tau)

    theta_shadow_boundary_flux(_...) = (return 0.0)
    bl_prop = (bl_prop..., theta_shadow_boundary_flux = theta_shadow_boundary_flux)

    coupling_lambda(_...) = (return 0.0)
    bl_prop = (bl_prop..., coupling_lambda = coupling_lambda)

    init_u(_...) = (return 0.0)
    bl_prop = (bl_prop..., init_u = init_u)

    bl_prop = (bl_prop..., LAW = CplMainBL)
end

"""
  Helper functions to communicate between components before and after timestepping
  - preatmos(csolver)
  - postatmos(csolver)
  - preocean(csolver)
  - postocean(csolver)
"""

function preatmos(csolver)
    mA = csolver.component_list.atmosphere.component_model
    mO = csolver.component_list.ocean.component_model
    # Set boundary SST used in atmos to SST of ocean surface at start of coupling cycle.
    mA.discretization.state_auxiliary.θ_secondary[mA.boundary] .=
        coupler_get(csolver.coupler, :Ocean_SST, mA.grid, DateTime(0), u"°C")
    # Set atmos boundary flux accumulator to 0.
    mA.state.F_accum .= 0

    @info(
        "preatmos",
        time = csolver.t,
        total_θ_atmos = weightedsum(mA.state, 1),
        total_θ_ocean = weightedsum(mO.state, 1),
        total_θ = weightedsum(mA.state, 1) + weightedsum(mO.state, 1),
        atmos_θ_surface_max = maximum(mA.state.θ[mA.boundary]),
        ocean_θ_surface_max = maximum(mO.state.θ[mO.boundary]),
    )
end

function postatmos(csolver)
    mA = csolver.component_list.atmosphere.component_model
    mO = csolver.component_list.ocean.component_model
    # Pass atmos exports to "coupler" namespace
    # 1. Save mean θ flux at the Atmos boundary during the coupling period
    coupler_put!(
        csolver.coupler,
        :Atmos_MeanAirSeaθFlux,
        mA.state.F_accum[mA.boundary] ./ csolver.dt,
        mA.grid,
        DateTime(0),
        u"°C",
    )

    @info(
        "postatmos",
        time = time = csolver.t + csolver.dt,
        total_θ_atmos = weightedsum(mA.state, 1),
        total_θ_ocean = weightedsum(mO.state, 1),
        total_F_accum = mean(mA.state.F_accum[mA.boundary]) * 1e6 * 1e6,
        total_θ = weightedsum(mA.state, 1) + weightedsum(mO.state, 1) + mean(mA.state.F_accum[mA.boundary]) * 1e6 * 1e6,
        F_accum_max = maximum(mA.state.F_accum[mA.boundary]),
        F_avg_max = maximum(mA.state.F_accum[mA.boundary] ./ csolver.dt),
        atmos_θ_surface_max = maximum(mA.state.θ[mA.boundary]),
        ocean_θ_surface_max = maximum(mO.state.θ[mO.boundary]),
    )
end

function preocean(csolver)
    mA = csolver.component_list.atmosphere.component_model
    mO = csolver.component_list.ocean.component_model
    # Set mean air-sea theta flux
    mO.discretization.state_auxiliary.F_prescribed[mO.boundary] .=
        coupler_get(csolver.coupler, :Atmos_MeanAirSeaθFlux, mO.grid, DateTime(0), u"°C")
    # Set ocean boundary flux accumulator to 0. (this isn't used)
    mO.state.F_accum .= 0

    @info(
        "preocean",
        time = csolver.t,
        F_prescribed_max = maximum(mO.discretization.state_auxiliary.F_prescribed[mO.boundary]),
        F_prescribed_min = maximum(mO.discretization.state_auxiliary.F_prescribed[mO.boundary]),
        ocean_θ_surface_max = maximum(mO.state.θ[mO.boundary]),
        ocean_θ_surface_min = maximum(mO.state.θ[mO.boundary]),
    )
end

function postocean(csolver)
    mA = csolver.component_list.atmosphere.component_model
    mO = csolver.component_list.ocean.component_model
    @info(
        "postocean",
        time = csolver.t + csolver.dt,
        ocean_θ_surface_max = maximum(mO.state.θ[mO.boundary]),
        ocean_θ_surface_min = maximum(mO.state.θ[mO.boundary]),
    )

    # Pass ocean exports to "coupler" namespace
    #  1. Ocean SST (value of θ at z=0)
    coupler_put!(csolver.coupler, :Ocean_SST, mO.state.θ[mO.boundary], mO.grid, DateTime(0), u"°C")
end

#end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
@testset "Correctness of conversion operations" begin

    @testset "Vector to permutation" begin
        # TODO: include cases when it fails
    end

    @testset "Permutation to matrix" begin
        # TODO
    end

    @testset "Matrix to permutation" begin
        # TODO: include validation that matrix represents a permutation
    end

end


text:
using DifferentialDynamicProgramming
using Base.Test

info("Compile time is high for this package, this is expected and is not an error.")
# write your own tests here
demoQP()
demo_linear()


text:
using Gadfly
using Distributions
using Colors
using DataFrames
using Query

Gadfly.push_theme(Theme(major_label_font="Droid Sans",
                        minor_label_font="Droid Sans",
                        major_label_font_size=18pt, minor_label_font_size=16pt,
                        line_width=0.8mm, key_label_font="Droid Sans",
                        lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2),
                        key_label_font_size=14pt,
                        default_color=colorant"#000000"))

colors = [colorant"#e41a1c", colorant"#377eb8", colorant"#4daf4a",
          colorant"#984ea3", colorant"#ff7f00", colorant"#ffff33"]

names = [:type, :problem, :id, :seed, :gen, :eval, :epoch, :total_epochs, :loss, :acc, :fit]

function get_res(log::String)
    readtable(log, header=false, separator=',', names=names)
end

function reducef(df, xmax)
    r = 1:length(df[:eval])
    if df[:eval][end] != xmax
        r = 1:(length(df[:eval])+1)
    end
    map(i->mapf(i, df, xmax), r)
end

function mapf(i::Int64, df, xmax::Int64)
    if i > length(df[:eval])
        if df[:eval][end] <= xmax
            return df[:fit][end] * ones(xmax - df[:eval][end])
        else
            return []
        end
    end
    lower = 0
    if i >= 2
        lower = df[:eval][i-1]
        return df[:fit][i] * ones(df[:eval][i] - lower)
    else
        return df[:fit][i] * ones(df[:eval][i])
    end
    return []
end

function get_stats(res::DataFrame; xmax::Int64 = maximum(res[:eval]))
    limited = res
    if xmax < maximum(res[:eval])
        limited = @from i in res begin
            @where i.eval <= xmax
            @select i
            @collect DataFrame
        end
    end
    filled = by(limited, [:id, :seed],
                df->reduce(vcat, reducef(df, xmax)))
    filled[:xs] = repeat(1:xmax, outer=Int64(size(filled,1)/xmax))
    stats = by(filled, [:xs, :id],
               df->DataFrame(stds=std(df[:x1]), means=mean(df[:x1]), mins=minimum(df[:x1]),
                             maxs=maximum(df[:x1])))
    stats[:stds][isnan.(stats[:stds])] = 0;
    stats[:lower] = stats[:means]-0.5*stats[:stds]
    stats[:upper] = stats[:means]+0.5*stats[:stds]
    stats
end

function plot_evo(stats::DataFrame, title::String, filename::String="training";
                  xmin=minimum(stats[:xs]),
                  xmax=maximum(stats[:xs]),
                  ymin=minimum(stats[:lower]),
                  ymax=maximum(stats[:upper]),
                  ylabel="Accuracy",
                  key_position=:right)
    plt = plot(stats, x="xs", y="means", ymin="lower", ymax="upper", color="id",
               Geom.line, Geom.ribbon,
               Scale.color_discrete_manual(colors...),
               Guide.title(title),
               Guide.xlabel("Evaluations"),
               Guide.ylabel(ylabel),
               Coord.cartesian(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),
               style(key_position=key_position))
    draw(PDF(string(filename, ".pdf"), 8inch, 6inch), plt)
    plt
end


text:
import Core: svec

Base.LambdaInfo() = ccall(:jl_new_lambda_info_uninit, Ref{LambdaInfo}, (Ptr{Void},), C_NULL)

function copyil(l::LambdaInfo)
  l′ = LambdaInfo()
  for f in fieldnames(l)
    isdefined(l, f) || continue
    val = getfield(l, f)
    isa(val, Array) && (val = copy(val))
    setfield!(l′, f, val)
  end
  return l′
end

function uncompress!(l::LambdaInfo)
  l.code = Base.uncompressed_ast(l)
  return l
end

typealias TupleType{T<:Tuple} Type{T}

jl_method_def(argdata::SimpleVector, f::LambdaInfo, isstaged::Bool = false) =
  ccall(:jl_method_def, Void, (Any, Any, Any), argdata, f, isstaged)

jl_method_def(types::TupleType, tvars::SimpleVector, f::LambdaInfo, isstaged::Bool = false) =
  jl_method_def(svec(types, tvars), f, isstaged)

tovec(sv::SimpleVector) = sv
tovec(x) = svec(x)

jl_method_def(f::LambdaInfo, isstaged = false) =
  jl_method_def(f.specTypes, tovec(f.def.tvars), f, isstaged)

function switchfunc!(l::LambdaInfo, f)
  l.specTypes = Tuple{typeof(f), l.specTypes.parameters[2:end]...}
  return l
end

function bump_slots!(l::LambdaInfo, start::Integer)
  map!(uncompress!(l).code) do ex
    MacroTools.prewalk(ex) do x
      isa(x, SlotNumber) && x.id ≥ start ? SlotNumber(x.id+1) :
      isa(x, NewvarNode) && x.slot.id ≥ start ? NewvarNode(SlotNumber(x.slot.id + 1)) :
      x
    end
  end
  return l
end

function insertarg!(l::LambdaInfo, n::Integer, v::Tuple{Symbol,Type})
  name, T = v
  insert!(l.slotnames, n, name)
  insert!(l.slotflags, n, 0)
  l.specTypes = Tuple{insert!(collect(l.specTypes.parameters), n, T)...}
  l.nargs += 1
  bump_slots!(l, n)
  return l
end


text:
# See https://github.com/TuringLang/Turing.jl/issues/1199
ChainRulesCore.@non_differentiable push!(
    vi::VarInfo,
    vn::VarName,
    r,
    dist::Distribution,
    gidset::Set{Selector}
)

ChainRulesCore.@non_differentiable updategid!(
    vi::AbstractVarInfo,
    vn::VarName,
    spl::Sampler,
)


text:
include("lightweight.jl")
include("incremental.jl")
include("involution_mh.jl")
include("involution_mh_tree.jl")