text:
function queryeaction(;parID::AbstractString="")

    if parID != ""

        allepar = readdlm(eraparametercopy(),',',comments=true);
        @info "$(Dates.now()) - Printing available actions for the given parID $(parID) ..."

        parIDs = allepar[:,2]; modIDs = allepar[:,1];
        if any(parIDs .== parID); modID = modIDs[parIDs .== parID]
            print(reset(),"$(parID) is a valid parID.  You can")
            if any(["csfc","cpre"] .== modID);
                  print(bold()," only analyse ")
            else; print(bold()," either download or analyse ")
            end
            print(reset(),"data for this variable.  If you are unsure as to how to proceed, please call the function",bold()," queryaction() ",reset(),"for more information.")
        end

    else

        @info "$(Dates.now()) - Printing basic help information for ClimateERA.jl ..."

        ftext = joinpath(@__DIR__,"../extra/actiontext.txt");
        lines = readlines(ftext); count = 0; nl = length(lines);
        for l in lines; count += 1;
           if count == nl; print(reset(),"$l\n\n");
           else;           print(reset(),"$l\n");
           end
        end

    end

end

function queryedataset()

    @info "$(Dates.now()) - The following datasets can be downloaded / manipulated with ClimateERA.jl ..."

    ftext = joinpath(@__DIR__,"../extra/datasettext.txt");
    lines = readlines(ftext); count = 0; nl = length(lines);
    for l in lines; count += 1;
       if any(count .== [1,13]); print(bold(),"$l\n");
       elseif count == nl;       print(reset(),"$l\n\n");
       else;                     print(reset(),"$l\n");
       end
    end

end

function queryemod(;modID::AbstractString="")

    mset = ["dsfc","dpre","msfc","mpre","csfc","cpre"]

    if modID != ""

        if any(mset .== modID)

            @info "$(Dates.now()) - $(uppercase(modID)) is a valid module in ClimateERA.jl.  Printing module information ...\n"

            ftext = joinpath(@__DIR__,"../extra/moduletext/$(modID).txt");
            lines = readlines(ftext); count = 0; nl = length(lines);
            for l in lines; count += 1;
               if count == 1;      print(bold(),"$l\n");
               elseif count == nl; print(reset(),"$l\n\n");
               else;               print(reset(),"$l\n");
               end
            end

        else
            @warn "$(Dates.now()) - $(uppercase(modID)) is not a valid module in ClimateERA.jl.  Please query a valid module."
        end

    else

        @info "$(Dates.now()) - The following are the modules available in ClimateERA.jl ...\n"

        for modID in mset
            ftext = joinpath(@__DIR__,"../extra/moduletext/$(modID).txt");
            lines = readlines(ftext); count = 0; nl = length(lines);
            for l in lines; count += 1;
               if count == 1;      print(bold(),"$l\n");
               elseif count == nl; print(reset(),"$l\n\n");
               else;               print(reset(),"$l\n");
               end
            end
        end

    end

end

function queryepar(;parID::AbstractString="",modID::AbstractString="")

    if parID=="" & modID==""; queryeparlist()
    elseif parID=="";         queryeparlist(modID);
    elseif modID=="";         queryeparinfo(parID);
    else;                     queryeparmod(parID,modID);
    end
end

function queryeparinfo(parID::AbstractString)

    allepar = readdlm(eraparametercopy(),',',comments=true);

    parIDs = allepar[:,2]
    if any(parIDs .== parID); parinfo = allepar[parIDs .== parID,:];

        @info "$(Dates.now()) - $(parID) is defined.  Printing variable information ...\n"

        print("    - ",bold(),"Full Name: ",reset(),"$(parinfo[:,6][1])\n");
        print("    - ",bold(),"Units: ",reset(),"$(parinfo[:,7][1])\n");
        print("    - ",bold(),"Module ID: ",reset(),"$(parinfo[:,1][1])\n");
        print("    - ",bold(),"ERA5 ID: ",reset(),"$(parinfo[:,4][1])\n");
        print("    - ",bold(),"ERA-Interim ID: ",reset(),"$(parinfo[:,5][1])\n\n");

    else

        @warn "$(Dates.now()) - $(parID) is not currently defined.  You may add it using the function addepar()."

    end

end

function queryeparmod(parID::AbstractString, modID::AbstractString)

    allepar = readdlm(eraparametercopy(),',',comments=true);

    parIDs = allepar[:,2]
    if any(parIDs .== parID); pmodID = allepar[parIDs .== parID,1][1];

        if pmodID == modID
              @info "$(Dates.now()) - $(parID) is a variable found in $(uppercase(modID))."
        else; @warn "$(Dates.now()) - $(uppercase(modID)) does not contain $(parID)."
        end

    else

        @warn "$(Dates.now()) - $(parID) is not currently defined.  You may add it using the function addepar()."

    end

end

function queryeparlist(modID::AbstractString)

    allepar = readdlm(eraparametercopy(),',',comments=true);
    mset = ["dsfc","dpre","msfc","mpre","csfc","cpre"]
    if any(mset .== modID)

        @info "$(Dates.now()) - $(uppercase(modID)) is a valid module in ClimateERA.jl.  Retrieving parameter information table ...\n"

    else
        @warn "$(Dates.now()) - $(uppercase(modID)) is not a valid module in ClimateERA.jl.  Please query a valid module."
    end

end

## Query ERA Dataset Timestep

hrstep(emod::Dict) = if emod["datasetID"] == 1; return 1; else; return 6 end
hrindy(emod::Dict) = if emod["datasetID"] == 1; return 24; else; return 4 end


text:
using Flux, Flux.Tracker
using LinearAlgebra
using Flux.Tracker: update!
using Word2Vec

sentences = tokenize_document("./data/shakespeare_all.txt")

vocab_hash, vocab, pq, data = process_input(sentences, 3, 5, :cbow)

data = minibatch(data, 500)

wv = WordVectors(25, vocab, vocab_hash)

hs = HeirarchicalSoftmax(25, pq, vocab_hash, vocab)

loss(z) = -1*log(hs(wv(z...)))

ps = params(hs.tree_vectors, wv.vectors)

using ProgressMeter

@showprogress 1 for i in data
        l = loss(i)
        #println(l)
        grad1 = Tracker.gradient(()->l, ps)
        a = grad1[hs.tree_vectors]
        b = grad1[wv.vectors]
        update!(hs.tree_vectors, -0.1*a)
        update!(wv.vectors, -0.1*b)
end

ext: jl


text:
@testset "Launch Configuration" begin
    kernel() = nothing

    device = AMDGPU.default_device()
    queue = AMDGPU.default_queue(device)

    # Group/grid size selection and aliases
    eval(:(@roc groupsize=2 $kernel()))
    eval(:(@roc groupsize=2 gridsize=4 $kernel()))
    eval(:(@roc gridsize=2 $kernel()))
    #eval(:(@roc threads=2 $kernel(1)))
    #eval(:(@roc blocks=4 $kernel(1)))
    #eval(:(@roc threads=2 gridsize=4 $kernel(1)))

    # Device/queue selection and aliases
    # FIXME: Test that device/queue are used!
    eval(:(@roc device=$device $kernel()))
    eval(:(@roc device=$device queue=$queue $kernel()))
    eval(:(@roc queue=$queue $kernel()))
    eval(:(@roc agent=$device $kernel()))
    eval(:(@roc stream=$queue $kernel()))
    eval(:(@roc agent=$device queue=$queue $kernel()))

    # Group size validity
    @test_throws ArgumentError eval(:(@roc groupsize=0 $kernel()))
    eval(:(@roc groupsize=1024 $kernel()))
    @test_throws ArgumentError eval(:(@roc groupsize=1025 $kernel()))

    # No-launch
    @test eval(:(@roc launch=true $kernel())) isa AMDGPU.RuntimeEvent
    @test eval(:(@roc launch=false $kernel())) isa AMDGPU.HostKernel
    @test_throws Exception eval(:(@roc launch=1 $kernel())) # TODO: ArgumentError
end

@testset "No-argument kernel" begin
    kernel() = nothing

    wait(@roc kernel())
end

@testset "Kernel argument alignment" begin
    function kernel(x, y)
        if Int64(x) != y
            error("Fail!")
        end
        nothing
    end
    x = rand(UInt32)
    y = Int64(x)
    wait(@roc kernel(x, y))
end

@testset "Signal waiting" begin
    kernel() = nothing

    sig = @roc kernel()
    wait(sig)
    wait(sig.event)
    wait(sig.event.signal)
    wait(sig.event.signal.signal[])
end


text:
# Assume homogeneous coordinates
function perturb(noise::GaussianNoise, σ::Real, 𝒟::Tuple{AbstractArray, Vararg{AbstractArray}})
    𝓔 = deepcopy(𝒟)
    S = length(𝓔)
    for s = 1:S
        ℳ = 𝓔[s]
        N = length(ℳ)
        for n = 1:N
            𝐦 = ℳ[n]
            D = length(𝐦)
            ℳ[n] = 𝐦 + SVector(randn((D,1))...) #𝐦 + push(σ*SVector(randn((D-1,1))...),0.0)
        end
    end
    𝓔
end


text:
# load packages/modules
using DataFrames
using AutoMLPipeline

# get data
profbdata = getprofb()
X = profbdata[:,2:end];
Y = profbdata[:,1] |> Vector;
topdf(x)=first(x,5)
topdf(profbdata)

# load filters
#### Decomposition
pca = SKPreprocessor("PCA");
fa = SKPreprocessor("FactorAnalysis");
ica = SKPreprocessor("FastICA");

#### Scaler
rb = SKPreprocessor("RobustScaler");
pt = SKPreprocessor("PowerTransformer");
norm = SKPreprocessor("Normalizer");
mx = SKPreprocessor("MinMaxScaler");

#### categorical preprocessing
ohe = OneHotEncoder()

#### Column selector
catf = CatFeatureSelector();
numf = NumFeatureSelector();

#### Learners
rf = SKLearner("RandomForestClassifier",Dict(:impl_args=>Dict(:n_estimators => 100)));
gb = SKLearner("GradientBoostingClassifier");
lsvc = SKLearner("LinearSVC");
mlp = SKLearner("MLPClassifier");
jrf = RandomForest();
stack = StackEnsemble();
svc = SKLearner("SVC");
ada = SKLearner("AdaBoostClassifier");
vote = VoteEnsemble();
best = BestLearner();

# filter categories and hotbit encode
pohe = @pipeline catf |> ohe;
tr = fit_transform!(pohe,X,Y)

pohe = @pipeline numf |> pca;
tr = fit_transform!(pohe,X,Y)

# filter numeric and apply pca and ica
pdec = @pipeline (numf |> pca) + (numf |> ica);
tr = fit_transform!(pdec,X,Y)

# filter numeric, apply rb/pt transforms and ica/pca extractions
ppt = @pipeline (numf |> rb |> ica) + (numf |> pt |> pca)
tr = fit_transform!(ppt,X,Y)

# rf learn baseline
rfp1 = @pipeline ( (catf |> ohe) + (catf |> ohe) + (numf) ) |> rf;
crossvalidate(rfp1, X,Y)

# rf learn: bigrams, 2-blocks
rfp2 = @pipeline ((catf |> ohe) + numf |> rb |> ica) + (numf |> pt |> pca) |> ada;
crossvalidate(rfp2, X,Y)

# lsvc learn: bigrams, 3-blocks
plsvc = @pipeline ((numf |> rb |> pca)+(numf |> rb |> fa)+(numf |> rb |> ica)+(catf |> ohe )) |> lsvc;
pred = fit_transform!(plsvc,X,Y)
crossvalidate(plsvc,X,Y)

@pipelinex ((numf |> rb |> pca)+(numf |> rb |> fa)+(numf |> rb |> ica)+(catf |> ohe )) |> lsvc


text:
N = 10
folders = ["definedtypes", "factors", "newlines", "quoting",
           "scaling", "separators", "typeinference"]

cd(Pkg.dir("DataFrames", "test", "data")) do
    for folder in folders
        for filename in readdir(folder)
            relpath = joinpath(folder, filename)
            df = benchmark(() -> readtable(relpath),
                           "DataFrame I/O",
                           relpath,
                           N)
            # TODO: Keep permanent record
            printtable(df, header=false)
        end
    end
end


text:
fsort([1,2,3])
fsort(rand(Float64, 100))


text:
function nlp_cvx_105_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - e and log expressions
    # - convex constraint, only in given domains
    # Variants
    #   010 - intersection point
    #   011 - intersection point
    #   012 - one binding constraint
    #   013 - one binding constraint

    model = Model(optimizer)

    @variable(model, x, start = 0.1)
    @variable(model, y)

    @objective(model, Min, -x - y)
    @NLconstraint(model, exp(x - 2.0) - 0.5 <= y)
    @NLconstraint(model, log(x) + 0.5 >= y)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, -4.176004405036646, tol = objective_tol)
    return check_solution(
        [x, y],
        [2.687422019398147, 1.488582385638499],
        tol = primal_tol,
    )
end


text:
using SymPy

x, y, t = @syms x y t
# u1 = exp(-t)*(0.5-x)^2*(2-y)^2
# u2 = exp(-t)*(0.5-x^2)^2*(2-sin(y))^2

u1 = exp(-t)*(2-x)*x*(1-y)
u2 = exp(-t)*(2-x)*x*(1-y)
u11 = diff(u1, x)
u12 = diff(u1, y)
u21 = diff(u2, x)
u22 = diff(u2, y)

ε = [u11; u22; u12+u21]

E = 1.0
ν = 0.35
H = E/(1+ν)/(1-2ν)*[
  1-ν ν 0
  ν 1-ν 0
  0 0 (1-2ν)/2
]

σ = H * ε

σ = [σ[1] σ[3];σ[3] σ[2]]

f1 = u1 - (diff(σ[1,1], x) + diff(σ[1,2], y))
f2 = u2 - (diff(σ[2,1], x) + diff(σ[2,2], y))

s1 = sympy.julia_code(f1)
s1 = replace(replace(s1, ".^"=>"^"), ".*"=>"*")

s2 = sympy.julia_code(f2)
s2 = replace(replace(s2, ".^"=>"^"), ".*"=>"*")

println(s1)

println(s2)

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
__precompile__()

module ApproxFun
    using Base, Compat, Plots, FastGaussQuadrature, FastTransforms, DualNumbers, BandedMatrices
    import FixedSizeArrays, ToeplitzMatrices

import Base.LinAlg: BlasFloat

import Base: values,getindex,setindex!,*,.*,+,.+,-,.-,==,<,<=,>,|,
                >=,./,/,.^,^,\,∪,transpose, size

