text:
"""
Function to shrink acceptbw by weighted quantile of previously accepted distances

# Arguments
* `distances::Array{Float64, 1}`: vector of distances or vector of distance vectors
* `weights::WeightVec`: weight vector must be in WeightVec format
* `acceptbwQuantile::Float64`: quantile (0,1) to shrink to

# Value
acceptbw acceptance acceptbws
"""
function quantile_acceptbw(d::Array{<:AbstractFloat, 1}, w::StatsBase.AnalyticWeights, q::AbstractFloat)
    return StatsBase.quantile(d, w, q)
end

function quantile_acceptbw{A <: SingleMeasureAbcPmc}(abcpmc::A, q::Float64)
    return StatsBase.quantile(abcpmc.distances, abcpmc.weights, q)
end

function independentquantile_acceptbw(d::Array{<:AbstractFloat, 2}, w::StatsBase.AnalyticWeights, q::AbstractFloat)
    m = size(d, 2)
    acceptbws = Array{Float64}(m)
    for jj in 1:m
        acceptbws[jj] = quantile(d[:, jj], w, q)
    end
    return acceptbws
end

function independentquantile_acceptbw(abcpmc::A, q::Float64) where A <: MultiMeasureAbcPmc
    return independentquantile_acceptbw(abcpmc.distances, abcpmc.weights, q)
end

function jointquantile_acceptbw(d::Array{<:AbstractFloat, 2}, w::StatsBase.AnalyticWeights, q::AbstractFloat)
    m = size(d, 2)
    return independentquantile_acceptbw(d, w, 1.0 - (1.0 - q) / m)
end

function jointquantile_acceptbw(abcpmc::A, q::Float64) where A <: MultiMeasureAbcPmc
    m = size(abcpmc.distances, 2)
    return independentquantile_acceptbw(abcpmc, 1.0 - (1.0 - q) / m)
end


text:
export getchar, set_homedir, get_homedir

charset = 'a':'z'

getchar(idx::Int64) = charset[idx]::Char

const homedir = Ref{String}("wrong")

function set_homedir(path)
    global homedir[] = path
end
get_homedir() = homedir[]

ext: jl


