text:
import CUDA.cu
export cu,array

function cu(mpo::MPO{FloatType,Array{FloatType,4}}) where FloatType<:Number
    return MPO{FloatType,CuArray{FloatType,4}}(mpo.L,mpo.S,mpo.bdim,map(CuArray,mpo.tensor))
end

function cu(mps::MPS{FloatType,Array{FloatType,3}}) where FloatType<:Number
    return MPS{FloatType,CuArray{FloatType,3}}(mps.L,mps.S,mps.bdim,map(CuArray,mps.tensor))
end

function array(mpo::MPO{FloatType,CuArray{FloatType,4}}) where FloatType<:Number
    return MPO{FloatType,Array{FloatType,4}}(mpo.L,mpo.S,mpo.bdim,map(Array,mpo.tensor))
end
function array(mps::MPS{FloatType,CuArray{FloatType,3}}) where FloatType<:Number
    return MPS{FloatType,Array{FloatType,3}}(mps.L,mps.S,mps.bdim,map(Array,mps.tensor))
end

max_stars_repo_name: qiyang-ustc/JMPS

max_issues_repo_name: qiyang-ustc/JMPS

max_forks_repo_name: qiyang-ustc/JMPS


text:
using PyPlot

function ts(t, alpha)
  if t > 0
    a = ts(t-1, alpha)
    last_element = a -> alpha * a[end] + randn()
    return push!(a, last_element(a))
  else
    return [0.0]
  end
end

plot(ts(200, 0.98))
plot(ts(200, 0.8))
plot(ts(200, 0))
legend()

ext: jl


text:
@timed_testset "turing_model" begin
    DATA_DIR = joinpath("..", "data")
    kidiq = CSV.read(joinpath(DATA_DIR, "kidiq.csv"), DataFrame)
    wells = CSV.read(joinpath(DATA_DIR, "wells.csv"), DataFrame)
    roaches = CSV.read(joinpath(DATA_DIR, "roaches.csv"), DataFrame)
    cheese = CSV.read(joinpath(DATA_DIR, "cheese.csv"), DataFrame)
    @timed_testset "Gaussian Model" begin
        f = @formula(kid_score ~ mom_iq * mom_hs)
        @testset "standardize=false" begin
            m = turing_model(f, kidiq)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 29.30 atol = 2.0
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.533 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.593 atol = 0.2
        end

        @testset "standardize=true" begin
            m = turing_model(f, kidiq; standardize=true)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 0.000 atol = 0.2
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.648 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.849 atol = 0.2
        end

        @testset "custom_priors" begin
            priors = CustomPrior(Normal(), Normal(28, 5), nothing)
            m = turing_model(f, kidiq; priors)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 28.758 atol = 2.0
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.539 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.3863 atol = 0.2
        end
        @testset "explicit calling Normal" begin
            m = turing_model(f, kidiq; model=Normal)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 29.30 atol = 2.0
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.533 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.593 atol = 0.2
        end
    end
    @timed_testset "TDist Model" begin
        f = @formula(kid_score ~ mom_iq * mom_hs)
        @testset "standardize=false" begin
            m = turing_model(f, kidiq; model=TDist)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 40.380 atol = 2.0
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.478 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.736 atol = 0.2
            @test quantile(chn)[:ν, Symbol("50.0%")] ≈ 1.039 atol = 0.5
        end

        @testset "custom_priors" begin
            priors = CustomPrior(Normal(), Normal(28, 5), Exponential(2))
            m = turing_model(f, kidiq; model=TDist, priors)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 35.506 atol = 2.0
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.522 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.628 atol = 0.2
            @test quantile(chn)[:ν, Symbol("50.0%")] ≈ 1.178 atol = 0.5
        end
    end
    @timed_testset "Bernoulli Model" begin
        f = @formula(switch ~ arsenic + dist + assoc + educ)
        @testset "standardize=false" begin
            m = turing_model(f, wells; model=Bernoulli)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ -0.153 atol = 0.2
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.467 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.009 atol = 0.2
        end

        @testset "custom_priors" begin
            priors = CustomPrior(Normal(), Normal(), nothing)
            m = turing_model(f, wells; model=Bernoulli, priors)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ -0.155 atol = 0.2
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.468 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.009 atol = 0.2
        end
    end
    @timed_testset "Poisson Model" begin
        f = @formula(y ~ roach1 + treatment + senior + exposure2)
        @testset "standardize=false" begin
            m = turing_model(f, roaches; model=Poisson)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 2.969 atol = 0.5
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.006 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.5145 atol = 0.2
        end

        @testset "custom_priors" begin
            priors = CustomPrior(Normal(2, 5), Normal(), nothing)
            m = turing_model(f, roaches; model=Poisson, priors)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 2.963 atol = 0.5
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.006 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.5145 atol = 0.2
        end
    end
    @timed_testset "NegativeBinomial Model" begin
        f = @formula(y ~ roach1 + treatment + senior + exposure2)
        @testset "standardize=false" begin
            m = turing_model(f, roaches; model=NegativeBinomial)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 2.448 atol = 0.5
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.013 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.734 atol = 0.2
            @test quantile(chn)[:ϕ⁻, Symbol("50.0%")] ≈ 1.401 atol = 0.2
        end

        @testset "custom_priors" begin
            priors = CustomPrior(Normal(), Normal(2, 5), Exponential(0.5))
            m = turing_model(f, roaches; model=NegativeBinomial, priors)
            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
            @test summarystats(chn)[:α, :mean] ≈ 2.422 atol = 0.5
            @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.013 atol = 0.2
            @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.732 atol = 0.2
            @test quantile(chn)[:ϕ⁻, Symbol("50.0%")] ≈ 3.56 atol = 0.2
        end
    end
    @timed_testset "Hierarchical Model" begin
        f = @formula(y ~ (1 | cheese) + background)
        m = turing_model(f, cheese)
        chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
        @test summarystats(chn)[:α, :mean] ≈ 68.33 atol = 2.0
        @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 6.928 atol = 0.2
        @test summarystats(chn)[Symbol("zⱼ[1]"), :mean] ≈ 0.306 atol = 0.2
        @test quantile(chn)[Symbol("zⱼ[2]"), Symbol("50.0%")] ≈ -1.422 atol = 0.5
    end
    @testset "Unsupported Model Likelihoods" begin
        @test_throws ArgumentError turing_model(@formula(y ~ x), nt_str; model=Binomial)
    end
    @testset "NegativeBinomial2" begin
        @test T.NegativeBinomial2(0, 1) == T.NegativeBinomial(1, 1)
        @test T.NegativeBinomial2(2, 8) == T.NegativeBinomial(8, 0.8)
        @test_throws DomainError T.NegativeBinomial2(0, -1)
    end
end


text:
using Test
using MaBLAS
using LinearAlgebra
using TimerOutputs
using VectorizationBase: pick_vector_width

noreturn_mul!(args...; kwargs...) = (MaBLAS.mul!(args...; kwargs...); nothing)

@testset "Size check" begin
    C = randn(12, 11); A = rand(2, 3); B = rand(2, 3)
    @test_throws DimensionMismatch MaBLAS.mul!(C, A, B)
end

@testset "Kernel sized matmul" begin
    _m, _k, _n = 8*4*4, 8, 6*4*4*3
    # lower cache_params to check multiple of cache sizes more easily
    cache_params = (_m, _k, _n)
    m, k, n = _m*2, _k*3, _n*5
    C = rand(m, n)
    A = rand(m, k)
    B = rand(k, n)
    T = Float64
    N = pick_vector_width(T)
    cache_params = cache_params .* 4N
    for α in (1, 2, 3, false, true), β in (1, 2, 3, false, true), packa in (true, false), packb in (true, false)
        for kernel_params in [(Val(N), Val(6)), (Val(N), Val(4)), (Val(2N), Val(3))]
            @test MaBLAS.mul!((copy(C)), A, B, α, β; cache_params=cache_params, packing=(packa, packb), kernel_params=kernel_params) ≈ LinearAlgebra.mul!((copy(C)), A, B, α, β)
        end
    end
    packa, packb = true, false
    α, β = randn(2)
    @test MaBLAS.mul!((copy(C)), A, B, α, β; cache_params=cache_params, packing=(packa, packb)) == MaBLAS.mul!((copy(C)), A, B, α, β; cache_params=cache_params, packing=(Val(packa), Val(packb)))