import BandedMatrices: bzeros, bandinds, bandrange, PrintShow, eachbandedindex, bandshift,
                        unsafe_getindex, unsafe_setindex!, bandwidth, AbstractBandedMatrix,
                        dot, dotu, normalize!, flipsign,
                        colstart, colstop, colrange, rowstart, rowstop, rowrange,
                        bandwidths, gbmv!, gbmm!

import Compat: view

import FixedSizeArrays: Vec

export pad!,pad,sample,chop!,complexroots,roots,svfft, reverseorientation

##Testing
export bisectioninv

include("LinearAlgebra/LinearAlgebra.jl")

include("Fun/Fun.jl")
include("Multivariate/Multivariate.jl")
include("Operators/Operator.jl")

include("Domains/Domains.jl")
include("Spaces/Spaces.jl")

## Further extra features

include("PDE/PDE.jl")
include("Extras/Extras.jl")
include("Plot/Plot.jl")

include("docs.jl")

## precompile
function _precompile_()
    precompile(Chebyshev,tuple())
    precompile(Fun,tuple())
    precompile(Base.exp,(Fun{Chebyshev{Interval{Float64}},Float64},))
    precompile(linsolve,(Vector{Operator{Float64}},Vector{Float64}))
    precompile(\,(Vector{Operator{Float64}},Vector{Float64}))
    precompile(+,(Int,Fun{Chebyshev{Interval{Float64}},Float64}))
    precompile(+,(Fun{Chebyshev{Interval{Float64}},Float64},Fun{Chebyshev{Interval{Float64}},Float64}))
    precompile(+,(Fun{ConstantSpace{AnyDomain},Float64},Fun{Chebyshev{Interval{Float64}},Float64}))
    precompile(union,(ConstantSpace{AnyDomain},Chebyshev{Interval{Float64}}))
    precompile(union,(Chebyshev{Interval{Float64}},ConstantSpace{AnyDomain}))
    precompile(Fun,(Fun{ConstantSpace{AnyDomain},Float64},Chebyshev{Interval{Float64}}))
end

_precompile_()

end #module


text:
module TestMatchingPursuit
using Test
using LinearAlgebra
using CompressedSensing: MP, mp, OMP, omp, SP, sp, rmp, sparse_data,
                        ols, ompr, srr, foba, perturb, gomp
using SparseArrays
# NOTE: these tests may rarely fail, because the randomly generated data is "hard"
# a small failure probability is theoretically expected
@testset "Matching Pursuits" begin
    n, m, k = 32, 48, 3
    A, x, b = sparse_data(n = n, m = m, k = k)
    δ = 1e-2 # slightly noisy
    y = perturb(b, δ/2)

    @testset "Matching Pursuit" begin
        xmp = mp(A, b, 10k) # giving more iterations to optimize
        @test isapprox(A*xmp, b, atol = 3δ)
        @test isapprox(xmp, x, atol = 3δ)
    end

    @testset "Orthogonal Matching Pursuit" begin
        # noiseless
        xomp = omp(A, b, k)
        @test xomp.nzind == x.nzind
        @test xomp.nzval ≈ x.nzval
        # noisy
        xomp = omp(A, y, k)
        @test xomp.nzind == x.nzind
        @test isapprox(xomp.nzval, x.nzval, atol = 2δ)
    end

    @testset "Generalized Orthogonal Matching Pursuit" begin
        n, m, k = 32, 48, 3
        A, x, b = sparse_data(n = n, m = m, k = k)
        y = perturb(b, δ/2)
        l = 2 # identify two atoms at a time
        # noiseless
        xgomp = gomp(A, b, l, k)
        @test xgomp.nzind == x.nzind
        @test xgomp.nzval ≈ x.nzval
        # noisy
        xgomp = gomp(A, y, l, k)
        @test xgomp.nzind == x.nzind
        @test isapprox(xgomp.nzval, x.nzval, atol = 2δ)
    end
end

end # TestMatchingPursuit


text:
module DolangTests

using Dolang
using DataStructures
using StaticArrays
using SparseArrays
using Test
#
# tests = length(ARGS) > 0 ? ARGS : [
#                                    "symbolic",
#                                    "incidence",
#                                    "factory",
#                                    "compiler",
#                                    "compiler_new",
#                                    "util",
#                                    "printing",
#                                    ]
# for t in tests
#     include("$(t).jl")
# end

include("symbolic.jl")
include("incidence.jl")
# include("compiler_new.jl")
include("util.jl")
include("printing.jl")
include("minilanguage.jl")

end

max_issues_repo_name: EconForge/Dolang


text:
@testset "Periods" begin
    s = 1.0seconds
    m = 1.0minutes
    h = 1.0hours
    d = 1.0days
    y = 1.0years
    c = 1.0centuries
    @test s == Period(seconds, 1.0)
    @test m == Period(minutes, 1.0)
    @test h == Period(hours, 1.0)
    @test d == Period(days, 1.0)
    @test y == Period(years, 1.0)
    @test c == Period(centuries, 1.0)

    @test seconds(s) == 1.0seconds
    @test seconds(m) == 60.0seconds
    @test seconds(h) == 3600.0seconds
    @test seconds(d) == 86400.0seconds
    @test seconds(y) == 3.15576e7seconds
    @test seconds(c) == 3.15576e9seconds

    @test minutes(s) == (1.0 / 60.0)minutes
    @test minutes(m) == 1.0minutes
    @test minutes(h) == 60.0minutes
    @test minutes(d) == 1440.0minutes
    @test minutes(y) == 525960.0minutes
    @test minutes(c) == 5.2596e7minutes

    @test hours(s) == (1.0 / 3600.0)hours
    @test hours(m) == (1.0 / 60.0)hours
    @test hours(h) == 1.0hours
    @test hours(d) == 24.0hours
    @test hours(y) == 8766.0hours
    @test hours(c) == 876600.0hours

    @test days(s) == (1.0 / 86400.0)days
    @test days(m) == (1.0 / 1440.0)days
    @test days(h) == (1.0 / 24.0)days
    @test days(d) == 1.0days
    @test days(y) == 365.25days
    @test days(c) == 36525.0days

    @test years(s) == (1.0 / 3.15576e7)years
    @test years(m) == (1.0 / 525960.0)years
    @test years(h) == (1.0 / 8766.0)years
    @test years(d) == (1.0 / 365.25)years
    @test years(y) == 1.0years
    @test years(c) == 100.0years

    @test centuries(s) == (1.0 / 3.15576e9)centuries
    @test centuries(m) == (1.0 / 5.2596e7)centuries
    @test centuries(h) == (1.0 / 876600.0)centuries
    @test centuries(d) == (1.0 / 36525.0)centuries
    @test centuries(y) == (1.0 / 100.0)centuries
    @test centuries(c) == 1.0centuries

    @test zero(Period{Year}) == 0.0years
    @test zero(1years) == 0years
    @test zero(1.0years) == 0.0years

    @test eltype(1.0years) == Float64
    @test eltype(typeof(1.0years)) == Float64

    @test 1.0seconds ≈ 1.0seconds
    @test 1.0 * seconds == seconds * 1.0
    @test 3 * 3.0seconds == 3.0seconds * 3
    @test 4.0seconds / 2 == 2.0seconds
    p = 1.0seconds
    @test -p == -1.0seconds
    @test 3.0seconds - 2.0seconds == 1.0seconds

    a = [1.0, 2.0, 3.0]
    @test a * seconds == [1.0seconds, 2.0seconds, 3.0seconds]
    @test a .* seconds == [1.0seconds, 2.0seconds, 3.0seconds]
    @test seconds * a == [1.0seconds, 2.0seconds, 3.0seconds]
    @test seconds .* a == [1.0seconds, 2.0seconds, 3.0seconds]
    @test Base.broadcastable(seconds) isa typeof(Ref(seconds))

    @test unit(1years) == years
    int_rng = 1seconds:3seconds
    @test step(int_rng) == 1seconds
    @test collect(int_rng) == [1seconds, 2seconds, 3seconds]
    float_rng = 1.0seconds:3.0seconds
    @test step(float_rng) == 1.0seconds
    @test collect(float_rng) == [1.0seconds, 2.0seconds, 3.0seconds]
    @test Period{Second,Float64}(1.0seconds) == 1.0seconds

    @test Periods.name(seconds, 1) == "second"
    @test Periods.name(seconds, 2) == "seconds"
    @test Periods.name(seconds, 1.0) == "seconds"
    @test Periods.name(minutes, 1) == "minute"
    @test Periods.name(minutes, 2) == "minutes"
    @test Periods.name(minutes, 1.0) == "minutes"
    @test Periods.name(hours, 1) == "hour"
    @test Periods.name(hours, 2) == "hours"
    @test Periods.name(hours, 1.0) == "hours"
    @test Periods.name(days, 1) == "day"
    @test Periods.name(days, 2) == "days"
    @test Periods.name(days, 1.0) == "days"
    @test Periods.name(years, 1) == "year"
    @test Periods.name(years, 2) == "years"
    @test Periods.name(years, 1.0) == "years"
    @test Periods.name(centuries, 1) == "century"
    @test Periods.name(centuries, 2) == "centuries"
    @test Periods.name(centuries, 1.0) == "centuries"

    @test string(1seconds) == "1 second"
    @test string(Int32(1) * seconds) == "1 second"
    @test string(1.0seconds) == "1.0 seconds"
    @test string(2seconds) == "2 seconds"
    @test string(1minutes) == "1 minute"
    @test string(Int32(1) * minutes) == "1 minute"
    @test string(1.0minutes) == "1.0 minutes"
    @test string(2minutes) == "2 minutes"
    @test string(1hours) == "1 hour"
    @test string(Int32(1) * hours) == "1 hour"
    @test string(1.0hours) == "1.0 hours"
    @test string(2hours) == "2 hours"
    @test string(1days) == "1 day"
    @test string(Int32(1) * days) == "1 day"
    @test string(1.0days) == "1.0 days"
    @test string(2days) == "2 days"
    @test string(1years) == "1 year"
    @test string(Int32(1) * years) == "1 year"
    @test string(1.0years) == "1.0 years"
    @test string(2years) == "2 years"
    @test string(1centuries) == "1 century"
    @test string(Int32(1) * centuries) == "1 century"
    @test string(1.0centuries) == "1.0 centuries"
    @test string(2centuries) == "2 centuries"
end


text:
####################################################################
# In this example we demonstrate using a generic solver
# (our version of BFGS) by sending in the alphafunc and alphagrad!
# functions as the objective and gradient.
#
# Both l2 and huber penalties are demonstrated

using RobustDMD

T = Float32

#--------------------------------------------------------------------
# Generate DMD Synthetic Data
#--------------------------------------------------------------------
println("generate synthetic data...")
# dimensions
m = Int32(100);    # temporal dimension
n = Int16(100);    # spatial dimension
h = 100;
k = 3;      # number of modes
# generate data
sigma = T(1e-4); # size of background noise
mu = T(0.1); # size of spikes
p = T(0.1); # frequency of spikes
X, Xclean, t, alphat, betat = genDMD(m, n, k, sigma, mu; seed=123456, p=p);

#--------------------------------------------------------------------
# Initial Guess
#--------------------------------------------------------------------
println("initialize alpha and B...")
# create initial guess using trap. rule with exact dmd
alpha0, B0 = dmdexactestimate(m,n,k,X,t);

#--------------------------------------------------------------------
# l2 Trial
#--------------------------------------------------------------------
println("l2 loss experiment...")
# DMD variables and parameters
# loss functions
lossf = (z) -> l2_func(z);
lossg = (z) -> l2_grad!(z);
l2Params = DMDParams(k, X, t, lossf, lossg);
l2DMD = DMDVariables(alpha0, B0, l2Params);
# solver variables
l2Svars = DMDVPSolverVariablesLSQ(l2Params);
# define outer problem
func = (alphar) -> alphafunc(alphar, l2DMD, l2Params, l2Svars);
grad! = (galphar, alphar) -> alphagrad!(alphar, galphar, l2DMD, l2Params, l2Svars);
# apply solver
alpha_BFGS_opts = BFGS_options(100, T(1e-6), true, true, true, 1);
alpha_BFGS_vars = BFGS_vars(2*k,T)
@time my_res = My_BFGS(func,grad!,l2DMD.alphar,alpha_BFGS_opts,alpha_BFGS_vars);

#--------------------------------------------------------------------
# Huber Trial
#--------------------------------------------------------------------
println("Huber loss experiment...")
# DMD variables
# loss functions
kappa = T(10*sigma);
lossf = (z) -> huber_func(z,kappa);
lossg = (z) -> huber_grad!(z,kappa);
huberParams = DMDParams(k, X, t, lossf, lossg);
huberDMD = DMDVariables(alpha0, B0, huberParams);
# solver variables
huberSvars = DMDVPSolverVariablesBFGS(huberParams,tol=T(1e-5));
# define outer problem
func = (alphar) -> alphafunc(alphar, huberDMD, huberParams, huberSvars);
grad! = (galphar, alphar) -> alphagrad!(alphar, galphar, huberDMD, huberParams, huberSvars);
# apply solver
alpha_BFGS_opts = BFGS_options(100, T(1e-6), true, true, true, 1);
alpha_BFGS_vars = BFGS_vars(2*k,T)
@time my_res = My_BFGS(func,grad!,huberDMD.alphar,alpha_BFGS_opts,alpha_BFGS_vars);

#--------------------------------------------------------------------
# Compare Result
#--------------------------------------------------------------------
println("compare result...")
err0 = besterrperm(alpha0, alphat);
err1 = besterrperm(l2DMD.alpha, alphat);
err2 = besterrperm(huberDMD.alpha, alphat);

