text:
#using MetidaNCA
using Test
using DataFrames, CSV, Plots
import TypedTables: Table

path     = dirname(@__FILE__)
io       = IOBuffer();
pkdata2  = CSV.File(joinpath(path, "csv", "pkdata2.csv")) |> DataFrame
multtimepk  = CSV.File(joinpath(path, "csv", "multtime.csv")) |> DataFrame
missingpk  = CSV.File(joinpath(path, "csv", "missingpk.csv")) |> DataFrame
aucallpk  = CSV.File(joinpath(path, "csv", "aucalltest.csv")) |> DataFrame
upkdata  = CSV.File(joinpath(path, "csv", "upkdata.csv")) |> DataFrame
pddata  = CSV.File(joinpath(path, "csv", "pddata.csv")) |> DataFrame
include("refdicts.jl")
# Cmax
# Tmax
# Cdose
# Tlag
# Clast
# AUClast
# AUMClast / AUMCtau
# AUCall
# Rsq
# Adjusted Rsq
# Kel
# HL
# LZint
# Clast_pred
# AUCinf
# AUCinf_pred
# AUMCinf
# AUMCinf_pred
# AUCpct
# MRTlast
# MRTinf / MRTtauinf
# MRTinf_pred
# Cllast
# Clinf / Cltau
# Vzlast
# Vzinf / Vztau
# Vssinf

# AUCtau
# Ctau
# Cavg
# Ctaumin
# Accind
# Fluc
# Fluctau
# Swing
# Swingtau
@testset "   Simple test                                             " begin
    # Basic dataset scenario
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 100, time = 0))
    sort!(ds, :Subject)
    show(io, ds)
    dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :lint)
    @test  MetidaNCA.getid(dsnca, :, :Subject) == collect(1:10)
    show(io, dsnca)

    # Export to tables
    mtds  = MetidaNCA.metida_table(ds)
    mtdst = Table(ds)
    @test size(mtds, 1) == size(pkdata2, 1)

    dsncafromds = MetidaNCA.nca(pkdata2, :Time, :Concentration, [:Subject, :Formulation])
    sort!(dsncafromds, :Subject)
    @test dsnca[:, :AUClast] == dsncafromds[:, :AUClast]

    # Plotting
    @test_nowarn MetidaNCA.pkplot(ds; typesort = :Subject, pagesort = nothing, sort = Dict(:Formulation => "R"))
    @test_nowarn MetidaNCA.pkplot(ds; typesort = :Formulation, pagesort = nothing, legend = true)
    @test_nowarn pl = MetidaNCA.pkplot(ds; elim = true, ls = true)
    @test_nowarn pl = MetidaNCA.pkplot(ds; typesort = :Subject, pagesort = :Formulation, elim = true, ls = true, title = "Plots")
    @test_nowarn pl = MetidaNCA.pkplot(ds; typesort = :Formulation, pagesort = :Subject, xticksn = 8, yticksn = 10)
    @test_nowarn pl = MetidaNCA.pkplot(ds; pagesort = [:Subject, :Formulation], legend = false)
    @test_nowarn pl = MetidaNCA.pkplot(ds[1]; ylims = (0, 10), yscale = :log10, legend = false)
    @test_nowarn MetidaNCA.plotstyle(40)

    # setdosetime!
    MetidaNCA.setdosetime!(ds, MetidaNCA.DoseTime(dose = 100, time = 0.25))
    @test first(MetidaNCA.nca!(ds)[:, :Cdose]) == 0

    # Single subject scenario
    tdat = pkdata2[1:16, :Time]
    cdat = pkdata2[1:16, :Concentration]
    ds = MetidaNCA.pkimport(tdat, cdat)
    show(io, ds)
    show(io, MetidaNCA.getdosetime(ds))
    show(io, MetidaNCA.getkelrange(ds))
    sbj = MetidaNCA.nca!(ds)
    show(io, sbj)
    show(io, MetidaNCA.getkeldata(sbj))
    ct = MetidaNCA.ctmax(ds)
    @test  sbj[:Cmax] == ct[1]
    @test  sbj[:Tmax] == ct[2]

    dsncafromds =  MetidaNCA.nca(pkdata2[1:16, :], :Time, :Concentration)
    @test  sbj[:AUClast]  ≈ dsncafromds[:AUClast]

    dsncafromds = MetidaNCA.nca(tdat, cdat)
    @test  sbj[:AUClast]  ≈ dsncafromds[:AUClast]

    # Missing NaN
    dsncafromds =  MetidaNCA.nca(missingpk, :Time, :Concentration)
    @test  sbj[:AUClast]  ≈ dsncafromds[:AUClast]

    dsncafromds =  MetidaNCA.nca(missingpk, :Time, :Concentration;
    limitrule = MetidaNCA.LimitRule(;lloq = 0, btmax = 0, atmax = NaN, nan = NaN, rm = true))
    @test  sbj[:AUClast]  ≈ dsncafromds[:AUClast]

    # Apply modify! function
    function newparam(data)
        data.result[:AUChalf] = data.result[:AUClast] / 2
    end

    dsncafromds =  MetidaNCA.nca(missingpk, :Time, :Concentration;
    limitrule = MetidaNCA.LimitRule(;lloq = 0, btmax = 0, atmax = NaN, nan = NaN, rm = true), modify! = newparam)
    dsncafromds[:AUChalf] ≈ dsncafromds[:AUClast] / 2

    #redirect_stderr(Base.DevNull())
    missingpk.ConcentrationStr = string.(missingpk.Concentration)
    @test_logs (:warn, "Some concentration values not a number, try to fix.") pkiw = MetidaNCA.pkimport(missingpk, :Time, :ConcentrationStr)

end

@testset "  #1 Linear trapezoidal, Dose 100, Dosetime 0, no tau      " begin

ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 100, time = 0))
sort!(ds, :Subject)
dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :lint)

    # Cmax
    @test dsnca[:, :Cmax] == refdict[:Cmax]

    # Tmax
    @test dsnca[:, :Tmax] == refdict[:Tmax]

    # Cdose
    @test round.(dsnca[:, :Cdose], sigdigits = 6) == round.(refdict[:Cdose], sigdigits = 6)

    # Tlag
    #=
    @test round.(dsnca[:, :Tlag], sigdigits = 6) == round.([0
    0
    0
    0
    0.5
    0
    0
    0
    0
    0], sigdigits = 6)
    =#

    # Clast
    @test dsnca[:, :Clast] == refdict[:Clast]

    # AUClast
    @test round.(dsnca[:, :AUClast], sigdigits = 6) == round.(refdict[:AUClast], sigdigits = 6)

    # AUMClast
    @test round.(dsnca[:, :AUMClast], sigdigits = 6) == round.(refdict[:AUMClast], sigdigits = 6)

    # AUCall
    @test round.(dsnca[:, :AUCall], sigdigits = 6) == round.(refdict[:AUCall], sigdigits = 6)

    # Rsq
    @test round.(dsnca[:, :Rsq], digits = 6) == round.(refdict[:Rsq], digits = 6)

    # Adjusted Rsq
    @test round.(dsnca[:, :ARsq], digits = 6) == round.(refdict[:ARsq], digits = 6)

    # Kel
    @test round.(dsnca[:, :Kel], sigdigits = 6) == round.(refdict[:Kel], sigdigits = 6)

    # HL
    @test round.(dsnca[:, :HL], sigdigits = 5) == round.(refdict[:HL], sigdigits = 5)

    # Clast_pred
    @test round.(dsnca[:, :Clast_pred], sigdigits = 6) == round.(refdict[:Clast_pred], sigdigits = 6)

    # AUCinf
    @test round.(dsnca[:, :AUCinf], sigdigits = 6) == round.(refdict[:AUCinf], sigdigits = 6)

    # AUCinf_pred

    # AUMCinf

    # AUMCinf_pred

    # AUCpct
    @test round.(dsnca[:, :AUCpct], sigdigits = 5) == round.(refdict[:AUCpct], sigdigits = 5)

    # MRTlast
    @test round.(dsnca[:, :MRTlast], digits = 6) == round.(refdict[:MRTlast], digits = 6)

    # MRTinf
    @test round.(dsnca[:, :MRTinf], digits = 5) == round.(refdict[:MRTinf], digits = 5)

    # MRTinf_pred

    # Cllast

    # Clinf
    @test round.(dsnca[:, :Clinf], sigdigits = 6) == round.(refdict[:Clinf], sigdigits = 6)

    # Vzlast

    # Vzinf
    @test round.(dsnca[:, :Vzinf], sigdigits = 6) == round.(refdict[:Vzinf], sigdigits = 6)

    # Vssinf

end

@testset "  #2 Linear up Log down, Dose 100, Dosetime 0.25, tau 9    " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 100, time = 0.25, tau = 9))
    sort!(ds, :Subject)
    dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luld)
    # Cmax
    @test dsnca[:, :Cmax] == refdict2[:Cmax]

    # Tmax
    @test dsnca[:, :Tmax] == refdict2[:Tmax]

    # Cdose
    @test round.(dsnca[:, :Cdose], sigdigits = 6) == round.(refdict2[:Cdose], sigdigits = 6)

    # Tlag
    #=
    @test round.(dsnca[:, :Tlag], sigdigits = 6) == round.([0
    0
    0
    0
    0.25
    0
    0
    0
    0
    0], sigdigits = 6)
    =#

    # Clast
    @test dsnca[:, :Clast] == refdict2[:Clast]

    # AUClast
    @test round.(dsnca[:, :AUClast], digits = 6) == round.(refdict2[:AUClast], digits = 6)
    # AUMClast / AUMCtau

    @test round.(dsnca[:, :AUMCtau], digits = 6) ==  round.(refdict2[:AUMCtau], digits = 6)

# AUCall
    @test round.(dsnca[:, :AUCall], digits = 6) == round.(refdict2[:AUCall], digits = 6)
# Rsq
# Adjusted Rsq
    # LZint
    @test round.(dsnca[:, :LZint], digits = 6) == round.(refdict2[:LZint], digits = 6)
# Kel
    @test round.(dsnca[:, :Kel], digits = 6) == round.(refdict2[:Kel], digits = 6)
# HL
    @test round.(dsnca[:, :HL], digits = 6) == round.(refdict2[:HL], digits = 6)
# Clast_pred
    @test round.(dsnca[:, :Clast_pred], digits = 6) == round.(refdict2[:Clast_pred], digits = 6)
# AUCinf
    @test round.(dsnca[:, :AUCinf], digits = 6) == round.(refdict2[:AUCinf], digits = 6)
# AUCinf_pred
# AUMCinf
# AUMCinf_pred
# AUCpct
    @test round.(dsnca[:, :AUCpct], digits = 6) == round.(refdict2[:AUCpct], digits = 6)
# MRTlast

# MRTinf / MRTtauinf
    @test round.(dsnca[:, :MRTtauinf], digits = 6) == round.(refdict2[:MRTtauinf], digits = 6)
# MRTinf_pred
# Cllast
# Clinf / Cltau
    @test round.(dsnca[:, :Cltau], digits = 6) == round.(refdict2[:Cltau], digits = 6)
# Vzlast
# Vzinf / Vztau
    @test round.(dsnca[:, :Vztau], digits = 6) == round.(refdict2[:Vztau], digits = 6)
# Vssinf

    # AUCtau
    @test round.(dsnca[:, :AUCtau], digits = 6) == round.(refdict2[:AUCtau], digits = 6)

end