end

@testset "Block size tests" begin
    @test MaBLAS.partition_k(2000, 532) == 500
    @test MaBLAS.partition_k(401, 400) == 201
    @test MaBLAS.partition_k(900, 400) == 300
end

@testset "Clean up loop tests" begin
    _m, _k, _n = 8*4*4, 8*2, 6*4*3
    # lower cache_params to check clean up loops more easily
    cache_params = (_m, _k, _n)
    α, β = randn(2)
    for T in (Float64, Float32)
        N = pick_vector_width(T)
        for packa in (true, false), packb in (true, false), (m, k, n) in [(73, 131, 257), (101, 103, 107), (109, 113, 127), (131, 137, 139), (149, 151, 157), (163, 167, 173), (179, 181, 191), (193, 197, 199)] # all primes
            for kernel_params in [(Val(N), Val(6)), (Val(N), Val(4)), (Val(2N), Val(3))]
                C = rand(T, m, n)
                A = rand(T, m, k)
                B = rand(T, k, n)
                @test MaBLAS.mul!((copy(C)), A, B, α, β; cache_params=cache_params, packing=(packa, packb), kernel_params=kernel_params) ≈ LinearAlgebra.mul!((copy(C)), A, B, α, β)
                @allocated noreturn_mul!((copy(C)), A, B, α, β; cache_params=cache_params, packing=(packa, packb), kernel_params=kernel_params)
                @allocated(noreturn_mul!((copy(C)), A, B, α, β; cache_params=cache_params, packing=(packa, packb), kernel_params=kernel_params)) <= (packa + packb) * 176 # minor allocation when packing
            end
        end
    end
end

@testset "View of matrix" begin
    V = randn(200, 200)
    A = @view V[1:100,  1:100]
    B = @view V[101:end, 1:100]
    C = @view V[1:100, 101:end]
    _m, _k, _n = 8*4*4, 8*2, 6*4*3
    cache_params = (cache_m=_m, cache_k=_k, cache_n=_n)
    @test LinearAlgebra.mul!((copy(C)), A, B) ≈ MaBLAS.mul!(C, A, B; cache_params=cache_params)
    for packa in (true, false), packb in (true, false)
        @test LinearAlgebra.mul!((copy(C)), A, B) ≈ MaBLAS.mul!(C, A, B; cache_params=cache_params, packing=(packa, packb))
    end

    A = @view V[1:2:end, 1:100]
    @test_throws ArgumentError MaBLAS.mul!(C, A, B)
    @test_throws ArgumentError MaBLAS.mul!(C, A, B; packing=(false, true))
    for packb in (true, false)
        @test LinearAlgebra.mul!((copy(C)), A, B) ≈ MaBLAS.mul!(C, A, B; cache_params=cache_params, packing=(true, packb))
    end

    for packa in (true, false), packb in (true, false)
        @test LinearAlgebra.mul!((copy(C)), B, A) ≈ MaBLAS.mul!(C, B, A; cache_params=cache_params, packing=(packa, packb))
    end

    C = @view V[1:2:end, 101:end]
    for packa in (true, false), packb in (true, false)
        @test_throws ArgumentError MaBLAS.mul!(C, A, B; cache_params=cache_params, packing=(packa, packb))
    end
end

@testset "Timer & Transposes/Adjoints" begin
    m, k, n = 73, 131, 257 # all prime numbers
    C = rand(m, n)
    timer = MaBLAS.get_timer()
    MaBLAS.enable_timer()
    MaBLAS.reset_timer!()
    for α in (1, 2, 3.0, false, true), β in (1, 2, 3.0, false, true), A in (rand(m, k), rand(k, m)'), B in (rand(k, n), rand(n, k)'), packa in (true, false), packb in (true, false)
        !packa && A isa Adjoint && continue
        AB = MaBLAS.mul!((copy(C)), A, B, α, β; packing=(packa, packb))
        @test AB ≈ LinearAlgebra.mul!((copy(C)), A, B, α, β)
    end
    tim, alloc = TimerOutputs.totmeasured(timer)
    @test tim > 0
    @test alloc == 0
    display(timer)

    MaBLAS.reset_timer!()
    MaBLAS.disable_timer()

    α, β = rand(2)
    for A in (rand(m, k), rand(k, m)', PermutedDimsArray(rand(k, m), (2, 1))), B in (rand(k, n), rand(n, k)', PermutedDimsArray(rand(n, k), (2, 1)))
        MaBLAS.mul!((copy(C)), A, B, α, β; packing=(true, true)) ≈ LinearAlgebra.mul!((copy(C)), A, B, α, β)
    end
    display(timer)
    @test all(iszero, TimerOutputs.totmeasured(timer))
end


text:
@inline function Phi_e(CellData,s,z)
   """
   Electrolyte Potential Transfer Function
   # Add Ins and Outs
       # Cell Data
       # Frequency Vector
       # Discretisation Locations
       # Electrode Definition
   """

CC_A = CellData.Const.CC_A   # Current-collector area [m^2]
De = CellData.Const.De # Electrolyte Diffusivity
κ_eff_Neg = CellData.Const.κ*CellData.Neg.ϵ_e^CellData.Neg.κ_brug
κ_eff_Sep = CellData.Const.κ*CellData.Sep.ϵ_e^CellData.Sep.κ_brug
κ_eff_Pos = CellData.Const.κ*CellData.Pos.ϵ_e^CellData.Pos.κ_brug
σ_eff_Neg = CellData.Neg.σ*CellData.Neg.ϵ_s^CellData.Neg.σ_brug #Effective Conductivity Neg
σ_eff_Pos = CellData.Pos.σ*CellData.Pos.ϵ_s^CellData.Pos.σ_brug #Effective Conductivity Pos
#dln = CellData.Const.dln  #Electrolyte activity coefficient term (Rod. 17)
#κ_D_eff = (2*R*CellData.Const.T/F)*κ_eff_Neg*(1-CellData.Const.t_plus)*(1+dln) #Diffision Effective Electrolyte Conductivity

#Defining SOC
θ_neg = CellData.Const.SOC * (CellData.Neg.θ_100-CellData.Neg.θ_0) + CellData.Neg.θ_0
θ_pos = CellData.Const.SOC * (CellData.Pos.θ_100-CellData.Pos.θ_0) + CellData.Pos.θ_0

#Beta's
βn = @. CellData.Neg.Rs*sqrt(s/CellData.Neg.Ds)
βp = @. CellData.Pos.Rs*sqrt(s/CellData.Pos.Ds)

#Prepare for j0
ce0 = CellData.Const.ce0
cs_max_neg = CellData.Neg.cs_max
cs0_neg = cs_max_neg * θ_neg
cs_max_pos = CellData.Pos.cs_max
cs0_pos = cs_max_pos * θ_pos
α_neg = CellData.Neg.α
α_pos = CellData.Pos.α

#Current Flux Density
if CellData.Const.CellTyp == "Doyle_94"
   κ_neg = CellData.Neg.k_norm/CellData.Neg.cs_max/ce0^(1-α_neg)
   κ_pos = CellData.Pos.k_norm/CellData.Pos.cs_max/ce0^(1-α_pos)
   j0_neg = κ_neg*(ce0*(cs_max_neg-cs0_neg))^(1-α_neg)*cs0_neg^α_neg
   j0_pos = κ_pos*(ce0*(cs_max_pos-cs0_pos))^(1-α_pos)*cs0_pos^α_pos
else
   j0_neg = CellData.Neg.k_norm*(ce0*(cs0_neg/cs_max_neg*(1-cs0_neg/cs_max_neg)))^(1-CellData.Neg.α)
   j0_pos = CellData.Pos.k_norm*(ce0*(cs0_pos/cs_max_pos*(1-cs0_pos/cs_max_pos)))^(1-CellData.Pos.α)
end

#Resistances
Rtot_neg = R*CellData.Const.T/(j0_neg*F^2) + CellData.Neg.RFilm
Rtot_pos = R*CellData.Const.T/(j0_pos*F^2) + CellData.Pos.RFilm

#OCP derivative
∂Uocp_neg = CellData.Const.∂Uocp("Neg",θ_neg)/cs_max_neg
∂Uocp_pos = CellData.Const.∂Uocp("Pos",θ_pos)/cs_max_pos

