text:
"""
Spherical cap, transient vibration.

The structure is loaded with a suddenly applied pressure.
The downward deflection at the centre is monitored.
One quarter of the entire shell is modeled.
"""
module spherical_cap_expl_examples

using LinearAlgebra
using SparseArrays
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.AssemblyModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using SymRCM
using VisualStructures: default_layout_3d, plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using PlotlyJS
# using Gnuplot; @gp "clear"
using PGFPlotsX
using CSV
using FinEtools.MeshExportModule.VTKWrite: vtkwritecollection, vtkwrite
using ThreadedSparseCSR
using UnicodePlots
using InteractiveUtils
using BenchmarkTools

const E = 10.5e6*phun("psi")
const nu = 0.3;
const rho = 2.45e-4*phun("lbf*s^2/in^4")
const R = 22.27*phun("in")
const h = 0.41*phun("in")
const Rxy = (R)*sin(26.67/180*pi)
const q = 600*phun("psi")
const ksi = 0.0
const omegad = 2*pi*1000
const color = "red"
const tend = 0.001 #6e-4
const drilling_stiffness_scale = 1.0
const visualize = true

function parloop_csr!(M, K, ksi, U0, V0, tend, dt, force!, peek, nthr)
    U = deepcopy(U0)
    V = deepcopy(U0)
    A = deepcopy(U0)
    F = deepcopy(U0)
    E = deepcopy(U0)
    C = deepcopy(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invMC = deepcopy(U0)
    invMC .= 1.0 ./ (vec(diag(M)) .+ (dt/2) .* C)
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end

    nth = (nthr == 0 ? Base.Threads.nthreads() : nthr)
    @info "$nth threads used"

    t = 0.0
    # Initial Conditions
    @. U = U0; @. V = V0
    A .= invMC .* force!(F, t);
    peek(0, U, V, t)
    @time for step in 1:nsteps
        # Displacement update
        @. U += dt*V + ((dt^2)/2)*A;
        # External loading
        force!(F, t);
        # Add elastic restoring forces
        F .-= ThreadedSparseCSR.bmul!(E, K, U)
        @inbounds @simd for i in eachindex(U)
            _Fi = F[i]; _Ai = A[i]; _Vi = V[i]
            # Add damping forces
            _Fi -= C[i] * (_Vi + (dt/2) * _Ai)
            # Compute the new acceleration.
            _A1i = invMC[i] * _Fi
            A[i] = _A1i
            # Update the velocity
            V[i] += (dt/2)* (_Ai + _A1i);
        end
        t = t + dt
        peek(step, U, V, t)
    end
end

function _execute_parallel_csr(n = 64, nthr = 0)
    color = "red"
    tolerance = min(R, h) / n  / 100

    tolerance = R/n/1000
    fens, fes = Q4circlen(Rxy, n)
    fens, fes = Q4toT3(fens, fes)

    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        x=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (x, y, sqrt(R^2-x^2-y^2))
    end

    let
        fens1 = FENodeSet(fens.xyz ./ phun("in"))
        vtkwrite("spherical_cap_expl_$n.vtu", fens1, fes)
    end

    # bfes = meshboundary(fes)
    # candidates = connectednodes(bfes)
    # fens, fes = mergenodes(fens, fes, tolerance, candidates)
    bfes = meshboundary(fes)
    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    # Renumber the nodes
    femm = FEMMBase(IntegDomain(fes, TriRule(1)))
    C = connectionmatrix(femm, count(fens))
    perm = symrcm(C)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), h), mater)
    # Set up
    femm.drilling_stiffness_scale = drilling_stiffness_scale

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    lx = selectnode(fens; box = Float64[0 0 -Inf Rxy-10*tolerance -Inf Inf], inflate = tolerance)
    ly = selectnode(fens; box = Float64[-Inf Rxy-10*tolerance 0 0 -Inf Inf], inflate = tolerance)
    lo = setdiff(connectednodes(meshboundary(fes)), vcat(lx, ly))
    for i in [2, 4, 6]
        setebc!(dchi, ly, true, i)
    end
    for i in [1, 5, 6]
        setebc!(dchi, lx, true, i)
    end
    for i in 1:6
        setebc!(dchi, lo, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi, perm);
    # numberdofs!(dchi);

    # Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    SM = FinEtoolsFlexStructures.AssemblyModule
    K = FEMMShellT3FFModule.stiffness(femm, SM.SysmatAssemblerSparseCSRSymm(0.0), geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    # Solve
    function pwr(K, M)
        invM = fill(0.0, size(M, 1))
        invM .= 1.0 ./ (vec(diag(M)))
        v = rand(size(M, 1))
        w = fill(0.0, size(M, 1))
        for i in 1:30
            ThreadedSparseCSR.bmul!(w, K, v)
            wn = norm(w)
            w .*= (1.0/wn)
            v .= invM .* w
            vn = norm(v)
            v .*= (1.0/vn)
        end
        sqrt((v' * (K * v)) / (v' * M * v))
    end
    @time omega_max = pwr(K, M)
    @show omega_max

    # @time evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=1, which=:LM, explicittransform=:none)
    # @show sqrt(evals)
    # @show omega_max = sqrt(evals[1])

    @show dt = Float64(0.9* 2/omega_max) * (sqrt(1+ksi^2) - ksi)

    U0 = gathersysvec(dchi)
    V0 = deepcopy(U0)

    cpoint = selectnode(fens; nearestto=[0 0 R])[1]

    cpointdof = dchi.dofnums[cpoint, 3]

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        r = vec(XYZ);
        r .= vec(r)/norm(vec(r))
        forceout[1:3] .= -r*q
        forceout[4:6] .= 0.0
        return forceout
    end

    # Distributed loading on the surface of the shell
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt, 6, computetrac!);
    Fmag = distribloads(lfemm, geom0, dchi, fi, 2);

    function force!(F, t)
        F .= Fmag
        return F
    end

    nsteps = Int(round(tend/dt))
    cdeflections = fill(0.0, nsteps+1)
    kinetic_energy = fill(0.0, nsteps+1)
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, V, t) = begin
        cdeflections[step+1] = U[cpointdof]
        kinetic_energy[step+1] = 1/2 * V' * M * V
        if rem(step+1, nbtw) == 0
            push!(displacements, deepcopy(U))
        end
        nothing
    end

    @info "$nsteps steps"
    parloop_csr!(M, K, ksi, U0, V0, nsteps*dt, dt, force!, peek, nthr)

    savecsv("spherical_cap_expl_$n-KE-$(drilling_stiffness_scale)", t = collect(0.0:dt:(nsteps*dt))./phun("milli*s"), v = kinetic_energy/phun("lbf*in"))

    if visualize
        reference_Bathe = [   0.5310   -0.0008
         10.6195    0.0069
         16.9912    0.0167
         23.8938    0.0302
         30.2655    0.0445
         32.3894    0.0503
         35.5752    0.0513
         41.9469    0.0471
         46.7257    0.0434
         51.5044    0.0471
         55.2212    0.0548
         60.0000    0.0667
         63.1858    0.0754
         67.4336    0.0855
         71.1504    0.0902
         78.5841    0.0855
         82.8319    0.0780
         85.4867    0.0630
         89.2035    0.0460
         94.5133    0.0241
        101.4159   -0.0066
        107.2566   -0.0293
        109.9115   -0.0439
        112.5664   -0.0497
        112.5664   -0.0500
        115.2212   -0.0515
        120.0000   -0.0470
        1.20000e+02  -3.98675e-02
        1.24194e+02  -1.69095e-02
        1.28387e+02   4.63576e-03
        1.34194e+02   3.43046e-02
        1.38065e+02   5.44371e-02
        1.42903e+02   6.92715e-02
        1.46452e+02   7.98675e-02
        1.50323e+02   8.76380e-02
        1.60323e+02   9.47020e-02
        1.70323e+02   9.82340e-02
        1.75806e+02   9.36424e-02
        1.78387e+02   8.62252e-02
        1.79677e+02   7.49227e-02
        1.83548e+02   5.93819e-02
        1.88387e+02   4.94923e-02
        1.93226e+02   3.25386e-02
        1.99355e+02   8.52097e-03  ]
        reference_Bathe[:, 1] .*= 6e-4/120/phun("Milli*s")
        reference_Bathe[:, 2] .*= -1
        reference_Belytschko = [
        0   0.001290322580645
        0.000036507936508  -0.003225806451613
        0.000065079365079  -0.010967741935484
        0.000095238095238  -0.022580645161290
        0.000123809523810  -0.038387096774194
        0.000155555555556  -0.045806451612903
        0.000193650793651  -0.039354838709677
        0.000217460317460  -0.034838709677419
        0.000255555555556  -0.041935483870968
        0.000287301587302  -0.054193548387097
        0.000315873015873  -0.070000000000000
        0.000355555555556  -0.089677419354839
        0.000401587301587  -0.092903225806452
        0.000439682539683  -0.082258064516129
        0.000466666666667  -0.058064516129032
        0.000482539682540  -0.037096774193548
        0.000498412698413  -0.012903225806452
        0.000526984126984   0.017096774193548
        0.000549206349206   0.034193548387097
        0.000582539682540   0.041612903225806
        0.000598412698413   0.038387096774194]
        reference_Belytschko[:, 1] .*= 1/phun("Milli*s")
        reference_Tabiei = [
        0.000003183023873   0.000649350649351
        0.000049336870027  -0.005844155844156
        0.000084350132626  -0.017207792207792
        0.000084350132626  -0.017532467532468
        0.000114588859416  -0.033766233766234
        0.000144827586207  -0.046103896103896
        0.000186206896552  -0.043181818181818
        0.000206896551724  -0.041883116883117
        0.000233952254642  -0.047077922077922
        0.000267374005305  -0.062987012987013
        0.000297612732095  -0.079220779220779
        0.000332625994695  -0.089610389610390
        0.000367639257294  -0.078246753246753
        0.000401061007958  -0.055844155844156
        0.000434482758621  -0.028896103896104
        0.000458355437666  -0.005194805194805
        0.000479045092838   0.017207792207792
        0.000507692307692   0.037012987012987
        0.000539522546419   0.047402597402597
        0.000572944297082   0.035389610389610
        0.000585676392573   0.025649350649351
        0.000600000000000   0.015909090909091
        ]
        reference_Tabiei[:, 1] .*= 1/phun("Milli*s")

        # reference[:, 2] .*= phun("in")
        # @gp  "reset" :-

        # @gp  :- "set terminal cairolatex standalone pdf size 16cm,10.5cm " :-
        # @gp  :- "set output 'fire_severity1.tex' "  :-
        # @gp "clear"
        # @gp "set key top left box  font \"Times-Roman\"" :-
        # @gp  :- reference_Bathe[:, 1] reference_Bathe[:, 2] " lw 2 lc rgb 'black' with points title 'Bathe' "  :-
        # @gp  :- reference_Belytschko[:, 1] reference_Belytschko[:, 2] " lw 2 lc rgb 'black' with points title 'Belytschko' "  :-
        # @gp  :- reference_Tabiei[:, 1] reference_Tabiei[:, 2] " lw 2 lc rgb 'black' with points title 'Tabiei' "  :-
        # @gp  :- reference_S4_impl[1]/phun("Milli*s")  reference_S4_impl[2] " lw 2 lc rgb 'black' with points title 'S4 impl' "  :-
        # @gp  :- reference_S4_expl[1]/phun("Milli*s")  reference_S4_expl[2] " lw 2 lc rgb 'black' with points title 'S4 expl' "  :-
        # @gp  :- collect(0.0:dt:(nsteps*dt))/phun("Milli*s")  cdeflections/phun("in") " lw 2 lc rgb '$color' with lines title 'Present' "  :-

        # @gp  :- "set tics font \"Times-Roman\"" :-
        # @gp  :- "set xlabel 'Time [ms]' font \"Times-Roman\"" :-
        # @gp  :- "set ylabel 'Deflection [in]' font \"Times-Roman\"" :-
        # @gp  :- "set title 'Spherical cap, deflection at the center' font \"Times-Roman\""
        # # @gp  :- "set output "
        # save(term="cairolatex pdf input color dashed size 5in,3.3in", output="spherical_cap_expl_$n.tex")

        objects = []

        xdata, ydata = collect(0.0:dt:(nsteps*dt))/phun("Milli*s"),  cdeflections/phun("in")
        @pgf p = PGFPlotsX.Plot(
        {
        color = "black",
        line_width  = 1.0,
        },
        Coordinates([v for v in  zip(xdata, ydata)])
        )
        push!(objects, p)
        push!(objects, LegendEntry("Present"))

        xdata, ydata = reference_Bathe[:, 1], reference_Bathe[:, 2]
        @pgf p = PGFPlotsX.Plot(
        {
        only_marks,
        color = "black",
        line_width  = 1.0,
        mark = "diamond",
        },
        Coordinates([v for v in  zip(xdata, ydata)])
        )
        push!(objects, p)
        push!(objects, LegendEntry("Bathe"))

        ft = CSV.File("$(@__DIR__)/S4_impl-t.csv")
        fw = CSV.File("$(@__DIR__)/S4_impl-w.csv")
        xdata, ydata = ft["t"][1:47:end]/phun("Milli*s"), fw["w"][1:47:end]
        @pgf p = PGFPlotsX.Plot(
        {
        only_marks,
        color = "black",
        line_width  = 1.0,
        mark = "square",
        },
        Coordinates([v for v in  zip(xdata, ydata)])
        )
        push!(objects, p)
        push!(objects, LegendEntry("A/impl"))

        ft = CSV.File("$(@__DIR__)/S4_expl-t.csv")
        fw = CSV.File("$(@__DIR__)/S4_expl-w.csv")
        xdata, ydata = ft["t"][1:43:end]/phun("Milli*s"), fw["w"][1:43:end]
        @pgf p = PGFPlotsX.Plot(
        {
        only_marks,
        color = "black",
        line_width  = 1.0,
        mark = "star",
        },
        Coordinates([v for v in  zip(xdata, ydata)])
        )
        push!(objects, p)
        push!(objects, LegendEntry("A/expl"))

        @pgf ax = Axis(
        {
        xlabel = "Time [ms]",
        ylabel = "Apex deflection [in]",
        # ymin = -1.1e-5,
        # ymax = 1.1e-5,
        xmin = 0.0,
        xmax = 1.0,
        xmode = "linear",
        ymode = "linear",
        yminorgrids = "true",
        grid = "both",
        legend_style = {
        at = Coordinate(0.5, 1.05),
        anchor = "south",
        legend_columns = -1
        },
        },
        objects...
        )

        display(ax)
        pgfsave("spherical_cap_expl_$n.pdf", ax)

        # Visualization
        @info "Dumping visualization"
        times = Float64[]
        vectors = []
        for i in 1:length(displacements)
            scattersysvec!(dchi, displacements[i])
            push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
            push!(times, i*dt*nbtw)
        end
        vtkwritecollection("spherical_cap_expl_$n", fens, fes, times; vectors = vectors)
    end