@testset "  #3 Linear trapezoidal, IV, Dose 120, Dosetime 0.0, tau 12" begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 120, time = 0, tau = 12))
    sort!(ds, :Subject)
    dsnca = MetidaNCA.nca!(ds, adm = :iv, calcm = :lint)

    # Cmax
    @test dsnca[:, :Cmax] == refdict3[:Cmax]

    # Tmax
    @test dsnca[:, :Tmax] == refdict3[:Tmax]

    # Cdose
    @test round.(dsnca[:, :Cdose], sigdigits = 6) == round.(refdict3[:Cdose], sigdigits = 6)

    # Tlag
    #@test round.(dsnca[:, :Tlag], sigdigits = 6) == round.(refdict3[:Tlag], sigdigits = 6)

    # Clast
    @test dsnca[:, :Clast] == refdict3[:Clast]

    # AUClast
    @test round.(dsnca[:, :AUClast], sigdigits = 6) == round.(refdict3[:AUClast], sigdigits = 6)
    # AUMClast / AUMCtau

    @test round.(dsnca[:, :AUMCtau], sigdigits = 6) ==  round.(refdict3[:AUMCtau], sigdigits = 6)

# AUCall
    @test round.(dsnca[:, :AUCall], sigdigits = 6) == round.(refdict3[:AUCall], sigdigits = 6)
# Rsq
# Adjusted Rsq
    # LZint
    @test round.(dsnca[:, :LZint], sigdigits = 6) == round.(refdict3[:LZint], sigdigits = 6)
# Kel
    @test round.(dsnca[:, :Kel], sigdigits = 6) == round.(refdict3[:Kel], sigdigits = 6)
# HL
    @test round.(dsnca[:, :HL], sigdigits = 6) == round.(refdict3[:HL], sigdigits = 6)
# Clast_pred
    @test round.(dsnca[:, :Clast_pred], sigdigits = 6) == round.(refdict3[:Clast_pred], sigdigits = 6)
# AUCinf
    @test round.(dsnca[:, :AUCinf], sigdigits = 6) == round.(refdict3[:AUCinf], sigdigits = 6)
# AUCinf_pred
# AUMCinf
# AUMCinf_pred
# AUCpct
    @test round.(dsnca[:, :AUCpct], sigdigits = 6) == round.(refdict3[:AUCpct], sigdigits = 6)
# MRTlast

# MRTinf / MRTtauinf
    @test round.(dsnca[:, :MRTtauinf], sigdigits = 6) == round.(refdict3[:MRTtauinf], sigdigits = 6)
# MRTinf_pred
# Cllast
# Clinf / Cltau
    @test round.(dsnca[:, :Cltau], sigdigits = 6) == round.(refdict3[:Cltau], sigdigits = 6)
# Vzlast
# Vzinf / Vztau
    @test round.(dsnca[:, :Vztau], sigdigits = 6) == round.(refdict3[:Vztau], sigdigits = 6)
# Vssinf

    # AUCtau
    @test round.(dsnca[:, :AUCtau], sigdigits = 6) == round.(refdict3[:AUCtau], sigdigits = 6)

    # AUClast
    # AUMClast / AUMCtau
    # AUCall
    # Rsq
    # Adjusted Rsq
    # Kel
    # HL
    # LZint
    # Clast_pred
    # AUCinf
    # AUCinf_pred
    # AUMCinf
    # AUMCinf_pred
    # AUCpct
    # MRTlast
    # MRTinf / MRTtauinf
    # MRTinf_pred
    # Cllast
    # Clinf / Cltau
    # Vzlast
    # Vzinf / Vztau
    # Vssinf

    # AUCtau
    # Ctau
    # Cavg
    # Ctaumin
    # Accind
    # Fluc
    # Fluctau
    # Swing
    # Swingtau

end

@testset "  #4 Log trapezoidal ATM, Dose 120, Dosetime 0, tau 12     " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 120, time = 0, tau = 12))
    sort!(ds, :Subject)
    dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :logt)

    # Cmax
    @test dsnca[:, :Cmax] == refdict4[:Cmax]

    # Tmax
    @test dsnca[:, :Tmax] == refdict4[:Tmax]

    # Cdose
    @test round.(dsnca[:, :Cdose], sigdigits = 4) == round.(refdict4[:Cdose], sigdigits = 4)

    # Tlag

    # Clast
    @test dsnca[:, :Clast] == refdict4[:Clast]

    # AUClast
    @test round.(dsnca[:, :AUClast], digits = 4) == round.(refdict4[:AUClast], digits = 4)
    # AUMClast / AUMCtau

    @test round.(dsnca[:, :AUMCtau], digits = 3) ==  round.(refdict4[:AUMCtau], digits = 3)

# AUCall
    @test round.(dsnca[:, :AUCall], digits = 4) == round.(refdict4[:AUCall], digits = 4)
# Rsq
# Adjusted Rsq
    # LZint
    @test round.(dsnca[:, :LZint], digits = 6) == round.(refdict4[:LZint], digits = 6)
# Kel
    @test round.(dsnca[:, :Kel], digits = 6) == round.(refdict4[:Kel], digits = 6)
# HL
    @test round.(dsnca[:, :HL], digits = 3) == round.(refdict4[:HL], digits = 3)
# Clast_pred
    @test round.(dsnca[:, :Clast_pred], digits = 4) == round.(refdict4[:Clast_pred], digits = 4)
# AUCinf
    @test round.(dsnca[:, :AUCinf], digits = 3) == round.(refdict4[:AUCinf], digits = 3)
# AUCinf_pred
# AUMCinf
# AUMCinf_pred
# AUCpct
    @test round.(dsnca[:, :AUCpct], digits = 4) == round.(refdict4[:AUCpct], digits = 4)
# MRTlast

# MRTinf / MRTtauinf
    @test round.(dsnca[:, :MRTtauinf], digits = 4) == round.(refdict4[:MRTtauinf], digits = 4)
# MRTinf_pred
# Cllast
# Clinf / Cltau
    @test round.(dsnca[:, :Cltau], digits = 6) == round.(refdict4[:Cltau], digits = 6)
# Vzlast
# Vzinf / Vztau
    @test round.(dsnca[:, :Vztau], digits = 6) == round.(refdict4[:Vztau], digits = 6)
# Vssinf

    # AUCtau
    @test round.(dsnca[:, :AUCtau], digits = 4) == round.(refdict4[:AUCtau], digits = 4)
    #=
    @test dsnca[:, :Cmax] == [190.869
    261.177
    105.345
    208.542
    169.334
    154.648
    153.254
    138.327
    167.347
    125.482]

    # Tmax
    @test dsnca[:, :Tmax] == [1
    1
    1.5
    1
    4
    2.5
    2.5
    4
    3
    2]

    # Cdose
    @test round.(dsnca[:, :Cdose], sigdigits = 6) == round.([0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0], sigdigits = 6)

    # Tlag
    #=
    @test round.(dsnca[:, :Tlag], sigdigits = 6) == round.([0
    0
    0
    0
    0.5
    0
    0
    0
    0
    0], sigdigits = 6)
    =#

    # Clast
    @test dsnca[:, :Clast] == [112.846
    85.241
    67.901
    97.625
    110.778
    69.501
    58.051
    74.437
    93.44
    42.191]

    # AUClast
    @test round.(dsnca[:, :AUClast], sigdigits = 6) == round.([9572.8582
    10054.0370
    5391.5322
    9296.2179
    9518.6531
    6948.5757
    6987.0645
    7064.7816
    8298.9634
    5485.6538], sigdigits = 6)
    =#
end

@testset "  Linear up Log down ATM, Dose 120, Dosetime 0, tau 12     " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 120, time = 0, tau = 12))
    sort!(ds, :Subject)
    dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luldt)

    # Cmax
    @test dsnca[:, :Cmax] == [190.869
    261.177
    105.345
    208.542
    169.334
    154.648
    153.254
    138.327
    167.347
    125.482]

    # Tmax
    @test dsnca[:, :Tmax] == [1
    1
    1.5
    1
    4
    2.5
    2.5
    4
    3
    2]

    # Cdose
    @test round.(dsnca[:, :Cdose], sigdigits = 6) == round.([0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0], sigdigits = 6)

    # Tlag
    #=
    @test round.(dsnca[:, :Tlag], sigdigits = 6) == round.([0
    0
    0
    0
    0.5
    0
    0
    0
    0
    0], sigdigits = 6)
    =#

    # Clast
    @test dsnca[:, :Clast] == [112.846
    85.241
    67.901
    97.625
    110.778
    69.501
    58.051
    74.437
    93.44
    42.191]

    @test round.(dsnca[:, :AUClast], sigdigits = 6) == round.([9573.810558691312
 10054.286478059563
  5392.457219413793
  9297.096334450325
  9519.181808199797
  6948.985621117448
  6988.960344867885
  7073.306755718137
  8303.373085532965
  5486.838889441992], sigdigits = 6)

  @test round.(dsnca[:, :AUCinf], sigdigits = 6) == round.([42913.407881300096
 16097.04111262767
 26022.090028134975
 21983.33845321829
 25778.19670343543
 15984.14736468632
 11648.33951823218
 15408.84256127436
 24853.53879891218
 7941.268553852982], sigdigits = 6)

end

@testset "  Linear up Log down, Dose 120, Dosetime 0, tau 12         " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 120, time = 0, tau = 12))
    sort!(ds, :Subject)
    dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luld)
    # Cmax
    @test dsnca[:, :Cmax] == [190.869
    261.177
    105.345
    208.542
    169.334
    154.648
    153.254
    138.327
    167.347
    125.482]

    # Tmax
    @test dsnca[:, :Tmax] == [1
    1
    1.5
    1
    4
    2.5
    2.5
    4
    3
    2]

    # Cdose
    @test round.(dsnca[:, :Cdose], sigdigits = 6) == round.([0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0], sigdigits = 6)

    # Tlag
    #=
    @test round.(dsnca[:, :Tlag], sigdigits = 6) == round.([0
    0
    0
    0
    0.5
    0
    0
    0
    0
    0], sigdigits = 6)
    =#

    # Clast
    @test dsnca[:, :Clast] == [112.846
    85.241
    67.901
    97.625
    110.778
    69.501
    58.051
    74.437
    93.44
    42.191]

    # AUClast
    @test round.(dsnca[:, :AUClast], sigdigits = 6) == round.([9573.8106
    10054.286
    5392.4572
    9297.0963
    9519.1809
    6948.9856
    6988.7726
    7073.0922
    8303.3586
    5486.8389], sigdigits = 6)

    # AUMCtau
    @test round.(dsnca[:, :AUMCtau], sigdigits = 6) == round.([9995.4521
14645.001
6030.6644
10314.429
11475.22
8486.2318
9023.7903
6403.6719
10120.542
8376.0054], sigdigits = 6)

    # AUCall
    @test round.(dsnca[:, :AUCall], sigdigits = 6) == round.([9573.8106
10054.286
5392.4572
9297.0963
9519.1809
6948.9856
6988.7726
7073.0922
8303.3586
5486.8389
], sigdigits = 6)

    # Rsq
    @test round.(dsnca[:, :Rsq], digits = 6) == round.([0.78607696
0.99276359
0.81358898
0.91885869
0.85335995
0.95011904
0.97031231
0.94796904
0.94753789
0.88092269], digits = 6)

    # Adjusted Rsq
    @test round.(dsnca[:, :ARsq], digits = 6) == round.([0.71476928
    0.99035145
    0.77630678
    0.83771737
    0.82891994
    0.92517856
    0.96041642
    0.92195356
    0.92130684
    0.86391165], digits = 6)

    # Kel
    @test round.(dsnca[:, :Kel], sigdigits = 6) == round.([0.0033847439
    0.014106315
    0.0032914304
    0.0076953442
    0.0068133279
    0.0076922807
    0.012458956
    0.0089300798
    0.0056458649
    0.017189737], sigdigits = 6)

    # HL
    @test round.(dsnca[:, :HL], sigdigits = 5) == round.([204.78571
    49.137367
    210.59148
    90.073577
    101.73401
    90.10945
    55.634451
    77.619371
    122.77077
    40.323315], sigdigits = 5)

    # Clast_pred
    @test round.(dsnca[:, :Clast_pred], sigdigits = 6) == round.([117.30578
    82.53669
    66.931057
    100.76793
    105.29832
    71.939942
    61.172702
    75.604277
    93.761762
    38.810857], sigdigits = 6)

    # AUCinf
    @test round.(dsnca[:, :AUCinf], sigdigits = 6) == round.([42913.408
16097.041
26022.09
21983.338
25778.196
15984.147
11648.152
15408.628
24853.524
7941.2686], sigdigits = 6)

    # AUCinf_pred

    # AUMCinf

    # AUMCinf_pred

    # AUCpct
    @test round.(dsnca[:, :AUCpct], sigdigits = 5) == round.([77.690398
37.539537
79.277386
57.708442
63.072742
56.525766
40.001017
54.096548
66.59082
30.907274], sigdigits = 5)

    # MRTlast
    #=
    @test round.(dsnca[:, :MRTlast], digits = 6) == round.([34.801023
    29.538786
    34.472406
    33.69408
    32.964438
    32.58076
    31.267574
    33.826053
    33.386807
    27.556657], digits = 6)