text:
function Sim_Model(CellData,Iapp,Tk,A0,B0,C0,D0)
    """
    Simulation of generated reduced-order models
    # Add License
    # Add Ins and Outs
        # Locations to be computed
        # Sampling Frequency
    """
    #Determine time span and allocate arrays
    tlength = size(Iapp,1)

    CellV_= Array{Float64}(undef,tlength,0)
    Ce_= Array{Float64}(undef,tlength,0)
    j0_ = Array{Float64}(undef,tlength,0)
    RtotNeg_ = Array{Float64}(undef,tlength,0)
    RtotPos_ = Array{Float64}(undef,tlength,0)
    ϕ_ẽ1_ = Array{Float64}(undef,tlength,0)
    ϕ_ẽ2_ = Array{Float64}(undef,tlength,0)
    Uocp_Neg_ = Array{Float64}(undef,tlength,0)
    Uocp_Pos_ = Array{Float64}(undef,tlength,0)
    η0_ = Array{Float64}(undef,tlength,0)
    η_neg_ = Array{Float64}(undef,tlength,0)
    ηL_ = Array{Float64}(undef,tlength,0)
    η_pos_ = Array{Float64}(undef,tlength,0)
    ϕ_e_ = Array{Float64}(undef,tlength,0)
    jNeg_ = Array{Float64}(undef,tlength,0)
    jPos_ = Array{Float64}(undef,tlength,0)
    Cse_Neg_ = Array{Float64}(undef,tlength,0)
    Cse_Pos_ = Array{Float64}(undef,tlength,0)

    #Selecting SS Models
    for γ in 1:9
        A = A0[γ]
        B = B0[γ]
        C = C0[γ]
        D = D0[γ]

        #Capturing Indices
        tfstr = Array{String}(undef,0,1)
        for i in 1:size(CellData.Transfer.tfs[:,1],1)
            t1 = t2 = Array{String}(undef,0,1)
            for j in 1:size(CellData.Transfer.tfs[i,3],1)
                t1 ="$(CellData.Transfer.tfs[i,1])_$(CellData.Transfer.tfs[i,2])"
                t2 = [t2; t1]
            end
            tfstr = [tfstr; t2]
        end

        CeInd = findall(isequal("C_e_Na"), tfstr)
        ϕ_ẽInd = findall(isequal("Phi_e_Na"), tfstr)
        CsePosInd = findall(isequal("C_se_Pos"), tfstr)
        CseNegInd = findall(isequal("C_se_Neg"), tfstr)
        ϕ_sNegInd = findall(isequal("Phi_s_Neg"), tfstr)
        ϕ_sPosInd = findall(isequal("Phi_s_Pos"), tfstr)
        ϕ_sePosInd = findall(isequal("Phi_se_Pos"), tfstr)
        ϕ_seNegInd = findall(isequal("Phi_se_Neg"), tfstr)
        FluxNegInd = findall(isequal("Flux_Neg"), tfstr)
        FluxPosInd = findall(isequal("Flux_Pos"), tfstr)

        CeNeg = CellData.Transfer.tfs[1,3].<=CellData.Neg.L
        CeNegOffset = CellData.Transfer.tfs[1,3].<CellData.Neg.L
        CeSep = CellData.Transfer.tfs[1,3].<=CellData.Neg.L+CellData.Sep.L
        CeSepOffset = CellData.Transfer.tfs[1,3].<CellData.Neg.L+CellData.Sep.L
        CePos = CellData.Transfer.tfs[1,3].<=CellData.Const.Ltot
        CeNegInd = findall(CeNeg .== 1)
        CeSepInd = findall(CeSep.-CeNegOffset .== 1)
        CePosInd = findall(CePos.-CeSepOffset .== 1)

        csegain_neg = C[CseNegInd[1][1],end]
        csegain_pos = C[CsePosInd[1][1],end]

        #Memory Allocation
        θ_neg = Array{Float64}(undef,tlength,1) .= 0.
        θ_pos = Array{Float64}(undef,tlength,1) .= 0.
        x = Array{Float64}(undef,tlength,size(A,1)) .= 0.
        y = Array{Float64}(undef,tlength,size(C,1)) .= 0.
        Cse_Neg = Array{Float64}(undef,tlength,size(CseNegInd,1)) .= 0.
        Cse_Pos = Array{Float64}(undef,tlength,size(CsePosInd,1)) .= 0.
        Ce = Array{Float64}(undef,tlength,size(CeInd,1)) .= CellData.Const.ce0
        η0 = Array{Float64}(undef,tlength,1) .= 0.
        η_neg = Array{Float64}(undef,tlength,size(FluxNegInd,1)) .= 0.
        ηL = Array{Float64}(undef,tlength,1) .= 0.
        η_pos = Array{Float64}(undef,tlength,size(FluxPosInd,1)) .= 0.
        ϕ_ẽ1 = Array{Float64}(undef,tlength,size(ϕ_ẽInd,1)) .= 0.
        ϕ_ẽ2 = Array{Float64}(undef,tlength,size(CeInd,1)) .= 0.
        ϕ_se_neg_0 = Array{Float64}(undef,tlength,1) .= 0. #Replace with length of ϕ_seNegInd @ zero
        jNeg = Array{Float64}(undef,tlength,size(FluxNegInd,1)) .= 0.
        j0 = Array{Float64}(undef,tlength,1) .= 0.
        jPos = Array{Float64}(undef,tlength,size(FluxPosInd,1)) .= 0.
        jL = Array{Float64}(undef,tlength,1) .= 0.
        j0_CC_neg = Array{Float64}(undef,tlength,1) .= 0.
        Rtot_neg = Array{Float64}(undef,tlength,1) .= 0.
        Rtot_pos = Array{Float64}(undef,tlength,1) .= 0.
        j0_CC_pos = Array{Float64}(undef,tlength,1) .= 0.
        Uocp_Neg = Array{Float64}(undef,tlength,1) .= 0.
        Uocp_Pos = Array{Float64}(undef,tlength,1) .= 0.
        Cell_V = Array{Float64}(undef,tlength,1) .= 0.
        ϕ_e = Array{Float64}(undef,tlength,size(CeInd,1)) .= 0.

        #Defining SOC
        SOC_Neg = CellData.Const.SOC * (CellData.Neg.θ_100-CellData.Neg.θ_0) + CellData.Neg.θ_0
        SOC_Pos = CellData.Const.SOC * (CellData.Pos.θ_100-CellData.Pos.θ_0) + CellData.Pos.θ_0
        θ_neg[1] = SOC_Neg
        θ_pos[1] = SOC_Pos

        #Loop through time
        #Compute dependent variables (voltage, flux, etc.)
        for i in 1:tlength-1
            cs_neg_avg = x[i,end] * csegain_neg + SOC_Neg * CellData.Neg.cs_max < 0. ? 0. : x[i,end] * csegain_neg + SOC_Neg * CellData.Neg.cs_max #Zero if < 0
            cs_pos_avg = x[i,end] * csegain_pos + SOC_Pos * CellData.Pos.cs_max < 0. ? 0. : x[i,end] * csegain_pos + SOC_Pos * CellData.Pos.cs_max #Zero if < 0

            if cs_neg_avg > CellData.Neg.cs_max
                cs_neg_avg = CellData.Neg.cs_max
            end

            if cs_pos_avg > CellData.Pos.cs_max
                cs_pos_avg = CellData.Pos.cs_max
            end

            #Reaction Rates
            if CellData.Const.CellTyp == "Doyle_94"
                k_neg = CellData.Neg.k_norm/CellData.Neg.cs_max/CellData.Const.ce0^(1-CellData.Neg.α)
                k_pos = CellData.Pos.k_norm/CellData.Pos.cs_max/CellData.Const.ce0^(1-CellData.Pos.α)
                jeq_neg = k_neg*sqrt(cs_neg_avg*CellData.Const.ce0*(CellData.Neg.cs_max-cs_neg_avg))
                jeq_pos = k_pos*sqrt(cs_pos_avg*CellData.Const.ce0*(CellData.Pos.cs_max-cs_pos_avg))
            else
                jeq_neg = CellData.Neg.k_norm*CellData.Neg.cs_max*(CellData.Const.ce0*(cs_neg_avg/CellData.Neg.cs_max*(1-cs_neg_avg/CellData.Neg.cs_max)))^(1-CellData.Neg.α)
                jeq_pos = CellData.Pos.k_norm*CellData.Pos.cs_max*(CellData.Const.ce0*(cs_pos_avg/CellData.Pos.cs_max*(1-cs_pos_avg/CellData.Pos.cs_max)))^(1-CellData.Pos.α)
            end

            θ_neg[i+1] = cs_neg_avg/CellData.Neg.cs_max
            θ_pos[i+1] = cs_pos_avg/CellData.Pos.cs_max
            Cell_SOC = (θ_neg[i]-CellData.Neg.θ_0)/(CellData.Neg.θ_100-CellData.Neg.θ_0)

            javg_neg = Iapp[i]/(CellData.Neg.as*F*CellData.Neg.L*CellData.Const.CC_A)
            javg_pos = Iapp[i]/(CellData.Pos.as*F*CellData.Pos.L*CellData.Const.CC_A)

            Arr_Factor = ((1/CellData.Const.T_ref)-(1/Tk[i]))/R
            κneg = CellData.Const.κf(mean(Ce[i,CeNegInd]))*exp(CellData.Const.Ea_κ*Arr_Factor)
            κpos = CellData.Const.κf(mean(Ce[i,CePosInd]))*exp(CellData.Const.Ea_κ*Arr_Factor) # Parameterise
            κsep = CellData.Const.κf(mean(Ce[i,CeSepInd]))*exp(CellData.Const.Ea_κ*Arr_Factor)
            σ_neg = CellData.Neg.σ*exp(CellData.Const.Ea_κ*Arr_Factor)
            σ_pos = CellData.Pos.σ*exp(CellData.Const.Ea_κ*Arr_Factor)
            κ_eff_Neg = κneg*(CellData.Neg.ϵ_e^(CellData.Neg.κ_brug))
            κ_eff_Sep = κsep*(CellData.Sep.ϵ_e^(CellData.Sep.κ_brug))
            κ_eff_Pos = κpos*(CellData.Pos.ϵ_e^(CellData.Pos.κ_brug))
            σ_eff_Neg = σ_neg*CellData.Neg.ϵ_s^CellData.Neg.σ_brug #Effective Conductivity Neg
            σ_eff_Pos = σ_pos*CellData.Pos.ϵ_s^CellData.Pos.σ_brug #Effective Conductivity Pos

            #Resistances
            Rtot_neg[i] = (Tk[i]*R)/(F^2*mean([abs(jeq_neg) abs(javg_neg)]))+CellData.Neg.RFilm
            Rtot_pos[i] = (Tk[i]*R)/(F^2*mean([abs(jeq_pos) abs(javg_pos)]))+CellData.Pos.RFilm

            #Condensing Variable
            ν_neg = @. CellData.Neg.L*sqrt((CellData.Neg.as*(1/κ_eff_Neg+1/σ_eff_Neg))/Rtot_neg[i])
            ν_pos = @. CellData.Pos.L*sqrt((CellData.Pos.as*(1/κ_eff_Pos+1/σ_eff_Pos))/Rtot_pos[i])

            #Relinearise dependent on ν, σ, κ
            #Call from CellData? List of functions composed from ROM creation?
            D = D_Linear(CellData, ν_neg, ν_pos, σ_eff_Neg, κ_eff_Neg, σ_eff_Pos, κ_eff_Pos, κ_eff_Sep)

            #SS Output
            y[i,:] = C*x[i,:] + D*Iapp[i]

            #Concentrations & Force electrode concentration maximum
            Cse_Neg[i,:] = (SOC_Neg.*CellData.Neg.cs_max .+ y[i,CseNegInd]) > ones(size(Cse_Neg,2))*CellData.Neg.cs_max ? ones(size(Cse_Neg,2))*CellData.Neg.cs_max : (SOC_Neg.*CellData.Neg.cs_max .+ y[i,CseNegInd])
            Cse_Pos[i,:] = (SOC_Pos.*CellData.Pos.cs_max .+ y[i,CsePosInd]) > ones(size(Cse_Neg,2))*CellData.Neg.cs_max ? ones(size(Cse_Neg,2))*CellData.Neg.cs_max : (SOC_Pos.*CellData.Pos.cs_max .+ y[i,CsePosInd])
            Ce[i,:] = @. CellData.Const.ce0 + y[i,CeInd]

            #Potentials
            Uocp_Neg[i] = CellData.Const.Uocp("Neg",Cse_Neg[i,1]/CellData.Neg.cs_max)
            Uocp_Pos[i] = CellData.Const.Uocp("Pos",Cse_Pos[i,1]/CellData.Pos.cs_max)
            ϕ_se_neg_0[i] = y[i, ϕ_seNegInd[1]] + Uocp_Neg[i] #Location 0
            ϕ_ẽ1[i,:] = y[i,ϕ_ẽInd]
            ϕ_ẽ2[i,:] = @. ((Tk[i]*2*R*(1-CellData.Const.tpf(Ce[i,:])))/F)*(log(Ce[i,:]/Ce[i,1]))
            ϕ_e[i,:] = @. [0; ϕ_ẽ1[i,:]]+ϕ_ẽ2[i,:]-ϕ_se_neg_0[i]

            #Flux
            jNeg[i,:] = y[i,FluxNegInd]
            j0[i] = y[i,FluxNegInd[1]]
            jPos[i,:] = y[i,FluxPosInd]
            jL[i] = y[i,FluxPosInd[1]]

            if CellData.Const.CellTyp == "Doyle_94"
                j0_CC_neg[i] = findmax([eps(); (k_neg*(CellData.Neg.cs_max-Cse_Neg[i,1])^(1-CellData.Neg.α))*((Cse_Neg[i,1]^CellData.Neg.α)*(Ce[i,1]^(1-CellData.Neg.α)))])[1]
                j0_neg = findmax([ones(size(Cse_Neg,2))*eps() (k_neg.*(Cse_Neg[i,:].^CellData.Neg.α).*(Ce[i,1].^(1-CellData.Neg.α))).*(CellData.Neg.cs_max.-Cse_Neg[i,:]).^(1-CellData.Neg.α)], dims=2)[1]
                η0[i] = Tk[i]*2*R/F*asinh(j0[i]/(2*j0_CC_neg[i]))
                η_neg[i,:] = @. (Tk[i]*2*R)/F*asinh((jNeg[i,:])/(2*j0_neg))

                j0_CC_pos[i] = findmax([eps(); (k_pos*(CellData.Pos.cs_max-Cse_Pos[i,1])^(1-CellData.Pos.α))*((Cse_Pos[i,1]^CellData.Pos.α)*(Ce[i,end]^(1-CellData.Pos.α)))])[1]
                j0_pos = findmax([ones(size(Cse_Pos,2))*eps() (k_pos.*(Cse_Pos[i,:].^CellData.Pos.α).*(Ce[i,1].^(1-CellData.Pos.α))).*(CellData.Pos.cs_max.-Cse_Pos[i,:]).^(1-CellData.Pos.α)], dims=2)[1]
                ηL[i] = (Tk[i]*2*R)/F*asinh(jL[i]/(2*j0_CC_pos[i]))
                η_pos[i,:] = @. (Tk[i]*2*R)/F*asinh(jPos[i,:]/(2*j0_pos))
            else
                j0_CC_neg[i] = findmax([eps(); CellData.Neg.k_norm*CellData.Neg.cs_max*(Ce[i,1]*(Cse_Neg[i,1]/CellData.Neg.cs_max*(1-Cse_Neg[i,1]/CellData.Neg.cs_max)))^(1-CellData.Neg.α)])[1]
                j0_neg = @. findmax([eps(); (CellData.Neg.k_norm*(Cse_Neg[i,:]^CellData.Neg.α)*(Ce[i,1]^(1-CellData.Neg.α)))*(CellData.Neg.cs_max-Cse_Neg[i,:])^(1-CellData.Neg.α)])[1]
                η0[i] = Tk[i]*2*R/F*asinh(j0[i])/(2*j0_CC_neg[i])
                η_neg[i,:] = @. (Tk[i]*2*R)/F*asinh(jNeg[i,:]/(2*j0_neg))

                j0_CC_pos[i] = findmax([eps(); CellData.Pos.k_norm*CellData.Pos.cs_max*(Ce[i,end]*(Cse_Pos[i,1]/CellData.Pos.cs_max*(1-Cse_Pos[i,1]/CellData.Pos.cs_max)))^(1-CellData.Pos.α)])[1]
                j0_pos = @. findmax([eps(); (CellData.Pos.k_norm*(Cse_Pos[i,:]^CellData.Pos.α)*(Ce[i,1]^(1-CellData.Pos.α)))*(CellData.Pos.cs_max-Cse_Pos[i,:])^(1-CellData.Pos.α)])[1]
                ηL[i] = (Tk[i]*2*R)/F*asinh(jL[i])/(2*j0_CC_pos[i])
                η_pos[i,:] = @. (Tk[i]*2*R)/F*asinh(jPos[i,:]/(2*j0_pos))
            end

            #Cell Voltage
            Cell_V[i] = @. (Uocp_Pos[i]-Uocp_Neg[i]) + (ηL[i]-η0[i]) + (ϕ_ẽ1[i,end]+ϕ_ẽ2[i,end]) + (CellData.Pos.RFilm*jL[i]-CellData.Neg.RFilm*j0[i])*F

            #ϕ_s
            ϕ_s_neg = y[i,ϕ_sNegInd]
            ϕ_s_pos = @. y[i,ϕ_sPosInd] + Cell_V[i]

            #Update States
            x[i+1,:] = A*x[i,:] + B*Iapp[i]
        end
        CellV_ = [CellV_ Cell_V]
        Ce_ = [Ce_ Ce]
        jNeg_ = jNeg
        jPos_ = jPos
        RtotNeg_ = [RtotNeg_ Rtot_neg]
        RtotPos_ = [RtotPos_ Rtot_pos]
        η0_ = [η0_ η0]
        η_neg_ = η_neg
        ηL_ = [ηL_ ηL]
        η_pos_ = η_pos
        ϕ_ẽ1_ = ϕ_ẽ1
        ϕ_ẽ2_ = ϕ_ẽ2
        ϕ_e_ = ϕ_e
        Uocp_Neg_ = [Uocp_Neg_ Uocp_Neg]
        Uocp_Pos_ = [Uocp_Pos_ Uocp_Pos]
        Cse_Neg_ = [Cse_Neg_ Cse_Neg]
        Cse_Pos_ = [Cse_Pos_ Cse_Pos]
    end
    return CellV_, Ce_, jNeg_, jPos_, RtotNeg_, RtotPos_, η0_, ηL_, η_neg_, η_pos_, ϕ_ẽ1_, ϕ_ẽ2_, Uocp_Neg_, Uocp_Pos_, ϕ_e_, Cse_Neg_, Cse_Pos_
