text:
# TODO: allow for different shapes and mirror sizes

using Plots
using LinearAlgebra

pyplot()

mutable struct Ray
    # Direction vector
    l::Vector{Float64}

    # Position vector
    p::Vector{Float64}

    # Current Index of Refraction
    ior::Float64
end

mutable struct Lens
    # Lens position
    p::Vector{Float64}

    # Lens radius
    r::Float64

    # Lens function
    f
end

function constant(p)
    return 1.5
end

function mag_index(p)
    return cos(sum(p[:].^2))
end

function inverse(p)
    mag = sum(p[:].^2)
    return 1/(mag)
end

function lens_normal_at(ray, lens)
    n = normalize(ray.p .- lens.p)
    if dot(-n, ray.l) < 0
        n *= -1
    end

    return n
end

function inside_of(ray::Ray, lens)
    return inside_of(ray.p, lens)
end

function inside_of(pos, lens)
    x = lens.p[1] - pos[1]
    y = lens.p[2] - pos[2]
    if (x^2 + y^2 <= lens.r^2)
        return true
    else
        return false
    end
end

function refract!(ray, n)
    c = dot(-n, ray.l);
    d = 1.0 - ior^2 * (1.0 - c^2);

    if (d < 0.0)
        return zeros(2);
    end

    ray.l = ior * ray.l + (ior * c - sqrt(d)) * n;
end

# note: light moves as a particular speed with respect to the medium it is
#       moving through, so...
#       n_2*v = n_1*l + (n_1*cos(theta_1) - n_2*cos(theta_2))*n
#       Other approximations: ior = n_1/n_2, c = -n*l
function refract!(ray, lens, ior)
    n = lens_normal_at(ray, lens)
    c = dot(-n, ray.l);
    d = 1.0 - ior^2 * (1.0 - c^2);

    if (d < 0.0)
        return zeros(2);
    end

    ray.l = ior * ray.l + (ior * c - sqrt(d)) * n;
end

mutable struct Mirror
    # Normal vector
    n::Vector{Float64}

    # Position of mirror
    p::Vector{Float64}

    # Mirror size
    scale::Float64

    Mirror(in_n, in_p) = new(in_n, in_p, 2.5)
end

function is_behind(ray, mirror)
    if dot(ray.p.-mirror.p, mirror.n) >= 0
        return true
    else
        return false
    end
end

# note: for reflection, l_x -> l_x, but l_y -> -l_y
#       In this case, the y component of l = cos(theta)*n
#       so new vector: v = l + 2cos(theta)*n
function reflect!(ray, n)
    ray.l .-= 2*dot(ray.l, n).*n
end

function draw_circle(p, r, res)
    return [x .+ (p[1], p[2]) for x in Plots.partialcircle(0, 2pi, res, r)]
end

function plot_rays(positions, mirrors::Vector{Mirror},
                   lenses::Vector{Lens}, filename)
    plt = plot(background_color=:black, aspect_ratio=:equal)
    for i = 1:size(positions)[1]
        plot!(plt, (positions[i,:,1], positions[i,:,2]); label = "ray",
              linecolor=:white)
    end

    for mirror in mirrors
        dir = [-mirror.n[2], mirror.n[1]]
        extents = zeros(2,2)
        extents[1,:] .= mirror.p .- mirror.scale * dir
        extents[2,:] .= mirror.p .+ mirror.scale * dir
        plot!(plt, (extents[:,1], extents[:,2]), label = "mirror")
    end

    for lens in lenses
        circle = draw_circle(lens.p, lens.r, 100)
        plot!(circle; label="lens", linecolor=:lightblue)
    end

    savefig(plt, filename)
end

function step!(ray::Ray, dt)
    ray.p .+= .+ ray.l.*dt/ray.ior
end

function propagate!(rays::Vector{Ray}, mirrors::Vector{Mirror},
                    lenses::Vector{Lens}, positions, n, dt)

    ior = 1.0
    for i = 1:n
        for j = 1:length(rays)
            for mirror in mirrors
                if is_behind(rays[j], mirror)
                    reflect!(rays[j], mirror.n)
                end
            end
            for lens in lenses
                if inside_of(rays[j], lens) && i > 2
                   ior = rays[j].ior/lens.f(rays[j].p .- lens.p)
                   rays[j].ior = lens.f(rays[j].p .- lens.p)
                   refract!(rays[j], lens, ior)
                elseif !inside_of(rays[j], lens) && i > 2 &&
                       rays[j].ior != 1
                    ior = rays[j].ior
                    rays[j].ior = 1
                    refract!(rays[j], lens, ior)
                end
                if rays[j].l == zeros(2)
                    reflect!(rays[j], lens_normal_at(rays[j], lens))
                end
            end
            step!(rays[j], dt)
            positions[j, i,:] .= rays[j].p
        end
    end
end

function parallel_propagate(ray_num, n; filename="check.png", dt = 0.1)

    rays = [Ray([1, 0],
            [0.1, float(i)],
            1.0) for i = 1:ray_num]

    #lenses = [Lens([10, 6], 6, constant)]
    lenses = [Lens([10, 6], 6, constant)]
    #lenses = [Lens([-15, -10], 5, 1.5)]
    mirrors = [Mirror([0.0, -1.0],[5.0, 0.0])]
#=
    mirrors = [Mirror([0.0, 1.0],[5.0, 5.0]),
               Mirror([0.0, -1.0],[5.0, 0.0]),
               Mirror([1.0, 0.0],[0.0, 2.5]),
               Mirror([-1.0, 0.0],[10, 2.5])]
=#

    positions = zeros(ray_num, n, 2)

    propagate!(rays, mirrors, lenses, positions, n, dt)

    plot_rays(positions, mirrors, lenses, filename)

    return rays
end


text:
abstract type AbstractIndividual end

variables(i::AbstractIndividual) = error("Function `variables` not implemented")

objectives(i::AbstractIndividual) = error("Function `objectives` not implemented")

constraints(i::AbstractIndividual) = error("Function `constraints` not implemented")

mutable struct Individual <: AbstractIndividual
    variables::AbstractVector
    objectives::Union{AbstractVector, Number, Nothing}
    constraints::Union{AbstractVector, Nothing}

    Individual(var) = new(var, nothing, nothing)
    Individual(var, obj) = new(var, obj, nothing)
    Individual(var, obj, cst) = new(var, obj, cst)
end

variables(i::Individual) = i.variables

objectives(i::Individual) = i.objectives

constraints(i::Individual) = i.constraints

copy(n::Nothing) = nothing

copy(i::Individual) = Individual(copy(i.variables), copy(i.objectives), copy(i.constraints))

ext: jl


text:
@testset "$TEST $G" begin

N = 10
d = 2
points = rand(d, N)
g, w = euclidean_graph(points)
@test nv(g) == N
@test ne(g) == N*(N-1) ÷ 2
@test maximum(x->x[2], w) <= sqrt(d)
@test minimum(x->x[2], w) >= 0
@test maximum(points) <= 1
@test minimum(points) >= 0.

g, w = euclidean_graph(points, bc=:periodic)
@test maximum(x->x[2], w) <= sqrt(d/2)
@test minimum(x->x[2], w) >= 0.

@test_throws ErrorException euclidean_graph(points, L=0.01,  bc=:periodic)
@test_throws ErrorException euclidean_graph(points, bc=:ciao)

end # testset


text:
# ERROR: ch_A not defined
global ch_A::BigInt = 13591409
global ch_B::BigInt = 545140134
global ch_C::BigInt = 640320
global ch_C3::BigInt = ch_C ^ 3
global ch_D::BigInt = 12

function ch_split(a::BigInt, b::BigInt)
    if b - a == 1
        g::BigInt = (6 * b - 5) * (2 * b - 1) * (6 * b - 1)
        g, div(ch_C3 * b ^ 3, 24), (-1^b) * g * (b * ch_B + ch_A)
    else
        mid = div(a + b, 2)
        g1, p1, q1 = ch_split(a, mid)
        g2, p2, q2 = ch_split(mid, b)
        g1 * g2, p1 * p2, q1 * p2 + q2 * p1
    end
end

function pi_again(digits)
    num_terms::BigInt = iround(2::BigInt + digits / 14.181647462)
    sqrt_C::BigInt = isqrt(ch_C * 100 ^ digits)
    g, p, q = ch_split(0, num_terms)
    div(p * ch_C * sqrt_C, ch_D * (q + p * ch_A))
end

printlin(map(pi_again, 0:10))


text:
using Test
using Random
using Statistics
using AdaptiveMCMC

Random.seed!(1234)

# Standard normal test target in 2d
log_p(x) = -(x[1]^2+x[2]^2)/2
x0 = zeros(2)
n = 100_000
test_means(X, atol=0.05) = all(abs.(mapslices(mean, X, dims=2)) .<= atol)
test_vars(X, atol=0.05) = all(abs.(mapslices(var, X, dims=2) .- 1.0) .<= atol)
test_stats(X, atol=0.05) = test_vars(X, atol) && test_means(X, atol)
test_acc(acc, acc_opt=0.234, atol=0.05) = all(abs.(acc .- acc_opt) .<= atol)

o_ram = adaptive_rwm(x0, log_p, n; algorithm=:ram)
@test test_acc(o_ram.accRWM)
@test test_stats(o_ram.X)
o_ram_32 = adaptive_rwm(Float32.(x0), log_p, n; algorithm=:ram)
@test test_acc(o_ram_32.accRWM)
@test test_stats(o_ram_32.X)

o_am = adaptive_rwm(x0, log_p, n; algorithm=:am)
@test o_am.S[1].L.L[1] ≈ 1.0 atol=0.05
@test test_stats(o_am.X)
o_am_32 = adaptive_rwm(x0, log_p, n; algorithm=:am)
@test o_am_32.S[1].L.L[1] ≈ 1.0 atol=0.05
@test test_stats(o_am_32.X)

o_asm = adaptive_rwm(x0, log_p, n; algorithm=:asm)
@test test_acc(o_asm.accRWM)
@test test_stats(o_asm.X)
o_asm_32 = adaptive_rwm(Float32.(x0), log_p, n; algorithm=:asm)
@test test_acc(o_asm_32.accRWM)
@test test_stats(o_asm_32.X)

o_aswam = adaptive_rwm(x0, log_p, n; algorithm=:aswam)
@test o_aswam.S[1].AM.L.L[1] ≈ 1.0 atol=0.05
@test test_acc(o_aswam.accRWM)
@test test_stats(o_aswam.X)
o_aswam_32 = adaptive_rwm(Float32.(x0), log_p, n; algorithm=:aswam)
@test o_aswam_32.S[1].AM.L.L[1] ≈ 1.0 atol=0.10
@test test_acc(o_aswam_32.accRWM)
@test test_stats(o_aswam_32.X)

o_apt = adaptive_rwm(x0, log_p, n; L=3, all_levels=true)
@test test_acc(o_apt.accRWM)
@test test_acc(o_apt.accSW)
@test test_stats(o_apt.X)
@test test_means(o_apt.allX[2], 0.10)
@test test_means(o_apt.allX[3], 0.15)
o_apt_32 = adaptive_rwm(Float32.(x0), log_p, n; L=3, all_levels=true)
@test test_acc(o_apt_32.accRWM)
@test test_acc(o_apt_32.accSW)
@test test_stats(o_apt_32.X)
@test test_means(o_apt_32.allX[2], 0.10)
@test test_means(o_apt_32.allX[3], 0.15)


text:
function interpolate_level_timeseries(t_o, level, dt)

    n_t = size(t_o)[1]
    n_ip = Int64(ceil((t_o[end]-t_o[1])/dt))

    T = eltype(level)
    t_interpolated = zeros(T, n_ip)
    for i in range(1, n_ip, step=1)
        t_interpolated[i] = t_o[1] + (i-1)*dt
    end

    # Interpolate the data
    f_level_interp = PCHIPInterpolation.Interpolator(t_o, reshape(level, (n_t, )))
    level_interpolated = f_level_interp.(t_interpolated)

    return level_interpolated
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

bt = backtrace()
have_backtrace = false
for l in bt
    lkup = ccall(:jl_lookup_code_address, Any, (Ptr{Void}, Cint), l, true)
    if lkup[1][1] == :backtrace
        @test lkup[1][5] == false # fromC
        global have_backtrace = true
        break
    end
end

@test have_backtrace

# Test location information for inlined code (ref issues #1334 #12544)
module test_inline_bt
using Test

function get_bt_frames(functionname, bt)
    for i = 1:length(bt)
        lkup = Base.StackTraces.lookup(bt[i])
        lkup[end].func == functionname && (return lkup)
    end
end

# same-file inline
eval(Expr(:function, Expr(:call, :test_inline_1),
          Expr(:block, Expr(:line, 99, Symbol("backtrace.jl")),
                       Expr(:block, Expr(:line, 42),
                                    Expr(:meta, :push_loc, Symbol("backtrace.jl"), :inlfunc),
                                    Expr(:line, 37),
                                    Expr(:call, :throw, "foo"),
                                    Expr(:meta, :pop_loc),
                                    Expr(:line, 99)))))

@test functionloc(test_inline_1) == ("backtrace.jl", 99)
try
    test_inline_1()
    error("unexpected")
catch err
    lkup = get_bt_frames(:test_inline_1, catch_backtrace())
    @test length(lkup) == 2
    @test endswith(string(lkup[2].file), "backtrace.jl")
    @test lkup[2].line == 42
    @test lkup[1].func == :inlfunc
    @test endswith(string(lkup[1].file), "backtrace.jl")
    @test lkup[1].line == 37
end

# different-file inline
const absfilepath = Sys.iswindows() ? "C:\\foo\\bar\\baz.jl" : "/foo/bar/baz.jl"
eval(Expr(:function, Expr(:call, :test_inline_2),
          Expr(:block, Expr(:line, 81, Symbol("backtrace.jl")),
                       Expr(:block, Expr(:meta, :push_loc, Symbol(absfilepath)),
                                    Expr(:line, 111),
                                    Expr(:call, :throw, "foo"),
                                    Expr(:meta, :pop_loc),
                                    Expr(:line, 99)))))