=#
    # MRTtauinf
    @test round.(dsnca[:, :MRTtauinf], digits = 4) == round.([302.522490
75.321653
312.895820
148.293830
172.561490
130.198390
91.786365
163.779880
176.558310
68.172066], digits = 4 )

    # MRTinf_pred

    # Cllast

    # Cltau
    @test round.(dsnca[:, :Cltau], sigdigits = 6) == round.([0.071896833
0.050424060
0.122488280
0.070172356
0.069035042
0.085192949
0.083697775
0.110260280
0.073377834
0.092799594], sigdigits = 6)

    # Vzlast

    # Vztau
    @test round.(dsnca[:, :Vztau], sigdigits = 6) == round.([21.241439
3.574574
37.214301
9.118807
10.132353
11.075122
6.717880
12.347065
12.996739
5.398547], sigdigits = 6)

    # Vssinf
end

@testset "  Linear trapezoidal, Dose 100, Dosetime 2.0, tau 10       " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 100, time = 2.0, tau = 10))
    sort!(ds, :Subject)
    dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :lint)

    # Cmax
    @test dsnca[:, :Cmax] == [150.854
    234.091
    100.943
    165.177
    169.334
    154.648
    153.254
    138.327
    167.347
    125.482]

    # Tmax
    @test dsnca[:, :Tmax] == [4
    2
    2
    2
    4
    2.5
    2.5
    4
    3
    2]

    # Cdose
    # Tlag
    # Clast
    # AUClast

    @test round.(dsnca[:, :AUClast], digits = 3) == round.([9283.0588
    9774.9218
    5255.0335
    9051.5985
    9441.4205
    6781.2608
    6858.7153
    6885.9150
    8094.8320
    5495.0990], digits = 3)

    # AUMClast / AUMCtau
    @test round.(dsnca[:, :AUMCtau], digits = 3) == round.([6915.4913
    10105.3000
    4166.2103
    7068.4668
    8032.4660
    5775.6840
    6243.2670
    4442.4953
    6999.9440
    5849.5573], digits = 3)
    # AUCall
    # Rsq
    # Adjusted Rsq
    # Kel
    # HL
    # LZint
    @test round.(dsnca[:, :LZint], sigdigits = 6) == round.([5.0084856
    5.4288976
    4.4406461
    5.166885
    5.1473571
    4.8296758
    5.0107459
    4.9684786
    4.9472594
    4.8963611], sigdigits = 6)
    # Clast_pred
    @test round.(dsnca[:, :Clast_pred], digits = 3) == round.([117.30578
    82.53669
    66.931057
    100.76793
    105.29832
    71.939942
    61.172702
    75.604277
    93.761762
    38.810857], digits = 3)
    # AUCinf
    @test round.(dsnca[:, :AUCinf], digits = 3) == round.([42622.6560
    15817.6760
    25884.6660
    21737.8410
    25700.4350
    15816.4220
    11518.0940
    15221.4510
    24644.9980
    7949.5287], digits = 3)
    # AUCinf_pred
    # AUMCinf
    # AUMCinf_pred
    # AUCpct
    # MRTlast
    # MRTinf / MRTtauinf
    @test round.(dsnca[:, :MRTtauinf], digits = 3) == round.([306.68478
    72.36944
    303.54748
    145.34822
    153.74241
    123.86557
    85.934338
    167.95643
    168.74491
    63.074784], digits = 3)

    # MRTinf_pred
    # Cllast
    # Clinf / Cltau
    # Vzlast
    # Vzinf / Vztau
    # Vssinf

    # AUCtau
    @test round.(dsnca[:, :AUCtau], digits = 3) == round.([1367.7388
    2043.0158
    838.8295
    1444.7985
    1618.6205
    1224.6608
    1265.7013
    880.311
    1417.942
    1167.911], digits = 3)
    # Ctau
    @test round.(dsnca[:, :Ctau], sigdigits = 6) == round.([144.964
    196.035
    76.027
    132.257
    154.066
    113.751
    123.37
    134.133
    135.58
    106.476], sigdigits = 6)
    # Cavg
    @test round.(dsnca[:, :Cavg], digits = 3) == round.([136.77388
    204.30158
    83.88295
    144.47985
    161.86205
    122.46608
    126.57013
    88.0311
    141.7942
    116.7911], digits = 3)
    # Ctaumin
    @test round.(dsnca[:, :Ctaumin], sigdigits = 6) == round.([121.239
    196.026
    69.985
    128.398
    151.452
    102.16
    105.513
    22.724
    105.438
    106.476], sigdigits = 6)
    # Accind
    @test round.(dsnca[:, :Accind], digits = 3) == round.([30.047153
    7.600775
    30.884671
    13.501282
    15.182793
    13.506455
    8.5367345
    11.705549
    18.216783
    6.3317425], digits = 3)
    # Fluc
    @test round.(dsnca[:, :Fluc], sigdigits = 6) == round.([21.652527
    18.631770
    36.906189
    25.456145
    11.047679
    42.859216
    37.719011
    131.320640
    43.661165
    16.273500], sigdigits = 6)

    # Swing
    @test round.(dsnca[:, :Swing], sigdigits = 6) == round.([0.24426958
    0.19418342
    0.44235193
    0.28644527
    0.11807041
    0.51378230
    0.45246557
    5.08726460
    0.58716023
    0.17850032], sigdigits = 6)

end

@testset "  Linear trapezoidal, Dose 100, Dosetime 0.0, tau 100      " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 100, time = 0.0, tau = 100))
    sort!(ds, :Subject)
    dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luld)

    # Cmax
    @test dsnca[:, :Cmax] == [190.869
    261.177
    105.345
    208.542
    169.334
    154.648
    153.254
    138.327
    167.347
    125.482]

    # Tmax
    @test dsnca[:, :Tmax] == [1
    1
    1.5
    1
    4
    2.5
    2.5
    4
    3
    2]

    # Cdose
    # Tlag
    # Clast
    # AUClast

    @test round.(dsnca[:, :AUClast], digits = 3) == round.([9573.810559
    10054.28648
    5392.457219
    9297.096334
    9519.180874
    6948.985621
    6988.772632
    7073.092214
    8303.358586
    5486.838889], digits = 3)

    # AUMClast / AUMCtau
    @test round.(dsnca[:, :AUMCtau], digits = 3) == round.([599491.4507
    469656.5157
    341841.895
    530699.316
    554228.2569
    381802.405
    343144.4367
    402118.7541
    487552.4608
    233710.5713], digits = 3)
    # AUCall
    # Rsq
    # Adjusted Rsq
    # Kel
    # HL
    # LZint
    @test round.(dsnca[:, :LZint], sigdigits = 6) == round.([5.0084856
    5.4288976
    4.4406461
    5.166885
    5.1473571
    4.8296758
    5.0107459
    4.9684786
    4.9472594
    4.8963611], sigdigits = 6)

    # Clast_pred
    @test round.(dsnca[:, :Clast_pred], digits = 3) == round.([117.30578
    82.53669
    66.931057
    100.76793
    105.29832
    71.939942
    61.172702
    75.604277
    93.761762
    38.810857], digits = 3)
    # AUCinf
    @test round.(dsnca[:, :AUCinf], digits = 3) == round.([42913.4078813004
    16097.0411126277
    26022.0900281352
    21983.3384532182
    25778.1957695968
    15984.1473646863
    11648.1518057779
    15408.6280190766
    24853.5242997128
    7941.2685538530], digits = 3)
    # AUCinf_pred
    # AUMCinf
    # AUMCinf_pred
    # AUCpct
    # MRTlast
    # MRTinf / MRTtauinf
    @test round.(dsnca[:, :MRTtauinf], digits = 3) == round.([286.7298872
    73.32812486
    309.1287691
    131.3893596
    155.1600979
    126.8510359
    79.61736185
    117.5547609
    177.1315649
    60.86551733], digits = 3)

    # MRTinf_pred
    # Cllast
    # Clinf / Cltau
    # Vzlast
    # Vzinf / Vztau
    # Vssinf

    # AUCtau
    @test round.(dsnca[:, :AUCtau], digits = 3) == round.([12646.63632
    11996.6718
    7195.902904
    11794.11692
    12274.83395
    8729.151856
    8395.400098
    8930.999936
    10727.4135
    6389.420453], digits = 3)
    # Ctau
    @test round.(dsnca[:, :Ctau], sigdigits = 6) == round.([106.6989373
    55.60460942
    61.0383917
    81.23535402
    87.01010737
    58.00027625
    43.15724396
    58.8781831
    80.05171762
    23.98401112], sigdigits = 6)
    # Cavg
    @test round.(dsnca[:, :Cavg], digits = 3) == round.([126.4663632
    119.966718
    71.95902904
    117.9411692
    122.7483395
    87.29151856
    83.95400098
    89.30999936
    107.274135
    63.89420453], digits = 3)
    # Ctaumin
    @test round.(dsnca[:, :Ctaumin], sigdigits = 6) == round.([0
    0
    0
    0
    0
    0
    0
    0
    0
    0], sigdigits = 6)
    # Accind
    @test round.(dsnca[:, :Accind], digits = 3) == round.([3.482585727
1.322732351
3.565571971
1.862990767
2.024054788
1.863483514
1.403869629
1.693257481
2.318008606
1.218397838], digits = 3)
    # Fluc
    @test round.(dsnca[:, :Fluc], sigdigits = 6) == round.([150.9247164
217.7078813
146.3958052
176.8186643
137.9521716
177.1626872
182.5452012
154.8841126
155.9993935
196.3902687], sigdigits = 6)

    # Swing

    # Swing tau
    @test round.(dsnca[:, :Swingtau], sigdigits = 6) == round.([0.788855679643705
3.697038657747540
0.725880991766796
1.567133516128610
0.946141719805552
1.666332127921700
2.551060863661350
1.349376164748240
1.090486062900490
4.231902178181850], sigdigits = 6)

end