end

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text: println("Hello minoJiro!")

max_stars_repo_name: minojiro/Hello-World

max_issues_repo_name: minojiro/Hello-World

max_forks_repo_name: minojiro/Hello-World


text:
"""
    ProofLine(line::Int, statement::AbstractExpression, argument::A, references::Vector{ProofLine}) where {A}

Represents a single line in a step-by-step proof, made up by a statement, an argument, and a set of references necessary to
make the provided argument. A line number is also required to serve as a unique human-readable identifier for each line.
"""
struct ProofLine{A}
    linenum::Int
    statement::AbstractExpression
    argument::A
    references::Vector{ProofLine{A}}
end
"""
    ProofLine(line::Int, statement::AbstractExpression, argument::A="N/A") where {A}

Creates a [`ProofLine`](@ref) without any references to other lines and without an argument. Implicitly this means in most proofs that
this line depends only on the previous one, but can also be taken to mean "is a conclusion of all statements above".
"""
function ProofLine(line::Int, statement::AbstractExpression, argument::A="N/A") where {A}
    ProofLine(line, statement, argument, ProofLine{A}[])
end
"""
    ProofLine(line::Int, statement::AbstractExpression, argument::A, reference::ProofLine{A}) where {A}

Convenience method. Creates a [`ProofLine`](@ref) with a single provided reference instead of a list of them.
"""
function ProofLine(line::Int, statement::AbstractExpression, argument::A, reference::ProofLine{A}) where {A}
    ProofLine(line, statement, argument, [reference])
end
# in most cases this shouldn't get shown since we also override Base.show(::IO, ::Vector{ProofLine})
# put another way, this is the best we can print out a proof line without context from the proof itself
function Base.show(io::IO, line::ProofLine)
    print(io, line.linenum)
    print(io, "\t")
    print(io, replace(string(line.statement), "Deductive." => ""))
    print(io, "\t")
    print(io, line.argument)

    if length(line.references) > 0
        print(io, "\t")
        print(io, "(" * join([string(ref.linenum) for ref ∈ line.references], ", ") * ")")
    end
end

function find_proof_line_by_statement(proof::Vector{ProofLine}, statement::AbstractExpression)
    for line ∈ proof
        if isequal(line.statement, statement)
            return line
        end
    end

    nothing
end

function Base.show(io::IO, m::MIME"text/plain", proof::Vector{ProofLine})
    proof_table = DataFrame("Line Number" => Int[], "Statement" => String[], "Argument" => String[], "References" => String[])

    for line ∈ proof
        push!(proof_table, Dict(
            "Line Number" => line.linenum,
            "Statement" => replace(string(line.statement), "Deductive." => ""),
            "Argument" => line.argument,
            "References" => join([string(ref.linenum) for ref ∈ line.references], ", ")
        ))
    end

    # display_size=(-1, -1) forces pretty_table to print all rows and columns regardless of screen size
    pretty_table(io, proof_table; display_size=(-1, -1))
end


text:
"
A pseudo-marginal Metropolis-Hastings algorithm which uses RE-ABC likelihood estimates.

## Arguments

- `a` Details of the application.
- `schedule` Sequence of epsilon values. Use this option for non-adaptive RE-ABC.
- `targetε` Target ε value. Use this option for adaptive RE-ABC.
- `N` Number of iterations to do.
- `nparticles` Number of particles to use.
- `prior` Prior distribution.
- `θinit` Initial state of chain.
- `inc` Distribution of proposal increments
- `k` Number of particles to accept at each iteration of adaptive RE-ABC.
"

function doPMMH(a::ApplicationDetails, schedule::Array{Float64, 1}, N::Int, nparticles::Int, prior::UnivariateDistribution, θinit::Float64, inc::UnivariateDistribution)
    states = zeros(Float64, N)
    loglikelihoods = zeros(Float64, N)
    states[1] = θ = θinit
    loglikelihoods[1] = ll = doSMC(a, [θ], schedule, nparticles=nparticles, quiet=true).llest
    logprior = logpdf(prior, θ)
    @showprogress for i in 2:N
        θprop = θ + rand(inc)
        logpriorProp = logpdf(prior, θprop)
        llBound = log(rand()) + ll + logprior - logpriorProp
        if llBound < Inf
            p = doSMC(a, [θprop], schedule, nparticles=nparticles, llBound=llBound, quiet=true)
            if !p.reachedBound ##Acceptance
                θ, ll = θprop, p.llest
            end
        end
        states[i], loglikelihoods[i] = θ, ll
    end
    states, loglikelihoods
end

function doPMMH(a::ApplicationDetails, schedule::Array{Float64, 1}, N::Int, nparticles::Int, prior::MultivariateDistribution, θinit::Array{Float64, 1}, inc::MultivariateDistribution)
    p = length(θinit)
    states = zeros(Float64, (N,p))
    loglikelihoods = zeros(Float64, N)
    states[1,:] = θ = copy(θinit)
    loglikelihoods[1] = ll = doSMC(a, θ, schedule, nparticles=nparticles, quiet=true).llest
    logprior = logpdf(prior, θ)
    @showprogress for i in 2:N
        θprop = θ + rand(inc)
        logpriorProp = logpdf(prior, θprop)
        llBound = log(rand()) + ll + logprior - logpriorProp
        if llBound < Inf
            p = doSMC(a, θprop, schedule, nparticles=nparticles, llBound=llBound, quiet=true)
            if !p.reachedBound ##Acceptance
                θ, ll = θprop, p.llest
            end
        end
        states[i,:], loglikelihoods[i] = θ, ll
    end
    states, loglikelihoods
end

function doPMMH(a::ApplicationDetails, targetε::Float64, N::Int, nparticles::Int, prior::UnivariateDistribution, θinit::Float64, inc::UnivariateDistribution; k=ceil(Int, nparticles/2))
    states = zeros(Float64, N)
    loglikelihoods = zeros(Float64, N)
    states[1] = θ = θinit
    loglikelihoods[1] = ll = doSMC(a, [θ], targetε, nparticles=nparticles, quiet=true, k=k).llest
    logprior = logpdf(prior, θ)
    @showprogress for i in 2:N
        θprop = θ + rand(inc)
        logpriorProp = logpdf(prior, θprop)
        llBound = log(rand()) + ll + logprior - logpriorProp
        if llBound < Inf
            p = doSMC(a, [θprop], targetε, nparticles=nparticles, llBound=llBound, quiet=true, k=k)
            if !p.reachedBound ##Acceptance
                θ, ll = θprop, p.llest
            end
        end
        states[i], loglikelihoods[i] = θ, ll
    end
    states, loglikelihoods
end

function doPMMH(a::ApplicationDetails, targetε::Float64, N::Int, nparticles::Int, prior::MultivariateDistribution, θinit::Array{Float64, 1}, inc::MultivariateDistribution; k=ceil(Int, nparticles/2))
    p = length(θinit)
    states = zeros(Float64, (N,p))
    loglikelihoods = zeros(Float64, N)
    states[1,:] = θ = copy(θinit)
    loglikelihoods[1] = ll = doSMC(a, θ, targetε, nparticles=nparticles, quiet=true, k=k).llest
    logprior = logpdf(prior, θ)
    @showprogress for i in 2:N
        θprop = θ + rand(inc)
        logpriorProp = logpdf(prior, θprop)
        llBound = log(rand()) + ll + logprior - logpriorProp
        if llBound < Inf
            p = doSMC(a, θprop, targetε, nparticles=nparticles, llBound=llBound, quiet=true, k=k)
            if !p.reachedBound ##Acceptance
                θ, ll = θprop, p.llest
            end
        end
        states[i,:], loglikelihoods[i] = θ, ll
    end
    states, loglikelihoods
end


text:
using EffectiveWaves
using LinearAlgebra