ν_neg = @. CellData.Neg.L*sqrt((CellData.Neg.as/σ_eff_Neg+CellData.Neg.as/κ_eff_Neg)/(Rtot_neg+∂Uocp_neg*(CellData.Neg.Rs/(F*CellData.Neg.Ds))*(tanh(βn)/(tanh(βn)-βn)))) #Condensing Variable - eq. 4.13
ν_neg_∞ = @. CellData.Neg.L*sqrt(CellData.Neg.as*((1/κ_eff_Neg)+(1/σ_eff_Neg))/(Rtot_neg))
ν_pos = @. CellData.Pos.L*sqrt((CellData.Pos.as/σ_eff_Pos+CellData.Pos.as/κ_eff_Pos)/(Rtot_pos+∂Uocp_pos*(CellData.Pos.Rs/(F*CellData.Pos.Ds ))*(tanh(βp)/(tanh(βp)-βp)))) #Condensing Variable - eq. 4.13
ν_pos_∞ = @. CellData.Pos.L*sqrt(CellData.Pos.as*((1/κ_eff_Pos)+(1/σ_eff_Pos))/(Rtot_pos))

# if tf = Rod17 #Consider splitting into different function calls.

#    for i < length(z)
#    pt = z[i]

#       if pt <= CellData.Neg.L
#          ϕ_tf = @. CellData.Neg.L*(σ_eff*(1-cosh(ν_neg*pt)/(as*κ_eff*(κ_eff+σ_eff)*sinh(ν_neg)*ν_neg) + (CellData.Neg.L*(κ_eff*cosh(ν_neg)-cosh(ν_neg*(1-pt))-pt*sinh(ν_neg)*ν_neg))/(A*κ_eff*(κ_eff+σ_eff)*sinh(ν_neg)*ν_neg)
#          zero_tf = @. -(CellData.Neg.L*pt^2)/(2*as*κ_eff) - (κ_d*CellData.Neg.L)(t_plus-1)*pt^2)/(2*A*ce0*De*F*κ_eff)
#          #D  =  @. CellData.Neg.L*(σ_eff*(1-cosh(ν_∞*pt)/(as*κ_eff*(κ_eff+σ_eff)*sinh(ν_∞)*ν_∞) + (CellData.Neg.L*(κ_eff*cosh(ν_∞)-cosh(ν_∞*(1-pt))-pt*sinh(ν_∞)*ν_∞))/(A*κ_eff*(κ_eff+σ_eff)*sinh(ν_∞)*ν_∞)
#          ϕ_tf[findall(s.==0),:] .= zero_tf[findall(s.==0),:]

#       elseif pt <= CellData.Neg.L + CellData.Sep.L
#          ϕ_tf1 = @. -pt*CellData.Sep.L/(as*κ_eff)
#          ϕ_tf2 = @. -(κ_D_eff/κ_eff_Pos*ce0)*(cS1*(e^ΛS1*pt-1)+cS2*(e^ΛS1*pt-1))
#          ϕ_n_tf = @. CellData.Neg.L*(σ_eff*(1-cosh(ν_neg*CellData.Neg.L)/(as*κ_eff*(κ_eff+σ_eff)*sinh(ν_neg)*ν_neg) + (CellData.Neg.L*(κ_eff*cosh(ν_neg)-cosh(ν_neg*(1-CellData.Neg.L))-CellData.Neg.L*sinh(ν_neg)*ν_neg))/(A*κ_eff*(κ_eff+σ_eff)*sinh(ν_neg)*ν_neg)
#          ϕ_tf = @. ϕ_n_tf + ϕ_tf1 + ϕ_tf2
#          zero_tf = @. -(CellData.Neg.L*pt^2)/(2*as*κ_eff) - (κ_d*CellData.Neg.L)(t_plus-1)*pt^2)/(2*A*ce0*De*F*κ_eff)
#          #D  = @. CellData.Sep.L*(σ_eff*(1-cosh(ν_∞*pt)/(as*κ_eff*(κ_eff+σ_eff)*sinh(ν_∞)*ν_∞) + (CellData.Sep.L*(κ_eff*cosh(ν_∞)-cosh(ν_∞*(1-pt))-pt*sinh(ν_∞)*ν_∞))/(A*κ_eff*(κ_eff+σ_eff)*sinh(ν_∞)*ν_∞) - pt*CellData.Sep.L/(as*κ_eff)
#          ϕ_tf[findall(s.==0),:] .= zero_tf[findall(s.==0),:]

#       else
#          ϕ_s_tf1 = @. -pt*CellData.Sep.L/(as*κ_eff)
#          ϕ_s_tf2 = @. -(κ_D_eff/κ_eff_Pos*ce0)*(cS1*(e^ΛS1*pt-1)+cS2*(e^ΛS1*pt-1))
#          ϕ_n_tf = @. CellData.Neg.L*(σ_eff*(1-cosh(ν_neg*CellData.Neg.L)/(as*κ_eff*(κ_eff+σ_eff)*sinh(ν_neg)*ν_neg) + (CellData.Neg.L*(κ_eff*cosh(ν_neg)-cosh(ν_neg*(1-CellData.Neg.L))-CellData.Neg.L*sinh(ν_neg)*ν_neg))/(A*κ_eff*(κ_eff+σ_eff)*sinh(ν_neg)*ν_neg)
#          ϕ_p_tf1 = ((pt-1)*CellData.Pos.L/κ_eff_Pos*CC_A)+(CellData.Pos.as*F*CellData.Pos.L^2/κ_eff_Pos) * (j_Pos_1*(e^(ΛP1*pt)+(1-pt)*ΛP1*e^ΛP1-e^ΛP1)/ΛP1^2 + j_Pos_2*(e^(-ΛP1*pt)+(pt-1)*ΛP1*e^-ΛP1-e^-ΛP1)/ΛP1^2 + j_Pos_1*(e^(ΛP1*pt)+(1-pt)*ΛP1*e^ΛP1-e^ΛP1)/ΛP1^2 + j_Pos_3*(e^(ΛP2*pt)+(1-pt)*ΛP2*e^ΛP2-e^ΛP2)/ΛP2^2 + j_Pos_4*(e^(-ΛP2*pt)+(pt-1)*ΛP2*e^-ΛP2-e^-ΛP2)/ΛP2^2)
#          ϕ_p_tf2 = -(κ_D_eff/κ_eff_Pos*ce0)*(cp1*(e^ΛP1*pt-e^Λp1)+cp2*(e^ΛP1*pt-e^Λp1)+cp3*(e^ΛP2*pt-e^Λp2)+cp4*(e^ΛP2*pt-e^Λp2))
#          ϕ_tf = @. ϕ_n_tf + ϕ_s_tf1 + ϕ_s_tf2 + ϕ_p_tf1 + ϕ_p_tf2

#          #D  = @.
#    end