@testset "  Linear up Log down, Dose 100, Dosetime 0.25, tau 9 IV    " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 100, time = 0.25, tau = 9))
    sort!(ds, :Subject)
    dsnca = MetidaNCA.nca!(ds, adm = :iv, calcm = :luld)

    @test dsnca[:, :Cdose] == [178.949
    62.222
    49.849
    52.421
    0.0
    57.882
    19.95
    142.34985100539438
    113.362
    13.634]

end

@testset "  Linear trapezoidal, Dose 100, Dosetime 0, no tau AUCall  " begin

    dsnca = MetidaNCA.nca(aucallpk, :Time, :Concentration; dosetime = MetidaNCA.DoseTime(dose = 100, time = 0), adm = :ev, calcm = :lint)
    @test dsnca[:AUClast] ≈ 9585.4218
    aucl =  MetidaNCA.linauc(72, 96, 112.846, 0)
    dsnca[:AUClast] + aucl ≈ dsnca[:AUCall]

end

@testset "  set-get*! tests                                          " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation])
    sort!(ds, :Subject)
    #=
    @testset "  #1 setkelauto!                                            " begin
        ka = MetidaNCA.setkelauto!(ds[1], false)
        @test MetidaNCA.getkelauto(ka) == true

        dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luldt)
    end
    =#
    @testset "  setdosetime!                                             " begin
        dt = MetidaNCA.DoseTime(dose = 110, time = 2.1, tau = 10)
        MetidaNCA.setdosetime!(ds[1], dt)
        dts = MetidaNCA.getdosetime(ds[1])
        @test dts.dose == 110
        @test dts.time == 2.1
        @test dts.tau == 10
        dt2 = MetidaNCA.DoseTime(dose = 100, time = 2.2, tau = 9)
        MetidaNCA.setdosetime!(ds, dt2, 4)
        MetidaNCA.setdosetime!(ds, dt2, [1,2,3])
        MetidaNCA.setdosetime!(ds, dt2, Dict(:Formulation => "R"))
        MetidaNCA.setdosetime!(ds, dt2)

        dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luldt)
    end
    @testset "  #2 setkelauto!                                            " begin
        kr =  MetidaNCA.ElimRange(kelstart = 4, kelend = 12)
        MetidaNCA.setkelrange!(ds, kr; kelauto = true)
        MetidaNCA.setkelauto!(ds, false, 4)
        MetidaNCA.setkelauto!(ds, false, [1,2,3])
        MetidaNCA.setkelauto!(ds, false, Dict(:Formulation => "R"))
        MetidaNCA.setkelauto!(ds, false)
        dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luldt)
        @test MetidaNCA.getkelauto(ds[1]) == false
    end
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation])
    sort!(ds, :Subject)
    @testset "  setkelrange!                                             " begin
        dsnca1 = deepcopy(MetidaNCA.nca!(ds[1], adm = :ev, calcm = :luldt))
        kr =  MetidaNCA.ElimRange(kelstart = 12, kelend = 16)
        MetidaNCA.setkelrange!(ds[1], kr)
        dsnca2 = deepcopy(MetidaNCA.nca!(ds[1], adm = :ev, calcm = :luldt))
        @test dsnca1.data.keldata.ar[3] ≈ 0.7147692761075757
        @test dsnca1.data.keldata.ar[3] ≈ dsnca2.data.keldata.ar[1]
        @test dsnca1.data.keldata.a[3] ≈ dsnca2.data.keldata.a[1]
        @test dsnca1.data.keldata.b[3] ≈ dsnca2.data.keldata.b[1]

        kr =  MetidaNCA.ElimRange(kelstart = 12, kelend = 16, kelexcl = Int[5,6])
        MetidaNCA.setkelrange!(ds[1], kr)
        dsnca3 = deepcopy(MetidaNCA.nca!(ds[1], adm = :ev, calcm = :luldt))
        @test dsnca1.data.keldata.ar[3] ≈ dsnca3.data.keldata.ar[1]
        @test dsnca1.data.keldata.a[3] ≈ dsnca3.data.keldata.a[1]
        @test dsnca1.data.keldata.b[3] ≈ dsnca3.data.keldata.b[1]

        kr =  MetidaNCA.ElimRange(kelexcl = Int[5,6])
        MetidaNCA.setkelrange!(ds[1], kr; kelauto = true)
        dsnca4 = deepcopy(MetidaNCA.nca!(ds[1], adm = :ev, calcm = :luldt))
        @test dsnca1.data.keldata.ar[3] ≈ dsnca4.data.keldata.ar[3]
        @test dsnca1.data.keldata.a[3] ≈ dsnca4.data.keldata.a[3]
        @test dsnca1.data.keldata.b[3] ≈ dsnca4.data.keldata.b[3]

        kr =  MetidaNCA.ElimRange(kelstart = 4, kelend = 12, kelexcl = Int[5,6])
        MetidaNCA.setkelrange!(ds[1], kr)
        krs = MetidaNCA.getkelrange(ds[1])
        @test krs.kelstart == 4
        @test krs.kelend == 12
        @test krs.kelexcl == [5,6]
        kr2 =  MetidaNCA.ElimRange(kelstart = 3, kelend = 12, kelexcl = Int[7])
        MetidaNCA.setkelrange!(ds, kr2, 4)
        MetidaNCA.setkelrange!(ds, kr2, [1,2,3])
        MetidaNCA.setkelrange!(ds, kr2, Dict(:Formulation => "R"))
        MetidaNCA.setkelrange!(ds, kr2)

        dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luldt)
    end
end
@testset "  applylimitrule!                                          " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation])
    sort!(ds, :Subject)
    lr = MetidaNCA.LimitRule(;lloq = 0.5, btmax = 0.0, atmax = NaN, nan = NaN, rm = true)
    MetidaNCA.applylimitrule!(ds[1], lr)
    dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luldt)

    lr = MetidaNCA.LimitRule(;lloq = 0.5, btmax = 0.5, atmax = NaN, nan = NaN, rm = true)
    function af(sbj)
        sbj.id[:Subject] == 1
    end
    MetidaNCA.applylimitrule!(af, ds, lr)
    @test ds[1].obs[1] ≈ 0.5
    MetidaNCA.applylimitrule!(ds, lr, 2)
    @test ds[2].obs[1] ≈ 0.5
    MetidaNCA.applylimitrule!(ds, lr, 3:4)
    @test ds[3].obs[1] ≈ 0.5
    @test ds[4].obs[1] ≈ 0.5
    MetidaNCA.applylimitrule!(ds, lr, Dict(:Formulation => "R"))
    @test ds[7].obs[1] ≈ 0.5
    MetidaNCA.applylimitrule!(ds, lr)
    @test ds[6].obs[1] ≈ 0.5

    ds = MetidaNCA.pkimport(missingpk, :Time, :Concentration)

    @test ismissing(ds.obs[13])
    @test isnan(ds.obs[15])
    @test length(ds) == 18
    MetidaNCA.applylimitrule!(ds, lr)
    @test length(ds) == 16

    ds = MetidaNCA.pkimport(missingpk, :Time, :Concentration)
    lr = MetidaNCA.LimitRule(;lloq = 180, btmax = 0.0, atmax = 0.5, nan = 1000, rm = false)
    MetidaNCA.applylimitrule!(ds, lr)
    @test ds.obs ==  [0.0
    0.0
  190.869
    0.5
    0.5
    0.5
    0.5
    0.5
    0.5
    0.5
    0.5
    0.5
 1000.0
    0.5
 1000.0
    0.5
    0.5
    0.5]
end

@testset "  kel                                                      " begin
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation])
    sort!(ds, :Subject)
    kr1 = MetidaNCA.ElimRange(kelstart = 4, kelend = 12, kelexcl = Int[5,6])
    kr2 = MetidaNCA.ElimRange(kelstart = 3, kelend = 12, kelexcl = Int[7])

    MetidaNCA.setkelrange!(ds, kr1, Dict(:Formulation => "T"))
    MetidaNCA.setkelrange!(ds, kr2, Dict(:Formulation => "R"))

    sub1 = MetidaNCA.subset(ds, Dict(:Formulation => "T"))
    @test MetidaNCA.getkelrange(sub1[1]) == kr1
    sub2 = MetidaNCA.subset(ds, Dict(:Formulation => "R"))
    @test MetidaNCA.getkelrange(sub2[1]) == kr2
    dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luldt)
end

@testset "  Output                                                   " begin
    io = IOBuffer();
    ds = MetidaNCA.pkimport(pkdata2, :Time, :Concentration, [:Subject, :Formulation]; dosetime = MetidaNCA.DoseTime(dose = 100, time = 0))
    sort!(ds, :Subject)
    @test_nowarn dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :lint, verbose = 1, io = io)
    show(io, ds[1])
    dt = MetidaNCA.DoseTime(dose = 100, time = 0.25, tau = 9)
    MetidaNCA.setdosetime!(ds, dt)
    @test_nowarn dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :luld, verbose = 1, io = io)
    show(io, ds[1])
    kr =  MetidaNCA.ElimRange(kelstart = 10, kelend = 16, kelexcl = Int[13,14])
    MetidaNCA.setkelrange!(ds, kr; kelauto = false)
    show(io, ds[1])
    @test_nowarn dsnca = MetidaNCA.nca!(ds, adm = :ev, calcm = :lint, verbose = 2, io = io)
end

@testset "  Multiple time                                            " begin
    io = IOBuffer();
    @test_logs (:warn,"Not all time values is unique, last observation used! ((1,))") ds = MetidaNCA.pkimport(multtimepk, :Time, :Concentration, :Subject)
    #@test ds[1].obs[6] == 129.59

end

include("upktest.jl")
include("pdtest.jl")


text:
"""
    OlsSureFn(y,x,m=1)

LS of y on x; for one n dependent variables, same regressors

# Usage
(b,res,yhat,Covb,R2a) = OlsSureFn(y,x,m)

# Input
- `y::Array`:     Txn, the n dependent variables
- `x::Array`:     Txk matrix of regressors (including deterministic ones)
- `m::Int`:       scalar, bandwidth in Newey-West

# Output
- `b::Array`:     n*kx1, regression coefficients
- `u::Array`:     Txn, residuals y - yhat
- `yhat::Array`:  Txn, fitted values x*b
- `V::Array`:     matrix, covariance matrix of sqrt(T)vec(b)
- `R2a::Number`:  n vector, R2 value

"""
function OlsSureFn(y,x,m=0)
    (T,n) = (size(y,1),size(y,2))
    k     = size(x,2)
    b     = x\y
    yhat  = x*b
    u     = y - yhat
    g     = zeros(T,n*k)
    for i = 1:n
      vv      = (1+(i-1)*k):(i*k)   #1:k,(1+k):2k,...
      g[:,vv] = x.*u[:,i]           #moment conditions for y[:,i] regression
    end
    S0    = NWFn(g,m)            #Newey-West covariance matrix
    Sxxi  = -x'x/T
    Sxx_1 = kron(Matrix(1.0I,n,n),inv(Sxxi))    #Matrix(1.0I,n,n) is identity matrix(n)
    V     = Sxx_1 * S0 * Sxx_1
    R2a   = 1.0 .- Compat.var(u,dims=1)./Compat.var(y,dims=1)  #0.7 syntax
    return b,u,yhat,V,R2a
end

ext: jl

max_stars_repo_name: snowdj/EmpiricalFinancePhD

max_issues_repo_name: AlTaken/EmpiricalFinancePhD

max_forks_repo_name: AlTaken/EmpiricalFinancePhD


text:
struct GradientNull <: Gradient end
struct Gradient0 <: Gradient
    index::Int
    offset::Int