include("test_discrete_solution.jl")
# Set parameters
particle_medium = Acoustic(3; ρ=0.01, c=0.01);
particle_medium = Acoustic(3; ρ=10.0, c=10.0);
medium = Acoustic(3; ρ=1.0, c=1.0);

R = 8.0
r = 1.0

separation_ratio = 1.02

ka = 0.4
k = ka / r

vol_fraction = 0.1

basis_order = Int(round(4. * ka)) .+ 1
basis_field_order = Int(round(2.0 * k * R)) .+ 1
# basis_field_order = 6

# basis_order = 1;
# basis_field_order = 6;
# basis_field_order = 10;

ω = k * real(medium.c)

s1 = Specie(
    particle_medium, Sphere(r);
    number_density = vol_fraction / volume(Sphere(r)),
    exclusion_distance = separation_ratio
);
species = [s1]

tol = 1e-7

# eff_medium = effective_medium(medium, species)

psource = PlaneSource(medium, [0.0,0.0,1.0]);
source = plane_source(medium; direction = [0.0,0.0,1.0])
sourceradial =  regular_spherical_source(medium, [1.0+0.0im];
   position = [0.0,0.0,0.0]
);

region_shape = Sphere([0.0,0.0,0.0], R)
material = Material(Sphere(R),species);

eff_medium = effective_medium(medium, species; numberofparticles = material.numberofparticles)
k_low = ω / eff_medium.c

keffs = wavenumbers(ω, medium, species;
    # num_wavenumbers = 2,
    mesh_points = 10,
    mesh_size = 2.0,
    num_wavenumbers = 3,
    basis_order = basis_order,
    tol = tol,
    symmetry = PlanarAzimuthalSymmetry(),
    numberofparticles = material.numberofparticles
)

pwavemode = WaveMode(ω, k_low, psource, material;
    basis_order = basis_order,
    basis_field_order = basis_field_order
);

pwavemode = WaveMode(ω, keffs[1], psource, material;
       basis_order = basis_order,
       basis_field_order = basis_field_order,
       source_basis_field_order = basis_field_order
   );

wavemode = WaveMode(ω, keffs[1], source, material;
    basis_order = basis_order,
    basis_field_order = basis_field_order,
    source_basis_field_order = basis_field_order
);

scat_field = scattering_field(wavemode)
pscat_field = scattering_field(pwavemode)

# res = discrete_system_residue(pscat_field, ω, source, material, AzimuthalSymmetry{3}();
#     basis_order = basis_order, mesh_points = 5,
#     rtol = 1e-2, maxevals = Int(1e4)
# )

# using Plots
# using Statistics

rtol = 1e-2; maxevals = Int(1e4);

field_order = max(basis_order,Int(round(basis_order/2 + basis_field_order/2)))
legendre_order = field_order + 1;

pair_corr = hole_correction_pair_correlation;

discrete_scat = discrete_system(ω, psource, material;
    basis_order = basis_order,
    basis_field_order = field_order,
    legendre_order = legendre_order,
    rtol = rtol, maxevals = maxevals
);

# xs = [ [cos(θ), sin(θ), 1.1] for θ in 0.0:0.6:(2π), r in LinRange(0.0,R,)]
xs = [ radial_to_cartesian_coordinates([R - 4*outer_radius(s1),θ,0.0]) for θ in 0.0:0.1:(π)];
xs = [ radial_to_cartesian_coordinates([R - 4*outer_radius(s1),0.0,φ]) for φ in 0.0:0.1:(π)];
xs = [ radial_to_cartesian_coordinates([r,0.0,0.0]) for r in 0.0:0.1:(R - 1.1 * outer_radius(s1))];
xs = [ radial_to_cartesian_coordinates([(R - outer_radius(s1)) * rand(),pi * rand(), pi * rand()]) for i in 1:1000];
# xs = [xs; [ radial_to_cartesian_coordinates([2.1*r,θ,0.0]) for θ in 0.0:0.1:(π)]];
# xs = [ radial_to_cartesian_coordinates([r,0.0,0.0]) for θ in 0.0:0.1:(R-r)];

pdata = pscat_field.(xs);
data = scat_field.(xs);

scatdata = discrete_scat.(xs);

maximum(norm.(data - pdata) ./ norm.(pdata))
maximum(norm.(scatdata - pdata) ./ norm.(pdata))
mean(norm.(scatdata - pdata) ./ norm.(pdata))

# pdata0 = [d[1] for d in pdata];
# pdata1 = [d[3] for d in pdata];
#
# pdata11 = [d[4] for d in pdata];
# pdata1m1 = [d[2] for d in pdata];
#
# data0 = [d[1] for d in data];
# data1 = [d[3] for d in data];

mat_coefs_pwave = material_scattering_coefficients(pwavemode);

# inds = findall(abs.(mat_coefs_pwave) .> 1e-6)

# material_scattering_coefficients(wavemode)

mat_coefs_pdisc = material_scattering_coefficients_discrete(pscat_field, ω, source, material;
    basis_order = basis_order,
    basis_field_order = basis_field_order,
    rtol = rtol,
    maxevals = maxevals
);

mat_coefs_disc = material_scattering_coefficients_discrete(discrete_scat, ω, source, material;
    basis_order = basis_order,
    basis_field_order = basis_field_order,
    rtol = rtol,
    maxevals = maxevals
);

L = length(mat_coefs_pdisc)

# should be very small
norm(mat_coefs_pwave[1:L] - mat_coefs_pdisc) / norm(mat_coefs_pwave[1:L])

norm(mat_coefs_disc - mat_coefs_pdisc) / norm(mat_coefs_disc)
norm(mat_coefs_pwave[1:L] - mat_coefs_disc) / norm(mat_coefs_pwave[1:L])

# [mat_coefs_pwave[inds] mat_coefs_disc2[inds] mat_coefs_disc[inds]]

# Calculate low frequency scattering
    Linc = basis_order + basis_field_order
    source_coefficients = regular_spherical_coefficients(source)(Linc,zeros(3),ω);

    material_low = Material(
        Sphere(outer_radius(material.shape) - outer_radius(s1)),
        species
    );

    effective_sphere = Particle(eff_medium, material_low.shape);
    Tmat = MultipleScattering.t_matrix(effective_sphere, medium, ω, Linc);
    scat_coef_low = Tmat * source_coefficients;

# 0.747
scatdata0 = [d[1] for d in scatdata];
scatdata1 = [d[3] for d in scatdata];
scatdata11 = [d[4] for d in scatdata];
scatdata1m1 = [d[2] for d in scatdata];

# plot([real.(scatdata0),real.(pdata0)])
# plot([imag.(scatdata0),imag.(pdata0)])
#
# plot([real.(scatdata1),real.(pdata1)])
# plot([imag.(scatdata1),imag.(pdata1)])