end

function test_parallel_csr(ns = [4*64], nthr = 0)
    @info "Clamped cylinder: parallel CSR"
    for n in ns
        _execute_parallel_csr(n, nthr)
    end
    return true
end

function allrun(ns = [22])
    println("#####################################################")
    println("# test_parallel_csr ")
    test_parallel_csr(ns)
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
# Types

const D4 = NTuple{4, Cdouble}

export Float256
struct Float256 <: AbstractFloat
    d4::D4
end

# Functions

function Base.tryparse(::Type{Float256}, str::String)
    r = tryparse(BigFloat, str)
    r === nothing && return r
    Float256(r)
end

export @f256_str
macro f256_str(str::String)
    :(Float256(BigFloat($str)))
end

# function Base.show(io::IO, x::Float256)
#     len = 80
#     buf = Vector{Cchar}(undef, len)
#     prec = 64
#     ccall((:c_qd_swrite, libqd), Cvoid, (Ref{D4}, Cint, Ptr{Cchar}, Cint),
#           x.d4, prec, buf, len)
#     print(io, unsafe_string(pointer(buf)))
# end

Base.show(io::IO, x::Float256) = show(io, BigFloat(x))

# function Float256(x::Float64)
#     r = Ref{D4}()
#     ccall((:c_qd_copy_d, libqd), Cvoid, (Cdouble, Ref{D4}), x, r)
#     Float256(r[])
# end
#
# function Float256(x::Float128)
#     r = Ref{D4}()
#     ccall((:c_qd_copy_dd, libqd), Cvoid, (Ref{D2}, Ref{D4}), x.d2, r)
#     Float256(r[])
# end

Float256(x::Float64) = Float256((x, 0.0, 0.0, 0.0))

Float256(x::Float128) = Float256((x.d2[1], x.d2[2], 0.0, 0.0))

function Float256(x::Union{Bool, Int8, Int16, Int32, UInt8, UInt16, UInt32,
                           Float16, Float32})
    Float256(Float64(x))
end

function Float256(x::Union{Int64, UInt64})
    hi = x >> 32 << 32
    lo = x - hi
    Float256(Float64(hi)) + Float256(Float64(lo))
end

function Float256(x::Rational)
    Float256(numerator(x)) / Float256(denominator(x))
end

function Float256(x::BigFloat)
    r1 = Float256(Float64(x))
    x -= BigFloat(r1)
    r2 = Float256(Float64(x))
    x -= BigFloat(r2)
    r3 = Float256(Float64(x))
    x -= BigFloat(r3)
    r4 = Float256(Float64(x))
    # x -= BigFloat(r4)
    r1 + r2 + r3 + r4
end

Base.Float64(x::Float256) = x.d4[1]

Float128(x::Float256) = Float128((x.d4[1], x.d4[2]))

Base.BigFloat(x::Float256) =
    big(x.d4[1]) + big(x.d4[2]) + big(x.d4[3]) + big(x.d4[4])

Base.big(::Type{Float256}) = BigFloat

Base.promote_rule(::Type{Float256}, ::Type{Float16}) = Float256
Base.promote_rule(::Type{Float256}, ::Type{Float32}) = Float256
Base.promote_rule(::Type{Float256}, ::Type{Float64}) = Float256
Base.promote_rule(::Type{Float256}, ::Type{Float128}) = Float256
Base.promote_rule(::Type{Float256}, ::Type{BigFloat}) = BigFloat
Base.promote_rule(::Type{Float256}, ::Type{I}) where {I <: Integer} = Float256
Base.promote_rule(::Type{Float256}, ::Type{BigInt}) = BigFloat

# Base.zero(::Type{Float256}) = c256_0
# Base.one(::Type{Float256}) = c256_1
# Base.eps(::Type{Float256}) = c256_eps
Base.zero(::Type{Float256}) = Float256(0.0)
Base.one(::Type{Float256}) = Float256(1.0)
Base.eps(::Type{Float256}) = Float256(2^-208)

# function Float256(::Irrational{:π})
#     r = Ref{D4}()
#     ccall((:c_qd_pi, libqd), Cvoid, (Ref{D4},), r)
#     Float256(r[])
# end
Float256(::Irrational{:ℯ}) = c256_e
Float256(::Irrational{:π}) = c256_pi

Base. +(x::Float256) = x

# function Base. -(x::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_neg, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
#     Float256(r[])
# end

function Base. -(a::Float256)
    Float256((- a.d4[1], - a.d4[2], - a.d4[3], - a.d4[4]))
end

# function Base.abs(x::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_abs, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
#     Float256(r[])
# end
Base.abs(a::Float256) = a.d4[1] < 0.0 ? -a : a

Base.inv(x::Float256) = one(Float64) / x

# function Base. +(x::Float256, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_add, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}), x.d4, y.d4, r)
#     Float256(r[])
# end
#
# function Base. +(x::Float128, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_add_dd_qd, libqd), Cvoid, (Ref{D2}, Ref{D4}, Ref{D4}),
#           x.d2, y.d4, r)
#     Float256(r[])
# end
#
# function Base. +(x::Float64, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_add_d_qd, libqd), Cvoid, (Cdouble, Ref{D4}, Ref{D4}),
#           x, y.d4, r)
#     Float256(r[])
# end
#
# function Base. +(x::Float256, y::Float128)
#     r = Ref{D4}()
#     ccall((:c_qd_add_qd_dd, libqd), Cvoid, (Ref{D4}, Ref{D2}, Ref{D4}),
#           x.d4, y.d2, r)
#     Float256(r[])
# end
#
# function Base. +(x::Float256, y::Float64)
#     r = Ref{D4}()
#     ccall((:c_qd_add_qd_d, libqd), Cvoid, (Ref{D4}, Cdouble, Ref{D4}),
#           x.d4, y, r)
#     Float256(r[])
# end

function Base. +(a::Float256, b::Float256)
    s0, t0 = two_sum(a.d4[1], b.d4[1])
    s1, t1 = two_sum(a.d4[2], b.d4[2])
    s2, t2 = two_sum(a.d4[3], b.d4[3])
    s3, t3 = two_sum(a.d4[4], b.d4[4])

    s1, t0 = two_sum(s1, t0)
    s2, t0, t1 = three_sum(s2, t0, t1)
    s3, t0, t2 = three_sum2(s3, t0, t2)
    t0 = t0 + t1 + t3

    s0, s1, s2, s3, t0 = renorm(s0, s1, s2, s3, t0)
    Float256((s0, s1, s2, s3))
end

function Base. +(a::Float256, b::Float64)
    c0, e = two_sum(a.d4[1], b)
    c1, e = two_sum(a.d4[2], e)
    c2, e = two_sum(a.d4[3], e)
    c3, e = two_sum(a.d4[4], e)

    c0, c1, c2, c3, e = renorm(c0, c1, c2, c3, e)

    Float256((c0, c1, c2, c3))
end

function Base. +(a::Float256, b::Float128)
    s0, t0 = two_sum(a.d4[1], b.d2[1])
    s1, t1 = two_sum(a.d4[2], b.d2[2])

    s1, t0 = two_sum(s1, t0)

    s2 = a.d4[3]
    s2, t0, t1 = three_sum(s2, t0, t1)

    s3, t0 = two_sum(t0, a.d4[4])
    t0 += t1

    s0, s1, s2, s3, t0 = renorm(s0, s1, s2, s3, t0)
    Float256((s0, s1, s2, s3))
end

Base. +(a::Float128, b::Float256) = b + a
Base. +(a::Float64, b::Float256) = b + a

# function Base. -(x::Float256, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_sub, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}), x.d4, y.d4, r)
#     Float256(r[])
# end
#
# function Base. -(x::Float128, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_sub_dd_qd, libqd), Cvoid, (Ref{D2}, Ref{D4}, Ref{D4}),
#           x.d2, y.d4, r)
#     Float256(r[])
# end
#
# function Base. -(x::Float64, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_sub_d_qd, libqd), Cvoid, (Cdouble, Ref{D4}, Ref{D4}),
#           x, y.d4, r)
#     Float256(r[])
# end
#
# function Base. -(x::Float256, y::Float128)
#     r = Ref{D4}()
#     ccall((:c_qd_sub_qd_dd, libqd), Cvoid, (Ref{D4}, Ref{D2}, Ref{D4}),
#           x.d4, y.d2, r)
#     Float256(r[])
# end
#
# function Base. -(x::Float256, y::Float64)
#     r = Ref{D4}()
#     ccall((:c_qd_sub_qd_d, libqd), Cvoid, (Ref{D4}, Cdouble, Ref{D4}),
#           x.d4, y, r)
#     Float256(r[])
# end

Base. -(a::Float256, b::Float256) = a + (- b)
Base. -(a::Float256, b::Float64) = a + (- b)
Base. -(a::Float256, b::Float128) = a + (- b)
Base. -(a::Float128, b::Float256) = a + (- b)
Base. -(a::Float64, b::Float256) = a + (- b)

# function Base. *(x::Float256, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_mul, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}), x.d4, y.d4, r)
#     Float256(r[])
# end
#
# # This function as well as its counterpart below are much less
# # accurate than expected. Their accuracy is that of Float128, not
# # Float256. We therefore disable them.
# # function Base. *(x::Float128, y::Float256)
# #     r = Ref{D4}()
# #     ccall((:c_qd_mul_dd_qd, libqd), Cvoid, (Ref{D2}, Ref{D4}, Ref{D4}),
# #           x.d2, y.d4, r)
# #     Float256(r[])
# # end
#
# function Base. *(x::Float64, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_mul_d_qd, libqd), Cvoid, (Cdouble, Ref{D4}, Ref{D4}),
#           x, y.d4, r)
#     Float256(r[])
# end
#
# # function Base. *(x::Float256, y::Float128)
# #     r = Ref{D4}()
# #     ccall((:c_qd_mul_qd_dd, libqd), Cvoid, (Ref{D4}, Ref{D2}, Ref{D4}),
# #           x.d4, y.d2, r)
# #     Float256(r[])
# # end
#
# function Base. *(x::Float256, y::Float64)
#     r = Ref{D4}()
#     ccall((:c_qd_mul_qd_d, libqd), Cvoid, (Ref{D4}, Cdouble, Ref{D4}),
#           x.d4, y, r)
#     Float256(r[])
# end

function Base. *(a::Float256, b::Float256)
    # sloppy_mul
    # quad-double * quad-double
    # a0 * b0                    0
    #      a0 * b1               1
    #      a1 * b0               2
    #           a0 * b2          3
    #           a1 * b1          4
    #           a2 * b0          5
    #                a0 * b3     6
    #                a1 * b2     7
    #                a2 * b1     8
    #                a3 * b0     9

    p0, q0 = two_prod(a.d4[1], b.d4[1])

    p1, q1 = two_prod(a.d4[1], b.d4[2])
    p2, q2 = two_prod(a.d4[2], b.d4[1])

    p3, q3 = two_prod(a.d4[1], b.d4[3])
    p4, q4 = two_prod(a.d4[2], b.d4[2])
    p5, q5 = two_prod(a.d4[3], b.d4[1])

    # Start Accumulation
    p1, p2, q0 = three_sum(p1, p2, q0)

    # Six-Three Sum  of p2, q1, q2, p3, p4, p5.
    p2, q1, q2 = three_sum(p2, q1, q2)
    p3, p4, p5 = three_sum(p3, p4, p5)
    # compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5).
    s0, t0 = two_sum(p2, p3)
    s1, t1 = two_sum(q1, p4)
    s2 = q2 + p5
    s1, t0 = two_sum(s1, t0)
    s2 += t0 + t1

    # O(eps^3) order terms
    s1 += a.d4[1] * b.d4[4] + a.d4[2] * b.d4[3] + a.d4[3] * b.d4[2] +
        a.d4[4] * b.d4[1] + q0 + q3 + q4 + q5
    p0, p1, s0, s1, s2 = renorm(p0, p1, s0, s1, s2)
    return Float256((p0, p1, s0, s1))
end

# This function is much less accurate than expected. Their accuracy is
# that of Float128, not Float256. We therefore disable it.
# function Base. *(a::Float256, b::Float128)
#     # quad-double * double-double
#     # a0 * b0                        0
#     #      a0 * b1                   1
#     #      a1 * b0                   2
#     #           a1 * b1              3
#     #           a2 * b0              4
#     #                a2 * b1         5
#     #                a3 * b0         6
#     #                     a3 * b1    7
#
#     p0, q0 = two_prod(a.d4[1], b.d2[1])
#     p1, q1 = two_prod(a.d4[1], b.d2[2])
#     p2, q2 = two_prod(a.d4[2], b.d2[1])
#     p3, q3 = two_prod(a.d4[2], b.d2[2])
#     p4, q4 = two_prod(a.d4[3], b.d2[1])
#
#     p1, p2, q0 = three_sum(p1, p2, q0)
#
#     # Five-Three-Sum
#     p2, p3, p4 = three_sum(p2, p3, p4)
#     q1, q2 = two_sum(q1, q2)
#     s0, t0 = two_sum(p2, q1)
#     s1, t1 = two_sum(p3, q2)
#     s1, t0 = two_sum(s1, t0)
#     s2 = t0 + t1 + p4
#     p2 = s0
#
#     p3 = a.d4[3] * b.d2[1] + a.d4[4] * b.d2[2] + q3 + q4
#     p3, q0, s1 = three_sum2(p3, q0, s1)
#     p4 = q0 + s2
#
#     p0, p1, p2, p3, p4 = renorm(p0, p1, p2, p3, p4)
#     Float256((p0, p1, p2, p3))
# end