ϕ_tf = Array{ComplexF64}(undef,length(z),length(s))
D_term = Array{String}(undef,length(z),1)
D = fill(0.0,length(z))
i=Int64(1)
# Loop Tf's
   for pt in z
         if pt <= CellData.Neg.L+eps()
            ϕ_tf[i,:] = @. (CellData.Neg.L*(σ_eff_Neg/κ_eff_Neg)*(1-cosh(ν_neg*pt/CellData.Neg.L)) - pt*ν_neg*sinh(ν_neg) + CellData.Neg.L*(cosh(ν_neg)-cosh(ν_neg*(CellData.Neg.L-pt)/CellData.Neg.L)))/(CC_A*(κ_eff_Neg+σ_eff_Neg)*sinh(ν_neg)*ν_neg) #Lee. Eqn. 4.22
            zero_tf = @. -(pt^2)/(2*CC_A*κ_eff_Neg*CellData.Neg.L)
            D[i]  =  @. (CellData.Neg.L*(σ_eff_Neg/κ_eff_Neg)*(1-cosh(ν_neg_∞*pt/CellData.Neg.L)) - pt*ν_neg_∞*sinh(ν_neg_∞) + CellData.Neg.L*(cosh(ν_neg_∞)-cosh(ν_neg_∞*(CellData.Neg.L-pt)/CellData.Neg.L)))/(CC_A*(κ_eff_Neg+σ_eff_Neg)*sinh(ν_neg_∞)*ν_neg_∞) #Lee. Eqn. 4.22 @ ∞
            ϕ_tf[i,findall(s.==0)] .= zero_tf
            D_term[i] = "@. ($(CellData.Neg.L)*($σ_eff_Neg/$κ_eff_Neg)*(1-cosh(ν_neg*$pt/$(CellData.Neg.L))) - $pt*ν_neg*sinh(ν_neg))/($CC_A*($κ_eff_Neg+$σ_eff_Neg)*sinh(ν_neg)*ν_neg) + ($(CellData.Neg.L)*(cosh(ν_neg)-cosh(ν_neg*($(CellData.Neg.L)-$pt)/$(CellData.Neg.L))/($CC_A*$κ_eff_Neg*($κ_eff_Neg+$σ_eff_Neg)*sinh(ν_neg)*ν_neg)))" #Lee. Eqn. 4.22 @ ∞

         elseif pt <= CellData.Neg.L + CellData.Sep.L + eps()
            ϕ_tf[i,:] = @. (CellData.Neg.L - pt)/(CC_A*κ_eff_Sep) + (CellData.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg/2)-ν_neg))/(CC_A*(κ_eff_Neg+σ_eff_Neg)*ν_neg) #Lee. Eqn. 4.23
            zero_tf = @. (2*κ_eff_Neg*CellData.Neg.L-κ_eff_Sep*CellData.Neg.L-2*κ_eff_Neg*pt)/(2*CC_A*κ_eff_Neg*κ_eff_Sep)
            D[i] = @. (CellData.Neg.L - pt)/(CC_A*κ_eff_Sep) + (CellData.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg_∞/2)-ν_neg_∞))/(CC_A*(κ_eff_Neg+σ_eff_Neg)*ν_neg_∞) #Lee. Eqn. 4.23 @ ∞
            ϕ_tf[i,findall(s.==0)] .= zero_tf
            D_term[i] = "@. ($(CellData.Neg.L - pt))/($CC_A*$κ_eff_Sep) + ($(CellData.Neg.L)*((1-$σ_eff_Neg/$κ_eff_Neg)*tanh(ν_neg/2)-ν_neg))/($CC_A*($κ_eff_Neg+$σ_eff_Neg)*ν_neg)" #Lee. Eqn. 4.23 @ ∞

         else
            #ϕ_tf[i,:] = @. -CellData.Sep.L/(CC_A*κ_eff_Sep) + (CellData.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg/2)-ν_neg))/(CC_A*(κ_eff_Neg+σ_eff_Neg)*ν_neg) - (CellData.Pos.L*(1+(σ_eff_Pos/κ_eff_Pos)*cosh(ν_pos)) + CellData.Pos.L*cosh((CellData.Neg.L+CellData.Sep.L-pt)*ν_pos/CellData.Pos.L) + CellData.Pos.L*((σ_eff_Pos/κ_eff_Pos)*cosh((Ltot-pt)*ν_pos/CellData.Pos.L)))/(CC_A*(κ_eff_Pos+σ_eff_Pos)*sinh(ν_pos)*ν_pos) + (CellData.Neg.L + CellData.Sep.L - pt)/(CC_A*((κ_eff_Pos+σ_eff_Pos)))  #Lee. Eqn. 4.24
            #ϕ_tf[i,:] = @. (-CellData.Neg.L*((σ_eff_Neg-κ_eff_Neg)*tanh(ν_neg/2)))/(CC_A*(κ_eff_Neg+σ_eff_Neg)*ν_neg) - CellData.Neg.L/(CC_A*(κ_eff_Neg+σ_eff_Neg)) - (CellData.Sep.L/(CC_A*κ_eff_Sep)) - CellData.Pos.L*(1-cosh((pt-1)*ν_pos))/((CC_A*(κ_eff_Pos+σ_eff_Pos)*sinh(ν_pos)*ν_pos)) - CellData.Pos.L*σ_eff_Pos*(cosh(ν_pos)-cosh(pt*ν_pos))/((CC_A*κ_eff_Pos*(κ_eff_Pos+σ_eff_Pos)*sinh(ν_pos)*ν_pos)) - CellData.Pos.L*(pt-1)/(CC_A*(κ_eff_Pos+σ_eff_Pos))
            ϕ_tf[i,:] = @. -CellData.Sep.L/(CC_A*κ_eff_Sep) + CellData.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg/2)-ν_neg)/(CC_A*(σ_eff_Neg+ κ_eff_Neg)*ν_neg) + (CellData.Pos.L*(-σ_eff_Pos*cosh(ν_pos) + σ_eff_Pos*cosh((CellData.Const.Ltot-pt)*ν_pos/CellData.Pos.L) +  κ_eff_Pos*(cosh((pt-CellData.Neg.L-CellData.Sep.L)*ν_pos/CellData.Pos.L)-1)) - (pt-CellData.Neg.L-CellData.Sep.L)*κ_eff_Pos*sinh(ν_pos)*ν_pos)/(CC_A*κ_eff_Pos*(κ_eff_Pos+σ_eff_Pos)*ν_pos*sinh(ν_pos))
            #zero_tf = @. -CellData.Neg.L/(2*CC_A*κ_eff_Neg) - CellData.Sep.L/(2*CC_A*κ_eff_Sep) + (CellData.Neg.L+CellData.Sep.L-pt)*(CellData.Neg.L+2*CellData.Pos.L+CellData.Sep.L-pt)/(2*CC_A*κ_eff_Pos*CellData.Pos.L)
            zero_tf = @. - CellData.Neg.L/(2*CC_A*κ_eff_Neg) - CellData.Sep.L/(CC_A*κ_eff_Sep) - (CellData.Pos.L-(CellData.Const.Ltot-pt)^2/CellData.Pos.L)/(2*CC_A*κ_eff_Pos)
            D[i] = @. -CellData.Sep.L/(CC_A*κ_eff_Sep) + CellData.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg_∞/2)-ν_neg_∞)/(CC_A*(σ_eff_Neg+ κ_eff_Neg)*ν_neg_∞) + (CellData.Pos.L*(-σ_eff_Pos*cosh(ν_pos_∞) + σ_eff_Pos*cosh((CellData.Const.Ltot-pt)*ν_pos_∞/CellData.Pos.L) +  κ_eff_Pos*(cosh((pt-CellData.Neg.L-CellData.Sep.L)*ν_pos_∞/CellData.Pos.L)-1)) - (pt-CellData.Neg.L-CellData.Sep.L)*κ_eff_Pos*sinh(ν_pos_∞)*ν_pos_∞)/(CC_A*κ_eff_Pos*(κ_eff_Pos+σ_eff_Pos)*ν_pos_∞*sinh(ν_pos_∞))
            #D[i]  = @. -CellData.Sep.L/(CC_A*κ_eff_Sep) + (CellData.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg_∞/2)-ν_neg_∞))/(CC_A*(κ_eff_Neg+σ_eff_Neg)*ν_neg_∞) - (CellData.Pos.L*(1+(σ_eff_Pos/κ_eff_Pos)*cosh(ν_pos_∞)) + CellData.Pos.L*cosh((CellData.Neg.L+CellData.Sep.L-pt)*ν_pos_∞/CellData.Pos.L) + CellData.Pos.L*((σ_eff_Pos/κ_eff_Pos)*cosh((Ltot-pt)*ν_pos_∞/CellData.Pos.L)))/(CC_A*(κ_eff_Pos+σ_eff_Pos)*sinh(ν_pos_∞)*ν_pos_∞) + (CellData.Neg.L + CellData.Sep.L - pt)/(CC_A*((κ_eff_Pos+σ_eff_Pos))) #Lee. Eqn. 4.24 @ ∞
            ϕ_tf[i,findall(s.==0)] .= zero_tf
            D_term[i] = "@. -$(CellData.Sep.L)/($CC_A*$κ_eff_Sep) + ($(CellData.Neg.L)*((1-$σ_eff_Neg/$κ_eff_Neg)*tanh(ν_neg/2)-ν_neg))/($CC_A*($κ_eff_Neg+$σ_eff_Neg)*ν_neg) - ($(CellData.Pos.L)*(1+($σ_eff_Pos/$κ_eff_Pos)*cosh(ν_pos)) + $(CellData.Pos.L)*cosh(($(CellData.Neg.L)+$(CellData.Sep.L)-$pt)*ν_pos/$(CellData.Pos.L)) + $(CellData.Pos.L)*(($σ_eff_Pos/$κ_eff_Pos)*cosh(($CellData.Const.Ltot-$pt)*ν_pos/$(CellData.Pos.L))))/($CC_A*($κ_eff_Pos+$σ_eff_Pos)*sinh(ν_pos)*ν_pos) + ($(CellData.Neg.L) + $(CellData.Sep.L) - $pt)/($CC_A*(($κ_eff_Pos+$σ_eff_Pos)))" #Lee. Eqn. 4.24 @ ∞
         end
      i=i+1
 end
 if Debug == 1
      println("D:Phi_e:",size(D))
      println("ϕ_tf:Phi_e:",ϕ_tf[:,2])
      println("D:Phi_e:Pos",D)
      println("z:Phi_e:Pos",z)
      println("ν_∞:Phi_e:Pos",ν_neg_∞)
      println("j0:Phi_e:Pos",j0_neg)
      println("Rtot:Phi_e:Pos",Rtot_neg)
      println("σ_eff:Phi_e:Pos",σ_eff_Neg)
      println("∂Uocp_elc:Phi_e:Pos",∂Uocp_neg)
      println("L:Phi_e:Pos",CellData.Neg.L)
      println("D_term:", D_term)
 end