# discrete_scat = discrete_system(ω, psource, material;
#     basis_order = basis_order, basis_field_order = basis_field_order,
#     rtol = 1e-3, atol=1e-3, maxevals=Int(10e4)
# );
# d2 = discrete_scat([1.0,2.0,1.1])
#
# discrete_scat = discrete_system(ω, psource, material;
#     basis_order = basis_order, basis_field_order = basis_field_order,
#     rtol = 1e-4, atol=1e-4, maxevals=Int(20e4)
# );
# d3 = discrete_scat([1.0,2.0,1.1])
# [d1,d2,d3]
#
# plot!(fun.(scatdata0 ))
#
# fun = imag
# fun = real
# plot([fun.(scatdata0),fun.(pdata0)])
#
#
# pdata[1]
# data[1]
#
# # Test chosen basis for field
#
#     R = outer_radius(material.shape)
#     legendre_order = basis_field_order
#     mesh_points = Int(round(1.5*legendre_order)) + 1
#
#     r1s = LinRange(0,R-r, mesh_points)
#     θ1s = LinRange(0,π, mesh_points)
#
#     T = Float64
#     using ClassicalOrthogonalPolynomials
#     function field_basis(rθ::AbstractVector{T},legendre_order::Int)
#         P = Legendre{T}()
#
#         # [P_0(cos(θ)), …, P_(legendre_order-1)(cos(θ))]
#         Pθs = P[cos(rθ[2]), 1:legendre_order]
#         Prs = P[2 * rθ[1] / (R - r) - one(T), 1:legendre_order]
#
#         # [Pr * Pθ for Pr in Prs, Pθ in Pθs][:]
#         return (Prs * transpose(Pθs))[:]
#     end
#
#     field_basis([1.0,0.2])
#
#     len = basisorder_to_basislength(Acoustic{T,3}, basis_order)
#     len_p = legendre_order^2
#     #
#     # function δφj(rθ1::AbstractVector{T}, legendre_order = 5)
#     #     basis1 = field_basis(rθ1, legendre_order)
#     #     return [b1 for b1 in basis1, n in 1:len][:]
#     # end
#
#     function δφj(rθ1::AbstractVector{T}, legendre_order::Int)
#         basis1 = field_basis(rθ1, legendre_order)
#         return reshape(
#             [
#                 (nd == n) ? b1 : zero(Complex{T})
#                 # (n,nd,b1)
#             # for n in -1:1, nd in 1:3, b1 in 0:2],
#             for n in 1:len, nd in 1:len, b1 in basis1],
#         (len, len * len_p))
#     end
#
#     As = [δφj([r1,θ1], legendre_order) for r1 in r1s, θ1 in θ1s][:];
#     A = vcat(As...)
#
#     rθφ2xyz = radial_to_cartesian_coordinates
#     b = [pscat_field(rθφ2xyz([r,θ,0.0])) for r in r1s, θ in θ1s][:]
#     b = vcat(b...)[:]
#
#     as = A \ b;
#     as = reshape(as,(len,:));
#
#     function scattered_field(xs::Vector{T})
#         rθφ = cartesian_to_radial_coordinates(xs)
#         return as * field_basis(rθφ[1:2],legendre_order)
#     end
#
#     rθs = [ [(R-r) * rand(), π * rand(),0.0] for i = 1:2000 ];
#     xs = rθφ2xyz.(rθs);
#
#     scatdata = scattered_field.(xs);
#     pdata = pscat_field.(xs);
#
#     maximum(norm.(scatdata - pdata) ./ norm.(pdata))
#     mean(norm.(scatdata - pdata) ./ norm.(pdata))
#
# # Test reshape used in the discrete_system
#     # ps = -1:1
#     # ns = 1:3
#     # nds = 1:2
#     # js = -1:2
#     #
#     # bs = [
#     #         [
#     #             (n,j)
#     #             # sum((n,j))
#     #         for n in ns]
#     # for j in js];
#     #
#     # bs = vcat(bs...);
#     #
#     # a = [
#     #         (nd,p)
#     #         # nd + p*1.0im
#     # for nd in nds, p in ps][:]
#     #
#     # function Kmat(j)
#     #     reshape(
#     #         [
#     #             (n,j,nd,p)
#     #             # Float64(n)^j + Float64(nd)^p * 1.0im + sum((n,j,nd,p)) + ((n == nd && p == j) ? 10.0 + 2im : 0.0)
#     #         for n in ns, nd in nds, p in ps],
#     #     (length(ns), length(nds) * length(ps)))
#     # end
#     #
#     # data = [
#     #         Kmat(j)
#     # for j in js];
#     #
#     # K = vcat(data...);
#     # bs = K * a
#     #
#     # a2 = inv(transpose(conj.(K)) * K) * transpose(conj.(K)) * bs;
#     # a3 = K \ bs;
#     # maximum(abs.(a - a2))
#     # maximum(abs.(a2 - a3))
#     #
#     # f(x) = cos(x[1]+x[2])
#     # hcubature(f, SVector(0.0,0.0), SVector(30pi,10pi);
#     #     rtol=rtol, atol=atol, maxevals=100000
#     # )
#
#     function δφ(j)
#         return reshape(
#             [
#                 (n,j,nd,p)
#             for n in ns, nd in nds, p in ps],
#         (length(ns), length(nds) * length(ps)))
#     end
#
#     data = [
#             δφ(j)
#     for j in js];
#
#     Kδφ = vcat(data...);
#
#
#
#     legendre_order = 15
#     legendre_order = 20
#     function field_basis(rθ::AbstractVector{T})
#         P = Legendre{T}()
#         # [P_0(cos(θ)), …, P_(legendre_order-1)(cos(θ))]
#         Pθs = P[cos(rθ[2]), 1:legendre_order]
#         Prs = P[2 * rθ[1] / (R - outer_radius(s1)) - one(T), 1:legendre_order]
#
#         # [Pr * Pθ for Pr in Prs, Pθ in Pθs][:]
#         return (Prs * transpose(Pθs))[:]
#     end
#
#
#     rs = 0.0:0.1:(R - outer_radius(s1))
#     θs = 0.0:0.1:π
#
#     associatedlegendre(m) = ((-1)^m*prod(1:2:(2m-1)))*(UltrasphericalWeight((m+1)/2).*Ultraspherical(m+1/2))
#     f(rθ) = sbesselj(3, 0.1 * rθ[1]) * associatedlegendre(2)[cos(rθ[2]),3] * exp(2im) + sbesselj(1,0.1 * rθ[1]) * associatedlegendre(1)[cos(rθ[2]),2] * exp(1im)
#
#     # dot(as,field_basis(rθ)) = f(rθ)
#
#     A = vcat([transpose(field_basis([r,θ])) for r in rs, θ in θs][:]...)
#     b = [f([r,θ]) for r in rs, θ in θs][:]
#
#     x = A \ b;
#
#     maximum(abs.(A * x - b)) / mean(abs.(b))
#
#     f1 = [ f([rs[10],θ])  for θ in 0.0:0.03:pi]
#     f2 = [ sum(field_basis([rs[10],θ]) .* x) for θ in 0.0:0.03:pi]
#
#     sum(field_basis([0.2,0.3]) .* x)


text:
using KernSmooth, RDatasets, Base.Test

expected = readcsv(joinpath(Pkg.dir("KernSmooth"), "test", "geyserfit.csv"), Float64, has_header=true)[1]

dat = dataset("MASS", "geyser")

x = array(dat[:Duration])
y = convert(Vector{Float64}, array(dat[:Waiting]))

resx, resy= locpoly(x, y, 0.25)

@test_approx_eq_eps(maximum(abs(resx - expected[:,1])), 0.0, sqrt(eps(Float64)))
@test_approx_eq_eps(maximum(abs(resy - expected[:,2])), 0.0, sqrt(eps(Float64)))

@test !any(isnan(resy))

#example from R function
@test_approx_eq dpill(x, y) 0.238350165926485


text:
module Solution
  export rentalcarcost
  rentalcarcost(d) = 40d - (d≥7 ? 50 : d≥3 ? 20 : 0)
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using Statistics
include("dsp.jl")

function sdnn(nn_interval::Vector{<:Real})

    sdnn = std(nn_interval)
    return sdnn
end

function sdann_asdnn(nn_interval::Vector{<:Real}, index_array::Vector{<:Integer},
    duration::Integer, min_number_of_interval::Integer = 3)

    if length(nn_interval) < min_number_of_interval
        return nothing, nothing
    end

    segment_index = get_segment_index(index_array, duration)

    average_array::Vector{Float64} = []
    std_array::Vector{Float64} = []

    start_idx = 1
    for end_idx in segment_index
        landmark = @view nn_interval[start_idx:end_idx]

        if length(landmark) >= min_number_of_interval
            push!(average_array, mean(landmark))
            push!(std_array, std(landmark))
        end
        start_idx = end_idx
    end

    if length(average_array) >= min_number_of_interval

        sdann = std(average_array)
        asdnn = mean(std_array)
        return sdann, asdnn
    end

    return nothing, nothing
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE

using FEMBase, Test, AsterReader
using AsterReader: aster_parse_nodes, aster_read_mesh, MEDFile, get_element_sets

datadir = first(splitext(basename(@__FILE__)))

@testset "parse nodes from ascii file" begin
    section = """
    N9  2.0 3.0 4.0
    COOR_3D
    N1          0.0 0.0 0.0
    N2          1.0 0.0 0.0
    N3          1.0 1.0 0.0
    N4          0.0 1.0 0.0
    N5          0.0 0.0 1.0
    N6          1.0 0.0 1.0
    N7          1.0 1.0 1.0
    N8          0.0 1.0 1.0
    FINSF
    absdflasdf
    N12 3.0 4.0 5.0 6.0
    N13 3.0 4.0 5.0
    """
    nodes = aster_parse_nodes(section)
    @test nodes[1] == Float64[0.0, 0.0, 0.0]
    @test nodes[8] == Float64[0.0, 1.0, 1.0]
    @test length(nodes) == 8
end

@testset "read mesh med file" begin
    meshfile = joinpath(datadir, "quad4.med")
    mesh = aster_read_mesh(meshfile)
    @test length(mesh["element_sets"]) == 5
    @test length(mesh["node_sets"]) == 4
    @test length(mesh["elements"]) == 5
    @test length(mesh["nodes"]) == 4
    for elset in ["BLOCK", "TOP", "BOTTOM", "LEFT", "RIGHT"]
        @test haskey(mesh["element_sets"], elset)
        @test length(mesh["element_sets"][elset]) == 1
    end
    for nset in ["TOP_LEFT", "TOP_RIGHT", "BOTTOM_LEFT", "BOTTOM_RIGHT"]
        @test haskey(mesh["node_sets"], nset)
        @test length(mesh["node_sets"][nset]) == 1
    end
end

@testset "test read element sets from med file, issue #111 at JuliaFEM.jl" begin
    meshfile = joinpath(datadir, "hexmeshOverlappingGroups.med")
    med = MEDFile(meshfile)
    element_sets = get_element_sets(med, "Mesh_1")
    @test element_sets[-10] == ["halfhex", "mosthex"]
    @test element_sets[-11] == ["halfhex"]
end

@testset "test read overlapping ets, issue #111 at JuliaFEM.jl" begin
    mesh_file = joinpath(datadir, "hexmeshOverlappingGroups.med")
    mesh = aster_read_mesh(mesh_file, "Mesh_1")
    @test length(mesh["element_sets"]["mosthex"]) == 273
    @test length(mesh["element_sets"]["halfhex"]) == 147
end