end
struct Gradient1{F,D1 <: Gradient} <: Gradient
    d1::D1
    fref1::MyRef{Float64}
    ref::MyRef{Float64}
    function Gradient1(e::Expression1{F,E1}, indexer = nothing) where {F,E1}
        d1 = Gradient(e.e1,indexer)
        return new{F,typeof(d1)}(d1,ref(e.e1),MyRef(1.))
    end
end
struct Gradient2F1{F,D1 <: Gradient, R<: Real} <: Gradient
    a::R
    d1::D1
    fref1::MyRef{Float64}
    ref::MyRef{Float64}
    function Gradient2F1(e::Expression2{F,E1,E2}, indexer = nothing) where {F,E1<:Real,E2}
        g1 = Gradient(e.e2,indexer)
        return new{F,typeof(g1),typeof(e.e1)}(e.e1,g1,ref(e.e2),MyRef(0.))
    end
end
struct Gradient2F2{F,D1 <: Gradient, R<: Real} <: Gradient
    a::R
    d1::D1
    fref1::MyRef{Float64}
    ref::MyRef{Float64}
    function Gradient2F2(e::Expression2{F,E1,E2}, indexer = nothing) where {F,E1,E2<:Real}
        g1 = Gradient(e.e1,indexer)
        return new{F,typeof(g1),typeof(e.e2)}(e.e2,g1,ref(e.e1),MyRef(0.))
    end
end
struct Gradient2{F,D1 <: Gradient,D2 <: Gradient} <: Gradient
    d1::D1
    d2::D2
    fref1::MyRef{Float64}
    fref2::MyRef{Float64}
    ref1::MyRef{Float64}
    ref2::MyRef{Float64}
    function Gradient2(e::Expression2{F,E1,E2}, indexer = nothing) where {F,E1,E2}
        g1 = Gradient(e.e1,indexer)
        g2 = Gradient(e.e2,indexer)
        return new{F,typeof(g1),typeof(g2)}(g1,g2,ref(e.e1),ref(e.e2),MyRef(0.),MyRef(0.))
    end
end
struct GradientSum{D,I} <: Gradient
    inner::I
    ds::Vector{D}
end

const GRADIENT_NULL = GradientNull()

@inline (::GradientNull)(z,x,p=nothing,d0=1) = nothing
@inline (d::Gradient0)(y,x,p=nothing,d0=1) = (@inbounds y[d.offset] += d0; nothing)
@inline (d::Gradient0)((j,y)::Tuple{Int,M},x,p=nothing,d0=1) where M <: AbstractMatrix = (@inbounds y[j,d.offset] += d0; nothing)
@inline (d::Gradient0)((j,y)::Tuple{Int,M},x,p=nothing,d0=1) where M <: AbstractVector = (@inbounds y[d.offset] += d0; nothing)

@inline function (d::GradientSum{D,I})(y,x,p=nothing,d0=1) where {D,I}
    inner(d)(y,x,p,d0)
    @simd for i in eachindex(d.ds)
        @inbounds d.ds[i](y,x,p,d0)
    end
end
@inline function (d::GradientSum{D,Nothing})(y,x,p=nothing,d0=1) where D
    @simd for i in eachindex(d.ds)
        @inbounds d.ds[i](y,x,p,d0)
    end
end

Gradient(e::V,(row,indexer)::Tuple{Int,Dict{Tuple{Int,Int},Int}}) where V <: Variable = Gradient0(index(e),set_indexer!(indexer,row,index(e)))
Gradient(e::V,indexer=nothing) where V <: Variable = Gradient0(index(e),set_indexer!(indexer,index(e)) )
Gradient(e::V,::Nothing) where V <: Variable = Gradient0(index(e),index(e))
Gradient(e::V,::Tuple{Int,Nothing}) where V <: Variable = Gradient0(index(e),index(e))
Gradient(e::P,::T) where {P<:Parameter,T} = GRADIENT_NULL
Gradient(e::C, indexer = nothing ) where {C<:Constant,T} = GRADIENT_NULL
Gradient(e::ExpressionSum{E,I},indexer = nothing) where {E,I} = GradientSum(Gradient(inner(e),indexer),[Gradient(ee,indexer) for ee in e.es])
Gradient(e::ExpressionSum{E,Nothing},indexer = nothing) where {E,I} = GradientSum(nothing,[Gradient(ee,indexer) for ee in e.es])
Gradient(e::Expression1{F,E}, indexer = nothing) where {F,E} = Gradient1(e,indexer)
Gradient(e::Expression2{F,E1,E2}, indexer = nothing) where {F,E1,E2} = Gradient2(e,indexer)
Gradient(e::Expression2{F,E1,E2}, indexer = nothing) where {F,E1<:Real,E2} = Gradient2F1(e,indexer)
Gradient(e::Expression2{F,E1,E2}, indexer = nothing) where {F,E1,E2<:Real} = Gradient2F2(e,indexer)


text:
module ShallowArrays

using Swizzles.WrapperArrays

using Base.Broadcast: BroadcastStyle
using Base: dataids, unaliascopy, unalias

export ShallowArray

"""
    ShallowArray

A convenience type for constructing simple wrapper arrays that behave almost
exactly like their parent. Provides default implementations of `AbstractArray`
methods. Subtypes of `ShallowArray` must define [`parent`](@ref) and
[`adopt`](@ref)

See also: [`parent`](@ref), [`adopt`](@ref)
"""
abstract type ShallowArray{T, N, Arg} <: AbstractArray{T, N} end

Base.parent(arr::ShallowArray) = throw(MethodError(parent, (arr)))
iswrapper(arr::ShallowArray) = true

IndexStyle(arr::ShallowArray) = IndexStyle(parent(arr))

Base.dataids(arr::ShallowArray) = dataids(parent(arr))
Base.unaliascopy(arr::A) where {A <:ShallowArray} = adopt(arr, unaliascopy(parent(arr)))::A
Base.unalias(dest, arr::A) where {A <:ShallowArray} = adopt(arr, unalias(dest, parent(arr)))::A

Base.eltype(::Type{<:ShallowArray{T}}) where {T} = T
Base.eltype(::ShallowArray{T}) where {T} = T

Base.ndims(::Type{<:ShallowArray{<:Any, N}}) where {N} = N
Base.ndims(::ShallowArray{<:Any, N}) where {N} = N

Base.size(arr::ShallowArray{<:Any, <:Any, <:AbstractArray}) = size(parent(arr))

Base.axes(arr::ShallowArray{<:Any, <:Any, <:AbstractArray}) = axes(parent(arr))

Base.@propagate_inbounds Base.getindex(arr::ShallowArray, inds...) = getindex(parent(arr), inds...)

Base.@propagate_inbounds Base.setindex!(arr::ShallowArray, val, inds...) = setindex!(parent(arr), val, inds...)

Base.eachindex(arr::ShallowArray) = eachindex(parent(arr))

@inline WrapperArrays.childstyle(::Type{<:ShallowArray}, S::BroadcastStyle) = S

@inline function Broadcast.BroadcastStyle(Arr::Type{<:ShallowArray{<:Any, <:Any, Arg}}) where {Arg}
    childstyle(Arr, BroadcastStyle(Arg))
end

function Base.show(io::IO, arr::ShallowArray{T, N, Arg}) where {T, N, Arg}
    print(io, typeof(arr))
    print(io, '(', parent(arr), ')')
    nothing
end

end


text:
# ===========
print("Test de l'analyse des arguments passées en ligne de commande... ")

# Si le premier paramètre est inconnu (action "xxx"), on doit lever une erreur
@test_throws ArgumentError args = Args.parse_commandline("xxx")

# Si le fichier d'instance n'existe pas, on doit lever une erreur
@test_throws ArgumentError args = Args.parse_commandline("test --infile fichier_bidon.alp")

# Test de l'analyse de la chaine passée en arguments
cli = "
    explore
    --infile $APPDIR/test/data/alp_01.alpx
    --loglevel 2
    --presort shuffle
    -L0
    -n100
    -d/tmp
"
args = Args.parse_commandline(cli)

# On vérifier que to_s_dict(args) retourne bien le dict des arguments
# Args.show_dict(args)
# println(Args.to_s_dict(args))
@test occursin(r"\spresort\s*=>\s*shuffle\s", Args.to_s_dict(args))

# Test d'accès aux arguments via le dict
@test  args[:infile] == "$APPDIR/test/data/alp_01.alpx"
@test  args[:level] == 0

# Test d'accès aux arguments via l'accesseur
@test  Args.get(:level) == args[:level]
@test  Args.get("level") == args[:level]

# Test de la modification d'un arguments
Args.set(:level, 3)
@test  args[:level] == 3
@test  Args.get(:level) == args[:level]
@test  Args.get("level") == args[:level]

Args.set(:level, 0) # Pour rendre les tests suivants silencieux !

print("ok\n")


text:
const PFCUT = [6:5:101]//100
const PFVAL = [10:8:26, 32:6:50, 54:4:98, 100]//100

function pricefraction{T<:FloatingPoint}(a::T)
    zero(T) <= a || error("a = ", a, ", but it must be >= 0.")
    a <= one(T) || error("a = ", a, ", but it must be <= 1.")
    convert(T, PFVAL[findfirst(a .< PFCUT)])
end

test = [0.:0.05:1., 0.51, 0.56, 0.61, rand(), rand(), rand(), rand()]

println("Testing the price fraction function")
for t in test
    println(@sprintf "    %.4f -> %.4f" t pricefraction(t))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module ShapleyValues

export shapleyvalues

include("utils.jl")
include("coalitions.jl")
include("samplespaces.jl")

"The core method that updates the Shapley value estimates."
function update_estimates!(deltas, csets, x, f, X, sampleWeights,
                           g, featureGroups::Array{Array{Int64,1}},
                           sampleCounts::Array{Int64,1}, sampleSpace)
    P,N = size(X)
    M = length(featureGroups)
    @assert length(sampleCounts) == M "sampleCounts should be an array of counts for each feature group!"

    # build the synthentic samples
    inds = collect(1:M)
    s1 = Array(eltype(X), P)
    s2 = Array(eltype(X), P)
    weights = Array(eltype(X), sum(sampleCounts))
    reset!(sampleSpace)
    pos = 1
    wpos = 1
    for i in 1:M
        for j in 1:sampleCounts[i]
            cset,w = next(csets[i])
            weights[wpos] = w
            wpos += 1

            for k in 1:N

                # save two synthetic samples with and without the current group replaced
                copy!(s2, X[:,k])
                for l in cset
                    for m in featureGroups[l]
                        s2[m] = x[m]
                    end
                end
                copy!(s1, s2)
                for m in featureGroups[i]
                    s1[m] = x[m]
                end
                addsample!(sampleSpace, s1)
                addsample!(sampleSpace, s2)
                pos += 2
            end
        end
    end

    # run the provided function
    y::Array{eltype(X),1} = vec(f(data(sampleSpace)))

    # sum the totals and keep an estimate of the variance differences
    pos = 1
    wpos = 1
    withiExp = zero(eltype(X))
    for i in 1:M
        for j in 1:sampleCounts[i]
            withiExp = zero(eltype(X))
            withoutiExp = zero(eltype(X))
            sumw = zero(eltype(X))
            for k in 1:N
                w = sampleWeights[k]
                withiExp += w*y[pos]
                withoutiExp += w*y[pos+1]
                sumw += w
                pos += 2
            end
            withiExp /= sumw
            withoutiExp /= sumw
            observe!(deltas[i], g(withiExp) - g(withoutiExp), weights[wpos])
            wpos += 1
        end
    end
end