@test functionloc(test_inline_2) == ("backtrace.jl", 81)
try
    test_inline_2()
    error("unexpected")
catch err
    lkup = get_bt_frames(:test_inline_2, catch_backtrace())
    @test length(lkup) == 2
    @test endswith(string(lkup[2].file), "backtrace.jl")
    @test lkup[2].line == 81
    @test string(lkup[1].file) == absfilepath
    @test lkup[1].line == 111
end

end # module

#issue 12977: line numbers for kwarg methods.
linenum = @__LINE__; f12977(; args...) = ()
loc = functionloc(f12977)
@test endswith(loc[1], "backtrace.jl")
@test loc[2] == linenum

@noinline function test_throw_commoning(x)
    if x==1; throw(AssertionError()); end
    if x==2; throw(AssertionError()); end
end

let
    local b1, b2
    try
        test_throw_commoning(1)
    catch
        b1 = stacktrace(catch_backtrace())
    end
    try
        test_throw_commoning(2)
    catch
        b2 = stacktrace(catch_backtrace())
    end
    i1 = findfirst(frame -> frame.func === :test_throw_commoning, b1)
    i2 = findfirst(frame -> frame.func === :test_throw_commoning, b2)
    @test i1 > 0 && i2 > 0
    @test b1[i1].line != b2[i2].line
end

module BackTraceTesting

using Test

@inline bt2() = backtrace()
@inline bt1() = bt2()
bt() = bt1()

lkup = map(StackTraces.lookup, bt())
hasbt = hasbt2 = false
for sfs in lkup
    for sf in sfs
        if sf.func == :bt
            global hasbt = true
        end
        if sf.func == :bt2
            global hasbt2 = true
        end
    end
end
@test hasbt
if Base.JLOptions().can_inline != 0
    @test_broken hasbt2
else
    @test hasbt2
end

function btmacro()
    ret = @timed backtrace()
    ret[1]
end
lkup = map(StackTraces.lookup, btmacro())
hasme = hasbtmacro = false
for sfs in lkup
    for sf in sfs
        if sf.func == Symbol("macro expansion")
            global hasme = true
        end
        if sf.func == :btmacro
            global hasbtmacro = true
        end
    end
end
@test hasme
@test hasbtmacro

end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
#= convex.jl
Philip Barrett, pobarrett@gmail.com
27may2016, Chicago

Calculates the convex hull of a set of points
=#

"""
    acw( p1, p2, p3 )
Positive if p1 -> p2 -> p3 is anti-clockwise
"""
function acw( p1, p2, p3 )
  return (p2[1] - p1[1]) * (p3[2] - p1[2]) - (p2[2] - p1[2]) * (p3[1] - p1[1])
end

"""
    grahamScan( pts::Matrix )
Computes a convex hull using the Graham Scan algortithm NB: Returns a polygon object even when a matrix is submitted
"""

function grahamScan( pts::Matrix )

  N = size(pts)[1]
      # Number of points
  lex = sortrows( pts[ :, [2, 1 ] ] )[ :, [ 2, 1 ] ]
      # Sorts the points lexicographically by y then x
  P = lex[ 1, : ]
      # The initial point
  otherpts = lex[ 2:end, : ]
      # The other points
  orderpts = [ P; otherpts ; P ]
      # Initialize the anti-clockwise ordered points
  cah = zeros( N - 1 )
      # The cosine of the angle between each point and P.  Will be in
      # [-1,1] because of choice of P

  ## Order the points anti-clockwise ##
  for( i in 1:(N-1) )
    diff = otherpts[i,:] - P
    cah[i] = diff[1] / norm(diff)
        # Cosine is adjacent over hypotenuse
  end
  orderpts[2:N,:] = otherpts[ sortperm(cah, rev=true), : ]
      # Order the points by the angle measure

  ## Create the output ##
  out = similar( pts )
  out[ 1:2, : ] = orderpts[1:2,:]
      # Initialize
  M = 2
      # Counts number of rows in convex hull
  for( i in 3:N )
    while( acw( out[ M-1, : ], out[ M, : ], orderpts[ i, : ] ) < 0 )
      if( M > 2 )
        M -= 1
      elseif( i == N )
        break
      else
        i += 1
      end
    end
    M += 1
        # Increment counter
    out[ M, : ] = orderpts[ i, : ]
        # Add to the output if we have an acw angle
  end
  return Polygon( pts=out[ 1:M, : ] )
end

function grahamScan( poly::Polygon )
    return grahamScan( poly.pts )
end

"""
    convexhull( poly::Polygon )
Computes a convex hull using the Graham Scan algortithm from the CHull2D package
"""
function chull( poly::Polygon )
  return Polygon( pts = chull(poly.pts) )
end

function chull( pts::Matrix )
  N = size(pts)[1]
      # Number of points
  oldpts = Point{2, Float64}[Point(pts[i, 1], pts[i, 2])
                            for i=1:N]
      # Convert poly.pts to an array of points
  ch = CHull2d.convexhull(oldpts, _at=false)
      # The convex hull object
  N_newpts = length(ch.extremepoints)
  newpts = Array(Float64, N_newpts, 2)
      # Create the container for the new points
  for i=1:N_newpts
      curr_pt = ch.extremepoints[i]
      newpts[i, 1] = curr_pt[1]
      newpts[i, 2] = curr_pt[2]
  end
  return newpts
end


text:
# SPDX-License-Identifier: X11
# 2020-12-09
# Day 5, Part 1

function seatid(str::AbstractString)::Integer
  row = -1
  s = 0
  e = 127
  for (i, c) ∈ enumerate(str[1:7])
    c == 'F' ? e -= 2^(7 - i) : s += 2^(7 - i)
  end
  @assert s == e
  row = s

  col = -1
  s = 0
  e = 7
  for (i, c) ∈ enumerate(str[8:10])
    c == 'L' ? e -= 2^(3 - i) : s += 2^(3 - i)
  end
  @assert s == e
  col = s

  return row * 8 + col
end

(println ∘ maximum)(seatid, readlines())

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
export filter

struct Filter{I,F}
    rows::I
    predicate::F
end

"""
    filter(rows, predicate)
    rows |> filter(predicate)

Filter `rows` to include only rows that fulfil `predicate`.
"""
filter(rows, predicate) = Filter{typeof(rows), typeof(predicate)}(rows, predicate)
filter(predicate) = rows -> filter(rows, predicate)

Base.IteratorEltype(::Type{Filter{I,F}}) where {I,F} = Base.HasEltype()
Base.eltype(::Type{Filter{I,F}}) where {I,F} = Base.eltype(I)
Base.IteratorSize(::Type{Filter{I,F}}) where {I,F} = Base.SizeUnknown()

Tables.istable(m::Filter{I,F}) where {I,F} = true
Tables.rowaccess(m::Filter{I,F}) where {I,F} = true
Tables.schema(m::Filter{I,F}) where {I,F} = Tables.schema(m.rows)

function Base.iterate(m::Filter{I,F}) where {I,F}
    elt_s = iterate(m.rows)
    state = elt_s
    iterate(m, state)
end

function Base.iterate(m::Filter{I,F}, state) where {I,F}
    elt_s = state

    while true
        # done
        elt_s == nothing && return nothing

        # predicate fulfilled
        m.predicate(elt_s[1]) && return elt_s[1], iterate(m.rows, elt_s[2])

        # forward to next row
        elt_s = iterate(m.rows, elt_s[2])
    end
end


text:
using Ft2k

# this is from watchmanfxns.jl
success = Ft2k.check_server_updating(Ft2k.conf) # 3.222670 seconds on macbook

tdict, fdict, ids, sqlite_db = Ft2k.get_dbvalues()

# create summary for the day to go out after 9

Ft2k.dailycheck(tdict, fdict, ids, sqlite_db)

# to get assigned feedings from sqlite_db:
# id = "f2b00d00"
# sqlite_db[1][sqlite_db[1][:,1] .== id, :]


text:
using Dates
using Printf

mutable struct TimeInfo{T<:Real}
    it      :: Int
    t       :: T
    dt      :: T
    runtime :: T
end
TimeInfo(it::Int, t::T, dt::T, runtime::T) where {T} = TimeInfo{T}(it, t, dt, runtime)
TimeInfo() = TimeInfo(0, 0.0, 0.0, 0.0)

mutable struct Field{A,G}
    name  :: String
    data  :: A
    chart :: G
end
Field(name::String, data, chart::Chart) = Field{typeof(data),typeof(chart)}(name, data, chart)

function out_info(it::Integer, t::Real, time_per_hour::Real, f, label::String, info_every::Integer,
                  header_every::Integer)

    if it % header_every == 0
        println("-------------------------------------------------------------")
        println("Iteration      Time | Time per hour  |           $label")
        println("                    |                |   minimum      maximum")
        println("-------------------------------------------------------------")
    end

    if it % info_every == 0
        @printf "%9d %9.3f | %9.3f      | %9.4g    %9.4g\n" it t time_per_hour minimum(f) maximum(f)
    end

    nothing
end

struct Output{T}
    dir              :: String
    prefix           :: String
    software         :: String
    software_version :: String
    tinfo            :: TimeInfo{T}

    function Output{T}(dir::String, prefix::String,
                       software::String, software_version::String,
                       tinfo::TimeInfo{T}) where {T}
        # if no name specified error out
        if dir == ""
            error("Output folder cannot be empty string.")
        end

        # create folder if it doesn't exist already
        if !isdir(dir)
            mkdir(dir)
        end

        if prefix == ""
            prefix = "data_"
        end

        new(dir, prefix, software, software_version, tinfo)
    end
end
function Output(dir::String, prefix::String, tinfo::TimeInfo{T}) where {T<:Real}
    software         = "Jecco"
    software_version = "0.8.0"
    Output{T}(dir, prefix, software, software_version, tinfo)
end

# make Output a callable struct

function (out::Output)(fields::Union{Vector, Tuple}; params::Union{Tuple, NamedTuple, Dict}=())
    it   = out.tinfo.it
    time = out.tinfo.t
    dt   = out.tinfo.dt

    filename = "$(out.prefix)$(lpad(string(it), 8, string(0))).h5"
    fullpath = abspath(out.dir, filename)

    # check if file already exists
    if isfile(fullpath)
        mode = "r+"
        firsttime = false
    else
        mode = "cw"
        firsttime = true
    end

    # open file
    fid = h5open(fullpath, mode)

    # create openPMD structure
    if firsttime
        setup_openpmd_file(out, fid)
    end

    if haskey(fid, "data/$it")
        grp = fid["data/$it/fields"]
    else
        grp = create_group(fid, "data/$it")
        attributes(grp)["time"] = time
        attributes(grp)["dt"] = dt
        grp = create_group(grp, "fields")
    end

    # write given parameters as attributes
    if firsttime
        for key in keys(params)
            name = String(key)
            val  = params[key]
            attributes(grp)[name] = val
        end
    end

    for field in fields
        write_hdf5(out, grp, field)
    end

    # close file
    close(fid)

    nothing
end

(out::Output)(fields::Vararg{Field,N}; params=()) where {N} =
    (out::Output)((fields...,), params=params)

write_hdf5(out::Output, grp::HDF5.Group, field::Field) =
    write_hdf5(out, grp, field.name, field.data, field.chart)

function write_hdf5(out::Output, grp::HDF5.Group, fieldname::String, data::AbstractArray,
                    chart::Chart)
    dset, dtype = create_dataset(grp, fieldname, data)
    # write actual data
    write_dataset(dset, dtype, data)
    setup_openpmd_mesh(dset, chart)

    nothing
end

function setup_openpmd_file(out::Output, fid::HDF5.File)
    timenow = now()
    date    = Dates.format(timenow, "yyyy-mm-dd HH:MM:SS")

    attributes(fid)["software"] = out.software
    attributes(fid)["softwareVersion"] = out.software_version
    attributes(fid)["openPMD"] = "1.1.0"
    attributes(fid)["openPMDextension"] = 0
    attributes(fid)["author"] = try
        ENV["USER"]
    catch e
        if isa(e, KeyError) # probably on windows
            splitdir(homedir())[end]
        else
            throw(e)
        end
    end
    attributes(fid)["date"] = string(date)
    attributes(fid)["iterationEncoding"] = "fileBased"
    attributes(fid)["iterationFormat"] = "$(out.prefix)%T.h5"

    attributes(fid)["basePath"] = "/data/%T/"
    attributes(fid)["meshesPath"] = "fields/"
end

openpmd_geometry(coord::CartesianCoord) = "cartesian"
openpmd_geometry(coord::GaussLobattoCoord) = "other"

function openpmd_geometry(chart::Chart)
    geometries = [openpmd_geometry(coord) for coord in chart.coords]
    if all(geometries .== "cartesian")
        geometry = "cartesian"
    else
        geometry = "other"
    end
    geometry
end

function setup_openpmd_mesh(dset::HDF5.Dataset, coord::AbstractCoord)
    attributes(dset)["geometry"]         = openpmd_geometry(coord)
    attributes(dset)["gridGlobalOffset"] = coord.min
    attributes(dset)["gridSpacing"]      = Jecco.delta(coord)
    attributes(dset)["gridMax"]          = coord.max
    attributes(dset)["gridType"]         = Jecco.coord_type(coord)
    attributes(dset)["axisLabels"]       = coord.name
    nothing
end
# there is no method to write BigFloat attributes in HDF5, so convert them to Float64
function setup_openpmd_mesh(dset::HDF5.Dataset, coord::AbstractCoord{N,BigFloat}) where {N}
    T = Float64
    attributes(dset)["geometry"]         = openpmd_geometry(coord)
    attributes(dset)["gridGlobalOffset"] = T(coord.min)
    attributes(dset)["gridSpacing"]      = T(Jecco.delta(coord))
    attributes(dset)["gridMax"]          = T(coord.max)
    attributes(dset)["gridType"]         = T(Jecco.coord_type(coord))
    attributes(dset)["axisLabels"]       = coord.name
    nothing
end