# This is an adaptation from the *(::Float256, ::Float256) case,
# manually partially evaluated to avoid operations involving zeros.
function Base. *(a::Float256, b::Float128)
    # sloppy_mul
    # quad-double * quad-double
    # a0 * b0                    0
    #      a0 * b1               1
    #      a1 * b0               2
    #           a1 * b1          4
    #           a2 * b0          5
    #                a2 * b1     8
    #                a3 * b0     9

    p0, q0 = two_prod(a.d4[1], b.d2[1])

    p1, q1 = two_prod(a.d4[1], b.d2[2])
    p2, q2 = two_prod(a.d4[2], b.d2[1])

    p4, q4 = two_prod(a.d4[2], b.d2[2])
    p5, q5 = two_prod(a.d4[3], b.d2[1])

    # Start Accumulation
    p1, p2, q0 = three_sum(p1, p2, q0)

    # Six-Three Sum  of p2, q1, q2, p3, p4, p5.
    p2, q1, q2 = three_sum(p2, q1, q2)
    p3, p4 = two_sum(p4, p5)
    # compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5).
    s0, t0 = two_sum(p2, p3)
    s1, t1 = two_sum(q1, p4)
    s2 = q2
    s1, t0 = two_sum(s1, t0)
    s2 += t0 + t1

    # O(eps^3) order terms
    s1 += a.d4[3] * b.d2[2] + a.d4[4] * b.d2[1] + q0 + q4 + q5
    p0, p1, s0, s1, s2 = renorm(p0, p1, s0, s1, s2)
    return Float256((p0, p1, s0, s1))
end

function Base. *(a::Float256, b::Float64)
    p0, q0 = two_prod(a.d4[1], b)
    p1, q1 = two_prod(a.d4[2], b)
    p2, q2 = two_prod(a.d4[3], b)
    p3 = a.d4[4] * b

    s0 = p0

    s1, s2 = two_sum(q0, p1)

    s2, q1, p2 = three_sum(s2, q1, p2)

    q1, q2, p3 = three_sum2(q1, q2, p3)
    s3 = q1

    s4 = q2 + p2

    s0, s1, s2, s3, s4 = renorm(s0, s1, s2, s3, s4)
    Float256((s0, s1, s2, s3))
end

Base. *(a::Float64, b::Float256) = b * a
Base. *(a::Float128, b::Float256) = b * a

function Base. /(x::Float256, y::Float256)
    r = Ref{D4}()
    ccall((:c_qd_div, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}), x.d4, y.d4, r)
    Float256(r[])
end

# function Base. /(x::Float128, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_div_dd_qd, libqd), Cvoid, (Ref{D2}, Ref{D4}, Ref{D4}),
#           x.d2, y.d4, r)
#     Float256(r[])
# end
#
# function Base. /(x::Float64, y::Float256)
#     r = Ref{D4}()
#     ccall((:c_qd_div_d_qd, libqd), Cvoid, (Cdouble, Ref{D4}, Ref{D4}),
#           x, y.d4, r)
#     Float256(r[])
# end
#
# function Base. /(x::Float256, y::Float128)
#     r = Ref{D4}()
#     ccall((:c_qd_div_qd_dd, libqd), Cvoid, (Ref{D4}, Ref{D2}, Ref{D4}),
#           x.d4, y.d2, r)
#     Float256(r[])
# end
#
# function Base. /(x::Float256, y::Float64)
#     r = Ref{D4}()
#     ccall((:c_qd_div_qd_d, libqd), Cvoid, (Ref{D4}, Cdouble, Ref{D4}),
#           x.d4, y, r)
#     Float256(r[])
# end

Base. \(x::Float256, y::Float256) = y / x
# Base. \(x::Float128, y::Float256) = y / x
# Base. \(x::Float64, y::Float256) = y / x
# Base. \(x::Float256, y::Float128) = y / x
# Base. \(x::Float256, y::Float64) = y / x

function Base.cmp(x::Float256, y::Float256)
    r = Ref{Cint}()
    ccall((:c_qd_comp, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{Cint}),
          x.d4, y.d4, r)
    Int(r[])
end

function Base.cmp(x::Float64, y::Float256)
    r = Ref{Cint}()
    ccall((:c_qd_comp_d_qd, libqd), Cvoid, (Cdouble, Ref{D4}, Ref{Cint}),
          x, y.d4, r)
    Int(r[])
end

function Base.cmp(x::Float256, y::Float64)
    r = Ref{Cint}()
    ccall((:c_qd_comp_qd_d, libqd), Cvoid, (Ref{D4}, Cdouble, Ref{Cint}),
          x.d4, y, r)
    Int(r[])
end

# Base. ==(x::Float256, y::Float256) = cmp(x, y) == 0
# Base. ==(x::Float64, y::Float256) = cmp(x, y) == 0
# Base. ==(x::Float256, y::Float64) = cmp(x, y) == 0
#
# Base. !=(x::Float256, y::Float256) = cmp(x, y) != 0
# Base. !=(x::Float64, y::Float256) = cmp(x, y) != 0
# Base. !=(x::Float256, y::Float64) = cmp(x, y) != 0
#
# Base. <(x::Float256, y::Float256) = cmp(x, y) < 0
# Base. <(x::Float64, y::Float256) = cmp(x, y) < 0
# Base. <(x::Float256, y::Float64) = cmp(x, y) < 0
#
# Base. <=(x::Float256, y::Float256) = cmp(x, y) <= 0
# Base. <=(x::Float64, y::Float256) = cmp(x, y) <= 0
# Base. <=(x::Float256, y::Float64) = cmp(x, y) <= 0
#
# Base. >(x::Float256, y::Float256) = cmp(x, y) > 0
# Base. >(x::Float64, y::Float256) = cmp(x, y) > 0
# Base. >(x::Float256, y::Float64) = cmp(x, y) > 0
#
# Base. >=(x::Float256, y::Float256) = cmp(x, y) >= 0
# Base. >=(x::Float64, y::Float256) = cmp(x, y) >= 0
# Base. >=(x::Float256, y::Float64) = cmp(x, y) >= 0

Base. ==(a::Float256, b::Float256) =
    a.d4[1] == b.d4[1] && a.d4[2] == b.d4[2] &&
    a.d4[3] == b.d4[3] && a.d4[4] == b.d4[4]
Base. ==(a::Float128, b::Float256) = b == a
Base. ==(a::Float64, b::Float256) = b == a
Base. ==(a::Float256, b::Float128) =
    a.d4[1] == b.d2[1] && a.d4[2] == b.d2[2] && a.d4[3] == 0.0 && a.d4[4] == 0.0
Base. ==(a::Float256, b::Float64) =
    a.d4[1] == b && a.d4[2] == 0.0 && a.d4[3] == 0.0 && a.d4[4] == 0.0

Base. !=(a::Float256, b::Float256) = !(a == b)
Base. !=(a::Float128, b::Float256) = !(a == b)
Base. !=(a::Float64, b::Float256) = !(a == b)
Base. !=(a::Float256, b::Float128) = !(a == b)
Base. !=(a::Float256, b::Float64) = !(a == b)

Base. <(a::Float256, b::Float256) =
    a.d4[1] < b.d4[1] ||
    (a.d4[1] == b.d4[1] && (a.d4[2] < b.d4[2] ||
                            (a.d4[2] == b.d4[2] && (a.d4[3] < b.d4[3] ||
                                                    (a.d4[3] == b.d4[3] &&
                                                     a.d4[4] < b.d4[4])))))
Base. <(a::Float128, b::Float256) = b > a
Base. <(a::Float64, b::Float256) = b > a
Base. <(a::Float256, b::Float128) =
    a.d4[1] < b.d2[1] ||
    (a.d4[1] == b.d2[1] && (a.d4[2] < b.d2[2] ||
                            (a.d4[2] == b.d2[2] && a.d4[3] < 0.0)))
Base. <(a::Float256, b::Float64) =
    a.d4[1] < b || (a.d4[1] == b && a.d4[2] < 0.0)

Base. <=(a::Float256, b::Float256) =
    a.d4[1] < b.d4[1] ||
    (a.d4[1] == b.d4[1] && (a.d4[2] < b.d4[2] ||
                            (a.d4[2] == b.d4[2] && (a.d4[3] < b.d4[3] ||
                                                    (a.d4[3] == b.d4[3] &&
                                                     a.d4[4] <= b.d4[4])))))
Base. <=(a::Float128, b::Float256) = b >= a
Base. <=(a::Float64, b::Float256) = b >= a
Base. <=(a::Float256, b::Float128) =
    a.d4[1] < b.d2[1] ||
    (a.d4[1] == b.d2[1] && (a.d4[2] < b.d2[2] ||
                            (a.d4[2] == b.d2[2] && a.d4[3] <= 0.0)))
Base. <=(a::Float256, b::Float64) =
    a.d4[1] < b || (a.d4[1] == b && a.d4[2] <= 0.0)

Base. >(a::Float256, b::Float256) =
    a.d4[1] > b.d4[1] ||
    (a.d4[1] == b.d4[1] && (a.d4[2] > b.d4[2] ||
                            (a.d4[2] == b.d4[2] && (a.d4[3] > b.d4[3] ||
                                                    (a.d4[3] == b.d4[3] &&
                                                     a.d4[4] > b.d4[4])))))
Base. >(a::Float128, b::Float256) = b < a
Base. >(a::Float64, b::Float256) = b < a
Base. >(a::Float256, b::Float128) =
    a.d4[1] > b.d2[1] ||
    (a.d4[1] == b.d2[1] && (a.d4[2] > b.d2[2] ||
                            (a.d4[2] == b.d2[2] && a.d4[3] > 0.0)))
Base. >(a::Float256, b::Float64) =
    a.d4[1] > b || (a.d4[1] == b && a.d4[2] > 0.0)

Base. >=(a::Float256, b::Float256) =
    a.d4[1] > b.d4[1] ||
    (a.d4[1] == b.d4[1] && (a.d4[2] > b.d4[2] ||
                            (a.d4[2] == b.d4[2] && (a.d4[3] > b.d4[3] ||
                                                    (a.d4[3] == b.d4[3] &&
                                                     a.d4[4] >= b.d4[4])))))
Base. >=(a::Float128, b::Float256) = b <= a
Base. >=(a::Float64, b::Float256) = b <= a
Base. >=(a::Float256, b::Float128) =
    a.d4[1] > b.d2[1] ||
    (a.d4[1] == b.d2[1] && (a.d4[2] > b.d2[2] ||
                            (a.d4[2] == b.d2[2] && a.d4[3] >= 0.0)))
Base. >=(a::Float256, b::Float64) =
    a.d4[1] > b || (a.d4[1] == b && a.d4[2] >= 0.0)