"Designed to determine the Shapley values (importance) of each feature for f(x)."
function shapleyvalues(x, f::Function, X, g::Function=identity; featureGroups=nothing, synthSampleSpace=nothing,
                       weights=nothing, nsamples=nothing)
    P,N = size(X)
    @assert length(x) == P "Provided 'x' length must match the data matrix features count ($(length(x)) != $P)!"

    # give default values to omitted arguments
    weights != nothing || (weights = ones(N))
    weights ./= sum(weights)
    featureGroups != nothing || (featureGroups = Array{Int64,1}[Int64[i] for i in 1:length(x)])
    featureGroups = convert(Array{Array{Int64,1},1}, featureGroups)
    @assert length(weights) == N "Provided 'weights' must match the number of representative data points (size(X)[2])!"

    # find the feature groups we will test. If a feature does not change from its
    # current value then we know it doesn't impact the model
    varyingInds = varying_groups(x, X, featureGroups)
    varyingFeatureGroups = featureGroups[varyingInds]
    M = length(varyingFeatureGroups)

    # more default values
    nsamples != nothing || (nsamples = 2M+1000)
    synthSampleSpace != nothing || (synthSampleSpace = samplespace(X, nsamples))
    @assert nsamples >= 2M "'nsamples' must be at least 2 times the number of varying feature groups!"

    # loop through the estimation process focusing samples on groups with high variance
    deltas = [MeanVarianceAccumulator() for i in 1:M]
    csets = [statefuliterator(coalitions(setdiff(collect(1:M), [i]))) for i in 1:M]
    numSets = length(csets[1])
    nsamples = min(nsamples, numSets*M)
    totalSamples = zeros(Int64, M)
    nextSamples = allocate_samples(ones(M), min(2M, nsamples), numSets*ones(M))
    while true

        # update our estimates for a block of samples
        update_estimates!(deltas, csets, x, f, X, weights, g, varyingFeatureGroups, nextSamples, synthSampleSpace)

        # keep track of our samples and optimize their allocation to minimize variance (Neyman allocation)
        totalSamples .+= nextSamples
        samplesLeft = nsamples-sum(totalSamples)
        if samplesLeft > 0
            # (N-n)/(N-1) is the finite sample variance correction without weighting
            # so we approximate weighted as (1-sumw)/1 since weights sum to 1 over the whole population
            vs = [(1-a.sumw) * var(a)/length(a) for a in deltas]
            sum(vs) > 0 || break
            nextSamples = allocate_samples(vs, min(round(Int, nsamples/3), samplesLeft), numSets .- totalSamples)
        else break end
    end

    # compute the Shapley values along with estimated variances of the estimates
    φ = zeros(length(featureGroups))
    φ[varyingInds] = [mean(a) for a in deltas]
    φVar = zeros(length(featureGroups))
    φVar[varyingInds] = [(1-a.sumw) * var(a)/length(a) for a in deltas]

    # find f(x) and E_x[f(x)]
    reset!(synthSampleSpace)
    addsample!(synthSampleSpace, x)
    fx = f(data(synthSampleSpace))[1]
    reset!(synthSampleSpace)
    for i in 1:N addsample!(synthSampleSpace, X[:,i]) end
    fnull = sum(f(data(synthSampleSpace)).*weights)

    # We ensure that the total of all features equals g(f(x)) - g(E_x[f(x)])
    trueSum = g(fx) - g(fnull)
    tmp = inv(φ*φ' + I*(trueSum*1e-8))
    β = inv(φ*φ' + I*(trueSum*1e-8))*φ*(sum(φ) - trueSum)
    #println("max β = ", maximum(β))
    φ .-= β.*φ
    @assert all(β .<= 1) "Rescaling failed! (indicates poor Shapley value estimates)"

    # return the Shapley values along with variances of the estimates
    fnull,φ,φVar
end

end # module


text:
using ShapeReconstructionPaLS.Utils
A = getDiamondModel([33;33;33]);
A = getBallModel([33;33;33]);

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
module TestExpressionsTupleForm

import Test: @testset, @test
import Volcanito: tuple_form, ColumnName
import MacroTools: postwalk, rmlines

@testset "tuple_form" begin
    input = tuple_form(
        :(x + y),
        (ColumnName(:x, false), ColumnName(:y, false)),
    )

    output = :(t -> $(Expr(:escape, :+))(t[1], t[2]))

    @test postwalk(rmlines, input) == postwalk(rmlines, output)

    input = tuple_form(
        :(y + x + $(Expr(:$, :z))),
        (ColumnName(:x, false), ColumnName(:y, false)),
    )

    output = :(t -> $(Expr(:escape, :+))(t[2], t[1], $(Expr(:escape, :z))))

    @test postwalk(rmlines, input) == postwalk(rmlines, output)
end

end


text:
include("linalg/old_linalg.jl")

# Calculate Ul, Dl, Tl =B(stop) ... B(start)
"""
Calculate effective(!) Green's function (direct, i.e. without stack) using QR DECOMPOSITION
"""
function calculate_slice_matrix_chain(mc::DQMC, start::Int, stop::Int, safe_mult::Int=mc.p.safe_mult)
  @assert 0 < start <= mc.p.slices
  @assert 0 < stop <= mc.p.slices
  @assert start <= stop

  flv = mc.model.flv
  N = mc.model.l.sites
  GreensType = geltype(mc)

  U = Matrix{GreensType}(I, flv*N, flv*N)
  D = ones(Float64, flv*N)
  T = Matrix{GreensType}(I, flv*N, flv*N)
  Tnew = Matrix{GreensType}(I, flv*N, flv*N)

  svs = zeros(flv*N,length(start:stop))
  svc = 1
  for k in start:stop
    if mod(k,safe_mult) == 0
      multiply_slice_matrix_left!(mc, mc.model, k, U)
      U *= spdiagm(0 => D)
      U, D, Tnew = decompose_udt(U)
      T =  Tnew * T
      svs[:,svc] = log.(D)
      svc += 1
    else
      multiply_slice_matrix_left!(mc, mc.model, k, U)
    end
  end
  return (U,D,T,svs)
end

# Calculate (Ur, Dr, Tr)' = B(stop) ... B(start)  => Ur,Dr, Tr = B(start)' ... B(stop)'
function calculate_slice_matrix_chain_dagger(mc::DQMC, start::Int, stop::Int, safe_mult::Int=mc.p.safe_mult)
  @assert 0 < start <= mc.p.slices
  @assert 0 < stop <= mc.p.slices
  @assert start <= stop

  flv = mc.model.flv
  N = mc.model.l.sites
  GreensType = geltype(mc)

  U = Matrix{GreensType}(I, flv*N, flv*N)
  D = ones(Float64, flv*N)
  T = Matrix{GreensType}(I, flv*N, flv*N)
  Tnew = Matrix{GreensType}(I, flv*N, flv*N)

  svs = zeros(flv*N,length(start:stop))
  svc = 1
  for k in reverse(start:stop)
    if mod(k,safe_mult) == 0
      multiply_daggered_slice_matrix_left!(mc, mc.model, k, U)
      U *= spdiagm(0 => D)
      U, D, Tnew = decompose_udt(U)
      T =  Tnew * T
      svs[:,svc] = log.(D)
      svc += 1
    else
      multiply_daggered_slice_matrix_left!(mc, mc.model, k, U)
    end
  end
  return (U,D,T,svs)
end

# Calculate G(slice) = [1+B(slice-1)...B(1)B(M) ... B(slice)]^(-1) and its singular values in a stable manner
function calculate_greens_and_logdet(mc::DQMC, slice::Int, safe_mult::Int=mc.p.safe_mult)
  GreensType = geltype(mc)
  flv = mc.model.flv
  N = mc.model.l.sites

  # Calculate Ur,Dr,Tr=B(slice)' ... B(M)'
  if slice <= mc.p.slices
    Ur, Dr, Tr = MonteCarlo.calculate_slice_matrix_chain_dagger(mc,slice,mc.p.slices, safe_mult)
  else
    Ur = Matrix{GreensType}(I, flv * N, flv * N)
    Dr = ones(Float64, flv * N)
    Tr = Matrix{GreensType}(I, flv * N, flv * N)
  end

  # Calculate Ul,Dl,Tl=B(slice-1) ... B(1)
  if slice-1 >= 1
    Ul, Dl, Tl = calculate_slice_matrix_chain(mc,1,slice-1, safe_mult)
  else
    Ul = eye(GreensType, flv * N)
    Dl = ones(Float64, flv * N)
    Tl = eye(GreensType, flv * N)
  end

  tmp = Tl * adjoint(Tr)
  U, D, T = decompose_udt(Diagonal(Dl) * tmp * Diagonal(Dr))
  U = Ul * U
  T *= adjoint(Ur)

  u, d, t = decompose_udt(adjoint(U) * inv(T) + Diagonal(D))

  T = inv(t * T)
  U *= u
  U = adjoint(U)
  d = 1. ./ d

  ldet = real(log(complex(det(U))) + sum(log.(d)) + log(complex(det(T))))

  return T * Diagonal(d) * U, ldet
end


text:
using Viznet
using Test

@testset "brush" begin
    include("brush.jl")
end

@testset "lattices" begin
    include("lattices/lattices.jl")
end

@testset "intersection" begin
    ring = [(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)]
    theta = π / 4.
    res = intersection(ring, theta, (-0.5, 0.0))
    @test all(res .≈ (0, 0.5))
end


text:
using Random

mutable struct Pixel{T}
    periods::Int
    val::T
    reliability::Float64
    groupsize::Int
    head::Pixel{T}
    last::Pixel{T}
    next::Union{Nothing, Pixel{T}}
    function Pixel{T}(periods, val, rel, gs) where T
        pixel = new(periods, val, rel, gs)
        pixel.head = pixel
        pixel.last = pixel
        pixel.next = nothing
        return pixel
    end
end
Pixel(v) = Pixel{typeof(v)}(0, v, rand(), 1)
@inline Base.length(p::Pixel) = p.head.groupsize

struct Edge{N}
    reliability::Float64
    pixel_1::CartesianIndex{N}
    pixel_2::CartesianIndex{N}
    periods::Int
end
function Edge{N}(pixel_image, ind1, ind2) where N
    @inbounds rel = pixel_image[ind1].reliability + pixel_image[ind2].reliability
    @inbounds periods = find_period(pixel_image[ind1].val, pixel_image[ind2].val)
    return Edge{N}(rel, ind1, ind2, periods)
end
Edge(g1, g2) = Edge(g1.reliability + g2.reliability,
                    g1, g2,
                    find_period(g1.val, g2.val))
Base.isless(e1::Edge, e2::Edge) = isless(e1.reliability, e2.reliability)

function unwrap!(wrapped_image::AbstractArray{T, N},
                 wrap_around::NTuple{N, Bool}=tuple(fill(false, N)...),
                 seed::Int=-1) where {T, N}

    if seed != -1
        Random.seed!(seed)
    end

    pixel_image = init_pixels(wrapped_image)
    calculate_reliability(pixel_image, wrap_around)
    edges = Edge{N}[]
    num_edges = _predict_num_edges(size(wrapped_image), wrap_around)
    sizehint!(edges, num_edges)
    for idx_dim=1:N
        populate_edges!(edges, pixel_image, idx_dim, wrap_around[idx_dim])
    end

    sort!(edges, alg=MergeSort)
    gather_pixels!(pixel_image, edges)
    unwrap_image!(wrapped_image, pixel_image)

    return wrapped_image
end

function _predict_num_edges(size_img, wrap_around)
    num_edges = 0
    for (size_dim, wrap_dim) in zip(size_img, wrap_around)
        num_edges += prod(size_img) * (size_dim-1) ÷ size_dim + wrap_dim * prod(size_img) ÷ size_dim
    end
    return num_edges
end

# function to broadcast
function init_pixels(wrapped_image)
    pixel_image = similar(wrapped_image, Pixel{eltype(wrapped_image)})
    @Threads.threads for i in eachindex(wrapped_image)
        @inbounds pixel_image[i] = Pixel(wrapped_image[i])
    end
    return pixel_image
end

function gather_pixels!(pixel_image, edges)
    for edge in edges
        @inbounds p1 = pixel_image[edge.pixel_1]
        @inbounds p2 = pixel_image[edge.pixel_2]
        merge_groups!(edge, p1, p2)
    end
end

function unwrap_image!(image, pixel_image)
    T = fieldtype(eltype(pixel_image), :val)
    this_pi = convert(T, π)
    @Threads.threads for i in eachindex(image)
        @inbounds image[i] = 2 * this_pi * pixel_image[i].periods + pixel_image[i].val
    end
end

function wrap_val(val)
    wrapped_val  = val
    wrapped_val += ifelse(val >  π, -2 * convert(typeof(val), π), zero(val))
    wrapped_val += ifelse(val < -π,  2 * convert(typeof(val), π), zero(val))
    return wrapped_val
end

function find_period(val_left, val_right)
    difference = val_left - val_right
    period  = 0
    period += ifelse(difference >  π, -1, 0)
    period += ifelse(difference < -π,  1, 0)
    return period
end

function merge_groups!(edge, pixel_1, pixel_2)
    if is_differentgroup(pixel_1, pixel_2)
        # pixel 2 is alone in group
        if is_pixelalone(pixel_2)
            merge_pixels!(pixel_1, pixel_2, -edge.periods)
        elseif is_pixelalone(pixel_1)
            merge_pixels!(pixel_2, pixel_1, edge.periods)
        else
            if is_bigger(pixel_1, pixel_2)
                merge_into_group!(pixel_1, pixel_2, -edge.periods)
            else
                merge_into_group!(pixel_2, pixel_1, edge.periods)
            end
        end
    end
end

@inline function is_differentgroup(p1::Pixel, p2::Pixel)
    return p1.head !== p2.head
end
@inline function is_pixelalone(pixel::Pixel)
    return pixel.head === pixel.last
end
@inline function is_bigger(p1::Pixel, p2::Pixel)
    return length(p1) ≥ length(p2)
end

function merge_pixels!(pixel_base::Pixel, pixel_target::Pixel, periods)
    pixel_base.head.groupsize += pixel_target.head.groupsize
    pixel_base.head.last.next = pixel_target.head
    pixel_base.head.last = pixel_target.head.last
    pixel_target.head = pixel_base.head
    pixel_target.periods = pixel_base.periods + periods
end

function merge_into_group!(pixel_base::Pixel, pixel_target::Pixel, periods)
    add_periods = pixel_base.periods + periods - pixel_target.periods
    pixel = pixel_target.head
    while pixel ≠ nothing
        # merge all pixels in pixel_target's group to pixel_base's group
        if pixel !== pixel_target
            pixel.periods += add_periods
            pixel.head = pixel_base.head
        end
        pixel = pixel.next
    end
    # assign pixel_target to pixel_base's group last
    merge_pixels!(pixel_base, pixel_target, periods)
end

function populate_edges!(edges, pixel_image::Array{T, N}, dim, connected) where {T, N}
    size_img       = collect(size(pixel_image))
    size_img[dim] -= 1
    idx_step       = fill(0, N)
    idx_step[dim] += 1
    idx_step_cart  = CartesianIndex{N}(idx_step...)
    idx_size       = CartesianIndex{N}(size_img...)
    for i in CartesianIndices((:).(Tuple(CartesianIndex{N}()),Tuple(idx_size)))
        push!(edges, Edge{N}(pixel_image, i, i+idx_step_cart))
    end
    if connected
        idx_step        = fill(0, N)
        idx_step[dim]   = -size_img[dim]
        idx_step_cart   = CartesianIndex{N}(idx_step...)
        edge_begin      = fill(1, N)
        edge_begin[dim] = size(pixel_image)[dim]
        edge_begin_cart = CartesianIndex{N}(edge_begin...)
        for i in CartesianIndices((:).(Tuple(edge_begin_cart), Tuple(CartesianIndex(size(pixel_image)))))
            push!(edges, Edge{N}(pixel_image, i, i+idx_step_cart))
        end
    end
end

function calculate_reliability(pixel_image::AbstractArray{T, N}, wrap_around) where {T, N}
    # get the shifted pixel indices in CartesinanIndex form
    # This gets all the nearest neighbors (CartesionIndex{N}() = one(CartesianIndex{N}))
    pixel_shifts = collect(CartesianIndices(tuple(repeat([-1:1],N)...)))
    size_img = size(pixel_image)
    # inner loop
    for i in CartesianIndices((:).(tuple(fill(2,length(size_img))...),size_img.-1))
        @inbounds pixel_image[i].reliability = calculate_pixel_reliability(pixel_image, i, pixel_shifts)
    end

    for (idx_dim, connected) in enumerate(wrap_around)
        if connected
            # first border
            pixel_shifts_border = copy(pixel_shifts)
            for (idx_ps, ps) in enumerate(pixel_shifts_border)
                # if the pixel shift goes out of bounds, we make the shift wrap
                if ps[idx_dim] == 1
                    new_ps          = fill(0, N)
                    new_ps[idx_dim] = -size_img[idx_dim]+1
                    pixel_shifts_border[idx_ps] = CartesianIndex{N}(new_ps...)
                end
            end
            border_begin          = fill(2, N)
            border_begin[idx_dim] = size_img[idx_dim]
            border_begin          = CartesianIndex{N}(border_begin...)
            border_end            = collect(size_img).-1
            border_end[idx_dim]   = size_img[idx_dim]
            border_end = CartesianIndex{N}(border_end...)
            for i in CartesianIndices((:).(Tuple(border_begin), Tuple(border_end)))
                @inbounds pixel_image[i].reliability = calculate_pixel_reliability(pixel_image, i, pixel_shifts_border)
            end
            # second border
            pixel_shifts_border = copyto!(pixel_shifts_border, pixel_shifts)
            for (idx_ps, ps) in enumerate(pixel_shifts_border)
                # if the pixel shift goes out of bounds, we make the shift wrap, this time to the other side
                if ps[idx_dim] == -1
                    new_ps          = fill(0, N)
                    new_ps[idx_dim] = size_img[idx_dim]-1
                    pixel_shifts_border[idx_ps] = CartesianIndex{N}(new_ps...)
                end
            end
            border_begin          = fill(2, N)
            border_begin[idx_dim] = 1
            border_begin          = CartesianIndex{N}(border_begin...)
            border_end            = collect(size_img).-1
            border_end[idx_dim]   = 1
            border_end = CartesianIndex{N}(border_end...)
            for i in CartesianIndices((:).(Tuple(border_begin), Tuple(border_end)))
                @inbounds pixel_image[i].reliability = calculate_pixel_reliability(pixel_image, i, pixel_shifts_border)
            end
        end
    end
end

function calculate_pixel_reliability(pixel_image::AbstractArray{T, N}, pixel_index, pixel_shifts) where {T, N}
    sum_val = zero(fieldtype(T, :val))
    for pixel_shift in pixel_shifts
        @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shift].val - pixel_image[pixel_index].val))^2
    end
    return sum_val