println("err exact DMD ..... ", err0)
println("err l2 fit ........ ", err1)
println("err huber fit ..... ", err2)


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^3*(a+b*exp(c+d*x))^-1, x) == :(1//4*a^-1*x^4+-6*a^-1*d^-4*Polylogarithms.polylog(4,-1*b*a^-1*exp(c+d*x))+-1*a^-1*d^-1*x^3*log(1+b*a^-1*exp(c+d*x))+-3*a^-1*d^-2*x^2*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+6*x*a^-1*d^-3*Polylogarithms.polylog(3,-1*b*a^-1*exp(c+d*x)))
@test integrate(x^2*(a+b*exp(c+d*x))^-1, x) == :(1//3*a^-1*x^3+2*a^-1*d^-3*Polylogarithms.polylog(3,-1*b*a^-1*exp(c+d*x))+-1*a^-1*d^-1*x^2*log(1+b*a^-1*exp(c+d*x))+-2*x*a^-1*d^-2*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x)))
@test integrate(x*(a+b*exp(c+d*x))^-1, x) == :((1/2)*a^-1*x^2+-1*a^-1*d^-2*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+-1*x*a^-1*d^-1*log(1+b*a^-1*exp(c+d*x)))
@test integrate((a+b*exp(c+d*x))^-1, x) == :(x*a^-1+-1*a^-1*d^-1*log(a+b*exp(c+d*x)))
@test integrate((a+b*exp(c+-1*d*x))^-1, x) == :(x*a^-1+a^-1*d^-1*log(a+b*exp(c+-1*d*x)))
@test integrate((a+b*exp(-1c+-1*d*x))^-1, x) == :(x*a^-1+a^-1*d^-1*log(a+b*exp(-1c+-1*d*x)))
@test integrate(x^3*(a+b*exp(c+d*x))^-2, x) == :(1//4*a^-2*x^4+-1*a^-2*d^-1*x^3+-6*a^-2*d^-4*Polylogarithms.polylog(3,-1*b*a^-1*exp(c+d*x))+-6*a^-2*d^-4*Polylogarithms.polylog(4,-1*b*a^-1*exp(c+d*x))+a^-1*d^-1*x^3*(a+b*exp(c+d*x))^-1+-1*a^-2*d^-1*x^3*log(1+b*a^-1*exp(c+d*x))+-3*a^-2*d^-2*x^2*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+3*a^-2*d^-2*x^2*log(1+b*a^-1*exp(c+d*x))+6*x*a^-2*d^-3*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+6*x*a^-2*d^-3*Polylogarithms.polylog(3,-1*b*a^-1*exp(c+d*x)))
@test integrate(x^2*(a+b*exp(c+d*x))^-2, x) == :(1//3*a^-2*x^3+-1*a^-2*d^-1*x^2+2*a^-2*d^-3*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+2*a^-2*d^-3*Polylogarithms.polylog(3,-1*b*a^-1*exp(c+d*x))+a^-1*d^-1*x^2*(a+b*exp(c+d*x))^-1+-1*a^-2*d^-1*x^2*log(1+b*a^-1*exp(c+d*x))+-2*x*a^-2*d^-2*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+2*x*a^-2*d^-2*log(1+b*a^-1*exp(c+d*x)))
@test integrate(x*(a+b*exp(c+d*x))^-2, x) == :((1/2)*a^-2*x^2+a^-2*d^-2*log(a+b*exp(c+d*x))+-1*x*a^-2*d^-1+-1*a^-2*d^-2*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+x*a^-1*d^-1*(a+b*exp(c+d*x))^-1+-1*x*a^-2*d^-1*log(1+b*a^-1*exp(c+d*x)))
@test integrate((a+b*exp(c+d*x))^-2, x) == :(x*a^-2+a^-1*d^-1*(a+b*exp(c+d*x))^-1+-1*a^-2*d^-1*log(a+b*exp(c+d*x)))
@test integrate((a+b*exp(c+-1*d*x))^-2, x) == :(x*a^-2+a^-2*d^-1*log(a+b*exp(c+-1*d*x))+-1*a^-1*d^-1*(a+b*exp(c+-1*d*x))^-1)
@test integrate((a+b*exp(-1c+-1*d*x))^-2, x) == :(x*a^-2+a^-2*d^-1*log(a+b*exp(-1c+-1*d*x))+-1*a^-1*d^-1*(a+b*exp(-1c+-1*d*x))^-1)
@test integrate(x^3*(a+b*exp(c+d*x))^-3, x) == :(1//4*a^-3*x^4+-9*a^-3*d^-4*Polylogarithms.polylog(3,-1*b*a^-1*exp(c+d*x))+-6*a^-3*d^-4*Polylogarithms.polylog(4,-1*b*a^-1*exp(c+d*x))+-3*a^-3*d^-4*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+-3//2*a^-3*d^-1*x^3+3//2*a^-3*d^-2*x^2+a^-2*d^-1*x^3*(a+b*exp(c+d*x))^-1+(1/2)*a^-1*d^-1*x^3*(a+b*exp(c+d*x))^-2+-1*a^-3*d^-1*x^3*log(1+b*a^-1*exp(c+d*x))+-3*x*a^-3*d^-3*log(1+b*a^-1*exp(c+d*x))+-3*a^-3*d^-2*x^2*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+6*x*a^-3*d^-3*Polylogarithms.polylog(3,-1*b*a^-1*exp(c+d*x))+9*x*a^-3*d^-3*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+-3//2*a^-2*d^-2*x^2*(a+b*exp(c+d*x))^-1+9//2*a^-3*d^-2*x^2*log(1+b*a^-1*exp(c+d*x)))
@test integrate(x^2*(a+b*exp(c+d*x))^-3, x) == :(1//3*a^-3*x^3+x*a^-3*d^-2+-1*a^-3*d^-3*log(a+b*exp(c+d*x))+2*a^-3*d^-3*Polylogarithms.polylog(3,-1*b*a^-1*exp(c+d*x))+3*a^-3*d^-3*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+-3//2*a^-3*d^-1*x^2+a^-2*d^-1*x^2*(a+b*exp(c+d*x))^-1+(1/2)*a^-1*d^-1*x^2*(a+b*exp(c+d*x))^-2+-1*x*a^-2*d^-2*(a+b*exp(c+d*x))^-1+-1*a^-3*d^-1*x^2*log(1+b*a^-1*exp(c+d*x))+-2*x*a^-3*d^-2*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+3*x*a^-3*d^-2*log(1+b*a^-1*exp(c+d*x)))
@test integrate(x*(a+b*exp(c+d*x))^-3, x) == :((1/2)*a^-3*x^2+-1*a^-3*d^-2*Polylogarithms.polylog(2,-1*b*a^-1*exp(c+d*x))+-3//2*x*a^-3*d^-1+-1//2*a^-2*d^-2*(a+b*exp(c+d*x))^-1+3//2*a^-3*d^-2*log(a+b*exp(c+d*x))+x*a^-2*d^-1*(a+b*exp(c+d*x))^-1+(1/2)*x*a^-1*d^-1*(a+b*exp(c+d*x))^-2+-1*x*a^-3*d^-1*log(1+b*a^-1*exp(c+d*x)))
@test integrate((a+b*exp(c+d*x))^-3, x) == :(x*a^-3+a^-2*d^-1*(a+b*exp(c+d*x))^-1+(1/2)*a^-1*d^-1*(a+b*exp(c+d*x))^-2+-1*a^-3*d^-1*log(a+b*exp(c+d*x)))
@test integrate((a+b*exp(c+-1*d*x))^-3, x) == :(x*a^-3+a^-3*d^-1*log(a+b*exp(c+-1*d*x))+-1*a^-2*d^-1*(a+b*exp(c+-1*d*x))^-1+-1//2*a^-1*d^-1*(a+b*exp(c+-1*d*x))^-2)
@test integrate((a+b*exp(-1c+-1*d*x))^-3, x) == :(x*a^-3+a^-3*d^-1*log(a+b*exp(-1c+-1*d*x))+-1*a^-2*d^-1*(a+b*exp(-1c+-1*d*x))^-1+-1//2*a^-1*d^-1*(a+b*exp(-1c+-1*d*x))^-2)
@test integrate((c+d*x)^3*(a+b*(F^(g*(e+f*x)))^n), x) == :(1//4*a*d^-1*(c+d*x)^4+b*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*(c+d*x)^3*log(F)^-1+-6*b*d^3*f^-4*g^-4*n^-4*(F^(e*g+f*g*x))^n*log(F)^-4+-3*b*d*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^n*(c+d*x)^2*log(F)^-2+6*b*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^n*log(F)^-3*(c+d*x))
@test integrate((c+d*x)^2*(a+b*(F^(g*(e+f*x)))^n), x) == :(1//3*a*d^-1*(c+d*x)^3+b*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*(c+d*x)^2*log(F)^-1+2*b*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^n*log(F)^-3+-2*b*d*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^n*log(F)^-2*(c+d*x))
@test integrate((a+b*(F^(g*(e+f*x)))^n)*(c+d*x), x) == :((1/2)*a*d^-1*(c+d*x)^2+b*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*log(F)^-1*(c+d*x)+-1*b*d*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^n*log(F)^-2)
@test integrate(a+b*(F^(g*(e+f*x)))^n, x) == :(a*x+b*f^-1*g^-1*n^-1*(F^(g*(e+f*x)))^n*log(F)^-1)
@test integrate((c+d*x)^-1*(a+b*(F^(g*(e+f*x)))^n), x) == :(a*d^-1*log(c+d*x)+b*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*d^-1*(F^(e*g+f*g*x))^n*expintegralei(f*g*n*d^-1*(c+d*x)*log(F)))
@test integrate((c+d*x)^-2*(a+b*(F^(g*(e+f*x)))^n), x) == :(-1*a*d^-1*(c+d*x)^-1+-1*b*d^-1*(F^(e*g+f*g*x))^n*(c+d*x)^-1+b*f*g*n*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*d^-2*(F^(e*g+f*g*x))^n*expintegralei(f*g*n*d^-1*(c+d*x)*log(F))*log(F))
@test integrate((c+d*x)^-3*(a+b*(F^(g*(e+f*x)))^n), x) == :(-1//2*a*d^-1*(c+d*x)^-2+-1//2*b*d^-1*(F^(e*g+f*g*x))^n*(c+d*x)^-2+-1//2*b*f*g*n*d^-2*(F^(e*g+f*g*x))^n*(c+d*x)^-1*log(F)+(1/2)*b*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*d^-3*f^2*g^2*n^2*(F^(e*g+f*g*x))^n*log(F)^2*expintegralei(f*g*n*d^-1*(c+d*x)*log(F)))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^2*(c+d*x)^3, x) == :(1//4*a^2*d^-1*(c+d*x)^4+(1/2)*b^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(2n)*(c+d*x)^3*log(F)^-1+-3//8*b^2*d^3*f^-4*g^-4*n^-4*(F^(e*g+f*g*x))^(2n)*log(F)^-4+-12*a*b*d^3*f^-4*g^-4*n^-4*(F^(e*g+f*g*x))^n*log(F)^-4+2*a*b*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*(c+d*x)^3*log(F)^-1+-3//4*d*b^2*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^(2n)*(c+d*x)^2*log(F)^-2+3//4*b^2*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^(2n)*log(F)^-3*(c+d*x)+-6*a*b*d*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^n*(c+d*x)^2*log(F)^-2+12*a*b*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^n*log(F)^-3*(c+d*x))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^2*(c+d*x)^2, x) == :(1//3*a^2*d^-1*(c+d*x)^3+(1/2)*b^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(2n)*(c+d*x)^2*log(F)^-1+1//4*b^2*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^(2n)*log(F)^-3+2*a*b*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*(c+d*x)^2*log(F)^-1+4*a*b*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^n*log(F)^-3+-1//2*d*b^2*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^(2n)*log(F)^-2*(c+d*x)+-4*a*b*d*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^n*log(F)^-2*(c+d*x))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^2*(c+d*x), x) == :((1/2)*a^2*d^-1*(c+d*x)^2+(1/2)*b^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(2n)*log(F)^-1*(c+d*x)+-1//4*d*b^2*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^(2n)*log(F)^-2+-2*a*b*d*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^n*log(F)^-2+2*a*b*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*log(F)^-1*(c+d*x))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^2, x) == :(x*a^2+(1/2)*b^2*f^-1*g^-1*n^-1*(F^(g*(e+f*x)))^(2n)*log(F)^-1+2*a*b*f^-1*g^-1*n^-1*(F^(g*(e+f*x)))^n*log(F)^-1)
@test integrate((a+b*(F^(g*(e+f*x)))^n)^2*(c+d*x)^-1, x) == :(a^2*d^-1*log(c+d*x)+F^(g*n*(2*e+-2*c*f*d^-1)+-2*g*n*(e+f*x))*b^2*d^-1*(F^(e*g+f*g*x))^(2n)*expintegralei(2*f*g*n*d^-1*(c+d*x)*log(F))+2*a*b*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*d^-1*(F^(e*g+f*g*x))^n*expintegralei(f*g*n*d^-1*(c+d*x)*log(F)))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^2*(c+d*x)^-2, x) == :(-1*a^2*d^-1*(c+d*x)^-1+-1*b^2*d^-1*(F^(e*g+f*g*x))^(2n)*(c+d*x)^-1+-2*a*b*d^-1*(F^(e*g+f*g*x))^n*(c+d*x)^-1+2*f*g*n*F^(g*n*(2*e+-2*c*f*d^-1)+-2*g*n*(e+f*x))*b^2*d^-2*(F^(e*g+f*g*x))^(2n)*expintegralei(2*f*g*n*d^-1*(c+d*x)*log(F))*log(F)+2*a*b*f*g*n*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*d^-2*(F^(e*g+f*g*x))^n*expintegralei(f*g*n*d^-1*(c+d*x)*log(F))*log(F))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^2*(c+d*x)^-3, x) == :(-1//2*a^2*d^-1*(c+d*x)^-2+-1//2*b^2*d^-1*(F^(e*g+f*g*x))^(2n)*(c+d*x)^-2+-1*a*b*d^-1*(F^(e*g+f*g*x))^n*(c+d*x)^-2+-1*f*g*n*b^2*d^-2*(F^(e*g+f*g*x))^(2n)*(c+d*x)^-1*log(F)+-1*a*b*f*g*n*d^-2*(F^(e*g+f*g*x))^n*(c+d*x)^-1*log(F)+2*F^(g*n*(2*e+-2*c*f*d^-1)+-2*g*n*(e+f*x))*b^2*d^-3*f^2*g^2*n^2*(F^(e*g+f*g*x))^(2n)*log(F)^2*expintegralei(2*f*g*n*d^-1*(c+d*x)*log(F))+a*b*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*d^-3*f^2*g^2*n^2*(F^(e*g+f*g*x))^n*log(F)^2*expintegralei(f*g*n*d^-1*(c+d*x)*log(F)))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^3*(c+d*x)^3, x) == :(1//4*a^3*d^-1*(c+d*x)^4+-2//27*b^3*d^3*f^-4*g^-4*n^-4*(F^(e*g+f*g*x))^(3n)*log(F)^-4+1//3*b^3*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(3n)*(c+d*x)^3*log(F)^-1+-18*b*a^2*d^3*f^-4*g^-4*n^-4*(F^(e*g+f*g*x))^n*log(F)^-4+3*b*a^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*(c+d*x)^3*log(F)^-1+-9//8*a*b^2*d^3*f^-4*g^-4*n^-4*(F^(e*g+f*g*x))^(2n)*log(F)^-4+-1//3*d*b^3*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^(3n)*(c+d*x)^2*log(F)^-2+2//9*b^3*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^(3n)*log(F)^-3*(c+d*x)+3//2*a*b^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(2n)*(c+d*x)^3*log(F)^-1+-9*b*d*a^2*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^n*(c+d*x)^2*log(F)^-2+18*b*a^2*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^n*log(F)^-3*(c+d*x)+-9//4*a*d*b^2*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^(2n)*(c+d*x)^2*log(F)^-2+9//4*a*b^2*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^(2n)*log(F)^-3*(c+d*x))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^3*(c+d*x)^2, x) == :(1//3*a^3*d^-1*(c+d*x)^3+1//3*b^3*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(3n)*(c+d*x)^2*log(F)^-1+2//27*b^3*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^(3n)*log(F)^-3+3*b*a^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*(c+d*x)^2*log(F)^-1+6*b*a^2*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^n*log(F)^-3+-2//9*d*b^3*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^(3n)*log(F)^-2*(c+d*x)+3//2*a*b^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(2n)*(c+d*x)^2*log(F)^-1+3//4*a*b^2*d^2*f^-3*g^-3*n^-3*(F^(e*g+f*g*x))^(2n)*log(F)^-3+-6*b*d*a^2*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^n*log(F)^-2*(c+d*x)+-3//2*a*d*b^2*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^(2n)*log(F)^-2*(c+d*x))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^3*(c+d*x), x) == :((1/2)*a^3*d^-1*(c+d*x)^2+-1//9*d*b^3*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^(3n)*log(F)^-2+1//3*b^3*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(3n)*log(F)^-1*(c+d*x)+-3*b*d*a^2*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^n*log(F)^-2+3*b*a^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*log(F)^-1*(c+d*x)+-3//4*a*d*b^2*f^-2*g^-2*n^-2*(F^(e*g+f*g*x))^(2n)*log(F)^-2+3//2*a*b^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(2n)*log(F)^-1*(c+d*x))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^3, x) == :(x*a^3+1//3*b^3*f^-1*g^-1*n^-1*(F^(g*(e+f*x)))^(3n)*log(F)^-1+3*b*a^2*f^-1*g^-1*n^-1*(F^(g*(e+f*x)))^n*log(F)^-1+3//2*a*b^2*f^-1*g^-1*n^-1*(F^(g*(e+f*x)))^(2n)*log(F)^-1)
@test integrate((a+b*(F^(g*(e+f*x)))^n)^3*(c+d*x)^-1, x) == :(a^3*d^-1*log(c+d*x)+F^(g*n*(3*e+-3*c*f*d^-1)+-3*g*n*(e+f*x))*b^3*d^-1*(F^(e*g+f*g*x))^(3n)*expintegralei(3*f*g*n*d^-1*(c+d*x)*log(F))+3*a*F^(g*n*(2*e+-2*c*f*d^-1)+-2*g*n*(e+f*x))*b^2*d^-1*(F^(e*g+f*g*x))^(2n)*expintegralei(2*f*g*n*d^-1*(c+d*x)*log(F))+3*b*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*a^2*d^-1*(F^(e*g+f*g*x))^n*expintegralei(f*g*n*d^-1*(c+d*x)*log(F)))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^3*(c+d*x)^-2, x) == :(-1*a^3*d^-1*(c+d*x)^-1+-1*b^3*d^-1*(F^(e*g+f*g*x))^(3n)*(c+d*x)^-1+-3*a*b^2*d^-1*(F^(e*g+f*g*x))^(2n)*(c+d*x)^-1+-3*b*a^2*d^-1*(F^(e*g+f*g*x))^n*(c+d*x)^-1+3*f*g*n*F^(g*n*(3*e+-3*c*f*d^-1)+-3*g*n*(e+f*x))*b^3*d^-2*(F^(e*g+f*g*x))^(3n)*expintegralei(3*f*g*n*d^-1*(c+d*x)*log(F))*log(F)+3*b*f*g*n*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*a^2*d^-2*(F^(e*g+f*g*x))^n*expintegralei(f*g*n*d^-1*(c+d*x)*log(F))*log(F)+6*a*f*g*n*F^(g*n*(2*e+-2*c*f*d^-1)+-2*g*n*(e+f*x))*b^2*d^-2*(F^(e*g+f*g*x))^(2n)*expintegralei(2*f*g*n*d^-1*(c+d*x)*log(F))*log(F))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^3*(c+d*x)^-3, x) == :(-1//2*a^3*d^-1*(c+d*x)^-2+-1//2*b^3*d^-1*(F^(e*g+f*g*x))^(3n)*(c+d*x)^-2+-3//2*a*b^2*d^-1*(F^(e*g+f*g*x))^(2n)*(c+d*x)^-2+-3//2*b*a^2*d^-1*(F^(e*g+f*g*x))^n*(c+d*x)^-2+-3//2*f*g*n*b^3*d^-2*(F^(e*g+f*g*x))^(3n)*(c+d*x)^-1*log(F)+-3*a*f*g*n*b^2*d^-2*(F^(e*g+f*g*x))^(2n)*(c+d*x)^-1*log(F)+-3//2*b*f*g*n*a^2*d^-2*(F^(e*g+f*g*x))^n*(c+d*x)^-1*log(F)+9//2*F^(g*n*(3*e+-3*c*f*d^-1)+-3*g*n*(e+f*x))*b^3*d^-3*f^2*g^2*n^2*(F^(e*g+f*g*x))^(3n)*log(F)^2*expintegralei(3*f*g*n*d^-1*(c+d*x)*log(F))+6*a*F^(g*n*(2*e+-2*c*f*d^-1)+-2*g*n*(e+f*x))*b^2*d^-3*f^2*g^2*n^2*(F^(e*g+f*g*x))^(2n)*log(F)^2*expintegralei(2*f*g*n*d^-1*(c+d*x)*log(F))+3//2*b*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*a^2*d^-3*f^2*g^2*n^2*(F^(e*g+f*g*x))^n*log(F)^2*expintegralei(f*g*n*d^-1*(c+d*x)*log(F)))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-1*(c+d*x)^3, x) == :(1//4*a^-1*d^-1*(c+d*x)^4+-1*a^-1*f^-1*g^-1*n^-1*(c+d*x)^3*log(F)^-1*log(1+b*a^-1*(F^(g*(e+f*x)))^n)+-6*a^-1*d^3*f^-4*g^-4*n^-4*log(F)^-4*Polylogarithms.polylog(4,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-3*d*a^-1*f^-2*g^-2*n^-2*(c+d*x)^2*log(F)^-2*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+6*a^-1*d^2*f^-3*g^-3*n^-3*log(F)^-3*(c+d*x)*Polylogarithms.polylog(3,-1*b*a^-1*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-1*(c+d*x)^2, x) == :(1//3*a^-1*d^-1*(c+d*x)^3+-1*a^-1*f^-1*g^-1*n^-1*(c+d*x)^2*log(F)^-1*log(1+b*a^-1*(F^(g*(e+f*x)))^n)+2*a^-1*d^2*f^-3*g^-3*n^-3*log(F)^-3*Polylogarithms.polylog(3,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-2*d*a^-1*f^-2*g^-2*n^-2*log(F)^-2*(c+d*x)*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-1*(c+d*x), x) == :((1/2)*a^-1*d^-1*(c+d*x)^2+-1*d*a^-1*f^-2*g^-2*n^-2*log(F)^-2*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-1*a^-1*f^-1*g^-1*n^-1*log(F)^-1*(c+d*x)*log(1+b*a^-1*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-1, x) == :(x*a^-1+-1*a^-1*f^-1*g^-1*n^-1*log(F)^-1*log(a+b*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-2*(c+d*x)^3, x) == :(1//4*a^-2*d^-1*(c+d*x)^4+-1*a^-2*f^-1*g^-1*n^-1*(c+d*x)^3*log(F)^-1+a^-1*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-1*(c+d*x)^3*log(F)^-1+-1*a^-2*f^-1*g^-1*n^-1*(c+d*x)^3*log(F)^-1*log(1+b*a^-1*(F^(g*(e+f*x)))^n)+-6*a^-2*d^3*f^-4*g^-4*n^-4*log(F)^-4*Polylogarithms.polylog(3,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-6*a^-2*d^3*f^-4*g^-4*n^-4*log(F)^-4*Polylogarithms.polylog(4,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-3*d*a^-2*f^-2*g^-2*n^-2*(c+d*x)^2*log(F)^-2*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+3*d*a^-2*f^-2*g^-2*n^-2*(c+d*x)^2*log(F)^-2*log(1+b*a^-1*(F^(g*(e+f*x)))^n)+6*a^-2*d^2*f^-3*g^-3*n^-3*log(F)^-3*(c+d*x)*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+6*a^-2*d^2*f^-3*g^-3*n^-3*log(F)^-3*(c+d*x)*Polylogarithms.polylog(3,-1*b*a^-1*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-2*(c+d*x)^2, x) == :(1//3*a^-2*d^-1*(c+d*x)^3+-1*a^-2*f^-1*g^-1*n^-1*(c+d*x)^2*log(F)^-1+a^-1*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-1*(c+d*x)^2*log(F)^-1+-1*a^-2*f^-1*g^-1*n^-1*(c+d*x)^2*log(F)^-1*log(1+b*a^-1*(F^(g*(e+f*x)))^n)+2*a^-2*d^2*f^-3*g^-3*n^-3*log(F)^-3*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+2*a^-2*d^2*f^-3*g^-3*n^-3*log(F)^-3*Polylogarithms.polylog(3,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-2*d*a^-2*f^-2*g^-2*n^-2*log(F)^-2*(c+d*x)*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+2*d*a^-2*f^-2*g^-2*n^-2*log(F)^-2*(c+d*x)*log(1+b*a^-1*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-2*(c+d*x), x) == :((1/2)*a^-2*d^-1*(c+d*x)^2+d*a^-2*f^-2*g^-2*n^-2*log(F)^-2*log(a+b*(F^(g*(e+f*x)))^n)+a^-1*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-1*log(F)^-1*(c+d*x)+-1*d*x*a^-2*f^-1*g^-1*n^-1*log(F)^-1+-1*d*a^-2*f^-2*g^-2*n^-2*log(F)^-2*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-1*a^-2*f^-1*g^-1*n^-1*log(F)^-1*(c+d*x)*log(1+b*a^-1*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-2, x) == :(x*a^-2+a^-1*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-1*log(F)^-1+-1*a^-2*f^-1*g^-1*n^-1*log(F)^-1*log(a+b*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-3*(c+d*x)^3, x) == :(1//4*a^-3*d^-1*(c+d*x)^4+-3//2*a^-3*f^-1*g^-1*n^-1*(c+d*x)^3*log(F)^-1+a^-2*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-1*(c+d*x)^3*log(F)^-1+(1/2)*a^-1*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-2*(c+d*x)^3*log(F)^-1+-1*a^-3*f^-1*g^-1*n^-1*(c+d*x)^3*log(F)^-1*log(1+b*a^-1*(F^(g*(e+f*x)))^n)+-9*a^-3*d^3*f^-4*g^-4*n^-4*log(F)^-4*Polylogarithms.polylog(3,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-6*a^-3*d^3*f^-4*g^-4*n^-4*log(F)^-4*Polylogarithms.polylog(4,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-3*a^-3*d^3*f^-4*g^-4*n^-4*log(F)^-4*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+3//2*d*a^-3*f^-2*g^-2*n^-2*(c+d*x)^2*log(F)^-2+-3*d*a^-3*f^-2*g^-2*n^-2*(c+d*x)^2*log(F)^-2*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-3*a^-3*d^2*f^-3*g^-3*n^-3*log(F)^-3*(c+d*x)*log(1+b*a^-1*(F^(g*(e+f*x)))^n)+6*a^-3*d^2*f^-3*g^-3*n^-3*log(F)^-3*(c+d*x)*Polylogarithms.polylog(3,-1*b*a^-1*(F^(g*(e+f*x)))^n)+9*a^-3*d^2*f^-3*g^-3*n^-3*log(F)^-3*(c+d*x)*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-3//2*d*a^-2*f^-2*g^-2*n^-2*(a+b*(F^(g*(e+f*x)))^n)^-1*(c+d*x)^2*log(F)^-2+9//2*d*a^-3*f^-2*g^-2*n^-2*(c+d*x)^2*log(F)^-2*log(1+b*a^-1*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-3*(c+d*x)^2, x) == :(1//3*a^-3*d^-1*(c+d*x)^3+-3//2*a^-3*f^-1*g^-1*n^-1*(c+d*x)^2*log(F)^-1+x*a^-3*d^2*f^-2*g^-2*n^-2*log(F)^-2+a^-2*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-1*(c+d*x)^2*log(F)^-1+(1/2)*a^-1*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-2*(c+d*x)^2*log(F)^-1+-1*a^-3*d^2*f^-3*g^-3*n^-3*log(F)^-3*log(a+b*(F^(g*(e+f*x)))^n)+-1*a^-3*f^-1*g^-1*n^-1*(c+d*x)^2*log(F)^-1*log(1+b*a^-1*(F^(g*(e+f*x)))^n)+2*a^-3*d^2*f^-3*g^-3*n^-3*log(F)^-3*Polylogarithms.polylog(3,-1*b*a^-1*(F^(g*(e+f*x)))^n)+3*a^-3*d^2*f^-3*g^-3*n^-3*log(F)^-3*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-1*d*a^-2*f^-2*g^-2*n^-2*(a+b*(F^(g*(e+f*x)))^n)^-1*log(F)^-2*(c+d*x)+-2*d*a^-3*f^-2*g^-2*n^-2*log(F)^-2*(c+d*x)*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+3*d*a^-3*f^-2*g^-2*n^-2*log(F)^-2*(c+d*x)*log(1+b*a^-1*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-3*(c+d*x), x) == :((1/2)*a^-3*d^-1*(c+d*x)^2+a^-2*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-1*log(F)^-1*(c+d*x)+(1/2)*a^-1*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-2*log(F)^-1*(c+d*x)+-1*d*a^-3*f^-2*g^-2*n^-2*log(F)^-2*Polylogarithms.polylog(2,-1*b*a^-1*(F^(g*(e+f*x)))^n)+-1*a^-3*f^-1*g^-1*n^-1*log(F)^-1*(c+d*x)*log(1+b*a^-1*(F^(g*(e+f*x)))^n)+-3//2*d*x*a^-3*f^-1*g^-1*n^-1*log(F)^-1+-1//2*d*a^-2*f^-2*g^-2*n^-2*(a+b*(F^(g*(e+f*x)))^n)^-1*log(F)^-2+3//2*d*a^-3*f^-2*g^-2*n^-2*log(F)^-2*log(a+b*(F^(g*(e+f*x)))^n))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^-3, x) == :(x*a^-3+a^-2*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-1*log(F)^-1+(1/2)*a^-1*f^-1*g^-1*n^-1*(a+b*(F^(g*(e+f*x)))^n)^-2*log(F)^-1+-1*a^-3*f^-1*g^-1*n^-1*log(F)^-1*log(a+b*(F^(g*(e+f*x)))^n))
@test integrate((c+d*x)^(1/2)*(a+b*exp(x)), x) == :(b*(c+d*x)^(1/2)*exp(x)+2//3*a*d^-1*(c+d*x)^3//2+-1//2*b*pi^(1/2)*d^(1/2)*erfi(d^-1//2*(c+d*x)^(1/2))*exp(-1*c*d^-1))
@test integrate((a+b*exp(x))^2*(c+d*x)^(1/2), x) == :((1/2)*b^2*(c+d*x)^(1/2)*exp(2x)+2//3*a^2*d^-1*(c+d*x)^3//2+2*a*b*(c+d*x)^(1/2)*exp(x)+-1*a*b*pi^(1/2)*d^(1/2)*erfi(d^-1//2*(c+d*x)^(1/2))*exp(-1*c*d^-1)+-1//8*2^(1/2)*pi^(1/2)*b^2*d^(1/2)*erfi(2^(1/2)*d^-1//2*(c+d*x)^(1/2))*exp(-2*c*d^-1))
@test integrate((a+b*exp(x))^3*(c+d*x)^(1/2), x) == :(1//3*b^3*(c+d*x)^(1/2)*exp(3x)+2//3*a^3*d^-1*(c+d*x)^3//2+3*b*a^2*(c+d*x)^(1/2)*exp(x)+3//2*a*b^2*(c+d*x)^(1/2)*exp(2x)+-3//2*b*pi^(1/2)*a^2*d^(1/2)*erfi(d^-1//2*(c+d*x)^(1/2))*exp(-1*c*d^-1)+-1//18*3^(1/2)*pi^(1/2)*b^3*d^(1/2)*erfi(3^(1/2)*d^-1//2*(c+d*x)^(1/2))*exp(-3*c*d^-1)+-3//8*a*2^(1/2)*pi^(1/2)*b^2*d^(1/2)*erfi(2^(1/2)*d^-1//2*(c+d*x)^(1/2))*exp(-2*c*d^-1))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^3*(c+d*x)^m, x) == :(a^3*d^-1*(1+m)^-1*(c+d*x)^(1+m)+3^(-1+-1m)*F^(g*n*(3*e+-3*c*f*d^-1)+-3*g*n*(e+f*x))*b^3*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(3n)*(-1*f*g*n*d^-1*(c+d*x)*log(F))^(-1m)*(c+d*x)^m*log(F)^-1*SpecialFunctions.gamma(1+m,-3*f*g*n*d^-1*(c+d*x)*log(F))+3*b*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*a^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*(-1*f*g*n*d^-1*(c+d*x)*log(F))^(-1m)*(c+d*x)^m*log(F)^-1*SpecialFunctions.gamma(1+m,-1*f*g*n*d^-1*(c+d*x)*log(F))+3*a*2^(-1+-1m)*F^(g*n*(2*e+-2*c*f*d^-1)+-2*g*n*(e+f*x))*b^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(2n)*(-1*f*g*n*d^-1*(c+d*x)*log(F))^(-1m)*(c+d*x)^m*log(F)^-1*SpecialFunctions.gamma(1+m,-2*f*g*n*d^-1*(c+d*x)*log(F)))
@test integrate((a+b*(F^(g*(e+f*x)))^n)^2*(c+d*x)^m, x) == :(a^2*d^-1*(1+m)^-1*(c+d*x)^(1+m)+2^(-1+-1m)*F^(g*n*(2*e+-2*c*f*d^-1)+-2*g*n*(e+f*x))*b^2*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^(2n)*(-1*f*g*n*d^-1*(c+d*x)*log(F))^(-1m)*(c+d*x)^m*log(F)^-1*SpecialFunctions.gamma(1+m,-2*f*g*n*d^-1*(c+d*x)*log(F))+2*a*b*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*(-1*f*g*n*d^-1*(c+d*x)*log(F))^(-1m)*(c+d*x)^m*log(F)^-1*SpecialFunctions.gamma(1+m,-1*f*g*n*d^-1*(c+d*x)*log(F)))
@test integrate((c+d*x)^m*(a+b*(F^(g*(e+f*x)))^n), x) == :(a*d^-1*(1+m)^-1*(c+d*x)^(1+m)+b*F^(g*n*(e+-1*c*f*d^-1)+-1*g*n*(e+f*x))*f^-1*g^-1*n^-1*(F^(e*g+f*g*x))^n*(-1*f*g*n*d^-1*(c+d*x)*log(F))^(-1m)*(c+d*x)^m*log(F)^-1*SpecialFunctions.gamma(1+m,-1*f*g*n*d^-1*(c+d*x)*log(F)))
@test integrate(F^(c+d*x)*x^3*(a+b*F^(c+d*x))^-1, x) == :(6*b^-1*d^-4*log(F)^-4*Polylogarithms.polylog(4,-1*b*F^(c+d*x)*a^-1)+b^-1*d^-1*x^3*log(F)^-1*log(1+b*F^(c+d*x)*a^-1)+-6*x*b^-1*d^-3*log(F)^-3*Polylogarithms.polylog(3,-1*b*F^(c+d*x)*a^-1)+3*b^-1*d^-2*x^2*log(F)^-2*Polylogarithms.polylog(2,-1*b*F^(c+d*x)*a^-1))
@test integrate(F^(c+d*x)*x^2*(a+b*F^(c+d*x))^-1, x) == :(-2*b^-1*d^-3*log(F)^-3*Polylogarithms.polylog(3,-1*b*F^(c+d*x)*a^-1)+b^-1*d^-1*x^2*log(F)^-1*log(1+b*F^(c+d*x)*a^-1)+2*x*b^-1*d^-2*log(F)^-2*Polylogarithms.polylog(2,-1*b*F^(c+d*x)*a^-1))
@test integrate(x*F^(c+d*x)*(a+b*F^(c+d*x))^-1, x) == :(b^-1*d^-2*log(F)^-2*Polylogarithms.polylog(2,-1*b*F^(c+d*x)*a^-1)+x*b^-1*d^-1*log(F)^-1*log(1+b*F^(c+d*x)*a^-1))
@test integrate(F^(c+d*x)*(a+b*F^(c+d*x))^-1, x) == :(b^-1*d^-1*log(F)^-1*log(a+b*F^(c+d*x)))
@test integrate(F^(c+d*x)*x^3*(a+b*F^(c+d*x))^-2, x) == :(a^-1*b^-1*d^-1*x^3*log(F)^-1+-1*b^-1*d^-1*x^3*(a+b*F^(c+d*x))^-1*log(F)^-1+6*a^-1*b^-1*d^-4*log(F)^-4*Polylogarithms.polylog(3,-1*b*F^(c+d*x)*a^-1)+-6*x*a^-1*b^-1*d^-3*log(F)^-3*Polylogarithms.polylog(2,-1*b*F^(c+d*x)*a^-1)+-3*a^-1*b^-1*d^-2*x^2*log(F)^-2*log(1+b*F^(c+d*x)*a^-1))
@test integrate(F^(c+d*x)*x^2*(a+b*F^(c+d*x))^-2, x) == :(a^-1*b^-1*d^-1*x^2*log(F)^-1+-1*b^-1*d^-1*x^2*(a+b*F^(c+d*x))^-1*log(F)^-1+-2*a^-1*b^-1*d^-3*log(F)^-3*Polylogarithms.polylog(2,-1*b*F^(c+d*x)*a^-1)+-2*x*a^-1*b^-1*d^-2*log(F)^-2*log(1+b*F^(c+d*x)*a^-1))
@test integrate(x*F^(c+d*x)*(a+b*F^(c+d*x))^-2, x) == :(x*a^-1*b^-1*d^-1*log(F)^-1+-1*x*b^-1*d^-1*(a+b*F^(c+d*x))^-1*log(F)^-1+-1*a^-1*b^-1*d^-2*log(F)^-2*log(a+b*F^(c+d*x)))
@test integrate(F^(c+d*x)*(a+b*F^(c+d*x))^-2, x) == :(-1*b^-1*d^-1*(a+b*F^(c+d*x))^-1*log(F)^-1)
@test integrate(F^(c+d*x)*x^3*(a+b*F^(c+d*x))^-3, x) == :((1/2)*a^-2*b^-1*d^-1*x^3*log(F)^-1+3*a^-2*b^-1*d^-4*log(F)^-4*Polylogarithms.polylog(2,-1*b*F^(c+d*x)*a^-1)+3*a^-2*b^-1*d^-4*log(F)^-4*Polylogarithms.polylog(3,-1*b*F^(c+d*x)*a^-1)+-3//2*a^-2*b^-1*d^-2*x^2*log(F)^-2+-1//2*b^-1*d^-1*x^3*(a+b*F^(c+d*x))^-2*log(F)^-1+-3*x*a^-2*b^-1*d^-3*log(F)^-3*Polylogarithms.polylog(2,-1*b*F^(c+d*x)*a^-1)+3*x*a^-2*b^-1*d^-3*log(F)^-3*log(1+b*F^(c+d*x)*a^-1)+-3//2*a^-2*b^-1*d^-2*x^2*log(F)^-2*log(1+b*F^(c+d*x)*a^-1)+3//2*a^-1*b^-1*d^-2*x^2*(a+b*F^(c+d*x))^-1*log(F)^-2)
@test integrate(F^(c+d*x)*x^2*(a+b*F^(c+d*x))^-3, x) == :(a^-2*b^-1*d^-3*log(F)^-3*log(a+b*F^(c+d*x))+(1/2)*a^-2*b^-1*d^-1*x^2*log(F)^-1+-1*x*a^-2*b^-1*d^-2*log(F)^-2+-1*a^-2*b^-1*d^-3*log(F)^-3*Polylogarithms.polylog(2,-1*b*F^(c+d*x)*a^-1)+-1//2*b^-1*d^-1*x^2*(a+b*F^(c+d*x))^-2*log(F)^-1+x*a^-1*b^-1*d^-2*(a+b*F^(c+d*x))^-1*log(F)^-2+-1*x*a^-2*b^-1*d^-2*log(F)^-2*log(1+b*F^(c+d*x)*a^-1))
@test integrate(x*F^(c+d*x)*(a+b*F^(c+d*x))^-3, x) == :((1/2)*x*a^-2*b^-1*d^-1*log(F)^-1+(1/2)*a^-1*b^-1*d^-2*(a+b*F^(c+d*x))^-1*log(F)^-2+-1//2*x*b^-1*d^-1*(a+b*F^(c+d*x))^-2*log(F)^-1+-1//2*a^-2*b^-1*d^-2*log(F)^-2*log(a+b*F^(c+d*x)))
@test integrate(F^(c+d*x)*(a+b*F^(c+d*x))^-3, x) == :(-1//2*b^-1*d^-1*(a+b*F^(c+d*x))^-2*log(F)^-1)

max_stars_repo_path: test/2 Exponentials/2.2 (c+d x)^m (F^(g (e+f x)))^n (a+b (F^(g (e+f x)))^n)^p.jl


text:
"""
	grdvolume(cmd0::String="", arg1=nothing, kwargs...)

Reads one 2-D grid and returns xyz-triplets.

Full option list at [`grdvolume`]($(GMTdoc)grdvolume.html)

Parameters
----------

- **C** | **cont** | **contour** :: [Type => Str | List]   ``Arg = cval or low/high/delta or rlow/high or rcval``

    Find area, volume and mean height (volume/area) inside the cval contour.
    ($(GMTdoc)grdvolume.html#c)
- **L** | **base_level** :: [Type => Number]          ``Arg = base``

    Also add in the volume from the level of the contour down to base [Default base is contour].
    ($(GMTdoc)grdvolume.html#l)
- $(GMT.opt_R)
- **S** | **unit** :: [Type => Str]              ``Arg = e|f|k|M|n|u``

    For geographical grids, append a unit from e|f|k|M|n|u [Default is meter (e)].
    ($(GMTdoc)grdvolume.html#s)
- **T** :: [Type => Str]                        ``Arg = [c|h]``

    Determine the single contour that maximized the average height (= volume/area).
    ($(GMTdoc)grdvolume.html#t)
- $(GMT.opt_V)
- **Z** | **scale** :: [Type => Str or List]     ``Arg = fact[/shift]``

    Optionally subtract shift before scaling data by fact. [Default is no scaling].
    ($(GMTdoc)grdvolume.html#z)
- $(GMT.opt_f)
- $(GMT.opt_o)
"""
function grdvolume(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdvolume", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:R :V_params :f :o])
	cmd  = parse_these_opts(cmd, d, [[:C :cont :contour], [:L :base_level], [:S :unit], [:T], [:Z :scale]])
	common_grd(d, cmd0, cmd, "grdvolume ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdvolume(arg1, cmd0::String=""; kw...) = grdvolume(cmd0, arg1; kw...)


text:
function change_form(::Type{LPForm{T, AT, VT}}, lp::LPForm) where {T, AT, VT}
    return LPForm{T, AT, VT}(
        lp.direction,
        lp.c,
        lp.A,
        lp.c_lb,
        lp.c_ub,
        lp.v_lb,
        lp.v_ub
    )
end

function change_form(::Type{LPForm{T, AT, VT}}, lp::LPStandardForm{T}) where {T, AT, VT}
    return LPForm{T, AT, VT}(
        lp.direction,
        lp.c,
        lp.A,
        lp.b,
        lp.b,
        fill(zero(T), length(lp.c)),
        fill(typemax(T), length(lp.c)),
    )
end
function change_form(::Type{LPForm{T, AT, VT}}, lp::LPGeometricForm{T}) where {T, AT, VT}
    return LPForm{T, AT, VT}(
        lp.direction,
        lp.c,
        lp.A,
        fill(typemin(T), length(lp.b)),
        lp.b,
        fill(typemin(T), length(lp.c)),
        fill(typemax(T), length(lp.c)),
    )
end
function change_form(::Type{LPForm{T, AT, VT}}, lp::LPSolverForm{T}) where {T, AT, VT}
    c_lb = fill(typemin(T), length(lp.b))
    c_ub = fill(typemax(T), length(lp.b))
    for i in eachindex(lp.b)
        if lp.senses[i] == LESS_THAN
            c_ub[i] = lp.b[i]
        elseif lp.senses[i] == GREATER_THAN
            c_lb[i] = lp.b[i]
        elseif lp.senses[i] == EQUAL_TO
            c_lb[i] = lp.b[i]
            c_ub[i] = lp.b[i]
        else
            error("invalid sign $(lp.senses[i])")
        end
    end
    return LPForm{T, AT, VT}(
        lp.direction,
        lp.c,
        lp.A,
        c_lb,
        c_ub,
        lp.v_lb,
        lp.v_ub,
    )
end

function change_form(::Type{LPGeometricForm{T, AT, VT}}, lp::LPGeometricForm) where {T, AT, VT}
    return LPGeometricForm{T, AT, VT}(
        lp.direction,
        lp.c,
        lp.A,
        lp.b
    )
end
function change_form(::Type{LPGeometricForm{T, AT, VT}}, lp::LPForm{T}) where {T, AT, VT}
    has_c_upper = Int[]
    has_c_lower = Int[]
    sizehint!(has_c_upper, length(lp.c_ub))
    sizehint!(has_c_lower, length(lp.c_ub))
    for i in eachindex(lp.c_ub)
        if _no_upper(lp.c_ub[i])
            push!(has_c_upper, i)
        end
        if _no_lower(lp.c_lb[i])
            push!(has_c_lower, i)
        end
    end
    has_v_upper = Int[]
    has_v_lower = Int[]
    sizehint!(has_v_upper, length(lp.v_ub))
    sizehint!(has_v_lower, length(lp.v_ub))
    for i in eachindex(lp.v_ub)
        if _no_upper(lp.v_ub[i])
            push!(has_v_upper, i)
        end
        if _no_lower(lp.v_lb[i])
            push!(has_v_lower, i)
        end
    end
    Id = Matrix{T}(I, length(lp.c), length(lp.c))
    new_A = vcat(
        lp.A[has_c_upper,:],
        -lp.A[has_c_lower,:],
        Id[has_v_upper,:],
        -Id[has_v_lower,:],
                )
    new_b = vcat(
        lp.c_ub[has_c_upper],
        -lp.c_lb[has_c_lower],
        lp.v_ub[has_v_upper],
        -lp.v_lb[has_v_lower],
    )
    return LPGeometricForm{T, AT, VT}(
        lp.direction,
        lp.c,
        new_A,
        new_b
    )
end
function change_form(::Type{LPGeometricForm{T, AT, VT}}, lp::F) where {T, AT, VT, F <: AbstractLPForm{T}}
    temp_lp = change_form(LPForm{T, AT, VT}, lp)
    return change_form(LPGeometricForm{T, AT, VT}, temp_lp)
end

function change_form(::Type{LPStandardForm{T, AT, VT}}, lp::LPStandardForm) where {T, AT, VT}
    return LPStandardForm(
        lp.direction,
        lp.c,
        lp.A,
        lp.b
    )
end
function change_form(::Type{LPStandardForm{T, AT, VT}}, lp::LPGeometricForm{T}) where {T, AT, VT}
    new_A = hcat(
        lp.A,
        -lp.A,
        AT(I, length(lp.b), length(lp.b))
    )
    new_c = vcat(
        lp.c,
        -lp.c,
        fill(0.0, length(lp.b))
    )
    return LPStandardForm{T, AT, VT}(
        lp.direction,
        new_c,
        new_A,
        copy(lp.b)
    )
end
function change_form(::Type{LPStandardForm{T, AT, VT}}, lp::F) where {T, AT, VT, F <: AbstractLPForm{T}}
    temp_lp = change_form(LPForm{T, AT, VT}, lp)
    new_lp = change_form(LPGeometricForm{T, AT, VT}, temp_lp)
    change_form(LPStandardForm{T, AT, VT}, new_lp)
end

function change_form(::Type{LPSolverForm{T, AT, VT}}, lp::LPSolverForm) where {T, AT, VT}
    return LPSolverForm{T, AT, VT}(
        lp.direction,
        lp.c,
        lp.A,
        lp.b,
        lp.senses,
        lp.v_lb,
        lp.v_ub
    )
end
function change_form(::Type{LPSolverForm{T, AT, VT}}, lp::LPForm{T}) where {T, AT, VT}
    new_A = copy(lp.A)
    senses = fill(LESS_THAN, length(lp.c_lb))
    new_b = fill(NaN, length(lp.c_lb))
    for i in eachindex(lp.c_lb)
        if lp.c_lb[i] == lp.c_ub[i]
            senses[i] = EQUAL_TO
            new_b[i] = lp.c_lb[i]
        elseif lp.c_lb[i] > -Inf && lp.c_ub[i] < Inf
            senses[i] = GREATER_THAN
            new_b[i] = lp.c_lb[i]
            push!(new_b, lp.c_ub[i])
            push!(sense, LESS_THAN)
            new_A = vcat(new_A, lp.A[i,:])
        elseif lp.c_lb[i] > -Inf
            senses[i] = GREATER_THAN
            new_b[i] = lp.c_lb[i]
        elseif lp.c_ub[i] < Inf
            senses[i] = LESS_THAN
            new_b[i] = lp.c_ub[i]
        end
    end
    return LPSolverForm{T, AT, VT}(
        lp.direction,
        lp.c,
        new_A,
        new_b,
        senses,
        lp.v_lb,
        lp.v_ub,
    )
end
function change_form(::Type{LPSolverForm{T, AT, VT}}, lp::F) where {T, AT, VT, F <: AbstractLPForm{T}}
    temp_lp = change_form(LPForm{T, AT, VT}, lp)
    change_form(LPSolverForm{T, AT, VT}, temp_lp)
end


text:
module EDKit

using LinearAlgebra, SparseArrays

import Base: +, -, *, /, Array, Vector, Matrix, size, length, eltype, digits
import LinearAlgebra: norm, mul!
import SparseArrays: sparse

export AbstractBasis, content, base, index, change!

"""
    AbstractBasis

Abstract type for all concrete `Basis` types.
A typical concrete `Basis` type usually has the field:

- `dgt`: Digits representation of a state. For a symmetric basis like translational-invariant basis, `dgt` is just a representation of the state.
- `I`  : List of indices in this basis.
- `R`  : List of normalization of each state in the basis.
- `B`  : Base of the basis.

Correspondingly, there are property functions for a `Basis` type:

- `digits` : Return digits representation.
- `content`: Return list of indices in this basis.
- `norm`   : Return list of norm.
- `base`   : Return the base of the basis.

In addition, we have the following property functions:

- `eltype` : Used for type promotion when constructing matrix representation of an `Operator`.
- `length` : Return the length of the digits of the basis.
- `size`   : Return the size of the matrix that can be constructed for an `Operator`.

The most inportant function for a `Basis` type is the index of a given digits and the digits representation of its ith content:

- `index`  : Return the norm and the index of the digits of the given basis.
- `change!`: Change the digits to the target index, and return the new normalization.

If we want to calculate the schmidt decomposition (in real space), we can use the `schmidt` function.
"""
abstract type AbstractBasis end

include("ToolKit.jl")
include("Schmidt.jl")
include("Symmetries.jl")
include("Basis.jl")
include("Operator.jl")

for file in readdir("$(@__DIR__)/algorithms/")
    if file[end-2:end] == ".jl"
        include("$(@__DIR__)/algorithms/$file")
    end
end

# Default function definitions:
@inline digits(b::AbstractBasis) = b.dgt
@inline content(b::AbstractBasis) = b.I
@inline content(b::AbstractBasis, i::Integer) = b.I[i]
@inline norm(b::AbstractBasis) = b.R
@inline norm(b::AbstractBasis, i::Integer) = b.R[i]
@inline base(b::AbstractBasis) = b.B
@inline eltype(::AbstractBasis) = ComplexF64
@inline length(b::AbstractBasis) = length(b.dgt)
@inline size(b::AbstractBasis, i::Integer) = isone(i) || isequal(i, 2) ? length(b.I) : 1
@inline size(b::AbstractBasis) = (size(b, 1), size(b, 2))
@inline change!(b::AbstractBasis, i::Integer) = (change!(b.dgt, content(b, i), base=b.B); norm(b, i))

end # module


text:
module ReadmeTests

println("Running code from README.md")

include_string(open(Pkg.dir("Fixtures","README.md")) do f
  join(filter(readlines(f)) do line
    (length(line)>0 && line[1]=='\t') || (length(line)>3 && line[1:4]=="    ")
  end, '\n')
end)

println("Done")

end


text:
using GraphQLClient

function verify_ipfshashes(xs::AbstractVector{<:AbstractString})
    return all(verify_ipfshash.(xs))
end

function ipfshash_in(
    whitelist::AbstractVector{T}, pinnedlist::AbstractVector{T}
) where {T<:AbstractString}
    return whitelist ∪ pinnedlist
end

function ipfshash_not_in(
    blacklist::AbstractVector{T}, frozenlist::AbstractVector{T}
) where {T<:AbstractString}
    return blacklist ∪ frozenlist
end

function frozen_stake(
    client::Client, id::AbstractString, frozenlist::AbstractVector{<:AbstractString}
)
    allocations_where_query = Dict("status" => "Active", "indexer" => id)

    allocations_query = GQLQuery(
        Dict("where" => allocations_where_query),
        ["allocatedTokens", "subgraphDeployment{ ipfsHash }"],
    )
    allocations_data = query(client, "allocations"; query_args=allocations_query.args, output_fields=allocations_query.fields).data["allocations"]
    return sum(
        map(
            a -> togrt(a["allocatedTokens"]),
            filter(
                a -> (a["subgraphDeployment"]["ipfsHash"] in frozenlist), allocations_data
            ),
        ),
    )
end

function query_subgraphs(
    client::Client, ipfsin::AbstractVector{T}, ipfs_notin::AbstractVector{T}
) where {T<:AbstractString}
    subgraph_where_query::Dict{AbstractString,Union{AbstractString,AbstractVector{<:AbstractString}}} = Dict(
        "signalledTokens_gte" => "1000000000000000000000"
    )
    if !isempty(ipfsin)
        subgraph_where_query["ipfsHash_in"] = ipfsin
    end
    if !isempty(ipfs_notin)
        subgraph_where_query["ipfsHash_not_in"] = ipfs_notin
    end
    subgraph_query = GQLQuery(
        Dict(
            "first" => 1000,
            "orderBy" => "signalledTokens",
            "orderDirection" => "desc",
            "where" => subgraph_where_query,
        ),
        ["id", "ipfsHash", "signalledTokens"],
    )
    subgraphs_data = query(client, "subgraphDeployments"; query_args=subgraph_query.args, output_fields=subgraph_query.fields).data["subgraphDeployments"]
    subgraphs = map(
        x -> SubgraphDeployment(x["id"], x["ipfsHash"], x["signalledTokens"]),
        subgraphs_data,
    )
    return subgraphs
end

function query_indexer_allocations(client::Client, indexer_id::AbstractString)
    indexer_query = GQLQuery(
        Dict(
            "first" => 1000,
            "where" =>
                Dict("stakedTokens_gte" => "100000000000000000000000", "id" => indexer_id),
        ),
        ["allocations{id,subgraphDeployment{ipfsHash}}"],
    )
    indexer_data = query(client, "indexers"; query_args=indexer_query.args, output_fields=indexer_query.fields).data["indexers"]
    indexer = Indexer(indexer_data[1]["allocations"])
    return indexer.allocations
end

function query_indexers(client::Client, subgraphs::AbstractVector{SubgraphDeployment})
    subgraph_ids = id.(subgraphs)
    allocations_where_query = Dict(
        "status" => "Active", "subgraphDeployment_in" => subgraph_ids
    )
    allocations_query = GraphQLClient.directly_write_query_args(
        Dict("where" => allocations_where_query)
    )
    indexer_query = GQLQuery(
        Dict(
            "first" => 1000,
            "where" => Dict("stakedTokens_gte" => "100000000000000000000000"),
        ),
        [
            "id",
            "delegatedTokens",
            "stakedTokens",
            "allocations($allocations_query){allocatedTokens,createdAtEpoch,subgraphDeployment{ipfsHash}}",
        ],
    )
    indexers_data = query(client, "indexers"; query_args=indexer_query.args, output_fields=indexer_query.fields).data["indexers"]
    indexers = map(
        x -> Indexer(x["id"], x["delegatedTokens"], x["stakedTokens"], x["allocations"]),
        indexers_data,
    )
    return indexers
end

function snapshot(
    client::Client, ipfsin::AbstractVector{T}, ipfs_notin::AbstractVector{T}
) where {T<:AbstractString}
    subgraphs = query_subgraphs(client, ipfsin, ipfs_notin)
    indexers = query_indexers(client, subgraphs)

    # Make repository
    repository = Repository(indexers, subgraphs)
    return repository
end


text:
module TestSomething

using Maybe
using Maybe: @something
using Maybe.Implementations: statements_in_braces
using Test

include("utils.jl")

@testset "statements_in_braces" begin
    stmts(ex::Expr) = statements_in_braces(ex.args[end])
    jl15"""
    @test stmts(:(@_{})) == Any[]
    @test stmts(:(@_{a})) == Any[:a]
    @test stmts(:(@_{a; b})) == Any[:a, :b]
    @test stmts(:(@_{a; b; c})) == Any[:a, :b, :c]
    @test stmts(:(@_{a; b; c;})) == Any[:a, :b, :c]
    @test_throws ErrorException stmts(:(@_{a; b, c;}))
    @test_throws ErrorException stmts(:(@_{a; b, c}))
    @test_throws ErrorException stmts(:(@_{a; b, c; d}))
    @test_throws ErrorException stmts(:(@_{a, b; c; d}))
    """
end

@testset "@something" begin
    hist = []
    rec(x) = push!(hist, x)
    @test @something(
        (rec(1); nothing),  # commas required
        (rec(2); nothing),
        (rec(3); :a),
        (rec(4); :b)
    ) === :a
    @test hist == [1, 2, 3]

    jl15"""
    hist = []
    @test @something{
        (rec(1); nothing);
        (rec(2); nothing);
        (rec(3); :a);
        (rec(4); :b);
    } === :a
    @test hist == [1, 2, 3]
    """

    hist = []
    @test (@something {
        (rec(1); nothing);
        (rec(2); nothing);
        (rec(3); :a);
        (rec(4); :b);
    }) === :a
    @test hist == [1, 2, 3]
end

@testset "with @?" begin
    d = Dict(:a => Dict(:b => 1, :c => nothing))
    @test (@something {
        @? d[:A][:B];
        @? d[:A][:b];
        @? d[:a][:b];
        @? d[:c][:d];
    }) === 1
    @test (@something {@? d[:a][:c]}) === nothing
    @test (@something (@? d[:a][:c])) === nothing
    jl15"""
    @test (@something {
        @? d[:A][:B];
        @? d[:A][:b];
        @? d[:a][:b];
        @? d[:c][:d];
    }) === 1
    @test (@something {@? d[:a][:c]}) === nothing
    @test (@something (@? d[:a][:c])) === nothing
    """
end

end  # module


text:
# Contains code for reconnecting a WebSocket.
#
# Retrying is done using "backoff", so that it doesn't retry too often. The backoffs are done in two
# layers. The first calculates a fixed backoff based on the atan trigonometric function. This is
# defined so that it will reach about 90% of the max backoff after 12 attempts.
# The other layer is a randomization of the previous backoff, jiggling the timeout a few seconds
# longer or shorter.

import Base: reset, retry

export AbstractBackoff, Backoff, RandomizedBackoff, reset, backoff_min, backoff_max
export AbstractRetry, Retry, retry, set_function

abstract type AbstractBackoff end

"A backoff that follows a atan curve, and reaches about 90% of max backoff in 12 attempts."
mutable struct Backoff <: AbstractBackoff
    min::Float64
    max::Float64
    state::Int

    Backoff(min::Float64, max::Float64) = new(min, max, 0)
end

"Reset the backoff to its initial state."
reset(b::Backoff) = b.state = 0

backoff_min(b::Backoff) = b.min
backoff_max(b::Backoff) = b.max

# make Backoff callable
"Get the next backoff value."
function (b::Backoff)()
    v = b.min + atan(b.state*b.state/32) * 2 / pi * (b.max - b.min)
    b.state += 1
    v
end

"Randomize another backoff by adding some randomness to the backoff time."
mutable struct RandomizedBackoff <: AbstractBackoff
    backoff::AbstractBackoff
    rng::AbstractRNG
    interval::Float64
end

show(io::IO, r::RandomizedBackoff) =
    show(io, "RandomizedBackoff($(r.backoff), $(r.interval))")

reset(b::RandomizedBackoff) = reset(b.backoff)
backoff_min(b::RandomizedBackoff) = backoff_min(b.backoff)
backoff_max(b::RandomizedBackoff) = backoff_max(b.backoff)

# make RandomizedBackoff callable
function (b::RandomizedBackoff)()
    (r,) = (rand(b.rng, Float64, 1) - 0.5) * 2 * b.interval
    v = b.backoff()
    max(backoff_min(b), min(backoff_max(b), v + r))
end

abstract type AbstractRetry end

default_timer_factory = (f, d) -> Timer(f, d)

"Start a timer for some function, based on a backoff."
mutable struct Retry <: AbstractRetry
    backoff::AbstractBackoff
    fun::Function
    timer_fun::Function

    Retry(backoff::AbstractBackoff, fun::Function;
          timer_fun::Function=default_timer_factory) = new(backoff, (t) -> fun(), timer_fun)
end

function retry(r::Retry)
    backoff_time = r.backoff()
    r.timer_fun(r.fun, backoff_time)
end

reset(r::Retry) = reset(r.backoff)

set_function(r::Retry, f::Function) = r.fun = (t) -> f()


text:
using DataFrames
using Dates
using Test

using Stonks: to_dataframe
using Stonks.Models: AssetInfo, AssetPrice

include("test_utils.jl")

@testset "Conversions" begin

  prices = test_price_data()
  info = test_info_data()

  @testset "Vector{AssetInfo} to DataFrame" begin
    df = to_dataframe(info)
    model_types = [
      (String(name), T) for (name, T) in zip(fieldnames(AssetInfo), AssetInfo.types)
    ]
    df_types = [(name, eltype(x)) for (name, x) in zip(names(df), eachcol(df))]
    @test length(info) == nrow(df)
    @test model_types == df_types
  end

  @testset "Vector{AssetPrice} to DataFrame" begin
    df = to_dataframe(prices)
    model_types = [
      (String(name), T) for (name, T) in zip(fieldnames(AssetPrice), AssetPrice.types)
    ]
    df_types = [(name, eltype(x)) for (name, x) in zip(names(df), eachcol(df))]
    @test length(prices) == nrow(df)
    @test model_types == df_types
  end
end


text:
using RNGPool
using CoupledConditionalSMC
import Statistics: mean, var, quantile
using JLD2

include("simpleModel.jl")

if length(ARGS) < 3
    n_rep = 100
    ofname = "simpleDemoCompareCost.jld2"
    task_id = 1
elseif length(ARGS) == 3
    n_rep = parse(Int64, ARGS[1])
    ofname = ARGS[2]
    task_id = parse(Int64, ARGS[3])
end

model, lM = setupSimpleModel(50000, 5.0)
setRNGs(12345 + task_id*n_rep)

function runDemo(algorithm, n::Int64, N::Int64, m::Int64;
    maxcost=1_000_000_000, maxmem::Int64=2^23)
    #N_ = Int(floor(min(n*N, maxmem)/n))
    #if N_<N
    #    println("Warning: Using N=", N_, " instead of N=", N)
    #end
    N_ = N
    maxiter = Int(floor(maxcost/N_/n))
    if algorithm == :BS || algorithm == :AS
        results = CoupledConditionalSMC.couplingTimes(model, lM, N_, n, m,
                                                    algorithm, true, true, maxiter)
    else
        results = CoupledConditionalSMC.couplingTimes(model, N_, n, m,
                                                      true, true, maxiter)
    end
    if all(ismissing.(results))
        results_ = ones(m)*Inf
    else
        results_ = results
    end
    #results.*N_/n
    results_, maxiter
end

# Time series lengths
ts = collect(1000:1000:5000)
Ns = 2 .^ (8:13)

BS_cost = zeros(length(ts), length(Ns), n_rep)
AS_cost = deepcopy(BS_cost)
maxIter = zeros(length(ts), length(Ns))

t0 = time()
#Threads.@threads
for i in length(ts):-1:1
  t = ts[i]
  for j in 1:length(Ns)
    N = Ns[j]
    println("BS: T=", t, ", N=", N)
    # BS estimator with given number of particles:
    BS_cost[i,j,:], _ = runDemo(:BS, t, N, n_rep)
    GC.gc()
    # AT estimator with t*N/10 particles:
    println("AS: T=", t, ", N=", N)
    AS_cost[i,j,:], maxIter[i,j] = runDemo(:AS, t, N, n_rep)
    GC.gc()
  end
end

elapsedTime = time()-t0

@save ofname BS_cost AS_cost ts n_rep Ns elapsedTime maxIter


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Geometry{Dim,T}

A geometry embedded in a `Dim`-dimensional space with coordinates of type `T`.
"""
abstract type Geometry{Dim,T} end

"""
    embeddim(geometry)

Return the number of dimensions of the space where the `geometry` is embedded.
"""
embeddim(::Type{<:Geometry{Dim,T}}) where {Dim,T} = Dim
embeddim(g::Geometry) = embeddim(typeof(g))

"""
    paramdim(geometry)

Return the number of parametric dimensions of the `geometry`. For example, a
sphere embedded in 3D has 2 parametric dimension (polar and azimuthal angles).
"""
paramdim(g::Geometry) = paramdim(typeof(g))

"""
    coordtype(geometry)

Return the machine type of each coordinate used to describe the `geometry`.
"""
coordtype(::Type{<:Geometry{Dim,T}}) where {Dim,T} = T
coordtype(g::Geometry) = coordtype(typeof(g))

"""
    centroid(geometry)

Return the centroid of the `geometry`.
"""
centroid(g::Geometry) = center(g)

"""
    measure(geometry)

Return the measure of the `geometry`, i.e. the length, area, or volume.
"""
function measure end

"""
    boundary(geometry)

Return the boundary of the `geometry`.
"""
function boundary end

# ----------------
# IMPLEMENTATIONS
# ----------------
include("primitives.jl")
include("polytopes.jl")


text:
#AuxiliaryData/{...}/{...}/{TAG}
function asdf_getaux( hdf_out::String )
  if isfile(hdf_out)
    io = h5open(hdf_out, "r+")
    fmt = read(attrs(io)["file_format"])
    (fmt == "ASDF") || (close(io); error("invalid ASDF file!"))
    if has(io, "AuxiliaryData")
      aux = io["AuxiliaryData"]
    else
      aux = g_create(io, "AuxiliaryData")
    end
  else
    io = h5open(hdf_out, "cw")
    attrs(io)["file_format"] = "ASDF"
    attrs(io)["file_format_version"] = "1.0.2"
    aux = g_create(io, "AuxiliaryData")
  end
  return io, aux
end

"""
    asdf_waux(hdf_out, path, X)

Write `X` to AuxiliaryData/path in `hdf_out`. If an object already exists at
AuxiliaryData/path, it will be deleted and overwritten with `X`.
"""
function asdf_waux(hdf_out::String, path::String, X::Union{HDF5Type,HDF5Array})
  # Correct leading /
  while startswith(path, "/")
    path = path[nextind(path, 1):lastindex(path)]
  end

  # Correct paths that start with AuxiliaryData
  startswith(path, "AuxiliaryData") && (path = String(split(path, "/", limit=2, keepempty=true)[2]))

  io, aux = asdf_getaux(hdf_out)

  # Remove existing object if it exists
  has(aux, path) && o_delete(aux, path)
  aux[path] = X
  close(io)
  return nothing
end


text:
# This is the Module for Code output in different style
# Specified for Genshin Impact, Types are Float32

module FileCodeOutput
using Dates

export Print2DMatrix,Print3DMatrix,PrintScriptVersion

# syntax dictionary
dictLanguage=Dict(
    "C"=>Dict(
    "termsuffix"=> "f",
    "termSplitter" => ", ",
    "lineStart" => "{",
    "lineEnd" => "},",
    "matrixStart" => "{",
    "matrixEnd" => "}",
    "comment" => "//",
    "terminator" => ";",
    ),
    "Julia"=>Dict(
    "termsuffix" => "f0",
    "termSplitter" => " ",
    "lineStart" => "",
    "lineEnd" => ";",
    "matrixStart" => "[",
    "matrixEnd" => "]",
    "comment" => "#",
    "terminator" => ";",
    ),
);

# Function to print Basic infomation
function PrintScriptVersion(f,s::String,v::String;codeStyle="C")
    print(f,dictLanguage[codeStyle]["comment"]*" This is a piece of code generated by hypergraph\\GenshinDataConverter\n");
    print(f,dictLanguage[codeStyle]["comment"]*" Code Style: "*codeStyle*"\n");
    print(f,dictLanguage[codeStyle]["comment"]*" Script Version: "*s*"\n");
    print(f,dictLanguage[codeStyle]["comment"]*" Your Version: "*v*"\n");
    print(f,dictLanguage[codeStyle]["comment"]*" Time of Code Generation: "*Dates.format(now(), "yyyy-mm-dd HH:MM")*"\n");
end

# Function to print 2D Data
function Print2DMatrix(f,M;tabLen=0,varName="",codeStyle="C")
    global dictLanguage

    # print variableName if needed
    if varName!=""
        if codeStyle=="C"
            print(f, "float "*varName);
            for i in 1:length(size(M))
                print(f,"["*string(size(M)[i])*"]");
            end
            print(f,"=");

        elseif codeStyle=="Julia"
            print(f,varName*"=");

        else
            error("Keyword codeStyle unsupported!")
        end
    end

    # print the mainpart of 2D Matrix
    print(f, repeat(" ",tabLen*4) * dictLanguage[codeStyle]["matrixStart"] * "\n");
    for i in 1:size(M)[1]
        print(f, repeat(" ", (tabLen+1)*4) * dictLanguage[codeStyle]["lineStart"]);
        for j in 1:size(M)[2]
            print(f, string(M[i,j]) * dictLanguage[codeStyle]["termsuffix"] * dictLanguage[codeStyle]["termSplitter"]);
        end
        print(f, dictLanguage[codeStyle]["lineEnd"] * "\n");
    end
    print(f, repeat(" ",tabLen*4) * dictLanguage[codeStyle]["matrixEnd"]);

    # print Terminator
    (varName!="")&&print(f, dictLanguage[codeStyle]["terminator"]);
    (varName=="")&&print(f, dictLanguage[codeStyle]["termSplitter"]);
    print(f,"\n");

end

# Function to print 3D Data
function Print3DMatrix(f,M;tabLen=0,varName="",codeStyle="C")
    global dictlanguage

    # print variableName if needed
    if varName!=""
        if codeStyle=="C"
            print(f, "float "*varName);
            for i in 1:length(size(M))
                print(f,"["*string(size(M)[i])*"]");
            end
            print(f,"=");

        elseif codeStyle=="Julia"
            print(f,varName*"=");

        else
            error("Keyword codeStyle unsupported!")
        end
    end

    # print the main part in different codeStyle
    if codeStyle=="C"
        print(f,"\n");
        print(f,repeat(" ",tabLen*4) * "{\n");
            for i in 1:size(M)[1]
                Print2DMatrix(f,M[i,:,:];tabLen=tabLen+1,codeStyle="C")
            end
        print(f,"}");

    elseif codeStyle=="Julia"
        print(f,"fill(0f0,"*string(size(M)[1])*","*string(size(M)[2])*","*string(size(M)[3])*");\n");
        for i in 1:size(M)[1]
            Print2DMatrix(f,M[i,:,:];varName=varName*"["*string(i)*",:,:]",codeStyle="Julia");
        end

    else
        error("Keyword codeStyle unsupported!")
    end

    # print Terminator
    (varName!="")&&(codeStyle=="C")&&print(f, dictLanguage[codeStyle]["terminator"]);
    print(f,"\n");
end

# Function to print a Sturcture/Dictionary or such
# function PrintDict()
# end

end


text:
module Simple

using Purple
using SymbolicUtils: Sym

b(x::Int) = 10 + x
function fn(x::Int, y::Int)
    z = 0
    z += x + y
    return b(z) + 10
end

println("Before:")
src = lift(fn,
           Tuple{Int, Int};
           clean = true,
           jit = false,
           opt = false)
display(src)

println("After:")
src = lift(fn,
           Tuple{Int, Int};
           jit = false,
           opt = false)
display(src)

f = lift(fn,
         Tuple{Int, Int};
         jit = true,
         opt = true)
ret = f(Sym{Int}(:x), Sym{Int}(:y))
display(ret)

end # module


text:
macro check_ptrs(args...)
    ex = Expr(:block)
    for x in args
        @assert x isa Symbol "$x should be a symbol."
        cond = :($(esc(x)).ptr != C_NULL)
        info = :(string(typeof($(esc(x)))) * " has a NULL pointer")
        push!(ex.args, Expr(:macrocall, Symbol("@assert"), nothing, cond, info))
    end
    return ex
end


text:
#Function to return valid intersections of lines with image boundary

function get_valid_intersections(intersections::Vector{Tuple{T,U}}, indsx::AbstractUnitRange, indsy::AbstractUnitRange) where {T<:Real, U<:Real}
    valid_intersections = Vector{Tuple{T,U}}(undef, 0)
    for intersection in intersections
        if first(indsx) <= intersection[1] <= last(indsx) && first(indsy) <= intersection[2] <= last(indsy)
            push!(valid_intersections,intersection)
        end
    end
    valid_intersections
end

# LineTwoPoints methods

LineTwoPoints(x0::Int, y0::Int, x1::Int, y1::Int) = LineTwoPoints(Point(x0, y0), Point(x1,y1))
LineTwoPoints(p1::CartesianIndex{2}, p2::CartesianIndex{2}) = LineTwoPoints(Point(p1), Point(p2))

draw!(img::AbstractArray{T,2}, line::LineTwoPoints, method::Function = bresenham) where {T<:Colorant} =
    draw!(img, line, oneunit(T), method)

function draw!(img::AbstractArray{T,2}, line::LineTwoPoints, color::T, method::Function = bresenham) where T<:Colorant
    indsy, indsx = axes(img)
    x1 = line.p1.x; y1 = line.p1.y
    x2 = line.p2.x; y2 = line.p2.y
    m = (y2 - y1) / (x2 - x1)
    intersections_x = [(x, y1 + m*(x-x1)) for x in (first(indsx), last(indsx))]
    intersections_y = [(x1 + (y-y1)/m, y) for y in (first(indsy), last(indsy))]
    valid_intersections = get_valid_intersections(vcat(intersections_x, intersections_y), indsx, indsy)
    if length(valid_intersections) > 0
        method(img, round(Int,valid_intersections[1][2]), round(Int,valid_intersections[1][1]), round(Int,valid_intersections[2][2]), round(Int,valid_intersections[2][1]), color)
    else
        img
    end
end

# LineNormal methods

LineNormal(τ::Tuple{T,U}) where {T<:Real, U<:Real} = LineNormal(τ...)

draw!(img::AbstractArray{T,2}, line::LineNormal, method::Function = bresenham) where {T<:Colorant} =
    draw!(img, line, oneunit(T), method)

function draw!(img::AbstractArray{T, 2}, line::LineNormal, color::T, method::Function = bresenham) where T<:Colorant
    indsy, indsx = axes(img)
    cosθ = cos(line.θ)
    sinθ = sin(line.θ)
    intersections_x = [(x, (line.ρ - x*cosθ)/sinθ) for x in (first(indsx), last(indsx))]
    intersections_y = [((line.ρ - y*sinθ)/cosθ, y) for y in (first(indsy), last(indsy))]
    valid_intersections = get_valid_intersections(vcat(intersections_x, intersections_y), indsx, indsy)
    if length(valid_intersections) > 0
        method(img, round(Int,valid_intersections[1][2]), round(Int,valid_intersections[1][1]), round(Int,valid_intersections[2][2]), round(Int,valid_intersections[2][1]), color)
    else
        img
    end
end

# LineSegment methods

LineSegment(x0::Int, y0::Int, x1::Int, y1::Int) = LineSegment(Point(x0, y0), Point(x1,y1))
LineSegment(p1::CartesianIndex, p2::CartesianIndex) = LineSegment(Point(p1), Point(p2))

draw!(img::AbstractArray{T,2}, line::LineSegment, method::Function = bresenham) where {T<:Colorant} =
    draw!(img, line, oneunit(T), method)

draw!(img::AbstractArray{T,2}, line::LineSegment, color::T, method::Function = bresenham) where {T<:Colorant} =
    method(img, line.p1.y, line.p1.x, line.p2.y, line.p2.x, color)

# Methods to draw lines

function bresenham(img::AbstractArray{T, 2}, y0::Int, x0::Int, y1::Int, x1::Int, color::T) where T<:Colorant
    dx = abs(x1 - x0)
    dy = abs(y1 - y0)

    sx = x0 < x1 ? 1 : -1
    sy = y0 < y1 ? 1 : -1;

    err = (dx > dy ? dx : -dy) / 2

    while true
        drawifinbounds!(img, y0, x0, color)
        (x0 != x1 || y0 != y1) || break
        e2 = err
        if e2 > -dx
            err -= dy
            x0 += sx
        end
        if e2 < dy
            err += dx
            y0 += sy
        end
    end

    img
end

fpart(pixel::T) where {T} = pixel - T(trunc(pixel))
rfpart(pixel::T) where {T} = oneunit(T) - fpart(pixel)

function swap(x, y)
    y, x
end

function xiaolin_wu(img::AbstractArray{T, 2}, y0::Int, x0::Int, y1::Int, x1::Int, color::T) where T<:Gray
    dx = x1 - x0
    dy = y1 - y0

    swapped=false
    if abs(dx) < abs(dy)
        x0, y0 = swap(x0, y0)
        x1, y1 = swap(x1, y1)
        dx, dy = swap(dx, dy)
        swapped=true
    end
    if x1 < x0
        x0, x1 = swap(x0, x1)
        y0, y1 = swap(y0, y1)
    end
    gradient = dy / dx

    xend = round(Int, x0)
    yend = y0 + gradient * (xend - x0)
    xgap = rfpart(x0 + 0.5)
    xpxl0 = xend
    ypxl0 = trunc(Int, yend)
    index = swapped ? CartesianIndex(xpxl0, ypxl0) : CartesianIndex(ypxl0, xpxl0)
    drawifinbounds!(img, index, T(rfpart(yend) * xgap))
    index = swapped ? CartesianIndex(xpxl0, ypxl0 + 1) : CartesianIndex(ypxl0 + 1, xpxl0)
    drawifinbounds!(img, index, T(fpart(yend) * xgap))
    intery = yend + gradient

    xend = round(Int, x1)
    yend = y1 + gradient * (xend - x1)
    xgap = fpart(x1 + 0.5)
    xpxl1 = xend
    ypxl1 = trunc(Int, yend)
    index = swapped ? CartesianIndex(xpxl1, ypxl1) : CartesianIndex(ypxl1, xpxl1)
    drawifinbounds!(img, index, T(rfpart(yend) * xgap))
    index = swapped ? CartesianIndex(xpxl1, ypxl1 + 1) : CartesianIndex(ypxl1 + 1, xpxl1)
    drawifinbounds!(img, index, T(fpart(yend) * xgap))

    for i in (xpxl0 + 1):(xpxl1 - 1)
        index = swapped ? CartesianIndex(i, trunc(Int, intery)) : CartesianIndex(trunc(Int, intery), i)
        drawifinbounds!(img, index, T(rfpart(intery)))
        index = swapped ? CartesianIndex(i, trunc(Int, intery) + 1) : CartesianIndex(trunc(Int, intery) + 1, i)
        drawifinbounds!(img, index, T(fpart(intery)))
        intery += gradient
    end
    img
end


text:
ds = DisplaySimulator(max_steps=10,
                      extra_initial=true,
                      extra_final=true,
                      rng=MersenneTwister(4))
m = BabyPOMDP()
@test simulate(ds, m, Starve()) ≈ 0.0

ds = DisplaySimulator(max_steps=1,
                      extra_initial=true,
                      extra_final=true,
                      rng=MersenneTwister(4))
m = SimpleGridWorld()
@test simulate(ds, m, FunctionPolicy(s->first(actions(m)))) ≈ 0.0


text:
# can we write a function adjusting the states so that process becomes monotone?
# resulting tree should be very close (in nested distance) to original process

#Example 1
trr1 = Tree([1,2,2,2]); trr2 = Tree([1,2,2,2]);
trr1.state = [0,0,0,0,0,0,0,1,3,20,40,2,4,21,41];
trr2.state = [0,0,0,0,0,0,0,1,21,3,41,2,20,4,40];

trr1.probability = [1.0,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5];
trr2.probability = trr1.probability;

nestedWasserstein(trr1,trr2,1.0)      # supposed to be 8.75

#Example 2
trr1 = Tree([1,2,4,1]); trr2 = Tree([1,2,4,1]);
trr1.state = [0,0,0,0,0,0,0,0,0,0,0,1,3,20,40,2,4,21,41];
trr2.state = [0,0,0,0,0,0,0,0,0,0,0,1,21,3,41,2,20,4,40];

trr1.probability = [1.0,0.5,0.5,0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0];
trr2.probability = trr1.probability;

nestedWasserstein(trr1,trr2,1.0)

# Example 3

trr1 = Tree([1,1,2]); trr2 = Tree([1,2,1]);
trr1.state = [2,2,3,1];
trr2.state = [2,2,2,3,1];

trr1.probability = [1.0,1.0,0.3,0.7];
trr2.probability = [1.0,0.3,0.7,1.0,1.0];

nestedWasserstein(trr1,trr2,1.0)


text:
export
    HOBO

import Base: *

"""
HOBO classes. Current HOBO object should not contain zero(T) values, which may result
in unnecessary memory blow.
"""
struct HOBO{T} <: BinaryOptimizationModel{T}
    vars::Variables
    values::Dict{Vector{Vector{Int}},T}
    function HOBO{T}(vars::Variables,
                     values::Dict{Vector{Vector{Int}},T}) where T<:Real
        final_values = Dict{Vector{Vector{Int}},T}()

        for (key, val) = values
            final_key = clean_key(HOBO, key)
            @assert all(length(i) == dim(vars) for i = key)
            @assert all(all(0 .< ind .<= indexcapacities(vars)) for ind=final_key)
            _dict_add!(final_values, final_key, val)
        end
        new{T}(vars, final_values)
    end
end

function HOBO(::Type{T}, vars::Variables) where T
    HOBO{T}(vars, Dict{Vector{Vector{Int}},T}())
end

HOBO(vars::Variables) = HOBO(Float64, vars::Variables)
HOBO{T}(vars::Variables) where T<:Real = HOBO(T, vars::Variables)

function HOBO(num::T, vars::Variables) where T<:Real
    HOBO{T}(vars, Dict(Vector{Int}[] => num))
end

clean_key(::Type{HOBO}, key::Vector{Vector{Int}}) = sort(unique(key))
clean_key(::Type{HOBO{T}}, key::Vector{Vector{Int}}) where T<:Real = clean_key(HOBO, key)

## may be improved - smart tree walker
function (hobo::HOBO{T})(keys_list::Set{Vector{Int}}) where T<:Real
    result = zero(T)
    for (key, val) = values(hobo)
        if all(k ∈ keys_list for k = key)
            result += val
        end
    end
    result
end

#(hobo::HOBO{<:Real})(keys_list) = (hobo::HOBO)(Set(keys_list...))
#(hobo::HOBO{<:Real})(keys_list::Vector{Vector{Int}}) = (hobo::HOBO)([keys_list])

## logical functions

or(hobo1::HOBO, hobo2::HOBO) = hobo1 + hobo2 - hobo1*hobo2
xor(hobo1::HOBO{T}, hobo2::HOBO{T}) where T<:Real = hobo1 + hobo2 - T(2)*hobo1*hobo2
and(hobo1::HOBO, hobo2::HOBO) = hobo1 * hobo2
neg(hobo::HOBO{T}) where T = one(T) - hobo


text:
function poly(x, c)
    out = x .* 0. .+ c[1]
    for i in 2:length(c)
        out .+= c[i] .* x.^(i-1)
    end
    return out
end

function ccm_unred(wave, ebv; R_V=3.1)
    x = 10000 ./ wave    # Convert to inverse microns
    npts = length( x )
    a = fill(0., npts)
    b = fill(0., npts)

    good = findall((x .>  0.3)  .&   (x .< 1.1)) #Infrared
    a[good] =  0.574 .* x[good].^(1.61)
    b[good] = -0.527 .* x[good].^(1.61)

    good = findall((x .>= 1.1)  .&   (x .< 3.3)) #Optical/NIR
    #Use new constants from O'Donnell (1994)
    y = x[good] .- 1.82
#     c1 = [ 1. , 0.17699, -0.50447, -0.02427,  0.72085,      #Original
#                 0.01979, -0.77530,  0.32999 ]               #coefficients
#     c2 = [ 0.,  1.41338,  2.28305,  1.07233, -5.38434,      #from CCM89
#                -0.62251,  5.30260, -2.09002 ]
      c1 = [ 1. , 0.104,   -0.609,    0.701,  1.137,         #New coefficients
                 -1.718,   -0.827,    1.647, -0.505 ]        #from O'Donnell
      c2 = [ 0.,  1.952,    2.908,   -3.989, -7.985,         #(1994)
                 11.102,    5.491,  -10.805,  3.347 ]
     a[good] = poly(y, c1)
     b[good] = poly(y, c2)

    good = findall((x .>= 3.3)  .&   (x .< 8.0)) #Mid-UV
    y = x[good]
    F_a = fill(0., length(good))
    F_b = fill(0., length(good))
    good1 = findall(y .> 5.9)
    y1 = y[good1] .- 5.9
    F_a[good1] = -0.04473 .* y1.^2 .- 0.009779 .* y1.^3
    F_b[good1] =   0.2130 .* y1.^2 .+   0.1207 .* y1.^3

    a[good] =  1.752 .- 0.316.*y .- (0.104 ./ ( (y.-4.67).^2 .+ 0.341 )) .+ F_a
    b[good] = -3.090 .+ 1.825.*y .+ (1.206 ./ ( (y.-4.62).^2 .+ 0.263 )) .+ F_b

    good = findall((x .>= 8.0)  .&   (x .< 11.)) #Far-UV
    y = x[good] .- 8.
    c1 = [ -1.073, -0.628,  0.137, -0.070 ]
    c2 = [ 13.670,  4.257, -0.420,  0.374 ]
    a[good] = poly(y, c1)
    b[good] = poly(y, c2)

    # Now apply extinction correction to input flux vector
    A_V = R_V * ebv
    A_lambda = A_V .* (a .+ b./R_V)
    return 10 .^(0.4.*A_lambda)
 end