text:
using Plots
include("bitcoin_btg.jl")
include("../../computation/tdist_mle.jl")

path = "../../datasets/"
data, target = bitcoin(path)

btg0 = load_bitcoin_btg()
(pdf0_raw, cdf0_raw, dpdf0_raw, quantInfo0_raw) = solve(btg0);
f = theta -> tdist_mle(btg0, theta, 1.0);

####################################
######## Optional Plotting #########
####################################
if false
    Plots.plot()
    plt(y -> fzero(x -> cdf0_raw(reshape([y], 1, 1), reshape([1.0], 1, 1), x) - 0.5, 0, 2), 0, 1)
    plt!(y -> fzero(x -> cdf0_raw(reshape([y], 1, 1), reshape([1.0], 1, 1), x) - 0.9, 0, 2), 0, 1)
    plt!(y -> fzero(x -> cdf0_raw(reshape([y], 1, 1), reshape([1.0], 1, 1), x) - 0.1, 0, 2), 0, 1)
end

####################################
############### Test ###############
####################################
if true #regular GP Test
    println("Start Vanilla GP Test")
    x = x_train #reshape(x_train, length(x_train) ,1)
    y = y_train #reshape(y_train, length(y_train), 1)
    # build and fit a GP
    #mymean = MeanLin(zeros(d))
    mymean = MeanZero()
    # mymean = MeanZero()
    kern = SE(0.0, 0.0)
    gp = GP(x, y, mymean, kern)
    optimize!(gp)
    # predict
    μ, σ² = predict_y(gp, dropdims(x_test, dims=2)); stdv = sqrt.(σ²)
    error_GP = abs.(μ .- y_test)
    error_abs_GP = mean(error_GP)
    error_sq_GP = mean(error_GP.^2)
    CI_test_GP = hcat(-1.96.*stdv .+ μ, 1.96.*stdv .+ μ)
    count_test_GP = sum((y_test .>= CI_test_GP[:, 1]) .* (y_test .<= CI_test_GP[:,2]))/n_test
    nlpd_GP = -mean(log.(pdf.(Normal(), (y_test.-μ)./stdv)./stdv))
    println("MSE:")
    display(error_sq_GP)
    println("Conf interval accuracy:")
    display(count_test_GP)
end

#
# MSE:
# 0.0017970008992554704
# Conf interval accuracy:
# 0.941033434650456
#
#

if true
    @info "Start BTG Test"
    before = Dates.now()
    count_test = 0
    error_abs = 0.
    error_sq = 0.
    nlpd = 0.
    id_fail = []
    id_nonproper = []
    CI_set = zeros(n_test, 2)
    median_set = zeros(n_test)
    for i in 1:n_test
        global error_abs, error_sq, nlpd, count_test
        # mod(i, 20) == 0 ? (@info i) : nothing
        # @info "i" i
        x_test_i = reshape(x_test[i, :], 1, d)
        Fx_test_i = reshape(Fx_test[i, :], 1, p)
        try
            #pdf_test_i, cdf_test_i, dpdf_test_i, quantbound_test_i, support_test_i = pre_process(x_test_i, Fx_test_i, pdf0_raw, cdf0_raw, dpdf0_raw, quantInfo0_raw)
            y_test_i = y_test[i]

            #median_test_i = max_train * quantile(cdf_test_i, quantbound_test_i, support_test_i)[1]
            median_test_i = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.5, 0, 10)
            # @info "True, median " y_test_i_true, median_test_i
            median_set[i] = median_test_i
            try
                #CI_test_i = max_train .* credible_interval(cdf_test_i, quantbound_test_i, support_test_i; mode=:equal, wp=.95)[1]
                aa = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.975, 0, 20)
                bb = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.025, 0, 20)
                CI_test_i = [bb, aa]
                count_test += (y_test_i >= CI_test_i[1])&&(y_test_i <= CI_test_i[2]) ? 1 : 0
                # @info "CI" CI_test_i
                #@info "95% confidence interval:", CI_test_i
                CI_set[i, :] = CI_test_i
            catch err
                append!(id_fail, i)
            end
            error_abs += abs(y_test_i - median_test_i)
            error_sq += (y_test_i - median_test_i)^2
            nlpd += log(pdf0_raw(y_test_i))
        # @info "Count, id_fail" count_test, id_fail
        catch err
            append!(id_nonproper, i)
        end
    end
        count_test /= n_test
        error_abs  /= n_test
        error_sq   /= n_test
        nlpd       /= -n_test
        after = Dates.now()
        elapsedmin = round(((after - before) / Millisecond(1000))/60, digits=5)
        println("error_sq")
        display(error_sq)
        println("count_test")
        display(count_test)
        println("elapsedmin")
        display(elapsedmin)
        # Plot
        PyPlot.close("all") #close existing windows
        PyPlot.plot(x_test, median_set, label = "BTG median")
        PyPlot.plot(xs, label, label = "true")
        PyPlot.fill_between(dropdims(x_test; dims = 2), CI_set[:, 1], CI_set[:, 2], alpha = 0.3, label = "95% confidence interval")
        PyPlot.scatter(x_train, y_train, s = 10, c = "k", marker = "*")
        #PyPlot.scatter(x_test, y_test, s = 10, c = "k", marker = "*")
        PyPlot.legend(fontsize=8)
        PyPlot.grid()
        PyPlot.title("BTG $myquadtype", fontsize=10)
        PyPlot.savefig("exp_synthetic_mle_btg4.pdf")
    end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
function acv2cov(acv)
# Construct a stationary covariance matrix from a vector of autovariances
 r = size(acv,1);
 covmat = NaN*zeros(r,r);
 covmat[:, 1] = acv;
 covmat[:, end] = acv[end:-1:1]
 for i in 2:r-1
    covmat[:,i] = [acv[i:-1:1]; acv[2:(r-i+1)]];
 end
return covmat
end

ext: jl


text:
function accumarray1d(A, subs, val)
    for i = 1:length(val)
        @inbounds A[subs[i]] += val[i]
    end
    return A
end

function accumarray2d(A, subs, val)
    for i = 1:length(val)
        A[:,subs[i]] = A[:,subs[i]] .+ val[:,i]
    end
    return A
end

function accumarray3d(A, subs, val)
    for i = 1:length(val)
        A[:,subs[i],:] = A[:,subs[i],:] .+ val[:,i,:]
    end
    return A
end

function accumarray(subs, val)
    n = maximum(subs)
    A = zeros(eltype(val), n)
    for i = 1:length(val)
        @inbounds A[subs[i]] += val[i]
    end
    A
end

function accumarray2d(subs, val)
    n = maximum(subs)
    m = size(val,1)
    A = zeros(eltype(val), (m, n))
    for i = 1:length(val)
        A[:,subs[i]] = A[:,subs[i]] .+ val[:,i]
    end
    return A
end

function accumarray3d(subs, val)
    n = maximum(subs)
    m = size(val,1)
    p = size(val,3)
    A = zeros(eltype(val), (m, n, p))
    for i = 1:length(val)
        A[:,subs[i],:] = A[:,subs[i],:] .+ val[:,i,:]
    end
    return A
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
function check_clock(cur_player::AbstractPlayer, cur_clock_symbol::Union{Void, Symbol}=nothing, expected_uuid::Union{Void, Base.Random.UUID}=nothing)
  cur_state = cur_player.state

  cur_state.is_playing || return false

  cur_state.is_done && return false

  ( cur_clock_symbol == nothing ) && return true

  actual_uuid = getfield(cur_player.clock, cur_clock_symbol)

  actual_uuid == expected_uuid
end


text:
const B23_n = [0.348_051_856_289_69E3,
                -0.116_718_598_799_75E1,
                0.101_929_700_393_26E-2,
                0.572_544_598_627_46E3,
                0.139_188_397_788_70E2]
"""
    B23(::Union{ThermoState.Types.Pressure,ThermoState.Types.Temperature}, InputValue)

Returns the boundary between regions 2 and 3.
InputType is either :T or :P to indicate that InputValue is temperature [K] or pressure [MPa]. The complimentary value is returned.
"""

function B23(InputType::Pressure, InputValue)
    n = B23_n
    return n[4] + √((InputValue - n[5])/n[3])
end

function B23(InputType::Temperature, InputValue)
    n = B23_n
    return n[1] + n[2] * InputValue + n[3] * InputValue^2
end

const B2bc_n = [0.905_842_785_147_23E3,
                -0.679_557_863_992_41,
                0.128_090_027_301_36E-3,
                0.265_265_719_084_28E4,
                0.452_575_789_059_48E1]
"""
    B2bc

    Returns the boundary between regions 2b and 2c, approx s=5.85kJ/kgK
    InputType is either :h or :P to indicate that InputValue is enthalpy [kJ/kg]
    or pressure [MPa]. The complimentary value is returned.
    Valid from saturation line at 554.485K & 6.54670MPa to 1019.32K & 100MPa
"""
function B2bc(InputType::Enthalpy, InputValue)
    n = B2bc_n
    return n[1] + n[2] * InputValue + n[3] * InputValue^2
end

function B2bc(InputType::Pressure, InputValue)
    n = B2bc_n
    return n[1] + n[2] * InputValue + n[3] * InputValue^2
end

function temperature_impl(mt::SingleSatP,model::IndustrialWater,p)
P = normalize_units(p)
Pc = pressure(model,CriticalPoint())
Pmin = 611.213 #different than triple point pressure
    if Pmin <= P <= Pc
        return temperature_impl(mt,IF97Region{:r4}(),P)
    else
        return DomainError(p, "Pressure not between triple and critical points")
    end
end

#temperature_impl(mt::SingleSatP,model::IndustrialWater,p)
function pressure_impl(mt::SingleSatT,model::IndustrialWater,t)
    T = normalize_units(t)
    Tc = temperature(model,CriticalPoint())
    T3= temperature(model,TriplePoint())
    if T3 <= T <= Tc
        return pressure_impl(mt,IF97Region{:r4}(),T)
    else
        return DomainError(t, "temperature not between triple and critical points")
    end
end

pressure_impl(mt::SingleΦT,model::IndustrialWater,t) = pressure_impl(QuickStates.sat_p(),model,t)
temperature_impl(mt::SingleΦP,model::IndustrialWater,p) = temperature_impl(QuickStates.sat_t(),model,p)

function pressure(model::IndustrialWater,st::ThermodynamicState,unit = u"Pa")
    return pressure(state_type(st),model,st,unit)
end

function pressure(mt::SingleSatT,model::IndustrialWater,st::ThermodynamicState,unit)
    t = temperature(FromState(),st)
    res = pressure_impl(mt,model,t)
    return convert_unit(u"Pa",unit,res)
end

function temperature(model::IndustrialWater,st::ThermodynamicState,unit = u"K")
    return pressure(state_type(st),model,st,unit)
end

function temperature(mt::SingleSatP,model::IndustrialWater,st::ThermodynamicState,unit)
    t = pressure(FromState(),st)
    res = temperature_impl(mt,model,t)
    return convert_unit(u"K",unit,res)
end

function region_id(mt::SinglePT,model::IndustrialWater,p, t)::Symbol
    function Psat(Tx)
        p = pressure_impl(QuickStates.sat_t(),model,Tx)
        return convert_unit(u"Pa",u"MPa",p)
    end
    TSat(Px) =  temperature_impl(QuickStates.sat_t(),model,Px)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    #=
    Region 1:
        273.15K ≤ T ≤ 623.15K  Psat(T) ≤ P ≤ 100MPa

    Region 2:
    273.15K ≤ T ≤ 623.15K  0 ≤ P ≤ Psat(T)
    623.15K <  T ≤ 863.15K  0 <  P ≤ P(T) from B23-model
    863.15K <  T ≤ 1073.15K 0 <  P ≤ 100MPa

    Region 2meta:
    From the saturated-vapour line to the 5% equlibrium moisture line, a.k.a the
    practical Wilson line:
        %equilibrium moisture = [h - h_liq(P)] / [h_vap(P) - h_liq(P)]
    for pressures from the triple point (273.16K, 611.657MPa) up to 10MPa.
    These equations only used on demand.

    Region 3:
        623.15K ≤ T ≤ T(P) from B23-model P(T) from B23-model ≤ P ≤ 100MPa

    Region 4:
    Returns the vapour-liquid phase boundary (saturation line).
    Valid from triple point to critical point
    273.15K ≤ T ≤ 647.096K
    Since this is a single floating point number, Region 4 is never returned and only used on demand.

    Region 5:emperature in [K]
    =#

    #Check region 5 first:
    if 1073.15 ≤ T ≤ 2273.15 && 0 ≤ P ≤ 50
        return :r5
    elseif T < 273.15 || T > 1073.15 || P < 0 || P > 100
        throw(DomainError((p, t), "Pressure/Temperature outside valid ranges.")) #Outside the regions where equations are available.
    elseif 273.15 ≤ T ≤ 623.15
        if Psat(T) ≤ P ≤ 100
            return :r1
        else
            return :r2
        end
    elseif 623.15 <  T ≤ 863.15
        if 0 <  P ≤ B23(Temperature(), T)
            return :r2
        else
            return :r3
        end
    else
        return :r2
    end
end
function region_id(mt::SinglePH,model::IndustrialWater,p, _h)::Symbol

    function Psat(Tx)
        p = pressure_impl(QuickStates.sat_t(),model,Tx)
        return convert_unit(u"Pa",u"MPa",p)
    end
    TSat(Px) =  temperature_impl(QuickStates.sat_t(),model,Px)
    P = convert_unit(u"Pa",u"MPa",p)
    h = convert_unit(u"kJ/g",u"kJ/kg",_h)
    #=
        Region 1:
        273.15K ≤ T ≤ 623.15K  Psat(T) ≤ P ≤ 100MPa

        Region 2:
        273.15K ≤ T ≤ 623.15K  0 ≤ P ≤ Psat(T)
        623.15K <  T ≤ 863.15K  0 <  P ≤ P(T) from B23-model
        863.15K <  T ≤ 1073.15K 0 <  P ≤ 100MPa
        For reverse:
            2a: P ≤ 4Mpa
            2b: s ≥ 5.85 kJ/kgK (or use function B2bc if P,h specificed)
            2c: s ≥ 5.85 kJ/kgK (or use function B2bc if P,h specificed)
    =#

    # Check overall region first
    if P > 100.0
        throw(DomainError(P, "Pressure not in valid ranges."))
    end

    T = temperature_impl(mt,IF97Region{:r1}(),p, _h)
    if 273.15 ≤ T ≤ 623.15
        # could be Region 1
        if P ≥ Psat(T)
            return :r1 #, Region1(:SpecificH, P, T) # Return forward h for consistency check
        end
    end

    if P ≤ 4.0
        # could be Region 2a
        T = temperature_impl(mt,IF97Region{:r2a}(),p, _h)
        if Tsat(P) ≤ T ≤ 1073.15
            return :r2a #, Region2(:SpecificH, P, T) # Return forward h for consistency check
        else
            throw(DomainError((P, T), "Pressure/Temperature outside valid ranges.")) # Only other region with backwards mdoels for P ≤ 4.0 is Region 1, which is already eliminated
        end
    else
        htest = B2bc(Pressure(), P)
        if h < htest
            # could be Region 2c
            T = temperature_impl(mt,IF97Region{:r2c}(),p, _h)
            if P ≤ Psat(623.15)
                if Tsat(P) ≤ T ≤ 1073.15
                    return :r2c #, Region2(:SpecificH, P, T) # Return forward h for consistency check
                end
            elseif B23(Pressure(), P) ≤ T ≤ 1073.15
                return :r2c #, Region2(:SpecificH, P, T) # Return forward h for consistency check
            else
                throw(DomainError((P, T), "Pressure/Temperature outside valid ranges."))
            end
        else
            # could be Region 2b
            T = temperature_impl(mt,IF97Region{:r2b}(),p, _h)
            if P ≤ Psat(623.15)
                if Tsat(P) ≤ T ≤ 1073.15
                    return :r2b #, Region2(:SpecificH, P, T) # Return forward h for consistency check
                end
            elseif B23(Pressure(), P) ≤ T ≤ 1073.15
                return :r2b #, Region2(:SpecificH, P, T) # Return forward h for consistency check
            else
                throw(DomainError((P, T), "Pressure/Temperature outside valid ranges."))
            end
        end
    end
    throw(DomainError((P, T), "Pressure/Temperature outside valid ranges."))