end

max_forks_repo_name: platawiec/Unwrap


text:
## Util.jl##

tupletypelength(a)=length(a.parameters)

macro replaceBase(fs...)
    b = Expr(:block)
    for f in fs
        if f == :/
            fa = :div_float
        elseif f == :^
            fa = :pow
        elseif f == :+
            fa = :vadd
        elseif f == :-
            fa = :vsub
        elseif f == :*
            fa = :vmul
        else
            fa = f
        end
        arg_consumed = false
        if fa in (:ceil,:floor,:sqrt,:rsqrt,:rec,
                  :exp,:exp2,:expm1,:log,:log1p,:log2,:log10,
                  :sin,:sinpi,:cos,:cospi,:tan,:tanpi,:asin,:acos,:atan,
                  :sinh,:cosh,:tanh,:asinh,:acosh,:atanh,
                  :trunc,:round,:exponent,:abs,:sincos,:cis)
            e = quote
                (Base.$f)(X::Array{T}) where {T <: Union{Float64,Float32}} = ($fa)(X)
                (Base.$f)(X::Union{Float64,Float32}) = ($fa)([X])[1]
            end
            arg_consumed = true
        end
        if fa in (:copysign,:atan,:pow,:rem,:div_float, :vadd, :vsub, :vmul)
            e = quote
                (Base.$f)(X::Array{T},Y::Array{T}) where {T <: Union{Float32,Float64}} = ($fa)(X,Y)
                (Base.$f)(X::T,Y::T) where {T <: Union{Float32,Float64}} = ($fa)([X],[Y])[1]
            end
            arg_consumed = true
        end
        if !arg_consumed
            error("Function $f not defined by AppleAccelerate.jl")
        end
        push!(b.args,e)
    end
    b
end


text:
# Create summary for Summary_value (default constructor by
# protobuf is broken).
const SummaryCollection(;kwargs...) = Summary(value=Base.Vector{Summary_Value}(); kwargs...)

function make_event(logger::Logger, summary::Summary;
                    step::Union{Nothing, Int}=nothing)
    # If the step is not set explicitly, get it from the logger
    if typeof(step) == Nothing
        step = iteration(logger)
    end
    Event(wall_time=time(), summary=summary, step=step)
end

function write_event(file::IOStream, event::Event)
    data = PipeBuffer();
    writeproto(data, event)

    #header
    header     = collect(reinterpret(UInt8, [data.size]))
    crc_header = reinterpret(UInt8, UInt32[masked_crc32c(header)])
    crc_data   = reinterpret(UInt8, UInt32[masked_crc32c(data.data)])

    write(file, header)
    write(file, crc_header)
    write(file, data.data)
    write(file, crc_data)
    flush(file)
end

write_event(logger::Logger, event::Event) = write_event(logger.file, event)


text:
function FESpace(
  model::DiscreteModel,
  cell_fe::CellFE;
  conformity=nothing,
  labels = get_face_labeling(model),
  dirichlet_tags=Int[],
  dirichlet_masks=nothing,
  constraint=nothing,
  vector_type::Union{Nothing,Type}=nothing)

  @assert num_cells(cell_fe) == num_cells(model) """\n
  The number of cells provided in the `cell_fe` argument ($(cell_fe.num_cells) cells)
  does not match the number of cells ($(num_cells(model)) cells) in the provided DiscreteModel.
  """

  trian = Triangulation(model)

  if vector_type == nothing
    cell_shapefuns, cell_dof_basis = compute_cell_space(cell_fe,trian)
    T = get_dof_value_type(cell_shapefuns,cell_dof_basis)
    _vector_type = Vector{T}
  else
    @assert vector_type <: AbstractVector """\n
    The (optional) argument vector_type has to be <: AbstractVector.
    """
    _vector_type = vector_type
  end

  if conformity in (L2Conformity(),:L2) && dirichlet_tags == Int[]
    F = _DiscontinuousFESpace(_vector_type,trian,cell_fe)
  else
    cell_conformity = CellConformity(cell_fe,conformity)
    F = _ConformingFESpace(
      _vector_type,
      model,
      labels,
      cell_fe,
      cell_conformity,
      dirichlet_tags,
      dirichlet_masks)
  end

  V = _add_constraint(F,cell_fe.max_order,constraint)

  V
end