res0 = zeros(length(z))

return ϕ_tf, D, res0, D_term

end

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text:
function evaluator_tests(make_eval::Function)
    @testset "Basic evaluation with a single-objective function to be minimized" begin
        e = make_eval()
        @test BlackBoxOptim.num_evals(e) == 0

        fit1 = fitness([0.0, 1.0], e)
        @test BlackBoxOptim.num_evals(e) == 1
        @test BlackBoxOptim.last_fitness(e) == fit1

        fit2 = fitness([2.0, 1.0], e)
        @test BlackBoxOptim.num_evals(e) == 2
        @test BlackBoxOptim.last_fitness(e) == fit2

        @test BlackBoxOptim.numdims(e) == 2

        @test BlackBoxOptim.is_better(0.0, 1.0, e)
        @test BlackBoxOptim.is_better(0.0, -1.0, e) == false

        a, b = [0.0, 0.5], [1.0, 1.0]

        @test BlackBoxOptim.is_better(a, 1.0, e)
        @test BlackBoxOptim.last_fitness(e) == 0.25

        @test BlackBoxOptim.is_better(b, 1.0, e) == false
        @test BlackBoxOptim.last_fitness(e) == 2.0

        @test BlackBoxOptim.best_of(a, b, e) == (a, 0.25)
        @test BlackBoxOptim.best_of(b, a, e) == (a, 0.25)
        BlackBoxOptim.shutdown!(e)
    end

    @testset "update_fitness!()" begin
        e = make_eval()

        candidates = [random_candidate(2, -1.0, 1.0) for i in 1:10]
        BlackBoxOptim.update_fitness!(e, candidates)
        @test BlackBoxOptim.num_evals(e) == length(candidates)
        @test all(c -> isfinite(c.fitness), candidates)
        BlackBoxOptim.shutdown!(e)
    end

    @testset "rank_by_fitness!()" begin
        e = make_eval()

        candidates = [random_candidate(2, -1.0, 1.0) for i in 1:10]
        # partially evaluate fitness
        BlackBoxOptim.update_fitness!(e, candidates[1:5])
        @test BlackBoxOptim.num_evals(e) == 5
        # complete fitness evaluation and sort by it
        BlackBoxOptim.rank_by_fitness!(e, candidates)
        @test BlackBoxOptim.num_evals(e) == 10
        @test sortperm(candidates, by = fitness) == collect(1:10)
        BlackBoxOptim.shutdown!(e)
    end
end

function multiobj_evaluator_tests(make_eval::Function)
    @testset "multi-objective problem" begin
        schaffer1(x) = (sum(abs2, x), sum(xx -> abs2(xx - 2.0), x))
        p = BlackBoxOptim.FunctionBasedProblem(schaffer1, "Schaffer1", ParetoFitnessScheme{2}(is_minimizing=true),
                                               RectSearchSpace(5, (-10.0, 10.0)))
        a = EpsBoxArchive(EpsBoxDominanceFitnessScheme(fitness_scheme(p)), max_size=100)

        e = make_eval(p, a)
        fit1 = fitness([0.0, 1.0, 2.0, 3.0, 4.0], e)
        @test BlackBoxOptim.num_evals(e) == 1
        @test BlackBoxOptim.last_fitness(e) == fit1.orig

        fit2 = fitness([0.0, -1.0, -2.0, -3.0, -4.0], e)
        @test BlackBoxOptim.num_evals(e) == 2
        @test BlackBoxOptim.last_fitness(e) == fit2.orig

        BlackBoxOptim.shutdown!(e)
    end
end

@testset "Evaluator" begin
    # Set up a small example problem
    f(x) = sum(abs2, x)
    p = minimization_problem(f, "", (-1.0, 1.0), 2)
    @testset "ProblemEvaluator" begin
        evaluator_tests(() -> BlackBoxOptim.ProblemEvaluator(p))
    end

    @testset "rank_by_fitness!()" begin
        e = BlackBoxOptim.ProblemEvaluator(p)

        candidates = [random_candidate(2, -1.0, 1.0) for i in 1:10]
        # partially evaluate fitness
        BlackBoxOptim.update_fitness!(e, candidates[1:5])
        @test BlackBoxOptim.num_evals(e) == 5
        # complete fitness evaluation and sort by it
        BlackBoxOptim.rank_by_fitness!(e, candidates)
        @test BlackBoxOptim.num_evals(e) == 10
        @test sortperm(candidates, by = fitness) == collect(1:10)
    end

    @testset "ParallelEvaluator" begin
        using Distributed

        evaluator_tests(() -> BlackBoxOptim.ParallelEvaluator(p, pids=workers()))
        multiobj_evaluator_tests((p, a) -> BlackBoxOptim.ParallelEvaluator(p, a))
    end

    @testset "MultithreadEvaluator" begin
        if Threads.nthreads() > 1
            evaluator_tests(() -> BlackBoxOptim.MultithreadEvaluator(p))
            multiobj_evaluator_tests((p, a) -> BlackBoxOptim.MultithreadEvaluator(p, a))
        else
            @warn "MultithreadEvaluator tests require >1 threads, $(Threads.nthreads()) found, use JULIA_NUM_THREADS"
            @test_skip Threads.nthreads() > 1
        end
    end
end


text:
# range, description, condition
const SHAPE_DIAGNOSTIC_CATEGORIES = (
    ("(-Inf, 0.5]", "good", x -> !ismissing(x) && x ≤ 0.5),
    ("(0.5, 0.7]", "okay", x -> !ismissing(x) && 0.5 < x ≤ 0.7),
    ("(0.7, 1]", "bad", x -> !ismissing(x) && 0.7 < x ≤ 1),
    ("(1, Inf)", "very bad", x -> !ismissing(x) && x > 1),
    ("——", "missing", ismissing),
)
const BAD_SHAPE_SUMMARY = "Resulting importance sampling estimates are likely to be unstable."
const VERY_BAD_SHAPE_SUMMARY = "Corresponding importance sampling estimates are likely to be unstable and are unlikely to converge with additional samples."

