text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8]);F = [81*x1^2*x3^2+306*x1^2*x3*x7+5545*x1^2*x4^2+8836*x1^2*x5^2+289*x1^2*x7^2+x1^2*x8^2+13348*x1*x2^2*x5+262*x1*x2*x4^2-2698*x1*x2*x4*x5-6768*x1*x2*x4*x8+124*x1*x2*x6*x8-1098*x1*x3^2*x6+54*x1*x3^2*x8+950*x1*x3*x4*x6-2074*x1*x3*x6*x7+102*x1*x3*x7*x8-50*x1*x3*x8^2+5904*x1*x4^2*x5-2204*x1*x4*x5*x8+12032*x1*x5*x6^2+13608*x1*x5*x6*x8-132*x1*x5*x8^2+5041*x2^4+1061*x2^2*x4^2+4402*x2^2*x4*x5-940*x2^2*x4*x8+5041*x2^2*x5^2+12932*x2^2*x6^2+10366*x2^2*x6*x8+2209*x2^2*x8^2-1550*x2*x3*x4*x6-3550*x2*x3*x5*x6-3100*x2*x3*x6*x8+820*x2*x4^2*x5-258*x2*x4*x5*x8+8236*x2*x5^2*x8-7192*x2*x5*x6^2-8184*x2*x5*x6*x8+8836*x3^4+10904*x3^3*x6-8836*x3^2*x4*x5-18048*x3^2*x5*x8+7710*x3^2*x6^2-366*x3^2*x6*x8+634*x3^2*x8^2-5452*x3*x4*x5*x6-11136*x3*x5*x6*x8+3300*x3*x5*x8^2+3890*x4^2*x5^2+9024*x4*x5^2*x8+3364*x5^2*x6^2+7656*x5^2*x6*x8+16936*x5^2*x8^2+4096*x6^4+9344*x6^3*x8+5329*x6^2*x8^2-558*x1^2*x3-3420*x1^2*x4-1054*x1^2*x7-5076*x1*x2*x4+31016*x1*x2*x5+234*x1*x3*x5-718*x1*x3*x6-184*x1*x3*x8-3882*x1*x5*x7+10440*x1*x5*x8+13774*x2^3-1480*x2^2*x4-3266*x2^2*x7+6956*x2^2*x8+124*x2*x3*x6-6068*x2*x4*x5+12416*x2*x6^2+14162*x2*x6*x8+2068*x3^2*x4-50*x3^2*x8+1276*x3*x4*x6-1702*x3*x5*x6-54*x3*x5*x8-1034*x4^2*x5-2112*x4*x5*x8-2944*x6^2*x7-3358*x6*x7*x8+9061*x1^2+4608*x1*x4-806*x1*x5+14885*x2^2+640*x2*x4-4462*x2*x7-3008*x2*x8-9775*x3^2-6032*x3*x6+121*x4^2+7512*x4*x5+169*x5^2+9984*x5*x8+529*x7^2-4736*x2-1144*x4+3728, 26696*x1*x2*x5+262*x1*x4^2-2698*x1*x4*x5-6768*x1*x4*x8+124*x1*x6*x8+20164*x2^3+2122*x2*x4^2+8804*x2*x4*x5-1880*x2*x4*x8+10082*x2*x5^2+25864*x2*x6^2+20732*x2*x6*x8+4418*x2*x8^2-1550*x3*x4*x6-3550*x3*x5*x6-3100*x3*x6*x8+820*x4^2*x5-258*x4*x5*x8+8236*x5^2*x8-7192*x5*x6^2-8184*x5*x6*x8-5076*x1*x4+31016*x1*x5+41322*x2^2-2960*x2*x4-6532*x2*x7+13912*x2*x8+124*x3*x6-6068*x4*x5+12416*x6^2+14162*x6*x8+29770*x2+640*x4-4462*x7-3008*x8-4736, 162*x1^2*x3+306*x1^2*x7-2196*x1*x3*x6+108*x1*x3*x8+950*x1*x4*x6-2074*x1*x6*x7+102*x1*x7*x8-50*x1*x8^2-1550*x2*x4*x6-3550*x2*x5*x6-3100*x2*x6*x8+35344*x3^3+32712*x3^2*x6-17672*x3*x4*x5-36096*x3*x5*x8+15420*x3*x6^2-732*x3*x6*x8+1268*x3*x8^2-5452*x4*x5*x6-11136*x5*x6*x8+3300*x5*x8^2-558*x1^2+234*x1*x5-718*x1*x6-184*x1*x8+124*x2*x6+4136*x3*x4-100*x3*x8+1276*x4*x6-1702*x5*x6-54*x5*x8-19550*x3-6032*x6, 11090*x1^2*x4+524*x1*x2*x4-2698*x1*x2*x5-6768*x1*x2*x8+950*x1*x3*x6+11808*x1*x4*x5-2204*x1*x5*x8+2122*x2^2*x4+4402*x2^2*x5-940*x2^2*x8-1550*x2*x3*x6+1640*x2*x4*x5-258*x2*x5*x8-8836*x3^2*x5-5452*x3*x5*x6+7780*x4*x5^2+9024*x5^2*x8-3420*x1^2-5076*x1*x2-1480*x2^2-6068*x2*x5+2068*x3^2+1276*x3*x6-2068*x4*x5-2112*x5*x8+4608*x1+640*x2+242*x4+7512*x5-1144, 17672*x1^2*x5+13348*x1*x2^2-2698*x1*x2*x4+5904*x1*x4^2-2204*x1*x4*x8+12032*x1*x6^2+13608*x1*x6*x8-132*x1*x8^2+4402*x2^2*x4+10082*x2^2*x5-3550*x2*x3*x6+820*x2*x4^2-258*x2*x4*x8+16472*x2*x5*x8-7192*x2*x6^2-8184*x2*x6*x8-8836*x3^2*x4-18048*x3^2*x8-5452*x3*x4*x6-11136*x3*x6*x8+3300*x3*x8^2+7780*x4^2*x5+18048*x4*x5*x8+6728*x5*x6^2+15312*x5*x6*x8+33872*x5*x8^2+31016*x1*x2+234*x1*x3-3882*x1*x7+10440*x1*x8-6068*x2*x4-1702*x3*x6-54*x3*x8-1034*x4^2-2112*x4*x8-806*x1+7512*x4+338*x5+9984*x8, 124*x1*x2*x8-1098*x1*x3^2+950*x1*x3*x4-2074*x1*x3*x7+24064*x1*x5*x6+13608*x1*x5*x8+25864*x2^2*x6+10366*x2^2*x8-1550*x2*x3*x4-3550*x2*x3*x5-3100*x2*x3*x8-14384*x2*x5*x6-8184*x2*x5*x8+10904*x3^3+15420*x3^2*x6-366*x3^2*x8-5452*x3*x4*x5-11136*x3*x5*x8+6728*x5^2*x6+7656*x5^2*x8+16384*x6^3+28032*x6^2*x8+10658*x6*x8^2-718*x1*x3+124*x2*x3+24832*x2*x6+14162*x2*x8+1276*x3*x4-1702*x3*x5-5888*x6*x7-3358*x7*x8-6032*x3, 306*x1^2*x3+578*x1^2*x7-2074*x1*x3*x6+102*x1*x3*x8-1054*x1^2-3882*x1*x5-3266*x2^2-2944*x6^2-3358*x6*x8-4462*x2+1058*x7, 2*x1^2*x8-6768*x1*x2*x4+124*x1*x2*x6+54*x1*x3^2+102*x1*x3*x7-100*x1*x3*x8-2204*x1*x4*x5+13608*x1*x5*x6-264*x1*x5*x8-940*x2^2*x4+10366*x2^2*x6+4418*x2^2*x8-3100*x2*x3*x6-258*x2*x4*x5+8236*x2*x5^2-8184*x2*x5*x6-18048*x3^2*x5-366*x3^2*x6+1268*x3^2*x8-11136*x3*x5*x6+6600*x3*x5*x8+9024*x4*x5^2+7656*x5^2*x6+33872*x5^2*x8+9344*x6^3+10658*x6^2*x8-184*x1*x3+10440*x1*x5+6956*x2^2+14162*x2*x6-50*x3^2-54*x3*x5-2112*x4*x5-3358*x6*x7-3008*x2+9984*x5];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
imu = -6 ± 0.4

u = complex(5 ± 0.3, imu)
v = complex(imu, 8 ± 0.9)

@testset "Addition and subtraction" begin
    @test @inferred(2u - u - 3v + 2v + pi) ≈ @inferred(u - v + pi)
    @test @inferred(u + v) == @inferred(v + u)
    for a in (u, v, z)
        @test @inferred(a + 0) == @inferred(0 + a) == a
        @test @inferred(a - 0) == @inferred(-(0 - a)) == a
        @test @inferred(a + float(im)) == @inferred(float(im) + a)
    end
end

@testset "Multiplication and power" begin
    @test @inferred(v*v*u + u*u*v) ≈ @inferred(u^2*v + v^2*u)
    @test @inferred(u * v) == @inferred(v * u)
    for a in (u, v, z)
        @test @inferred(a * 1) == @inferred(1 * a) == a
        @test @inferred(a ^ 1) == a
        @test @inferred(a ^ 0) == one(a)
        @test @inferred(a * im) == @inferred(im * a)
    end
end

@testset "Division" begin
    @test @inferred(v/v + u/u) ≈ complex(2 ± 0) atol = 5e-18
    for a in (u, v, z)
        @test @inferred(a / 1) == a
        @test @inferred(a / im) == @inferred(-im * a)
    end
end

@testset "Inverse" begin
    @test @inferred(inv(v) + inv(u)) ≈ @inferred((v + u)/(v*u))
end

@testset "Trig functions" begin
    for a in (u, v)
        @test @inferred(cos(acos(a)))   ≈ a
        @test @inferred(cosh(acosh(a))) ≈ a
        @test @inferred(sin(asin(a)))   ≈ a
        @test @inferred(sinh(asinh(a))) ≈ a
        @test @inferred(tan(a))         ≈ @inferred(sin(a)/cos(a))
        @test @inferred(tan(atan(a)))   ≈ a
        @test @inferred(tanh(a))        ≈ @inferred(sinh(a)/cosh(a))
        @test @inferred(tanh(atanh(a))) ≈ a
        @test @inferred(sinpi(a))       ≈ sin(pi * a)
        @test @inferred(cospi(a))       ≈ cos(pi * a)
        @test @inferred(sinc(a))        ≈ (sin(pi * a) / (pi * a))
        @test @inferred(cosc(a))        ≈ ((cos(pi * a) - sin(pi * a) / (pi * a)) / a)
    end
end

@testset "Exponential and logarithm" begin
    @test @inferred(u ^ v) ≈ @inferred(exp(v * log(u)))
    for a in (u, v, z)
        @test @inferred(exp(a))          ≈ ℯ^a
        @test @inferred(expm1(a))        ≈ exp(a) - 1
        @test @inferred(exp10(a))        ≈ 10^a
        @test @inferred(exp2(a))         ≈ 2^a
        @test @inferred(exp(log(a)))     ≈ a
        @test @inferred(log(ℯ, a))       ≈ log(a)
        @test @inferred(log(10, a))      ≈ log10(a)
        @test @inferred(log(2, a))       ≈ log2(a)
        @test @inferred(log1p(a))        ≈ log(1 + a)
        @test @inferred(exp2(log2(a)))   ≈ a
        @test @inferred(exp10(log10(a))) ≈ a
        @test @inferred(expm1(log1p(a))) ≈ a
    end
end

@testset "Square root" begin
    for a in (u, v, z)
        @test @inferred(sqrt(a))         ≈ @inferred(a^0.5)
        @test @inferred(sqrt(a)*sqrt(a)) ≈ a
        @test @inferred(sqrt(big(a)))    ≈ sqrt(a)
    end
end

@testset "Absolute value" begin
    for a in (u, v, z); @test @inferred(abs(a)) ≈ @inferred(sqrt(abs2(a))); end
end

@testset "Representation" begin
    # Make sure the printed representation of a Measurement object is correctly parsed as
    # the same object (well, the tag will be different, but that's not important here).
    for a in (u, v, z); @test eval(Meta.parse(repr(a))) == a; end
end

@testset "Other" begin
    ⊕(a::Complex{<:Measurement}, b::Integer) =
        Measurements.result(Measurements.value(a) + b, (1, 0, 0, 1), a)
    ⊗(a::Complex{<:Measurement}, b::Integer) =
        Measurements.result(Measurements.value(a) * b, (b, 0, 0, b), a)
    for a in (u, v, z)
        @test a ⊕ 5 == a + 5
        @test a ⊗ 5 == a * 5
    end
end

# TODO: Define a method for lgamma with complex argument.
# # Factorial and gamma
# for a in (u, v); test_approx_eq(gamma(a), factorial(a - 1)); end
# for a in (u, v); test_approx_eq(gamma(a + 1), factorial(a)); end
# for a in (u, v); test_approx_eq(lgamma(a), log(gamma(a))); end

ext: jl


text:
import Pkg
Pkg.add("PyPlot")
using PyPlot

max_iters = 100
width = 1000
height = 1000
inf = 4
z_start = complex(0.5)