text:
export iauEpv00
"""
Earth position and velocity, heliocentric and barycentric, with
respect to the Barycentric Celestial Reference System.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
      date1,date2  double        TDB date (Note 1)

Returned:
      pvh          double[2][3]  heliocentric Earth position/velocity
      pvb          double[2][3]  barycentric Earth position/velocity

      Returned (function value):
                  int           status: 0 = OK
                                       +1 = warning: date outside
                                                the range 1900-2100 AD

Notes:

      1. The TDB date date1+date2 is a Julian Date, apportioned in any
      convenient way between the two arguments.  For example,
      JD(TDB)=2450123.7 could be expressed in any of these ways, among
      others:

            date1          date2

            2450123.7           0.0       (JD method)
            2451545.0       -1421.3       (J2000 method)
            2400000.5       50123.2       (MJD method)
            2450123.5           0.2       (date & time method)

      The JD method is the most natural and convenient to use in cases
      where the loss of several decimal digits of resolution is
      acceptable.  The J2000 method is best matched to the way the
      argument is handled internally and will deliver the optimum
      resolution.  The MJD method and the date & time methods are both
      good compromises between resolution and convenience.  However,
      the accuracy of the result is more likely to be limited by the
      algorithm itself than the way the date has been expressed.

      n.b. TT can be used instead of TDB in most applications.

      2. On return, the arrays pvh and pvb contain the following:

            pvh[0][0]  x       }
            pvh[0][1]  y       } heliocentric position, au
            pvh[0][2]  z       }

            pvh[1][0]  xdot    }
            pvh[1][1]  ydot    } heliocentric velocity, au/d
            pvh[1][2]  zdot    }

            pvb[0][0]  x       }
            pvb[0][1]  y       } barycentric position, au
            pvb[0][2]  z       }

            pvb[1][0]  xdot    }
            pvb[1][1]  ydot    } barycentric velocity, au/d
            pvb[1][2]  zdot    }

      The vectors are with respect to the Barycentric Celestial
      Reference System.  The time unit is one day in TDB.

      3. The function is a SIMPLIFIED SOLUTION from the planetary theory
      VSOP2000 (X. Moisson, P. Bretagnon, 2001, Celes. Mechanics &
      Dyn. Astron., 80, 3/4, 205-213) and is an adaptation of original
      Fortran code supplied by P. Bretagnon (private comm., 2000).

      4. Comparisons over the time span 1900-2100 with this simplified
      solution and the JPL DE405 ephemeris give the following results:

                                    RMS    max
            Heliocentric:
                  position error    3.7   11.2   km
                  velocity error    1.4    5.0   mm/s

            Barycentric:
                  position error    4.6   13.4   km
                  velocity error    1.4    4.9   mm/s

      Comparisons with the JPL DE406 ephemeris show that by 1800 and
      2200 the position errors are approximately double their 1900-2100
      size.  By 1500 and 2500 the deterioration is a factor of 10 and
      by 1000 and 3000 a factor of 60.  The velocity accuracy falls off
      at about half that rate.

      5. It is permissible to use the same array for pvh and pvb, which
      will receive the barycentric values.

This revision:  2017 March 16

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauEpv00(double date1, double date2,
#              double pvh[2][3], double pvb[2][3])

function iauEpv00(date1::Real, date2::Real)
      # Initialize function return variables
      pvh = zeros(Float64, 3, 2)
      pvb = zeros(Float64, 3, 2)

      status = ccall((:iauEpv00, libsofa_c), Cint,
            (Cdouble, Cdouble,
            Ref{Cdouble}, Ref{Cdouble}),
            convert(Float64, date1), convert(Float64, date2),
            pvh, pvb)

      return status, SMatrix{2,3}(pvh'), SMatrix{2,3}(pvb')
   end


text:
function auto_vecjac!(du, f, x, v, cache1 = nothing, cache2 = nothing)
    !hasmethod(f, (typeof(x),)) &&
        error("For inplace function use autodiff = false")
    du .= reshape(auto_vecjac(f, x, v), size(du))
end

function auto_vecjac(f, x, v)
    vv, back = Zygote.pullback(f, x)
    return vec(back(reshape(v, size(vv)))[1])
end


text:
module MultiFloats

export MultiFloat, renormalize,
       Float64x1, Float64x2, Float64x3, Float64x4,
       Float64x5, Float64x6, Float64x7, Float64x8,
       use_clean_multifloat_arithmetic,
       use_sloppy_multifloat_arithmetic,
       use_very_sloppy_multifloat_arithmetic

include("./MultiFloatsCodeGen.jl")
using .MultiFloatsCodeGen

####################################################### DEFINITION OF MULTIFLOAT

struct MultiFloat{T<:AbstractFloat,N} <: AbstractFloat
    _limbs::NTuple{N,T}
end

const Float16x{N} = MultiFloat{Float16,N}
const Float32x{N} = MultiFloat{Float32,N}
const Float64x{N} = MultiFloat{Float64,N}

const AF = Core.AbstractFloat
const MF = MultiFloat

const Float64x1 = Float64x{1}
const Float64x2 = Float64x{2}
const Float64x3 = Float64x{3}
const Float64x4 = Float64x{4}
const Float64x5 = Float64x{5}
const Float64x6 = Float64x{6}
const Float64x7 = Float64x{7}
const Float64x8 = Float64x{8}

################################################ CONVERSION FROM PRIMITIVE TYPES

@inline MultiFloat{T,N}(x::MultiFloat{T,N}) where {T<:AF,N} = x

@inline MultiFloat{T,N}(x::T) where {T<:AF,N} =
    MultiFloat{T,N}((x, ntuple(_ -> zero(T), N - 1)...))

@inline MultiFloat{T,N}(x::MultiFloat{T,M}) where {T<:AbstractFloat,M,N} =
    MultiFloat{T,N}((
        ntuple(i -> x._limbs[i], min(M, N))...,
        ntuple(_ -> zero(T), max(N - M, 0))...))

# Values of the types Bool, Int8, UInt8, Int16, UInt16, Float16, Int32, UInt32,
# and Float32 can be converted losslessly to a single Float64, which has 53
# bits of integer precision.

@inline Float64x{N}(x::Bool   ) where {N} = Float64x{N}(Float64(x))
@inline Float64x{N}(x::Int8   ) where {N} = Float64x{N}(Float64(x))
@inline Float64x{N}(x::UInt8  ) where {N} = Float64x{N}(Float64(x))
@inline Float64x{N}(x::Int16  ) where {N} = Float64x{N}(Float64(x))
@inline Float64x{N}(x::UInt16 ) where {N} = Float64x{N}(Float64(x))
@inline Float64x{N}(x::Float16) where {N} = Float64x{N}(Float64(x))
@inline Float64x{N}(x::Int32  ) where {N} = Float64x{N}(Float64(x))
@inline Float64x{N}(x::UInt32 ) where {N} = Float64x{N}(Float64(x))
@inline Float64x{N}(x::Float32) where {N} = Float64x{N}(Float64(x))

# Values of the types Int64, UInt64, Int128, and UInt128 cannot be converted
# losslessly to a single Float64 and must be split into multiple components.

@inline Float64x1(x::Int64  ) = Float64x1(Float64(x))
@inline Float64x1(x::UInt64 ) = Float64x1(Float64(x))
@inline Float64x1(x::Int128 ) = Float64x1(Float64(x))
@inline Float64x1(x::UInt128) = Float64x1(Float64(x))

@inline function Float64x2(x::Int128)
    x0 = Float64(x)
    x1 = Float64(x - Int128(x0))
    Float64x2((x0, x1))
end

@inline function Float64x2(x::UInt128)
    x0 = Float64(x)
    x1 = Float64(reinterpret(Int128, x - UInt128(x0)))
    Float64x2((x0, x1))
end

@inline function Float64x{N}(x::Int64) where {N}
    x0 = Float64(x)
    x1 = Float64(x - Int64(x0))
    Float64x{N}((x0, x1, ntuple(_ -> 0.0, N - 2)...))
end

@inline function Float64x{N}(x::UInt64) where {N}
    x0 = Float64(x)
    x1 = Float64(reinterpret(Int64, x - UInt64(x0)))
    Float64x{N}((x0, x1, ntuple(_ -> 0.0, N - 2)...))
end

@inline function Float64x{N}(x::Int128) where {N}
    x0 = Float64(x)
    r1 = x - Int128(x0)
    x1 = Float64(r1)
    x2 = Float64(r1 - Int128(x1))
    Float64x{N}((x0, x1, x2, ntuple(_ -> 0.0, N - 3)...))
end

@inline function Float64x{N}(x::UInt128) where {N}
    x0 = Float64(x)
    r1 = reinterpret(Int128, x - UInt128(x0))
    x1 = Float64(r1)
    x2 = Float64(r1 - Int128(x1))
    Float64x{N}((x0, x1, x2, ntuple(_ -> 0.0, N - 3)...))
end

################################################## CONVERSION TO PRIMITIVE TYPES

@inline Base.Float16(x::Float64x{N}) where {N} = Float16(x._limbs[1])
@inline Base.Float32(x::Float64x{N}) where {N} = Float32(x._limbs[1])

@inline Base.Float16(x::Float16x{N}) where {N} = x._limbs[1]
@inline Base.Float32(x::Float32x{N}) where {N} = x._limbs[1]
@inline Base.Float64(x::Float64x{N}) where {N} = x._limbs[1]

###################################################### CONVERSION FROM BIG TYPES

function MultiFloat{T,N}(x::BigFloat) where {T<:AF,N}
    setprecision(Int(precision(x))) do
        r = Vector{BigFloat}(undef, N)
        y = Vector{T}(undef, N)
        r[1] = x
        y[1] = T(r[1])
        for i = 2 : N
            r[i] = r[i-1] - y[i-1]
            y[i] = T(r[i])
        end
        MultiFloat{T,N}((y...,))
    end
end

function MultiFloat{T,N}(x::BigInt) where {T<:AF,N}
    y = Vector{T}(undef, N)
    for i = 1 : N
        y[i] = T(x)
        x -= BigInt(y[i])
    end
    MultiFloat{T,N}((y...,))
end

MultiFloat{T,N}(x::Rational{U}) where {T<:AF,N,U} =
    MF{T,N}(numerator(x)) / MF{T,N}(denominator(x))

######################################################## CONVERSION TO BIG TYPES

Base.BigFloat(x::MultiFloat{T,N}) where {T<:AF,N} =
    +(ntuple(i -> BigFloat(x._limbs[N-i+1]), N)...)

################################################################ PROMOTION RULES

Base.promote_rule(::Type{MF{T,N}}, ::Type{T      }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{Int8   }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{Int16  }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{Int32  }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{Int64  }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{Int128 }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{Bool   }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{UInt8  }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{UInt16 }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{UInt32 }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{UInt64 }) where {T<:AF,N} = MF{T,N}
Base.promote_rule(::Type{MF{T,N}}, ::Type{UInt128}) where {T<:AF,N} = MF{T,N}

Base.promote_rule(::Type{MF{T,N}}, ::Type{BigFloat}) where {T<:AF,N} = BigFloat

Base.promote_rule(::Type{Float32x{N}}, ::Type{Float16}) where {N} = Float32x{N}
Base.promote_rule(::Type{Float64x{N}}, ::Type{Float16}) where {N} = Float64x{N}
Base.promote_rule(::Type{Float64x{N}}, ::Type{Float32}) where {N} = Float64x{N}

####################################################################### PRINTING

@inline function renormalize(x::MF{T,N}) where {T<:AF,N}
    total = +(x._limbs...)
    if isfinite(total)
        while true
            x0::MF{T,N} = x + zero(T)
            if !(x0._limbs != x._limbs); break; end
            x = x0
        end
        x
    else
        MultiFloat{T,N}(ntuple(_ -> total, N))
    end
end

@inline renormalize(x::T) where {T<:Number} = x

function call_normalized(callback, x::MultiFloat{T,N}) where {T<:AF,N}
    x = renormalize(x)
    if !isfinite(x._limbs[1])
        callback(x._limbs[1])
    else
        i = N
        while (i > 0) && iszero(x._limbs[i])
            i -= 1
        end
        if iszero(i)
            callback(zero(T))
        else
            setprecision(() -> callback(BigFloat(x)),
                precision(T) + exponent(x._limbs[1]) - exponent(x._limbs[i]))
        end
    end
end

function Base.show(io::IO, x::MultiFloat{T,N}) where {T<:AF,N}
    call_normalized(y -> show(io, y), x)
end

################################################################################

# Thanks to Greg Plowman (https://github.com/GregPlowman) for suggesting
# implementations of Printf.fix_dec and Printf.ini_dec for @printf support.

import Printf: fix_dec, ini_dec

if VERSION < v"1.1"

    fix_dec(out, x::MultiFloat{T,N}, flags::String, width::Int, precision::Int, c::Char) where {T<:AF,N} =
        call_normalized(d -> fix_dec(out, BigFloat(d), flags, width, precision, c), x)

    ini_dec(out, x::MultiFloat{T,N}, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) where {T<:AF,N} =
        call_normalized(d -> ini_dec(out, BigFloat(d), ndigits, flags, width, precision, c), x)

else

    fix_dec(out, x::MultiFloat{T,N}, flags::String, width::Int, precision::Int, c::Char, digits) where {T<:AF,N} =
        call_normalized(d -> fix_dec(out, BigFloat(d), flags, width, precision, c, digits), x)

    ini_dec(out, x::MultiFloat{T,N}, ndigits::Int, flags::String, width::Int, precision::Int, c::Char, digits) where {T<:AF,N} =
        call_normalized(d -> ini_dec(out, BigFloat(d), ndigits, flags, width, precision, c, digits), x)

end

################################################################################

@inline Base.:(==)(x::MF{T,N}, y::MF{T,N}) where {T<:AF,N} = _eq(renormalize(x), renormalize(y))
@inline Base.:(!=)(x::MF{T,N}, y::MF{T,N}) where {T<:AF,N} = _ne(renormalize(x), renormalize(y))
@inline Base.:(< )(x::MF{T,N}, y::MF{T,N}) where {T<:AF,N} = _lt(renormalize(x), renormalize(y))
@inline Base.:(> )(x::MF{T,N}, y::MF{T,N}) where {T<:AF,N} = _gt(renormalize(x), renormalize(y))
@inline Base.:(<=)(x::MF{T,N}, y::MF{T,N}) where {T<:AF,N} = _le(renormalize(x), renormalize(y))
@inline Base.:(>=)(x::MF{T,N}, y::MF{T,N}) where {T<:AF,N} = _ge(renormalize(x), renormalize(y))

@inline scale(a::T, x::MultiFloat{T,N}) where {T<:AF,N} =
    MultiFloat{T,N}(ntuple(i -> a * x._limbs[i], N))

@inline function Base.ldexp(x::MF{T,N}, n::U) where {T<:AF,N,U<:Integer}
    x = renormalize(x)
    MultiFloat{T,N}(ntuple(i -> ldexp(x._limbs[i], n), N))
end

################################################################################

@inline Base.zero(::Type{MF{T,N}}) where {T<:AF,N} = MF{T,N}(zero(T)  )
@inline Base.one( ::Type{MF{T,N}}) where {T<:AF,N} = MF{T,N}(one( T)  )
@inline Base.eps( ::Type{MF{T,N}}) where {T<:AF,N} = MF{T,N}(eps( T)^N)

@inline _iszero(x::MF{T,N}) where {T<:AF,N} =
    (&)(ntuple(i -> iszero(x._limbs[i]), N)...)
@inline _isone( x::MF{T,N}) where {T<:AF,N} =
    isone(x._limbs[1]) & (&)(ntuple(i -> iszero(x._limbs[i + 1]), N - 1)...)

@inline Base.iszero(x::MF{T,1}) where {T<:AF  } =  iszero(x._limbs[1])
@inline Base.isone( x::MF{T,1}) where {T<:AF  } =  isone( x._limbs[1])
@inline Base.iszero(x::MF{T,N}) where {T<:AF,N} = _iszero(renormalize(x))
@inline Base.isone( x::MF{T,N}) where {T<:AF,N} = _isone( renormalize(x))

################################################################################

@inline Base.precision(::Type{MF{T,N}}) where {T<:AF,N} = N * precision(T)

@inline Base.floatmin(::Type{MF{T,N}}) where {T<:AF,N} = MF{T,N}(floatmin(T))
@inline Base.floatmax(::Type{MF{T,N}}) where {T<:AF,N} = MF{T,N}(floatmax(T))

@inline Base.typemin(::Type{MF{T,N}}) where {T<:AF,N} = MF{T,N}(ntuple(_ -> typemin(T), N))
@inline Base.typemax(::Type{MF{T,N}}) where {T<:AF,N} = MF{T,N}(ntuple(_ -> typemax(T), N))

@inline Base.exponent(   x::MF{T,N}) where {T<:AF,N} = exponent(   renormalize(x)._limbs[1])
@inline Base.signbit(    x::MF{T,N}) where {T<:AF,N} = signbit(    renormalize(x)._limbs[1])
@inline Base.issubnormal(x::MF{T,N}) where {T<:AF,N} = issubnormal(renormalize(x)._limbs[1])
@inline Base.isfinite(   x::MF{T,N}) where {T<:AF,N} = isfinite(   renormalize(x)._limbs[1])
@inline Base.isinf(      x::MF{T,N}) where {T<:AF,N} = isinf(      renormalize(x)._limbs[1])
@inline Base.isnan(      x::MF{T,N}) where {T<:AF,N} = isnan(      renormalize(x)._limbs[1])

@inline function Base.nextfloat(x::MF{T,N}) where {T<:AF,N}
    y = renormalize(x)
    MF{T,N}((ntuple(i -> y._limbs[i], N - 1)..., nextfloat(y._limbs[N])))
end

@inline function Base.prevfloat(x::MF{T,N}) where {T<:AF,N}
    y = renormalize(x)
    MF{T,N}((ntuple(i -> y._limbs[i], N - 1)..., prevfloat(y._limbs[N])))
end

import LinearAlgebra: floatmin2
@inline floatmin2(::Type{MF{T,N}}) where {T<:AF,N} =
    MF{T,N}(ldexp(one(T), div(exponent(floatmin(T)) - N * exponent(eps(T)), 2)))

################################################################################

@inline Base.inv(x::MF{T,N}) where {T<:AF,N} = one(MF{T,N}) / x

@inline function Base.abs(x::MF{T,N}) where {T<:AF,N}
    x = renormalize(x)
    ifelse(signbit(x._limbs[1]), -x, x)
end

@inline function Base.abs2(x::MF{T,N}) where {T<:AF,N}
    x = renormalize(x)
    renormalize(x * x)
end

################################################################################

# To-do list of transcendental math functions to be implemented

Base.exp(  ::MF{T,N}) where {T<:AF,N} = error(  "exp(MultiFloat) not yet implemented")
Base.expm1(::MF{T,N}) where {T<:AF,N} = error("expm1(MultiFloat) not yet implemented")
Base.log(  ::MF{T,N}) where {T<:AF,N} = error(  "log(MultiFloat) not yet implemented")
Base.log2( ::MF{T,N}) where {T<:AF,N} = error( "log2(MultiFloat) not yet implemented")
Base.log10(::MF{T,N}) where {T<:AF,N} = error("log10(MultiFloat) not yet implemented")
Base.log1p(::MF{T,N}) where {T<:AF,N} = error("log1p(MultiFloat) not yet implemented")

Base.sin(::MF{T,N}) where {T<:AF,N} = error("sin(MultiFloat) not yet implemented")
Base.cos(::MF{T,N}) where {T<:AF,N} = error("cos(MultiFloat) not yet implemented")
Base.tan(::MF{T,N}) where {T<:AF,N} = error("tan(MultiFloat) not yet implemented")
Base.sec(::MF{T,N}) where {T<:AF,N} = error("sec(MultiFloat) not yet implemented")
Base.csc(::MF{T,N}) where {T<:AF,N} = error("csc(MultiFloat) not yet implemented")
Base.cot(::MF{T,N}) where {T<:AF,N} = error("cot(MultiFloat) not yet implemented")

Base.sinh(::MF{T,N}) where {T<:AF,N} = error("sinh(MultiFloat) not yet implemented")
Base.cosh(::MF{T,N}) where {T<:AF,N} = error("cosh(MultiFloat) not yet implemented")
Base.tanh(::MF{T,N}) where {T<:AF,N} = error("tanh(MultiFloat) not yet implemented")
Base.sech(::MF{T,N}) where {T<:AF,N} = error("sech(MultiFloat) not yet implemented")
Base.csch(::MF{T,N}) where {T<:AF,N} = error("csch(MultiFloat) not yet implemented")
Base.coth(::MF{T,N}) where {T<:AF,N} = error("coth(MultiFloat) not yet implemented")

Base.sind(::MF{T,N}) where {T<:AF,N} = error("sind(MultiFloat) not yet implemented")
Base.cosd(::MF{T,N}) where {T<:AF,N} = error("cosd(MultiFloat) not yet implemented")
Base.tand(::MF{T,N}) where {T<:AF,N} = error("tand(MultiFloat) not yet implemented")
Base.secd(::MF{T,N}) where {T<:AF,N} = error("secd(MultiFloat) not yet implemented")
Base.cscd(::MF{T,N}) where {T<:AF,N} = error("cscd(MultiFloat) not yet implemented")
Base.cotd(::MF{T,N}) where {T<:AF,N} = error("cotd(MultiFloat) not yet implemented")

Base.asin(::MF{T,N}) where {T<:AF,N} = error("asin(MultiFloat) not yet implemented")
Base.acos(::MF{T,N}) where {T<:AF,N} = error("acos(MultiFloat) not yet implemented")
Base.atan(::MF{T,N}) where {T<:AF,N} = error("atan(MultiFloat) not yet implemented")
Base.asec(::MF{T,N}) where {T<:AF,N} = error("asec(MultiFloat) not yet implemented")
Base.acsc(::MF{T,N}) where {T<:AF,N} = error("acsc(MultiFloat) not yet implemented")
Base.acot(::MF{T,N}) where {T<:AF,N} = error("acot(MultiFloat) not yet implemented")

Base.asinh(::MF{T,N}) where {T<:AF,N} = error("asinh(MultiFloat) not yet implemented")
Base.acosh(::MF{T,N}) where {T<:AF,N} = error("acosh(MultiFloat) not yet implemented")
Base.atanh(::MF{T,N}) where {T<:AF,N} = error("atanh(MultiFloat) not yet implemented")
Base.asech(::MF{T,N}) where {T<:AF,N} = error("asech(MultiFloat) not yet implemented")
Base.acsch(::MF{T,N}) where {T<:AF,N} = error("acsch(MultiFloat) not yet implemented")
Base.acoth(::MF{T,N}) where {T<:AF,N} = error("acoth(MultiFloat) not yet implemented")

Base.asind(::MF{T,N}) where {T<:AF,N} = error("asind(MultiFloat) not yet implemented")
Base.acosd(::MF{T,N}) where {T<:AF,N} = error("acosd(MultiFloat) not yet implemented")
Base.atand(::MF{T,N}) where {T<:AF,N} = error("atand(MultiFloat) not yet implemented")
Base.asecd(::MF{T,N}) where {T<:AF,N} = error("asecd(MultiFloat) not yet implemented")
Base.acscd(::MF{T,N}) where {T<:AF,N} = error("acscd(MultiFloat) not yet implemented")
Base.acotd(::MF{T,N}) where {T<:AF,N} = error("acotd(MultiFloat) not yet implemented")

################################################################################

@inline function two_sum(a::T, b::T) where {T<:AF}
    s = a + b
    v = s - a
    s, (a - (s - v)) + (b - v)
end

@inline function quick_two_sum(a::T, b::T) where {T<:AF}
    s = a + b
    s, b - (s - a)
end

@inline function two_prod(a::T, b::T) where {T<:AF}
    p = a * b
    p, fma(a, b, -p)
end

################################################################################

@inline unsafe_sqrt(x::Float32) = Base.sqrt_llvm(x)
@inline unsafe_sqrt(x::Float64) = Base.sqrt_llvm(x)
@inline unsafe_sqrt(x::T) where {T <: Real} = sqrt(x)

@inline Base.:(==)(x::MF{T,1}, y::MF{T,1}) where {T<:AF} = (x._limbs[1] == y._limbs[1])
@inline Base.:(!=)(x::MF{T,1}, y::MF{T,1}) where {T<:AF} = (x._limbs[1] != y._limbs[1])
@inline Base.:(< )(x::MF{T,1}, y::MF{T,1}) where {T<:AF} = (x._limbs[1] <  y._limbs[1])
@inline Base.:(> )(x::MF{T,1}, y::MF{T,1}) where {T<:AF} = (x._limbs[1] >  y._limbs[1])
@inline Base.:(<=)(x::MF{T,1}, y::MF{T,1}) where {T<:AF} = (x._limbs[1] <= y._limbs[1])
@inline Base.:(>=)(x::MF{T,1}, y::MF{T,1}) where {T<:AF} = (x._limbs[1] >= y._limbs[1])

@inline Base.:+(a::MF{T,1}, b::MF{T,1}) where {T<:AF} = MF{T,1}(a._limbs[1] + b._limbs[1])
@inline Base.:+(a::MF{T,1}, b::T      ) where {T<:AF} = MF{T,1}(a._limbs[1] + b          )
@inline Base.:*(a::MF{T,1}, b::MF{T,1}) where {T<:AF} = MF{T,1}(a._limbs[1] * b._limbs[1])
@inline Base.:*(a::MF{T,1}, b::T      ) where {T<:AF} = MF{T,1}(a._limbs[1] * b          )
@inline Base.:/(a::MF{T,1}, b::MF{T,1}) where {T<:AF} = MF{T,1}(a._limbs[1] / b._limbs[1])
@inline _sqrt(  x::MF{T,1}            ) where {T<:AF} = MF{T,1}(unsafe_sqrt(x._limbs[1]) )

@inline function Base.sqrt(x::MF{T,N}) where {T<:AF,N}
    x = renormalize(x)
    if iszero(x)
        return x
    else
        return _sqrt(x)
    end
end

function use_clean_multifloat_arithmetic(n::Integer=8)
    for i = 2 : n
        eval(multifloat_eq_func(       i              ))
        eval(multifloat_ne_func(       i              ))
        eval(multifloat_lt_func(       i              ))
        eval(multifloat_gt_func(       i              ))
        eval(multifloat_le_func(       i              ))
        eval(multifloat_ge_func(       i              ))
        eval(two_pass_renorm_func(     i, sloppy=false))
        eval(multifloat_add_func(      i, sloppy=false))
        eval(multifloat_float_add_func(i, sloppy=false))
        eval(multifloat_mul_func(      i, sloppy=false))
        eval(multifloat_float_mul_func(i, sloppy=false))
        eval(multifloat_div_func(      i, sloppy=false))
        eval(multifloat_sqrt_func(     i, sloppy=false))
    end
    for (_, v) in MultiFloatsCodeGen.MPADD_CACHE
        eval(v)
    end
end

function use_sloppy_multifloat_arithmetic(n::Integer=8)
    for i = 2 : n
        eval(multifloat_eq_func(       i             ))
        eval(multifloat_ne_func(       i             ))
        eval(multifloat_lt_func(       i             ))
        eval(multifloat_gt_func(       i             ))
        eval(multifloat_le_func(       i             ))
        eval(multifloat_ge_func(       i             ))
        eval(two_pass_renorm_func(     i, sloppy=true))
        eval(multifloat_add_func(      i, sloppy=true))
        eval(multifloat_float_add_func(i, sloppy=true))
        eval(multifloat_mul_func(      i, sloppy=true))
        eval(multifloat_float_mul_func(i, sloppy=true))
        eval(multifloat_div_func(      i, sloppy=true))
        eval(multifloat_sqrt_func(     i, sloppy=true))
    end
    for (_, v) in MultiFloatsCodeGen.MPADD_CACHE
        eval(v)
    end
end

function use_very_sloppy_multifloat_arithmetic(n::Integer=8)
    for i = 2 : n
        eval(multifloat_eq_func(       i             ))
        eval(multifloat_ne_func(       i             ))
        eval(multifloat_lt_func(       i             ))
        eval(multifloat_gt_func(       i             ))
        eval(multifloat_le_func(       i             ))
        eval(multifloat_ge_func(       i             ))
        eval(one_pass_renorm_func(     i, sloppy=true))
        eval(multifloat_add_func(      i, sloppy=true))
        eval(multifloat_float_add_func(i, sloppy=true))
        eval(multifloat_mul_func(      i, sloppy=true))
        eval(multifloat_float_mul_func(i, sloppy=true))
        eval(multifloat_div_func(      i, sloppy=true))
        eval(multifloat_sqrt_func(     i, sloppy=true))
    end
    for (_, v) in MultiFloatsCodeGen.MPADD_CACHE
        eval(v)
    end
end

################################################################################

@inline Base.:+(x::T, y::MF{T,N}) where {T<:AF,N} = y + x
@inline Base.:-(x::T, y::MF{T,N}) where {T<:AF,N} = -(y + (-x))
@inline Base.:*(x::T, y::MF{T,N}) where {T<:AF,N} = y * x

@inline Base.:-(x::MF{T,N}) where {T<:AF,N} = MF{T,N}(ntuple(i -> -x._limbs[i], N))
@inline Base.:-(x::MF{T,N}, y::MF{T,N}) where {T<:AF,N} = x + (-y)
@inline Base.:-(x::MF{T,N}, y::T      ) where {T<:AF,N} = x + (-y)

@inline Base.hypot(x::MF{T,N}, y::MF{T,N}) where {T<:AF,N} = sqrt(x*x + y*y)

################################################################################

use_sloppy_multifloat_arithmetic()

end # module MultiFloats