end

function region_id(mt::SinglePS,model::IndustrialWater,p, _s)::Symbol
    function Psat(Tx)
        p = pressure_impl(QuickStates.sat_t(),model,Tx)
        return convert_unit(u"Pa",u"MPa",p)
    end
    TSat(Px) =  temperature_impl(QuickStates.sat_t(),model,Px)

    P = convert_unit(u"Pa",u"MPa",p)
    s = convert_unit(u"J/(kg*K)",u"kJ/(kg*K)",_s)
    #=
        Region 1:
        273.15K ≤ T ≤ 623.15K  Psat(T) ≤ P ≤ 100MPa

        Region 2:
        273.15K ≤ T ≤ 623.15K  0 ≤ P ≤ Psat(T)
        623.15K <  T ≤ 863.15K  0 <  P ≤ P(T) from B23-model
        863.15K <  T ≤ 1073.15K 0 <  P ≤ 100MPa
        For reverse:
            2a: P ≤ 4Mpa
            2c: s ≥ 5.85 kJ/kgK
            2b: s ≥ 5.85 kJ/kgK
    =#

    # Check overall region first
    if P > 100
        throw(DomainError((P, s), "Pressure/entropy not in valid ranges."))
    end

    T = temperature_impl(mt,IF97Region{:r1}(),p, _s)

    if 273.15 ≤ T ≤ 623.15
        # could be Region 1
        if P ≥ Psat(T)
            return :r1 #, Region1(:SpecificS, P, T) # Return forward s for consistency check
        end
    end

    if P ≤ 4.0
        # could be Region 2a
        T = temperature_impl(mt,IF97Region{:r2a}(),p, _s)
        if Tsat(P) ≤ T ≤ 1073.15
            return :r2a #, Region2(:SpecificS, P, T) # Return forward s for consistency check
        else
            throw(DomainError((P, s), "Pressure/entropy not in valid ranges.")) # Only other region with backwards mdoels for P ≤ 4.0 is Region 1, which is already eliminated
        end
    else
        if s < 5.85
            # could be Region 2c
            T = temperature_impl(mt,IF97Region{:r2c}(),p, _s)
            if P ≤ Psat(623.15)
                if Tsat(P) ≤ T ≤ 1073.15
                    return :r2 #c, Region2(:SpecificS, P, T) # Return forward s for consistency check
                end
            elseif B23(Pressure(),P) ≤ T ≤ 1073.15
                return :r2c #, Region2(:SpecificS, P, T) # Return forward s for consistency check
            else
                throw(DomainError((P, s), "Pressure/entropy not in valid ranges."))
            end
        else
            # could be Region 2b
            T = temperature_impl(mt,IF97Region{:r2b}(),p, _s)
            if P ≤ Psat(623.15)
                if Tsat(P) ≤ T ≤ 1073.15
                    return :r2b #, Region2(:SpecificS, P, T) # Return forward s for consistency check
                end
            elseif B23(Pressure(), P) ≤ T ≤ 1073.15
                return :r2b #, Region2(:SpecificS, P, T) # Return forward s for consistency check
            else
                throw(DomainError((P, s), "Pressure/entropy not in valid ranges."))

            end
        end
    end
    throw(DomainError((P, s), "Pressure/entropy not in valid ranges."))
end


text:
using Test
using HyperbolicPlane

p = HPoint(0)
q = HPoint(.2)
r = HPoint(0-0.2im)
@test dist(p,q) == dist(r)

p[:color] = "red"
@test p[:color] == "red"

P = RandomHPoint()
Q = RandomHPoint()
S = P+Q
L = HLine(P,Q)
@test issubset(S,L)
M = midpoint(S)
@test in(P,L)
@test in(P,S)
@test issubset(P+M,S)

@test length(S) == dist(P,Q)

x1 = length(S)
x2 = length(P+M) + length(M+Q)
@test abs(x1-x2) < 100*eps(1.)

T = P+Q+M
@test area(T) < 100*eps(1.0)

L = HLine(1,3)
LL = HLine(2,pi)
P = meet(L,LL)
@test in(P,L)
@test in(P,LL)

TT = reflect_across(T,L)
TT = reflect_across(TT,L)
@test T == TT

S = RandomHSegment()
L = HLine(S)
SS = S'
LL = HLine(SS)
@test L == LL'

SS = -S
LL = HLine(SS)
@test L == -LL

a,b,c,d = [RandomHPoint() for t=1:4]
X = HPolygon(a,b,c,d,HPoint(im*eps(1.0)))
Y = HPolygon(d,c,b,a,HPoint(0))
@test X==Y

C = HContainer()
P = RandomHPolygon(5)
for k=1:10
    add_object!(C,P)
end