function setup_openpmd_mesh(dset::HDF5.Dataset, chart::Chart)
    mins      = Jecco.min(chart)
    deltas    = Jecco.delta(chart)
    maxs      = Jecco.max(chart)
    gridtypes = Jecco.coord_type(chart)
    names     = Jecco.name(chart)

    attributes(dset)["geometry"]         = openpmd_geometry(chart)

    # Julia, like Fortran and Matlab, stores arrays in column-major order. HDF5
    # uses C's row-major order, and consequently every array's dimensions are
    # inverted compared to what is seen with tools like h5dump. This is the same
    # convention as for the Fortran and Matlab HDF5 interfaces. For consistency,
    # then, we here flip the order of the grid arrays (and remember to flip back
    # when reading data in). The advantage is that no data rearrangement takes
    # place when reading or writing the data itself, which is more intensive.
    attributes(dset)["gridGlobalOffset"] = mins[end:-1:1]
    attributes(dset)["gridSpacing"]      = deltas[end:-1:1]
    attributes(dset)["gridMax"]          = maxs[end:-1:1]
    attributes(dset)["gridType"]         = gridtypes[end:-1:1]
    attributes(dset)["axisLabels"]       = names[end:-1:1]
    nothing
end
function setup_openpmd_mesh(dset::HDF5.Dataset, chart::Chart{N,BigFloat}) where{N}
    T = Float64
    mins      = T.(Jecco.min(chart))
    deltas    = T.(Jecco.delta(chart))
    maxs      = T.(Jecco.max(chart))
    gridtypes = Jecco.coord_type(chart)
    names     = Jecco.name(chart)

    attributes(dset)["geometry"]         = openpmd_geometry(chart)
    attributes(dset)["gridGlobalOffset"] = mins[end:-1:1]
    attributes(dset)["gridSpacing"]      = deltas[end:-1:1]
    attributes(dset)["gridMax"]          = maxs[end:-1:1]
    attributes(dset)["gridType"]         = gridtypes[end:-1:1]
    attributes(dset)["axisLabels"]       = names[end:-1:1]
    nothing
end


text:
include("resources/HTTPMessages.jl")

import Base.==
using ..HTTP
using ..HTTPMessages
using HTTP.Messages
using HTTP.Parsers
using HTTP.Strings: tocameldash
using ..Test

import Base.==

const strict = false

==(a::Request,b::Request) = (a.method         == b.method)    &&
                            (a.version        == b.version)   &&
                            (a.headers        == b.headers)   &&
                            (a.body           == b.body)

macro errmsg(expr)
    esc(quote
        try
            $expr
        catch e
            sprint(show, e)
        end
    end)
end