"""
    PSISResult

Result of Pareto-smoothed importance sampling (PSIS) using [`psis`](@ref).

# Properties

  - `log_weights`: un-normalized Pareto-smoothed log weights
  - `weights`: normalized Pareto-smoothed weights (allocates a copy)
  - `pareto_shape`: Pareto ``k=ξ`` shape parameter
  - `nparams`: number of parameters in `log_weights`
  - `ndraws`: number of draws in `log_weights`
  - `nchains`: number of chains in `log_weights`
  - `reff`: the ratio of the effective sample size of the unsmoothed importance ratios and
    the actual sample size.
  - `ess`: estimated effective sample size of estimate of mean using smoothed importance
    samples (see [`ess_is`](@ref))
  - `log_weights_norm`: the logarithm of the normalization constant of `log_weights`
  - `tail_length`: length of the upper tail of `log_weights` that was smoothed
  - `tail_dist`: the generalized Pareto distribution that was fit to the tail of
    `log_weights`. Note that the tail weights are scaled to have a maximum of 1, so
    `tail_dist * exp(maximum(log_ratios))` is the corresponding fit directly to the tail of
    `log_ratios`.

# Diagnostic

The `pareto_shape` parameter ``k=ξ`` of the generalized Pareto distribution `tail_dist` can
be used to diagnose reliability and convergence of estimates using the importance weights
[^VehtariSimpson2021].

  - if ``k < \\frac{1}{3}``, importance sampling is stable, and importance sampling (IS) and
    PSIS both are reliable.
  - if ``k ≤ \\frac{1}{2}``, then the importance ratio distributon has finite variance, and
    the central limit theorem holds. As ``k`` approaches the upper bound, IS becomes less
    reliable, while PSIS still works well but with a higher RMSE.
  - if ``\\frac{1}{2} < k ≤ 0.7``, then the variance is infinite, and IS can behave quite
    poorly. However, PSIS works well in this regime.
  - if ``0.7 < k ≤ 1``, then it quickly becomes impractical to collect enough importance
    weights to reliably compute estimates, and importance sampling is not recommended.
  - if ``k > 1``, then neither the variance nor the mean of the raw importance ratios
    exists. The convergence rate is close to zero, and bias can be large with practical
    sample sizes.

See [`paretoshapeplot`](@ref) for a diagnostic plot.

[^VehtariSimpson2021]: Vehtari A, Simpson D, Gelman A, Yao Y, Gabry J. (2021).
    Pareto smoothed importance sampling.
    [arXiv:1507.02646v7](https://arxiv.org/abs/1507.02646v7) [stat.CO]
"""
struct PSISResult{T,W<:AbstractArray{T},N,R,L,D}
    log_weights::W
    log_weights_norm::N
    reff::R
    tail_length::L
    tail_dist::D
end

function Base.propertynames(r::PSISResult)
    return [fieldnames(typeof(r))..., :weights, :nparams, :ndraws, :nchains, :pareto_shape]
end

function Base.getproperty(r::PSISResult, k::Symbol)
    if k === :weights
        return exp.(getfield(r, :log_weights) .- getfield(r, :log_weights_norm))
    elseif k === :nparams
        log_weights = getfield(r, :log_weights)
        return ndims(log_weights) == 1 ? 1 : size(log_weights, 1)
    elseif k === :ndraws
        log_weights = getfield(r, :log_weights)
        return ndims(log_weights) == 1 ? length(log_weights) : size(log_weights, 2)
    elseif k === :nchains
        log_weights = getfield(r, :log_weights)
        return size(log_weights, 3)
    end
    k === :pareto_shape && return pareto_shape(r)
    k === :ess && return ess_is(r)
    return getfield(r, k)
end

function Base.show(io::IO, ::MIME"text/plain", r::PSISResult)
    npoints = r.nparams
    nchains = r.nchains
    println(
        io, "PSISResult with $npoints parameters, $(r.ndraws) draws, and $nchains chains"
    )
    return _print_pareto_shape_summary(io, r; newline_at_end=false)
end

function pareto_shape_summary(r::PSISResult; kwargs...)
    return _print_pareto_shape_summary(stdout, r; kwargs...)
end

function _print_pareto_shape_summary(io::IO, r::PSISResult; kwargs...)
    k = as_array(pareto_shape(r))
    ess = as_array(ess_is(r))
    npoints = r.nparams
    rows = map(SHAPE_DIAGNOSTIC_CATEGORIES) do (range, desc, cond)
        inds = findall(cond, k)
        count = length(inds)
        perc = 100 * count / npoints
        ess_min = if count == 0 || desc == "too few draws"
            missing
        else
            minimum(view(ess, inds))
        end
        return (range=range, desc=desc, count_perc=(count, perc), ess_min=ess_min)
    end
    rows = filter(r -> r.count_perc[1] > 0, rows)
    formats = Dict(
        "good" => (),
        "okay" => (; color=:yellow),
        "bad" => (bold=true, color=:light_red),
        "very bad" => (bold=true, color=:red),
        "missing" => (),
    )

    col_padding = " "
    col_delim = ""
    col_delim_tot = col_padding * col_delim * col_padding
    col_widths = [
        maximum(r -> length(r.range), rows),
        maximum(r -> length(r.desc), rows),
        maximum(r -> ndigits(r.count_perc[1]), rows),
        floor(Int, log10(maximum(r -> r.count_perc[2], rows))) + 6,
    ]

    println(io, "Pareto shape (k) diagnostic values:")
    printstyled(
        io,
        col_padding,
        " "^col_widths[1],
        col_delim_tot,
        " "^col_widths[2],
        col_delim_tot,
        _pad_right("Count", col_widths[3] + col_widths[4] + 1),
        col_delim_tot,
        "Min. ESS";
        bold=true,
    )
    for r in rows
        count, perc = r.count_perc
        perc_str = "($(round(perc; digits=1))%)"
        println(io)
        print(io, col_padding, _pad_left(r.range, col_widths[1]), col_delim_tot)
        print(io, _pad_right(r.desc, col_widths[2]), col_delim_tot)
        format = formats[r.desc]
        printstyled(io, _pad_left(count, col_widths[3]); format...)
        printstyled(io, " ", _pad_right(perc_str, col_widths[4]); format...)
        print(io, col_delim_tot, r.ess_min === missing ? "——" : floor(Int, r.ess_min))
    end
    return nothing
end

_pad_left(s, nchars) = " "^(nchars - length("$s")) * "$s"
_pad_right(s, nchars) = "$s" * " "^(nchars - length("$s"))

function _promote_result_type(::Type{PSISResult{T,W,N,R,L,D}}) where {T,W,N,R,L,D}
    return PSISResult{T,W,N,R,L,D2} where {D2<:Union{D,Missing,GeneralizedPareto{T}}}
end

"""
    psis(log_ratios, reff = 1.0; kwargs...) -> PSISResult
    psis!(log_ratios, reff = 1.0; kwargs...) -> PSISResult

Compute Pareto smoothed importance sampling (PSIS) log weights [^VehtariSimpson2021].

While `psis` computes smoothed log weights out-of-place, `psis!` smooths them in-place.

# Arguments

  - `log_ratios`: an array of logarithms of importance ratios, with one of the following
    sizes:

      + `(ndraws,)`: a vector of draws for a single parameter from a single chain
      + `(nparams, ndraws)`: a matrix of draws for a multiple parameter from a single chain
      + `(nparams, ndraws, nchains)`: an array of draws for multiple parameters from
        multiple chains, e.g. as might be generated with Markov chain Monte Carlo.

  - `reff::Union{Real,AbstractVector}`: the ratio(s) of effective sample size of
    `log_ratios` and the actual sample size `reff = ess/(ndraws * nchains)`, used to account
    for autocorrelation, e.g. due to Markov chain Monte Carlo.

# Keywords

  - `improved=false`: If `true`, use the adaptive empirical prior of [^Zhang2010].
    If `false`, use the simpler prior of [^ZhangStephens2009], which is also used in
    [^VehtariSimpson2021].
  - `warn=true`: If `true`, warning messages are delivered

# Returns

  - `result`: a [`PSISResult`](@ref) object containing the results of the Pareto-smoothing.

A warning is raised if the Pareto shape parameter ``k ≥ 0.7``. See [`PSISResult`](@ref) for
details and [`paretoshapeplot`](@ref) for a diagnostic plot.

[^VehtariSimpson2021]: Vehtari A, Simpson D, Gelman A, Yao Y, Gabry J. (2021).
    Pareto smoothed importance sampling.
    [arXiv:1507.02646v7](https://arxiv.org/abs/1507.02646v7) [stat.CO]
[^ZhangStephens2009]: Jin Zhang & Michael A. Stephens (2009)
    A New and Efficient Estimation Method for the Generalized Pareto Distribution,
    Technometrics, 51:3, 316-325,
    DOI: [10.1198/tech.2009.08017](https://doi.org/10.1198/tech.2009.08017)
[^Zhang2010]: Jin Zhang (2010) Improving on Estimation for the Generalized Pareto Distribution,
    Technometrics, 52:3, 335-339,
    DOI: [10.1198/TECH.2010.09206](https://doi.org/10.1198/TECH.2010.09206)
"""
psis, psis!

function psis(logr, reff=1; kwargs...)
    T = float(eltype(logr))
    logw = similar(logr, T)
    copyto!(logw, logr)
    return psis!(logw, reff; kwargs...)
end