function Base.ceil(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_ceil, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.floor(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_floor, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

# Note: This breaks ties the wrong way
function Base.round(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_nint, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.trunc(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_aint, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base. ^(x::Float256, i::Integer)
    r = Ref{D4}()
    ccall((:c_qd_npwr, libqd), Cvoid, (Ref{D4}, Cint, Ref{D4}), x.d4, i, r)
    Float256(r[])
end

# function Base.literal_pow(::typeof(^), x::Float256, ::Val{2})
#     r = Ref{D4}()
#     ccall((:c_qd_sqr, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
#     Float256(r[])
# end

function Base.literal_pow(::typeof(^), a::Float256, ::Val{2})
    # quad-double ^ 2  = (x0 + x1 + x2 + x3) ^ 2
    #                  = x0 ^ 2 + 2 x0 * x1 + (2 x0 * x2 + x1 ^ 2)
    #                             + (2 x0 * x3 + 2 x1 * x2)

    p0, q0 = two_sqr(a.d4[1])
    p1, q1 = two_prod(2.0 * a.d4[1], a.d4[2])
    p2, q2 = two_prod(2.0 * a.d4[1], a.d4[3])
    p3, q3 = two_sqr(a.d4[2])

    p1, q0 = two_sum(q0, p1)

    q0, q1 = two_sum(q0, q1)
    p2, p3 = two_sum(p2, p3)

    s0, t0 = two_sum(q0, p2)
    s1, t1 = two_sum(q1, p3)

    s1, t0 = two_sum(s1, t0)
    t0 += t1

    s1, t0 = quick_two_sum(s1, t0)
    p2, t1 = quick_two_sum(s0, s1)
    p3, q0 = quick_two_sum(t1, t0)

    p4 = 2.0 * a.d4[1] * a.d4[4]
    p5 = 2.0 * a.d4[2] * a.d4[3]

    p4, p5 = two_sum(p4, p5)
    q2, q3 = two_sum(q2, q3)

    t0, t1 = two_sum(p4, q2)
    t1 = t1 + p5 + q3

    p3, p4 = two_sum(p3, t0)
    p4 = p4 + q0 + t1

    p0, p1, p2, p3, p4 = renorm(p0, p1, p2, p3, p4)
    Float256((p0, p1, p2, p3))
end

function Base.sqrt(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_sqrt, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

export nroot
function nroot(x::Float256, i::Integer)
    r = Ref{D4}()
    ccall((:c_qd_nroot, libqd), Cvoid, (Ref{D4}, Cint, Ref{D4}), x.d4, i, r)
    Float256(r[])
end

function Base.exp(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_exp, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

Base.exp10(x::Float256) = exp(c256_log10 * x)

Base.exp2(x::Float256) = exp(c256_log2 * x)

function Base.log(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_log, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.log10(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_log10, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

Base.log2(x::Float256) = c256_1_log2 * log(x)

function Base.sin(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_sin, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.cos(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_cos, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.tan(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_tan, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.asin(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_asin, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.acos(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_acos, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.atan(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_atan, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.atan(x::Float256, y::Float256)
    r = Ref{D4}()
    ccall((:c_qd_atan2, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}),
          x.d4, y.d4, r)
    Float256(r[])
end

Base.hypot(x::Float256, y::Float256) = sqrt(x^2 + y^2)

function Base.sinh(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_sinh, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.cosh(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_cosh, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.tanh(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_tanh, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.asinh(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_asinh, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.acosh(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_acosh, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.atanh(x::Float256)
    r = Ref{D4}()
    ccall((:c_qd_atanh, libqd), Cvoid, (Ref{D4}, Ref{D4}), x.d4, r)
    Float256(r[])
end

function Base.sincos(x::Float256)
    r = Ref{D4}()
    s = Ref{D4}()
    ccall((:c_qd_sincos, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}), x.d4, r, s)
    Float256(r[]), Float256(s[])
end

export sincosh
function sincosh(x::Float256)
    r = Ref{D4}()
    s = Ref{D4}()
    ccall((:c_qd_sincosh, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}),
          x.d4, r, s)
    Float256(r[]), Float256(s[])
end

function Base.rand(::Type{Float256})
    r = Ref{D4}()
    ccall((:c_qd_rand, libqd), Cvoid, (Ref{D4},), r)
    Float256(r[])
end

function Base.rand(::Type{Float256}, dims::Dims)
    r = Array{Float256}(undef, dims)
    for i in LinearIndices(r)
        r[i] = rand(Float256)
    end
    r
end

Base.rand(::Type{Float256}, dims::Integer...) = rand(Float256, Dims(dims))

# Constants

# const c256_0 = Float256(0.0)
# const c256_1 = Float256(1.0)
# const c256_eps = Float256(2^-208)

const c256_pi = begin
    r = Ref{D4}()
    ccall((:c_qd_pi, libqd), Cvoid, (Ref{D4},), r)
    Float256(r[])
end
const c256_e = exp(Float256(1.0))
const c256_log10 = log(Float256(10.0))
const c256_log2 = log(Float256(2.0))
const c256_1_log2 = inv(c256_log2)

ext: jl


text:
import UncertainData.ConstrainedIndexValueResampling

"""
    ConstrainedTest(test::CausalityTest, constraints::ConstrainedResampling, n::Int)
    ConstrainedTest(test::CausalityTest, constraints::ConstrainedResampling)

A causality test where the supports of the uncertain values furnishing the
uncertain values in the datasets are truncated (according to the provided
`constraints`) prior to performing the `test`.

`n` controls the number of independent resamplings over which the test
is performed (if not provided, `n` is set to `1` by default).

## Examples

Assume we want to apply a causality test to two uncertain datasets `X` and `Y`,
but restricting the ranges of values their elements can take while resampling.

```julia
# Constraints for X.indices and X.values
cx = (TruncateQuantiles(0.3, 0.7), TruncateStd(1.5))

# Constraints for Y.indices and Y.values
cy = (TruncateQuantiles(0.3, 0.7), TruncateStd(1.5))

# Need to gather to feed to the ConstrainedTest constructor.
cs = ConstrainedIndexValueResampling(cx, cy)

# A cross mapping test applied to 100 independent realisations
# of `X` and `Y`, resampled after constraining the data.
ctest = ConstrainedTest(CrossMappingTest(), cs, 100)
```
"""
struct ConstrainedTest{
        CT <: CausalityTest,
        CR <: ConstrainedIndexValueResampling{N_VARIABLES, N_DATASETS}
            where {N_VARIABLES, N_DATASETS}} <: PreprocessedDataCausalityTest{CT}
    test::CT
    constrained_resampling::CR
    n::Int
end

function ConstrainedTest(test, constrained_resampling)
    ConstrainedTest(test, constrained_resampling, 1)
end

export ConstrainedTest


text:
using FileIO
using Images

field = [ "s3", "s8", "s14", "s15", "s20","s37", "s40"]
data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200316";
output_dir = "$(data_dir)_visualization"
try
    mkdir("$output_dir")
catch
end

function max_projection(pos)
    println("Processing $pos")
    raw_img = load(File(format"TIFF", "$data_dir/HE7-11-1-80uw-PWM_1_$pos.ome.tiff"));
    img_size = size(raw_img);
    img_max = zeros(N0f16, img_size[1], img_size[2], Int(img_size[3]/20));
    for i in 1:Int(img_size[3]/20)
        img_max[:, :, i] = maximum(raw_img[:, :, (i-1)*20+1:i*20], dims=3);
    end
    save("$output_dir/$pos.tiff", img_max);
end

for pos in field
    max_projection(pos)
end


text:
immutable MVHistory{H<:UnivalueHistory} <: MultivalueHistory
    storage::Dict{Symbol, H}
end

function MVHistory{H<:UnivalueHistory}(::Type{H} = History)
    MVHistory{H}(Dict{Symbol, H}())
end

# ==========================================================================
# Functions

Base.length(history::MVHistory, key::Symbol) = length(history.storage[key])
Base.enumerate(history::MVHistory, key::Symbol) = enumerate(history.storage[key])
Base.first(history::MVHistory, key::Symbol) = first(history.storage[key])
Base.last(history::MVHistory, key::Symbol) = last(history.storage[key])

function Base.push!{I,H<:UnivalueHistory,V}(
        history::MVHistory{H},
        key::Symbol,
        iteration::I,
        value::V)
    if !haskey(history.storage, key)
        _hist = H(V, I)
        push!(_hist, iteration, value)
        history.storage[key] = _hist
    else
        push!(history.storage[key], iteration, value)
    end
    value
end

function Base.push!{H<:UnivalueHistory,V}(
        history::MVHistory{H},
        key::Symbol,
        value::V)
    if !haskey(history.storage, key)
        _hist = H(V, Int)
        push!(_hist, value)
        history.storage[key] = _hist
    else
        push!(history.storage[key], value)
    end
    value
end

function Base.getindex(history::MVHistory, key::Symbol)
    history.storage[key]
end

Base.haskey(history::MVHistory, key::Symbol) = haskey(history.storage, key)

function Base.get(history::MVHistory, key::Symbol)
    l = length(history, key)
    k, v = first(history.storage[key])
    karray = zeros(typeof(k), l)
    varray = Array(typeof(v), l)
    i = 1
    for (k, v) in enumerate(history, key)
        karray[i] = k
        varray[i] = v
        i += 1
    end
    karray, varray
end

function Base.show{H}(io::IO, history::MVHistory{H})
    print(io, "MVHistory{$H}")
    for (key, val) in history.storage
        print(io, "\n", "  :$(key) => $(val)")
    end
end

using Base.Meta

"""
Easily add to a MVHistory object `tr`.

Example:

```julia
using ValueHistories, OnlineStats
v = Variance(BoundedEqualWeight(30))
tr = MVHistory()
for i=1:100
    r = rand()
    fit!(v,r)
    μ,σ = mean(v),std(v)

    # add entries for :r, :μ, and :σ using their current values
    @trace tr i r μ σ
end
"""
macro trace(tr, i, vars...)
    block = Expr(:block)
    for v in vars
        push!(block.args, :(push!($(esc(tr)), $(quot(Symbol(v))), $(esc(i)), $(esc(v)))))
    end
    block
end


text:
const INFILE = "SOCATv5_all.tsv"::String
const OUTDIR = "datasets"::String

function run()

    # Open input file
    inchan = open(INFILE)

    # Tracking variables
    currentdataset = Nullable{String}()
    outchan = Nullable{Any}()

    currentline = readline(inchan)
    linecount = 1

    # Loop until the end of the file
    while length(currentline) > 0

        fields = split(currentline, "\t")
        dataset = fields[1]

        # See if the dataset has changed
        if dataset != currentdataset

            # Close the existing output file
            if !isnull(outchan)
                close(outchan)
            end

            currentdataset = dataset

            # Open a new output file
            outchan = open(OUTDIR * "/" * currentdataset * ".tsv", "w")
            #write(outchan, "date\tlon\tlat\tfco2\n")
            write(outchan, "lon\tlat\tfco2\n")
            print("\033[2K\r $linecount $currentdataset")
        end

        # Write the line to the output file, only using the bits we need
        #year = fields[5]
        #month = fields[6]
        #day = fields[7]
        #hour = fields[8]
        #minute = fields[9]
        #second = fields[10][1:2]
        lon = fields[11]
        lat = fields[12]
        fco2 = fields[24]

#       write(outchan, "$year/$month/$day $hour:$minute:$second\t$lon\t$lat\t$fco2\n")
        write(outchan, "$lon\t$lat\t$fco2\n")

        currentline = readline(inchan)
        linecount += 1
    end

    if !isnull(outchan)
        close(outchan)
    end

    # Close input file
    close(inchan)

    print("\n")
end

@time run()


text:
@testset "thermalization.jl" begin

    @testset "_weights_canonical" begin
        # check weighting
        range = collect(1:15)
        @test XXZNumerics._weights_canonical(log.(range), -1) ≈ range./sum(range)
        @test XXZNumerics._weights_canonical(log.(range), -2) ≈ range.^2 ./ sum(range .^ 2)

        # check extreme cases
        @test .!(any(isnan, XXZNumerics._weights_canonical(exp.(range), 100)))
        @test XXZNumerics._weights_canonical(exp.(range), 0) ≈ ones(length(range))/length(range)
        @test XXZNumerics._weights_canonical(exp.(range), 1) ≈ XXZNumerics._weights_canonical(sort!(-exp.(range)), -1)[end:-1:1]
    end

    @testset "findβ" begin
        @test findβ([1,2,3,4], 2.5) ≈ 0
        @test findβ([1,2,3,4], 3) < 0
        @test findβ([1,2,3,4], 2) > 0
        @test findβ([1,2,3,4], 2.5-0.7) ≈ -findβ([1,2,3,4], 2.5+0.7) # should be symmetric around the average energy as spectrum is symmetric
        # can't think of other good tests...
    end

    @testset "canonical_ensemble" begin
        l = 10
        states = sort!(rand(l)*10 .- 5)
        @test canonical_ensemble(states, states[1]) ≈ vcat(1, zeros(l-1))
        @test canonical_ensemble(states, states[end]) ≈ vcat(zeros(l-1), 1)
        @test canonical_ensemble(states, mean(states)) ≈ ones(l)/l
    end

    @testset "microcanonical_ensemble" begin
        @test microcanonical_ensemble([-2, -1, 0.5, 1, 2], 0, 1.1) ≈ [0,1,1,1,0]/3
        @test microcanonical_ensemble([-2, -1, 0.5, 1, 2], 0, 0.9) ≈ [0,0,1,0,0]
        @test microcanonical_ensemble([-2, -1, 0.5, 1, 2], 0, 0) ≈ [0,0,1,0,0] # choose closest state
        @test microcanonical_ensemble([-2, -1, 0.5, 1, 2], -0.6, 0) ≈ [0,1,0,0,0]

        @test microcanonical_ensemble([-40, -1, 0.5, 1, 40], 0.7) ≈ [0,0,1,1,0]/2 # default arg -> ΔE = 80*0.5% = 0.4
    end
end


text:
using Test
using BulletsMaggots

@test bm_count(1111, 2222) == (0, 0)
@test bm_count(1111, 1234) == (1, 0)
@test bm_count(1234, 4321) == (0, 4)
@test bm_solver(5678, false) > 0


text:
module Memorization

export Network

@enum TieBreak negative positive

include("./Table.jl")
include("./Network.jl")

function numberofneurons(net::Network)

    [sum([ isassigned(layer, tt) for tt in 1:length(layer) ])  for layer in net.layers]

end

end # module


text:
using Base.Test
using StochasticProcesses

let p = FirstTime(BrownianMotion(), (t,x)->x>.1), t=linspace(0, 1, 10)
  sim(p, t, 3)
end


text:
## TODO consolidate low-level part code. I thought newer code was better. Now, it seems the older code is useful

posnegi(x::Mxpr,n::Integer) = n >= 0 ? n : length(x) + n + 1
## WARNING: The following method is not the same as the others, it won't work in the general case (not for sequence specs)
## we may need to do something about this.
## In posnegi(obj,n) in general, obj is the object a part of which we want to get
## In posnegi below nmax is not this object
posnegi(nmax::Integer,n::Integer) = n >= 0 ? n : nmax + n + 1
posnegi(x,n::Integer) = n
posnegi(x,a) = a

# This *is* used. Eg 1[0] --> Integer
# FIXME: throw Symata error, not Julia
function get_part_one_ind(x, ind::Integer)
    ind != 0 && throw(BoundsError("part specification $ind is longer than the depth of the object."))
    return head(x)
end

function get_part_one_ind(texpr::V,ind::Integer) where V<:Union{Mxpr,Array}
    ind = posnegi(texpr,ind)
    texpr = ind == 0 ? mhead(texpr) : texpr[ind]
    return texpr
end
get_part_one_ind(texpr::Dict,tind) = texpr[tind]  # Part 0 can't return "Dict" because it could be a key.
get_part_one_ind(texpr::Tuple,tind) = texpr[tind]
function get_part_one_ind(texpr::Mxpr,tind::Mxpr{:Span})
    spanargs = margs(tind)
    lsp = length(spanargs)
    if lsp == 2
        nargs = view(margs(texpr),spanargs[1]:spanargs[2]) # need function to do this translation
    elseif lsp == 3
        nargs = view(margs(texpr),spanargs[1]:spanargs[3]:spanargs[2])
    end
    texpr = mxpr(mhead(texpr),nargs...) # we need splice to copy Int Array to Any Array
    return texpr
end

### Part

# Get part of expression. Julia :ref is mapped to :Part
# a[i] parses to Part(a,i), and a[i][j] to Part(Part(a,i),j)
# a[i,j] parses to Part(a,i,j).

@mkapprule Part  """
    Part(expr,n) or expr[n]

returns the `n`th element of expression `expr`.

    Part(expr,n1,n2,...) or expr[n1,n2,...]

return the subexpression at index `[n1,n2,...]`.

The same can be achieved less efficiently with `expr[n1][n2]...`.

    expr[n] = val

set the `n`th part of `expr` to `val`. `n` and `val` are evaluated
normally. `expr` is evaluated once.

`n` less than zero counts backward from the end of the expression or subexpression.

`expr[n]` also returns the `n`th element of instances of several
Julia types such as `Array`, or the element with key `n` for `Dict`'s.

`Part` also indexes into some sympy objects, such as sympy Tuples.
"""

@doap Part(texpr,tinds...) = getpart2(texpr,tinds...)

function getpart2(texpr,tinds...)
    for tind in tinds
        texpr = get_part_one_ind(texpr,tind)
    end
    return texpr
end

### Span

@sjdoc Span """
    Span(a,b) or a:b

represents elements `a` through `b`.

    Span(a,b,c) or a:b:c

represents elements `a` through `b` in steps of `c`.

`expr[a:b]` returns elements a through `b` of expr, with the same head as `expr`.
"""

### Position

@mkapprule Position """
    Position(expr,x)

return a list of part specifications (indices) of positions in
`expr` at which `x` occurs.
"""

@doap Position(expr,subx) = tolistoflists(find_positions(expr,subx))
@curry_second Position

"""
    find_positions(ex,subx)

return a list of part indices at which subexpression `subx` occurs in expression `ex`.
"""
function find_positions(ex,subx)
    posns = Array{Array{Int,1}}(undef, 0)
    lev = Array{Int}(undef, 100)
    clev::Int = 1
    lev[clev] = 0
    capt = capturealloc()
    psubx = patterntoBlank(subx)
    _find_positions(ex,psubx,lev,posns,clev,capt)
    return posns
end

# TODO: split out literal pattern code: This routine will be the main routine for
# finding positions, inlcuding when we know that the pattern is literal.
# The literal pattern code is commented out below
function _find_positions(ex::Mxpr,psubx,lev,posns,clev,capt)
    if clev > length(lev) push!(lev,0) end
    args = margs(ex)
    @inbounds for i in 1:length(args)
        lev[clev] = i
        _find_positions(args[i],psubx,lev,posns,clev+1,capt)
    end
    (gotmatch, capt) = match_and_capt(mhead(ex),psubx,capt)
    if gotmatch
#    if mhead(ex) == subx      # Use this for literal pattern
        lev[clev] = 0
        nlev = copy(lev)
        push!(posns,view(nlev,1:clev))
    end
    (gotmatch, capt) = match_and_capt(ex,psubx,capt)
    if gotmatch
#    if ex == subx
        nlev = copy(lev)
        push!(posns,view(nlev,1:clev-1))
    end
end

function _find_positions(ex,psubx,lev,posns,clev,capt)
    (gotmatch, capt) = match_and_capt(ex,psubx,capt)
    if gotmatch
#    if ex == subx
        nlev = copy(lev)
        push!(posns,view(nlev,1:clev-1))
    end
end

#### setpart

## setpart2! does the same, but I think it  is better and simpler and doesnt seem to be broken

## getpart is used only here. it is defined in mxpr.jl. It is probably not fast
## maybe getpart2 above is faster and handles more cases.
"""
    setpart!(mx::Mxpr,val,inds...)

set the part of `mx` indexed by `inds` to `val`.

The last index equal to zero means the head of the subexpression specified by the previous indices.
Because most heads are effectively immutable, we must replace the final subexpression in this case.
"""
function setpart!(mx::Mxpr,val,inds...)
    if length(inds) == 0
#        return mxpr(val,margs(mx))
        return mxprnewhead(mx, val)
    elseif length(inds) == 1
        if inds[1] == 0  # replacing head
#            mx = mxpr(val,margs(mx)...)
            mx = mxprnewhead(mx,val)
            # Function application is very slow if this branch is taken.
            return mx
        else
            margs(mx)[inds[1]] = val
        end
    elseif inds[end] == 0
        ninds = inds[1:end-2]
        if length(ninds) == 0
            symerror("Not implemented getting no part of $mx")
        end
        ex = getpart(mx, ninds...)
        ex1 = margs(ex)[inds[end-1]]
#        margs(ex)[inds[end-1]] = mxpr(val, margs(ex1))
        margs(ex)[inds[end-1]] = mxprnewhead(ex1,val)
    else
        ex = getpart(mx, inds[1:end-1]...)
        margs(ex)[inds[end]] = val
    end
    mx
end

## This is older, but seems better.
## This checks for negative indices. If calling setpart!, this checking must be done before calling
## setpart! is intended for dev code where negative indices won't appear. But, it is probably not
## necessary to have setpart!
##
## A complication arises because the head (Head) of an Mxpr cannot be changed because the head is
## also in the type Mxpr{sym}. If the head needs to be changed, a copy must be made. The exception
## is if the head is not a symbol. In this case, the type is Mxpr{GenHead}. The head is always stored
## in a field as well. It is possible, but undesirable for the type and field disagree.
## This is why the previous level expression is stored in setpart2!
function setpart2!(mx::Mxpr, val, inds...)
    ex = mx
    exlast = ex
    ninds = length(inds)
    for j in 1:ninds-1
        exlast = ex
        ex = get_part_one_ind(ex,inds[j])
    end
    ind = posnegi(ex,inds[ninds])
    if ind == 0
        exlast[inds[ninds-1]] = mxprnewhead(ex,val)
    else
        ex[ind] = val
    end
    val
end

## either setpart! or setpart2! works here.
Base.setindex!(mx::Mxpr, val, inds...) = setpart2!(mx,val,inds...)

#### localize variable

"""
    maybe_localize_variable(var::Symbol,expr)

If `var` is unbound (that is, evaluates to itself), return `(var,expr)`.
If `var` is bound, return `(var1,expr1)`, where `var1` is a new symbol (from `gensym`)
and `expr1` is a copy of `expr` with `var` replaced everywhere by `var1`.
"""
function maybe_localize_variable(var,expr::Mxpr)
    isbound(var) || return (var,expr)
    localize_variable(var,expr)
end

"""
    localize_variable(var::Symbol,expr)

return `(var1,expr1)`, where `var1` is a new symbol (from `gensym`)
and `expr1` is a copy of `expr` with `var` replaced everywhere by `var1`.
"""
function localize_variable(var,expr::Mxpr)
    sym = get_localized_symbol(var)  # not really necessary that named contains var
    positions = find_positions(expr,var)
    nexpr = deepcopy(expr)
    for pos in positions
        nexpr = setpart!(nexpr, sym, pos...)
    end
    (sym, nexpr)
end

function localize_variable(var,expr::SJSym)
    sym = get_localized_symbol(var)
    var == expr && return (sym,sym)
    return (sym, expr)
end

localize_variable(var,expr) = (get_localized_symbol(var), expr)

function localize_variables(vars::Array,expr::Mxpr)
    syms = map(get_localized_symbol, vars)
    nexpr = deepcopy(expr)
    for (sym,var) in zip(syms,vars)
        positions = find_positions(expr,var)
        for pos in positions
            setpart!(nexpr, sym, pos...)
        end
    end
    (syms, nexpr)
end

function localize_variables(vars::Array,expr::SJSym)
    syms = map(get_localized_symbol, vars)
    for (sym,var) in zip(syms,vars)
        var == expr && return (syms,sym)
    end
    return (syms, expr)
end

localize_variables(vars::Array,expr) = (map(get_localized_symbol, vars), expr)

### Take

@sjdoc Take """
    Take(expr,n)

returns `expr` keeping only the first `n` arguments.

    Take(expr,-n)

returns `expr` dropping all but the last `n` arguments.

    Take(expr,[n])

returns `expr` keeping only the `n`th arguments.

    Take(expr,[m,n])

returns `expr` keeping only the `m`th through `n`th arguments.

    Take(expr,[m,n,s])

returns `expr` keeping only the `m`th through `n`th arguments, with step `s`.

    Take(expr,spec1,spec2,...)

apply `spec1` at level `1`, `spec2` at level `2`...
"""
@mkapprule Take
@doap Take(x::Mxpr, rawspecs...) = take(x,map(sequencespec, rawspecs)...)
take(x, onespec, specs...) = mxpr(mhead(x),mapmargs(t -> take(t,specs...), margs(take(x,onespec)))...)
take(x,spec::SequenceN)     = mxpr(mhead(x), margs(x)[spec.n>=0 ? (1:spec.n) : (length(x)+spec.n+1:end)]...)
take(x,spec::SequenceUpToN) = mxpr(mhead(x), margs(x)[1:min(spec.n,length(x))]...)
take(x,spec::SequenceNOnly) = mxpr(mhead(x), margs(x)[posnegi(x,spec.n)])
take(x,spec::SequenceMN)    = mxpr(mhead(x), margs(x)[posnegi(x,spec.m):posnegi(x,spec.n)]...)
take(x,spec::SequenceMNS)   = mxpr(mhead(x), margs(x)[posnegi(x,spec.m):spec.s:posnegi(x,spec.n)]...)
take(x,spec::SequenceNone)  = mxpr(mhead(x))
take(x,spec::SequenceAll)   = mxpr(mhead(x), margs(x)...)

### Drop

# With a little work, we could raise an error if the user tries to drop more elements
# than exists. In fact, the behavior is as if we always use UpTo.
@sjdoc Drop """
    Drop(expr,n)
    Drop(expr,-n)
    Drop(expr,[m,n])

!!! note
    `Drop(expr,UpTo(n))` is not necessary because `Drop(expr,n)` already gives the
    desired behavior.
"""

@mkapprule Drop
@doap Drop(x::Mxpr, rawspecs...) = drop(x,map(sequencespec, rawspecs)...)
drop(x, onespec, specs...) = mxpr(mhead(x),map(t -> drop(t,specs...), margs(drop(x,onespec)))...)
drop(x,spec::SequenceN) = mxpr(mhead(x), margs(x)[spec.n>0 ? (((spec.n)+1):end) : (1:length(x)+spec.n)]...)
drop(x,spec::SequenceUpToN) = drop(x,SequenceN(spec.n))

function drop(x,spec::SequenceNOnly)
    n = posnegi(x,spec.n)
    ma = margs(x)
    mxpr(mhead(x), ma[1:n-1]..., ma[n+1:end]...)
end

function drop(x,spec::SequenceMN)
    m = posnegi(x,spec.m)
    n = posnegi(x,spec.n)
    ma = margs(x)
    mxpr(mhead(x), ma[1:m-1]..., ma[n+1:end]...)
end

### Extract

@sjdoc Extract """
    Extract(expr, list)

returns the subexpression specified by `list`.

    Extract(expr, [list1, list2, ...])

returns a list of subexpressions.

    Extract(expr, [list1, list2, ...], h)

wraps each subexpression with head `h`.

    Extract(list)

operator form of `Extract`.

`list` is a part specification as returned by `Part`.
"""
@mkapprule Extract

@doap function Extract(expr, p::ListT)
    isa(p[1],Integer) && return expr[margs(p)...]
    ! isa(p[1],ListT) && symerror("Extract: Bad part specification ", p)
    nargs = newargs(length(p))
    for i in 1:length(p)
        nargs[i] = expr[margs(p[i])...]
    end
    return MList(nargs)
end

@doap function Extract(expr, p::ListT, h)
    isa(p[1],Integer) && return mxpr(h,expr[margs(p)...])
    ! isa(p[1],ListT) && symerror("Extract: Bad part specification ", p)
    nargs = newargs(length(p))
    for i in 1:length(p)
        nargs[i] = mxpr(h,expr[margs(p[i])...])
    end
    return MList(nargs)
end

@curry_last Extract

### Replace

@sjdoc Replace """
    Replace(expr,rule)

replace parts in expr according to `rule`.

    Replace(expr,rule,n)

replace at levels less than or equal to `n`.

    Replace(expr,rule,[n])

replace at level `n` only.
"""

@sjseealso_group(Replace,ReplaceAll)
@sjexamp(Replace,
         ("Clear(a,b,c)",""),
         ("Replace( Cos(a+b)^2 + Sin(a+c)^2, Cos(x_)^2 + Sin(x_)^2 => 1)",
          "Cos(a + b)^2 + Sin(a + c)^2", "This expression does not match the pattern."),
         ("Replace( Cos(a+b)^2 + Sin(a+b)^2, Cos(x_)^2 + Sin(x_)^2 => 1)",
          "1", "This expression does match the pattern."))

@mkapprule Replace

# function apprules(mx::Mxpr{:Replace})
#     doreplace(mx,margs(mx)...)
# end

# function doreplace(mx,expr,r::Rules)
#     (success, result) = replace(expr,r)
#     result
# end

@doap function Replace(expr,r::Rules)
    (success, result) = replace(expr,r)
    result
end

@doap function Replace(expr,r::Rules,inlevelspec)
    levelspec = make_level_specification(expr,inlevelspec)
    (success, result) = replace(levelspec,expr,r)
    result
end

@curry_second Replace

#### ReplaceAll

@sjdoc ReplaceAll """
    ReplaceAll(expr,rule)

replace parts at all levels in expr according to `rule`. This includes
the `Head` of `expr`, and the entire `expr`.

    ReplaceAll(expr,List(rule1,rule2,...))

replace parts at all levels according to the
list of rules. If given explicitly, the rules should be given as `List(...)` rather than
`[...]` because of a parsing error.

    op = ReplaceAll(rule)

Define function `op(expr)` that returns `ReplaceAll(expr,rule)`.
"""

@mkapprule ReplaceAll

@doap ReplaceAll(expr,r::Rules) =replaceall(expr,r)

@doap function ReplaceAll(expr,rs::Mxpr{:List})
    if listoflistsq(rs)
        maplist( r ->  _doreplaceall(mx,expr,r), rs)
    else
        _doreplaceall(mx,expr,rs::Mxpr{:List})
    end
end

function _doreplaceall(mx,expr,rs::Mxpr{:List})
    rsa = Array{Any}(undef, 0)
    for i in 1:length(rs)
        if isa(rs[i],Rules)
            push!(rsa, rs[i])
        else
            @symwarn("ReplaceAll expected Rule  got ", rs[i])
            nothing  # do something better here, like return mx
        end
    end
    replaceall(expr,rsa)
end

@sjexamp( ReplaceAll,
         ("ClearAll(zz,b,c)",""),
         ("zz = 10 * b^2 * (c+d)","zz = 10 * b^2 * (c+d)"),
         ("ReplaceAll(zz, List(c => 3,d => 2) )", "50*b^2"))

@curry_second ReplaceAll

### ReplaceRepeated

@mkapprule ReplaceRepeated :options => Dict( :MaxIterations => 65536 )

@sjdoc ReplaceRepeated """
    ReplaceRepeated(expr,rules)

perform `ReplaceAll(expr,rules)` repeatedly until `expr` no longer changes.
"""

@doap ReplaceRepeated(expr,r::T; kws...) where {T<:Rules} = replacerepeated(expr,r; kws...)

## FIXME: use get to get keywords
#function do_ReplaceRepeated(mx::Mxpr{:ReplaceRepeated},expr,rs::Mxpr{:List}; kws...)
@doap function ReplaceRepeated(expr,rs::Mxpr{:List}; kws...)
    rsa = Array{Any}(undef, 0)
    for r in rs
        if isa(r,Rules)
            push!(rsa, r)
        else
            @symwarn("ReplaceRepeated expected Rule, got ", r)
            nothing  # do something better here, like return mx
        end
    end
    replacerepeated(expr,rsa; kws...)
end

@doap ReplaceRepeated(a,b; kws...) = mx

### MatchQ

@sjdoc MatchQ """
    MatchQ(expr,pattern)

return `True` if `expr` matches `pattern`.

    op = MatchQ(pattern)

Define a function `op(expr)` that returns `MatchQ(expr,pattern)`.
For example, `myintq = MatchQ(_Integer)`.
"""
@mkapprule MatchQ
@doap MatchQ(expr,pat) = matchq(expr,pat)
matchq(expr,pat) = match_no_capture(expr,patterntoBlank(pat))
@curry_last MatchQ

@sjexamp( MatchQ,
         ("MatchQ( 1, _Integer)", "true"),
         ("ClearAll(gg,xx,b)",""),
         ("MatchQ( gg(xx) , _gg)", "true"),
         ("MatchQ( b^2, _^2)","true"))


text:
import HTTP2
using HTTP2.Frame
using Base.Test

# A server example
HTTP2.serve(8000, b"
Hello, world!
")


text:
Ad = [
   [-1.0/6.0  3.0/6.0 -3.0/6.0 1.0/6.0];
   [ 3.0/6.0 -6.0/6.0  0.0/6.0 4.0/6.0];
   [-3.0/6.0  3.0/6.0  3.0/6.0 1.0/6.0];
   [ 1.0/6.0  0.0/6.0  0.0/6.0 0.0/6.0]
]

dAd = [
   [ 0.0 -0.5  1.0 -0.5];
   [ 0.0  1.5 -2.0  0.0];
   [ 0.0 -1.5  1.0  0.5];
   [ 0.0  0.5  0.0  0.0]
]

d2Ad = [
   [ 0.0 0.0 -1.0  1.0];
   [ 0.0 0.0  3.0 -2.0];
   [ 0.0 0.0 -3.0  1.0];
   [ 0.0 0.0  1.0  0.0]
]

U(s,i) = Symbol(string(s,i))
U(s,i,j) = Symbol(string(s,i,"_",j))

function create_Phi(d, extrap, diff)
    lines = []
    for i=1:d
        block = []
        rhs_1 = U("tp", i,1)
        rhs_2 = U("tp", i,2)
        rhs_4 = U("tp", i,4)
        rhs_3 = U("tp", i,3)
        if extrap == "none"
            for j=1:4
                eq = :($(U("Phi_",i,j)) = ($(Ad[j,1])*$rhs_1 + $(Ad[j,2])*$rhs_2 + $(Ad[j,3])*$rhs_3 + $(Ad[j,4])*$rhs_4) )
                push!(lines,eq)
            end
        elseif extrap == "natural"
            block = quote
                if $(U("t",i))<0
                    $( [ :($(U("Phi_",i,j)) = ($(dAd[j,4])*$rhs_3 + $(Ad[j,4])) ) for j=1:4 ]...)
                elseif $(U("t",i))>1
                    $( [ :($(U("Phi_",i,j)) = (3*$(Ad[j,1]) + 2*$(Ad[j,2]) + $(Ad[j,3]))*($rhs_3-1) + ($(Ad[j,1])+$(Ad[j,2])+$(Ad[j,3])+$(Ad[j,4])) ) for j=1:4 ]...)
                else
                    $( [ :($(U("Phi_",i,j)) = ($(Ad[j,1])*$rhs_1 + $(Ad[j,2])*$rhs_2 + $(Ad[j,3])*$rhs_3 + $(Ad[j,4])*$rhs_4) ) for j=1:4 ]...)
                end
            end
            # for l in block.args
            push!(lines,block.args...)
        end
        if diff
            for j=1:4
                eq = :($(U("dPhi_",i,j)) = ($(dAd[j,1])*$rhs_1 + $(dAd[j,2])*$rhs_2 + $(dAd[j,3])*$rhs_3 + $(dAd[j,3])*$rhs_4 )*$(U("dinv",i)) )
                push!(lines,eq)
            end
        end
    end
    return lines
end

function tensor_prod(symbs, inds, add_index=false)
    if length(symbs)==0
        subscripts = [:($(U("i",i))+$(inds[i])) for i=1:length(inds)]
        if add_index
            return :(C[k,$(subscripts...)])
        else
            return :(C[$(subscripts...)])
        end
    else
        h = symbs[1]
        if length(symbs)>1
            q = symbs[2:end]
        else
            q = []
        end
        exprs = []
        for i = 1:4
            e = Meta.parse( string(h,"_",i,".*(",tensor_prod(q,cat(inds,[i]; dims=1),add_index),")") )
            push!(exprs,e)
        end
        return :(+($(exprs...)))
    end
end

tensor_prod([], Int64[1,2,3,4],)            # C[i1+1,i2+2,i3+3,i4+4]
tensor_prod([], Int64[1,2,3,4], true)      # C[k,i1+1,i2+2,i3+3,i4+4]

tensor_prod(["Phi_1"], Int64[])            # Phi_1_1 * C[i1 + 1] + Phi_1_2 * C[i1 + 2] + Phi_1_3 * C[i1 + 3] + Phi_1_4 * C[i1 + 4]
tensor_prod(["Phi_1", "Phi_2"], Int64[])   # Phi_1_1 * (Phi_2_1 * C[i1 + 1,i2 + 1] + Phi_2_2 * C[i1 + 1,i2 + 2] + Phi_2_3 * C[i1 + 1,i2 + 3] + Phi_2_4 * C[i1 + 1,i2 + 4]) + Phi_1_2 * (Phi_2_1 * C[i1 + 2,i2 + 1] + Phi_2_2 * C[i1 + 2,i2 + 2] + Phi_2_3 * C[i1 + 2,i2 + 3] + Phi_2_4 * C[i1 + 2,i2 + 4]) + Phi_1_3 * (Phi_2_1 * C[i1 + 3,i2 + 1] + Phi_2_2 * C[i1 + 3,i2 + 2] + Phi_2_3 * C[i1 + 3,i2 + 3] + Phi_2_4 * C[i1 + 3,i2 + 4]) + Phi_1_4 * (Phi_2_1 * C[i1 + 4,i2 + 1] + Phi_2_2 * C[i1 + 4,i2 + 2] + Phi_2_3 * C[i1 + 4,i2 + 3] + Phi_2_4 * C[i1 + 4,i2 + 4])
tensor_prod(["Phi_1", "dPhi_2"], Int64[])  # Phi_1_1 * (dPhi_2_1 * C[i1 + 1,i2 + 1] + dPhi_2_2 * C[i1 + 1,i2 + 2] + dPhi_2_3 * C[i1 + 1,i2 + 3] + dPhi_2_4 * C[i1 + 1,i2 + 4]) + Phi_1_2 * (dPhi_2_1 * C[i1 + 2,i2 + 1] + dPhi_2_2 * C[i1 + 2,i2 + 2] + dPhi_2_3 * C[i1 + 2,i2 + 3] + dPhi_2_4 * C[i1 + 2,i2 + 4]) + Phi_1_3 * (dPhi_2_1 * C[i1 + 3,i2 + 1] + dPhi_2_2 * C[i1 + 3,i2 + 2] + dPhi_2_3 * C[i1 + 3,i2 + 3] + dPhi_2_4 * C[i1 + 3,i2 + 4]) + Phi_1_4 * (dPhi_2_1 * C[i1 + 4,i2 + 1] + dPhi_2_2 * C[i1 + 4,i2 + 2] + dPhi_2_3 * C[i1 + 4,i2 + 3] + dPhi_2_4 * C[i1 + 4,i2 + 4])

function create_parameters(d)
    lines = []
    for i=1:d
        block = quote
            $(U("M",i)) = orders[$i]
            $(U("start",i)) =  a[$i]
            $(U("dinv",i)) = (orders[$i]-1.0)/(b[$i]-a[$i])
        end
        for ll in block.args
            push!(lines, ll)
        end
    end
    return lines
end

function create_local_parameters(d)
    lines = []
    for i=1:d
        bl = quote
            $(U("x",i)) = S[n][$i]
            $(U("u",i)) = ($(U("x",i)) - $(U("start",i)))*$(U("dinv",i))
            $(U("i",i)) = (floor(Int,$(U("u",i)) ))
            $(U("i",i)) = max( min($(U("i",i)),$(U("M",i))-2), 0 )
            $(U("t",i)) = $(U("u",i))-$(U("i",i))
            $(U("tp",i,1)) = $(U("t",i))*$(U("t",i))*$(U("t",i))
            $(U("tp",i,2)) = $(U("t",i))*$(U("t",i))
            $(U("tp",i,3)) = $(U("t",i))
            $(U("tp",i,4)) = 1.0;
        end
        for ll in bl.args
            push!(lines, ll)
        end
    end
    return lines
end

function create_function(d,extrap="natural")
    expr = quote
        function $(Symbol(string("eval_UC_spline!")))( a, b, orders, C::Array{T,d}, S::Vector{SVector{d,Float64}}, V::Vector{T}) where d where T
            $(create_parameters(d)...)
            N = size(S,1)
            # @fastmath @inbounds @simd( # doesn't seem to make any difference
            for n=1:N
                $(create_local_parameters(d)...)
                $(create_Phi(d,extrap,false)...)
                V[n] = $( tensor_prod([string("Phi_",i) for i=1:d], Int64[]) )
            end
            # )
        end
    end
    return expr
end


text:
facts("Writemime") do
    workdir = joinpath(tempdir(), "Images")
    if !isdir(workdir)
        mkdir(workdir)
    end
    context("no compression or expansion") do
        A = U8[0.01 0.99; 0.25 0.75]
        fn = joinpath(workdir, "writemime.png")
        open(fn, "w") do file
            writemime(file, MIME("image/png"), grayim(A), minpixels=0, maxpixels=typemax(Int))
        end
        b = convert(Image{Gray{U8}}, load(fn))
        @fact data(b) --> A
    end
    context("small images (expansion)") do
        A = U8[0.01 0.99; 0.25 0.75]
        fn = joinpath(workdir, "writemime.png")
        open(fn, "w") do file
            writemime(file, MIME("image/png"), grayim(A), minpixels=5, maxpixels=typemax(Int))
        end
        b = convert(Image{Gray{U8}}, load(fn))
        @fact data(b) --> A[[1,1,2,2],[1,1,2,2]]
    end
    context("big images (use of restrict)") do
        A = U8[0.01 0.4 0.99; 0.25 0.8 0.75; 0.6 0.2 0.0]
        Ar = restrict(A)
        fn = joinpath(workdir, "writemime.png")
        open(fn, "w") do file
            writemime(file, MIME("image/png"), grayim(A), minpixels=0, maxpixels=5)
        end
        b = convert(Image{Gray{U8}}, load(fn))
        @fact data(b) --> convert(Array{U8}, Ar)
        # a genuinely big image (tests the defaults)
        abig = grayim(rand(UInt8, 1024, 1023))
        fn = joinpath(workdir, "big.png")
        open(fn, "w") do file
            writemime(file, MIME("image/png"), abig, maxpixels=10^6)
        end
        b = convert(Image{Gray{U8}}, load(fn))
        abigui = convert(Array{UFixed8,2}, data(restrict(abig, (1,2))))
        @fact data(b) --> abigui
    end
end


text:
import Base.Iterators: flatten

a = ones(5)
b = [1 3 5]

c = reinterpret(UInt8, a)
println(c)

z1 = map(x -> reinterpret(UInt8, [x]), a)
z2 = map(x -> reverse(reinterpret(UInt8, [x])), a)

println(z1)
println(z2)

y1 = collect(flatten(z1))
y2 = collect(flatten(z2))

println(y1)
println(y2)


text:
#### AuxiliaryFuncs

export ActiveThermoState

using LambertW
export update_dt!
export heaviside

"""
    heaviside(x_1, x_2)

Heaviside function
 - `x_1`
 - `x_2` value specified at edge case `x_1 = 0`
"""
heaviside(x_1, x_2) = x_1 == 0 ? x_2 : typeof(x_1)(x_1 > 0)

"""
    ActiveThermoState(q, aux, k, i)

Returns a `ThermodynamicState` using grid-mean
quantities at element `k`.
"""
@inline function ActiveThermoState(param_set, q, aux, k, i)
    return PhaseEquil_pθq(
        param_set,
        aux[:p_0, k],
        q[:θ_liq, k, i],
        q[:q_tot, k, i],
    )
end

"""
    ActiveThermoState(q, aux, k, i)

Returns a `ThermodynamicState` using grid-mean
quantities at element `k`.
"""
@inline function ActiveThermoState(param_set, q, aux, k::Dual, i)
    return PhaseEquil_pθq.(
        Ref(param_set),
        aux[:p_0, k],
        q[:θ_liq, k, i],
        q[:q_tot, k, i],
    )
end

function update_dt!(grid, params, q, t)
    gm, en, ud, sd, al = allcombinations(q)
    @unpack Δt, Δt_min = params
    u_max = max([q[:w, k, i] for i in ud for k in over_elems(grid)]...)
    Δt = [min(Δt_min, 0.5 * grid.Δz / max(u_max, 1e-10))]
    Δti = [1 / Δt[1]]
    params[:Δt] = Δt
    params[:Δti] = Δti

    # println("----------------------------")
    t[1] += Δt[1]
    percent_done = t[1] / params[:t_end] * 100.0
    # @show t[1], Δt[1], percent_done, params[:t_end]
end

function lamb_smooth_minimum(l, lower_bound, upper_bound)
    leng = size(l)
    xmin = minimum(l)
    lambda0 = max(
        xmin * lower_bound / real(LambertW.lambertw(2.0 / MathConstants.e)),
        upper_bound,
    )

    i = 1
    num = 0
    den = 0
    while (tuple(i) < leng)
        num += l[i] * exp(-(l[i] - xmin) / lambda0)
        den += exp(-(l[i] - xmin) / lambda0)
        i += 1
    end
    smin = num / den

    return smin
end


text:
@testset "evodist" begin
    @testset "estimate_distance" begin
        out = estimate_distance("GTCGATCGACTAGCT", "GACGCTCGACTAGGT")
        @test out == 0.23963672893825572
    end

end


text:
@testset "profile" begin

CUDAdrv.Profile.start()
CUDAdrv.Profile.stop()

CUDAdrv.@profile begin end

end


text:
#= MIT License

Copyright (c) 2020, 2021 Uwe Fechner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE. =#

# data structures for the flight state and the flight log
# functions for creating a demo flight state, demo flight log, loading and saving flight logs
# function se() for reading the settings
# the parameter P is the number of points of the tether, equal to segments+1
# in addition helper functions for working with rotations

"""
    Settings

Flat struct, defining the settings of the Simulator and the Viewer.

$(TYPEDFIELDS)
"""
@with_kw mutable struct Settings @deftype Float64
    project::String       = ""
    log_file::String      = ""
    "file name of the 3D model of the kite for the viewer"
    model::String         = ""
    "name of the kite model to use (KPS3 or KPS4)"
    physical_model::String = ""
    version::Int64 = 1
    "number of tether segments"
    segments::Int64       = 0
    sample_freq::Int64    = 0
    time_lapse            = 0
    zoom                  = 0
    kite_scale            = 1.0
    fixed_font::String    = ""
    abs_tol               = 0.0
    rel_tol               = 0.0
    max_iter::Int64       = 1
    v_reel_out            = 0
    c0                    = 0
    c_s                   = 0
    c2_cor                = 0
    k_ds                  = 0
    "projected kite area            [m^2]"
    area                  = 0
    "kite mass incl. sensor unit     [kg]"
    mass                  = 0
    "height of the kite               [m]"
    height_k              = 0
    alpha_cl::Vector{Float64} = []
    cl_list::Vector{Float64}  = []
    alpha_cd::Vector{Float64} = []
    cd_list::Vector{Float64}  = []
    "width of the kite                [m]"
    width                 = 0
    alpha_zero            = 0
    alpha_ztip            = 0
    "relative nose distance; increasing m_k increases C2 of the turn-rate law"
    m_k                   = 0
    rel_nose_mass         = 0
    "mass of the top particle relative to the sum of top and side particles"
    rel_top_mass          = 0
    "relative side area               [%]"
    rel_side_area         = 0
    "max depower angle              [deg]"
    alpha_d_max           = 0
    "mass of the kite control unit   [kg]"
    kcu_mass              = 0
    "power to steering line distance  [m]"
    power2steer_dist      = 0
    depower_drum_diameter = 0
    depower_offset        = 0
    tape_thickness        = 0
    v_depower             = 0
    v_steering            = 0
    depower_gain          = 0
    steering_gain         = 0
    v_wind                = 0
    v_wind_ref::Vector{Float64} = [] # wind speed vector at reference height
    h_ref                 = 0
    rho_0                 = 0
    z0                    = 0
    profile_law::Int64    = 0
    alpha                 = 0
    cd_tether             = 0
    d_tether              = 0
    d_line                = 0
    "height of the bridle                    [m]"
    h_bridle              = 0
    l_bridle              = 0
    l_tether              = 0
    damping               = 0
    c_spring              = 0
    "density of Dyneema                   [kg/m³]"
    rho_tether            = 0
    "axial tensile modulus of the tether     [Pa]"
    e_tether              = 0
    "initial elevation angle                [deg]"
    elevation             = 0
    "simulation time                   [sim only]"
    depower               = 0
    sim_time              = 0
    "temperature at reference height         [°C]"
    temp_ref              = 0
    "height of groundstation above see level  [m]"
    height_gnd            = 0
    "turbulence intensity relative to Cabau, NL"
    use_turbulence        = 0
    "wind speeds at ref height for calculating the turbulent wind field [m/s]"
    v_wind_gnds::Vector{Float64} = []
    "average height during reel out           [m]"
    avg_height            = 0
    "relative turbulence at the v_wind_gnds"
    rel_turbs::Vector{Float64} = []
    "the expected value of the turbulence intensity at 15 m/s"
    i_ref                 = 0
    "five times the average wind speed in m/s at hub height over the full year    [m/s]"
    v_ref                 = 0
    "grid resolution in z direction                                                 [m]"
    height_step           = 0
    "grid resolution in x and y direction                                           [m]"
    grid_step             = 0
end
const SETTINGS = Settings()

"""
    set_data_path(data_path="")

Set the directory for log and config files.

If called without argument, use the data path of the package to obtain the default settings
when calling se().
"""
function set_data_path(data_path="")
    if data_path==""
        data_path = joinpath(dirname(dirname(pathof(KiteUtils))), "data")
    end
    if data_path != DATA_PATH[1]
        DATA_PATH[1] = data_path
        SETTINGS.segments == 0 # enforce reloading of settings.yaml
    end
end

"""
    load_settings(project="")

Load the project with the given file name. The default project is determined by the content of the file system.yaml .

The project must include the path and the suffix .yaml .
"""
function load_settings(project="")
    SETTINGS.segments=0
    se(project)
end

"""
    copy_settings()

Copy the default settings.yaml and system.yaml files to the folder DATAPATH
(it will be created if it doesn't exist).
"""
function copy_settings()
    if ! isdir(DATA_PATH[1])
        mkdir(DATA_PATH[1])
    end
    src_path = joinpath(dirname(pathof(@__MODULE__)), "..", DATA_PATH[1])
    cp(joinpath(src_path, "settings.yaml"), joinpath(DATA_PATH[1], "settings.yaml"), force=true)
    cp(joinpath(src_path, "system.yaml"), joinpath(DATA_PATH[1], "system.yaml"), force=true)
    chmod(joinpath(DATA_PATH[1], "settings.yaml"), 0o664)
    chmod(joinpath(DATA_PATH[1], "system.yaml"), 0o664)
end

"""
    se()

Getter function for the [`Settings`](@ref) struct.

The default project is determined by the content of the file system.yaml .
"""
function se(project="")
    if SETTINGS.segments == 0
        if project == ""
            # determine which project to load
            dict = YAML.load_file(joinpath(DATA_PATH[1], "system.yaml"))
            SETTINGS.project = dict["system"]["project"]
        end
        # load project from YAML
        dict = YAML.load_file(joinpath(DATA_PATH[1], SETTINGS.project))
        tmp = split(dict["system"]["log_file"], "/")
        SETTINGS.log_file    = joinpath(tmp[1], tmp[2])
        SETTINGS.segments    = dict["system"]["segments"]
        SETTINGS.sample_freq = dict["system"]["sample_freq"]
        SETTINGS.time_lapse  = dict["system"]["time_lapse"]
        SETTINGS.sim_time    = dict["system"]["sim_time"]
        SETTINGS.zoom        = dict["system"]["zoom"]
        SETTINGS.kite_scale  = dict["system"]["kite_scale"]
        SETTINGS.fixed_font  = dict["system"]["fixed_font"]

        SETTINGS.l_tether    = dict["initial"]["l_tether"]
        SETTINGS.v_reel_out  = dict["initial"]["v_reel_out"]
        SETTINGS.elevation   = dict["initial"]["elevation"]
        SETTINGS.depower     = dict["initial"]["depower"]

        SETTINGS.abs_tol     = dict["solver"]["abs_tol"]
        SETTINGS.rel_tol     = dict["solver"]["rel_tol"]
        SETTINGS.max_iter    = dict["solver"]["max_iter"]

        SETTINGS.c0          = dict["steering"]["c0"]
        SETTINGS.c_s         = dict["steering"]["c_s"]
        SETTINGS.c2_cor      = dict["steering"]["c2_cor"]
        SETTINGS.k_ds        = dict["steering"]["k_ds"]

        SETTINGS.alpha_d_max    = dict["depower"]["alpha_d_max"]
        SETTINGS.depower_offset = dict["depower"]["depower_offset"]

        SETTINGS.model         = dict["kite"]["model"]
        SETTINGS.physical_model= dict["kite"]["physical_model"]
        SETTINGS.version       = dict["kite"]["version"]
        SETTINGS.area          = dict["kite"]["area"]
        SETTINGS.rel_side_area = dict["kite"]["rel_side_area"]
        SETTINGS.mass          = dict["kite"]["mass"]
        SETTINGS.height_k      = dict["kite"]["height"]
        SETTINGS.alpha_cl      = dict["kite"]["alpha_cl"]
        SETTINGS.cl_list       = dict["kite"]["cl_list"]
        SETTINGS.alpha_cd      = dict["kite"]["alpha_cd"]
        SETTINGS.cd_list       = dict["kite"]["cd_list"]

        SETTINGS.width         = dict["kps4"]["width"]
        SETTINGS.alpha_zero    = dict["kps4"]["alpha_zero"]
        SETTINGS.alpha_ztip    = dict["kps4"]["alpha_ztip"]
        SETTINGS.m_k           = dict["kps4"]["m_k"]
        SETTINGS.rel_nose_mass = dict["kps4"]["rel_nose_mass"]
        SETTINGS.rel_top_mass  = dict["kps4"]["rel_top_mass"]

        SETTINGS.l_bridle      = dict["bridle"]["l_bridle"]
        SETTINGS.h_bridle      = dict["bridle"]["h_bridle"]
        SETTINGS.d_line        = dict["bridle"]["d_line"]

        SETTINGS.kcu_mass         = dict["kcu"]["kcu_mass"]
        SETTINGS.power2steer_dist = dict["kcu"]["power2steer_dist"]
        SETTINGS.depower_drum_diameter = dict["kcu"]["depower_drum_diameter"]

        SETTINGS.tape_thickness   = dict["kcu"]["tape_thickness"]
        SETTINGS.v_depower        = dict["kcu"]["v_depower"]
        SETTINGS.v_steering       = dict["kcu"]["v_steering"]
        SETTINGS.depower_gain     = dict["kcu"]["depower_gain"]
        SETTINGS.steering_gain    = dict["kcu"]["steering_gain"]

        SETTINGS.cd_tether   = dict["tether"]["cd_tether"]
        SETTINGS.d_tether    = dict["tether"]["d_tether"]
        SETTINGS.damping     = dict["tether"]["damping"]
        SETTINGS.c_spring    = dict["tether"]["c_spring"]
        SETTINGS.rho_tether  = dict["tether"]["rho_tether"]
        SETTINGS.e_tether    = dict["tether"]["e_tether"]

        SETTINGS.v_wind      = dict["environment"]["v_wind"]
        SETTINGS.v_wind_ref  = dict["environment"]["v_wind_ref"]
        SETTINGS.h_ref       = dict["environment"]["h_ref"]
        SETTINGS.rho_0       = dict["environment"]["rho_0"]
        SETTINGS.z0          = dict["environment"]["z0"]
        SETTINGS.alpha       = dict["environment"]["alpha"]
        SETTINGS.profile_law = dict["environment"]["profile_law"]
        SETTINGS.temp_ref    = dict["environment"]["temp_ref"]            # temperature at reference height         [°C]
        SETTINGS.height_gnd  = dict["environment"]["height_gnd"]          # height of groundstation above see level [m]
        SETTINGS.use_turbulence = dict["environment"]["use_turbulence"]   # turbulence intensity relative to Cabau, NL
        SETTINGS.v_wind_gnds = dict["environment"]["v_wind_gnds"]         # wind speeds at ref height for calculating the turbulent wind field [m/s]
        SETTINGS.avg_height  = dict["environment"]["avg_height"]          # average height during reel out          [m]
        SETTINGS.rel_turbs   = dict["environment"]["rel_turbs"]           # relative turbulence at the v_wind_gnds
        SETTINGS.i_ref       = dict["environment"]["i_ref"]               # is the expected value of the turbulence intensity at 15 m/s.
        SETTINGS.v_ref       = dict["environment"]["v_ref"]               # five times the average wind speed in m/s at hub height over the full year    [m/s]
                                                                          # Cabau: 8.5863 m/s * 5.0 = 42.9 m/s
        SETTINGS.height_step = dict["environment"]["height_step"]         # use a grid with 2m resolution in z direction                                 [m]
        SETTINGS.grid_step   = dict["environment"]["grid_step"]           # grid resolution in x and y direction                                         [m]
    end
    return SETTINGS
end


text:
using FastTransforms, ToeplitzMatrices

@testset "BigFloat Toeplitz" begin
    T = Toeplitz(BigFloat[1,2,3,4,5], BigFloat[1,6,7,8,0])
    @test T*ones(BigFloat,5) ≈ [22,24,19,16,15]

    let n = 512
        r = map(BigFloat,rand(n))
        T = Toeplitz(r,[r[1];map(BigFloat,rand(n-1))])
        @test T*ones(BigFloat,n) ≈ Matrix(T)*ones(BigFloat,n)

        T = TriangularToeplitz(BigFloat[1,2,3,4,5],:L)
        @test T*ones(BigFloat,5) ≈ Matrix(T)*ones(BigFloat,5)

        r = map(BigFloat,rand(n))
        T = TriangularToeplitz(r,:L)
        @test T*ones(BigFloat,n) ≈ Matrix(T)*ones(BigFloat,n)

        T = TriangularToeplitz(BigFloat[1,2,3,4,5],:U)
        @test T*ones(BigFloat,5) ≈ Matrix(T)*ones(BigFloat,5)

        r = map(BigFloat,rand(n))
        T = TriangularToeplitz(r,:U)
        @test T*ones(BigFloat,n) ≈ Matrix(T)*ones(BigFloat,n)
    end
end


text:
#=
Names w/o apostrophe
function 1: read name lists, e.g. reaction names, compound names
  --> return an array of strings
function 2: read bounds, e.g. upper/lower bounds
  --> return an array of tuples (compound::String, bound::Float64)
function 3: read S matrix
  --> return an array of tuples
      (compound::String, reaction::String, coefficient::Float64)
Q: can be combined together? Yup
=#

function read_GAMS_format_file(filePath::String)
  returnArr = Array{Any, 1}()
  open(filePath, "r") do file
    lines = readlines(file)[2:(end-1)]
    if (length(split(lines[1])) <= 1)
      # return an array of strings
      read_names_list(lines, returnArr)
    else
      tokens = split(split(lines[1])[1], ".")
      if length(tokens) == 1
        read_bounds_list(lines, returnArr)
      else
        read_S_matrix_list(lines, returnArr)
      end
    end
  end

  return returnArr
end

function read_names_list(rawData::Array, newData::Array)
  for (id, val) in enumerate(rawData)
    val = split(val)
    if length(val) == 0
      # println("empty line found")
    else
      push!(newData, val[1][2: (end-1)])
    end
  end
end

function read_bounds_list(rawData::Array, newData::Array)
  for (id, val) in enumerate(rawData)
    tokens = split(val)
    if length(tokens) < 1
      # println("weird line found")
    else
      push!(newData, (tokens[1][2:end-1], parse(Float64, tokens[2])))
    end
  end
end

function read_S_matrix_list(rawData::Array, newData::Array)
  for (id, val) in enumerate(rawData)
    tk1 = split(val)
    if length(tk1) == 0
      # println("empty line found")
    else
      tk2 = split(tk1[1], ".")
      push!(newData, (tk2[1][2:end-1], tk2[2][2:end-1], parse(Float64, tk1[2])))
    end
  end
end


text:
module TestHarness

import DataFrames
import MacroTools
import Pkg
import Test: @testset

Pkg.activate("..")

anyerrors = false

println("Running tests:")

paths = (
    "query/expression_operations/validate.jl",
    "query/expression_operations/get_alias.jl",
    "query/expression_operations/find_column_names.jl",
    "query/expression_operations/index_column_names.jl",
    "query/expression_operations/rewrite_column_names.jl",
    "query/expression_operations/safe_tuple_name.jl",
    "query/expression_operations/predicates.jl",
    "query/expression_operations/fusion.jl",
    "query/expression_operations/passes/locals.jl",
    "query/expression_operations/passes/lift.jl",
    "query/expression_operations/passes/tvl.jl",
    "query/expression_operations/forms/tuple_form.jl",
    "query/expression_operations/forms/broadcast_form.jl",
    "query/expression_operations/expression/column_names_tuple_expr.jl",
    "query/expression_operations/expression/expression.jl",

    "query/macros/select.jl",
    "query/macros/where.jl",
    "query/macros/group_by.jl",
    "query/macros/aggregate_vector.jl",
    "query/macros/order_by.jl",
    "query/macros/limit.jl",
    "query/macros/inner_join.jl",
    "query/macros/left_join.jl",
    "query/macros/right_join.jl",
    "query/macros/outer_join.jl",

    "logical_nodes/projection.jl",
    "logical_nodes/selection.jl",
    "logical_nodes/group_by.jl",
    "logical_nodes/aggregate_vector.jl",
    "logical_nodes/order_by.jl",
    "logical_nodes/limit.jl",
    "logical_nodes/inner_join.jl",
    "logical_nodes/left_join.jl",
    "logical_nodes/right_join.jl",
    "logical_nodes/outer_join.jl",

    "physical_operations/dataframes/iterators.jl",
    "physical_operations/dataframes/projection.jl",
    "physical_operations/dataframes/selection.jl",
    "physical_operations/dataframes/group_by.jl",
    "physical_operations/dataframes/aggregate_vector.jl",
    "physical_operations/dataframes/order_by.jl",
    "physical_operations/dataframes/limit.jl",
    "physical_operations/dataframes/inner_join.jl",
    "physical_operations/dataframes/left_join.jl",
    "physical_operations/dataframes/right_join.jl",
    "physical_operations/dataframes/outer_join.jl",
)

@testset "All tests" begin
    for path in paths
        try
            include(path)
            println("\t\033[1m\033[32mPASSED\033[0m: $(path)")
        catch e
            global anyerrors = true
            println("\t\033[1m\033[31mFAILED\033[0m: $(path)")
            showerror(stdout, e, backtrace())
            println()
        end
    end
end

if anyerrors
    throw("Tests failed")
end

end


text:
# ---
# title: 1189. Maximum Number of Balloons
# id: problem1189
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table, String
# link: <https://leetcode.com/problems/maximum-number-of-balloons/description/>
# hidden: true
# ---
#
# Given a string `text`, you want to use the characters of `text` to form as
# many instances of the word **" balloon"** as possible.
#
# You can use each character in `text` **at most once**. Return the maximum
# number of instances that can be formed.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/09/05/1536_ex1_upd.JPG)**
#
#
#
#     Input: text = "nlaebolko"
#     Output: 1
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/09/05/1536_ex2_upd.JPG)**
#
#
#
#     Input: text = "loonbalxballpoon"
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: text = "leetcode"
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= text.length <= 10^4`
#   * `text` consists of lower case English letters only.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
UTIL_inf() = 2^31 - 1

function UTIL_concatenate_solutions_vectors(population::Vector{Sol}, offsprings::Vector{Sol})

    debug && DEBUG_feasible_solutions(population)
    debug && DEBUG_feasible_solutions(offsprings)

    result = vcat(population,offsprings)

    return result
end

function UTIL_remove_doublons_from_front(front::Vector{Sol})

    debug && DEBUG_is_front(front)

    new_front = Vector{Sol}(undef, 0)

    for sol in front
        if !(sol in new_front)
            push!(new_front, sol)
        else
            println("REMOVE DOUBLON")
        end
    end

    debug && DEBUG_is_front(new_front)
    debug && DEBUG_unicity(new_front)

    return new_front
end

function UTIL_termination_criteria_met(it::Int, nb_it::Int)
    return it > nb_it
end

function UTIL_termination_criteria_met(t::Float64, max_t::Float64)
    return t > max_t
end


text:
using BenchmarkTools
using LinearAlgebra
using Newton
using NLsolve
using StableRNGs

LinearAlgebra.BLAS.set_num_threads(1)

function create_rf!(dim, rng=StableRNG(0))
    A = -rand(rng, dim, dim)
    b = -rand(rng, dim)

    function rf!(r, x)
        r .= b
        mul!(r, A, x, 1, 1)
        @simd for i=1:length(x)
            @inbounds r[i] += exp(x[i]) - x[i]^2
        end
        return r
    end

    return rf!
end

dim = 20
rf = create_rf!(dim)
r, x = [zeros(dim) for _ in 1:2]

# Check that rf! is non-allocating
println("@btime rf!")
@btime rf(r_, x_) setup=(r_=copy(r); x_=copy(x)) evals=1;

# Get time for newtonsolve!
# First time setup
cache = NewtonCache(x, rf)
drdx = get_drdx(cache)
println("@btime newtonsolve!")
@btime newtonsolve!(x_, $drdx, rf, $cache) setup=(x_=copy(x)) evals=1;

# Get time for nlsolve
println("@btime nlsolve")
@btime nlsolve(rf, x_, autodiff=:forward) setup=(x_=copy(x)) evals=1;

println("Benchmark (dim=$dim) complete")


text:
# LIST OF OPERATORNAMES
functions_generate = Tuple{String, Int64, Int64}[
    ("Square", 2, 2),
    ("SquareOctagon", 2, 2),
    ("ShastrySutherland", 2, 2),
    ("Triangular", 2, 2),
    ("Honeycomb", 2, 2),
    ("Kagome", 2, 2),

    ("Cubic", 3, 3),
    ("FCC", 3, 3),
    ("Diamond", 3, 3),
    ("Pyrochlore", 3, 3),

    ("Hyperkagome", 3, 3),

    ("8_3_a", 3, 3),
    ("8_3_b", 3, 3),
    ("8_3_c", 3, 3),
    ("8_3_n", 3, 3),

    ("9_3_a", 3, 3),

    ("10_3_a", 3, 3),
    ("10_3_b", 3, 3),
    ("10_3_c", 3, 3),
    ("10_3_d", 3, 3)
]

# go through all elements and generate interface functions
for operatorpair in functions_generate

    # get D and N
    D = operatorpair[2]
    N = operatorpair[3]
    # generate a global name for getUnitcellXXX Syntax
    if isdigit(operatorpair[1][1])
        function_name = "getUnitcell_"*lowercase(filter(c -> isascii(c), operatorpair[1]))
        function_id   = lowercase(filter(c -> isascii(c), operatorpair[1]))
    else
        function_name = "getUnitcell" * filter(c -> isascii(c) && isletter(c), uppercase(operatorpair[1][1]) * operatorpair[1][2:end])
        function_id   = lowercase(filter(c -> isascii(c) && isletter(c), operatorpair[1]))
    end

    # evaluate global export
    eval( quote
    export $(Symbol(function_name))
    end )

    # global interface type (tmp)
    function_id_val = Val{Symbol(function_id)}
    # evaluate global interface scheme
    eval( quote
    function getUnitcell(
                identifier :: $(function_id_val),
                args...
            )

        # return the fitting honeycomb function
        return $(Symbol(function_name))(args...)
    end
    end )

    # evaluating type specific master functions

    # Referencing to individual functions by wrapping in Val()
    eval( quote
    function $(Symbol(function_name))(
                unitcell_type  :: Type{U},
                implementation :: Int64 = 1
            ) :: U where {LS,LB,S<:AbstractSite{LS,$D},B<:AbstractBond{LB,$N},U<:AbstractUnitcell{S,B}}

        # call the respective subfunction by converting to val type
        return $(Symbol(function_name))(unitcell_type, Val(implementation))
    end
    end )

    # Fallback for all implementations (if Val{V} is not found)
    eval( quote
    function $(Symbol(function_name))(
                unitcell_type  :: Type{U},
                implementation :: Val{V}
            ) :: U where {LS,LB,S<:AbstractSite{LS,$D},B<:AbstractBond{LB,$N},U<:AbstractUnitcell{S,B},V}

        # fallback / fail due to missing implementation
        error("Implementation " * string(V) * $(" of " * function_id * " unitcell (label types ") * string(LS) * " / " * string(LB) * ") not implemented yet")
    end
    end )

    # evaluating WRAPPER FUNCTIONS (for concrete Unitcell type) call general function

    # wrapper function for passing no label types (and implementation) (DEFAULT)
    eval( quote
    function $(Symbol(function_name))(
                implementation :: Int64 = 1
            ) :: Unitcell{Site{Int64,$D},Bond{Int64,$N}}

        # create a suitable unitcell of the Unitcell type
        return $(Symbol(function_name))(Unitcell{Site{Int64,$D},Bond{Int64,$N}}, implementation)
    end
    end )

    # wrapper function for passing common label type (and implementation)
    eval( quote
    function $(Symbol(function_name))(
                label_type     :: Type{L},
                implementation :: Int64 = 1
            ) :: Unitcell{Site{L,$D},Bond{L,$N}} where L

        # create a suitable unitcell of the Unitcell type
        return $(Symbol(function_name))(Unitcell{Site{L,$D},Bond{L,$N}}, implementation)
    end
    end )

    # wrapper function for passing site / bond label types (and implementation)
    eval( quote
    function $(Symbol(function_name))(
                label_type_site :: Type{LS},
                label_type_bond :: Type{LB},
                implementation  :: Int64 = 1
            ) :: Unitcell{Site{LS,$D},Bond{LB,$N}} where {LS,LB}

        # create a suitable unitcell of the Unitcell type
        return $(Symbol(function_name))(Unitcell{Site{LS,$D},Bond{LB,$N}}, implementation)
    end
    end )

end

# SYNONYMS
getUnitcellHyperhoneycomb = getUnitcell_10_3_b
getUnitcellHyperoctagon   = getUnitcell_10_3_a
getUnitcellSC             = getUnitcellCubic

export getUnitcellHyperhoneycomb, getUnitcellHyperoctagon, getUnitcellSC