function _add_constraint(F,order,constraint)
  if constraint == nothing
    V = F
  elseif constraint == :zeromean
    trian = get_triangulation(F)
    dΩ = Measure(trian,order)
    V = ZeroMeanFESpace(F,dΩ)
  else
    @unreachable """\n
    The passed option constraint=$constraint is not valid.
    Valid values for constraint: nothing, :zeromean
    """
  end
  V
end

function FESpace(
  model::RestrictedDiscreteModel, cell_fe::CellFE;constraint=nothing,kwargs...)
  model_portion = model.model
  V_portion = FESpace(model_portion,cell_fe;constraint=nothing,kwargs...)
  F = ExtendedFESpace(V_portion,model)
  V = _add_constraint(F,cell_fe.max_order,constraint)
  V
end

function FESpace(
  model::DiscreteModel,
  cell_reffe::AbstractArray{<:ReferenceFE};
  conformity=nothing,
  kwargs...)
  trian = Triangulation(model)
  cell_map = get_cell_map(trian)
  if conformity in (nothing, :default, :L2, L2Conformity())
    conf = conformity
  else
    conf = CellConformity(cell_reffe,conformity)
  end
  cell_fe = CellFE(cell_map,cell_reffe)
  FESpace(model,cell_fe;conformity=conf,kwargs...)
end

function FESpace(model::DiscreteModel, reffe::Tuple{<:ReferenceFEName,Any,Any}; kwargs...)
  basis, reffe_args,reffe_kwargs = reffe
  cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
  FESpace(model,cell_reffe;kwargs...)
end

function FESpace(model::DiscreteModel, reffe::ReferenceFE; kwargs...)
  cell_reffe = Fill(reffe,num_cells(model))
  FESpace(model,cell_reffe;kwargs...)
end

"""
"""
function TestFESpace(args...;kwargs...)
  FESpace(args...;kwargs...)
end


text:
# Autogenerated wrapper script for Libtiff_jll for aarch64-apple-darwin
export libtiff

using JpegTurbo_jll
using Zlib_jll
using Zstd_jll
JLLWrappers.@generate_wrapper_header("Libtiff")
JLLWrappers.@declare_library_product(libtiff, "@rpath/libtiff.5.dylib")
function __init__()
    JLLWrappers.@generate_init_header(JpegTurbo_jll, Zlib_jll, Zstd_jll)
    JLLWrappers.@init_library_product(
        libtiff,
        "lib/libtiff.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
Random.seed!(1)
a = mapweight(grid2(20),x->1/(rand(1)[1]));
la = lap(a)
n = size(la)[1]
b = randn(n)
b = b .- mean(b);

d = Array(diag(la))
pre = (x)->(x ./ d)

x = pcg(la,b,pre,maxits=2000,maxtime=1,verbose=true);
x = Laplacians.pcgBLAS(la,b,pre,maxits=2000,maxtime=1,verbose=true);
x = Laplacians.pcgSlow(la,b,pre,maxits=2000,maxtime=1,verbose=true);
x = pcg(la,b,pre,maxits=2000,verbose=true);
x = Laplacians.pcgBLAS(la,b,pre,maxits=2000,verbose=true);
x = Laplacians.pcgBLAS(la,b,pre,maxits=2000,verbose=true,maxtime=1e-10);
x = Laplacians.pcgSlow(la,b,pre,maxits=2000,verbose=true);
x = Laplacians.pcgSlow(la,b,pre,maxits=2000,verbose=true,maxtime=1e-10);
x = Laplacians.pcgBLAS(la,b,pre,tol=1e-2);
x = Laplacians.pcgSlow(la,b,pre,tol=1e-2);
x = pcg(la,b,pre,tol=1e-2);

x = cg(la,b,verbose=true,maxtime=1,maxits=1000);
x = Laplacians.cgBLAS(la,b,verbose=true,maxtime=1,maxits=1000);
x = Laplacians.cgBLAS(la,b,verbose=true,maxtime=1e-10,maxits=1000);
x = Laplacians.cgSlow(la,b,verbose=true,maxtime=1,maxits=1000);
x = Laplacians.cgSlow(la,b,verbose=true,maxtime=1e-10,maxits=1000);
x = cg(la,b,verbose=true,maxits=2000);
x = Laplacians.cgBLAS(la,b,verbose=true,maxits=2000);
x = Laplacians.cgSlow(la,b,verbose=true,maxits=2000);
x = cg(la,b,verbose=true,tol=0.5);
x = Laplacians.cgBLAS(la,b,verbose=true,tol=0.5);
x = Laplacians.cgSlow(la,b,verbose=true,tol=0.5);

bbig = convert(Array{BigFloat,1},b);
f = cgSolver(la)
x = f(bbig);

x = pcg(la,bbig,pre,maxits=200,maxtime=1,verbose=true);
x = cg(la,bbig,maxits=200,maxtime=1,verbose=true);

n = 100
a = wted_chimera(n,3)
la = lap(a)

f = cgSolver(la,verbose=false)
b = randn(n)
b = b .- mean(b)
x = f(b)
norm(la*x-b)

f(zeros(100))

x = f(b,maxits=15,verbose=true, maxtime=1)

d = Diagonal(1 ./ diag(la))
pre = x -> Array(d*x)

x = pcg(la, b, sparse(Diagonal(diag(la))),maxits=10)
x = pcg(la, b, pre, maxits=10, verbose=true)
f = pcgSolver(la, spdiagm(0=>diag(la)) ,maxits=10,verbose=true)
x = f(b)
f = pcgSolver(la,pre,maxits=10,verbose=false)
x = f(b,verbose=true, maxits=1000, maxtime = 2)

f(zeros(100))

t = akpw(a)
lt = lap(t)

f = pcgLapSolver(a, t)
pcgIts = [1]
x = f(b,verbose=true, pcgIts=pcgIts)
bbig = convert(Array{BigFloat,1},b);
x = f(bbig,verbose=true, pcgIts=pcgIts)

sdda = copy(la)
sdda[1,1] += 1
sddt = copy(lt)
sddt[1,1] += 1

f = pcgSolver(sdda, sddt)
x = f(b,verbose=true, pcgIts=pcgIts)
x = f(bbig,verbose=true, pcgIts=pcgIts)


text:
function qsat(water, P, T, cn::Constants)

    Tc = T - cn.Tm
    if (Tc > 0.0 || water)
        es = cn.e0 * exp(17.5043 * Tc / (241.3 + Tc))
    else
        es = cn.e0 * exp(22.4422 * Tc / (272.186 + Tc))
    end
    Qs = cn.eps * es / P

    return Qs

end


text:
module VarinnerBenchs

using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.Fields: MockField, MockBasis
using InplaceArrays.TensorValues
using FillArrays

@inline function loop(a,cache)
  for i in eachindex(a)
    ai = getindex!(cache,a,i)
  end
end

function bench1(n)
  p1 = Point(2,2)
  p2 = Point(4,2)
  p3 = Point(1,3)
  p4 = Point(5,2)
  x = [p1,p2,p3,p4]
  np = length(x)
  va = 3.0
  vb = 3.3
  d = 2
  fa = MockField{d}(va)
  fb = MockField{d}(vb)
  wa = 4.5
  wb = 4.5
  ndofa = 8
  ndofb = 5
  ba = MockBasis{d}(wa,ndofa)
  bb = MockBasis{d}(wb,ndofb)
  l = n
  ax = fill(x,l)
  afa = Fill(fa,l)
  afb = Fill(fb,l)
  aba = Fill(ba,l)
  abb = Fill(bb,l)

  ag_ff = varinner(afa,afb)
  agx_ff = evaluate(ag_ff,ax)
  cagx_ff = array_cache(agx_ff)
  @time loop(agx_ff,cagx_ff)

  ag_bf = varinner(aba,afb)
  agx_bf = evaluate(ag_bf,ax)
  cagx_bf = array_cache(agx_bf)
  @time loop(agx_bf,cagx_bf)

  ag_bb = varinner(aba,abb)
  agx_bb = evaluate(ag_bb,ax)
  cagx_bb = array_cache(agx_bb)
  @time loop(agx_bb,cagx_bb)

end

for n in (1,1,10,1000,100000)
  @eval begin
    println("+++ runing suite for n = $($n) +++")
    bench1($n)
  end
end

end # module


text:
import VTUFileHandler
import VTUFileHandler: set_uncompress_keywords, set_interpolation_keywords, VTUFile

set_uncompress_keywords(["xRamp"])
set_interpolation_keywords(["xRamp"])
vtu = VTUFile("vox8.vtu");
vtu += vtu/4;
vtu *= 4.0;
vtu -= 2.0;
vtu ^= 2.0;

write(vtu)

max_stars_repo_name: nagelt/VTUFileHandler

max_issues_repo_name: nagelt/VTUFileHandler

max_forks_repo_name: nagelt/VTUFileHandler


text:
# Base model struct
Base.@kwdef mutable struct guessing{T <: AbstractFloat} <: IRTmodel
    a::AbstractVector{T} = [0.0]
    d::AbstractVector{T} = [0.0]
    c::AbstractVector{T} = [0.0]
    group::Vector{Int64} = [1]
    fixed::Bool = false
    estc::Bool = true
end

function _distribute!(new, old::guessing)
    old.a = new[begin:length(old.a)]
    old.d = [new[length(old.a)+1]]
    old.c = [new[end]]
end

function _distribute!(new, old::guessing, q)
    old.a = new[begin:length(old.a)][q]
    old.d = [new[length(old.a)+1]]
    old.c = [new[end]]
end

function _distribute(new, old::guessing)
    out = copy(old)
    out.a = new[begin:length(out.a)]
    out.d = [new[length(out.a)+1]]
    out.c = [new[end]]
    return out
end

function _distribute(new, old::guessing, q)
    out = copy(old)
    out.a = new[begin:length(out.a)][q]
    out.d = [new[length(out.a)+1]]
    out.c = [new[end]]
    return out
end

function _take(par::guessing)
    if par.estc == true
        return [par.a; par.d; par.c]
    else
        return [par.a; par.d]
    end
end

function _take(par::guessing, q)
    if par.estc == true
        return [par.a[q]; par.d; par.c]
    else
        return [par.a[q]; par.d]
    end
end

function _copy(x::guessing)
    guessing(a = x.a, d = x.d, c = x.c, group = x.group, fixed = x.fixed, estc = x.estc)
end

function _rescale!(par::guessing, μ, Σ, q)
    if length(q) == 1
        # This item follows unidimensional model.
        A = 1/ sqrt(Σ[1]); K = -A*μ[1]
        par.a = par.a / A
        par.d = @. par.d - par.a*K/A
        par.c = [par.c[1]]
    else
        # This item seems to follow the multidimensional model.
        L = cholesky(Σ).L
        par.a = par.a * L # Dimension mismatch ?
        par.d = @. par.d + par.a * μ
        par.c = [par.c[1]]
    end
end

function _irf(p::guessing, θ, u, q)::Float64
    if u == 0
        return 1.0 - p.c[1] - (1.0 - p.c[1])*logistic(p.a[q]'θ[q] + p.d[1])
    else
        return p.c[1] + (1.0 - p.c[1])*logistic(p.a[q]'θ[q] + p.d[1])
    end
end

function _irf(p, m::guessing, θ, u, q)::eltype(p)
    a = p[1:length(m.a[q])]
    d = p[length(m.a[q])+1]
    c = m.estc ? p[end] : 0.0
    if u == 0
        return 1.0 - c[1] - (1.0 - c[1])*logistic(a[q]'θ[q] + d[1])
    else
        return c[1] + (1.0 - c[1])*logistic(a[q]'θ[q] + d[1])
    end
end