T = RandomHTriangle()
add_object!(C,T)
TT = T'
add_object!(C,TT')
@test length(C) == 2

a,b,c = endpoints(T)
@test a+c+b == T

P = RandomHPolygon(3)
T = HTriangle(P)
aP = sort(angles(P))
aT = angles(T)
@test sum(abs.(aP-aT)) < 100*eps(1.0)

T = HPoint(1.0,2.0) + HPoint(3.0,-1) + HPoint(2,0)
P = HPolygon(T)
@test perimeter(T) == perimeter(P)

# point on line stuff
L = RandomHLine()
P = point_on_line(L)
@test in(P,L)

PP = points_on_line(L,2)
P,Q = PP
@test in(P,L) && in(Q,L)

R = RandomHRay()
V = get_vertex(R)
P = point_on_ray(R)
@test in(P,R)
A = midpoint(V,P)
@test in(A,R)

# Polygon area
X = equilateral(6,1)
T = X.plist[1] + X.plist[2] + HPoint()
AX = area(X)
AT = area(T)
@test abs(AX-6AT) <= equality_threshold() * eps(1.0)


text:
export isbtcaddress

function isbtcaddress(str::AbstractString)::Bool
    bench32Reg = r"^(bc1)[a-z0-9]{25,39}$"
    base58Reg = r"^(1|3)[A-HJ-NP-Za-km-z1-9]{25,39}$"

    if startswith(str, "bc1")
        return contains(str, bench32Reg)
    end

    return contains(str, base58Reg)
end


text:
(\)(A::AbstractTriangular{<:CuSparseMatrix},B::CuMatrix)       = sm('N',A,B,'O')
(\)(transA::Transpose{<:Any, AbstractTriangular{<:CuSparseMatrix}}, B::CuMatrix) = sm('T',transA.parent,B,'O')
(\)(adjA::Adjoint{<:Any, AbstractTriangular{<:CuSparseMatrix}},B::CuMatrix) = sm('C',adjA.parent,B,'O')

(*)(A::CuSparseMatrix,B::CuVector)       = mv('N',A,B,'O')
(*)(transA::Transpose{<:Any, CuSparseMatrix},B::CuVector)  = mv('T',transA.parent,B,'O')
(*)(adjA::Adjoint{<:Any, CuSparseMatrix},B::CuVector)  = mv('C',transA.parent,B,'O')
(*)(A::HermOrSym{T,CuSparseMatrix{T}},B::CuVector{T}) where T = mv('N',A,B,'O')
(*)(transA::Transpose{<:Any, HermOrSym{T,CuSparseMatrix{T}}},B::CuVector{T}) where T = mv('T',transA.parent,B,'O')
(*)(adjA::Adjoint{<:Any, HermOrSym{T,CuSparseMatrix{T}}},B::CuVector{T}) where T = mv('C',adjA.parent,B,'O')

(*)(A::CuSparseMatrix{T},B::CuMatrix{T})       where T = mm2('N','N',A,B,'O')
(*)(A::CuSparseMatrix{T},transB::Transpose{<:Any, CuMatrix{T}})  where T = mm2('N','T',A,transB.parent,'O')
(*)(transA::Transpose{<:Any, CuSparseMatrix{T}},B::CuMatrix{T})  where T = mm2('T','N',transA.parent,B,'O')
(*)(transA::Transpose{<:Any, CuSparseMatrix{T}},transB::Transpose{<:Any, CuMatrix{T}}) where T = mm2('T','T',transA.parent,transB.parent,'O')
(*)(adjA::Adjoint{<:Any, CuSparseMatrix{T}},B::CuMatrix{T})  where T = mm2('C','N',adjA.parent,B,'O')

(*)(A::HermOrSym,B::CuMatrix) = mm('N',A,B,'O')
(*)(transA::Transpose{<:Any, <:HermOrSym{<:Number, <:CuSparseMatrix}},B::CuMatrix) = mm('T',transA.parent,B,'O')
(*)(adjA::Adjoint{<:Any, <:HermOrSym{<:Number, <:CuSparseMatrix}},B::CuMatrix) = mm('C',adjA.parent, B,'O')

(\)(A::AbstractTriangular{<:CuSparseMatrix},B::CuVector)       = sv2('N',A,B,'O')
(\)(transA::Transpose{<:Any, AbstractTriangular{<:CuSparseMatrix}},B::CuVector) = sv2('T',transA.parent,B,'O')
(\)(adjA::Adjoint{<:Any, AbstractTriangular{T,CuSparseMatrix{T}}},B::CuVector{T}) where T = sv2('C',adjA.parent,B,'O')
(\)(A::AbstractTriangular{T,CuSparseMatrixHYB{T}},B::CuVector{T})       where T = sv('N',A,B,'O')
(\)(transA::Transpose{<:Any, AbstractTriangular{T,CuSparseMatrixHYB{T}}},B::CuVector{T}) where T = sv('T',transA.parent,B,'O')
(\)(adjA::Adjoint{<:Any, AbstractTriangular{T,CuSparseMatrixHYB{T}}},B::CuVector{T}) where T = sv('C',adjA.parent,B,'O')

(+)(A::Union{CuSparseMatrixCSR,CuSparseMatrixCSC},B::Union{CuSparseMatrixCSR,CuSparseMatrixCSC}) = geam(A,B,'O','O','O')
(-)(A::Union{CuSparseMatrixCSR,CuSparseMatrixCSC},B::Union{CuSparseMatrixCSR,CuSparseMatrixCSC}) = geam(A,-one(eltype(A)),B,'O','O','O')


text:
function io_model(t::Trace, dim::Int64, data::Int64)
    z = sample(t, "z", Dirichlet(ones(dim)))
    observe(t, "x", Categorical(z), data)
end

@testset "io 1" begin
    t = trace()
    dim = 10
    data = 7
    io_model(t, dim, data)
    testpath = joinpath(@__DIR__, "TESTIO")
    t_table = to_table(t)
    @info "Created JuliaDB table: $t_table"

    mkpath(testpath)
    db_file = joinpath(testpath, "test.jdb")
    save(t, db_file)
    loaded_trace = load(db_file)
    @info "z node: $(loaded_trace["z"])"
    @info "x node: $(loaded_trace["x"])"
    @test loaded_trace["x"].address == "x"
    @test typeof(loaded_trace["x"].dist) == Categorical{Float64, Vector{Float64}}

    # test save / load of results
    results = mh(io_model; params = (dim, data), burn = 0, thin = 1, num_iterations=10)
    results_file = joinpath(testpath, "io_model.csm")
    save(results, results_file)
    loaded_results = load(joinpath(testpath, "io_model.csm"))

    # use saved results for prediction
    updated_model = update(io_model, loaded_results)
    new_t, _ = updated_model(trace(), dim, data)
    @info "Trace created using saved results: $new_t"
    @test new_t["z"].interpretation == EMPIRICAL
    rm(testpath, recursive=true)
end


text:
@doc doc"""Represents a NeoVim window.""" ->
type Window
    nvim::Nvim
    ext::MsgPack.Ext
end


text:
"""
    abstract_eval(ctx::LintContext, ex) :: Nullable

Like `eval`, but does it in the current context and without any dynamism.
Returns `Nullable()` if the result can't be evaluated.
"""
abstract_eval(ctx::LintContext, ex::Symbol) =
    flatten(BROADCAST(extractobject, lookup(ctx, ex)))

"""
    abstract_eval(ctx::LintContext, ex::Expr) :: Nullable

If the given expression is curly, and each component of the curly is a constant
object in the given `ctx`, construct the object `x` as would have been done in
the program itself, and return `Nullable(x)`.

Otherwise, if the given expression is `foo.bar`, and `foo` is a standard
library object with attribute `bar`, then construct `foo.bar` as would be done
in the program itself and return it.

Otherwise, return `Nullable()`.
"""
abstract_eval(ctx::LintContext, ex::Expr) = begin
    if isexpr(ex, :curly)
        # TODO: when 0.5 support dropped, remove [...] around ctx
        objs = abstract_eval.([ctx], ex.args)
        if all(!isnull, objs)
            try
                Nullable(Core.apply_type(get.(objs)...))
            catch
                Nullable()
            end
        else
            Nullable()
        end
    elseif isexpr(ex, :(.))
        head = ex.args[1]
        tail = ex.args[2].value
        obj = abstract_eval(ctx, head)
        if !isnull(obj)
            try
                Nullable(getfield(get(obj), tail))
            catch
                Nullable()
            end
        else
            Nullable()
        end
    else
        Nullable()
    end
end

"""
    abstract_eval(ctx::LintContext, ex)

Return the literal embedded within a `Nullable{Any}`.
"""
abstract_eval(ctx::LintContext, ex) = lexicalvalue(ex)


text:
const observables = [
    "Phosphorylated_Akt"
    "Phosphorylated_ERK"
    "Phosphorylated_cFos"
];

function observables_index(observable_name::String)::Int

    return findfirst(isequal(observable_name),observables)
end


text:
using Slack

"""
    sendtoslack(data::String, user_endpoint::String)
Takes in a string of data and the endpoint url and returns the HTTP response.
The response will be "ok" if the request is sent successfully.
"""
function sendtoslack(data::String, user_endpoint::String)
    HTTP.setuseragent!("HTTP.jl")

    params = Dict("text"=>"$(data)")
    base_url = "https://hooks.slack.com"

    endpoint = user_endpoint
    url = base_url * endpoint

    r = HTTP.post(url,
                 ["Content-Type" => "application/json"],
                 JSON.json(params))

    return (String(r.body))
end

"""
    sendattachmenttoslack(data::String, user_endpoint::String)
Takes in a string of data and the endpoint url and returns the HTTP response.
The response will be "ok" if the request is sent successfully.
"""
function sendattachmenttoslack(data, user_endpoint::String)
    HTTP.setuseragent!("HTTP.jl")

    params = data
    base_url = "https://hooks.slack.com"

    endpoint = user_endpoint
    url = base_url * endpoint

    r = HTTP.post(url,
                 ["Content-Type" => "application/json"],
                 JSON.json(params))

    return (String(r.body))
end

"""
    sendfiletoslack(filepath, message, channel, token)

Sends a file with a message to the specified channel.
https://api.slack.com/methods/files.upload
"""
function sendfiletoslack(filepath, message, channel, token)
    HTTP.setuseragent!("HTTP.jl")

    url = "https://slack.com/api/files.upload"
    headers = Dict("Authorization"=>"Bearer $token")
    data = Dict(
        "file" => open(filepath),
        "initial_comment" => message,
        "channels" => channel
    )
    body = HTTP.Form(data)

    r = HTTP.post(url, headers, body)

    return (String(r.body))
end


text:
ENV["CARGO_TARGET_DIR"] = @__DIR__

const rustprojname = "RustDylib"
const rustlibname = "rustdylib"
const juliapackage = "JuliaPackageWithRustDep"

const libname = Sys.iswindows() ? rustlibname : "lib" * rustlibname
# Windows .dlls do not have the "lib" prefix

function build_dylib()
    clean()

    run(Cmd(`cargo build --release`, dir=joinpath(@__DIR__, rustprojname)))

    release_dir = joinpath(@__DIR__, "release")
    dylib = dylib_filename()

    release_dylib_filepath = joinpath(release_dir, dylib)
    @assert isfile(release_dylib_filepath) "$release_dylib_filepath not found. Build may have failed."
    mv(release_dylib_filepath, joinpath(@__DIR__, dylib))
    rm(release_dir, recursive=true)

    write_deps_file(libname, dylib, juliapackage)
end

function dylib_filename()
    @static if Sys.isapple()
        "$libname.dylib"
    elseif Sys.islinux()
        "$libname.so"
    elseif Sys.iswindows()
        "$libname.dll"
    else
        error("Not supported: $(Sys.KERNEL)")
    end
end

function write_deps_file(libname, libfile, juliapackage)
    script = """
import Libdl

const $libname = joinpath(@__DIR__, "$libfile")

function check_deps()
    global $libname
    if !isfile($libname)
        error("\$$libname does not exist, Please re-run Pkg.build(\\"$juliapackage\\"), and restart Julia.")
    end

    if Libdl.dlopen_e($libname) == C_NULL
        error("\$$libname cannot be opened, Please re-run Pkg.build(\\"$juliapackage\\"), and restart Julia.")
    end
end
"""

    open(joinpath(@__DIR__, "deps.jl"), "w") do f
        write(f, script)
    end
end

function clean()
    deps_file = joinpath(@__DIR__, "deps.jl")
    isfile(deps_file) && rm(deps_file)

    release_dir = joinpath(@__DIR__, "release")
    isdir(release_dir) && rm(release_dir, recursive=true)

    dylib_file = joinpath(@__DIR__, dylib_filename())
    isfile(dylib_file) && rm(dylib_file)
end

build_dylib()


text:
function default_scaleselection(X::AbstractMatrix, ::ObsDim.Constant{1})
    collect(1:size(X, 2))
end

function default_scaleselection(X::AbstractMatrix, ::ObsDim.Constant{2})
    collect(1:size(X, 1))
end

function default_scaleselection(X::AbstractMatrix, ::ObsDim.Last)
    collect(1:size(X, 1))
end

function default_scaleselection(x::AbstractVector)
    collect(1:length(x))
end

function default_scaleselection(x::AbstractVector, ::ObsDim.Last)
    collect(1:length(x))
end

function default_scaleselection(x::AbstractVector, ::ObsDim.Constant{M}) where {M}
    collect(1:length(x))
end

function default_scaleselection(D::AbstractDataFrame)
    valid_columns(D)
end

function valid_columns(D::AbstractDataFrame)
    valid_colnames = Symbol[]
    for colname in names(D)
        if (eltype(D[colname]) <: Real) & !any(isna, D[colname])
            push!(valid_colnames, colname)
        else
            warn("Skipping \"$colname\" because it either contains NA or is not of type <: Real")
        end
    end
    valid_colnames
end

function valid_columns(D::AbstractDataFrame, colnames)
    valid_colnames = Symbol[]
    for colname in colnames
        if (eltype(D[colname]) <: Real) & !(any(isna, D[colname]))
            push!(valid_colnames, colname)
        else
            warn("Skipping \"$colname\" because it either contains NA or is not of type <: Real")
        end
    end
    valid_colnames
end


text:
import Base: log,exp,cos,sin,cosh,sinh
mutable struct Log <: Component
	x
	y
end
Log(ex::EX)=Log(ex,:e)
log(ex::Ex)=Log(ex,:e)
log(ex::Ex,p::EX)=Log(ex,p)
function simplify(l::Log)
	l=Log(simplify(l.x),simplify(l.y))
	if l.x==l.y
		return 1
	end
	if isa(l.x,Number)
		if l.y==:e
			return log(l.x)
		elseif isa(l.y,Number)
			return log(l.x,l.y)
		end
	end
	return l
end

mutable struct Exp <: Component
	x
end
exp(ex::Ex)=Exp(ex)
function simplify(l::Exp)
	x=simplify(l.x)
	if isa(x,Number)
		return exp(x)
	end
	return Exp(x)
end
mutable struct Cos <: Component
	x
end
cos(ex::Ex)=Cos(ex)
function simplify(l::Cos)
	x=simplify(l.x)
	if isa(x,Number)
		return cos(x)
	end
	return cos(x)
end
mutable struct Sin <: Component
	x
end
sin(ex::Ex)=Sin(ex)
function simplify(l::Sin)
	x=simplify(l.x)
	if isa(x,Number)
		return sin(x)
	end
	return Sin(x)
end
for c=((:Cosh,:cosh),(:Sinh,:sinh))
	eval(Meta.parse("mutable struct $(c[1])<:Component;x;end"))
	eval(Meta.parse("$(c[2])(ex::Ex)=$(c[1])(ex)"))
	eval(Meta.parse("function simplify(l::$(c[1]))
			x=simplify(l.x)
			if isa(x,Number)
				return $(c[2])(x)
			end
			return $(c[1])(x)
		end"))
end


text:
"""
    clusterPalette(n; alpha=1.0, h=0.0, hl=360.0, s::Array{Float64,1}=[0.7, 1.0], v=[1.0, 0.7])

Create an array of categorical colors suitable for plotting colored clusters.

# Arguments
- `n` number of colors
- `alpha` transparency
- `h`, `hs` starting hue and hue range
- `s`,`v` arrays of saturations and color values in the palette
   (the numbers will be used in a cyclic manner)
"""
function clusterPalette(n; alpha = 1.0, h = 0.0, hl = 360.0, s = [0.7, 1.0], v = [1.0, 0.7])
    [
        Colors.convert(
            Colors.RGBA,
            Colors.HSVA(h + hl * (i - 1) / n, s[1+i%length(s)], v[1+i%length(v)], alpha),
        ) for i = 1:n
    ]
end

"""
    linearPaletteBlend(c, r)

Internal function for picking the color blends from palette arrays
"""
function linearPaletteBlend(c, r)
    k = size(c, 1)
    kr = (k - 1) * r
    i = Int64(trunc(kr))

    if i < 0
        return c[1, :]
    end

    if i >= (k - 1)
        return c[k, :]
    end

    a = kr - i
    return (1 - a) * c[i+1, :] + a * c[i+2, :]
end

# data taken from a modified RbYlBu from EmbedSOM, which was originally taken
# from ColorBrewer and slightly modified
palRdYlBu = Matrix{Float64}(
    [
        33 38 149
        53 101 180
        100 157 209
        145 195 226
        184 217 200
        255 255 168
        254 224 144
        253 174 97
        244 109 67
        215 48 39
        165 0 38
    ],
)

"""
    expressionPalette(n; alpha=1.0, palette=palRdYlBu)

Create a palette of `n` colors with alpha `alpha` that is useful for displaying
expressions; by default a variation of the widely used RdYlBu palette is used.
Any supplied `palette` with colors ranging from 0 to `maxValue` can be expanded.
"""
function expressionPalette(
    n::Int64;
    alpha::Float64 = 1.0,
    palette::Matrix{Float64} = palRdYlBu,
    maxValue::Float64 = 255.0,
)
    [
        Colors.RGBA(x[1], x[2], x[3], alpha) for
        x in [linearPaletteBlend(palette, (i - 1) / (n - 1)) for i = 1:n] / maxValue
    ]
end

"""
    classColors(classes::Array{Int64, 1}, palette)::Matrix{Float64}

Take integer class IDs from `classes` and the corresponding colors in
`Palette`, and produce a corresponding color matrix usable in `rasterize`.

Optionally, `classCount` specifies the number of classes; by default taken as a
maximum of `classes` if `classCount` zero.
"""
function classColors(classes::Array{Int64,1}, palette)::Matrix{Float64}
    hcat([[palette[i].r, palette[i].g, palette[i].b, palette[i].alpha] for i in classes]...)
end

"""
    expressionColors(expressions::Array{Float64, 1}, palette)::Matrix{Float64}

Convert array of floating point `expressions` in range [0..1] to an array of
colors usable in `rasterize`, using the `palette` as a color lookup.
"""
function expressionColors(expressions::Array{Float64,1}, palette)::Matrix{Float64}
    nColors = length(palette)
    hcat(
        [
            [palette[c].r, palette[c].g, palette[c].b, palette[c].alpha] for
            c in Array{Int64,1}(1 .+ trunc.((nColors - 1) * expressions))
        ]...,
    )
end

"""
    scaleMinMax(expressions::Array{Float64, 1})::Array{Float64, 1}

Linearly scale `expressions` into the interval [0..1], usable in
`expressionColors`.
"""
function scaleMinMax(expressions::Array{Float64,1})::Array{Float64,1} #TODO quantiles?
    expressions .-= minimum(expressions)
    m = maximum(expressions)
    if m > 0
        expressions ./= m
    end
    expressions
end
"""
    scaleNorm(expressions::Array{Float64, 1})::Array{Float64, 1}

Scale `expressions` to quantile in their empirical normal distribution, i.e. to [0..1], usable in
`expressionColors`.
"""
function scaleNorm(expressions::Array{Float64,1})::Array{Float64,1}
    len = length(expressions)
    if len == 0
        return []
    end

    # compute the empirical distribution
    mu = sum(expressions) / len
    var = sum((expressions .- mu) .^ 2) / len

    if var == 0
        var = 1 #save the day
    end

    dist = Distributions.Normal(mu, sqrt(var))

    Distributions.cdf.(dist, expressions)
end


text:
int_rules_1_3_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.3.2*P(x)*Q(x)^p =#
    @apply_utils Antiderivative(sqrt(~v) / (~d + ~(e') * (~x) ^ 4), ~x) => With([a = Coeff(~v, ~x, 0), b = Coeff(~v, ~x, 2), c = Coeff(~v, ~x, 4)], (a / ~d) * Subst(Antiderivative(1 / ((1 - 2 * b * (~x) ^ 2) + (b ^ 2 - 4 * a * c) * (~x) ^ 4), ~x), ~x, ~x / sqrt(~v)) <-- EqQ(c * ~d + a * ~e, 0) && PosQ(a * c)) <-- FreeQ([~d, ~e], ~x) && PolyQ(~v, (~x) ^ 2, 2)
    @apply_utils Antiderivative(sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) / (~d + ~(e') * (~x) ^ 4), ~x) => With([q = sqrt((~b) ^ 2 - 4 * ~a * ~c)], ((-(~a) * sqrt(~b + q)) / (2 * sqrt(2) * Rt(-(~a) * ~c, 2) * ~d)) * atan((sqrt(~b + q) * ~x * ((~b - q) + 2 * ~c * (~x) ^ 2)) / (2 * sqrt(2) * Rt(-(~a) * ~c, 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) + ((~a * sqrt(-(~b) + q)) / (2 * sqrt(2) * Rt(-(~a) * ~c, 2) * ~d)) * atanh((sqrt(-(~b) + q) * ~x * (~b + q + 2 * ~c * (~x) ^ 2)) / (2 * sqrt(2) * Rt(-(~a) * ~c, 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~c * ~d + ~a * ~e, 0) && NegQ(~a * ~c))
    @apply_utils Antiderivative((~P) ^ ~p * (~Q) ^ ~(q'), ~x) => With([PP = Factor(ReplaceAll(~P, (~x->begin
                                        #= none:1 =#
                                        #= none:3 =#
                                        sqrt(~x)
                                    end)))], Antiderivative(ExpandIntegrand(ReplaceAll(PP, (~x->begin
                                                #= none:4 =#
                                                #= none:6 =#
                                                (~x) ^ 2
                                            end)) ^ ~p * (~Q) ^ ~q, ~x), ~x) <-- Not(SumQ(NonfreeFactors(PP, ~x)))) <-- FreeQ(~q, ~x) && (PolyQ(~P, (~x) ^ 2) && (PolyQ(~Q, ~x) && ILtQ(~p, 0)))
    @apply_utils Antiderivative((~P) ^ ~p * (~Q) ^ ~(q'), ~x) => With([PP = Factor(~P)], Antiderivative(ExpandIntegrand(PP ^ ~p * (~Q) ^ ~q, ~x), ~x) <-- Not(SumQ(NonfreeFactors(PP, ~x)))) <-- FreeQ(~q, ~x) && (PolyQ(~P, ~x) && (PolyQ(~Q, ~x) && (IntegerQ(~p) && NeQ(~P, ~x))))
    @apply_utils Antiderivative((~P) ^ ~p * ~Qm, ~x) => With([PP = Factor(~P)], Antiderivative(ExpandIntegrand(PP ^ ~p * ~Qm, ~x), ~x) <-- QuadraticProductQ(PP, ~x)) <-- PolyQ(~Qm, ~x) && (PolyQ(~P, ~x) && ILtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(d') * (~x) ^ 3) ^ ~(p'), ~x) => (1 / (3 ^ (3 * ~p) * (~a) ^ (2 * ~p))) * Antiderivative((~e + ~f * ~x) ^ ~m * (3 * ~a - ~b * ~x) ^ ~p * (3 * ~a + 2 * ~b * ~x) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m], ~x) && (EqQ(4 * (~b) ^ 3 + 27 * (~a) ^ 2 * ~d, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(d') * (~x) ^ 3) ^ ~p, ~x) => ((~a + ~b * ~x + ~d * (~x) ^ 3) ^ ~p / ((3 * ~a - ~b * ~x) ^ ~p * (3 * ~a + 2 * ~b * ~x) ^ (2 * ~p))) * Antiderivative((~e + ~f * ~x) ^ ~m * (3 * ~a - ~b * ~x) ^ ~p * (3 * ~a + 2 * ~b * ~x) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~p], ~x) && (EqQ(4 * (~b) ^ 3 + 27 * (~a) ^ 2 * ~d, 0) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(d') * (~x) ^ 3) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~e + ~f * ~x) ^ ~m * (~a + ~b * ~x + ~d * (~x) ^ 3) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m], ~x) && (NeQ(4 * (~b) ^ 3 + 27 * (~a) ^ 2 * ~d, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(d') * (~x) ^ 3) ^ ~p, ~x) => With([r = Rt(-9 * ~a * (~d) ^ 2 + sqrt(3) * ~d * sqrt(4 * (~b) ^ 3 * ~d + 27 * (~a) ^ 2 * (~d) ^ 2), 3)], (1 / (~d) ^ (2 * ~p)) * Antiderivative((~e + ~f * ~x) ^ ~m * Simp(((18 ^ (1 / 3) * ~b * ~d) / (3r) - r / 18 ^ (1 / 3)) + ~d * ~x, ~x) ^ ~p * Simp((((~b * ~d) / 3 + (12 ^ (1 / 3) * (~b) ^ 2 * (~d) ^ 2) / (3 * r ^ 2) + r ^ 2 / (3 * 12 ^ (1 / 3))) - ~d * ((2 ^ (1 / 3) * ~b * ~d) / (3 ^ (1 / 3) * r) - r / 18 ^ (1 / 3)) * ~x) + (~d) ^ 2 * (~x) ^ 2, ~x) ^ ~p, ~x)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m], ~x) && (NeQ(4 * (~b) ^ 3 + 27 * (~a) ^ 2 * ~d, 0) && ILtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(d') * (~x) ^ 3) ^ ~p, ~x) => With([r = Rt(-9 * ~a * (~d) ^ 2 + sqrt(3) * ~d * sqrt(4 * (~b) ^ 3 * ~d + 27 * (~a) ^ 2 * (~d) ^ 2), 3)], ((~a + ~b * ~x + ~d * (~x) ^ 3) ^ ~p / (Simp(((18 ^ (1 / 3) * ~b * ~d) / (3r) - r / 18 ^ (1 / 3)) + ~d * ~x, ~x) ^ ~p * Simp((((~b * ~d) / 3 + (12 ^ (1 / 3) * (~b) ^ 2 * (~d) ^ 2) / (3 * r ^ 2) + r ^ 2 / (3 * 12 ^ (1 / 3))) - ~d * ((2 ^ (1 / 3) * ~b * ~d) / (3 ^ (1 / 3) * r) - r / 18 ^ (1 / 3)) * ~x) + (~d) ^ 2 * (~x) ^ 2, ~x) ^ ~p)) * Antiderivative((~e + ~f * ~x) ^ ~m * Simp(((18 ^ (1 / 3) * ~b * ~d) / (3r) - r / 18 ^ (1 / 3)) + ~d * ~x, ~x) ^ ~p * Simp((((~b * ~d) / 3 + (12 ^ (1 / 3) * (~b) ^ 2 * (~d) ^ 2) / (3 * r ^ 2) + r ^ 2 / (3 * 12 ^ (1 / 3))) - ~d * ((2 ^ (1 / 3) * ~b * ~d) / (3 ^ (1 / 3) * r) - r / 18 ^ (1 / 3)) * ~x) + (~d) ^ 2 * (~x) ^ 2, ~x) ^ ~p, ~x)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~p], ~x) && (NeQ(4 * (~b) ^ 3 + 27 * (~a) ^ 2 * ~d, 0) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~P3) ^ ~(p'), ~x) => With([a = Coeff(~P3, ~x, 0), b = Coeff(~P3, ~x, 1), c = Coeff(~P3, ~x, 2), d = Coeff(~P3, ~x, 3)], Subst(Antiderivative(((3 * d * ~e - c * ~f) / (3d) + ~f * ~x) ^ ~m * Simp((((2 * c ^ 3 - 9 * b * c * d) + 27 * a * d ^ 2) / (27 * d ^ 2) - ((c ^ 2 - 3 * b * d) * ~x) / (3d)) + d * (~x) ^ 3, ~x) ^ ~p, ~x), ~x, ~x + c / (3d)) <-- NeQ(c, 0)) <-- FreeQ([~e, ~f, ~m, ~p], ~x) && PolyQ(~P3, ~x, 3)
    @apply_utils Antiderivative(~x / sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(e') * (~x) ^ 4), ~x) => With([Px = ((1 / 320) * (33 * (~b) ^ 2 * ~c + 6 * ~a * (~c) ^ 2 + 40 * (~a) ^ 2 * ~e) - (22 / 5) * ~a * ~c * ~e * (~x) ^ 2) + (22 / 15) * ~b * ~c * ~e * (~x) ^ 3 + (1 / 4) * ~e * (5 * (~c) ^ 2 + 4 * ~a * ~e) * (~x) ^ 4 + (4 / 3) * ~b * (~e) ^ 2 * (~x) ^ 5 + 2 * ~c * (~e) ^ 2 * (~x) ^ 6 + (~e) ^ 3 * (~x) ^ 8], (1 / (8 * Rt(~e, 2))) * log(Px + Dist(1 / (8 * Rt(~e, 2) * ~x), D(Px, ~x), ~x) * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2 + ~e * (~x) ^ 4))) <-- FreeQ([~a, ~b, ~c, ~e], ~x) && (EqQ(71 * (~c) ^ 2 + 100 * ~a * ~e, 0) && EqQ(1152 * (~c) ^ 3 - 125 * (~b) ^ 2 * ~e, 0))
    @apply_utils Antiderivative((~A + ~(B') * ~x) / sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(d') * (~x) ^ 3 + ~(e') * (~x) ^ 4), ~x) => ~B * Subst(Antiderivative(~x / sqrt((((((-3 * (~d) ^ 4 + 16 * ~c * (~d) ^ 2 * ~e) - 64 * ~b * ~d * (~e) ^ 2) + 256 * ~a * (~e) ^ 3) / (256 * (~e) ^ 3) + ((((~d) ^ 3 - 4 * ~c * ~d * ~e) + 8 * ~b * (~e) ^ 2) * ~x) / (8 * (~e) ^ 2)) - ((3 * (~d) ^ 2 - 8 * ~c * ~e) * (~x) ^ 2) / (8 * ~e)) + ~e * (~x) ^ 4), ~x), ~x, ~d / (4 * ~e) + ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (EqQ(~B * ~d - 4 * ~A * ~e, 0) && (EqQ(~d * ((141 * (~d) ^ 3 - 752 * ~c * ~d * ~e) - 400 * ~b * (~e) ^ 2) + 16 * (~e) ^ 2 * (71 * (~c) ^ 2 + 100 * ~a * ~e), 0) && EqQ(144 * (3 * (~d) ^ 2 - 8 * ~c * ~e) ^ 3 + 125 * (((~d) ^ 3 - 4 * ~c * ~d * ~e) + 8 * ~b * (~e) ^ 2) ^ 2, 0)))
    @apply_utils Antiderivative((~f + ~(g') * (~x) ^ 2) / ((~d + ~(e') * ~x + ~(d') * (~x) ^ 2) * sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(b') * (~x) ^ 3 + ~(a') * (~x) ^ 4)), ~x) => ((~a * ~f) / (~d * Rt((~a) ^ 2 * (2 * ~a - ~c), 2))) * atan((~a * ~b + ((4 * (~a) ^ 2 + (~b) ^ 2) - 2 * ~a * ~c) * ~x + ~a * ~b * (~x) ^ 2) / (2 * Rt((~a) ^ 2 * (2 * ~a - ~c), 2) * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2 + ~b * (~x) ^ 3 + ~a * (~x) ^ 4))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~b * ~d - ~a * ~e, 0) && (EqQ(~f + ~g, 0) && PosQ((~a) ^ 2 * (2 * ~a - ~c))))
    @apply_utils Antiderivative((~f + ~(g') * (~x) ^ 2) / ((~d + ~(e') * ~x + ~(d') * (~x) ^ 2) * sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(b') * (~x) ^ 3 + ~(a') * (~x) ^ 4)), ~x) => ((-(~a) * ~f) / (~d * Rt(-((~a) ^ 2) * (2 * ~a - ~c), 2))) * atanh((~a * ~b + ((4 * (~a) ^ 2 + (~b) ^ 2) - 2 * ~a * ~c) * ~x + ~a * ~b * (~x) ^ 2) / (2 * Rt(-((~a) ^ 2) * (2 * ~a - ~c), 2) * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2 + ~b * (~x) ^ 3 + ~a * (~x) ^ 4))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~b * ~d - ~a * ~e, 0) && (EqQ(~f + ~g, 0) && NegQ((~a) ^ 2 * (2 * ~a - ~c))))
    @apply_utils Antiderivative(~P3 / (~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(d') * (~x) ^ 3 + ~(e') * (~x) ^ 4), ~x) => With([q = sqrt((8 * (~a) ^ 2 + (~b) ^ 2) - 4 * ~a * ~c), A = Coeff(~P3, ~x, 0), B = Coeff(~P3, ~x, 1), C = Coeff(~P3, ~x, 2), D = Coeff(~P3, ~x, 3)], (1 / q) * Antiderivative(((~b * A - 2 * ~a * B) + 2 * ~a * D + A * q + ((2 * ~a * A - 2 * ~a * C) + ~b * D + D * q) * ~x) / (2 * ~a + (~b + q) * ~x + 2 * ~a * (~x) ^ 2), ~x) - (1 / q) * Antiderivative(((((~b * A - 2 * ~a * B) + 2 * ~a * D) - A * q) + (((2 * ~a * A - 2 * ~a * C) + ~b * D) - D * q) * ~x) / (2 * ~a + (~b - q) * ~x + 2 * ~a * (~x) ^ 2), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~P3, ~x, 3) && (EqQ(~a, ~e) && EqQ(~b, ~d)))
    @apply_utils Antiderivative(((~x) ^ ~(m') * ~P3) / (~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(d') * (~x) ^ 3 + ~(e') * (~x) ^ 4), ~x) => With([q = sqrt((8 * (~a) ^ 2 + (~b) ^ 2) - 4 * ~a * ~c), A = Coeff(~P3, ~x, 0), B = Coeff(~P3, ~x, 1), C = Coeff(~P3, ~x, 2), D = Coeff(~P3, ~x, 3)], (1 / q) * Antiderivative(((~x) ^ ~m * ((~b * A - 2 * ~a * B) + 2 * ~a * D + A * q + ((2 * ~a * A - 2 * ~a * C) + ~b * D + D * q) * ~x)) / (2 * ~a + (~b + q) * ~x + 2 * ~a * (~x) ^ 2), ~x) - (1 / q) * Antiderivative(((~x) ^ ~m * ((((~b * A - 2 * ~a * B) + 2 * ~a * D) - A * q) + (((2 * ~a * A - 2 * ~a * C) + ~b * D) - D * q) * ~x)) / (2 * ~a + (~b - q) * ~x + 2 * ~a * (~x) ^ 2), ~x)) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (PolyQ(~P3, ~x, 3) && (EqQ(~a, ~e) && EqQ(~b, ~d)))
    @apply_utils Antiderivative((~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2) / (~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(d') * (~x) ^ 3 + ~(e') * (~x) ^ 4), ~x) => With([q = Rt(~C * (2 * ~e * (~B * ~d - 4 * ~A * ~e) + ~C * ((~d) ^ 2 - 4 * ~c * ~e)), 2)], ((-2 * (~C) ^ 2) / q) * atanh(((~C * ~d - ~B * ~e) + 2 * ~C * ~e * ~x) / q) + ((2 * (~C) ^ 2) / q) * atanh((~C * (((4 * ~B * ~c * ~C - 3 * (~B) ^ 2 * ~d) - 4 * ~A * ~C * ~d) + 12 * ~A * ~B * ~e + 4 * ~C * ((2 * ~c * ~C - ~B * ~d) + 2 * ~A * ~e) * ~x + 4 * ~C * (2 * ~C * ~d - ~B * ~e) * (~x) ^ 2 + 8 * (~C) ^ 2 * ~e * (~x) ^ 3)) / (q * ((~B) ^ 2 - 4 * ~A * ~C)))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B, ~C], ~x) && (EqQ(((~B) ^ 2 * ~d + 2 * ~C * (~b * ~C + ~A * ~d)) - 2 * ~B * (~c * ~C + 2 * ~A * ~e), 0) && (EqQ((((2 * (~B) ^ 2 * ~c * ~C - 8 * ~a * (~C) ^ 3) - (~B) ^ 3 * ~d) - 4 * ~A * ~B * ~C * ~d) + 4 * ~A * ((~B) ^ 2 + 2 * ~A * ~C) * ~e, 0) && PosQ(~C * (2 * ~e * (~B * ~d - 4 * ~A * ~e) + ~C * ((~d) ^ 2 - 4 * ~c * ~e)))))
    @apply_utils Antiderivative((~(A') + ~(C') * (~x) ^ 2) / (~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(d') * (~x) ^ 3 + ~(e') * (~x) ^ 4), ~x) => With([q = Rt(~C * (-8 * ~A * (~e) ^ 2 + ~C * ((~d) ^ 2 - 4 * ~c * ~e)), 2)], ((-2 * (~C) ^ 2) / q) * atanh((~C * (~d + 2 * ~e * ~x)) / q) + ((2 * (~C) ^ 2) / q) * atanh((~C * (((~A * ~d - 2 * (~c * ~C + ~A * ~e) * ~x) - 2 * ~C * ~d * (~x) ^ 2) - 2 * ~C * ~e * (~x) ^ 3)) / (~A * q))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~C], ~x) && (EqQ(~b * ~C + ~A * ~d, 0) && (EqQ(~a * (~C) ^ 2 - (~A) ^ 2 * ~e, 0) && PosQ(~C * (-8 * ~A * (~e) ^ 2 + ~C * ((~d) ^ 2 - 4 * ~c * ~e)))))
    @apply_utils Antiderivative((~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2) / (~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(d') * (~x) ^ 3 + ~(e') * (~x) ^ 4), ~x) => With([q = Rt(-(~C) * (2 * ~e * (~B * ~d - 4 * ~A * ~e) + ~C * ((~d) ^ 2 - 4 * ~c * ~e)), 2)], ((2 * (~C) ^ 2) / q) * atan(((~C * ~d - ~B * ~e) + 2 * ~C * ~e * ~x) / q) - ((2 * (~C) ^ 2) / q) * atan((~C * (((4 * ~B * ~c * ~C - 3 * (~B) ^ 2 * ~d) - 4 * ~A * ~C * ~d) + 12 * ~A * ~B * ~e + 4 * ~C * ((2 * ~c * ~C - ~B * ~d) + 2 * ~A * ~e) * ~x + 4 * ~C * (2 * ~C * ~d - ~B * ~e) * (~x) ^ 2 + 8 * (~C) ^ 2 * ~e * (~x) ^ 3)) / (q * ((~B) ^ 2 - 4 * ~A * ~C)))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B, ~C], ~x) && (EqQ(((~B) ^ 2 * ~d + 2 * ~C * (~b * ~C + ~A * ~d)) - 2 * ~B * (~c * ~C + 2 * ~A * ~e), 0) && (EqQ((((2 * (~B) ^ 2 * ~c * ~C - 8 * ~a * (~C) ^ 3) - (~B) ^ 3 * ~d) - 4 * ~A * ~B * ~C * ~d) + 4 * ~A * ((~B) ^ 2 + 2 * ~A * ~C) * ~e, 0) && NegQ(~C * (2 * ~e * (~B * ~d - 4 * ~A * ~e) + ~C * ((~d) ^ 2 - 4 * ~c * ~e)))))
    @apply_utils Antiderivative((~(A') + ~(C') * (~x) ^ 2) / (~a + ~(b') * ~x + ~(c') * (~x) ^ 2 + ~(d') * (~x) ^ 3 + ~(e') * (~x) ^ 4), ~x) => With([q = Rt(-(~C) * (-8 * ~A * (~e) ^ 2 + ~C * ((~d) ^ 2 - 4 * ~c * ~e)), 2)], ((2 * (~C) ^ 2) / q) * atan((~C * ~d + 2 * ~C * ~e * ~x) / q) - ((2 * (~C) ^ 2) / q) * atan((-(~C) * (-(~A) * ~d + 2 * (~c * ~C + ~A * ~e) * ~x + 2 * ~C * ~d * (~x) ^ 2 + 2 * ~C * ~e * (~x) ^ 3)) / (~A * q))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~C], ~x) && (EqQ(~b * ~C + ~A * ~d, 0) && (EqQ(~a * (~C) ^ 2 - (~A) ^ 2 * ~e, 0) && NegQ(~C * (-8 * ~A * (~e) ^ 2 + ~C * ((~d) ^ 2 - 4 * ~c * ~e)))))
    @apply_utils Antiderivative(~Px * (~P4) ^ ~p, ~x) => With([a = Coeff(~P4, ~x, 0), b = Coeff(~P4, ~x, 1), c = Coeff(~P4, ~x, 2), d = Coeff(~P4, ~x, 3), e = Coeff(~P4, ~x, 4)], (1 / a ^ (3 * ~p)) * Antiderivative(ExpandIntegrand((~Px * (a - b * ~x) ^ -(~p)) / (a ^ 5 - b ^ 5 * (~x) ^ 5) ^ -(~p), ~x), ~x) <-- NeQ(a, 0) && (EqQ(c, b ^ 2 / a) && (EqQ(d, b ^ 3 / a ^ 2) && EqQ(e, b ^ 4 / a ^ 3)))) <-- FreeQ(~p, ~x) && (PolyQ(~P4, ~x, 4) && (PolyQ(~Px, ~x) && ILtQ(~p, 0)))
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~n) / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ ~n + ~(d') * (~x) ^ ~n2), ~x) => (~A) ^ 2 * (~n - 1) * Subst(Antiderivative(1 / (~a + (~A) ^ 2 * ~b * (~n - 1) ^ 2 * (~x) ^ 2), ~x), ~x, ~x / (~A * (~n - 1) - ~B * (~x) ^ ~n)) <-- FreeQ([~a, ~b, ~c, ~d, ~A, ~B, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~n, 2) && (EqQ(~a * (~B) ^ 2 - (~A) ^ 2 * ~d * (~n - 1) ^ 2, 0) && EqQ(~B * ~c + 2 * ~A * ~d * (~n - 1), 0))))
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~A + ~(B') * (~x) ^ ~(n'))) / (~a + ~(b') * (~x) ^ ~(k') + ~(c') * (~x) ^ ~(n') + ~(d') * (~x) ^ ~n2), ~x) => (((~A) ^ 2 * ((~m - ~n) + 1)) / (~m + 1)) * Subst(Antiderivative(1 / (~a + (~A) ^ 2 * ~b * ((~m - ~n) + 1) ^ 2 * (~x) ^ 2), ~x), ~x, (~x) ^ (~m + 1) / (~A * ((~m - ~n) + 1) + ~B * (~m + 1) * (~x) ^ ~n)) <-- FreeQ([~a, ~b, ~c, ~d, ~A, ~B, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~k, 2 * (~m + 1)) && (EqQ(~a * (~B) ^ 2 * (~m + 1) ^ 2 - (~A) ^ 2 * ~d * ((~m - ~n) + 1) ^ 2, 0) && EqQ(~B * ~c * (~m + 1) - 2 * ~A * ~d * ((~m - ~n) + 1), 0))))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) / (~d + ~(e') * (~x) ^ 2 + ~(f') * (~x) ^ 4 + ~(g') * (~x) ^ 6), ~x) => With([q = Rt((-(~a) * ~c * (~f) ^ 2 + 12 * (~a) ^ 2 * (~g) ^ 2 + ~f * (3 * (~c) ^ 2 * ~d - 2 * ~a * ~b * ~g)) / (~c * ~g * (3 * ~c * ~d - ~a * ~f)), 2), r = Rt(((~a * ~c * (~f) ^ 2 + 4 * ~g * (~b * ~c * ~d + (~a) ^ 2 * ~g)) - ~f * (3 * (~c) ^ 2 * ~d + 2 * ~a * ~b * ~g)) / (~c * ~g * (3 * ~c * ~d - ~a * ~f)), 2)], ((~c / (~g * q)) * atan((r + 2 * ~x) / q) - (~c / (~g * q)) * atan((r - 2 * ~x) / q)) - (~c / (~g * q)) * atan((((3 * ~c * ~d - ~a * ~f) * ~x) / (~g * q * (~b * ~c * ~d - 2 * (~a) ^ 2 * ~g) * ((~b * ~c * ~d - ~a * ~b * ~f) + 4 * (~a) ^ 2 * ~g))) * (((~b * (~c) ^ 2 * ~d * ~f - ~a * (~b) ^ 2 * ~f * ~g) - 2 * (~a) ^ 2 * ~c * ~f * ~g) + 6 * (~a) ^ 2 * ~b * (~g) ^ 2 + ~c * (((3 * (~c) ^ 2 * ~d * ~f - ~a * ~c * (~f) ^ 2) - ~b * ~c * ~d * ~g) + 2 * (~a) ^ 2 * (~g) ^ 2) * (~x) ^ 2 + (~c) ^ 2 * ~g * (3 * ~c * ~d - ~a * ~f) * (~x) ^ 4))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(((9 * (~c) ^ 3 * (~d) ^ 2 - ~c * ((~b) ^ 2 + 6 * ~a * ~c) * ~d * ~f) + (~a) ^ 2 * ~c * (~f) ^ 2 + 2 * ~a * ~b * (3 * ~c * ~d + ~a * ~f) * ~g) - 12 * (~a) ^ 3 * (~g) ^ 2, 0) && (EqQ(((3 * (~c) ^ 4 * (~d) ^ 2 * ~e - 3 * (~a) ^ 2 * (~c) ^ 2 * ~d * ~f * ~g) + (~a) ^ 3 * ~c * (~f) ^ 2 * ~g + 2 * (~a) ^ 3 * (~g) ^ 2 * (~b * ~f - 6 * ~a * ~g)) - (~c) ^ 3 * ~d * ((2 * ~b * ~d * ~f + ~a * ~e * ~f) - 12 * ~a * ~d * ~g), 0) && (NeQ(3 * ~c * ~d - ~a * ~f, 0) && (NeQ(~b * ~c * ~d - 2 * (~a) ^ 2 * ~g, 0) && (NeQ((~b * ~c * ~d - ~a * ~b * ~f) + 4 * (~a) ^ 2 * ~g, 0) && PosQ((-(~a) * ~c * (~f) ^ 2 + 12 * (~a) ^ 2 * (~g) ^ 2 + ~f * (3 * (~c) ^ 2 * ~d - 2 * ~a * ~b * ~g)) / (~c * ~g * (3 * ~c * ~d - ~a * ~f))))))))
    @apply_utils Antiderivative((~a + ~(c') * (~x) ^ 4) / (~d + ~(e') * (~x) ^ 2 + ~(f') * (~x) ^ 4 + ~(g') * (~x) ^ 6), ~x) => With([q = Rt((-(~a) * ~c * (~f) ^ 2 + 12 * (~a) ^ 2 * (~g) ^ 2 + 3 * ~f * (~c) ^ 2 * ~d) / (~c * ~g * (3 * ~c * ~d - ~a * ~f)), 2), r = Rt(((~a * ~c * (~f) ^ 2 + 4 * (~a) ^ 2 * (~g) ^ 2) - 3 * (~c) ^ 2 * ~d * ~f) / (~c * ~g * (3 * ~c * ~d - ~a * ~f)), 2)], ((~c / (~g * q)) * atan((r + 2 * ~x) / q) - (~c / (~g * q)) * atan((r - 2 * ~x) / q)) - (~c / (~g * q)) * atan((~c * (3 * ~c * ~d - ~a * ~f) * ~x * ((2 * (~a) ^ 2 * ~f * ~g - ((3 * (~c) ^ 2 * ~d * ~f - ~a * ~c * (~f) ^ 2) + 2 * (~a) ^ 2 * (~g) ^ 2) * (~x) ^ 2) - ~c * (3 * ~c * ~d - ~a * ~f) * ~g * (~x) ^ 4)) / (8 * (~a) ^ 4 * (~g) ^ 3 * q))) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(((9 * (~c) ^ 3 * (~d) ^ 2 - 6 * ~a * (~c) ^ 2 * ~d * ~f) + (~a) ^ 2 * ~c * (~f) ^ 2) - 12 * (~a) ^ 3 * (~g) ^ 2, 0) && (EqQ((((3 * (~c) ^ 4 * (~d) ^ 2 * ~e - 3 * (~a) ^ 2 * (~c) ^ 2 * ~d * ~f * ~g) + (~a) ^ 3 * ~c * (~f) ^ 2 * ~g) - 12 * (~a) ^ 4 * (~g) ^ 3) - ~a * (~c) ^ 3 * ~d * (~e * ~f - 12 * ~d * ~g), 0) && (NeQ(3 * ~c * ~d - ~a * ~f, 0) && PosQ((-(~a) * ~c * (~f) ^ 2 + 12 * (~a) ^ 2 * (~g) ^ 2 + 3 * (~c) ^ 2 * ~d * ~f) / (~c * ~g * (3 * ~c * ~d - ~a * ~f))))))
    @apply_utils Antiderivative(~u * (~Q6) ^ ~p, ~x) => With([a = Coeff(~Q6, ~x, 0), b = Coeff(~Q6, ~x, 2), c = Coeff(~Q6, ~x, 3), d = Coeff(~Q6, ~x, 4), e = Coeff(~Q6, ~x, 6)], (1 / (3 ^ (3 * ~p) * a ^ (2 * ~p))) * Antiderivative(ExpandIntegrand(~u * (3a + 3 * Rt(a, 3) ^ 2 * Rt(c, 3) * ~x + b * (~x) ^ 2) ^ ~p * ((3a - 3 * (-1) ^ (1 / 3) * Rt(a, 3) ^ 2 * Rt(c, 3) * ~x) + b * (~x) ^ 2) ^ ~p * (3a + 3 * (-1) ^ (2 / 3) * Rt(a, 3) ^ 2 * Rt(c, 3) * ~x + b * (~x) ^ 2) ^ ~p, ~x), ~x) <-- EqQ(b ^ 2 - 3 * a * d, 0) && EqQ(b ^ 3 - 27 * a ^ 2 * e, 0)) <-- ILtQ(~p, 0) && (PolyQ(~Q6, ~x, 6) && (EqQ(Coeff(~Q6, ~x, 1), 0) && (EqQ(Coeff(~Q6, ~x, 5), 0) && RationalFunctionQ(~u, ~x))))
    @apply_utils Antiderivative(~Pm / ~Qn, ~x) => With([m = Expon(~Pm, ~x), n = Expon(~Qn, ~x)], (Coeff(~Pm, ~x, m) * log(~Qn)) / (n * Coeff(~Qn, ~x, n)) + Simplify(~Pm - (Coeff(~Pm, ~x, m) * D(~Qn, ~x)) / (n * Coeff(~Qn, ~x, n))) * Antiderivative(1 / ~Qn, ~x) <-- EqQ(m, n - 1) && EqQ(D(Simplify(~Pm - (Coeff(~Pm, ~x, m) / (n * Coeff(~Qn, ~x, n))) * D(~Qn, ~x)), ~x), 0)) <-- PolyQ(~Pm, ~x) && PolyQ(~Qn, ~x)
    @apply_utils Antiderivative(~Pm * (~Qn) ^ ~p, ~x) => With([m = Expon(~Pm, ~x), n = Expon(~Qn, ~x)], (Coeff(~Pm, ~x, m) * (~Qn) ^ (~p + 1)) / (n * (~p + 1) * Coeff(~Qn, ~x, n)) + Simplify(~Pm - (Coeff(~Pm, ~x, m) * D(~Qn, ~x)) / (n * Coeff(~Qn, ~x, n))) * Antiderivative((~Qn) ^ ~p, ~x) <-- EqQ(m, n - 1) && EqQ(D(Simplify(~Pm - (Coeff(~Pm, ~x, m) / (n * Coeff(~Qn, ~x, n))) * D(~Qn, ~x)), ~x), 0)) <-- FreeQ(~p, ~x) && (PolyQ(~Pm, ~x) && (PolyQ(~Qn, ~x) && NeQ(~p, -1)))
    @apply_utils Antiderivative(~Pm / ~Qn, ~x) => With([m = Expon(~Pm, ~x), n = Expon(~Qn, ~x)], (Coeff(~Pm, ~x, m) * log(~Qn)) / (n * Coeff(~Qn, ~x, n)) + (1 / (n * Coeff(~Qn, ~x, n))) * Antiderivative(ExpandToSum(n * Coeff(~Qn, ~x, n) * ~Pm - Coeff(~Pm, ~x, m) * D(~Qn, ~x), ~x) / ~Qn, ~x) <-- EqQ(m, n - 1)) <-- PolyQ(~Pm, ~x) && PolyQ(~Qn, ~x)
    @apply_utils Antiderivative(~Pm * (~Qn) ^ ~p, ~x) => With([m = Expon(~Pm, ~x), n = Expon(~Qn, ~x)], (Coeff(~Pm, ~x, m) * (~Qn) ^ (~p + 1)) / (n * (~p + 1) * Coeff(~Qn, ~x, n)) + (1 / (n * Coeff(~Qn, ~x, n))) * Antiderivative(ExpandToSum(n * Coeff(~Qn, ~x, n) * ~Pm - Coeff(~Pm, ~x, m) * D(~Qn, ~x), ~x) * (~Qn) ^ ~p, ~x) <-- EqQ(m, n - 1)) <-- FreeQ(~p, ~x) && (PolyQ(~Pm, ~x) && (PolyQ(~Qn, ~x) && NeQ(~p, -1)))
    @apply_utils Antiderivative(~Pm * (~Qn) ^ ~(p'), ~x) => With([m = Expon(~Pm, ~x), n = Expon(~Qn, ~x)], (Coeff(~Pm, ~x, m) * (~x) ^ ((m - n) + 1) * (~Qn) ^ (~p + 1)) / ((m + n * ~p + 1) * Coeff(~Qn, ~x, n)) + (1 / ((m + n * ~p + 1) * Coeff(~Qn, ~x, n))) * Antiderivative(ExpandToSum((m + n * ~p + 1) * Coeff(~Qn, ~x, n) * ~Pm - Coeff(~Pm, ~x, m) * (~x) ^ (m - n) * (((m - n) + 1) * ~Qn + (~p + 1) * ~x * D(~Qn, ~x)), ~x) * (~Qn) ^ ~p, ~x) <-- LtQ(1, n, m + 1) && m + n * ~p + 1 < 0) <-- FreeQ(~p, ~x) && (PolyQ(~Pm, ~x) && (PolyQ(~Qn, ~x) && LtQ(~p, -1)))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.3 Miscellaneous/1.3.2 P(x) Q(x)^p.jl


text:
#=

Generate plots of value of employment and unemployment in the McCall model.

=#

using Plots, LaTeXStrings
pyplot()

mcm = McCallModel()
V, U = solve_mccall_model(mcm)
U_vec = U .* ones(length(mcm.w_vec))

plot(mcm.w_vec,
        [V U_vec],
        lw=2,
        alpha=0.7,
        label=[L"$V$" L"$U$"])


text:
# Plotting routines -- uses PyPlot

using PyPlot

# Bodemag with frsp and om given
function BodeMagPlot(freqResponse::Vector{ComplexF64},om::Vector{Float64})
    plot(om,abs.(freqResponse),linewidth=.75);
    xscale("log");
    yscale("log");
    xlabel("Frequency (rad/s)");
    ylabel("Magnitude");
    tight_layout();
end

# Bodemag directly from system
function BodeMagPlot(sys::StateSpace,no::Integer, ni::Integer, om::Vector{Float64}=Array{Float64}(undef,0))
    if (ni>sys.nu) error("Input index exceeds number of system inputs") end
    if (no>sys.ny) error("Output index exceeds number of system output") end

    if isempty(om)
        freqResponse,om = FrequencyResponse(sys,no,ni);
    else
        freqResponse = FrequencyResponse(sys,no,ni,om);
    end
    BodeMagPlot(freqResponse,om);
end

# Bode phase with frsp and om given -- this method should be sufficient. Other methods are not necessary.
function BodePhasePlot(freqResponse::Vector{ComplexF64},om::Vector{Float64})
    phase = 180*atan.(imag.(freqResponse),real.(freqResponse))/pi; # Convert it to [-180,180]
    plot(om,phase,linewidth=.75);
    xscale("log");
    #yscale("log");
    xlabel("Frequency (rad/s)");
    ylabel("Phase");
    tight_layout();
end

function BodePlot(freqResponse::Vector{ComplexF64},om::Vector{Float64})
    subplot(211); BodeMagPlot(freqResponse,om);
    subplot(212); BodePhasePlot(freqResponse,om);
end

# Bodemag directly from system
function BodePlot(sys::StateSpace,no::Integer, ni::Integer, om::Vector{Float64}=Array{Float64}(undef,0))
    if (ni>sys.nu) error("Input index exceeds number of system inputs") end
    if (no>sys.ny) error("Output index exceeds number of system output") end

    if isempty(om)
        freqResponse,om = FrequencyResponse(sys,no,ni);
    else
        freqResponse = FrequencyResponse(sys,no,ni,om);
    end
    BodePlot(freqResponse,om);
end

function NyquistPlot(sys::StateSpace,no::Integer, ni::Integer, om::Vector{Float64}=Array{Float64}(undef,0))
    if (ni>sys.nu) error("Input index exceeds number of system inputs") end
    if (no>sys.ny) error("Output index exceeds number of system output") end

    if isempty(om)
        f1,om = FrequencyResponse(sys,no,ni);
        f2 = FrequencyResponse(sys,no,ni,-om);
    else
        f1 = FrequencyResponse(sys,no,ni,om);
        f2 = FrequencyResponse(sys,no,ni,-om);
    end
    plot(real.(f1),imag.(f1),linewidth=0.5); # Positive ω
    plot(real.(f2),imag.(f2),linewidth=0.5); # Negative ω
    grid(which="both",color=[1,1,1]*0.9,linewidth=0.25)
    xlabel("Real");
    ylabel("Imaginary");
    tight_layout();
end

function PoleZeroMap()
    # Todo ...
end

max_stars_repo_name: isrlab/RobustControl

max_issues_repo_name: isrlab/RobustControl

max_forks_repo_name: isrlab/RobustControl


text:
using Oceananigans.Fields: XFaceField, YFaceField, ZFaceField, CenterField, PressureField
using Oceananigans.BuoyancyModels: BuoyancyField

has_buoyancy_tracer(model) = :b ∈ keys(model.tracers)

"""
    horizontally_averaged_velocities(model)

Returns a dictionary in which `:u` represents the horizontal average of
`model.velocities.u` and `:v` represents the horizontal average of
`model.velocities.v`.
"""
function horizontally_averaged_velocities(model)

    u, v, w = model.velocities

    averages = Dict(
                    :u => AveragedField(u, dims=(1, 2)),
                    :v => AveragedField(v, dims=(1, 2))
                   )

    return averages
end

"""
    horizontally_averaged_tracers(model)

Returns a dictionary containing horizontal averages of
every tracer in `model.tracers`. The dictionary keys
correspond to the names of the tracer.
"""
function horizontally_averaged_tracers(model)

    averages = Dict()

    for tracer in keys(model.tracers)
        c = getproperty(model.tracers, tracer)
        averages[tracer] = AveragedField(c, dims=(1, 2))
    end

    return averages
end

"""
    velocity_covariances(model; u_scratch = XFaceField(model.architecture, model.grid),
                                v_scratch = YFaceField(model.architecture, model.grid),
                                w_scratch = ZFaceField(model.architecture, model.grid),
                                c_scratch = CenterField(model.architecture, model.grid))

Returns a dictionary containing horizontal averages of the velocity covariances uᵢuⱼ.

The scratch fields `u_scratch`, `v_scratch`, `w_scratch` and `c_scratch` specify
scratch space for computations at `u`, `v`, `w`, and tracer ("`c`") locations, respectively.
"""
function velocity_covariances(model; u_scratch = XFaceField(model.architecture, model.grid),
                                     v_scratch = YFaceField(model.architecture, model.grid),
                                     w_scratch = ZFaceField(model.architecture, model.grid),
                                     c_scratch = CenterField(model.architecture, model.grid))

    u, v, w = model.velocities

    covariances = Dict(
                       :uu => AveragedField(u * u, dims=(1, 2)),
                       :vv => AveragedField(v * v, dims=(1, 2)),
                       :ww => AveragedField(w * w, dims=(1, 2)),
                       :uv => AveragedField(u * v, dims=(1, 2)),
                       :wv => AveragedField(w * v, dims=(1, 2)),
                       :wu => AveragedField(w * u, dims=(1, 2))
                       )

    return covariances
end

"""
    tracer_covariances(model; b = BuoyancyField(model),
                              u_scratch = XFaceField(model.architecture, model.grid),
                              v_scratch = YFaceField(model.architecture, model.grid),
                              w_scratch = ZFaceField(model.architecture, model.grid),
                              c_scratch = CenterField(model.architecture, model.grid))

Returns a dictionary containing horizontal averages of
tracer variance, tracer-velocity covariances, and tracer-buoyancy covariance.

The scratch fields `u_scratch`, `v_scratch`, `w_scratch` and `c_scratch` specify
scratch space for computations at `u`, `v`, `w`, and tracer ("`c`") locations, respectively.

`b` is the `BuoyancyField(model)` for model. If `b` is not specified, new memory is allocated to
store the computation of buoyancy in three-dimensions over `model.grid`.
"""
function tracer_covariances(model; b = BuoyancyField(model),
                                   u_scratch = XFaceField(model.architecture, model.grid),
                                   v_scratch = YFaceField(model.architecture, model.grid),
                                   w_scratch = ZFaceField(model.architecture, model.grid),
                                   c_scratch = CenterField(model.architecture, model.grid))

    u, v, w = model.velocities

    covariances = Dict()

    for tracer in keys(model.tracers)
        c = getproperty(model.tracers, tracer)

        # Keys
        cc = Symbol(tracer, tracer)
        uc = Symbol(:u, tracer)
        vc = Symbol(:v, tracer)
        wc = Symbol(:w, tracer)

        covariances[cc] = AveragedField(c * c, dims=(1, 2))
        covariances[uc] = AveragedField(u * c, dims=(1, 2))
        covariances[vc] = AveragedField(v * c, dims=(1, 2))
        covariances[wc] = AveragedField(w * c, dims=(1, 2))

        # Add covariance of tracer with buoyancy
        if tracer != :b && !has_buoyancy_tracer(model)
            bc = Symbol(:b, tracer)
            covariances[bc] = AveragedField(b * c, dims=(1, 2))
        end
    end

    return covariances
end

"""
    third_order_velocity_statistics(model, u_scratch = XFaceField(model.architecture, model.grid),
                                           v_scratch = YFaceField(model.architecture, model.grid),
                                           w_scratch = ZFaceField(model.architecture, model.grid),
                                           c_scratch = CenterField(model.architecture, model.grid),
                                           p = PressureField(model))

Returns a dictionary of functions that calculate horizontally-averaged third-order statistics
that involve the velocity field. Includes statistics associated with pressure.

The scratch fields `u_scratch`, `v_scratch`, `w_scratch` and `c_scratch` specify
scratch space for computations at `u`, `v`, `w`, and tracer ("`c`") locations, respectively.
"""
function third_order_velocity_statistics(model; u_scratch = XFaceField(model.architecture, model.grid),
                                                v_scratch = YFaceField(model.architecture, model.grid),
                                                w_scratch = ZFaceField(model.architecture, model.grid),
                                                c_scratch = CenterField(model.architecture, model.grid),
                                                p = PressureField(model))

    u, v, w = model.velocities
    T,S = model.tracers
    # Pressure-strain terms
    Σˣʸ = 0.5 * (∂y(u) + ∂x(v)) # FFC
    Σˣᶻ = 0.5 * (∂z(u) + ∂x(w)) # FCF
    Σʸᶻ = 0.5 * (∂z(v) + ∂y(w)) # CFF

    third_order_statistics = Dict(
                                  :uuu => AveragedField(u * u * u, dims=(1, 2)),
                                  :uuv => AveragedField(u * u * v, dims=(1, 2)),
                                  :uuw => AveragedField(u * u * w, dims=(1, 2)),
                                  :uvv => AveragedField(u * v * v, dims=(1, 2)),
                                  :uww => AveragedField(u * w * w, dims=(1, 2)),

                                  :vvv => AveragedField(v * v * v, dims=(1, 2)),
                                  :vvw => AveragedField(v * v * w, dims=(1, 2)),
                                  :vww => AveragedField(v * w * w, dims=(1, 2)),

                                  :www => AveragedField(w * w * w, dims=(1, 2)),
                                  :wwt => AveragedField(w * w * T, dims=(1, 2)),
                                  :wws => AveragedField(w * w * S, dims=(1, 2)),
                                  :wvu => AveragedField(w * v * u, dims=(1, 2)),

                                  :up => AveragedField(u * p,     dims=(1, 2)),
                                  :vp => AveragedField(v * p,     dims=(1, 2)),
                                  :wp => AveragedField(w * p,     dims=(1, 2)),

                                  :pux => AveragedField(p * ∂x(u), dims=(1, 2)),
                                  :puy => AveragedField(p * ∂y(u), dims=(1, 2)),
                                  :puz => AveragedField(p * ∂z(u), dims=(1, 2)),

                                  :pvx => AveragedField(p * ∂x(v), dims=(1, 2)),
                                  :pvy => AveragedField(p * ∂y(v), dims=(1, 2)),
                                  :pvz => AveragedField(p * ∂z(v), dims=(1, 2)),

                                  :pwx => AveragedField(p * ∂x(w), dims=(1, 2)),
                                  :pwy => AveragedField(p * ∂y(w), dims=(1, 2)),
                                  :pwz => AveragedField(p * ∂z(w), dims=(1, 2)),
                                 )

    return third_order_statistics
end

"""
    third_order_tracer_statistics(model; u_scratch = XFaceField(model.architecture, model.grid),
                                         v_scratch = YFaceField(model.architecture, model.grid),
                                         w_scratch = ZFaceField(model.architecture, model.grid),
                                         c_scratch = CenterField(model.architecture, model.grid))

Returns a dictionary of functions that calculate horizontally-averaged third-order statistics
that involve tracers.

The scratch fields `u_scratch`, `v_scratch`, `w_scratch` and `c_scratch` specify
scratch space for computations at `u`, `v`, `w`, and tracer ("`c`") locations, respectively.
"""
function third_order_tracer_statistics(model; u_scratch = XFaceField(model.architecture, model.grid),
                                              v_scratch = YFaceField(model.architecture, model.grid),
                                              w_scratch = ZFaceField(model.architecture, model.grid),
                                              c_scratch = CenterField(model.architecture, model.grid),
                                              p = PressureField(model))

    u, v, w = model.velocities

    third_order_statistics = Dict()

    for tracer in keys(model.tracers)
        c = getproperty(model.tracers, tracer)

        # Keys
        wcc = Symbol(:w, tracer, tracer)
        cwu = Symbol(tracer, :wu)
        cpx = Symbol(tracer, :px)
        cpy = Symbol(tracer, :py)
        cpz = Symbol(tracer, :pz)

        third_order_statistics[wcc] = AveragedField(w * c * c, dims=(1, 2))
        third_order_statistics[cwu] = AveragedField(c * w * u, dims=(1, 2))
        third_order_statistics[cpx] = AveragedField(c * ∂x(p), dims=(1, 2))
        third_order_statistics[cpy] = AveragedField(c * ∂y(p), dims=(1, 2))
        third_order_statistics[cpz] = AveragedField(c * ∂z(p), dims=(1, 2))
    end

    return third_order_statistics
end

function first_order_statistics(model; b = BuoyancyField(model),
                                       u_scratch = XFaceField(model.architecture, model.grid),
                                       v_scratch = YFaceField(model.architecture, model.grid),
                                       w_scratch = ZFaceField(model.architecture, model.grid),
                                       c_scratch = CenterField(model.architecture, model.grid),
                                       p = PressureField(model))

    output = merge(
                   horizontally_averaged_velocities(model),
                   horizontally_averaged_tracers(model),
                   )

    output[:p] = AveragedField(p, dims=(1, 2))

    if !has_buoyancy_tracer(model)
        output[:b] = AveragedField(b, dims=(1, 2))
    end

    return output
end

function second_order_statistics(model; b = BuoyancyField(model),
                                        u_scratch = XFaceField(model.architecture, model.grid),
                                        v_scratch = YFaceField(model.architecture, model.grid),
                                        w_scratch = ZFaceField(model.architecture, model.grid),
                                        c_scratch = CenterField(model.architecture, model.grid))

    output = merge(
                   velocity_covariances(model, u_scratch = u_scratch,
                                               v_scratch = v_scratch,
                                               w_scratch = w_scratch,
                                               c_scratch = c_scratch),

                   tracer_covariances(model, b = b,
                                             u_scratch = u_scratch,
                                             v_scratch = v_scratch,
                                             w_scratch = w_scratch,
                                             c_scratch = c_scratch)

                  )

    if !has_buoyancy_tracer(model)
        u, v, w = model.velocities

        output[:ub] = AveragedField(u * b, dims=(1, 2))
        output[:vb] = AveragedField(v * b, dims=(1, 2))
        output[:wb] = AveragedField(w * b, dims=(1, 2))
    end

    return output
end

function third_order_statistics(model; u_scratch = XFaceField(model.architecture, model.grid),
                                       v_scratch = YFaceField(model.architecture, model.grid),
                                       w_scratch = ZFaceField(model.architecture, model.grid),
                                       c_scratch = CenterField(model.architecture, model.grid),
                                       p = PressureField(model))

    output = merge(
                   third_order_velocity_statistics(model, p = p,
                                                   u_scratch = u_scratch,
                                                   v_scratch = v_scratch,
                                                   w_scratch = w_scratch,
                                                   c_scratch = c_scratch),

                   third_order_tracer_statistics(model, p = p,
                                                 u_scratch = u_scratch,
                                                 v_scratch = v_scratch,
                                                 w_scratch = w_scratch,
                                                 c_scratch = c_scratch)
                   )

    return output
end

function first_through_second_order(model; b = BuoyancyField(model), p = PressureField(model),
                                           u_scratch = XFaceField(model.architecture, model.grid),
                                           v_scratch = YFaceField(model.architecture, model.grid),
                                           w_scratch = ZFaceField(model.architecture, model.grid),
                                           c_scratch = CenterField(model.architecture, model.grid))

    output = merge(
                   first_order_statistics(model, b = b, p = p,
                                          u_scratch = u_scratch,
                                          v_scratch = v_scratch,
                                          w_scratch = w_scratch,
                                          c_scratch = c_scratch),

                   second_order_statistics(model,
                                           u_scratch = u_scratch,
                                           v_scratch = v_scratch,
                                           w_scratch = w_scratch,
                                           c_scratch = c_scratch,
                                                   b = b),

                  )

    return output
end

function first_through_third_order(model; b = BuoyancyField(model), p = PressureField(model),
                                          u_scratch = XFaceField(model.architecture, model.grid),
                                          v_scratch = YFaceField(model.architecture, model.grid),
                                          w_scratch = ZFaceField(model.architecture, model.grid),
                                          c_scratch = CenterField(model.architecture, model.grid))

    output = merge(
                   first_order_statistics(model, b = b, p = p,
                                          u_scratch = u_scratch,
                                          v_scratch = v_scratch,
                                          w_scratch = w_scratch,
                                          c_scratch = c_scratch),

                   second_order_statistics(model,
                                           u_scratch = u_scratch,
                                           v_scratch = v_scratch,
                                           w_scratch = w_scratch,
                                           c_scratch = c_scratch,
                                                   b = b),

                   third_order_statistics(model, p = p,
                                          u_scratch = u_scratch,
                                          v_scratch = v_scratch,
                                          w_scratch = w_scratch,
                                          c_scratch = c_scratch)

                  )

    return output
end


text:
"""
Copyright (C) 2019 Mo Zhou <lumin@debian.org>
MIT/Expat license.
"""
module DistroHelper

const JL_REGISTRY_URL = "https://github.com/JuliaRegistries/General/raw/master/"

using Pkg

"""
Return the slug string for the given package in specific version.
The latest version will be used when version is omitted.

```
slug("Package.toml", "Versions.toml", "0.20.0")
slug("Package.toml", "Versions.toml")
```
"""
function slug(package::AbstractString, versions::AbstractString, version::Any = nothing)
	if !any(map(isfile, (package, versions)))
		uuid = Base.UUID(package)
		sha1 = Base.SHA1(hex2bytes(versions))
	else
		Package  = Pkg.TOML.parsefile(package)
		Versions = Pkg.TOML.parsefile(versions)
		latest_ver = maximum(VersionNumber.(keys(Versions)))
		version = string((version == nothing) ? latest_ver : version)
		uuid = Base.UUID(Package["uuid"])
		sha1 = Base.SHA1(hex2bytes(Versions[version]["git-tree-sha1"]))
	end
	return Base.version_slug(uuid,sha1)
end

"""
Return the dependencies for a specific version

```
depends("Deps.toml", "0.20.0")
depends("Deps.toml")
```
"""
function depends(deps::AbstractString, version::Any = nothing)
	alldeps = Pkg.TOML.parsefile(deps)
	bounds = map(x -> VersionNumber.(x), split.(keys(alldeps), "-"))
	ver	= (version == nothing) ? maximum(maximum(bounds)) : VersionNumber(version)
	mask = [ver >= x[1] && ver <= x[2] for x in bounds]
	realdeps = Dict{String,Any}()
	for (k, v) in zip(mask, values(alldeps))
		if k for (i, j) in v realdeps[i] = j end end
	end
	return realdeps
end

"""
Gather manifest for stdlib
"""
function gather_stdlib_manifest()
	p = "/usr/share/julia/stdlib/v1.0"
	stdlibs = values(Pkg.Types.gather_stdlib_uuids())
	manifest = Dict{String,Any}()
	for lib in stdlibs
		d = Dict{String,Any}()
		project = Pkg.TOML.parsefile(joinpath(p, lib, "Project.toml"))
		name, uuid = project["name"], project["uuid"]
		d["uuid"] = uuid
		if haskey(project, "deps")
			d["deps"] = sort(string.(keys(project["deps"])))
		end
		manifest[name] = [d]
	end
	#Pkg.TOML.print(manifest)
	return manifest
end

"""
Generate Project.toml based on the given METADATA/Registry files.
Output will be written to pwd if dest is not specified.
"""
function distro_project(package::AbstractString, dest::Any = nothing)
	package = Pkg.TOML.parsefile(package)
	name, uuid = package["name"], package["uuid"]
	project = Dict{String,Any}()
	deps = Dict{String,Any}()
	deps[name] = uuid
	project["deps"] = deps
	open((dest==nothing) ? "dh_Project.toml" : dest, "w") do io
		Pkg.TOML.print(io, project)
	end
end

"""
Generate Manifest.toml entry for a specific package
"""
function distro_manifest_(package::AbstractString, versions::AbstractString, deps::AbstractString)
	d = Dict{String,Any}()
	package = Pkg.TOML.parsefile(package)
	versions = Pkg.TOML.parsefile(versions)
	latest_ver = string(maximum(VersionNumber.(keys(versions))))
	deps = depends(deps, latest_ver)  # XXX: version?
	d["deps"] = sort(collect(keys(deps)))
	d["git-tree-sha1"] = versions[latest_ver]["git-tree-sha1"]
	d["uuid"] = package["uuid"]
	d["version"] = latest_ver
	manifest = Dict{String,Any}()
	manifest[package["name"]] = [d]
	return manifest
end

"""
Generate Manifest.toml based on the given METADATA
FIXME: WIP
"""
function distro_manifest(package::AbstractString, versions::AbstractString, deps::AbstractString, dest::Any=nothing)
	manifest = gather_stdlib_manifest()
	merge!(manifest, distro_manifest_(package, versions, deps))
	open((dest==nothing) ? "dh_Manifest.toml" : dest, "w") do io
		Pkg.TOML.print(io, manifest)
	end
end

"""
Merge different Project.toml and Manifest.toml files
"""
function distro_merge(p::AbstractString, dest::Any=nothing)
	packages = filter(x -> isdir(joinpath(p, x)), readdir(p))
	println("Found $(length(packages)) packages: $packages")
	env = nothing
	for package in packages
		e = Pkg.Types.EnvCache(joinpath(p, package, "Project.toml"))
		if nothing == env
			env = e
		else
			merge!(env.project["deps"], e.project["deps"])
			merge!(env.manifest, e.manifest)
		end
	end
	mkpath(dest)
	open(joinpath(nothing==dest ? "." : dest, "Project.toml"), "w") do io
		Pkg.TOML.print(io, env.project)
	end
	open(joinpath(nothing==dest ? "." : dest, "Manifest.toml"), "w") do io
		Pkg.TOML.print(io, env.manifest)
	end
	#Pkg.Types.write_env(Pkg.Types.Context(env=env))
end

"""
Download metadata/registry of the specified package
"""
function metadata(name::AbstractString)
	baseurl = JL_REGISTRY_URL * "$(uppercase(first(name)))/$(name)"
	run(`wget -nv -c $(baseurl)/Compat.toml`)
	run(`wget -nv -c $(baseurl)/Deps.toml`)
	run(`wget -nv -c $(baseurl)/Package.toml`)
	run(`wget -nv -c $(baseurl)/Versions.toml`)
end

"""
Download not only the metadata, but also the source tarball.
https://github.com/JuliaIO/JSON.jl/archive/v0.20.0.tar.gz
"""
function download(name::AbstractString, version::Any=nothing)
	metadata(name)
	repo = Pkg.TOML.parsefile("Package.toml")["repo"]
	repo = replace(repo, r"\.git$" => "")
	if version == nothing
		versions = Pkg.TOML.parsefile("Versions.toml")
		version = string(maximum(VersionNumber.(keys(versions))))
	end
	dest = lowercase("julia-$(name)_$(version).orig.tar.gz")
	run(`wget -nv -c $repo/archive/v$version.tar.gz -O $dest`)
end

end # module


text:
using GLPlot, GLAbstraction
f = frag"""
float function(float x) {
    return sin(x*x*x)*sin(x);
}
"""

glplot(f)


text:
## CPU

function _scattering(p::NTuple{2,A}, xi::NTuple{2,B}, F::Array{C,2}) where {A,B,C}
    # unpack tuples
    x, y = xi
    dx, dy = x[2] - x[1], y[2] - y[1]
    px, py = p

    # check that the particle is inside the grid
    isinside(px, py, x, y)

    # indices of lowermost-left corner of the cell
    # containing the particle
    idx_x, idx_y = parent_cell(p, dxi)

    # normalize particle coordinates
    tx = (px - x[idx_x]) / dx
    ty = (py - y[idx_y]) / dy

    # Interpolate field F onto particle
    Fp = bilinear(
        tx,
        ty,
        F[idx_x, idx_y],
        F[idx_x + 1, idx_y],
        F[idx_x, idx_y + 1],
        F[idx_x + 1, idx_y + 1],
    )

    return Fp
end

function scattering(xi, F::Array{T,2}, particle_coords) where {T}
    np = length(particle_coords[1])
    Fp = zeros(T, np)

    Threads.@threads for i in 1:np
        @inbounds Fp[i] = _scattering((particle_coords[1][i], particle_coords[2][i]), xi, F)
    end

    return Fp
end

##  CUDA

function _scattering!(
    Fp::CuDeviceVector{T,1}, p::NTuple{2,A}, xi::NTuple{2,B}, F::CuDeviceMatrix{T,1}
) where {A,B,T}
    ix = (blockIdx().x - 1) * blockDim().x + threadIdx().x

    # unpack tuples
    dx, dy = x[2] - x[1], y[2] - y[1]
    px, py = p
    x, y = xi

    @inbounds if ix ≤ length(px)
        # check that the particle is inside the grid
        isinside(px[ix], py[ix], x, y)

        # indices of lowermost-left corner of the cell
        # containing the particle
        idx_x, idx_y = parent_cell((px[ix], py[ix]), dxi)

        # normalize particle coordinates
        tx = (px[ix] - x[idx_x]) / dx
        ty = (py[ix] - y[idx_y]) / dy

        # Interpolate field F onto particle
        Fp[ix] = bilinear(
            tx,
            ty,
            F[idx_x, idx_y],
            F[idx_x + 1, idx_y],
            F[idx_x, idx_y + 1],
            F[idx_x + 1, idx_y + 1],
        )
    end

    return nothing
end

function scattering(
    xi, Fd::CuArray{T,2}, particle_coords::NTuple{2,CuArray}; nt=512
) where {T}
    N = length(particle_coords[1])
    Fpd = CuArray{T}(undef, N)

    numblocks = ceil(Int, N / nt)
    CUDA.@sync begin
        @cuda threads = nt blocks = numblocks _scattering!(Fpd, particle_coords, xi, Fd)
    end
end

function scattering!(
    Fpd, xi, Fd::CuArray{T,2}, particle_coords::NTuple{2,CuArray}; nt=512
) where {T}
    N = length(particle_coords[1])
    numblocks = ceil(Int, N / nt)
    CUDA.@sync begin
        @cuda threads = nt blocks = numblocks _scattering!(Fpd, particle_coords, xi, Fd)
    end
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to simplified dipole model of the Earth geomagnetic field.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] http://helios.fmi.fi/~juusolal/geomagnetism/Lectures/Chapter3_dipole.pdf
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/earth/geomagnetic_field_models/dipole
# =================================================

# Function: geomag_dipole
# -----------------------

@testset "Function geomag_dipole" begin

    d2r = π/180

    # Test 1
    # --------------------------------------------------------------------------
    #
    # Position aligned with the dipole moment vector.

    year     = SatelliteToolbox._dipole_pole[32,1]
    pole_lat = SatelliteToolbox._dipole_pole[32,2] * d2r
    pole_lon = SatelliteToolbox._dipole_pole[32,3] * d2r
    m        = SatelliteToolbox._dipole_mag[32,2]  * 1e22

    # Distance from the Earth center.
    r = R0 + 196e3

    # Compute using the simplified equations [1] [nT].
    k0     = 1e-7m
    B_norm = 2k0 / r^3 * 1e9

    # Rotate to ECEF.
    Deg = angle_to_dcm(-(π / 2 - pole_lat), -pole_lon, 0, :YZX)
    B_g = SVector{3}(0, 0, -B_norm)
    B_e_expected = Deg * B_g

    # Compute using the function.
    r_g = SVector{3}(0, 0, r)
    r_e = Deg * r_g

    B_e_result_f64 = geomag_dipole(r_e, year)
    B_e_result_f32 = geomag_dipole(Float32.(r_e), year)

    @test B_e_expected ≈ B_e_result_f64 atol=1e-9
    @test eltype(B_e_result_f64) === Float64

    @test B_e_expected ≈ B_e_result_f32 atol=1e-1
    @test eltype(B_e_result_f32) === Float32

    # Test 2
    # --------------------------------------------------------------------------
    #
    # Position at the magnetic Equator, which must have half the magnitude of
    # that of the test 1.

    B_g = SVector{3}(0, 0, B_norm / 2)
    B_e_expected = Deg * B_g

    r_g = SVector{3}(r, 0, 0)
    r_e = Deg * r_g

    B_e_result_f64 = geomag_dipole(r_e, year)
    B_e_result_f32 = geomag_dipole(Float32.(r_e), year)

    @test B_e_expected ≈ B_e_result_f64 atol=1e-9
    @test eltype(B_e_result_f64) === Float64

    @test B_e_expected ≈ B_e_result_f32 atol=1e-1
    @test eltype(B_e_result_f32) === Float32
end


text:
module CIFAR10

import ..Metalhead
import ..Metalhead: testimgs
import ..Metalhead: ValidationImage, TrainingImage, ValData, TestData, TrainData, ObjectClass, labels

abstract type DataSet <: Metalhead.DataSet end
using Images
using ColorTypes
using ColorTypes: N0f8

const C10Labels = ["airplane", "automobile", "bird", "cat", "deer", "dog",
    "frog", "horse", "ship", "truck"]

struct C10Class <: ObjectClass
    class::Int
end
labels(::Type{C10Class}) = C10Labels

struct BinPackedFS <: DataSet
    folder::String
end

Base.size(v::ValData{<:DataSet}) = (10000,)

Base.size(v::TrainData{<:DataSet}) = (50000,)

testimgs(::DataSet) = error("CIFAR10 does not specify a test set (test_batch is considered the validation set)")

function bytes_to_image(bytes::Vector{UInt8})
    channelwise = reshape(bytes, (1024, 3))
    transpose(reshape(
        RGB.((reinterpret.(N0f8, channelwise[:, i]) for i = 1:3)...),
        (32, 32)))
end

function Base.getindex(v::ValData{BinPackedFS}, i::Integer)
    ValidationImage(DataSet, i, open(joinpath(v.set.folder, "test_batch.bin")) do f
        seek(f, (i-1)*3073)
        label = read(f, UInt8)
        bytes = read(f, 3072)
        bytes_to_image(bytes), C10Class(label+1)
    end...)
end

function Base.getindex(v::TrainData{BinPackedFS}, i::Integer)
    batch, num = divrem(i - 1, 10000)
    file = "data_batch_$(batch+1).bin"
    TrainingImage(DataSet, i, open(joinpath(v.set.folder, file)) do f
        seek(f, num * 3073)
        label = read(f, UInt8)
        bytes = read(f, 3072)
        bytes_to_image(bytes), C10Class(label+1)
    end...)
end

end


text:
using CUDA
using MPI
using Test
using LinearAlgebra
using Random
using StaticArrays
using KernelAbstractions: CPU, CUDADevice
using ClimateMachine
using ClimateMachine.SystemSolvers
using ClimateMachine.MPIStateArrays
using Random

Random.seed!(1235)

let
    ClimateMachine.init()
    mpicomm = MPI.COMM_WORLD
    ArrayType = ClimateMachine.array_type()
    n = 100
    T = Float64
    A = rand(n, n)
    scale = 2.0
    ϵ = 0.1
    # Matrix 1
    A = A' * A .* ϵ + scale * I

    err_thresh = sqrt(eps(T))

    # Matrix 2
    # A = Diagonal(collect(1:n) * 1.0)
    positive_definite = minimum(eigvals(A)) > eps(1.0)
    @test positive_definite

    b = ones(n) * 1.0
    mulbyA!(y, x) = (y .= A * x)

    tol = sqrt(eps(T))
    method(b, tol) = ConjugateGradient(b, max_iter = n)
    linearsolver = method(b, tol)

    x = ones(n) * 1.0
    x0 = copy(x)
    iters = linearsolve!(mulbyA!, nothing, linearsolver, x, b; max_iters = Inf)
    exact = A \ b
    x0 = copy(x)

    @testset "Array test" begin
        @test norm(x - exact) / norm(exact) < err_thresh
        @test norm(A * x - b) / norm(b) < err_thresh
    end

    # Testing for CUDA CuArrays
    if CUDA.has_cuda_gpu()
        at_A = ArrayType(A)
        at_b = ArrayType(b)
        at_x = ArrayType(ones(n) * 1.0)
        mulbyat_A!(y, x) = (y .= at_A * x)
        at_method(b, tol) = ConjugateGradient(b, max_iter = n)
        linearsolver = at_method(at_b, tol)
        iters = linearsolve!(
            mulbyat_A!,
            nothing,
            linearsolver,
            at_x,
            at_b;
            max_iters = n,
        )
        exact = at_A \ at_b

        @testset "CuArray test" begin
            @test norm(at_x - exact) / norm(exact) < err_thresh
            @test norm(at_A * at_x - at_b) / norm(at_b) < err_thresh
        end
    end

    mpi_b = MPIStateArray{T}(mpicomm, ArrayType, 4, 4, 4)
    mpi_x = MPIStateArray{T}(mpicomm, ArrayType, 4, 4, 4)
    mpi_A = ArrayType(randn(4^3, 4^3))
    mpi_A .= mpi_A' * mpi_A

    function mpi_mulby!(x, y)
        fy = y.data[:]
        fx = mpi_A * fy
        x.data[:] .= fx[:]
        return nothing
    end

    mpi_b.data[:] .= ArrayType(randn(4^3))
    mpi_x.data[:] .= ArrayType(randn(4^3))

    mpi_method(mpi_b, tol) = ConjugateGradient(mpi_b, max_iter = n)
    linearsolver = mpi_method(mpi_b, tol)
    iters = linearsolve!(
        mpi_mulby!,
        nothing,
        linearsolver,
        mpi_x,
        mpi_b;
        max_iters = n,
    )

    exact = mpi_A \ mpi_b[:]

    @testset "MPIStateArray test" begin
        @test norm(mpi_x.data[:] - exact) / norm(exact) < err_thresh
        mpi_Ax = MPIStateArray{T}(mpicomm, ArrayType, 4, 4, 4)
        mpi_mulby!(mpi_Ax, mpi_x)
        @test norm(mpi_Ax - mpi_b) / norm(mpi_b) < err_thresh
    end

    # ## More Complex Example
    function closure_linear_operator!(A, tup)
        function linear_operator!(y, x)
            alias_x = reshape(x, tup)
            alias_y = reshape(y, tup)
            for i6 in 1:tup[6]
                for i4 in 1:tup[4]
                    for i2 in 1:tup[2]
                        for i1 in 1:tup[1]
                            tmp = alias_x[i1, i2, :, i4, :, i6][:]
                            tmp2 = A[i1, i2, i4, i6] * tmp
                            alias_y[i1, i2, :, i4, :, i6] .=
                                reshape(tmp2, (tup[3], tup[5]))
                        end
                    end
                end
            end
        end
    end

    tup = (3, 4, 7, 2, 20, 2)

    B = [
        randn(tup[3] * tup[5], tup[3] * tup[5])
        for i1 in 1:tup[1], i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
    ]
    columnwise_A = [
        B[i1, i2, i4, i6] * B[i1, i2, i4, i6]' + 10I
        for i1 in 1:tup[1], i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
    ]
    columnwise_inv_A = [
        inv(columnwise_A[i1, i2, i4, i6])
        for i1 in 1:tup[1], i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
    ]
    columnwise_linear_operator! = closure_linear_operator!(columnwise_A, tup)
    columnwise_inverse_linear_operator! =
        closure_linear_operator!(columnwise_inv_A, tup)

    mpi_tup = (tup[1] * tup[2] * tup[3], tup[4], tup[5] * tup[6])
    b = randn(mpi_tup)
    x = randn(mpi_tup)

    linearsolver = ConjugateGradient(
        x,
        max_iter = tup[3] * tup[5],
        dims = (3, 5),
        reshape_tuple = tup,
    )

    iters =
        linearsolve!(columnwise_linear_operator!, nothing, linearsolver, x, b)
    x_exact = copy(x)
    columnwise_inverse_linear_operator!(x_exact, b)

    @testset "Columnwise test" begin
        @test norm(x - x_exact) / norm(x_exact) < err_thresh
    end
end

nothing


text:
# The API server serving functions from srvrfn.jl.
# Runs server in blocking mode when invoked directly with "--runsrvr" argument.
# Call `run_srvr` to start server otherwise.
using JuliaWebAPI
using Logging
using Sockets
using ZMQ
using JSON
using HTTP

include("srvrfn.jl")

const SRVR_ADDR = "tcp://127.0.0.1:9999"
const JSON_RESP_HDRS = Dict{String,String}("Content-Type" => "application/json; charset=utf-8")
const BINARY_RESP_HDRS = Dict{String,String}("Content-Type" => "application/octet-stream")

function run_srvr(fmt, tport, async=false, openaccess=false)
    global_logger(SimpleLogger(open("apisrvr_test.log", "a"), Logging.Info))
    @info("queue is at $SRVR_ADDR")

    api = APIResponder(tport, fmt, nothing, openaccess)
    @info("responding with: $api")

    register(api, testfn1; resp_json=true, resp_headers=JSON_RESP_HDRS)
    register(api, testfn2)
    register(api, testbinary; resp_headers=BINARY_RESP_HDRS)
    register(api, testArray)
    register(api, testFile; resp_json=true, resp_headers=JSON_RESP_HDRS)
    register(api, testException; resp_json=true, resp_headers=JSON_RESP_HDRS)

    process(api; async=async)
end

function test_preproc(req::HTTP.Request)
    respcode = HTTP.header(req, "juliawebapi")
    isempty(respcode) ? JuliaWebAPI.default_preproc(req) : HTTP.Response(parse(Int, respcode))
end

function run_httprpcsrvr(fmt, tport, async=false)
    run_srvr(fmt, tport, true, true)
    apiclnt = APIInvoker(ZMQTransport(SRVR_ADDR, REQ, false), fmt)
    if async
        @async run_http(apiclnt, 8888, test_preproc; reuseaddr=true)
    else
        run_http(apiclnt, 8888, test_preproc; reuseaddr=true)
    end
end

function wait_for_httpsrvr()
    while true
        try
            sock = connect("localhost", 8888)
            close(sock)
            return
        catch
            @info("waiting for httpserver to come up at port 8888...")
            sleep(5)
        end
    end
end

# run in blocking mode if invoked with run flag
!isempty(ARGS) && (ARGS[1] == "--runsrvr") && run_srvr(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPI.ZMQTransport(SRVR_ADDR, ZMQ.REP, true), false)

# run http rpc server if invoked with flag
!isempty(ARGS) && (ARGS[1] == "--runhttprpcsrvr") && run_httprpcsrvr(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPI.ZMQTransport(SRVR_ADDR, ZMQ.REP, true), false)


text:
@testset "Graph Utils" begin

    # define matrix
    M = SparseMatrixCSC([
        0 1 5 6
        1 0 4 10
        5 4 0 3
        6 10 3 0
    ])

    # test children of AdjacencyMatrix
    @test WeightedNetwork.children(M, 1) == [2, 3, 4]
    @test WeightedNetwork.children(M, 2) == [1, 3, 4]
    @test WeightedNetwork.children(M, 3) == [1, 2, 4]
    @test WeightedNetwork.children(M, 4) == [1, 2, 3]

    # test parents of AdjacencyMatrix
    @test WeightedNetwork.parents(M, 1) == [2, 3, 4]
    @test WeightedNetwork.parents(M, 2) == [1, 3, 4]
    @test WeightedNetwork.parents(M, 3) == [1, 2, 4]
    @test WeightedNetwork.parents(M, 4) == [1, 2, 3]

    SyM = zeros(10, 10)
    SyM[3, 1] = -0.9521
    SyM[10, 1] = 0.8540
    SyM[3, 2] = 0.1250
    SyM[8, 2] = 0.0397
    SyM[1, 3] = -0.9521
    SyM[2, 3] = 0.1250
    SyM[8, 3] = 0.5301
    SyM[6, 4] = 1.6742
    SyM[5, 5] = -0.9382
    SyM[8, 5] = 0.1092
    SyM[4, 6] = 1.6742
    SyM[8, 7] = 0.3891
    SyM[10, 7] = -1.6066
    SyM[2, 8] = 0.0397
    SyM[3, 8] = 0.5301
    SyM[5, 8] = 0.1092
    SyM[7, 8] = 0.3891
    SyM[1, 10] = 0.8540
    SyM[7, 10] = -1.6066

    SyM2 = deepcopy(SyM)
    SyM2[7, 3] = 1.9876

    # test symmetric of Matrix
    @test WeightedNetwork.is_symmetric(SparseMatrixCSC(SyM)) == true
    @test WeightedNetwork.is_symmetric(SparseMatrixCSC(SyM2)) == false

end

max_stars_repo_name: tragisch/NetworkOptimization

max_issues_repo_name: tragisch/NetworkOptimization

max_forks_repo_name: tragisch/NetworkOptimization


text:
"""
A Julia Package for finding the trace of large sparse matrices
using Trace Estimation Techniques
"""

module TraceEstimation

include("slq.jl")

end


text:
using TrackFiles
using Test

@testset "TrackFiles.jl" begin
    # Write your tests here.
end


text:
using PyCall
pygui(:tk)
using LinearAlgebra

include("wannierize.jl")
PyPlot.rc("font", family="serif")
PyPlot.rc("xtick", labelsize="x-small")
PyPlot.rc("ytick", labelsize="x-small")
PyPlot.rc("figure", figsize=(4,3))
#PyPlot.rc("text", usetex=false)

# Either run in batch mode, and take arguments from the CLI, or run in interactive mode and specify the file names manually
if(length(ARGS) >= 1)
    filename = ARGS[1]
    methodStr = ARGS[2]
    if methodStr == "parallelTransport"
	method = "parallel transport"
    elseif methodStr == "logInterpolation"
	method = "log interpolation"
    else
	println("ERROR: interpolation method not recognized.")
    end

    if length(ARGS) >= 3
	nbeg = parse(Int64,ARGS[5])
	nend = parse(Int64,ARGS[6])
    else
	nbeg = 1
	nend = Inf
    end
else
    filename = "tests/KaneMele/kanemele_0.0_25"
    #filename = "tests/silicon/silicon"
    method = "parallel transport"
    #method = "log interpolation"
end

read_amn = false #true
read_eig = false #true

p = read_system(filename,read_amn,read_eig)

if(length(ARGS) <1)
    nbeg = 1
    nend = p.nwannier
end

println("Method of interpolation: $method")

println("$(p.nband) bands, wannierizing bands $nbeg:$nend, $(p.N1) x $(p.N2) x $(p.N3) grid, $(p.nntot) neighbors")

#Build the gauge that makes the Bloch frame continuous on the Brillouin Zone.
#This is equivalent to building a set of algebraic decaying Wannier functions
p, interp = make_wannier(p,method)

#Plot and display the results
close("all")
plot_results(p,interp)
print_error(p)

max_stars_repo_name: antoine-levitt/wannier

max_issues_repo_name: antoine-levitt/wannier

max_forks_repo_name: antoine-levitt/wannier


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia}"
JULIA_CMD="${JULIA_CMD:-$JULIA --color=yes --startup-file=no}"
exec $JULIA_CMD -e 'include(popfirst!(ARGS))' "${BASH_SOURCE[0]}" "$@"
=#

module CLITest

using Test

const jlcli = joinpath(dirname(@__DIR__), "jlcli", "jlcli.py")

function test_cli(args)
    julia = Base.julia_cmd().exec[1]
    project = Base.active_project()
    cmd = `$jlcli --julia=$julia --project=$project $args`

    println()
    printstyled("_" ^ displaysize(stdout)[2]; color=:blue)
    @info "Executing a command:" cmd
    flush(stdout)
    flush(stderr)

    @test success(pipeline(cmd; stdout=stdout, stderr=stderr))
end

@testset begin
    test_cli(`--print-result --eval VERSION`)
    test_cli(`--print-result --callany='Statistics=10745b16-79ce-11e8-11f9-7d13ad32a3b2:std' '{"args":[[1,2,3]]}'`)
    test_cli(`--print-result --adhoccli='Unicode=4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5:normalize' -- -h`)
    test_cli(`--print-result --adhoccli='Unicode=4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5:normalize' -- --casefold JuLiA`)
end

end  # module