x = range(-2, stop=2, length=width)
y = range(-2, stop=2, length=height)
x = repeat(x',width,1)
y = repeat(y,1,height)

function julia(X, Y)
    z = z_start
    c = complex(X,Y)
    for i in 0:max_iters
        z = z * z + c
        if abs(z) > inf
            return i
        end
    end
    return 0
end
z = julia.(x, y)

contour(x, y, z)
show()

max_stars_repo_name: gaurapanasenko/experiments

max_issues_repo_name: gaurapanasenko/experiments

max_forks_repo_name: gaurapanasenko/experiments


text:
function linshrink(F::AbstractMatrix, S::AbstractMatrix, λ::Real)
    return Symmetric((1.0 .- λ).*S .+ λ.*F)
end

function linshrink(F::UniformScaling, S::AbstractMatrix, λ::Real)
    return Symmetric((1.0 .- λ).*S + λ.*F)
end

function linshrink!(F::AbstractMatrix, S::AbstractMatrix, λ::Real)
    F .= (1.0 .- λ).*S .+ λ.*F
    return Symmetric(F)
end


text:
module PkgDemo

export add_nums
export cat_str

include("add_nums.jl")
include("cat_str.jl")
end


text:
import DataStructures

"""
    strongly_connected_components(V::Array, E::Array{Tuple})

Strongly connected components of a directed graph in reverse topological
order. Translated from sympy/utilities/iteratbles.py/strongly_connected_components.

"""
function strongly_connected_components(V, E)
    Gmap = Dict(vi => [] for vi in V)
    for (v1, v2) in E
        push!(Gmap[v1], v2)
    end
    # Non-recursive Tarjan's algorithm:
    lowlink = Dict()
    indices = Dict()
    stack = DataStructures.OrderedDict() #########
    callstack = []
    components = []
    next_index = 1

    function start(v)
        index = length(stack) + 1
        indices[v] = lowlink[v] = index
        stack[v] = missing
        push!(callstack, (v, Gmap[v]))  #####
    end

    function finish(v1)
        # Finished a component?
        if lowlink[v1] == indices[v1]
            component = [pop!(stack)[1]]

            while component[end] != v1
                push!(component, pop!(stack)[1])
            end
            push!(components, reverse(component))
        end
        v2, _ = pop!(callstack)
        if length(callstack) > 0
            v1, _ = callstack[end]
            lowlink[v1] = min(lowlink[v1], lowlink[v2])
        end
    end

    for v in V
        if v in keys(indices)
            continue
        end
        start(v)
        while length(callstack) > 0
            v1, it1 = callstack[end]

            if length(it1) == 0
                v2 = missing
            else
                v2 = popfirst!(it1)  ######
            end
            # Finished children of v1?
            if v2 === missing
                finish(v1)
            # Recurse on v2
            elseif !(v2 in keys(indices))
                start(v2)
            elseif v2 in keys(stack)
                lowlink[v1] = min(lowlink[v1], indices[v2])
            end
        end
    end
    # Reverse topological sort order:
    return components
end

"""
Computes the period of an irreducible transition matrix
T must be a 2D array
"""
function breadth_first_search(T)
    n = size(T)[1]

    non_tree_edge_values = Set()
    visited = Set(1)

    overall_newly_visited = Set(1)
    level = Dict(1 => 0)
    current_level = 0
    done = false  # imitate a do-while loop
    while !done  # runs at most len(class_) times
        done = length(visited) == n
        current_level += 1
        new_overall_newly_visited = Set()

        # this loop and the while loop above run a combined len(class_) number of times.
        # so this triple nested loop runs through each of the n states once.
        for i in overall_newly_visited

            # the loop below runs len(class_) number of times
            # complexity is around about O(n * avg(len(class_)))
            newly_visited = Set(j for j in 1:n if T[i, j] != 0)
            new_overall_newly_visited = union(new_overall_newly_visited, newly_visited)

            new_tree_edges = setdiff(newly_visited, visited)
            for j in new_tree_edges
                level[j] = current_level
            end

            new_non_tree_edges = intersect(newly_visited, visited)
            new_non_tree_edge_values = Set(level[i]-level[j]+1 for j in new_non_tree_edges)

            non_tree_edge_values = union(non_tree_edge_values, new_non_tree_edge_values)
            visited = union(visited, new_tree_edges)
        end
        overall_newly_visited = new_overall_newly_visited
    end

    positive_ntev = Set(val_e for val_e in non_tree_edge_values if val_e > 0)
    if length(positive_ntev) == 0
        return n
    else
        return gcd(positive_ntev...)
    end
end


text:
export Fourier,Taylor,Hardy,CosSpace,SinSpace,Laurent

for T in (:CosSpace,:SinSpace)
    @eval begin
        struct $T{D<:PeriodicDomain,R} <: Space{D,R}
            domain::D
            $T{D,R}(d::PeriodicDomain) where {D,R} = new(convert(D,d))
            $T{D,R}(d::D) where{D,R} = new(d)
        end
        $T(d::PeriodicDomain) = $T{typeof(d),real(prectype(d))}(d)
        $T(d) = $T(convert(PeriodicDomain, d))
        $T() = $T(PeriodicSegment())
        spacescompatible(a::$T,b::$T) = domainscompatible(a,b)
        hasfasttransform(::$T) = true
        canonicalspace(S::$T) = Fourier(domain(S))
        setdomain(S::$T,d::PeriodicDomain) = $T(d)
    end
end

"""
`CosSpace()` is the space spanned by `[1,cos θ,cos 2θ,...]`
"""
CosSpace()

"""
`SinSpace()` is the space spanned by `[sin θ,sin 2θ,...]`
"""
SinSpace()

# s == true means analytic inside, taylor series
# s == false means anlytic outside and decaying at infinity
"""
`Hardy{false}()` is the space spanned by `[1/z,1/z^2,...]`.
`Hardy{true}()` is the space spanned by `[1,z,z^2,...]`.
"""
struct Hardy{s,D<:PeriodicDomain,R} <: Space{D,R}
    domain::D
    Hardy{s,D,R}(d) where {s,D,R} = new{s,D,R}(d)
    Hardy{s,D,R}() where {s,D,R}  = new{s,D,R}(D())
end

# The <: PeriodicDomain is crucial for matching Basecall overrides
"""
`Taylor()` is the space spanned by `[1,z,z^2,...]`.
This is a type alias for `Hardy{true}`.
"""
const Taylor{D<:PeriodicDomain,R} = Hardy{true,D,R}

@containsconstants CosSpace
@containsconstants Taylor

Base.promote_rule(::Type{Fun{S,V}},::Type{T}) where {T<:Number,S<:Union{Hardy{true},CosSpace},V} =
    Fun{S,promote_type(V,T)}
Base.promote_rule(::Type{Fun{S}},::Type{T}) where {T<:Number,S<:Union{Hardy{true},CosSpace}} =
    Fun{S,T}

(H::Type{Hardy{s}})(d::PeriodicDomain) where {s} = Hardy{s,typeof(d),complex(prectype(d))}(d)
(H::Type{Hardy{s}})() where {s} = Hardy{s}(Circle())

canonicalspace(S::Hardy) = S
setdomain(S::Hardy{s},d::PeriodicDomain) where {s} = Hardy{s}(d)

spacescompatible(a::Hardy{s},b::Hardy{s}) where {s} = domainscompatible(a,b)
hasfasttransform(::Hardy) = true

for (Typ,Plfft!,Plfft,Pltr!,Pltr) in ((:TransformPlan,:plan_fft!,:plan_fft,:plan_transform!,:plan_transform),
                           (:ITransformPlan,:plan_ifft!,:plan_ifft,:plan_itransform!,:plan_itransform))
    @eval begin
        $Pltr!(sp::Hardy,x::AbstractVector{T}) where {T<:Complex} = $Typ(sp,$Plfft!(x),Val{true})
        $Pltr!(::Hardy,x::AbstractVector{T}) where {T<:Real} =
            error("In place variants not possible with real data.")

        $Pltr(sp::Hardy,x::AbstractVector{T}) where {T<:Complex} = $Typ(sp,$Pltr!(sp,x),Val{false})
        function $Pltr(sp::Hardy,x::AbstractVector{T}) where T
            plan = $Pltr(sp,Array{Complex{T}}(undef,length(x))) # we can reuse vector in itransform
            $Typ{T,typeof(sp),false,typeof(plan)}(sp,plan)
        end

        *(P::$Typ{T,HS,false},vals::AbstractVector{T}) where {T<:Complex,HS<:Hardy} = P.plan*copy(vals)
        *(P::$Typ{T,HS,false},vals::AbstractVector{T}) where {T,HS<:Hardy} = P.plan*Vector{Complex{T}}(vals)
    end
end

*(P::TransformPlan{T,Hardy{true,DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR} =
    lmul!(one(T)/length(vals),P.plan*vals)
*(P::ITransformPlan{T,Hardy{true,DD,RR},true},cfs::AbstractVector{T}) where {T,DD,RR} =
    lmul!(length(cfs),P.plan*cfs)
*(P::TransformPlan{T,Hardy{false,DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR} =
    lmul!(one(T)/length(vals),reverse!(P.plan*vals))
*(P::ITransformPlan{T,Hardy{false,DD,RR},true},cfs::AbstractVector{T}) where {T,DD,RR} =
    lmul!(length(cfs),P.plan*reverse!(cfs))

transform(sp::Hardy,vals::AbstractVector,plan) = plan*vals
itransform(sp::Hardy,vals::AbstractVector,plan) = plan*vals

evaluate(f::AbstractVector,S::Taylor{D,R},z) where {D<:PeriodicDomain,R} = horner(f,fromcanonical(Circle(),tocanonical(S,z)))
function evaluate(f::AbstractVector,S::Taylor{D,R},z) where {D<:Circle,R}
    z=mappoint(S,𝕌,z)
    d=domain(S)
    horner(f,z)
end

function evaluate(f::AbstractVector,S::Hardy{false,D},z) where D<:PeriodicDomain
    z=mappoint(S,𝕌,z)
    z=1/z
    z*horner(f,z)
end
function evaluate(f::AbstractVector,S::Hardy{false,D},z) where D<:Circle
    z=mappoint(S,𝕌,z)
    z=1/z
    z*horner(f,z)
end

##TODO: fast routine

function horner(c::AbstractVector,kr::AbstractRange{Int},x)
    T = promote_type(eltype(c),eltype(x))
    if isempty(c)
        return zero(x)
    end

    ret = zero(T)
    @inbounds for k in reverse(kr)
        ret = muladd(x,ret,c[k])
    end

    ret
end

function horner(c::AbstractVector,kr::AbstractRange{Int},x::AbstractVector)
    n,T = length(x),promote_type(eltype(c),eltype(x))
    if isempty(c)
        return zero(x)
    end

    ret = zeros(T,n)
    @inbounds for k in reverse(kr)
        ck = c[k]
        @simd for i = 1:n
            ret[i] = muladd(x[i],ret[i],ck)
        end
    end

    ret
end

horner(c::AbstractVector,x) = horner(c,1:length(c),x)
horner(c::AbstractVector,x::AbstractArray) = horner(c,1:length(c),x)
horner(c::AbstractVector,kr::AbstractRange{Int},x::AbstractArray) = reshape(horner(c,kr,vec(x)),size(x))

## Cos and Sin space

points(sp::CosSpace,n) = points(domain(sp),2n-2)[1:n]

plan_transform(::CosSpace,x::AbstractVector) = plan_chebyshevtransform(x;kind=2)
plan_itransform(::CosSpace,x::AbstractVector) = plan_ichebyshevtransform(x;kind=2)
transform(::CosSpace,vals,plan) = plan*vals
itransform(::CosSpace,cfs,plan) = plan*cfs

function evaluate(f::AbstractVector,S::CosSpace,t)
    if t ∈ domain(S)
        clenshaw(Chebyshev(),f,cos(tocanonical(S,t)))
    else
        zero(cfstype(f))
    end
end

points(sp::SinSpace,n)=points(domain(sp),2n+2)[2:n+1]

for (Typ,Pltr!,Pltr) in ((:TransformPlan,:plan_transform!,:plan_transform),
                         (:ITransformPlan,:plan_itransform!,:plan_itransform))
    @eval begin
        $Pltr!(sp::SinSpace{DD},x::AbstractVector{T}) where {DD,T<:fftwNumber} =
            $Typ(sp,plan_r2r!(x,RODFT00),Val{true})
        $Pltr(sp::SinSpace{DD},x::AbstractVector{T}) where {DD,T<:fftwNumber} =
            $Typ(sp,$Pltr!(sp,x),Val{false})
        $Pltr!(sp::SinSpace{DD},x::AbstractVector{T}) where {DD,T} =
            error("transform for SinSpace only implemented for fftwNumbers")
        $Pltr(sp::SinSpace{DD},x::AbstractVector{T}) where {DD,T} =
            error("transform for SinSpace only implemented for fftwNumbers")

        *(P::$Typ{T,SinSpace{D,R},false},vals::AbstractVector{T}) where {T,D,R} = P.plan*copy(vals)
    end
end

*(P::TransformPlan{T,SinSpace{DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR} =
    lmul!(one(T)/(length(vals)+1),P.plan*vals)
*(P::ITransformPlan{T,SinSpace{DD,RR},true},cfs::AbstractVector{T}) where {T,DD,RR} =
    lmul!(one(T)/2,P.plan*cfs)

transform(sp::SinSpace,vals::AbstractVector,plan) = plan*vals
itransform(sp::SinSpace,vals::AbstractVector,plan) = plan*vals

evaluate(f::AbstractVector,S::SinSpace,t) = sineshaw(f,tocanonical(S,t))

## Laurent space
"""
`Laurent()` is the space spanned by the complex exponentials
```
    1,exp(-im*θ),exp(im*θ),exp(-2im*θ),…
```
See also `Fourier`.
"""
const Laurent{DD,RR} = SumSpace{Tuple{Hardy{true,DD,RR},Hardy{false,DD,RR}},DD,RR}

##FFT That interlaces coefficients

plan_transform!(sp::Laurent{DD,RR},x::AbstractVector{T}) where {DD,RR,T<:Complex} =
    TransformPlan(sp,plan_fft!(x),Val{true})
plan_itransform!(sp::Laurent{DD,RR},x::AbstractVector{T}) where {DD,RR,T<:Complex} =
    ITransformPlan(sp,plan_ifft!(x),Val{true})

plan_transform!(sp::Laurent{DD,RR},x::AbstractVector{T}) where {DD,RR,T} =
    error("In place variants not possible with real data.")
plan_itransform!(sp::Laurent{DD,RR},x::AbstractVector{T}) where {DD,RR,T} =
    error("In place variants not possible with real data.")

plan_transform(sp::Laurent{DD,RR},x::AbstractVector{T}) where {T<:Complex,DD,RR} =
    TransformPlan(sp,plan_transform!(sp,x),Val{false})
plan_itransform(sp::Laurent{DD,RR},x::AbstractVector{T}) where {T<:Complex,DD,RR} =
    ITransformPlan(sp,plan_itransform!(sp,x),Val{false})

function plan_transform(sp::Laurent{DD,RR},x::AbstractVector{T}) where {T,DD,RR}
    plan = plan_transform(sp,Array{Complex{T}}(undef, length(x))) # we can reuse vector in itransform
    TransformPlan{T,typeof(sp),false,typeof(plan)}(sp,plan)
end
function plan_itransform(sp::Laurent{DD,RR},x::AbstractVector{T}) where {T,DD,RR}
    plan = plan_itransform(sp,Array{Complex{T}}(undef, length(x))) # we can reuse vector in itransform
    ITransformPlan{T,typeof(sp),false,typeof(plan)}(sp,plan)
end

function *(P::TransformPlan{T,Laurent{DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR}
    n = length(vals)
    vals = lmul!(inv(convert(T,n)),P.plan*vals)
    reverseeven!(interlace!(vals,1))
end

function *(P::ITransformPlan{T,Laurent{DD,RR},true},cfs::AbstractVector{T}) where {T,DD,RR}
    n = length(cfs)
    reverseeven!(cfs)
    cfs[:]=[cfs[1:2:end];cfs[2:2:end]]  # TODO: deinterlace!
    lmul!(n,cfs)
    P.plan*cfs
end

*(P::TransformPlan{T,Laurent{DD,RR},false},vals::AbstractVector{T}) where {T<:Complex,DD,RR} =
    P.plan*copy(vals)
*(P::TransformPlan{T,Laurent{DD,RR},false},vals::AbstractVector{T}) where {T,DD,RR} =
    P.plan*Vector{Complex{T}}(vals)
*(P::ITransformPlan{T,Laurent{DD,RR},false},vals::AbstractVector{T}) where {T<:Complex,DD,RR} =
    P.plan*copy(vals)
*(P::ITransformPlan{T,Laurent{DD,RR},false},vals::AbstractVector{T}) where {T,DD,RR} =
    P.plan*Vector{Complex{T}}(vals)

transform(::Laurent{DD,RR},vals,plan) where {DD,RR} = plan*vals
itransform(::Laurent{DD,RR},cfs,plan) where {DD,RR} = plan*cfs

transform(sp::Laurent{DD,RR},vals::AbstractVector) where {DD,RR} = plan_transform(sp,vals)*vals
itransform(sp::Laurent{DD,RR},cfs::AbstractVector) where {DD,RR} = plan_itransform(sp,cfs)*cfs

function evaluate(f::AbstractVector,S::Laurent{DD,RR},z) where {DD,RR}
    z = mappoint(domain(S),Circle(),z)
    invz = 1/z
    horner(f,1:2:length(f),z) + horner(f,2:2:length(f),invz).*invz
end

function Base.conj(f::Fun{Laurent{DD,RR}}) where {DD,RR}
    ncoefficients(f) == 0 && return f

    cfs = Array{cfstype(f)}(undef, iseven(ncoefficients(f)) ? ncoefficients(f)+1 : ncoefficients(f))
    cfs[1] = conj(f.coefficients[1])
    cfs[ncoefficients(f)] = 0
    for k=2:2:ncoefficients(f)-1
        cfs[k] = conj(f.coefficients[k+1])
    end
    for k=3:2:ncoefficients(f)+1
        cfs[k] = conj(f.coefficients[k-1])
    end
    Fun(space(f),cfs)
end

## Fourier space

"""
`Fourier()` is the space spanned by the trigonemtric polynomials
```
    1,sin(θ),cos(θ),sin(2θ),cos(2θ),…
```
See also `Laurent`.
"""
const Fourier{DD,RR} = SumSpace{Tuple{CosSpace{DD,RR},SinSpace{DD,RR}},DD,RR}

Laurent(d::PeriodicDomain) = Laurent{typeof(d),complex(prectype(d))}(d)
Fourier(d::PeriodicDomain) = Fourier{typeof(d),real(prectype(d))}(d)

for Typ in (:Laurent,:Fourier)
    @eval begin
        $Typ() = $Typ(PeriodicSegment())
        $Typ(d) = $Typ(convert(PeriodicDomain, d))

        hasfasttransform(::$Typ{D,R}) where {D,R} = true
    end
end

Laurent(S::Fourier{DD,RR}) where {DD,RR} = Laurent(domain(S))
Fourier(S::Laurent{DD,RR}) where {DD,RR} = Fourier(domain(S))

for T in (:CosSpace,:SinSpace)
    @eval begin
        # override default as canonicalspace must be implemented
        maxspace(::$T,::Fourier{D,R}) where {D,R} = NoSpace()
        maxspace(::Fourier{D,R},::$T) where {D,R} = NoSpace()
    end
end

points(sp::Fourier{D,R},n) where {D,R}=points(domain(sp),n)

struct IFourierTransformPlan{T,SP,PL} <: AbstractTransformPlan{T}
    space::SP
    plan::PL
    work::Vector{T}
end

plan_transform!(sp::Fourier{D,R},x::AbstractVector{T}) where {T<:fftwNumber,D,R} =
    TransformPlan(sp,plan_r2r!(x, R2HC),Val{true})
plan_itransform!(sp::Fourier{D,R},x::AbstractVector{T}) where {T<:fftwNumber,D,R} =
    IFourierTransformPlan(sp, plan_r2r!(x, HC2R), copy(x)) # copy(x) is workspace data

for (Typ,Pltr!,Pltr) in ((:TransformPlan,:plan_transform!,:plan_transform),
                         (:ITransformPlan,:plan_itransform!,:plan_itransform))
    @eval begin
        $Pltr(sp::Fourier{DD,RR},x::AbstractVector{T}) where {T<:fftwNumber,DD,RR} =
            $Typ(sp,$Pltr!(sp,x),Val{false})
        $Pltr!(sp::Fourier{DD,RR},x::AbstractVector{T}) where {T,DD,RR} =
            error("transform for Fourier only implemented for fftwNumbers")
        $Pltr(sp::Fourier{DD,RR},x::AbstractVector{T}) where {T,DD,RR} =
            error("transform for Fourier only implemented for fftwNumbers")

        *(P::$Typ{T,Fourier{DD,RR},false},vals::AbstractVector{T}) where {T,DD,RR} = P.plan*copy(vals)
        *(P::$Typ{T,Fourier{DD,RR},false},vals::AbstractVector{<:Complex}) where {T,DD,RR} =
                P.plan*real(vals) .+ im.*(P.plan*imag(vals))
    end
end

mul!(cfs::AbstractVector{T}, P::TransformPlan{T,Fourier{DD,RR},true}, vals::AbstractVector{T}) where {T,DD,RR} =
    P*copyto!(cfs, vals)

mul!(cfs::AbstractVector{T}, P::TransformPlan{T,Fourier{DD,RR},false}, vals::AbstractVector{T}) where {T,DD,RR} =
    P.plan*copyto!(cfs, vals)

function *(P::TransformPlan{T,Fourier{DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR}
    n = length(vals)
    cfs = lmul!(convert(T,2)/n,P.plan*vals)
    cfs[1] /= 2
    if iseven(n)
        cfs[n÷2+1] /= 2
    end

    negateeven!(reverseeven!(interlace!(cfs,1)))
end

mul!(vals::AbstractVector{T}, P::IFourierTransformPlan{T,Fourier{DD,RR}}, cfs::AbstractVector{T}) where {T,DD,RR} =
    P*copyto!(vals, cfs)

mul!(vals::AbstractVector{T}, P::ITransformPlan{T,Fourier{DD,RR},false},cfs::AbstractVector{T}) where {T,DD,RR} =
    mul!(vals, P.plan, cfs)

function *(P::IFourierTransformPlan{T,Fourier{DD,RR}},cfs::AbstractVector{T}) where {T,DD,RR}
    n = length(cfs)
    reverseeven!(negateeven!(cfs))

    # allocation free version of
    #    cfs[:] = [view(cfs,1:2:n); view(cfs, 2:2:n)]
    P.work .= cfs
    j = 1
    for k=1:2:n
        cfs[j] = P.work[k]
        j+=1
    end
    for k=2:2:n
        cfs[j] = P.work[k]
        j+=1
    end

    if iseven(n)
        cfs[n÷2+1] *= 2
    end
    cfs[1] *= 2
    P.plan*lmul!(inv(convert(T,2)),cfs)
end

transform(sp::Fourier{DD,RR},vals::AbstractVector,plan) where {DD,RR} = plan*vals
itransform(sp::Fourier{DD,RR},cfs::AbstractVector,plan) where {DD,RR} = plan*cfs

transform(sp::Fourier{DD,RR},vals::AbstractVector) where {DD,RR} = plan_transform(sp,vals)*vals
itransform(sp::Fourier{DD,RR},cfs::AbstractVector) where {DD,RR} = plan_itransform(sp,cfs)*cfs

canonicalspace(S::Laurent{DD,RR}) where {DD<:PeriodicSegment,RR} = Fourier(domain(S))
canonicalspace(S::Fourier{DD,RR}) where {DD<:Circle,RR} = Laurent(domain(S))
canonicalspace(S::Laurent{DD,RR}) where {DD<:PeriodicLine,RR} = S

## Ones and zeros

for sp in (:Fourier,:CosSpace,:Laurent,:Taylor)
    @eval begin
        one(::Type{T},S::$sp{DD,RR}) where {T<:Number,DD,RR} = Fun(S,fill(one(T),1))
        one(S::$sp{DD,RR}) where {DD,RR} = Fun(S,fill(1.0,1))
    end
end

function Fun(::typeof(identity), S::Taylor{DD,RR}) where {DD<:Circle,RR}
    d=domain(S)
    if d.orientation
        Fun(S,[d.center,d.radius])
    else
        error("Cannot create identity on $S")
    end
end

Fun(::typeof(identity), S::Fourier{DD,RR}) where {DD<:Circle,RR} =
    Fun(Fun(identity, Laurent(domain(S))),S)

reverseorientation(f::Fun{Fourier{DD,RR}}) where {DD,RR} =
    Fun(Fourier(reverseorientation(domain(f))),alternatesign!(copy(f.coefficients)))
function reverseorientation(f::Fun{Laurent{DD,RR}}) where {DD,RR}
    # exp(im*k*x) -> exp(-im*k*x), or equivalentaly z -> 1/z
    n=ncoefficients(f)
    ret=Array{cfstype(f)}(undef, iseven(n) ? n+1 : n)  # since z -> 1/z we get one more coefficient
    ret[1]=f.coefficients[1]
    for k=2:2:length(ret)-1
        ret[k+1]=f.coefficients[k]
    end
    for k=2:2:n-1
        ret[k]=f.coefficients[k+1]
    end
    iseven(n) && (ret[n] = 0)

    Fun(Laurent(reverseorientation(domain(f))),ret)
end

include("calculus.jl")
include("specialfunctions.jl")
include("FourierOperators.jl")
include("LaurentOperators.jl")
include("LaurentDirichlet.jl")


text:
"""
data_populate(c=.5; half=false, level=5)
"""
function data_populate(c=.5; half=false, level=5)
fmtplot = ("k-", "k--", "k-.", "k-.", "k>:")
dimbase=1024
aymin=1.e-15
half ? (T=Float16) : (T=Float64)
c==1.0 ? (maxit=30) : (maxit=10)
if T == Float16
DataC=zeros(maxit+1,level,2)
BuildHist!(DataC,dimbase, c, level, maxit, T)
else
DataC=zeros(maxit+1,level,4)
for TI in (Float64, Float32)
    BuildHist!(DataC,dimbase, c, level, maxit, TI)
end
end
fname=Fname(c,T)
writempdata(fname,DataC)
end

function BuildHist!(DataC,dimbase, c, pmax, maxit, T)
if T==Float16
pstart=1
else
pstart = (T==Float32)*3 + (T==Float64)*1
end
for p=0:pmax-1
n=dimbase*(2^p)
dout=heqtest(n, c, "no"; jprecision=T, jmaxit=maxit)
DataC[:,p+1,pstart].=dout.exactout.history
DataC[:,p+1,pstart+1].=dout.fdout.history
end
end

max_stars_repo_name: ctkelley/MPResults


text:
using SeisIO.SUDS

vfname        = "SampleFiles/Restricted/10081701.WVP"
sac_fstr      = "SampleFiles/Restricted/20081701.sac-0*"
eq_wvm1_pat   = "SampleFiles/SUDS/*1991191072247.wvm1.sac"
eq_wvm1_file  = "SampleFiles/SUDS/eq_wvm1.sud"
eq_wvm2_pat   = "SampleFiles/SUDS/*1991191072247.wvm2.sac"
eq_wvm2_file  = "SampleFiles/SUDS/eq_wvm2.sud"
lsm_pat       = "SampleFiles/SUDS/*1992187065409.sac"
lsm_file      = "SampleFiles/SUDS/lsm.sud"
rot_pat       = "SampleFiles/SUDS/*1993258220247.sac"
rot_file      = "SampleFiles/SUDS/rotate.sud"

printstyled("  SUDS\n", color=:light_green)

printstyled("    readsudsevt\n", color=:light_green)
# read into SeisEvent
W = SUDS.readsudsevt(rot_file)
W = SUDS.readsudsevt(lsm_file)
W = SUDS.readsudsevt(eq_wvm1_file)
W = SUDS.readsudsevt(eq_wvm2_file)

printstyled("    in read_data\n", color=:light_green)
if safe_isfile(vfname)
  redirect_stdout(out) do
    SUDS.suds_support()
    S = read_data("suds", vfname, v=3, full=true)
    S = read_data("suds", eq_wvm1_file, v=3, full=true)
    S = read_data("suds", eq_wvm2_file, v=3, full=true)
    S = read_data("suds", lsm_file, v=3, full=true)
    S = read_data("suds", rot_file, v=3, full=true)
    S = read_data("suds", "SampleFiles/SUDS/eq_wvm*sud", v=3, full=true)
  end
end

printstyled("    equivalence to SAC readsuds\n", color=:light_green)
# Volcano-seismic event supplied by W. McCausland
if has_restricted
  S = SUDS.read_suds(vfname, full=true)
  S2 = read_data("sac", sac_fstr, full=true)
  for n = 1:S2.n
    id = S2.id[n]
    if startswith(id, ".")
      id = "OV"*id
    end
    i = findid(id, S)
    if i > 0
      @test isapprox(S.x[i], S2.x[n])
      @test isapprox(Float32(S2.fs[n]), Float32(S.fs[i]))
      @test abs(S.t[i][1,2] - S2.t[i][1,2]) < 2000 # SAC only has ~1 ms precision
    end
  end
else
  printstyled("      skipped volcano-seismic test file\n", color=:red)
end

# SUDS sample files
# from eq_wvm1.sud
S = SUDS.read_suds(eq_wvm1_file)
S1 = read_data("sac", eq_wvm1_pat)
for n = 1:S1.n
  i = findid(S1.id[n], S)
  if i > 0
    @test isapprox(S.x[i], S1.x[n])
  else
    @warn(string(S1.id[n], " not found in S; check id conversion!"))
  end
end

# from eq_wvm2.sud
S = SUDS.read_suds(eq_wvm2_file)
S1 = read_data("sac", eq_wvm2_pat)
for n = 1:S1.n
  i = findid(S1.id[n], S)
  if i > 0
    @test isapprox(S.x[i], S1.x[n])
  else
    @warn(string(S1.id[n], " not found in S; check id conversion!"))
  end
end

# from lsm.sud
S = SUDS.read_suds(lsm_file)
S1 = read_data("sac", lsm_pat)
for n = 1:S1.n
  i = findid(S1.id[n], S)
  if i > 0
    @test isapprox(S.x[i], S1.x[n])
  else
    @warn(string(S1.id[n], " not found in S; check id conversion!"))
  end
end

# from rot.sud
S = SUDS.read_suds(rot_file)
S1 = read_data("sac", rot_pat)
for n = 1:S1.n
  i = findid(S1.id[n], S)
  if i > 0
    @test isapprox(S.x[i], S1.x[n])
  else
    @warn(string(S1.id[n], " not found in S; check id conversion!"))
  end
end

nothing


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

#
# This implements Agha's example 3.2.1
#

using Actors, Printf
import Actors: spawn, newLink

mutable struct StackNode{T}
    content::T
    link::Link
end

struct Pop <: Msg
    customer::Link
end

struct Push{T} <: Msg
    content::T
end

struct Print <: Msg end

forward!(lk::L, msg::M) where {L<:Link, M<:Msg} = send(lk, msg)

# ----- this is essentially Agha's code example written in Julia/YAActL
function stack_node(sn::StackNode, msg::Pop)
    isnothing(sn.content) || become(forward!, sn.link)
    send(msg.customer, Response(sn.content))
end

function stack_node(sn::StackNode, msg::Push)
    P = spawn(Bhv(stack_node, sn))
    become(stack_node, StackNode(msg.content, P))
end

mystack = spawn(Bhv(stack_node, StackNode(nothing, newLink())))

response = newLink()

send(mystack, Pop(response))  # new stack
receive(response).y           # returns nothing
send(mystack, Push(1))        # push 1 on the stack
send(mystack, Pop(response))  # pop it
receive(response).y           # returns 1, 1st node now forwards messages
send(mystack, Pop(response)); # pop again
receive(response).y           # now returns nothing

for i ∈ 1:5
    send(mystack, Push(i))
end

# (send(mystack, Print()); sleep(0.1))

for i ∈ 1:5
    send(mystack, Pop(response))
    println(receive(response).y)
end


text:
struct MinChanges
    image::Matrix{UInt8}
    cover::Matrix{UInt8}
    distortion::Int
end

MinChanges(image::Array{UInt8}) = MinChanges(deepcopy(image),deepcopy(image),0)
getindex(image::MinChanges,i::Int64) = image.image[i]

function setindex!(a::MinChanges,v,i::Int64)
    a.distortion -= a.image[i] != a.cover[i]
    a.image[i] = v
    a.distortion += a.image[i] != a.cover[i]
end

function LSBcostfun(image::MinChanges,pidx)
    if image.image[pidx]&0x01==0x01
        return(image.distortion + 1,0)
    else
        return(0,image.distortion + 1)
    end
end

function copy!(dest::MinChanges,src::MinChanges)
    copy!(dest.image,src.image)
    dest.distortion = src.distortion
    return(dest)
end

readlsb(image::MinChanges,embpath) = readlsb(image.image,embpath)

function readlsb(img::Matrix{UInt8},embpath)
    y=zeros(UInt8,length(embpath))
    for (i,ii) in enumerate(embpath)
        y[i] = img[ii] & 0x1
    end
    return(y)
end

function matchlsb!(img::MinChanges,mbit,pidx::Int)
    original=img[pidx]
    if img[pidx]&0x1!=mbit
        if original==255
            img[pidx]=original-1
            return
        end
        if original==0
            img[pidx]=original+1
            return
        end
        if rand()>0.5
            img[pidx]=original-1
        else
            img[pidx]=original+1
        end
    end
    @assert(mod(img[pidx],2) == mbit)
end


text:
using PolynomialZeros
using Polynomials
using Test
using LinearAlgebra

include("test-amvw.jl")
include("test-multroot.jl")

## Some Polynomial Familes
function Wilkinson(n, T=Float64)
    x = variable(T)
    prod(x-i for i in 1:n)
end

function Chebyshev(n, T=Float64)
    x = variable(T)
    n == 0 && return one(x)
    n == 1 && return x
    2 * x * Chebyshev(n-1, T) - Chebyshev(n-2, T)
end

function Legendre(n, T=Float64)
    x = variable(T)
    n == 0 && return one(x)
    n == 1 && return x
    1/(n+1) * ( (2n+1)*x*Legendre(n-1, T) - n * Legendre(n-2, T))
end

## only a few
function Cyclotomic(n, T=Float64)
    x = variable(T)
    n == 1 && return x-1
    n == 2 && return x + 1
    n == 5 && return x^4 + x^3 + x^2 + x + 1
    n == 10 && return x^4 - x^3 + x^2 -x + 1
    n == 20 && return x^8 - x^6 + x^4 - x^2 + 1
    throw(DomainError)
end

@testset "Over.C" begin
    @test length(poly_roots(Wilkinson(10), Over.C)) == 10
    @test length(poly_roots(x -> x^5 - x -1, Over.C)) == 5
    @test maximum(norm.(poly_roots(Chebyshev(5), Over.C))) <= 1

    p = Cyclotomic(20)
    rts = poly_roots(p, Over.C)
    @test maximum(norm.(p.(rts) ./ polyder(p).(rts))) <= 1e-14

    ## we have different methods
    rts = 1.0:6; p = poly(rts)
    fn = x -> x^5 - x - 1
    for m in [:PolynomialRoots, :roots, :amvw]
        @test maximum(norm.(sort(poly_roots(p, Over.C, method=m), by=norm) .- rts)) <= 1e-10
        poly_roots(fn, Over.C, method=m)
    end

    # test deflation
    fn = x -> x^3 * (x-1) * (x^2 + 1)
    @test sum(iszero.(poly_roots(fn, Over.C))) == 3

end

@testset "Over.R" begin
    @test length(poly_roots(Wilkinson(10), Over.R)) == 10
    @test length(poly_roots(x -> x^5 - x -1, Over.R)) == 1
    @test norm(poly_roots(x -> x^5 - x -1, Over.R)[1] - 1.1673039782614187) <= 1e-14
    @test maximum(norm.(poly_roots(Chebyshev(5), Over.R))) <= 1

    p = Legendre(11)
    rts = poly_roots(p, Over.R, square_free=true) # needs help here
    @test maximum(p.(rts) ./ polyder(p).(rts)) <= 1e-14

    p = poly([1.0, 2, 2, 3, 3, 3])
    rts = poly_roots(p, Over.R)
    @test length(rts) == 3

    # test deflation
    fn = x -> x^3 * (x-1) * (x^2 + 1)
    ds = sort(poly_roots(fn, Over.R)) .- sort([1,0])
    @test norm(ds) <= 1e-14

end

@testset "Over.Q" begin
    @test length(poly_roots(Wilkinson(5, Int), Over.Q)) == 5
    @test length(poly_roots(x -> x^5 - x -1, Over.Q)) == 0

    # test deflation
    fn = x -> x^3 * (x-1) * (x^2 + 1)
    ds = sort(poly_roots(fn, Over.R)) .- sort([1,0])
    @test norm(ds) <= 1e-14

    # test for input
    poly_roots(x -> x^4 - 1, Over.Q)
    poly_roots(x -> x^4 - 1//1, Over.Q)
    @test_throws ArgumentError poly_roots(x -> x^4 - 1.0, Over.Q)
end

@testset "Over.Z" begin
    @test length(poly_roots(x -> (x-1) * (2x-1) * (x^2 + 1), Over.Z)) == 1

    # test for input
    poly_roots(x -> x^4 - 1, Over.Z)
    @test_throws MethodError poly_roots(x -> x^4 - 1//1, Over.Z)
    @test_throws MethodError poly_roots(x -> x^4 - 1.0, Over.Z)
end

@testset "Over.Zp{q}" begin
    p = x -> x^8 - 1
    @test length(poly_roots(p, Over.Zp{7})) == 2
    @test length(poly_roots(p, Over.Zp{17})) == 8
end

@testset "special cases" begin
    x = variable()
    p = (x-1)*(x^2 + 1)
    rts_c = poly_roots(p, Over.C)
    rts_r = poly_roots(p, Over.R)
    @test length(rts_c) == 3
    @test length(rts_r) == 1

    x = variable(Int)
    p = (x-1)*(2x-3)
    rts_q = poly_roots(p, Over.Q)
    rts_z = poly_roots(p, Over.Z)
    @test length(rts_q) == 2
    @test length(rts_z) == 1
end

@testset "Different Types" begin

    FTs = [Float16, Float32, Float64, BigFloat]
    fn = x -> x^5 - x - 1
    fn1 = x -> (x-1)*(x-2)*(x^5 - x - 1)

    # over.C
    ## all methods promote, this just checks for errors
    for T in [Float64, BigFloat]
       poly_roots(fn, Over.CC{T}, method=:PolynomialRoots)
    end
    for T in [Float64]
       poly_roots(fn, Over.CC{T}, method=:roots)
    end
    for T in FTs
        poly_roots(fn, Over.CC{T}, method=:amvw)
    end
    @test_throws MethodError poly_roots(fn, Over.CC{Int})

    # over.R
    [@test eltype(poly_roots(fn, Over.RR{T})) == T for T in FTs]
    @test_throws MethodError poly_roots(fn, Over.RR{Int})

    # over.Z
    [poly_roots(fn1, Over.ZZ{T})  for T in [Int16, Int32, Int64, Int128]]
    @test eltype(poly_roots(fn1, Over.ZZ{Int32})) == Int32
    @test_throws MethodError poly_roots(fn1, Over.ZZ{Float64})

    p = poly([3.0])^5
    @test_throws MethodError poly_roots(p, Over.Z) # p has Float64 coefficients.
    @test poly_roots(convert(Poly{Int}, p), Over.Z) == [-3.0]

end


text:
"""Reference for parameters update when present"""
struct UpdateRef{T}
    access_ref::Union{Symbol,String}
end

const DRDA = Dict{UpdateRef, JuMP.Containers.DenseAxisArray}


text:
export @nvrtc

type CuFunction
    m::CuModule # avoid CuModule gc-ed
    ptr::Ptr{Void}
end

function CuFunction(m::CuModule, name::String)
    p = CUfunction[0]
    cuModuleGetFunction(p, m, name)
    CuFunction(m, p[1])
end

Base.unsafe_convert(::Type{Ptr{Void}}, f::CuFunction) = f.ptr

box(x) = pointer_from_objref(x)

macro nvrtc(expr)
    expr.head == :string || throw("expr is not string")
    idx = findfirst(expr.args) do a
        typeof(a) <: String && match(r"__global__ void", a) != nothing
    end
    idx == 0 && throw("Cannot find \"__global__ void\".")
    dict = Dict()
    for i = idx+1:length(expr.args)
        arg = expr.args[i]
        typeof(arg) == Symbol && (dict[arg] = arg)
        typeof(arg) == String && contains(arg, ")") && break
    end
    syms = Expr(:tuple, keys(dict)...)
    dict = Dict()

    quote
        local dict = $dict
        local key = $(esc(syms))
        if haskey(dict, key)
            dict[key]
        else
            local code = $(esc(expr))
            local ptx = NVRTC.compile(code)
            f = load_ptx(ptx)
            dict[key] = f
            f
        end
    end
end

function load_ptx(ptx::String)
    p = Ptr{Void}[0]
    cuModuleLoadData(p, pointer(ptx))
    mod = CuModule(p[1])
    # TODO: multi-device
    for line in split(ptx,'\n')
        m = match(r".visible .entry (.+)\(", line) # find function name
        m == nothing && continue
        fname = Symbol(m[1])
        return CuFunction(mod, string(fname))
    end
end

function (f::CuFunction)(args...;
    dx=1, dy=1, dz=1, bx=128, by=1, bz=1, sharedmem=0, stream=C_NULL)

    argptrs = Ptr{Void}[box(a) for a in args]
    gx = ceil(dx / bx)
    gy = ceil(dy / by)
    gz = ceil(dz / bz)
    cuLaunchKernel(f, gx, gy, gz, bx, by, bz, sharedmem, stream, argptrs, C_NULL)
end


text:
#=
How to use:
```
$ julia script.jl a b c
a
b
c
```
=#

for arg in ARGS
    println(arg)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
@testset "acquisition functions" begin
    for ac in [ProbabilityOfImprovement(), ExpectedImprovement(),
               UpperConfidenceBound(), ThompsonSamplingSimple(), MutualInformation()]
        model = GPE(rand(3, 4), rand(4), MeanZero(), SEIso(0., 0.))
        acfunc = BayesianOptimization.acquisitionfunction(ac, model)
        x = rand(3, 2)
        acvector = acfunc(x)
        @test length(acvector) == 2
        if typeof(ac) != ThompsonSamplingSimple
            @test acvector[1] == acfunc(x[:, 1])
        end
    end
end


text:
function gencond( kappa, m, n )
  A = randn(m,n)
  (U,S,V) = svd(A)
  S = logspace(0, -log10(kappa), n)
  return U*diagm(S)*V'
end


text:
# # Langmuir turbulence example
#
# This example implements the Langmuir turbulence simulation reported in section
# 4 of
#
# > [McWilliams, J. C. et al., "Langmuir Turbulence in the ocean," Journal of Fluid Mechanics (1997)](https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/langmuir-turbulence-in-the-ocean/638FD0E368140E5972144348DB930A38).
#
# This example demonstrates
#
#   * How to run large eddy simulations with surface wave effects via the
#     Craik-Leibovich approximation.
#
#   * How to specify time- and horizontally-averaged output.

# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, Plots"
# ```

using Oceananigans
using Oceananigans.Units: minute, minutes, hours

# ## Model set-up
#
# To build the model, we specify the grid, Stokes drift, boundary conditions, and
# Coriolis parameter.
#
# ### Domain and numerical grid specification
#
# We create a grid with modest resolution. The grid extent is similar, but not
# exactly the same as that in McWilliams et al. (1997).

grid = RectilinearGrid(size=(32, 32, 48), extent=(128, 128, 96))

# ### The Stokes Drift profile
#
# The surface wave Stokes drift profile prescribed in McWilliams et al. (1997)
# corresponds to a 'monochromatic' (that is, single-frequency) wave field.
#
# A monochromatic wave field is characterized by its wavelength and amplitude
# (half the distance from wave crest to wave trough), which determine the wave
# frequency and the vertical scale of the Stokes drift profile.

using Oceananigans.BuoyancyModels: g_Earth

 amplitude = 0.8 # m
wavelength = 60 # m
wavenumber = 2π / wavelength # m⁻¹
 frequency = sqrt(g_Earth * wavenumber) # s⁻¹

## The vertical scale over which the Stokes drift of a monochromatic surface wave
## decays away from the surface is `1/2wavenumber`, or
const vertical_scale = wavelength / 4π

## Stokes drift velocity at the surface
const Uˢ = amplitude^2 * wavenumber * frequency # m s⁻¹

# The `const` declarations ensure that Stokes drift functions compile on the GPU.
# To run this example on the GPU, write `architecture = GPU()` in the constructor
# for `NonhydrostaticModel` below.
#
# The Stokes drift profile is

uˢ(z) = Uˢ * exp(z / vertical_scale)

# which we'll need for the initial condition.
#
# !!! info "The Craik-Leibovich equations in Oceananigans"
#     Oceananigans implements the Craik-Leibovich approximation for surface wave effects
#     using the _Lagrangian-mean_ velocity field as its prognostic momentum variable.
#     In other words, `model.velocities.u` is the Lagrangian-mean ``x``-velocity beneath surface
#     waves. This differs from models that use the _Eulerian-mean_ velocity field
#     as a prognostic variable, but has the advantage that ``u`` accounts for the total advection
#     of tracers and momentum, and that ``u = v = w = 0`` is a steady solution even when Coriolis
#     forces are present. See the
#     [physics documentation](https://clima.github.io/OceananigansDocumentation/stable/physics/surface_gravity_waves/)
#     for more information.
#
# The vertical derivative of the Stokes drift is

∂z_uˢ(z, t) = 1 / vertical_scale * Uˢ * exp(z / vertical_scale)

# Finally, we note that the time-derivative of the Stokes drift must be provided
# if the Stokes drift changes in time. In this example, the Stokes drift is constant
# and thus the time-derivative of the Stokes drift is 0.

# ### Boundary conditions
#
# At the surface at ``z=0``, McWilliams et al. (1997) impose a wind stress
# on ``u``,

Qᵘ = -3.72e-5 # m² s⁻², surface kinematic momentum flux

u_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

# McWilliams et al. (1997) impose a linear buoyancy gradient `N²` at the bottom
# along with a weak, destabilizing flux of buoyancy at the surface to faciliate
# spin-up from rest.

Qᵇ = 2.307e-9 # m³ s⁻², surface buoyancy flux
N² = 1.936e-5 # s⁻², initial and bottom buoyancy gradient

b_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ),
                                                bottom = GradientBoundaryCondition(N²))

# !!! info "The flux convention in Oceananigans"
#     Note that Oceananigans uses "positive upward" conventions for all fluxes. In consequence,
#     a negative flux at the surface drives positive velocities, and a positive flux of
#     buoyancy drives cooling.

# ### Coriolis parameter
#
# McWilliams et al. (1997) use

coriolis = FPlane(f=1e-4) # s⁻¹

# which is typical for mid-latitudes on Earth.

# ## Model instantiation
#
# We are ready to build the model. We use a fifth-order Weighted Essentially
# Non-Oscillatory (WENO) advection scheme and the `AnisotropicMinimumDissipation`
# model for large eddy simulation. Because our Stokes drift does not vary in ``x, y``,
# we use `UniformStokesDrift`, which expects Stokes drift functions of ``z, t`` only.

model = NonhydrostaticModel(
              advection = WENO5(),
            timestepper = :RungeKutta3,
                   grid = grid,
                tracers = :b,
               buoyancy = BuoyancyTracer(),
               coriolis = coriolis,
                closure = AnisotropicMinimumDissipation(),
           stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
    boundary_conditions = (u=u_boundary_conditions, b=b_boundary_conditions),
)

# ## Initial conditions
#
# We make use of random noise concentrated in the upper 4 meters
# for buoyancy and velocity initial conditions,

Ξ(z) = randn() * exp(z / 4)
nothing # hide

# Our initial condition for buoyancy consists of a surface mixed layer 33 m deep,
# a deep linear stratification, plus noise,

initial_mixed_layer_depth = 33 # m
stratification(z) = z < - initial_mixed_layer_depth ? N² * z : N² * (-initial_mixed_layer_depth)

bᵢ(x, y, z) = stratification(z) + 1e-1 * Ξ(z) * N² * model.grid.Lz

# The velocity initial condition in McWilliams et al. (1997) is zero *Eulerian-mean* velocity.
# This means that we must add the Stokes drift profile to the Lagrangian-mean ``u`` velocity field
# modeled by Oceananigans.jl. We also add noise scaled by the friction velocity to ``u`` and ``w``.

uᵢ(x, y, z) = uˢ(z) + sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)

wᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)

set!(model, u=uᵢ, w=wᵢ, b=bᵢ)

# ## Setting up the simulation

simulation = Simulation(model, Δt=45.0, stop_time=4hours)

# We use the `TimeStepWizard` for adaptive time-stepping
# with a Courant-Freidrichs-Lewy (CFL) number of 1.0,

wizard = TimeStepWizard(cfl=1.0, max_change=1.1, max_Δt=1minute)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

# ### Nice progress messaging
#
# We define a function that prints a helpful message with
# maximum absolute value of ``u, v, w`` and the current wall clock time.

using Printf

function print_progress(simulation)
    u, v, w = simulation.model.velocities

    ## Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
                   iteration(simulation),
                   prettytime(time(simulation)),
                   prettytime(simulation.Δt),
                   maximum(abs, u), maximum(abs, v), maximum(abs, w),
                   prettytime(simulation.run_wall_time))

    @info msg

    return nothing
end

simulation.callbacks[:progress] = Callback(print_progress, IterationInterval(10))

# ## Output
#
# ### A field writer
#
# We set up an output writer for the simulation that saves all velocity fields,
# tracer fields, and the subgrid turbulent diffusivity.

output_interval = 5minutes

fields_to_output = merge(model.velocities, model.tracers, (νₑ=model.diffusivity_fields.νₑ,))

simulation.output_writers[:fields] =
    JLD2OutputWriter(model, fields_to_output,
                     schedule = TimeInterval(output_interval),
                     prefix = "langmuir_turbulence_fields",
                     force = true)

# ### An "averages" writer
#
# We also set up output of time- and horizontally-averaged velocity field and
# momentum fluxes,

u, v, w = model.velocities

U = AveragedField(u, dims=(1, 2))
V = AveragedField(v, dims=(1, 2))
B = AveragedField(model.tracers.b, dims=(1, 2))

wu = AveragedField(w * u, dims=(1, 2))
wv = AveragedField(w * v, dims=(1, 2))

simulation.output_writers[:averages] =
    JLD2OutputWriter(model, (u=U, v=V, b=B, wu=wu, wv=wv),
                     schedule = AveragedTimeInterval(output_interval, window=2minutes),
                     prefix = "langmuir_turbulence_averages",
                     force = true)

# ## Running the simulation
#
# This part is easy,

run!(simulation)

# # Making a neat movie
#
# We look at the results by plotting vertical slices of ``u`` and ``w``, and a horizontal
# slice of ``w`` to look for Langmuir cells.

k = searchsortedfirst(grid.zᵃᵃᶠ[:], -8)
nothing # hide

# Making the coordinate arrays takes a few lines of code,

xw, yw, zw = nodes(model.velocities.w)
xu, yu, zu = nodes(model.velocities.u)
nothing # hide

# Next, we open the JLD2 file, and extract the iterations we ended up saving at,

using JLD2, Plots

fields_file = jldopen(simulation.output_writers[:fields].filepath)
averages_file = jldopen(simulation.output_writers[:averages].filepath)

iterations = parse.(Int, keys(fields_file["timeseries/t"]))

# This utility is handy for calculating nice contour intervals:

function nice_divergent_levels(c, clim; nlevels=20)
    levels = range(-clim, stop=clim, length=nlevels)
    cmax = maximum(abs, c)
    clim < cmax && (levels = vcat([-cmax], levels, [cmax]))
    return (-clim, clim), levels
end
nothing # hide

# Finally, we're ready to animate.

@info "Making an animation from the saved data..."

anim = @animate for (i, iter) in enumerate(iterations)

    @info "Drawing frame $i from iteration $iter \n"

    ## Load 3D fields from fields_file
    t = fields_file["timeseries/t/$iter"]
    w_snapshot = fields_file["timeseries/w/$iter"]
    u_snapshot = fields_file["timeseries/u/$iter"]

    B_snapshot = averages_file["timeseries/b/$iter"][1, 1, :]
    U_snapshot = averages_file["timeseries/u/$iter"][1, 1, :]
    V_snapshot = averages_file["timeseries/v/$iter"][1, 1, :]
    wu_snapshot = averages_file["timeseries/wu/$iter"][1, 1, :]
    wv_snapshot = averages_file["timeseries/wv/$iter"][1, 1, :]

    ## Extract slices
    wxy = w_snapshot[:, :, k]
    wxz = w_snapshot[:, 1, :]
    uxz = u_snapshot[:, 1, :]

    wlims, wlevels = nice_divergent_levels(w, 0.03)
    ulims, ulevels = nice_divergent_levels(w, 0.05)

    B_plot = plot(B_snapshot, zu,
                  label = nothing,
                  legend = :bottom,
                  xlabel = "Buoyancy (m s⁻²)",
                  ylabel = "z (m)")

    U_plot = plot([U_snapshot V_snapshot], zu,
                  label = ["\$ \\bar u \$" "\$ \\bar v \$"],
                  legend = :bottom,
                  xlabel = "Velocities (m s⁻¹)",
                  ylabel = "z (m)")

    wu_label = "\$ \\overline{wu} \$"
    wv_label = "\$ \\overline{wv} \$"

    fluxes_plot = plot([wu_snapshot, wv_snapshot], zw,
                       label = [wu_label wv_label],
                       legend = :bottom,
                       xlabel = "Momentum fluxes (m² s⁻²)",
                       ylabel = "z (m)")

    wxy_plot = contourf(xw, yw, wxy';
                        color = :balance,
                        linewidth = 0,
                        aspectratio = :equal,
                        clims = wlims,
                        levels = wlevels,
                        xlims = (0, grid.Lx),
                        ylims = (0, grid.Ly),
                        xlabel = "x (m)",
                        ylabel = "y (m)")

    wxz_plot = contourf(xw, zw, wxz';
                              color = :balance,
                          linewidth = 0,
                        aspectratio = :equal,
                              clims = wlims,
                             levels = wlevels,
                              xlims = (0, grid.Lx),
                              ylims = (-grid.Lz, 0),
                             xlabel = "x (m)",
                             ylabel = "z (m)")

    uxz_plot = contourf(xu, zu, uxz';
                              color = :balance,
                          linewidth = 0,
                        aspectratio = :equal,
                              clims = ulims,
                             levels = ulevels,
                              xlims = (0, grid.Lx),
                              ylims = (-grid.Lz, 0),
                             xlabel = "x (m)",
                             ylabel = "z (m)")

    wxy_title = @sprintf("w(x, y, t) (m s⁻¹) at z=-8 m and t = %s ", prettytime(t))
    wxz_title = @sprintf("w(x, z, t) (m s⁻¹) at y=0 m and t = %s", prettytime(t))
    uxz_title = @sprintf("u(x, z, t) (m s⁻¹) at y=0 m and t = %s", prettytime(t))

    plot(wxy_plot, B_plot, wxz_plot, U_plot, uxz_plot, fluxes_plot,
         layout = Plots.grid(3, 2, widths=(0.7, 0.3)), size = (900.5, 1000.5),
         title = [wxy_title "" wxz_title "" uxz_title ""])

    if iter == iterations[end]
        close(fields_file)
        close(averages_file)
    end
end

mp4(anim, "langmuir_turbulence.mp4", fps = 8) # hide


text:
Base.@kwdef mutable struct Args
	learning_rate::Float64 = 1e-2
    batchsize::Int = 50
    epochs::Int = 1000
	sd_λ::Float64 = 3e-3
    wd_λ::Float64 = 1e-2
end

struct Experiment
    optimizer::String
    offset::Float64
    sd::Bool
end


text:
using ThreadingExperiments, .Threads
using Test

@testset "ThreadingExperiments.jl" begin
    for i in 1:nthreads()
        @test onthread(threadid, i) == i
    end
    @test_throws AssertionError onthread(threadid, 0815)
end


text:
#Compute Threshold Functions

module ThresholdFunctions

using ..MathFuns

@doc raw"""
    Fb1(k,m,T)
Compute the one bosonic propagator threshold function without external momentum dependence.
```math
\tilde{\partial_k}\int \frac{d^3q}{2\pi^3}\sum_{q_0}G_b(q_0)
```
"""
function Fb1(k, m2, T)
    -1 / 24 * (
        k^4 * (
            2 * T * coth(Eb(k, m2) / (2 * T)) +
            Eb(k, m2) * csch(Eb(k, m2) / (2 * T))^2
        )
    ) / (Eb(k, m2)^3 * pi^2 * T)
end

@doc raw"""
    Fb2(k,m,T)
Compute the two bosonic propagator threshold function without external momentum dependence.
```math
\tilde{\partial_k}\int \frac{d^3q}{2\pi^3}\sum_{q_0}G_b(q_0)G_b(q_0)
```
If we want to start from the momentum dependent case and take the zero momentum limit, we need to
take the $p_0\to 0$ first then take the $p_s\to 0$, they are noncommutative, or we need to take the $p_0$
goes faster to $0$ than $p_s$
"""
function Fb2(k,m2,T)
    (
        k^4 * (
            -6 * T^2 * coth(Eb(k, m2) / (2 * T)) -
            csch(Eb(k, m2) / (2 * T))^2 *
            Eb(k, m2) *
            (3 * T + coth(Eb(k, m2) / (2 * T)) * Eb(k, m2))
        )
    ) / (48 * pi^2 * T^2 * Eb(k, m2)^5)
end

function ReFb2positive(p0, k, m2, T,ϵ=0.1)
    -1 / 24 * (
        k^4 *
        csch(Eb(k, m2) / (2 * T))^2 *
        (p0 - 2 * Eb(k, m2))^2 *
        sinh(Eb(k, m2) / T)
    ) / (pi^2 * (ϵ + (p0 - 2 * Eb(k, m2))^2)^2 * Eb(k, m2)^3) -
    (k^4 * csch(Eb(k, m2) / (2 * T))^2 * sinh(Eb(k, m2) / T)) /
    (24 * pi^2 * Eb(k, m2)^3 * (p0 + 2 * Eb(k, m2))^2) +
    (
        k^4 *
        csch(Eb(k, m2) / (2 * T))^2 *
        (p0 - 2 * Eb(k, m2)) *
        (Eb(k, m2) + 2 * T * sinh(Eb(k, m2) / T))
    ) / (48 * pi^2 * T * (ϵ + (p0 - 2 * Eb(k, m2))^2) * Eb(k, m2)^4) -
    (
        k^4 *
        csch(Eb(k, m2) / (2 * T))^2 *
        (Eb(k, m2) + 2 * T * sinh(Eb(k, m2) / T))
    ) / (48 * pi^2 * T * Eb(k, m2)^4 * (p0 + 2 * Eb(k, m2)))
end

#we add the ϵ by analitical continuation
function ReFb2positive_ac(ω, k, m2, T, ϵ = 0.1)
    -1 / 6 * (
        k^4 * (
            (ϵ^2 + (ω - 2 * Eb(k, m2))^2) *
            Eb(k, m2) *
            (ϵ^2 - ω^2 + 4 * Eb(k, m2)^2) *
            (ϵ^2 + (ω + 2 * Eb(k, m2))^2) +
            T *
            (
                (ϵ - ω) * (ϵ + ω) * (ϵ^2 + ω^2)^2 +
                4 * (5 * ϵ^4 - 14 * ϵ^2 * ω^2 + 5 * ω^4) * Eb(k, m2)^2 +
                112 * (ϵ - ω) * (ϵ + ω) * Eb(k, m2)^4 +
                192 * Eb(k, m2)^6
            ) *
            sinh(Eb(k, m2) / T)
        )
    ) / (
        pi^2 *
        T *
        (-1 + cosh(Eb(k, m2) / T)) *
        Eb(k, m2)^3 *
        (
            (ϵ^2 + ω^2)^2 +
            8 * (ϵ - ω) * (ϵ + ω) * Eb(k, m2)^2 +
            16 * Eb(k, m2)^4
        )^2
    )
end

function ReFb2positiveEu(p0, k, m2, T)
    -1 / 12 * (
        k^4 *
        csch(Eb(k, m2) / (2 * T))^2 *
        (
            Eb(k, m2) * (p0^2 + 4 * Eb(k, m2)^2) +
            T * (p0^2 + 12 * Eb(k, m2)^2) * sinh(Eb(k, m2) / T)
        )
    ) / (pi^2 * T * Eb(k, m2)^3 * (p0^2 + 4 * Eb(k, m2)^2)^2)
end

function ImFb2positive(p0, k, m2, T, ϵ = 0.1)
    -1 / 12 *
    (k^4 * coth(Eb(k, m2) / (2 * T)) * deltadxfun(p0 - 2 * Eb(k, m2), ϵ)) /
    (pi * Eb(k, m2)^3) -
    (
        k^4 *
        csch(Eb(k, m2) / (2 * T))^2 *
        deltafun(p0 - 2 * Eb(k, m2), ϵ) *
        (Eb(k, m2) + 2 * T * sinh(Eb(k, m2) / T))
    ) / (48 * pi * T * Eb(k, m2)^4)
end

function ImFb2(p0, k, m2, T, ϵ = 0.1)
    if p0>0.0
        ImFb2positive(p0, k, m2, T, ϵ)
    elseif p0<0.0
        -ImFb2positive(-p0, k, m2, T, ϵ)
    elseif p0==0.0
        0.0
    end
end

function ReFb2(p0, k, m2, T, ϵ = 0.1)
    if p0>=0.0
        ReFb2positive(p0, k, m2, T, ϵ)
    elseif p0<0.0
        ReFb2positive(-p0, k, m2, T, ϵ)
    end
end

function ReFb2_ac(p0, k, m2, T, ϵ = 0.1)
    if p0>=0.0
        ReFb2positive_ac(p0, k, m2, T, ϵ)
    elseif p0<0.0
        ReFb2positive_ac(-p0, k, m2, T, ϵ)
    end
end

function ReFb2Eu(p0, k, m2, T)
    if p0>=0.0
        ReFb2positiveEu(p0, k, m2, T)
    elseif p0<0.0
        ReFb2positiveEu(-p0, k, m2, T)
    end
end

end


text:
mutable struct Point <: GeoInterface.AbstractPoint
    ptr::GEOSGeom

    function Point(ptr::GEOSGeom)
        point = new(ptr)
        finalizer(point, destroyGeom)
        point
    end
    Point(coords::Vector{Float64}) = Point(createPoint(coords))
    Point(x::Real, y::Real) = Point(createPoint(x,y))
    Point(x::Real, y::Real, z::Real) = Point(createPoint(x,y,z))
    Point(obj::T) where {T<:GeoInterface.AbstractPoint} = Point(GeoInterface.coordinates(obj))
end

mutable struct MultiPoint <: GeoInterface.AbstractMultiPoint
    ptr::GEOSGeom

    function MultiPoint(ptr::GEOSGeom)
        multipoint = new(ptr)
        finalizer(multipoint, destroyGeom)
        multipoint
    end
    MultiPoint(multipoint::Vector{Vector{Float64}}) = MultiPoint(createCollection(GEOS_MULTIPOINT, GEOSGeom[createPoint(coords) for coords in multipoint]))
    MultiPoint(obj::T) where {T<:GeoInterface.AbstractMultiPoint} = MultiPoint(GeoInterface.coordinates(obj))
end

mutable struct LineString <: GeoInterface.AbstractLineString
    ptr::GEOSGeom

    function LineString(ptr::GEOSGeom)
        line = new(ptr)
        finalizer(line, destroyGeom)
        line
    end
    LineString(line::Vector{Vector{Float64}}) = LineString(createLineString(line))
    LineString(obj::T) where {T<:GeoInterface.AbstractLineString} = LineString(GeoInterface.coordinates(obj))
end

mutable struct MultiLineString <: GeoInterface.AbstractMultiLineString
    ptr::GEOSGeom

    function MultiLineString(ptr::GEOSGeom)
        multiline = new(ptr)
        finalizer(multiline, destroyGeom)
        multiline
    end
    MultiLineString(multiline::Vector{Vector{Vector{Float64}}}) = MultiLineString(createCollection(GEOS_MULTILINESTRING, GEOSGeom[createLineString(coords) for coords in multiline]))
    MultiLineString(obj::T) where {T<:GeoInterface.AbstractMultiLineString} = MultiLineString(GeoInterface.coordinates(obj))
end

mutable struct LinearRing <: GeoInterface.AbstractLineString
    ptr::GEOSGeom

    function LinearRing(ptr::GEOSGeom)
        ring = new(ptr)
        finalizer(ring, destroyGeom)
        ring
    end
    LinearRing(ring::Vector{Vector{Float64}}) = LinearRing(createLinearRing(ring))
    LinearRing(obj::T) where {T<:GeoInterface.AbstractLineString} = LinearRing(GeoInterface.coordinates(obj))
end

mutable struct Polygon <: GeoInterface.AbstractPolygon
    ptr::GEOSGeom

    function Polygon(ptr::GEOSGeom)
        polygon = new(ptr)
        finalizer(polygon, destroyGeom)
        polygon
    end
    function Polygon(coords::Vector{Vector{Vector{Float64}}})
        exterior = createLinearRing(coords[1])
        interiors = GEOSGeom[createLinearRing(lr) for lr in coords[2:end]]
        polygon = new(createPolygon(exterior,interiors))
        finalizer(polygon, destroyGeom)
        polygon
    end
    Polygon(obj::T) where {T<:GeoInterface.AbstractPolygon} = Polygon(GeoInterface.coordinates(obj))
end

mutable struct MultiPolygon <: GeoInterface.AbstractMultiPolygon
    ptr::GEOSGeom

    function MultiPolygon(ptr::GEOSGeom)
        multipolygon = new(ptr)
        finalizer(multipolygon, destroyGeom)
        multipolygon
    end
    MultiPolygon(multipolygon::Vector{Vector{Vector{Vector{Float64}}}}) =
        MultiPolygon(createCollection(GEOS_MULTIPOLYGON,
                                      GEOSGeom[createPolygon(createLinearRing(coords[1]),
                                                             GEOSGeom[createLinearRing(c) for c in coords[2:end]])
                                               for coords in multipolygon]))
    MultiPolygon(obj::T) where {T<:GeoInterface.AbstractMultiPolygon} = MultiPolygon(GeoInterface.coordinates(obj))
end

mutable struct GeometryCollection <: GeoInterface.AbstractGeometryCollection
    ptr::GEOSGeom

    function GeometryCollection(ptr::GEOSGeom)
        geometrycollection = new(ptr)
        finalizer(geometrycollection, destroyGeom)
        geometrycollection
    end
    GeometryCollection(collection::Vector{GEOSGeom}) = GeometryCollection(createCollection(GEOS_GEOMETRYCOLLECTION, collection))
end

for geom in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon, :GeometryCollection)
    @eval begin
        function destroyGeom(obj::$geom)
            destroyGeom(obj.ptr)
            obj.ptr = C_NULL
        end
    end
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module PGFPlots

export LaTeXString, @L_str, @L_mstr
export plot, Axis, Axes, PolarAxis, GroupPlot, Plots, ColorMaps, save, define_color
export pushPGFPlotsOptions, popPGFPlotsOptions, resetPGFPlotsOptions, pgfplotsoptions
export pushPGFPlotsPreamble, popPGFPlotsPreamble, resetPGFPlotsPreamble, pgfplotspreamble
export pushPGFPlots, popPGFPlots
import Colors: RGB
import Contour: contours, levels

using Compat
using Discretizers

include("colormaps.jl")
include("plots.jl")

import TikzPictures: TikzPicture, PDF, TEX, SVG, save, LaTeXString, @L_str, @L_mstr

_pgfplotsoptions = Any[""]

pushPGFPlotsOptions(options::AbstractString) = push!(_pgfplotsoptions, options)
function popPGFPlotsOptions()
    if length(_pgfplotsoptions) > 1
        pop!(_pgfplotsoptions)
    end
end
resetPGFPlotsOptions() = deleteat!(_pgfplotsoptions, 2:length(_pgfplotsoptions))

function pushPGFPlots(options::AbstractString)
    warn("Use pushPGFPlotsOptions")
    pushPGFPlotsOptions(options)
end

function popPGFPlots()
    warn("Use popPGFPlotsOptions")
    popPGFPlotsOptions()
end

pgfplotsoptions() = join(_pgfplotsoptions, ",\n")

_pgfplotspreamble = Any[readstring(joinpath(dirname(@__FILE__), "preamble.tex"))]

pushPGFPlotsPreamble(preamble::AbstractString) = push!(_pgfplotspreamble, preamble)
function popPGFPlotsPreamble()
    if length(_pgfplotspreamble) > 1
        pop!(_pgfplotspreamble)
    end
end
resetPGFPlotsPreamble() = deleteat!(_pgfplotspreamble,2:length(_pgfplotspreamble))

pgfplotspreamble() = join(_pgfplotspreamble, "\n")

function define_color(name::AbstractString, color::RGB)
    _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{rgb}{$(color.r), $(color.g), $(color.b)}"
end

define_color(name::AbstractString, color) = define_color(name, convert(RGB, color))

function define_color{T<:AbstractFloat}(name::AbstractString, color::Vector{T})
    if length(color) == 3
        r, g, b = color[1], color[2], color[3]
        _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{rgb}{$r, $g, $b}"
    elseif length(color) == 4
        c, m, y, k = color[1], color[2], color[3], color[4]
        _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{cmyk}{$c, $m, $y, $k}"
    else
        error("Color must have three or four components")
    end
end

function define_color{T<:Integer}(name::AbstractString, color::Vector{T})
    @assert(length(color) == 3, "Color must have three components")
    r, g, b = color[1], color[2], color[3]
    _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{RGB}{$r, $g, $b}"
end

function define_color(name::AbstractString, color::UInt32)
    _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{HTML}{$(uppercase(hex(color)))}"
end

function define_color(name::AbstractString, color::Real)
    _pgfplotspreamble[end] = _pgfplotspreamble[end] * "\n\\definecolor{$name}{gray}{$color}"
end

histogramMap = Dict(
    :bins => "bins",
    :density => "density",
    :cumulative => "cumulative"
    )

linearMap = Dict(
    :mark => "mark",
    :markSize => "mark size",
    :style => "",
    :onlyMarks => "only marks"
    )

scatterMap = Dict(
    :mark => "mark",
    :markSize => "mark size",
    :style => "",
    :onlyMarks => "only marks",
    :scatterClasses => "scatter/classes"
    )

errorbarsMap = Dict(
    :mark => "mark",
    :style => ""
    )

quiverMap = Dict(
    :style => ""
    )

contourMap = Dict(
    :number => "number",
    :levels => "levels",
    :style => "",
    :labels => "labels"
    )

using .Plots
using .ColorMaps

typealias IntegerRange @compat Tuple{Integer,Integer}
typealias RealRange @compat Tuple{Real,Real}

type Axis
    plots::Vector{Plot}
    title
    xlabel
    xlabelStyle
    ylabel
    ylabelStyle
    zlabel
    zlabelStyle
    xmin
    xmax
    ymin
    ymax
    axisEqual
    axisEqualImage
    enlargelimits
    axisOnTop
    view
    width
    height
    style
    legendPos
    legendStyle
    xmode
    ymode
    colorbar
    hideAxis
    axisLines
    axisKeyword

    Axis{P <: Plot}(plots::Vector{P};title=nothing, xlabel=nothing, xlabelStyle=nothing, ylabel=nothing, ylabelStyle=nothing, zlabel=nothing, zlabelStyle=nothing, xmin=nothing, xmax=nothing,
                    ymin=nothing, ymax=nothing, axisEqual=nothing, axisEqualImage=nothing, enlargelimits=nothing, axisOnTop=nothing, view=nothing, width=nothing,
                    height=nothing, style=nothing, legendPos=nothing, legendStyle=nothing, xmode=nothing, ymode=nothing, colorbar=nothing, hideAxis=nothing, axisLines=nothing, axisKeyword="axis") =
        new(plots, title, xlabel, xlabelStyle, ylabel, ylabelStyle, zlabel, zlabelStyle, xmin, xmax, ymin, ymax, axisEqual, axisEqualImage, enlargelimits, axisOnTop, view, width, height, style, legendPos, legendStyle, xmode, ymode, colorbar, hideAxis, axisLines, axisKeyword
            )

    Axis(;kwargs...) = Axis(Plot[]; kwargs...)
    Axis(plot::Plot; kwargs...) = Axis(Plot[plot]; kwargs...)

    Axis(plot::Contour; kwargs...) = Axis(Plot[plot]; kwargs..., view="{0}{90}")
    Axis(plots::Vector{Contour}; kwargs...) = Axis(Plot[plots...]; kwargs..., view="{0}{90}")
end

PolarAxis(args...; kwargs...) = Axis(args...; kwargs..., axisKeyword = "polaraxis")

typealias Axes Vector{Axis}

function Base.push!(g::Axis, p::Plot)
    push!(g.plots, p)
end

axisMap = Dict(
    :title => "title",
    :xlabel => "xlabel",
    :xlabelStyle => "xlabel style",
    :ylabel => "ylabel",
    :ylabelStyle => "ylabel style",
    :zlabel => "zlabel",
    :zlabelStyle => "zlabel style",
    :xmin => "xmin",
    :xmax => "xmax",
    :ymin => "ymin",
    :ymax => "ymax",
    :axisEqual => "axis equal",
    :axisEqualImage => "axis equal image",
    :enlargelimits => "enlargelimits",
    :axisOnTop => "axis on top",
    :view => "view",
    :width => "width",
    :height => "height",
    :style => "",
    :legendPos => "legend pos",
    :legendStyle => "legend style",
    :xmode => "xmode",
    :ymode => "ymode",
    :colorbar => "colorbar",
    :hideAxis => "hide axis",
    :axisLines => "axis lines"
    )

type GroupPlot
    axes::AbstractArray{Axis,1}
    dimensions::IntegerRange
    style
    groupStyle
    GroupPlot(;style=nothing,groupStyle=nothing) = new(AbstractArray{Axis,1}[], (0,0), style, groupStyle)
    GroupPlot(axes::AbstractArray{Axis,1};style=nothing,groupStyle=nothing) = new(axes, (length(axes), 1), style, groupStyle)
    GroupPlot(plots;style=nothing,groupStyle=nothing) = new([Axis(p) for p in plots], (length(plots), 1), style, groupStyle)
    GroupPlot(dimensions::IntegerRange;style=nothing,groupStyle=nothing) = new(AbstractArray{Axis,1}[], dimensions, style, groupStyle)
    GroupPlot(rows::Integer, columns::Integer;style=nothing,groupStyle=nothing) = new(AbstractArray{Axis,1}[], (rows, columns), style, groupStyle)
end

function Base.push!(g::GroupPlot, a::Axis)
    push!(g.axes, a)
    if length(g.axes) > prod(g.dimensions)
        g.dimensions = (length(g.axes), 1)
    end
    g
end

function Base.push!(g::GroupPlot, p::Plot)
    push!(g, Axis(p))
end

function printList{T}(o::IOBuffer, a::AbstractArray{T,1}; brackets=false)
    first = true
    for elem in a
        if first
            first = false
            if brackets
                print(o, "{")
            end
            print(o, "$elem")
        else
            print(o, ", $elem")
        end
    end
    if !first && brackets
        print(o, "}")
    end
end

function printObject(o::IOBuffer, object)
    print(o, "$(object)")
end

function printObject{T}(o::IOBuffer, object::AbstractArray{T,1})
    printList(o, object, brackets = true)
end

function optionHelper(o::IOBuffer, m, object; brackets=false, otherOptions=Dict{AbstractString,AbstractString}[], otherText=nothing)
    first = true
    for (sym, str) in m
        if getfield(object,sym) != nothing
            if first
                first = false
                if brackets
                    print(o, "[")
                end
            else
                print(o, ", ")
            end
            if length(str) > 0
                print(o, "$str = {")
            end
            printObject(o, getfield(object,sym))
            if length(str) > 0
                print(o, "}")
            end
        end
    end
    for (k, v) in otherOptions
        if first
            first = false
            if brackets
                print(o, "[")
            end
        else
            print(o, ", ")
        end
        print(o, "$k = $v")
    end
    if otherText != nothing
        for t in otherText
            if t != nothing
                if first
                    first = false
                    if brackets
                        print(o, "[")
                    end
                else
                    print(o, ", ")
                end
                print(o, "$t")
            end
        end
    end
    if !first && brackets
        print(o, "]")
    end
end

function plotHelper(o::IOBuffer, p::Histogram)
    if (p.discretization == :pgfplots && p.bins > 0) ||
       (p.discretization == :default && length(p.data) ≤ Plots.THRESHOLD_NSAMPLES_DISC_OURSELVES)

        print(o, "\\addplot+ [mark=none, $(p.style), hist={")
        optionHelper(o, histogramMap, p)
        print(o, "}] table [row sep=\\\\, y index = 0] {")
        println(o, "data\\\\")
        for d in p.data
            println(o, "$d \\\\ ")
        end
        println(o, "};")

    else
        # discretize using Discretizers.jl

        if p.discretization == :specified && p.bins > 0
            edges = binedges(DiscretizeUniformWidth(p.bins), convert(Vector{Float64}, p.data))
        else
            discretization = p.discretization == :default ? :auto : p.discretization
            edges = binedges(DiscretizeUniformWidth(discretization), convert(Vector{Float64}, p.data))
        end

        linear = Plots._construct_histogram_linear_data(p.data, edges, p.density, p.cumulative)
        if p.style != "fill=blue!10"
            linear.style = p.style
            if !contains(linear.style, "ybar interval")
                linear.style = "ybar interval," * linear.style
            end
        end
        plotHelper(o, linear)
    end
end

plotLegend(o::IOBuffer, entry) = nothing
plotLegend(o::IOBuffer, entry::AbstractString) = println(o, "\\addlegendentry{$entry}")
function plotLegend{T <: AbstractString}(o::IOBuffer, entries::Vector{T})
    for entry in entries
        plotLegend(o, entry)
    end
end

function plotHelper(o::IOBuffer, p::Linear)
    print(o, "\\addplot+ ")
    optionHelper(o, linearMap, p, brackets=true)
    println(o, "coordinates {")
    for i = 1:size(p.data,2)
        println(o, "($(p.data[1,i]), $(p.data[2,i]))")
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IOBuffer, p::Scatter)
    if p.scatterClasses == nothing
        print(o, "\\addplot+[scatter, scatter src=explicit, ")
    else
        print(o, "\\addplot+[scatter, scatter src=explicit symbolic, ")
    end
    optionHelper(o, scatterMap, p)
    println(o, "] coordinates {")
    for i = 1:size(p.data,2)
        println(o, "($(p.data[1,i]), $(p.data[2,i])) [$(p.data[3,i])]")
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

# Specific version for Linear3 type
# Changes are addplot3 (vs addplot) and iterate over all 3 columns
function plotHelper(o::IOBuffer, p::Linear3)
    print(o, "\\addplot3+ ")
    optionHelper(o, linearMap, p, brackets=true)
    println(o, "coordinates {")
    for i = 1:size(p.data,2)
        println(o, "($(p.data[1,i]), $(p.data[2,i]), $(p.data[3,i]))")
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IOBuffer, p::Node)

    axis = p.axis != nothing ? p.axis : "axis cs"

    if p.style != nothing
        println(o, "\\node at ($(axis):$(p.x), $(p.y)) [$(p.style)] {$(p.data)};")
    else
        println(o, "\\node at ($(axis):$(p.x), $(p.y)) {$(p.data)};")
    end
end

function plotHelper(o::IOBuffer, p::ErrorBars)
    print(o, "\\addplot+ [")
    optionHelper(o, errorbarsMap, p)
    print(o, ",error bars/.cd, x dir=both, x explicit, y dir=both, y explicit")
    println(o,"]")
    println(o, "coordinates {")
    for i = 1:size(p.data,2)
        println(o, "($(p.data[1,i]), $(p.data[2,i])) +=($(p.data[3,i]),$(p.data[4,i])) -=($(p.data[5,i]),$(p.data[6,i]))")
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IOBuffer, p::Quiver)
    print(o, "\\addplot+ ")
    optionHelper(o, quiverMap, p, brackets=true, otherOptions=Dict("quiver"=>"{u=\\thisrow{u},v=\\thisrow{v}}"))
    println(o, "table {")
    println(o, "x y u v")
    for i = 1:size(p.data,2)
        println(o, "$(p.data[1,i]) $(p.data[2,i]) $(p.data[3,i]) $(p.data[4,i])")
    end
    println(o, "};")
    plotLegend(o, p.legendentry)
end

function plotHelper(o::IOBuffer, p::Contour)
    arg = 5
    if p.number != nothing
        arg = p.number
    elseif p.levels != nothing
        arg = p.levels
    end
    C = contours(p.xbins, p.ybins, convert(Matrix{Float64}, p.data), arg)
    if p.labels == nothing
        p.labels = true
    end
    if p.labels
        if p.style != nothing
            print(o, "\\addplot3[contour prepared, $(p.style)] table {")
        else
            print(o, "\\addplot3[contour prepared] table {")
        end
    else
        if p.style != nothing
            print(o, "\\addplot3[contour prepared={labels=false}, $(p.style)] table {")
        else
            print(o, "\\addplot3[contour prepared={labels=false}] table {")
        end
    end
    for c in levels(C)
        level = c.level
        for l in c.lines
            for v in l.vertices
                println(o, "$(v[1]) $(v[2]) $level")
            end
            println(o)
        end
    end
    println(o, "};")
end

function plotHelper(o::IOBuffer, p::Circle)
    if p.style != nothing
        println(o, "\\draw[$(p.style)] (axis cs:$(p.xc), $(p.yc)) circle[radius=$(p.radius)];")
    else
        println(o, "\\draw (axis cs:$(p.xc), $(p.yc)) circle[radius=$(p.radius)];")
    end
end

function plotHelper(o::IOBuffer, p::Ellipse)
    if p.style != nothing
        println(o, "\\draw[$(p.style)] (axis cs:$(p.xc), $(p.yc)) ellipse[x radius=$(p.xradius), y radius=$(p.yradius)];")
    else
        println(o, "\\draw (axis cs:$(p.xc), $(p.yc)) ellipse[x radius=$(p.xradius), y radius=$(p.yradius)];")
    end
end

function plotHelper(o::IOBuffer, p::Command)
    println(o, p.cmd*";") #PGFPlots expects commands to be terminated with a ;
end

function plotHelper(o::IOBuffer, p::Image)
    if p.zmin == p.zmax
        error("Your colorbar range limits must not be equal to each other.")
    end
    if p.style != nothing
        println(o, "\\addplot [$(p.style), point meta min=$(p.zmin), point meta max=$(p.zmax)] graphics [xmin=$(p.xmin), xmax=$(p.xmax), ymin=$(p.ymin), ymax=$(p.ymax)] {$(p.filename)};")
    else
        println(o, "\\addplot [point meta min=$(p.zmin), point meta max=$(p.zmax)] graphics [xmin=$(p.xmin), xmax=$(p.xmax), ymin=$(p.ymin), ymax=$(p.ymax)] {$(p.filename)};")
    end
end

# plot option string and contents; no \begin{axis} or \nextgroupplot
function plotHelper(o::IOBuffer, axis::Axis)
    optionHelper(o, axisMap, axis, brackets=true, otherText=[axisOptions(p) for p in axis.plots])
    for p in axis.plots
        plotHelper(o, p)
    end
end

function plot(axis::Axis)
    o = IOBuffer()
    print(o, "\\begin{$(axis.axisKeyword)}")
    plotHelper(o, axis)
    println(o, "\\end{$(axis.axisKeyword)}")
    TikzPicture(takebuf_string(o), options=pgfplotsoptions(), preamble=pgfplotspreamble())
end

function plot(axes::Axes)
    o = IOBuffer()

    for axis in axes
        print(o, "\\begin{$(axis.axisKeyword)}")
        plotHelper(o, axis)
        println(o, "\\end{$(axis.axisKeyword)}")
    end
    TikzPicture(takebuf_string(o), options=pgfplotsoptions(), preamble=pgfplotspreamble())
end

function plot(p::GroupPlot)
    o = IOBuffer()
    style = ""
    if p.style != nothing
        style = p.style * ", "
    end
    groupStyle = ""
    if p.groupStyle != nothing
        groupStyle = p.groupStyle * ", "
    end
    println(o, "\\begin{groupplot}[$(style)group style={$(groupStyle)group size=$(p.dimensions[1]) by $(p.dimensions[2])}]")
    for a in p.axes
        print(o, "\\nextgroupplot ")
        plotHelper(o, a)
    end
    println(o, "\\end{groupplot}")
    mypreamble = pgfplotspreamble() * "\\usepgfplotslibrary{groupplots}"
    TikzPicture(takebuf_string(o), options=pgfplotsoptions(), preamble=mypreamble)
end

typealias Plottable Union{Plot,GroupPlot,Axis,Axes,TikzPicture}

plot(p::Plot) = plot(Axis(p))

function plot(p::Contour)
    plot(Axis(p, xmin=p.xbins[1], xmax=p.xbins[end], ymin=p.ybins[1], ymax=p.ybins[end]))
end

plot{A<:Real,B<:Real}(x::AbstractArray{A,1}, y::AbstractArray{B,1}) = plot(Linear(x, y))

plot{A<:Real,B<:Real,C<:Real}(x::AbstractVector{A}, y::AbstractVector{B}, z::AbstractVector{C}) = plot(Linear3(x, y, z))

function Plots.Linear(f::Function, range::RealRange; mark="none", style=nothing, legendentry=nothing)
    x = linspace(range[1], range[2])
    y = map(f, x)
    Linear(x, y, mark=mark, style=style, legendentry=legendentry)
end

plot(f::Function, range::RealRange) = plot(Linear(f, range))

plot(tkz::TikzPicture) = tkz # tikz pic doesn't need plot, here for convenience

Base.mimewritable(::MIME"image/svg+xml", p::Plottable) = true

cleanup(p::Axis) = map(cleanup, p.plots)
cleanup(axes::Axes) = map(cleanup, axes)
cleanup(p::GroupPlot) = map(cleanup, p.axes)
cleanup(p::Plot) = nothing
cleanup(p::Circle) = nothing
cleanup(p::Ellipse) = nothing
cleanup(p::Command) = nothing
cleanup(p::Image) = rm(p.filename)
cleanup(p::Contour) = nothing

axisOptions(p::Plot) = nothing

function colormapOptions(cm::ColorMaps.Gray)
    if cm.invert
        return "colormap={wb}{gray(0cm)=(1); gray(1cm)=(0)}"
    else
        return "colormap={wb}{gray(0cm)=(0); gray(1cm)=(1)}"
    end
end

function colormapOptions(cm::ColorMaps.RGBArray)
    o = IOBuffer()
    print(o, "colormap={mycolormap}{ ")
    n = length(cm.colors)
    for i = 1:n
        c = cm.colors[i]
        print(o, "rgb($(i-1)cm)=($(c.r),$(c.g),$(c.b)) ")
    end
    print(o, "}")
    takebuf_string(o)
end

function axisOptions(p::Image)
    if p.colorbar
        cmOpt = colormapOptions(p.colormap)
        return "enlargelimits = false, axis on top, $cmOpt, colorbar"
    else
        return "enlargelimits = false, axis on top"
    end
end

@compat function Base.show(f::IO, a::MIME"image/svg+xml", p::Plottable)
    r = Base.show(f, a, plot(p))
    cleanup(p)
    r
end

function save(filename::AbstractString, o::Plottable; include_preamble::Bool=true)
    _, ext = splitext(filename)
    ext = lowercase(ext)
    if ext == ".pdf"
        save(PDF(filename), plot(o))
    elseif ext == ".svg"
        save(SVG(filename), plot(o))
    elseif ext == ".tex"
        save(TEX(filename, include_preamble=include_preamble), plot(o))
    elseif ext == "." || ext == ""
        error("You must specify a file extension.")
    else
        error("Unsupported file extensions: $ext")
    end
end

end # module


text:
export

    ProgressDisplay,
    StoreSnapshots,
    PointMicrophone

# --------------------------------------------------------------------

mutable struct ProgressDisplay
    dt::Float64
    bar::Progress
    showvalues::Bool

    function ProgressDisplay(dt::Number = .1, showvalues = true)
        hook = new()
        hook.dt = Float64(dt)
        hook.showvalues = showvalues
        hook
    end
end

function hook_init!(hook::ProgressDisplay, state::State, sim::Simulator)
    hook.bar = Progress(sim.maxiter, dt = hook.dt)
end

function hook_update!(hook::ProgressDisplay, state::State, sim::Simulator)
    if hook.showvalues
        next!(hook.bar; showvalues = [(:time,state.t), (:iter,state.iter)])
    else
        next!(hook.bar)
    end
end

# --------------------------------------------------------------------

mutable struct StoreSnapshots
    dt::Float64
    last::Float64
    index::Int
    buffer
    snapshots

    function StoreSnapshots(dt::Number)
        hook = new()
        hook.dt = Float64(dt)
        hook
    end
end

function hook_init!(hook::StoreSnapshots, state::State{N,T}, sim::Simulator) where {N,T}
    count = floor(Int, sim.duration / hook.dt)
    hook.last = typemin(Float64)
    hook.index = 1
    hook.buffer = zeros(T, size(state.current))
    hook.snapshots = zeros(T, (size(state.current)..., count))
    hook
end

function hook_update!(hook::StoreSnapshots, state::State{N}, sim::Simulator) where N
    if state.t - (hook.last + hook.dt) >= 0
        if hook.index <= size(hook.snapshots, ndims(hook.snapshots))
            copy!(hook.buffer, state.current)
            copy!(view(hook.snapshots, ntuple(_->:,N)..., hook.index), hook.buffer)
            hook.index += 1
        end
        hook.last = state.t
    end
    hook
end

# --------------------------------------------------------------------

mutable struct PointMicrophone{N}
    position::NTuple{N,Float64}
    startiter::Int
    gridpos::NTuple{N,Int}
    signal

    function PointMicrophone(position::Vararg{Number,N}) where N
        hook = new{N}()
        hook.position = Float64.(position)
        hook.startiter = 1
        hook
    end
end

function hook_init!(hook::PointMicrophone, state::State{N,T}, sim::Simulator) where {N,T}
    hook.signal = zeros(T, sim.maxiter)
    hook.startiter = state.iter
    hook.gridpos = floor.(Int, hook.position ./ sim.wave.dx) .+ 1
    hook
end

function hook_update!(hook::PointMicrophone, state::State{N}, sim::Simulator) where N
    hook.signal[state.iter - hook.startiter] = state.current[hook.gridpos...]
    hook
end

# --------------------------------------------------------------------

mutable struct PointSpeaker{N,S<:AbstractVector}
    position::NTuple{N,Float64}
    gridpos::NTuple{N,Int}
    signal::S

    function PointSpeaker(position::NTuple{N,Number}, signal::S) where {N,S<:AbstractVector}
        hook = new{N,S}()
        hook.position = Float64.(position)
        hook.signal = signal
        hook
    end
end

function hook_init!(hook::PointSpeaker, state::State{N,T}, sim::Simulator) where {N,T}
    hook.gridpos = floor.(Int, hook.position ./ sim.wave.dx) .+ 1
    hook
end

function hook_update!(hook::PointSpeaker, state::State{N}, sim::Simulator) where N
    if state.iter <= length(hook.signal)
        state.current[hook.gridpos...] = hook.signal[state.iter]
    end
    hook
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using MXNet

#--------------------------------------------------------------------------------
# define MLP
# the following two ways are equivalent

#-- Option 1: explicit composition
# data = mx.Variable(:data)
# fc1  = mx.FullyConnected(data, name=:fc1, num_hidden=128)
# act1 = mx.Activation(fc1, name=:relu1, act_type=:relu)
# fc2  = mx.FullyConnected(act1, name=:fc2, num_hidden=64)
# act2 = mx.Activation(fc2, name=:relu2, act_type=:relu)
# fc3  = mx.FullyConnected(act2, name=:fc3, num_hidden=10)
# mlp  = mx.SoftmaxOutput(fc3, name=:softmax)

#-- Option 2: using the mx.chain macro
# mlp = @mx.chain mx.Variable(:data)             =>
#   mx.FullyConnected(name=:fc1, num_hidden=128) =>
#   mx.Activation(name=:relu1, act_type=:relu)   =>
#   mx.FullyConnected(name=:fc2, num_hidden=64)  =>
#   mx.Activation(name=:relu2, act_type=:relu)   =>
#   mx.FullyConnected(name=:fc3, num_hidden=10)  =>
#   mx.SoftmaxOutput(name=:softmax)

#-- Option 3: using nn-factory
mlp = @mx.chain mx.Variable(:data) =>
  mx.MLP([128, 64, 10])            =>
  mx.SoftmaxOutput(name=:softmax)

# data provider
batch_size = 100
include("mnist-data.jl")
train_provider, eval_provider = get_mnist_providers(batch_size)

# setup model
model = mx.FeedForward(mlp, context=mx.cpu())

# optimizer
optimizer = mx.SGD(η=0.1, μ=0.9, λ=0.00001)

# fit parameters
mx.fit(model, optimizer, train_provider, eval_data=eval_provider, n_epoch=20)

#--------------------------------------------------------------------------------
# Optional, demonstration of the predict API
probs = mx.predict(model, eval_provider)

# collect all labels from eval data
labels = reduce(
  vcat,
  copy(mx.get(eval_provider, batch, :softmax_label)) for batch ∈ eval_provider)
# labels are 0...9
labels .= labels .+ 1

# Now we use compute the accuracy
pred = map(i -> argmax(probs[1:10, i]), 1:size(probs, 2))
correct = sum(pred .== labels)
@printf "Accuracy on eval set: %.2f%%\n" 100correct/length(labels)


text:
module RegressionDiscontinuity

# Julia Standard Library packages and imports
import Base: size, getindex, getproperty, propertynames, show
using DelimitedFiles
using LinearAlgebra
using Statistics; import Statistics: var

# Julia Registry packages
using Distributions
using IntervalSets
using JuMP
using LaTeXStrings
using GLM
using QuadGK
using RecipesBase
using Requires
using Roots
using Setfield
using StatsBase; import StatsBase: fit, weights, nobs
using StatsDiscretizations
using Tables

include("types.jl")
include("running_variable.jl")
include("discretization.jl")
include("load_example_data.jl")
include("kernels.jl")
include("imbens_kalyanaraman.jl")
include("local_linear.jl")
include("minmax_optimal.jl")
include("density_test.jl")

function __init__()
    @require Empirikos="cab608d6-c565-4ea1-96d6-ce5441ba21b0" include("nir.jl")
end

export fit

export RunningVariable,
    RDData,
    Treated,
    Untreated,
    Rectangular,
    bandwidth,
    ImbensKalyanaraman,
    linearweights,
    EickerHuberWhite,
    Homoskedastic,
    NaiveLocalLinearRD,
    ImbensWagerOptRD,
    McCraryTest
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

# # Strong Landau Damping
#
# Electrostatic example of strong Landau damping
# $$
# f(x,v) =\frac{1}{2\pi\sigma^2}  \exp
# \Big( - \frac{v_1^2 + v_2^2}{2\sigma^2} \Big)
# ( 1+\alpha \cos(k x)􏰁),
# $$
#

using GEMPIC
using TimerOutputs

const to = TimerOutput()

# +
function run( nstep :: Int )

   σ, μ = 0.02, 0.0
   kx, α = 1.004355, 0.001
   xmin, xmax = 0, 2π/kx
   nx = 1024
   n_particles = 200000
   mesh = Mesh( xmin, xmax, nx)
   spline_degree = 3

   df = CosSumGaussianSpin([[kx]], [α], [[σ]], [[μ]] )

   particle_group2 = ParticleGroup{1,1}( n_particles, n_spin=3)
   sampler = ParticleSampler{1,1}( :sobol, n_particles)

   sample!(particle_group2, sampler, df, mesh)

   particle_group = ParticleGroup{1,1}( n_particles, n_spin=3)

   for  i_part = 1:n_particles

       x = zeros( 1 )
       v = zeros( 1 )
       s = zeros( 3 )
       w = zeros( 1 )
       x = GEMPIC.get_x(particle_group2, i_part)
       v = GEMPIC.get_v(particle_group2, i_part)

       s1 = GEMPIC.get_spin(particle_group2, i_part, 1)
       s2 = GEMPIC.get_spin(particle_group2, i_part, 2)
       s3 = GEMPIC.get_spin(particle_group2, i_part, 3)

       w = GEMPIC.get_weights(particle_group2, i_part)

       GEMPIC.set_x(particle_group, i_part, x[1])
       GEMPIC.set_v(particle_group, i_part, v[1])
       GEMPIC.set_spin(particle_group, i_part, 1, s1)
       GEMPIC.set_spin(particle_group, i_part, 2, s2)
       GEMPIC.set_spin(particle_group, i_part, 3, s3)
       GEMPIC.set_weights(particle_group, i_part, w[1])

   end

   kernel_smoother2 = ParticleMeshCoupling( mesh, n_particles, spline_degree-2, :galerkin)
   kernel_smoother1 = ParticleMeshCoupling( mesh, n_particles, spline_degree-1, :galerkin)
   kernel_smoother0 = ParticleMeshCoupling( mesh, n_particles, spline_degree, :galerkin)

   rho = zeros(Float64, nx)
   efield_poisson = zeros(Float64, nx)

   # Init!ialize the field solver
   maxwell_solver = Maxwell1DFEM(mesh, spline_degree)
   # efield by Poisson
   solve_poisson!( efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho )

   # Initialize the arrays for the spline coefficients of the fields
   k0 = 12.0523
   E0 = 10
   ww = 12.104827940833333

   Ey(x) = E0*cos(k0*x)
   Ez(x) = E0*sin(k0*x)
   Ay(x) = -E0/ww*sin(k0*x)
   Az(x) = E0/ww*cos(k0*x)

   efield_dofs = [efield_poisson, zeros(Float64, nx), zeros(Float64, nx)]
   afield_dofs = [zeros(Float64, nx), zeros(Float64, nx)]

   l2projection!( efield_dofs[2], maxwell_solver, Ey, spline_degree)
   l2projection!( efield_dofs[3], maxwell_solver, Ez, spline_degree)
   l2projection!( afield_dofs[1], maxwell_solver, Ay, spline_degree)
   l2projection!( afield_dofs[2], maxwell_solver, Az, spline_degree)

   propagator = HamiltonianSplittingSpin( maxwell_solver,
                                      kernel_smoother0,
                                      kernel_smoother1,
                                      kernel_smoother2,
                                      particle_group,
                                      efield_dofs,
                                      afield_dofs)

   efield_dofs_n = propagator.e_dofs

   thdiag = TimeHistoryDiagnosticsSpin( particle_group, maxwell_solver,
                           kernel_smoother0, kernel_smoother1 );

   Δt = 0.002

   for istep in 1:nstep

       @timeit to "Operator HE" GEMPIC.operatorHE( propagator, 0.5Δt)
       @timeit to "Operator Hp" GEMPIC.operatorHp( propagator, 0.5Δt)
       @timeit to "Operator HA" GEMPIC.operatorHA( propagator, 0.5Δt)
       @timeit to "Operator Hs" GEMPIC.operatorHs( propagator, 1.0Δt)
       @timeit to "Operator HA" GEMPIC.operatorHA( propagator, 0.5Δt)
       @timeit to "Operator Hp" GEMPIC.operatorHp( propagator, 0.5Δt)
       @timeit to "Operator HB" GEMPIC.operatorHE( propagator, 0.5Δt)

   end

end

run(10)

println(to)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SingleMass <: Shaft
        H::Float64
        D::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of single mass shaft model. Typically represents the rotor mass.

# Arguments
- `H::Float64`: Rotor inertia constant in MWs/MVA, validation range: `(0, nothing)`
- `D::Float64`: Rotor natural damping in pu, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	δ: rotor angle,
	ω: rotor speed
- `n_states::Int`: SingleMass has 1 state
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct SingleMass <: Shaft
    "Rotor inertia constant in MWs/MVA"
    H::Float64
    "Rotor natural damping in pu"
    D::Float64
    ext::Dict{String, Any}
    "The states are:
	δ: rotor angle,
	ω: rotor speed"
    states::Vector{Symbol}
    "SingleMass has 1 state"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function SingleMass(H, D, ext=Dict{String, Any}(), )
    SingleMass(H, D, ext, [:δ, :ω], 2, InfrastructureSystemsInternal(), )
end

function SingleMass(; H, D, ext=Dict{String, Any}(), states=[:δ, :ω], n_states=2, internal=InfrastructureSystemsInternal(), )
    SingleMass(H, D, ext, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function SingleMass(::Nothing)
    SingleMass(;
        H=0,
        D=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`SingleMass`](@ref) `H`."""
get_H(value::SingleMass) = value.H
"""Get [`SingleMass`](@ref) `D`."""
get_D(value::SingleMass) = value.D
"""Get [`SingleMass`](@ref) `ext`."""
get_ext(value::SingleMass) = value.ext
"""Get [`SingleMass`](@ref) `states`."""
get_states(value::SingleMass) = value.states
"""Get [`SingleMass`](@ref) `n_states`."""
get_n_states(value::SingleMass) = value.n_states
"""Get [`SingleMass`](@ref) `internal`."""
get_internal(value::SingleMass) = value.internal

"""Set [`SingleMass`](@ref) `H`."""
set_H!(value::SingleMass, val) = value.H = val
"""Set [`SingleMass`](@ref) `D`."""
set_D!(value::SingleMass, val) = value.D = val
"""Set [`SingleMass`](@ref) `ext`."""
set_ext!(value::SingleMass, val) = value.ext = val