function psis!(
    logw::AbstractVector,
    reff=1;
    sorted::Bool=false, # deprecated
    improved::Bool=false,
    warn::Bool=true,
)
    S = length(logw)
    reff_val = first(reff)
    M = tail_length(reff_val, S)
    if M < 5
        warn &&
            @warn "$M tail draws is insufficient to fit the generalized Pareto distribution. Total number of draws should in general exceed 25."
        return PSISResult(logw, LogExpFunctions.logsumexp(logw), reff_val, M, missing)
    end
    perm = partialsortperm(logw, (S - M):S)
    cutoff_ind = perm[1]
    tail_inds = @view perm[2:(M + 1)]
    logu = logw[cutoff_ind]
    logw_tail = @views logw[tail_inds]
    if !all(isfinite, logw_tail)
        warn &&
            @warn "Tail contains non-finite values. Generalized Pareto distribution cannot be reliably fit."
        return PSISResult(logw, LogExpFunctions.logsumexp(logw), reff_val, M, missing)
    end
    _, tail_dist = psis_tail!(logw_tail, logu, M, improved)
    warn && check_pareto_shape(tail_dist)
    return PSISResult(logw, LogExpFunctions.logsumexp(logw), reff_val, M, tail_dist)
end
function psis!(logw::AbstractArray, reff=1; warn::Bool=true, kwargs...)
    # allocate containers, calling psis! for first parameter to determine eltype
    logw_firstdraw = first_draw(logw)
    reffs = reff isa Number ? fill!(similar(logw_firstdraw), reff) : reff
    r1 = psis!(vec(param_draws(logw, 1)), reffs[1]; warn=false, kwargs...)
    results = similar(logw_firstdraw, _promote_result_type(typeof(r1)))
    i, inds = Iterators.peel(eachindex(results, reffs))
    results[i] = r1
    # call psis! for remaining parameters
    Threads.@threads for i in collect(inds)
        results[i] = psis!(vec(param_draws(logw, i)), reffs[i]; warn=false, kwargs...)
    end
    # combine results
    logw_norms = map(r -> r.log_weights_norm, results)
    tail_lengths = map(r -> r.tail_length, results)
    tail_dists = map(r -> r.tail_dist, results)
    result = PSISResult(logw, logw_norms, reffs, tail_lengths, tail_dists)
    # warn for bad shape
    warn && check_pareto_shape(result)
    return result
end

pareto_shape(::Missing) = missing
pareto_shape(dist::GeneralizedPareto) = dist.k
pareto_shape(r::PSISResult) = pareto_shape(getfield(r, :tail_dist))
pareto_shape(dists) = map(pareto_shape, dists)

check_pareto_shape(result::PSISResult) = check_pareto_shape(result.tail_dist)
function check_pareto_shape(dist::GeneralizedPareto)
    k = pareto_shape(dist)
    if k > 1
        @warn "Pareto shape k = $(@sprintf("%.2g", k)) > 1. $VERY_BAD_SHAPE_SUMMARY"
    elseif k > 0.7
        @warn "Pareto shape k = $(@sprintf("%.2g", k)) > 0.7. $BAD_SHAPE_SUMMARY"
    end
    return nothing
end
function check_pareto_shape(dists::AbstractVector{<:Union{Missing,GeneralizedPareto}})
    nmissing = count(ismissing, dists)
    ngt07 = count(x -> !(ismissing(x)) && pareto_shape(x) > 0.7, dists)
    ngt1 = iszero(ngt07) ? ngt07 : count(x -> !(ismissing(x)) && pareto_shape(x) > 1, dists)
    if ngt07 > ngt1
        @warn "$(ngt07 - ngt1) parameters had Pareto shape values 0.7 < k ≤ 1. $BAD_SHAPE_SUMMARY"
    end
    if ngt1 > 0
        @warn "$ngt1 parameters had Pareto shape values k > 1. $VERY_BAD_SHAPE_SUMMARY"
    end
    if nmissing > 0
        @warn "For $nmissing parameters, the generalized Pareto distribution could not be fit to the tail draws. Total number of draws should in general exceed 25, and the tail draws must be finite."
    end
    return nothing
end

tail_length(reff, S) = min(cld(S, 5), ceil(Int, 3 * sqrt(S / reff)))

function psis_tail!(logw, logμ, M=length(logw), improved=false)
    T = eltype(logw)
    logw_max = logw[M]
    # to improve numerical stability, we first shift the log-weights to have a maximum of 0,
    # equivalent to scaling the weights to have a maximum of 1.
    μ_scaled = exp(logμ - logw_max)
    w = (logw .= exp.(logw .- logw_max))
    tail_dist_scaled = fit_gpd(w; sorted=true, improved=improved, μ=μ_scaled)
    tail_dist_adjusted = prior_adjust_shape(tail_dist_scaled, M)
    # undo the scaling
    k = tail_dist_adjusted.k
    if isfinite(k)
        p = uniform_probabilities(T, M)
        @inbounds for i in eachindex(logw, p)
            # undo scaling in the log-weights
            logw[i] = min(log(quantile(tail_dist_adjusted, p[i])), 0) + logw_max
        end
    end
    return logw, tail_dist_adjusted
end


text:
""" Functorial data migration for C-sets.
"""
module DataMigration
export Delta, Sigma, migrate!

using ...Present, ...Theories
using ...Theories: SchemaDesc, ob, hom, dom, codom, attr, adom
using ..FinSets, ..CSets, ..Limits, ...Graphs, ..FreeDiagrams
using ..Categories, ..FinCats
import ...Present: Presentation
import ...CategoricalAlgebra.FreeDiagrams: FreeDiagram

const ACSetCat = TypeCat{ACSet,ACSetTransformation}

""" Abstract type for functor that migrates data from one ACSet to another.
"""
abstract type MigrationFunctor <: Functor{ACSetCat,ACSetCat} end

codom_type(F::MigrationFunctor) = F.codom_type

(F::MigrationFunctor)(X::ACSet) = F(codom_type(F)(), X)

# Delta migration
#################

""" Pullback functorial data migration from one ACSet to another.

Note that this operation is contravariant: the data is transferred from `Y` to
`X` but the functor, represented by two dictionaries, maps the schema for `X`
to the schema for `Y`.

When the functor is the identity, this operation is equivalent to
[`copy_parts!`](@ref).
"""
struct DeltaMigration <: MigrationFunctor
  F::FinFunctor # on the schemas
  codom_type::Type
end

"""   Delta(F::FinFunctor)

Given a functor ``F: \\mathcal{C} \\to \\mathcal{D}`` produces a `MigrationFunctor`
which maps a ``\\mathcal{D}``-set ``X`` to the ``\\mathcal{C}``-set
``\\Delta_F(X) := X \\circ F``.

See (Spivak, 2014, *Category Theory for the Sciences*) for details.
"""
Delta(F::FinFunctor, codom) = DeltaMigration(F, codom)
Delta(F::FinFunctor, dom, codom) = DeltaMigration(F, codom)

(ΔF::DeltaMigration)(X::ACSet, Y::ACSet) = migrate!(X, Y, ΔF.F)

"""   migrate!(X::ACSet, Y::ACSet, FOb, FHom)

Migrates the data from `Y` to `X` via the pullback data migration induced by the
functor defined on objects by `FOb` and on morphisms by `FHom`.
"""
function migrate!(X::StructACSet{S}, Y::ACSet, F::FinFunctor) where S
  partsX = Dict(c => add_parts!(X, c, nparts(Y, nameof(ob_map(F,c))))
                for c in ob(S))

  for (f,c,d) in zip(hom(S), dom(S), codom(S))
    set_subpart!(X, partsX[c], f, partsX[d][subpart(Y, hom_map(F,f))])
  end
  for (f,c) in zip(attr(S), adom(S))
    set_subpart!(X, partsX[c], f, subpart(Y, hom_map(F,f)))
  end

  return X
end

function migrate!(X::ACSet, Y::ACSet, FOb, FHom)
  F = FinFunctor(FOb, FHom, FinCat(Presentation(X)), FinCat(Presentation(Y)))
  migrate!(X, Y, F)
end

function (::Type{T})(Y::ACSet, FOb::AbstractDict, FHom::AbstractDict) where T <: ACSet
  X = T()
  migrate!(X, Y, FOb, FHom)
end

function (::Type{T})(Y::ACSet, F::FinFunctor) where T <: ACSet
  X = T()
  migrate!(X, Y, F)
end

# Sigma migration
#################