@testset "HTTP.parser" begin
    @testset "parse - Strings" begin
        @testset "Requests - $request" for request in requests
            r = parse(Request, request.raw)

            if r.method == "CONNECT"
                host, port = split(r.target, ":")
                @test host == request.host
                @test port == request.port
            else
                if r.target == "*"
                    @test r.target == request.request_path
                else
                    target = parse(HTTP.URI, r.target)
                    @test target.query == request.query_string
                    @test target.fragment == request.fragment
                    @test target.path == request.request_path
                    @test target.host == request.host
                    @test target.userinfo == request.userinfo
                    @test target.port in (request.port, "80", "443")
                    @test string(target) == request.request_url
                end
            end

            r_headers = [tocameldash(n) => String(v) for (n,v) in r.headers]

            @test r.version.major == request.http_major
            @test r.version.minor == request.http_minor
            @test r.method == string(request.method)
            @test length(r.headers) == request.num_headers
            @test r_headers == request.headers
            @test String(r.body) == request.body

            @test_broken HTTP.http_should_keep_alive(HTTP.DEFAULT_PARSER) == request.should_keep_alive
            @test_broken String(collect(upgrade[])) == request.upgrade
        end

        @testset "Request - Headers" begin
            reqstr = "GET http://www.techcrunch.com/ HTTP/1.1\r\n" *
                   "Host: www.techcrunch.com\r\n" *
                   "User-Agent: Fake\r\n" *
                   "Accept: text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8\r\n" *
                   "Accept-Language: en-us,en;q=0.5\r\n" *
                   "Accept-Encoding: gzip,deflate\r\n" *
                   "Accept-Charset: ISO-8859-1,utf-8;q=0.7,*;q=0.7\r\n" *
                   "Keep-Alive: 300\r\n" *
                   "Content-Length: 7\r\n" *
                   "Proxy-Connection: keep-alive\r\n\r\n1234567"
            req = Request("GET", "http://www.techcrunch.com/", ["Host"=>"www.techcrunch.com","User-Agent"=>"Fake","Accept"=>"text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8","Accept-Language"=>"en-us,en;q=0.5","Accept-Encoding"=>"gzip,deflate","Accept-Charset"=>"ISO-8859-1,utf-8;q=0.7,*;q=0.7","Keep-Alive"=>"300","Content-Length"=>"7","Proxy-Connection"=>"keep-alive"])
            req.body = HTTP.bytes("1234567")
            @test parse(Request,reqstr).headers == req.headers
            @test parse(Request,reqstr) == req
        end

        @testset "Request - Hostname - URL" begin
            reqstr = "GET / HTTP/1.1\r\n" *
                   "Host: foo.com\r\n\r\n"
            req = Request("GET", "/", ["Host"=>"foo.com"])
            @test parse(Request, reqstr) == req
        end

        @testset "Request - Hostname - Path" begin
            reqstr = "GET //user@host/is/actually/a/path/ HTTP/1.1\r\n" *
                   "Host: test\r\n\r\n"
            req = Request("GET", "//user@host/is/actually/a/path/",
                       ["Host"=>"test"])
            @test parse(Request, reqstr) == req
        end

        @testset "Request - Hostname - Path - ParseError" begin
            reqstr = "GET  HTTP/1.1\r\n" *
                    "Host: test\r\n\r\n"
            @test_throws HTTP.ParseError parse(Request, reqstr)
        end

        @testset "Request - Hostname - URL - ParseError" begin
            reqstr = "GET ../../../../etc/passwd HTTP/1.1\r\n" *
                "Host: test\r\n\r\n"
            @test_throws HTTP.ParseError HTTP.URI(parse(Request, reqstr).target)
        end

        @testset "Request - HTTP - Bytes" begin
            reqstr = "POST / HTTP/1.1\r\n" *
                   "Host: foo.com\r\n" *
                   "Transfer-Encoding: chunked\r\n\r\n" *
                   "3\r\nfoo\r\n" *
                   "3\r\nbar\r\n" *
                   "0\r\n" *
                   "Trailer-Key: Trailer-Value\r\n" *
                   "\r\n"
            req = Request("POST", "/",
                        ["Host"=>"foo.com", "Transfer-Encoding"=>"chunked", "Trailer-Key"=>"Trailer-Value"])
            req.body = HTTP.bytes("foobar")
            @test parse(Request, reqstr) == req
        end

        @test_skip @testset "Request - HTTP - ParseError" begin
            reqstr = "POST / HTTP/1.1\r\n" *
               "Host: foo.com\r\n" *
               "Transfer-Encoding: chunked\r\n" *
               "Content-Length: 9999\r\n\r\n" * # to be removed.
               "3\r\nfoo\r\n" *
               "3\r\nbar\r\n" *
               "0\r\n" *
               "\r\n"

            @test_throws HTTP.ParseError parse(Request, reqstr)
        end

        @testset "Request - URL" begin
            reqstr = "CONNECT www.google.com:443 HTTP/1.1\r\n\r\n"
            req = Request("CONNECT", "www.google.com:443")
            @test parse(Request, reqstr) == req
        end

        @testset "Request - Localhost" begin
            reqstr = "CONNECT 127.0.0.1:6060 HTTP/1.1\r\n\r\n"
            req = Request("CONNECT", "127.0.0.1:6060")
            @test parse(Request, reqstr) == req
        end

        @testset "Request - RPC" begin
            reqstr = "CONNECT /_goRPC_ HTTP/1.1\r\n\r\n"
            req = HTTP.Request("CONNECT", "/_goRPC_")
            @test parse(Request, reqstr) == req
        end

        @testset "Request - NOTIFY" begin
            reqstr = "NOTIFY * HTTP/1.1\r\nServer: foo\r\n\r\n"
            req = Request("NOTIFY", "*", ["Server"=>"foo"])
            @test parse(Request, reqstr) == req
        end

        @testset "Request - OPTIONS" begin
            reqstr = "OPTIONS * HTTP/1.1\r\nServer: foo\r\n\r\n"
            req = Request("OPTIONS", "*", ["Server"=>"foo"])
            @test parse(Request, reqstr) == req
        end

        @testset "Request - GET" begin
            reqstr = "GET / HTTP/1.1\r\nHost: issue8261.com\r\nConnection: close\r\n\r\n"
            req = Request("GET", "/", ["Host"=>"issue8261.com", "Connection"=>"close"])
            @test parse(Request, reqstr) == req
        end

        @testset "Request - HEAD" begin
            reqstr = "HEAD / HTTP/1.1\r\nHost: issue8261.com\r\nConnection: close\r\nContent-Length: 0\r\n\r\n"
            req = Request("HEAD", "/", ["Host"=>"issue8261.com", "Connection"=>"close", "Content-Length"=>"0"])
            @test parse(Request, reqstr) == req
        end

        @testset "Request - POST" begin
            reqstr = "POST /cgi-bin/process.cgi HTTP/1.1\r\n" *
                   "User-Agent: Mozilla/4.0 (compatible; MSIE5.01; Windows NT)\r\n" *
                   "Host: www.tutorialspoint.com\r\n" *
                   "Content-Type: text/xml; charset=utf-8\r\n" *
                   "Content-Length: 19\r\n" *
                   "Accept-Language: en-us\r\n" *
                   "Accept-Encoding: gzip, deflate\r\n" *
                   "Connection: Keep-Alive\r\n\r\n" *
                   "first=Zara&last=Ali\r\n\r\n"
            req = Request("POST", "/cgi-bin/process.cgi",
                        ["User-Agent"=>"Mozilla/4.0 (compatible; MSIE5.01; Windows NT)",
                         "Host"=>"www.tutorialspoint.com",
                         "Content-Type"=>"text/xml; charset=utf-8",
                         "Content-Length"=>"19",
                         "Accept-Language"=>"en-us",
                         "Accept-Encoding"=>"gzip, deflate",
                         "Connection"=>"Keep-Alive"])
            req.body = HTTP.bytes("first=Zara&last=Ali")
            @test parse(Request, reqstr) == req
        end

        @testset "Request - Pair{}" begin
            r = parse(Request,"GET / HTTP/1.1\r\n" * "Test: Düsseldorf\r\n\r\n")
            @test Pair{String,String}[r.headers...] == ["Test" => "Düsseldorf"]
        end

        @testset "Request - Methods - 1" begin
            for m in ["GET", "PUT", "M-SEARCH", "FOOMETHOD"]
                r = parse(Request,"$m / HTTP/1.1\r\n\r\n")
                @test r.method == string(m)
            end
        end

        @testset "Request - Methods - 2" begin
            for m in ("ASDF","C******","COLA","GEM","GETA","M****","MKCOLA","PROPPATCHA","PUN","PX","SA")
                @test parse(Request,"$m / HTTP/1.1\r\n\r\n").method == m
            end
        end

        @testset "Request - HTTPS" begin
            reqstr = "GET / HTTP/1.1\r\n" *
            "X-SSL-FoooBarr:   -----BEGIN CERTIFICATE-----\r\n" *
            "\tMIIFbTCCBFWgAwIBAgICH4cwDQYJKoZIhvcNAQEFBQAwcDELMAkGA1UEBhMCVUsx\r\n" *
            "\tETAPBgNVBAoTCGVTY2llbmNlMRIwEAYDVQQLEwlBdXRob3JpdHkxCzAJBgNVBAMT\r\n" *
            "\tAkNBMS0wKwYJKoZIhvcNAQkBFh5jYS1vcGVyYXRvckBncmlkLXN1cHBvcnQuYWMu\r\n" *
            "\tdWswHhcNMDYwNzI3MTQxMzI4WhcNMDcwNzI3MTQxMzI4WjBbMQswCQYDVQQGEwJV\r\n" *
            "\tSzERMA8GA1UEChMIZVNjaWVuY2UxEzARBgNVBAsTCk1hbmNoZXN0ZXIxCzAJBgNV\r\n" *
            "\tBAcTmrsogriqMWLAk1DMRcwFQYDVQQDEw5taWNoYWVsIHBhcmQYJKoZIhvcNAQEB\r\n" *
            "\tBQADggEPADCCAQoCggEBANPEQBgl1IaKdSS1TbhF3hEXSl72G9J+WC/1R64fAcEF\r\n" *
            "\tW51rEyFYiIeZGx/BVzwXbeBoNUK41OK65sxGuflMo5gLflbwJtHBRIEKAfVVp3YR\r\n" *
            "\tgW7cMA/s/XKgL1GEC7rQw8lIZT8RApukCGqOVHSi/F1SiFlPDxuDfmdiNzL31+sL\r\n" *
            "\t0iwHDdNkGjy5pyBSB8Y79dsSJtCW/iaLB0/n8Sj7HgvvZJ7x0fr+RQjYOUUfrePP\r\n" *
            "\tu2MSpFyf+9BbC/aXgaZuiCvSR+8Snv3xApQY+fULK/xY8h8Ua51iXoQ5jrgu2SqR\r\n" *
            "\twgA7BUi3G8LFzMBl8FRCDYGUDy7M6QaHXx1ZWIPWNKsCAwEAAaOCAiQwggIgMAwG\r\n" *
            "\tA1UdEwEB/wQCMAAwEQYJYIZIAYb4QgHTTPAQDAgWgMA4GA1UdDwEB/wQEAwID6DAs\r\n" *
            "\tBglghkgBhvhCAQ0EHxYdVUsgZS1TY2llbmNlIFVzZXIgQ2VydGlmaWNhdGUwHQYD\r\n" *
            "\tVR0OBBYEFDTt/sf9PeMaZDHkUIldrDYMNTBZMIGaBgNVHSMEgZIwgY+AFAI4qxGj\r\n" *
            "\tloCLDdMVKwiljjDastqooXSkcjBwMQswCQYDVQQGEwJVSzERMA8GA1UEChMIZVNj\r\n" *
            "\taWVuY2UxEjAQBgNVBAsTCUF1dGhvcml0eTELMAkGA1UEAxMCQ0ExLTArBgkqhkiG\r\n" *
            "\t9w0BCQEWHmNhLW9wZXJhdG9yQGdyaWQtc3VwcG9ydC5hYy51a4IBADApBgNVHRIE\r\n" *
            "\tIjAggR5jYS1vcGVyYXRvckBncmlkLXN1cHBvcnQuYWMudWswGQYDVR0gBBIwEDAO\r\n" *
            "\tBgwrBgEEAdkvAQEBAQYwPQYJYIZIAYb4QgEEBDAWLmh0dHA6Ly9jYS5ncmlkLXN1\r\n" *
            "\tcHBvcnQuYWMudmT4sopwqlBWsvcHViL2NybC9jYWNybC5jcmwwPQYJYIZIAYb4QgEDBDAWLmh0\r\n" *
            "\tdHA6Ly9jYS5ncmlkLXN1cHBvcnQuYWMudWsvcHViL2NybC9jYWNybC5jcmwwPwYD\r\n" *
            "\tVR0fBDgwNjA0oDKgMIYuaHR0cDovL2NhLmdyaWQt5hYy51ay9wdWIv\r\n" *
            "\tY3JsL2NhY3JsLmNybDANBgkqhkiG9w0BAQUFAAOCAQEAS/U4iiooBENGW/Hwmmd3\r\n" *
            "\tXCy6Zrt08YjKCzGNjorT98g8uGsqYjSxv/hmi0qlnlHs+k/3Iobc3LjS5AMYr5L8\r\n" *
            "\tUO7OSkgFFlLHQyC9JzPfmLCAugvzEbyv4Olnsr8hbxF1MbKZoQxUZtMVu29wjfXk\r\n" *
            "\thTeApBv7eaKCWpSp7MCbvgzm74izKhu3vlDk9w6qVrxePfGgpKPqfHiOoGhFnbTK\r\n" *
            "\twTC6o2xq5y0qZ03JonF7OJspEd3I5zKY3E+ov7/ZhW6DqT8UFvsAdjvQbXyhV8Eu\r\n" *
            "\tYhixw1aKEPzNjNowuIseVogKOLXxWI5vAi5HgXdS0/ES5gDGsABo4fqovUKlgop3\r\n" *
            "\tRA==\r\n" *
            "\t-----END CERTIFICATE-----\r\n" *
            "\r\n"

            r = parse(Request, reqstr)
            @test r.method == "GET"
            @test "GET / HTTP/1.1X-SSL-FoooBarr:   $(header(r, "X-SSL-FoooBarr"))" == replace(reqstr, "\r\n" => "")
            @test_throws HTTP.HTTP.ParseError HTTP.parse(Request, "GET / HTTP/1.1\r\nHost: www.example.com\r\nConnection\r\033\065\325eep-Alive\r\nAccept-Encoding: gzip\r\n\r\n")

            r = parse(Request,"GET /bad_get_no_headers_no_body/world HTTP/1.1\r\nAccept: */*\r\n\r\nHELLO")
            @test String(r.body) == ""
        end

        @testset "Response - readheaders() - 1" begin
            respstr = "HTTP/1.1 200 OK\r\n" * "Content-Length: " * "1844674407370955160" * "\r\n\r\n"
            r = Response()
            readheaders(IOBuffer(respstr), r)
            @test r.status == 200
            @test [r.headers...] == ["Content-Length"=>"1844674407370955160"]
        end

        @testset "Response - readheaders() - 2" begin
            respstr = "HTTP/1.1 200 OK\r\n" * "Transfer-Encoding: chunked\r\n\r\n" * "FFFFFFFFFFFFFFE" * "\r\n..."
            r = Response()
            readheaders(IOBuffer(respstr), r)
            @test r.status == 200
            @test [r.headers...] == ["Transfer-Encoding"=>"chunked"]
        end
    end

    @testset "parse - $response" for response in responses
        r = parse(Response, response.raw)
        r_headers = [tocameldash(n) => String(v) for (n,v) in r.headers]

        @test r.version.major == response.http_major
        @test r.version.minor == response.http_minor
        @test r.status == response.status_code
        @test HTTP.Messages.statustext(r.status) == response.response_status
        @test length(r.headers) == response.num_headers
        @test r_headers == response.headers
        @test String(r.body) == response.body

        @test_skip @test HTTP.http_should_keep_alive(HTTP.DEFAULT_PARSER) == response.should_keep_alive
    end

    @testset "Parse Errors" begin
        @testset "Requests" begin
            @testset "ArgumentError - Invalid Base 10 Digit" begin
                reqstr = "GET / HTTP/1.1\r\n" * "Content-Length: 0\r\nContent-Length: 1\r\n\r\n"
                e = try parse(Request, reqstr) catch e e end
                @test isa(e, ArgumentError)
            end

            @testset "EOFError - 1" begin
                reqstr = "GET / HTTP/1.1\r\n" * "Transfer-Encoding: chunked\r\nContent-Length: 1\r\n\r\n"
                e = try parse(Request, reqstr) catch e e end
                @test isa(e, EOFError)
            end

            @testset "EOFError - 2" begin
                reqstr = "GET / HTTP/1.1\r\nheader: value\nhdr: value\r\n"
                e = try parse(Request, reqstr) catch e e end
                @test isa(e, EOFError)
            end

            @testset "Invalid Request Line" begin
                reqstr = "GET / HTP/1.1\r\n\r\n"
                e = try parse(Request, reqstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code ==:INVALID_REQUEST_LINE
            end

            @testset "Invalid Header Field - 1" begin
                reqstr = "GET / HTTP/1.1\r\n" * "Fo@: Failure\r\n\r\n"
                e = try parse(Request, reqstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code ==:INVALID_HEADER_FIELD
            end

            @testset "Invalid Header Field - 2" begin
                reqstr = "GET / HTTP/1.1\r\n" * "Foo\01\test: Bar\r\n\r\n"
                e = try parse(Request, reqstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code ==:INVALID_HEADER_FIELD
            end

            @testset "Invalid Header Field - 3" begin
                reqstr = "GET / HTTP/1.1\r\n" * "Foo: 1\rBar: 1\r\n\r\n"
                e = try parse(Request, reqstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code ==:INVALID_HEADER_FIELD
            end

            @testset "Invalid Header Field - 4" begin
                reqstr = "GET / HTTP/1.1\r\n" * "name\r\n" * " : value\r\n\r\n"
                e = try parse(Request, reqstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code ==:INVALID_HEADER_FIELD
            end

            @testset "Invalid Request Line" begin
                for m in ("HTTP/1.1", "hello world")
                    reqstr = "$m / HTTP/1.1\r\n\r\n"
                    e = try parse(Request, reqstr) catch e e end
                    @test isa(e, HTTP.ParseError) && e.code ==:INVALID_REQUEST_LINE
                end
            end

            @testset "Strict Headers - 1" begin
                reqstr = "GET / HTTP/1.1\r\n" * "Foo: F\01ailure\r\n\r\n"
                strict && @test_throws HTTP.ParseError parse(Request,reqstr)

                if !strict
                    r = HTTP.parse(HTTP.Messages.Request, reqstr)
                    @test r.method == "GET"
                    @test r.target == "/"
                    @test length(r.headers) == 1
                end
            end

            @testset "Strict Headers - 2" begin
                reqstr = "GET / HTTP/1.1\r\n" * "Foo: B\02ar\r\n\r\n"
                strict && @test_throws HTTP.ParseError parse(Request, reqstr)

                if !strict
                    r = parse(HTTP.Messages.Request, reqstr)
                    @test r.method == "GET"
                    @test r.target == "/"
                    @test length(r.headers) == 1
                end
            end
        end

        @testset "Responses" begin
            @testset "ArgumentError - Invalid Base 10 Digit" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Content-Length: 0\r\nContent-Length: 1\r\n\r\n"
                e = try parse(Response, respstr) catch e e end
                @test isa(e, ArgumentError)
            end

            @testset "Chunk Size Exceeds Limit" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Transfer-Encoding: chunked\r\n\r\n" * "FFFFFFFFFFFFFFF" * "\r\n..."
                e = try parse(Response,respstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code == :CHUNK_SIZE_EXCEEDS_LIMIT
            end

            @testset "Chunk Size Exceeds Limit" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Transfer-Encoding: chunked\r\n\r\n" * "10000000000000000" * "\r\n..."
                e = try parse(Response,respstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code == :CHUNK_SIZE_EXCEEDS_LIMIT
            end

            @testset "EOF Error" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Transfer-Encoding: chunked\r\nContent-Length: 1\r\n\r\n"
                e = try parse(Response, respstr) catch e e end
                @test isa(e, EOFError)
            end

            @test_skip @testset "Invalid Content Length" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Content-Length: " * "18446744073709551615" * "\r\n\r\n"
                e = try parse(Response,respstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code == Parsers.HPE_INVALID_CONTENT_LENGTH
            end

            @testset "Invalid Header Field - 1" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Fo@: Failure\r\n\r\n"
                e = try parse(Response, respstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code ==:INVALID_HEADER_FIELD
            end

            @testset "Invalid Header Field - 2" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Foo\01\test: Bar\r\n\r\n"
                e = try parse(Response, respstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code ==:INVALID_HEADER_FIELD
            end

            @testset "Invalid Header Field - 3" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Foo: 1\rBar: 1\r\n\r\n"
                e = try parse(Response, respstr) catch e e end
                @test isa(e, HTTP.ParseError) && e.code ==:INVALID_HEADER_FIELD
            end

            @testset "Strict Headers - 1" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Foo: F\01ailure\r\n\r\n"
                strict && @test_throws HTTP.ParseError parse(Response,respstr)

                if !strict
                    r = parse(HTTP.Messages.Response, respstr)
                    @test r.status == 200
                    @test length(r.headers) == 1
                end
            end

            @testset "Strict Headers - 2" begin
                respstr = "HTTP/1.1 200 OK\r\n" * "Foo: B\02ar\r\n\r\n"
                strict && @test_throws HTTP.ParseError parse(Response,respstr)

                if !strict
                    r = parse(HTTP.Messages.Response, respstr)
                    @test r.status == 200
                    @test length(r.headers) == 1
                end
            end
        end
    end
end


text:
using PyPlot
using JLD
function plot_stable_sens()
	X = load("../data/obj_erg_avg/cos4y_s4.jld")
	s4_arr = X["s4"]
	J_arr = X["J"]
	fig, ax = subplots(1,1)
	ax.plot(s4_arr, J_arr, ".", ms=10.0)
    ax.xaxis.set_tick_params(labelsize=32)
    ax.yaxis.set_tick_params(labelsize=32)
    ax.set_xlabel(L"$s_4$",fontsize=32)
    ax.set_ylabel(L"$\langle J\rangle$",fontsize=32)
	ax.grid(true)

	X = load("../data/stable_sens/dJds4.jld")
	dJds = X["dJds"]
	s4 = X["s4"]

	X = load("../data/obj_erg_avg/cos4y_s4_sens.jld")
	J = X["J"]

	eps = 2.e-2
	n = size(dJds)[1]
	J_pts = reshape([J .- eps*dJds J .+ eps*dJds], n,
					2)'
	s_pts = reshape([s4 .- eps s4 .+ eps], n,
					2)'

	ax.plot(s_pts[:,1:3:end-1], J_pts[:,1:3:end-1], "k",lw=6.0)
	ax.plot(s_pts[:,end],J_pts[:,end],"k",lw=6.0,label=L"$d_s\langle J\rangle$ from S3")
	ax.legend(fontsize=32)

end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
#=
Given a string of round, curly, and square open and closing brackets, return whether the brackets are balanced (well-formed).

For example, given the string "([])[]({})", you should return true.

Given the string "([)]" or "((()", you should return false.
=#

using Test
include("Solutions/problem27_are_brackets_balanced.jl")

@test isbalanced("([])[]({})") == true
@test isbalanced("([)]") == false
@test isbalanced("((()") == false

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
function SymmLQ(arg0::IterationManager, arg1::jdouble, arg2::jboolean)
    return SymmLQ((IterationManager, jdouble, jboolean), arg0, arg1, arg2)
end

function SymmLQ(arg0::jint, arg1::jdouble, arg2::jboolean)
    return SymmLQ((jint, jdouble, jboolean), arg0, arg1, arg2)
end

function should_check(obj::SymmLQ)
    return jcall(obj, "shouldCheck", jboolean, ())
end

function solve(obj::SymmLQ, arg0::RealLinearOperator, arg1::RealLinearOperator, arg2::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealLinearOperator, RealVector), arg0, arg1, arg2)
end

function solve(obj::SymmLQ, arg0::RealLinearOperator, arg1::RealLinearOperator, arg2::RealVector, arg3::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2, arg3)
end

function solve(obj::SymmLQ, arg0::RealLinearOperator, arg1::RealLinearOperator, arg2::RealVector, arg3::jboolean, arg4::jdouble)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealLinearOperator, RealVector, jboolean, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function solve(obj::SymmLQ, arg0::RealLinearOperator, arg1::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealVector), arg0, arg1)
end

function solve(obj::SymmLQ, arg0::RealLinearOperator, arg1::RealVector, arg2::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2)
end

function solve(obj::SymmLQ, arg0::RealLinearOperator, arg1::RealVector, arg2::jboolean, arg3::jdouble)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealVector, jboolean, jdouble), arg0, arg1, arg2, arg3)
end

function solve_in_place(obj::SymmLQ, arg0::RealLinearOperator, arg1::RealLinearOperator, arg2::RealVector, arg3::RealVector)
    return jcall(obj, "solveInPlace", RealVector, (RealLinearOperator, RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2, arg3)
end

function solve_in_place(obj::SymmLQ, arg0::RealLinearOperator, arg1::RealLinearOperator, arg2::RealVector, arg3::RealVector, arg4::jboolean, arg5::jdouble)
    return jcall(obj, "solveInPlace", RealVector, (RealLinearOperator, RealLinearOperator, RealVector, RealVector, jboolean, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function solve_in_place(obj::SymmLQ, arg0::RealLinearOperator, arg1::RealVector, arg2::RealVector)
    return jcall(obj, "solveInPlace", RealVector, (RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2)
end


text:
import DataArrays.skipmissing

##### To define functions for calculating the distances:

function calc_distance(ss1::ExoplanetsSysSim.CatalogSummaryStatistics, ss2::ExoplanetsSysSim.CatalogSummaryStatistics, return_KS_or_AD::String ; AD_mod::Bool=false, all_dist::Bool=false, save_dist::Bool=true)
    #This function calculates the total KS distance between two summary statistics (simulated observed catalogs).
    #If 'all_dist=true', the function outputs the individual distances in the distance function.
    #If 'save_dist=true', the function also saves the distances (individual and total) to a file (assuming file 'f' is open for writing).

    M_cat_obs1 = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system for ss1
    M_cat_obs2 = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system for ss2
    for k in 1:get_int(sim_param,"max_tranets_in_sys")
        append!(M_cat_obs1, k*ones(Int64, ss1.stat["num n-tranet systems"][k]))
        append!(M_cat_obs2, k*ones(Int64, ss2.stat["num n-tranet systems"][k]))
    end

    #To handle empty arrays:
    if sum(M_cat_obs1 .>= 2) < 2 || sum(M_cat_obs2 .>= 2) < 2 #need at least 2 multi-systems per catalog in order to be able to compute AD distances for distributions of ratios of observables
        println("Not enough observed multi-planet systems in one of the catalogs to compute the AD distance.")
        d = ones(Int64,8)*1e6

        println("Distances: ", d, [sum(d)])
        if save_dist
            println(f, "Dist_KS: ", d, [sum(d)])
            println(f, "Dist_AD: ", d, [sum(d)])
        end

        if all_dist
            return d
        else
            return sum(d)
        end
    end

    const max_incl_sys = get_real(sim_param,"max_incl_sys")
    cos_factor = cos(max_incl_sys*pi/180) #factor to divide the number of targets in simulation by to get the actual number of targets needed (with an isotropic distribution of system inclinations) to produce as many transiting systems for a single observer

    #To compute the KS distances:
    d_KS = Array{Float64}(8)
    d_KS[1] = abs(ss1.stat["num_tranets"]/(ss1.stat["num targets"]/cos_factor) - ss2.stat["num_tranets"]/(ss2.stat["num targets"]))
    d_KS[2] = ksstats_ints(M_cat_obs1, M_cat_obs2)[5]
    d_KS[3] = ksstats(ss1.stat["P list"], ss2.stat["P list"])[5]
    d_KS[4] = ksstats(ss1.stat["period_ratio_list"], ss2.stat["period_ratio_list"])[5]
    d_KS[5] = ksstats(ss1.stat["duration list"], ss2.stat["duration list"])[5]
    d_KS[6] = ksstats(ss1.stat["duration_ratio_list"], ss2.stat["duration_ratio_list"])[5]
    d_KS[7] = ksstats(ss1.stat["depth list"], ss2.stat["depth list"])[5]
    d_KS[8] = ksstats(ss1.stat["radius_ratio_list"], ss2.stat["radius_ratio_list"])[5]

    #To compute the AD distances:
    if AD_mod
        ADdist = ADstats_mod
    else
        ADdist = ADstats
    end

    d_AD = Array{Float64}(8)
    d_AD[1] = abs(ss1.stat["num_tranets"]/(ss1.stat["num targets"]/cos_factor) - ss2.stat["num_tranets"]/(ss2.stat["num targets"]))
    d_AD[2] = ksstats_ints(M_cat_obs1, M_cat_obs2)[5]
    d_AD[3] = ADdist(ss1.stat["P list"], ss2.stat["P list"])
    d_AD[4] = ADdist(ss1.stat["period_ratio_list"], ss2.stat["period_ratio_list"])
    d_AD[5] = ADdist(ss1.stat["duration list"], ss2.stat["duration list"])
    d_AD[6] = ADdist(ss1.stat["duration_ratio_list"], ss2.stat["duration_ratio_list"])
    d_AD[7] = ADdist(ss1.stat["depth list"], ss2.stat["depth list"])
    d_AD[8] = ADdist(ss1.stat["radius_ratio_list"], ss2.stat["radius_ratio_list"])

    #To print and/or write the distances to file:
    println("KS Distances: ", d_KS, [sum(d_KS)])
    println("AD Distances: ", d_AD, [sum(d_AD)])
    if save_dist
        println(f, "Dist_KS: ", d_KS, [sum(d_KS)])
        println(f, "Dist_AD: ", d_AD, [sum(d_AD)])
    end

    #To return the distances or total distance:
    if all_dist
        if return_KS_or_AD == "KS"
            return d_KS
        elseif return_KS_or_AD == "AD"
            return d_AD
        end
    else
        if return_KS_or_AD == "KS"
            return sum(d_KS)
        elseif return_KS_or_AD == "AD"
            return sum(d_AD)
        end
    end
end

function calc_distance_Kepler(ss1::ExoplanetsSysSim.CatalogSummaryStatistics, return_KS_or_AD::String ; AD_mod::Bool=false, all_dist::Bool=false, save_dist::Bool=true)
    #This function calculates the total KS distance between a population generated by our model and the actual Kepler population (which must be loaded in).
    #If 'all_dist=true', the function outputs the individual distances in the distance function.
    #If 'save_dist=true', the function also saves the distances (individual and total) to a file (assuming file 'f' is open for writing).

    M_cat_obs = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system
    for k in 1:get_int(sim_param,"max_tranets_in_sys")
        append!(M_cat_obs, k*ones(Int64, ss1.stat["num n-tranet systems"][k]))
    end

    if sum(M_cat_obs .>= 2) < 2 #need at least 2 observed multi-systems in order to be able to compute AD distances for distributions of ratios of observables
        println("Not enough observed multi-planet systems in the simulated catalog.")
        d = ones(Int64,8)*1e6

        println("Distances: ", d, [sum(d)])
        if save_dist
            println(f, "Dist_KS: ", d, [sum(d)])
            println(f, "Dist_AD: ", d, [sum(d)])
        end

        if all_dist
            return d
        else
            return sum(d)
        end
    end

    const max_incl_sys = get_real(sim_param,"max_incl_sys")
    cos_factor = cos(max_incl_sys*pi/180) #factor to divide the number of targets in simulation by to get the actual number of targets needed (with an isotropic distribution of system inclinations) to produce as many transiting systems for a single observer

    #To compute the KS distances:
    d_KS = Array{Float64}(8)
    d_KS[1] = abs(ss1.stat["num_tranets"]/(ss1.stat["num targets"]/cos_factor) - length(P_confirmed)/N_Kepler_targets)
    d_KS[2] = ksstats_ints(M_cat_obs, M_confirmed)[5]
    d_KS[3] = ksstats(ss1.stat["P list"], P_confirmed)[5]
    d_KS[4] = ksstats(ss1.stat["period_ratio_list"], R_confirmed)[5]
    d_KS[5] = ksstats(ss1.stat["duration list"].*24, t_D_confirmed)[5] #transit durations in simulations are in days, while in the Kepler catalog are in hours
    d_KS[6] = ksstats(ss1.stat["duration_ratio_list"], xi_confirmed)[5]
    d_KS[7] = ksstats(ss1.stat["depth list"], D_confirmed)[5]
    d_KS[8] = ksstats(ss1.stat["radius_ratio_list"].^2, D_ratio_confirmed)[5] #simulations save radius ratios while we computed transit duration ratios from the Kepler catalog

    #To compute the AD distances:
    if AD_mod
        ADdist = ADstats_mod
    else
        ADdist = ADstats
    end

    d_AD = Array{Float64}(8)
    d_AD[1] = abs(ss1.stat["num_tranets"]/(ss1.stat["num targets"]/cos_factor) - length(P_confirmed)/N_Kepler_targets)
    d_AD[2] = ksstats_ints(M_cat_obs, M_confirmed)[5]
    d_AD[3] = ADdist(ss1.stat["P list"], P_confirmed)
    d_AD[4] = ADdist(ss1.stat["period_ratio_list"], R_confirmed)
    d_AD[5] = ADdist(ss1.stat["duration list"].*24, t_D_confirmed) #transit durations in simulations are in days, while in the Kepler catalog are in hours
    d_AD[6] = ADdist(ss1.stat["duration_ratio_list"], xi_confirmed)
    d_AD[7] = ADdist(ss1.stat["depth list"], D_confirmed)
    d_AD[8] = ADdist(ss1.stat["radius_ratio_list"].^2, D_ratio_confirmed) #simulations save radius ratios while we computed transit duration ratios from the Kepler catalog

    #To print and/or write the distances to file:
    println("KS Distances: ", d_KS, [sum(d_KS)])
    println("AD Distances: ", d_AD, [sum(d_AD)])
    if save_dist
        println(f, "Dist_KS: ", d_KS, [sum(d_KS)])
        println(f, "Dist_AD: ", d_AD, [sum(d_AD)])
    end

    #To return the distances or total distance:
    if all_dist
        if return_KS_or_AD == "KS"
            return d_KS
        elseif return_KS_or_AD == "AD"
            return d_AD
        end
    else
        if return_KS_or_AD == "KS"
            return sum(d_KS)
        elseif return_KS_or_AD == "AD"
            return sum(d_AD)
        end
    end
end

function target_function(active_param::Vector{Float64}, use_KS_or_AD::String, Kep_or_Sim::String ; AD_mod::Bool=false, weights::Vector{Float64}=ones(8), all_dist::Bool=false, save_dist::Bool=true)
    #This function takes in the values of the active model parameters, generates a simulated observed catalog, and computes the distance function.
    #If 'all_dist=true', the function outputs the individual distances in the distance function.
    #If 'save_dist=true', the function also saves the distances (unweighted and weighted, individual and total) to a file (assuming file 'f' is open for writing).

    println("Active parameter values: ", active_param)
    if save_dist
        println(f, "Active_params: ", active_param) #to write the params to file
    end

    global sim_param
    sim_param_here = deepcopy(sim_param)
    ExoplanetsSysSim.update_sim_param_from_vector!(active_param,sim_param_here)
    cat_phys = generate_kepler_physical_catalog(sim_param_here)
    cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param_here)
    cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param_here)
    summary_stat = calc_summary_stats_model(cat_obs,sim_param_here)

    if Kep_or_Sim == "Kep"
        dist = calc_distance_Kepler(summary_stat, use_KS_or_AD; AD_mod=AD_mod, all_dist=true, save_dist=save_dist)
    elseif Kep_or_Sim == "Sim"
        global summary_stat_ref
        dist = calc_distance(summary_stat, summary_stat_ref, use_KS_or_AD; AD_mod=AD_mod, all_dist=true, save_dist=save_dist)
    end

    weighted_dist = dist./weights
    #used_dist = weighted_dist[1:6] #choose a subset of the distances to pass into the optimizer

    if weights != ones(8)
        println("Weighted distances: ", weighted_dist, [sum(weighted_dist)])
        if save_dist
            println(f, "Dist_weighted: ", weighted_dist, [sum(weighted_dist)])
        end
    end

    if all_dist
        return weighted_dist
    else
        return sum(weighted_dist)
    end
end

function compute_weights_target_fitness_std_perfect_model(num_evals::Int64, use_KS_or_AD::String ; AD_mod::Bool=false, weight::Bool=true, save_dist::Bool=true)
    tic()
    active_param_true = make_vector_of_sim_param(sim_param)
    println("# True values: ", active_param_true)
    if save_dist
        println(f, "# Format: Dist: [distances][total distance]")
    end

    dists_true = zeros(num_evals,8)
    for i in 1:num_evals
        dists_true[i,:] = target_function(active_param_true, use_KS_or_AD, "Sim"; AD_mod=AD_mod, all_dist=true, save_dist=save_dist)
    end
    mean_dists = transpose(mean(dists_true,1))[:,] #array of mean distances for each individual distance
    mean_dist = mean(sum(dists_true, 2)) #mean total distance
    #std_dists = transpose(std(dists_true, 1))[:,] #array of std distances for each individual distance
    rms_dists = transpose(sqrt.(mean(dists_true.^2, 1)))[:,] #array of rms (std around 0)  distances for each individual distance
    std_dist = std(sum(dists_true, 2)) #std of total distance
    rms_dist = sqrt(mean(sum(dists_true, 2).^2)) #rms (std around 0) of total distance; should be similar to the mean total distance

    if weight
        weights = rms_dists #to use the array 'rms_dists' as the weights for the individual distances
    else
        weights = ones(8)
    end

    weighted_dists_true = zeros(num_evals,8)
    for i in 1:num_evals
        weighted_dists_true[i,:] = dists_true[i,:]./weights
    end
    mean_weighted_dists = transpose(mean(weighted_dists_true,1))[:,] #array of mean weighted distances for each individual distance
    mean_weighted_dist = mean(sum(weighted_dists_true,2)) #mean weighted total distance
    std_weighted_dist = std(sum(weighted_dists_true,2)) #std of weighted total distance
    t_elapsed = toc()

    println("Mean dists: ", mean_dists)
    println("Rms dists: ", rms_dists)
    println("Mean weighted dists: ", mean_weighted_dists)
    println("Distance using true values: ", mean_dist, " +/- ", std_dist)
    println("Weighted distance using true values: ", mean_weighted_dist, " +/- ", std_weighted_dist)
    if save_dist
        println(f, "Mean: ", mean_dists, [mean_dist])
        println(f, "Rms: ", rms_dists, [rms_dist])
        println(f, "Mean weighted dists: ", mean_weighted_dists, [mean_weighted_dist])
        println(f, "# Distance using true values (default parameter values): ", mean_dist, " +/- ", std_dist)
        println(f, "# Weighted distance using true values (default parameter values): ", mean_weighted_dist, " +/- ", std_weighted_dist)
        println(f, "# elapsed time: ", t_elapsed, " seconds")
        println(f, "#")
    end

    return (active_param_true, weights, mean_weighted_dist, std_weighted_dist)
end


text:
doc"""
    VonMises(μ, κ)

The *von Mises distribution* with mean `μ` and concentration `κ` has probability density function

$f(x; \mu, \kappa) = \frac{1}{2 \pi I_0(\kappa)} \exp \left( \kappa \cos (x - \mu) \right)$

```julia
VonMises()       # von Mises distribution with zero mean and unit concentration
VonMises(κ)      # von Mises distribution with zero mean and concentration κ
VonMises(μ, κ)   # von Mises distribution with mean μ and concentration κ
```

External links

* [von Mises distribution on Wikipedia](http://en.wikipedia.org/wiki/Von_Mises_distribution)

"""
immutable VonMises <: ContinuousUnivariateDistribution
    μ::Float64      # mean
    κ::Float64      # concentration
    I0κ::Float64    # I0(κ), where I0 is the modified Bessel function of order 0

    function VonMises(μ::Real, κ::Real)
        @check_args(VonMises, κ > zero(κ))
        new(μ, κ, besseli(0, κ))
    end
    VonMises(κ::Real) = VonMises(0.0, κ)
    VonMises() = new(0.0, 1.0)
end

show(io::IO, d::VonMises) = show(io, d, (:μ, :κ))

@distr_support VonMises d.μ - π d.μ + π

#### Parameters

params(d::VonMises) = (d.μ, d.κ)

#### Statistics

mean(d::VonMises) = d.μ
median(d::VonMises) = d.μ
mode(d::VonMises) = d.μ
circvar(d::VonMises) = 1.0 - besseli(1, d.κ) / d.I0κ
entropy(d::VonMises) = log(twoπ * d.I0κ) - d.κ * (besseli(1, d.κ) / d.I0κ)

cf(d::VonMises, t::Real) = (besseli(abs(t), d.κ) / d.I0κ) * cis(t * d.μ)

#### Evaluation

pdf(d::VonMises, x::Float64) = exp(d.κ * cos(x - d.μ)) / (twoπ * d.I0κ)
logpdf(d::VonMises, x::Float64) = d.κ * cos(x - d.μ) - log(d.I0κ) - log2π

cdf(d::VonMises, x::Float64) = _vmcdf(d.κ, d.I0κ, x - d.μ, 1.0e-15)

function _vmcdf(κ::Float64, I0κ::Float64, x::Float64, tol::Float64)
    j = 1
    cj = besseli(j, κ) / j
    s = cj * sin(j * x)
    while abs(cj) > tol
        j += 1
        cj = besseli(j, κ) / j
        s += cj * sin(j * x)
    end
    return (x + 2.0 * s / I0κ) / twoπ + 0.5
end

#### Sampling

sampler(d::VonMises) = VonMisesSampler(d.μ, d.κ)


text:
# Full Training Test (needs GPU for reasonable training time)using MuZero
using MuZero
import Zygote
using Grids
include("game.jl")
include("inference_tests.jl")
using Base.Threads: @spawn, @threads
config = getconfig(Tac)
config.batch_size = 16
config.root_dirichlet_alpha=0.3
config.num_actors = 3
config.weight_decay = 1e-4
config.lr_init=0.002f0
config.save_interval = 10000
config.num_simulations = 100
config.num_unroll_steps = 1
losses, dynamics_stats = [], []

policy_perfs = [];
p_tests = getpolicytests(Tac)
p_tests.period = 1000
testfn = gettestfn(p_tests)
function calc_policy_perf(l)
    acc, _ = testfn(l.test)
    push!(policy_perfs, acc)
end

test_cb = LearnerCallback(
    "test_cb",
    calc_policy_perf,
    period=1000,
)

learner = Learner(
    Tac;
    config=config,
    dir="output",
    cb=[test_cb],
    tests=[p_tests]
)

#selfplay!(learner; limit=config.batch_size)
# @spawn run_selfplay!(config, storage, replay_buffer; limit=33)
# @spawn run_selfplay!(config, storage, replay_buffer; limit=33)
# @spawn run_selfplay!(config, storage, replay_buffer; limit=34)
# rand(1:length(replay_buffer))

for i = 1:(config.num_actors)
    @spawn selfplay!(learner)
end
learner.replay_buffer |> length

#cb = [(1000, 1, [calc_policy_perf]),]

#config.training_steps = 3
#train_network!(learner, losses=losses, callbacks=cb)

config.training_steps = 200001
@show @time fit!(learner; steps=3)
@spawn fit!(learner)
#@spawn train_network!(learner, losses=losses, callbacks=cb)
#@spawn train_network!(config, storage, replay_buffer, losses=losses, callbacks=cb)

@show 1
learner.stats.losses |> length
@show 2
getpct.(policy_perfs)

using MuZero
using Base.Threads: @spawn

include("game.jl")
include("inference_tests.jl")
cd(@__DIR__)
config = getconfig(Tac)
config.checkpoint_interval = 1
config.save_interval = 1000
config.batch_size = 10
config.num_filters = 3
config.representation_resnet_layers = 1
config.dynamics_resnet_layers = 1
config.reward_fc_layers = 1
config.policy_fc_layers = 1
config.value_fc_layers = 1

pv_tests = get_policyvalue_tests(Tac)
pv_tests.period = 1
test_cb = LearnerCallback(
    "test_cb",
    (l)->println("hello from test_cb")
)

learner = Learner(
    Tac;
    dir="output",
    config=config,
    tests=[pv_tests],
    cb=[test_cb]
)

@spawn selfplay!(learner)
length(learner.replay_buffer)
fit!(learner, steps=10)

learner.stats.losses |> length
learner |> gettraindfall |> display
learner |> gettraindfall |> mean100 |> display
_, _, mse = gettestfn(pv_tests)(training_network(learner)) ; mse
_, _, mse = gettestfn(pv_tests)(test_network(learner)) ; mse

learner |> gettraindf |> display
learner |> gettraindfall |> mean100 |> display

using MuZero
include("game.jl")

show_policy_tests(learner)
policy_tests(learner)[3246] |> svg

g = Tac()
g |> display

step!(g, 4)
g |> display

legal_actions(g)

step!(g, 3)
g |> svg

step!(g, [7,9])
g |> svg

s = g |> get_state

isvalid(s)

ps = piece.(Tac, [:X, :O, :X, :O, :X, :O, :X, :O, :X])
s = TacState(ps, PLAYER_2)

ps2 = [piece.(Tac, [:X, :O, :X, :O, :X, :O, :X, :O])..., nothing]
s2 = TacState(ps2, PLAYER_2)

sp = getstatesplus(Tac) ; length(sp)
all_states = states(Tac) ; length(all_states)

getpolicytests(Tac)

target_s = "N0pMDQQAAAA0EAEIVGFjU3RhdGUfTgEETWFpbkQ0EAEKQm9hcmRTdGF0ZR8LnWwJrspjS75KQULJfWQ93gEGTXVaZXJvRAEAAADoFAkDAQMAAwIDAAMAAwADAAMAAwDg"
target_s |> decode |> display

s_plus = getstatesplus(Tac)

for i in 11:82
    @show i
    s_plus[i] |> display
end


text:
import Base.show, Base.==, Base.ctranspose, Base.*
import Base.getindex

export SimpleGraph, IntGraph, StringGraph
export show, NV, NE, has, vertex_type, fastN!
export vlist, elist, neighbors, getindex, deg, deg_hist

"""
The `SimpleGraph` type represents a simple graph; that is, an
undirected graph with no loops and no multiple edges.

Use `SimpleGraph()` to create a new graph in which the vertices may be
`Any` type. Use `SimpleGraph{T}()` to create a new graph in which the
vertices are of type `T`. See `IntGraph` and `StringGraph` as special
cases.
"""
type SimpleGraph{T} <: AbstractSimpleGraph
    V::Set{T}          # Vertex set
    E::Set{Tuple{T,T}}      # Edge set
    N::Dict{T,Set{T}}  # Optional neighbor sets
    Nflag::Bool        # Tells if N is used or not (default on)
    function SimpleGraph(Nflag::Bool=true)
        V = Set{T}()
        E = Set{Tuple{T,T}}()
        N = Dict{T,Set{T}}()
        G = new(V,E,N,Nflag)
    end
end

function show(io::IO, G::SimpleGraph)
    print(io,"SimpleGraph{$(vertex_type(G))} (n=$(NV(G)), m=$(NE(G)))")
end

# Default constructor uses Any type vertices
SimpleGraph(Nflag::Bool = true) = SimpleGraph{Any}(Nflag)

# A StringGraph has vertices of type ASCIIString.
"""
A `StringGraph` is a `SimpleGraph` whose vertices are of type
`ASCIIString`.

When constructed with `StringGraph()` creates an empty
`SimpleGraph{ASCIIString}`.

When invoked as `StringGraph(file::AbstractString)` opens the named
file for parsing as a graph. Each line of the file should contain one
or two tokens separated by white space. If the line contains a single
token, we add that token as a vertex. If the line contains two (or
more) tokens, then the first two tokens are taken as vertex names and
(assuming they are different) the corresponding edge is created. Any
extra tokens on the line are ignored. Lines that begin with a \# are
ignored.
"""
StringGraph() = SimpleGraph{ASCIIString}()

function StringGraph(file::AbstractString)
    G = StringGraph()
    load!(G,file)
    return G
end

# Helper function for StrinGraph(file), and can be used to add
# vertices and edges to a graph (assuming its vertex type can
# accomodate strings).
function load!(G::SimpleGraph, file::AbstractString)
    f = open(file, "r")
    while(~eof(f))
        line = chomp(readline(f))
        tokens = split(line)

        if (length(tokens) == 0)
            continue
        end

        if (tokens[1][1] == '#')
            continue
        end

        add!(G,tokens[1])
        if length(tokens) > 1
            add!(G,tokens[1],tokens[2])
        end
    end
end

"""
`IntGraph()` creates a new `SimpleGraph` whose vertices are of type
`Int`. Called as `IntGraph(n::Int)` prepopulates the vertex set with
vertices `1:n`.
"""

IntGraph() = SimpleGraph{Int}()

# With a postive integer argument, adds 1:n as vertex set, but no
# edges.
function IntGraph(n::Int)
    G = IntGraph()
    for v=1:n
        add!(G,v)
    end
    return G
end

"""
`vertex_type(G)` returns the data type of the vertices this graph may hold.
For example, if `G=IntGraph()` then this returns `Int64`.`
"""
vertex_type{T}(G::SimpleGraph{T}) = T

    # number of vertices and edges

"""
`NV(G)` returns the number of vertices in `G`.
"""
NV(G::AbstractSimpleGraph) = length(G.V)

"""
`NE(G)` returns the number of edges in `G`.
"""
NE(G::SimpleGraph) = length(G.E)

"""
`has(G,v)` returns `true` iff `v` is a vertex of `G`.

`has(G,v,w)` returns `true` iff `(v,w)` is an edge of `G`.
"""
has(G::AbstractSimpleGraph, v) = in(v,G.V)
has(G::SimpleGraph, v, w) = in((v,w), G.E) || in((w,v), G.E)

# fastN(G,true) creates an additional data structure to speed up
# neighborhood lookups.

"""
`fastN!(G,flg=true)` is used to turn on (or off) fast neighborhood
lookup in graphs. Switching this off decreases the size of the data
structure holding the graph, but slows down look up of edges.

**Note**: Fast neighborhood look up is on by default.
"""
function fastN!{T}(G::SimpleGraph{T},flg::Bool=true)
    # if no change, do nothing
    if flg == G.Nflag
        return flg
    end

    # if setting the flag to false, erase the G.N data structure;
    # otherwise build it.
    if flg == false
        G.N = Dict{T, Set{T}}()
        sizehint!(G.N, NV(G))
    else
        # build the G.N structure.
        # start with empty sets for each vertex
        for v in G.V
            G.N[v] = Set{T}()
        end
        # now iterate over the edge set and populate G.N sets
        for e in G.E
            v,w = e
            push!(G.N[v],w)
            push!(G.N[w],v)
        end
    end
    G.Nflag = flg
    return flg
end

# Create a mapping between G.V and 1:n. This is not exposed outside
# this module; it's a helper function used by other functions. This
# has been crafted to work with either SimpleGraph or SimpleDigraph
# arguments.
function vertex2idx(G::AbstractSimpleGraph)
    T = vertex_type(G)
    d = Dict{T,Int}()
    V = vlist(G)
    n = NV(G)

    for k=1:n
        d[V[k]] = k
    end

    return d
end

    # get the vertices as a (sorted if possible) list
"""
`vlist(G)` returns the vertices of `G` as a list (array).
"""
function vlist(G::AbstractSimpleGraph)
    result = collect(G.V)
    try
        sort!(result)
    end
    return result
end

"""
`elist(G)` returns the edges of `G` as a list (array).
"""
function elist(G::SimpleGraph)
    result = collect(G.E)
    try
        sort!(result)
    end
    return result
end

    # Get the neighbors of a vertex
"""
`neighbors(G,v)` returns a list of the neighbors of `v`.

May also be invoked as `G[v]`.
"""
function neighbors{T}(G::SimpleGraph{T}, v)
    if ~has(G,v)
        error("Graph does not contain requested vertex")
    end

    if G.Nflag
        N = collect(G.N[v])

    else
        N = T[]
        for w in G.V
            if has(G,v,w)
                append!(N,[w])
            end
        end
    end
    try
        sort!(N)
    end
    return N
end

# Here is another way to access the neighbors of a vertex: G[v]
getindex(G::SimpleGraph,v) = neighbors(G,v)

# And here's a getindex way to check for edges: G[u,v] is a shortcut
# for has(G,u,v).
getindex(G::SimpleGraph,v,w) = has(G,v,w)

# Degree of a vertex
"""
`deg(G,v)` gives the degree of `v` in `G`.

`deg(G)` gives the degree sequence (sorted).
"""
function deg(G::SimpleGraph, v)
    if ~has(G,v)
        error("Graph does not contain requested vertex")
    end
    if G.Nflag
        return length(G.N[v])
    end
    return length(G[v])
end

# Degree sequence
function deg{T}(G::SimpleGraph{T})
    if G.Nflag
        ds = [deg(G,v) for v in G.V]
    else
        dd = Dict{T,Int}()
        for v in G.V
            dd[v] = 0
        end
        for e in G.E
            v,w = e
            dd[v] += 1
            dd[w] += 1
        end
        ds = collect(values(dd))
    end
    sort!(ds, lt = >)
    return ds
end

# Report how many vertices we have each possible degree.
# If G has n vertices, this returns an n-long vector whose
    # k'th entry is the number of vertices of degree k-1.

"""
`deg_hist(G)` gives a tally of how many vertices of each degree are
present in the graph. Because Julia arrays are 1-based, the indexing
is a bit off. Specifically, entry `k` in the returned array is the
number of vertices of degree `k-1`.
"""
function deg_hist{T}(G::SimpleGraph{T})
    n = NV(G)
    degs = deg(G)
    result = zeros(Int,n)
    for d in degs
        result[d+1] += 1
    end
    return result
end

import Base.hash

function hash(G::SimpleGraph, h::UInt64 = UInt64(0))
    return hash(G.V,h) + hash(G.E,h)
end


text:
struct CohModel <: TimeSeriesModel{2,Real} end
function WhittleLikelihoodInference.add_sdf!(out, ::CohModel, ω)
    out[1] = 4
    out[2] = 3+4im
    out[3] = 9
end
@testset "coherancy" begin
    @test coherancy(CohModel(), 1.0) ≈ [1.0 (1/2-2im/3); (1/2+2im/3) 1.0]
    @test coherance(CohModel(), 1.0) ≈ [1.0 5/6; 5/6 1.0]
    @test groupdelay(CohModel(), 1.0) ≈ [0 -0.9272952180016122; 0.9272952180016122 0]
    @test_throws ErrorException coherancy(OU(1.0,1.0), 1.0)
    @test_throws ErrorException coherance(OU(1.0,1.0), 1.0)
    @test_throws ErrorException groupdelay(OU(1.0,1.0), 1.0)
end


text:
using LinearAlgebra

using Keldysh, HDF5

parse_param(::Type{T}, s::AbstractString) where T = parse(T, s)
parse_param(::Type{String}, s::AbstractString) = string(s)
function parse_params(args, param_def)
  N = length(param_def)

  names = ntuple(N) do i
    param_def[i][2]
  end

  vals = map(param_def) do (type, name, default)
    i = findfirst(isequal("--$name"), args)
    return isnothing(i) ? default : parse_param(type, args[i+1])
  end

  return NamedTuple{names}(vals)
end

@enum SpinEnum spin_up = UInt8(1) spin_down=UInt8(2)
Base.to_index(A, sp::SpinEnum) = Int(sp)
flip(sp::SpinEnum) = sp == spin_up ? spin_down : spin_up

# type representing state of anderson impurity
# occupation of up/down stored in first two bits
struct FockState
    state::UInt8
    function FockState(s)
        @assert s < 4
        new(s)
    end
end

# check whether ith spin is occupied
# note sp ∈ [1, 2] so don't need to do (1 << (sp - 1))
Base.getindex(st::FockState, sp::SpinEnum) = (st.state & UInt8(sp)) > 0

# flip occupation of ith spin
# note sp ∈ [1, 2] so don't need to do (sp << (i - 1))
flip(st::FockState, sp::SpinEnum) = FockState(st.state ⊻ UInt8(sp))

# convert internal state to index
Base.to_index(st::FockState) = Int(st.state + 1)

"""
Compute the populations (i.e. the diagonal components of the impurity density matrix) from a propagator
"""
function populations(P)
  nstates = length(P)
  ξ = [1.0, -1.0, -1.0, 1.0]
  P_lsr_diag = reduce(hcat, (1.0im * ξ[s] * diag(P[s][:lesser]) for s in 1:nstates))
  Zt = sum(P_lsr_diag, dims=2)
  ρt = P_lsr_diag ./ Zt
  return ρt, Zt
end

struct NCAParams
  dyson_rtol::Float64
  dyson_atol::Float64
  dyson_max_iter::Int
  max_order::Int
  function NCAParams(dyson_rtol, dyson_atol, dyson_max_iter, max_order)
    @assert 1 <= max_order <= 2
    new(dyson_rtol, dyson_atol, dyson_max_iter, max_order)
  end
end
NCAParams(; dyson_rtol = 1e-6, dyson_atol = 1e-10, dyson_max_iter = 100, max_order = 1) = NCAParams(dyson_rtol, dyson_atol, dyson_max_iter, max_order)

struct NCAData{T <: AbstractTimeGF, U <: AbstractTimeGrid}
  P0::Array{T,1} # bare propagator
  Δ::Array{T, 1} # hybridization function

  P::Array{T,1} # dressed propagator
  Σ::Array{T,1} # self-energy
  ΣxP::Array{T,1} # self-energy convolved with propagator
  G::Array{T,1} # green's function

  grid::U # time grid
  states::NTuple{4, FockState}
  spins::Tuple{SpinEnum, SpinEnum}
end

function NCAData(P0, Δ)
  states = ntuple(i -> FockState(i-1), 4)
  spins = instances(SpinEnum)

  statesize = length(states)
  indexsize = length(spins)

  @assert length(P0) == statesize
  @assert length(Δ) == indexsize

  grid = first(P0).grid

  X = P0[1]
  P = [zero(X) for _ in 1:statesize]
  Σ = [zero(X) for _ in 1:statesize]
  ΣxP = [zero(X) for _ in 1:statesize]
  G = [zero(X) for _ in 1:indexsize]

  NCAData(P0, Δ, P, Σ, ΣxP, G, grid, states, spins)
end

function Σnca(data::NCAData, t1::TimeGridPoint, t2::TimeGridPoint, st_sigma::FockState)
  sum(data.spins) do sp
    st_prop = flip(st_sigma, sp)
    1.0im * data.P[st_prop][t1, t2] * (st_sigma[sp] ? data.Δ[sp][t1, t2] : -data.Δ[sp][t2, t1, false])
  end
end

# collect all one-crossing terms
# define 4 points t1 > t2 > t3 > t4
#           _________
#           |        |
#      __________    |
# _____|    |    |___|___
#     t1   t2    t3  t4
# <------------------------
# t
# st_sigma <- st0 <- st1 <- st2 <- st_sigma
#
# delta[sp0](t1,t3) * delta[sp1](t2, t4)
function Σoca(data::NCAData, t1::TimeGridPoint, t4::TimeGridPoint, st_sigma::FockState)
  Δ, P, grid = data.Δ, data.P, data.grid

  sum(data.spins) do sp1
    sp0 = flip(sp1)

    st2 = flip(st_sigma, sp1)
    st1 = flip(st2, sp0)
    st0 = flip(st1, sp1)

    # delta[sp1](t2, t4)
    h1 = t2 -> 1.0im * (st_sigma[sp1] ? Δ[sp1][t2, t4] : -Δ[sp1][t4, t2, false])

    # delta[sp0](t1, t3)
    h0 = t3 -> 1.0im * (st_sigma[sp0] ? Δ[sp0][t1, t3] : -Δ[sp0][t3, t1, false])

    # integrate over t3
    f = t2 -> h1(t2) * integrate(t3 -> h0(t3) * P[st1][t2, t3] * P[st2][t3, t4], grid, t2, t4)

    # integrate over t2
    return integrate(t2 -> P[st0][t1, t2] * f(t2), grid, t1, t4)
  end
end

function dyson!(data::NCAData, t1::TimeGridPoint, t2::TimeGridPoint, params::NCAParams)
  @assert t1.cidx >= t2.cidx

  p_t1t2_cur = zeros(ComplexF64, length(data.states))
  p_t1t2_next = zeros(ComplexF64, length(data.states))

  for st in data.states
    p_t1t2_cur[st] = data.P0[st][t1,t2]
    data.P[st][t1,t2] = data.P0[st][t1,t2] # initial guess
  end

  ↻ = (A, B) -> integrate(t -> @inbounds(A[t1, t] * B[t, t2]), data.grid, t1, t2)

  done = false
  iter = 1
  diff = 0.0
  while iter <= params.dyson_max_iter && !done

    for st in data.states
      p_t1t2_next[st] = 0.0

      data.Σ[st][t1, t2] = 0.0
      params.max_order > 0 && (data.Σ[st][t1, t2] += Σnca(data, t1, t2, st))
      params.max_order > 1 && (data.Σ[st][t1, t2] += Σoca(data, t1, t2, st))

      # p = p₀ + p₀ ↻ Σ ↻ p
      data.ΣxP[st][t1, t2] = data.Σ[st] ↻ data.P[st]
      p_t1t2_next[st] += data.P0[st] ↻ data.ΣxP[st]

      p_t1t2_next[st] += data.P0[st][t1, t2]
    end

    diff = norm(p_t1t2_cur - p_t1t2_next)
    done = diff < max(params.dyson_atol, params.dyson_rtol * norm(p_t1t2_cur))
    for st in data.states
      data.P[st][t1,t2] = p_t1t2_next[st]
    end
    p_t1t2_cur .= p_t1t2_next
    iter += 1
  end
end

function nca!(data::NCAData, params::NCAParams)
  N = length(data.grid)
  for d in 0:(N-1) # solve diagonal by diagonal
    println("diagonal $(d+1)/$N")
    for j in 1:(N-d)
      i = j + d
      t1 = data.grid[i]
      t2 = data.grid[j]
      dyson!(data, t1, t2, params)
    end
  end

  return data
end

function make_bare_prop(grid::KeldyshTimeGrid, ρ, ϵ, U)
  E = [0.0, ϵ, ϵ, 2*ϵ + U]
  ξ = [1.0, -1.0, -1.0, 1.0]
  P = map(1:4) do s
    GenericTimeGF(grid, 1, true) do t1, t2
      t1.cidx < t2.cidx && return 0.0
      ϕ = integrate(t -> E[s], grid, t1, t2)
      heaviside(t1.bpoint, t2.bpoint) ? -im * exp(-im * ϕ) : -im * ξ[s] * ρ[s] * exp(-im * ϕ)
    end
  end
  return P
end

function make_bare_prop(grid::FullTimeGrid, ϵ, U)
  E = [0.0, ϵ, ϵ, 2*ϵ + U]
  ξ = [1.0, -1.0, -1.0, 1.0]
  P = map(1:4) do s
    GenericTimeGF(grid, 1, true) do t1, t2
      t1.cidx < t2.cidx && return 0.0
      ϕ = integrate(t -> E[s], grid, t1, t2)
      heaviside(t1.bpoint, t2.bpoint) ? -im * exp(-im * ϕ) : -im * ξ[s] * exp(-im * ϕ)
    end
  end
  return P
end

param_def = [(String, :contour, "keldysh"),
             (Float64, :tmax, 10.0),
             (Float64, :beta, 5.0),
             (Int, :nt, 201),
             (Int, :ntau, 101),
             (Float64, :rtol, 1e-6),
             (Float64, :atol, 1e-10),
             (Int, :max_iter, 200),
             (String, :mode, "nca"),
             (Float64, :D, 10.0),
             (Float64, :eps, -3.0),
             (Float64, :U, 8.0),
             (String, :output_file, "output.h5")]

function main()
  println("running nca...")

  p = parse_params(ARGS, param_def)

  for (k,v) in pairs(p)
    println("$k : $v")
  end

  data =
  if p.contour == "keldysh"
    c = twist(KeldyshContour(tmax=p.tmax))
    grid = KeldyshTimeGrid(c, p.nt)

    ρ = [1.0, 0.0, 0.0, 0.0]
    P0 = make_bare_prop(grid, ρ, p.eps, p.U)

    dos = Keldysh.flat_dos(ν=10.0, D=p.D)
    Δ = [GenericTimeGF(dos, p.beta, grid) for s in 1:2]

    NCAData(P0, Δ)
  elseif p.contour == "full"
    c = twist(FullContour(tmax=p.tmax, β=p.beta))
    grid = FullTimeGrid(c, p.nt, p.ntau)
    P0 = make_bare_prop(grid, p.eps, p.U)

    dos = Keldysh.flat_dos(ν=10.0, D=p.D)
    Δ = [GenericTimeGF(dos, grid) for s in 1:2]

    NCAData(P0, Δ)
  else
    error("unknown contour $(p.contour)")
  end

  max_order = Dict("nca"=>1, "oca"=>2)[p.mode]
  params = NCAParams(dyson_rtol = p.rtol, dyson_atol = p.atol, dyson_max_iter = p.max_iter, max_order = max_order)

  nca!(data, params)

  t = collect(realtimes(data.grid))
  ρt, Zt = populations(data.P)

  h5open(p.output_file, "w") do h5f
    for (k,v) in pairs(p)
      write(h5f, "/input/params/$k", v)
    end

    h5f["/output/obs/pop/rho"] = ρt
    h5f["/output/obs/pop/Z"] = Zt
    h5f["/output/obs/pop/t"] = t
  end
end

if abspath(PROGRAM_FILE) == @__FILE__
  main()
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _premat(equinox1::T, equinox2::T, FK4::Bool) where {T<:AbstractFloat}
    t = (equinox2 - equinox1) / 1000
    if FK4
        st = (equinox1 - 1900) / 1000
        #  Compute 3 rotation angles
        a = sec2rad(t*(23042.53 + st*(139.75 + 0.06*st)
                       + t*(30.23 - 0.27*st + 18.0*t)))
        b = sec2rad(t*t*(79.27 + 0.66*st + 0.32*t)) + a
        c = sec2rad(t*(20046.85 - st*(85.33 + 0.37*st)
                       + t*(-42.67 - 0.37*st -41.8*t)))
    else
        st = (equinox1 / 1000) - 2
        # Compute 3 rotation angles
        a = sec2rad(t *(23062.181 + st*(139.656 + 0.0139*st)
                        + t*(30.188 - 0.344*st + 17.998*t)))
        b = sec2rad(t*t*(79.280 + 0.410*st + 0.205*t)) + a
        c = sec2rad(t*(20043.109 - st*(85.33 + 0.217*st)
                       + t*(-42.665 - 0.217*st - 41.833*t)))
    end
    sa, ca = sincos(a)
    sb, cb = sincos(b)
    sc, cc = sincos(c)
    return @SMatrix [ca * cb * cc - sa * sb   -ca * sb - sa * cb * cc   -cb*sc;
                     sa * cb + ca * sb * cc    ca * cb - sa * sb * cc   -sb*sc;
                     ca * sc                  -sa * sc                   cc]
end

"""
    premat(equinox1, equinox2[, FK4=true]) -> precession_matrix

### Purpose ###

Return the precession matrix needed to go from `equinox1` to `equinox2`.

### Explanation ###

This matrix is used by `precess` and `baryvel` to precess astronomical
coordinates.

### Arguments ###

* `equinox1`: original equinox of coordinates.
* `equinox2`: equinox of precessed coordinates.
* `FK4` (optional boolean keyword): if this keyword is set to `true`, the FK4
  (B1950.0) system precession angles are used to compute the precession matrix.
  When it is `false`, the default, use FK5 (J2000.0) precession angles.

### Output ###

A 3×3 matrix, used to precess equatorial rectangular coordinates.

### Example ###

Return the precession matrix from 1950.0 to 1975.0 in the FK4 system

```jldoctest
julia> using AstroLib

julia> premat(1950,1975,FK4=true)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9} with indices SOneTo(3)×SOneTo(3):
 0.999981    -0.00558775  -0.00242909
 0.00558775   0.999984    -6.78691e-6
 0.00242909  -6.78633e-6   0.999997
```

### Method ###

FK4 constants from "Computational Spherical Astronomy" by Taff (1983),
p. 24. (FK4). FK5 constants from "Explanatory Supplement To The Astronomical
Almanac" 1992, page 104 Table 3.211.1
(https://archive.org/details/131123ExplanatorySupplementAstronomicalAlmanac).

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
premat(eq1::Real, eq2::Real; FK4::Bool=false) =
    _premat(promote(float(eq1), float(eq2))..., FK4)


text:
# ---
# title: 1446. Consecutive Characters
# id: problem1446
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/consecutive-characters/description/>
# hidden: true
# ---
#
# Given a string `s`, the power of the string is the maximum length of a non-
# empty substring that contains only one unique character.
#
# Return _the power_  of the string.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "leetcode"
#     Output: 2
#     Explanation: The substring "ee" is of length 2 with the character 'e' only.
#
#
# **Example 2:**
#
#
#
#     Input: s = "abbcccddddeeeeedcba"
#     Output: 5
#     Explanation: The substring "eeeee" is of length 5 with the character 'e' only.
#
#
# **Example 3:**
#
#
#
#     Input: s = "triplepillooooow"
#     Output: 5
#
#
# **Example 4:**
#
#
#
#     Input: s = "hooraaaaaaaaaaay"
#     Output: 11
#
#
# **Example 5:**
#
#
#
#     Input: s = "tourist"
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 500`
#   * `s` contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

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


text:
import Pkg

function maketempdir()::String
    tmpdir = mktempdir()
    atexit(() -> rm(tmpdir; force = true, recursive = true))
    return tmpdir
end

function allowed_intersection_predictmdfull_predictmd()::Vector{String}
    result::Vector{String} = String[
        "Distributed",
        "Pkg"
        ]
    return result
end

function allowed_intersection_predictmdfull_predictmdextra()::Vector{String}
    result::Vector{String} = String[
        "Distributed",
        "Pkg"
        ]
    return result
end

function main()::Nothing
    predictmdfull_dir = pwd()
    predictmdfull_toml_file_name = joinpath(predictmdfull_dir, "Project.toml")
    predictmdfull_toml = Pkg.TOML.parsefile(predictmdfull_toml_file_name)

    predictmd_toml_file_url = "https://raw.githubusercontent.com/bcbi/PredictMD.jl/master/Project.toml"
    predictmdextra_toml_file_url = "https://raw.githubusercontent.com/bcbi/PredictMDExtra.jl/master/Project.toml"
    tmpdir_predictmd = maketempdir()
    tmpdir_predictmdextra = maketempdir()
    predictmd_toml_file_name = joinpath(tmpdir_predictmd, "Project.toml")
    predictmdextra_toml_file_name = joinpath(tmpdir_predictmdextra, "Project.toml")
    Base.download(predictmd_toml_file_url, predictmd_toml_file_name)
    Base.download(predictmdextra_toml_file_url, predictmdextra_toml_file_name)
    predictmd_toml = Pkg.TOML.parsefile(predictmd_toml_file_name)
    predictmdextra_toml = Pkg.TOML.parsefile(predictmdextra_toml_file_name)
    rm(tmpdir_predictmd; force = true, recursive = true)
    rm(tmpdir_predictmdextra; force = true, recursive = true)

    predictmd_deps = sort(unique(keys(predictmd_toml["deps"])))
    predictmdextra_deps = sort(unique(keys(predictmdextra_toml["deps"])))
    predictmdfull_deps = sort(unique(keys(predictmdfull_toml["deps"])))

    intersection_predictmdfull_predictmd = intersect(predictmdfull_deps, predictmd_deps)
    intersection_predictmdfull_predictmdextra = intersect(predictmdfull_deps, predictmdextra_deps)

    @info("Actual intersection of PredictMDFull and PredictMD", intersection_predictmdfull_predictmd, repr(intersection_predictmdfull_predictmd))
    @info("Allowed intersection of PredictMDFull and PredictMD", allowed_intersection_predictmdfull_predictmd(), repr(allowed_intersection_predictmdfull_predictmd()))
    if !issubset(intersection_predictmdfull_predictmd, allowed_intersection_predictmdfull_predictmd())
        error("Intersection of PredictMDFull and PredictMD contains non-allowed packages")
    end

    @info("Actual intersection of PredictMDFull and PredictMDExtra", intersection_predictmdfull_predictmdextra, repr(intersection_predictmdfull_predictmdextra))
    @info("Allowed intersection of PredictMDFull and PredictMDExtra", allowed_intersection_predictmdfull_predictmdextra(), repr(allowed_intersection_predictmdfull_predictmdextra()))
    if !issubset(intersection_predictmdfull_predictmdextra, allowed_intersection_predictmdfull_predictmdextra())
        error("Intersection of PredictMDFull and PredictMDExtra contains non-allowed packages")
    end
end

main()


text:
function characteristic_matrix_ellipse(q0::Float64, q1::Float64, q2::Float64, q3::Float64, r1::Float64, r2::Float64)

	a11::Float64 = 4.0 * r2^2 * (q0 * q3 + q1 * q2)^2 + r1^2 * (q0^2 + q1^2 - q2^2 - q3^2)^2
	a12::Float64 = r2^2 * (2.0 * q0 * q3 + 2.0 * q1 * q2) * (q0^2 - q1^2 + q2^2 - q3^2) - r1^2 * (2.0 * q0 * q3 - 2.0 * q1 * q2) * (q0^2 + q1^2 - q2^2 - q3^2)
	a13::Float64 = r1^2 * (2.0 * q0 * q2 + 2.0 * q1 * q3) * (q0^2 + q1^2 - q2^2 - q3^2) - r2^2 * (2.0 * q0 * q1 - 2.0 * q2 * q3) * (2.0 * q0 * q3 + 2.0 * q1 * q2)
	a21::Float64 = a12
	a22::Float64 = 4.0 * r1^2 * (q0 * q3 - q1 * q2)^2 + r2^2 * (q0^2 - q1^2 + q2^2 - q3^2)^2
	a23::Float64 = - r2^2 * (2.0 * q0 * q1 - 2.0 * q2 * q3) * (q0^2 - q1^2 + q2^2 - q3^2) - r1^2 * (2.0 * q0 * q2 + 2.0 * q1 * q3) * (2.0 * q0 * q3 - 2.0 * q1 * q2)
	a31::Float64 = a13
	a32::Float64 = a23
	a33::Float64 = 4.0 * r1^2 * (q0 * q2 + q1 * q3)^2 + 4.0 * r2^2 * (q0 * q1 - q2 * q3)^2

	return (a11, a12, a13, a21, a22, a23, a31, a32, a33)
end


text:
Base.@kwdef mutable struct Options{T}
    line_search::Symbol=:armijo
    max_iterations::Int=100
    max_dual_updates::Int=10
    min_step_size::T=1.0e-5
    objective_tolerance::T=1.0e-3
    lagrangian_gradient_tolerance::T=1.0e-3
    constraint_tolerance::T=5.0e-3
    constraint_norm::T=Inf
    initial_constraint_penalty::T=1.0
    scaling_penalty::T=10.0
    max_penalty::T=1.0e8
    reset_cache::Bool=false
    verbose=true
end


text:
# Autogenerated wrapper script for WCS_jll for x86_64-w64-mingw32
export libwcs

JLLWrappers.@generate_wrapper_header("WCS")
JLLWrappers.@declare_library_product(libwcs, "libwcs.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libwcs,
        "bin\\libwcs.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module MemTransportTests

using Thrift
using Test

function testmemtransport()
    @testset "memory transport" begin
        t = TMemoryTransport()
        p = TCompactProtocol(t)

        @test open(t) === nothing
        @test close(t) === nothing
        @test flush(t) === nothing
        @test isa(Thrift.rawio(t), IOBuffer)
        @test isopen(t)

        s1 = "Hello World"
        write(p, s1)
        s2 = read(p, typeof(s1))
        @test s2 == s1

        t = TMemoryTransport()
        p = TCompactProtocol(t)
        write(p, s1)
        t = TMemoryTransport(take!(t.buff))
        p = TCompactProtocol(t)
        s2 = read(p, typeof(s1))
        @test s2 == s1
    end
end

testmemtransport()

end # module MemTransportTests


text:
## CacheElement
mutable struct CacheElement{T}
    data::T
    dirty::Bool
    CacheElement(data::T) where {T} = new{T}(data, true)
end

@inline setdirty!(element::CacheElement) = (element.dirty = true; nothing)
@inline isdirty(element::CacheElement) = element.dirty


text:
using Test, BenchmarkTools

using StaticArrays: SVector
using StaticKernels: Kernel
using NNlib: DenseConvDims, conv!
using ImageFiltering: centered, imfilter!, Inner
using LocalFilters: convolve!
import LocalFilters
using OffsetArrays: OffsetArray
using LoopVectorization: @avx

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5
pd = 19

kernel_axes(k) = map(n -> -n÷2:n÷2, size(k))
inner_axes(k, a) = map(axes(a), kernel_axes(k)) do ax, kx
        first(ax) - first(kx) : last(ax) - last(kx)
    end
inner_size(k, a) = length.(inner_axes(k, a))

function bench_linfilt(k, a, ::Val{:StaticKernels})
    print(rpad("  StaticKernels", pd))

    ktuple = Tuple(k)
    @inline function wf(w)
        return sum(SVector(Tuple(w)) .* SVector(ktuple))
    end
    @assert size(k, 1) == size(k, 2)
    K = size(k, 1)
    kern = Kernel{kernel_axes(k)}(wf)
    b1 = similar(a, size(kern, a))
    @btime map!($kern, $b1, $a)
    b1
end

function bench_linfilt(k, a, ::Val{:LoopVectorization})
    print(rpad("  LoopVectorization", pd))

    # from discourse forum:
    # https://discourse.julialang.org/t/ann-statickernels-jl-fast-kernel-operations-on-arrays/37658
    function filter2davx!(out::AbstractMatrix, A::AbstractMatrix, kern)
        @avx for J in CartesianIndices(out)
            tmp = zero(eltype(out))
            for I ∈ CartesianIndices(kern)
                tmp += A[I + J] * kern[I]
            end
            out[J] = tmp
        end
        out
    end

    b5 = similar(a, inner_axes(k, a))
    k5 = reshape(k, kernel_axes(k))
    @btime $filter2davx!($b5, $a, $k5)
    #filter2davx!(b5, a, k5)
    # weird OffsetArray broadcast problem
    b5 = parent(b5)
end

function bench_linfilt(k, a, ::Val{:NNlib})
    print(rpad("  NNlib", pd))

    ar = reshape(a, (size(a)..., 1, 1))
    b2 = similar(a, (inner_size(k, a)..., 1, 1))
    kr = reshape(k, (size(k)..., 1, 1))
    cdims = DenseConvDims(ar, kr; stride=1, padding=0, dilation=1, flipkernel=true)
    @btime conv!($b2, $ar, $kr, $cdims)
    b2
end

function bench_linfilt(k, a, ::Val{:ImageFiltering})
    print(rpad("  ImageFiltering", pd))

    b3 = similar(a, inner_axes(k, a))
    kc = centered(k)
    @btime imfilter!($b3, $a, $kc, Inner())
    # weird OffsetArray broadcast problem
    b3 = parent(b3)
end

function bench_linfilt(k, a, ::Val{:LocalFilters})
    print(rpad("  LocalFilters", pd))

    b4 = similar(a)
    # mathematical convolution
    kr = LocalFilters.Kernel(reverse(reverse(k, dims=1), dims=2))
    @btime convolve!($b4, $a, $kr)
    # no way to only compute inner part
    b4 = b4[inner_axes(k, a)...]
end

# Linear Filtering

for N in (100, 1000), K in (3, 5, 7)
    a = rand(N, N)
    k = rand(K, K)

    println("Array: $(size(a)), Kernel: $(size(k))")

    b = bench_linfilt(k, a, Val(:StaticKernels))

    # disabled: :LocalFilters
    for s in [:LoopVectorization, :NNlib, :ImageFiltering]
        @test b ≈ bench_linfilt(k, a, Val(s))
    end
end