""" Left pushforward functorial data migration from one ACSet to another.
"""
struct SigmaMigration{CC} <: MigrationFunctor
  F::FinFunctor # on the schemas
  codom_type::Type
  comma_cats::CC

  function SigmaMigration(F::FinFunctor, codom_type)
    cc = get_comma_cats(F)
    new{typeof(cc)}(F, codom_type, cc)
  end
end

"""   Sigma(F::FinFunctor)

Given a functor ``F: \\mathcal{C} \\to \\mathcal{D}`` produces a `MigrationFunctor`
which maps a ``\\mathcal{C}``-set ``X`` to the ``\\mathcal{D}``-set ``\\Sigma_F(X)``.
``\\Sigma_F`` is left-adjoint to ``\\Delta_F``.

For ``d \\in \\mathsf{ob}\\mathcal{D}``,
``\\Sigma_F(X)(d) := \\mathsf{colim}_{F \\downarrow d} X \\circ \\pi`` where
``\\pi: F\\downarrow d \\to \\mathcal{C}`` is the projection.

See (Spivak, 2014, *Category Theory for the Sciences*) for details.
"""
Sigma(F::FinFunctor, codom) = SigmaMigration(F, codom)
Sigma(F::FinFunctor, dom, codom) = SigmaMigration(F, codom)

function (ΣF::SigmaMigration)(Y::ACSet, X::ACSet)
  comma_cats = ΣF.comma_cats
  diagramD = FreeDiagram(presentation(codom(ΣF.F)))

  # define Y on objects by taking colimits
  colimX = map(parts(diagramD, :V)) do i
    F∇d = ob(comma_cats, i)
    Xobs = map(ob(F∇d)) do (c,_)
      FinSet{Int,Int}(nparts(X, nameof(c)))
    end
    Xhoms = map(parts(F∇d, :E)) do g
      FinFunction(X[nameof(hom(F∇d, g))], Xobs[src(F∇d, g)], Xobs[tgt(F∇d, g)])
    end
    colimit(FreeDiagram(Xobs, collect(zip(Xhoms, src(F∇d), tgt(F∇d)))))
  end

  for d in parts(diagramD, :V)
    add_parts!(Y, nameof(ob(diagramD, d)), length(apex(colimX[d])))
  end

  # Define Y on morphisms by the universal property
  for g in parts(diagramD, :E)
    if nparts(Y, nameof(dom(hom(diagramD, g)))) == 0
      continue
    end
    set_subpart!(Y, nameof(hom(diagramD, g)),
      collect(universal(colimX[src(diagramD, g)],
        Multicospan(legs(colimX[tgt(diagramD, g)])[collect(hom(comma_cats, g)[:V])])
      )))
  end
  return Y
end

""" add_hom!(d::FreeDiagram{Ob, Hom}, src_ob::Ob, tgt_ob::Ob, hom::Hom)

Adds a hom to `d` from the vertex with object `src_ob` to the vertex with object `tgt_ob`.
"""
function add_hom!(d::FreeDiagram, src_ob, tgt_ob, hom)
  src_idx = first(incident(d, src_ob, :ob))
  tgt_idx = first(incident(d, tgt_ob, :ob))
  return add_edge!(d, src_idx, tgt_idx, hom = hom)
end

"""   comma_cats(diagramD::FreeDiagram{FreeSchema.Ob, FreeSchema.Hom}, FOb, FHom)

Given a free category ``\\mathcal{D}`` with no cycles and
a functor represented by the pair `(FOb, FHom)`, returns a diagram
``\\mathcal{D} \\to \\mathsf{Cat}`` which on objects takes ``d \\in \\mathcal{D}`` to the
comma category ``F \\downarrow d``.
"""
function get_comma_cats(F::FinFunctor)
  diagramD = FreeDiagram(presentation(codom(F)))
  FObInv = y -> filter(x -> ob_map(F,x) == y, ob_generators(dom(F)))
  FHomInv = g -> filter(f -> hom_map(F,f) == g, hom_generators(dom(F)))
  comma_cats = FreeDiagram{FreeDiagram, ACSetTransformation}()
  add_vertices!(comma_cats,
    nparts(diagramD, :V),
    ob = map(parts(diagramD, :V)) do _
      FreeDiagram{Pair{FreeSchema.Ob, FreeSchema.Hom}, FreeSchema.Hom}()
    end
  )

  for d in topological_sort(diagramD)
    F∇d = ob(comma_cats, d)
    id_d = id(ob(diagramD, d))

    # If FOb[c] = d, add an objects (c, id(d)) to F∇d
    preimage_d = FObInv(ob(diagramD, d))
    id_obs = add_vertices!(F∇d, length(preimage_d), ob = map(c->c=>id_d, preimage_d))

    # if FHom[h] = id(d), add a hom h: (dom(h), id(d)) -> (codom(h), id(d)) to F∇d
    id_homs = map(FHomInv(id_d)) do h
      add_hom!(F∇d, dom(h) => id_d, codom(h) => id_d, h)
    end

    # for g: d' -> d in D (note that F∇d' must have already been constructed)
    #     populate F∇d with obs and homs coming from F∇d′
    for g in incident(diagramD, d, :tgt)
      d′ = src(diagramD, g)
      F∇g = comma_cat_hom!(F∇d, ob(comma_cats, d′), id_d, hom(diagramD, g), FHomInv)
      add_edge!(comma_cats, d′, d, hom=F∇g)
    end
  end

  return comma_cats
end

function comma_cat_hom!(F∇d, F∇d′, id_d, g, FHomInv)
  # for (c,f) in F∇d' add ob (c, compose(f,g)) to F∇d
  vs = add_vertices!(F∇d, nparts(F∇d′, :V), ob = map(ob(F∇d′)) do (c,f)
    c => compose(f, g)
  end)

  # for h: (c, f) -> (c', f') in diagram_d', add hom h to F∇d
  es = add_edges!(F∇d, vs[src(F∇d′)], vs[tgt(F∇d′)], hom = hom(F∇d′))

  # for every (c,f) in the recently added objects. If FHom[h] == f, then add the hom
  #       h: (c,f) -> (codom(h), idv)
  # relies on D being free
  for (c, f) in ob(F∇d, vs)
    for h in FHomInv(f)
      dom(h) == c && add_hom!(F∇d, c => f, codom(h) => id_d, h)
    end
  end

  # return the inclusion from F∇d′ to F∇d
  return ACSetTransformation((V = collect(vs), E = collect(es)), F∇d′, F∇d)
end

### Translate between ACSets, Presentations, and FreeDiagrams
#############################################################

# FIXME: These functions do not belong here.

""" Get the Schema from an ACSet
"""
function Presentation(::StructACSet{S}) where S
  return Presentation(S)
end

"""   FreeDiagram(pres::Presentation{FreeSchema, Symbol})

Returns a `FreeDiagram` whose objects are the generating objects of `pres` and
whose homs are the generating homs of `pres`.
"""
function FreeDiagram(pres::Presentation{Schema, Symbol}) where Schema
  obs = Array{FreeSchema.Ob}(generators(pres, :Ob))
  homs = Array{FreeSchema.Hom}(generators(pres, :Hom))
  doms = map(h -> generator_index(pres, nameof(dom(h))), homs)
  codoms = map(h -> generator_index(pres, nameof(codom(h))), homs)
  return FreeDiagram(obs, collect(zip(homs, doms, codoms)))
end

end


text:
using Test
using LightGraphs
using SimpleTensorNetworks
using Random

@testset "greedy" begin
    include("edgeinfo.jl")
end

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

@testset "other" begin
    for seed in 1:10
        Random.seed!(3)
        g = random_regular_graph(10, 3)
        tn = TensorNetwork([LabeledTensor(randn(2,2,2), [e for e in edges(g) if i ∈ (e.src, e.dst)]) for i=1:10])
        tcs, scs, order = SimpleTensorNetworks.order_greedy(SimpleTensorNetworks.layout(tn); strategy="min_reduce")
        trees = SimpleTensorNetworks.build_trees(10, order)
        order_new = SimpleTensorNetworks.build_order(trees)
        #@test order == order_new
        tc, sc = abstract_contract(SimpleTensorNetworks.layout(tn), order)
        tc2, sc2 = abstract_contract(tn, trees)
        @test sort(tc) ≈ sort(tc2)
    end

    @testset "log2sumexp2" begin
        x = randn(10)
        @test log2(sum(exp2.(x))) ≈ log2sumexp2(x)
    end
end