text:
# SimulationResult sould provide an interface for working with properties of a separate particle
# and with physical properties of the whole system.
struct SimulationResult{sType <: NBodySystem}
    solution::DiffEqBase.AbstractTimeseriesSolution
    simulation::NBodySimulation{sType}
end

function Base.show(stream::IO, sr::SimulationResult)
    print(stream, "N: ")
    show(stream, length(sr.simulation.system.bodies))
    println(stream)
    show(stream, sr.simulation)
    print(stream, "Time steps: ")
    show(stream, sr.solution.destats.naccept)
    println(stream)
    print(stream, "t: ", minimum(sr.solution.t), ", ", maximum(sr.solution.t))
    println(stream)
end

(sr::SimulationResult)(args...; kwargs...) = return sr.solution(args...; kwargs...)

# This iterator interface is implemented specifically for making animation.
# Probably, there will be a wrapper for this in the future.

function Base.iterate(sr::SimulationResult,state=1)
  state > length(sr.solution.t) && return nothing
  (sr, sr.solution.t[state]), state + 1
end

function get_velocity(sr::SimulationResult, time::Real, i::Integer=0)
    n = get_coordinate_vector_count(sr.simulation.system)

    if typeof(sr.solution[1]) <: RecursiveArrayTools.ArrayPartition
        velocities = sr(time).x[1]
        if i <= 0
            return velocities[:, 1:n]
        else
            return velocities[:, i]
        end
    else
        velocities = sr(time)
        if sr.simulation.thermostat isa LangevinThermostat
            velocities = reshape(velocities, 3, 2*n)
        end

        if i <= 0
            return @view velocities[:, n + 1:2*n]
        else
            return @view velocities[:, n + i]
        end
    end
end

function get_position(sr::SimulationResult, time::Real, i::Integer=0)
    if typeof(sr.solution[1]) <: RecursiveArrayTools.ArrayPartition
        positions = sr(time).x[2]
    else
        positions = sr(time)
    end

    n = get_coordinate_vector_count(sr.simulation.system)

    if i <= 0
        return @view positions[:, 1:n]
    else
        return @view positions[:, i]
    end
end

function get_masses(system::NBodySystem)
    n = length(system.bodies)
    masses = zeros(n)
    for i = 1:n
        masses[i] = system.bodies[i].m
    end
    return masses
end

function get_masses(system::WaterSPCFw)
    n = length(system.bodies)
    ms = zeros(Real, 3 * n)
    for i = 1:n
        ms[3 * (i - 1) + 1] = system.mO
        ms[3 * (i - 1) + 2] = system.mH
        ms[3 * (i - 1) + 3] = system.mH
    end
    return ms
end

function get_coordinate_vector_count(system::NBodySystem)
    return length(system.bodies)
end

function get_coordinate_vector_count(system::WaterSPCFw)
    return 3* length(system.bodies)
end

# n - number of particles
# nc - number of constraints
# ndf - number of degrees of freedom
function get_degrees_of_freedom(system::NBodySystem)
    n = length(system.bodies)
    nc = 0
    ndf = 3*n-nc
    (n, nc, ndf)
end

function get_degrees_of_freedom(system::WaterSPCFw)
    n = 3*length(system.bodies)
    nc = 2*length(system.bodies)
    ndf = 3*n-nc
    (n, nc, ndf)
end

function temperature(result::SimulationResult, time::Real)
    kb = result.simulation.kb
    (n, nc, ndf) = get_degrees_of_freedom(result.simulation.system)
    vs = get_velocity(result, time)
    ms = get_masses(result.simulation.system)
    return md_temperature(vs, ms, kb, n, nc)
end

function kinetic_energy(velocities, masses)
    ke = sum(dot(vec(sum(velocities.^2, dims=1)), masses / 2))
    return ke
end

function kinetic_energy(sr::SimulationResult, time::Real)
    vs = get_velocity(sr, time)
    ms = get_masses(sr.simulation.system)
    return kinetic_energy(vs, ms)
end

function potential_energy(coordinates, simulation::NBodySimulation)
    e_potential = 0
    system = simulation.system
    n = length(system.bodies)
    if :lennard_jones ∈ keys(system.potentials)
        p = system.potentials[:lennard_jones]
        (ms, indxs) = obtain_data_for_lennard_jones_interaction(system)
        e_potential += lennard_jones_potential(p, indxs, coordinates, simulation.boundary_conditions)
    end

    if :electrostatic ∈ keys(system.potentials)
        p = system.potentials[:electrostatic]
        (qs, ms, indxs, exclude) = obtain_data_for_electrostatic_interaction(simulation.system)
        e_potential += electrostatic_potential(p, indxs, exclude, qs, coordinates, simulation.boundary_conditions)
    end
    e_potential
end

function potential_energy(coordinates, simulation::NBodySimulation{<:WaterSPCFw})
    e_potential = 0
    system = simulation.system
    p = system.lj_parameters
    (ms, indxs) = obtain_data_for_lennard_jones_interaction(system)
    e_potential += lennard_jones_potential(p, indxs, coordinates, simulation.boundary_conditions)

    p = system.e_parameters
    (qs, ms, indxs, exclude) = obtain_data_for_electrostatic_interaction(simulation.system)
    e_potential += electrostatic_potential(p, indxs, exclude, qs, coordinates, simulation.boundary_conditions)

    p = system.scpfw_parameters
    (ms, neighbouhoods) = obtain_data_for_harmonic_bond_interaction(simulation.system, p)
    e_potential += harmonic_bonds_potential(p, coordinates, ms, neighbouhoods)

    p = system.scpfw_parameters
    (ms, bonds) = obtain_data_for_valence_angle_harmonic_interaction(simulation.system)
    e_potential += valence_angle_harmonic_potential(coordinates, bonds)
    e_potential
end

function lennard_jones_potential(p::LennardJonesParameters, indxs::Vector{<:Integer}, coordinates, pbc::BoundaryConditions)
    e_lj = 0
    n = length(indxs)
    @inbounds for ind_i = 1:n
        i = indxs[ind_i]
        ri = @SVector [coordinates[1, i], coordinates[2, i], coordinates[3, i]]
        for ind_j = ind_i + 1:n
            j = indxs[ind_j]
            rj = @SVector [coordinates[1, j], coordinates[2, j], coordinates[3, j]]

            (rij, r, r2) = get_interparticle_distance(ri, rj, pbc)

            if !(r2 < p.R2)
                r2 = p.R2
            end

            σ_rij_6 = (p.σ2 / r2)^3
            σ_rij_12 = σ_rij_6^2
            e_lj += (σ_rij_12 - σ_rij_6 )
        end
    end

    return 4 * p.ϵ * e_lj
end

function electrostatic_potential(p::ElectrostaticParameters, indxs::Vector{<:Integer}, exclude::Dict{Int,Vector{Int}}, qs, rs, pbc::BoundaryConditions)
    e_el = 0.0

    n = length(indxs)
    for ind_i = 1:n
        i = indxs[ind_i]
        ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
        e_el_i = 0.0
        for ind_j = ind_i + 1:n
            j = indxs[ind_j]
            if !in(j, exclude[i])
                rj = @SVector [rs[1, j], rs[2, j], rs[3, j]]

                (rij, r, r2) = get_interparticle_distance(ri, rj, pbc)

                if r2 < p.R2
                    e_el_i += qs[j] / r
                else
                    e_el_i += qs[j] / p.R
                end
            end
        end
        e_el += e_el_i * qs[i]
    end

    return e_el * p.k
end

function harmonic_bonds_potential(p::SPCFwParameters,
    rs,
    ms::Vector{<:Real},
    neighborhoods::Dict{Int,Vector{Tuple{Int,Float64}}})

    e_harmonic = 0

    @inbounds for (i, neighborhood) ∈ neighborhoods
        ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
        for (j, k) in neighborhood
            rj = @SVector [rs[1, j], rs[2, j], rs[3, j]]
            rij = ri - rj
            r = norm(rij)
            d = r - p.rOH
            e_harmonic += d^2 * k
        end
    end
    return e_harmonic / 4
end

function valence_angle_harmonic_potential(rs,
    bonds::Vector{Tuple{Int,Int,Int,Float64,Float64}})

    e_valence = 0

    for (a, b, c, valence_angle, k) ∈ bonds
        ra = @SVector [rs[1, a], rs[2, a], rs[3, a]]
        rb = @SVector [rs[1, b], rs[2, b], rs[3, b]]
        rc = @SVector [rs[1, c], rs[2, c], rs[3, c]]

        rba = ra - rb
        rbc = rc - rb

        currenct_angle = acos(dot(rba, rbc) / (norm(rba) * norm(rbc)))
        e_valence += k * (currenct_angle - valence_angle)^2
    end
    return e_valence / 2
end

function potential_energy(sr::SimulationResult, time::Real)
    e_potential = 0
    coordinates = get_position(sr, time)
    return potential_energy(coordinates, sr.simulation)
end

function total_energy(sr::SimulationResult, time::Real)
    e_kin = kinetic_energy(sr, time)
    e_pot = potential_energy(sr, time)
    e_kin + e_pot
end

function initial_energy(simulation::NBodySimulation)
    (u0, v0, n) = gather_bodies_initial_coordinates(simulation)
    ms = get_masses(simulation.system)
    return potential_energy(u0, simulation) + kinetic_energy(v0, ms)
end

function run_simulation(s::NBodySimulation, args...; kwargs...)
    if s.thermostat isa LangevinThermostat
        calculate_simulation_sde(s, args...; kwargs...)
    else
        calculate_simulation(s, args...; kwargs...)
    end
end

function calculate_simulation(s::NBodySimulation, alg_type=Tsit5(), args...; kwargs...)
    solution = solve(ODEProblem(s), alg_type, args...; kwargs...)
    return SimulationResult(solution, s)
end

# this should be a method for integrators designed for the SecondOrderODEProblem (It is worth somehow to sort them from other algorithms)
function calculate_simulation(s::NBodySimulation, alg_type::Union{VelocityVerlet,DPRKN6,Yoshida6}, args...; kwargs...)
    cb = obtain_callbacks_for_so_ode_problem(s)
    solution = solve(SecondOrderODEProblem(s), alg_type, args...; callback=cb, kwargs...)
    return SimulationResult(solution, s)
end

function calculate_simulation_sde(s::NBodySimulation, args...; kwargs...)
    solution = solve(SDEProblem(s), args...; kwargs...)
    return SimulationResult(solution, s)
end

function obtain_callbacks_for_so_ode_problem(s::NBodySimulation)
    callback_array = Vector{DiffEqBase.DECallback}()

    if s.thermostat isa AndersenThermostat
        push!(callback_array, get_andersen_thermostating_callback(s))
    end

    return CallbackSet(tuple(callback_array...)...)
end

function get_andersen_thermostating_callback(s::NBodySimulation)
    p = s.thermostat::AndersenThermostat
    n = length(s.system.bodies)

    condition = function (u, t, integrator)
        true
    end
    affect! = function (integrator)
        collision_prob = p.ν * (integrator.t - integrator.tprev)
        for i = 1:n
            if rand() < collision_prob
                apply_andersen_rescaling_velocity(integrator, i, s.kb, s.system, p)
            end
        end
    end
    cb = DiscreteCallback(condition, affect!)
end

function apply_andersen_rescaling_velocity(integrator, i, kb, system::PotentialNBodySystem, p::AndersenThermostat)
    v_dev = sqrt(kb * p.T / system.bodies[1].m)
    @. integrator.u.x[1][:,i] = v_dev * randn()
end

function apply_andersen_rescaling_velocity(integrator, i, kb, system::WaterSPCFw, p::AndersenThermostat)
    vO = sqrt(kb*p.T/system.mO)
    vH = sqrt(kb*p.T/system.mH)
    @. integrator.u.x[1][:,3*(i-1)+1] = vO * randn()
    @. integrator.u.x[1][:,3*(i-1)+2] = vH * randn()
    @. integrator.u.x[1][:,3*(i-1)+3] = vH * randn()
end

@recipe function generate_data_for_scatter(sr::SimulationResult{<:PotentialNBodySystem}, time::Real=0.0)
    solution = sr.solution
    n = length(sr.simulation.system.bodies)

    if :gravitational ∈ keys(sr.simulation.system.potentials)

        xlims --> 1.1 * [minimum(solution[1,1:n,:]), maximum(solution[1,1:n,:])]
        ylims --> 1.1 * [minimum(solution[2,1:n,:]), maximum(solution[2,1:n,:])]

        for i in 1:n
            @series begin
                label --> "Orbit $i"
                vars --> (3 * (i - 1) + 1, 3 * (i - 1) + 2)
                solution
            end
        end
    else
        positions = get_position(sr, time)
        borders = sr.simulation.boundary_conditions
        if borders isa PeriodicBoundaryConditions
            xlims --> 1.1 * [borders[1], borders[2]]
            ylims --> 1.1 * [borders[3], borders[4]]
            zlims --> 1.1 * [borders[5], borders[6]]
        elseif borders isa CubicPeriodicBoundaryConditions
            xlims --> 1.1 * [0, borders.L]
            ylims --> 1.1 * [0, borders.L]
            zlims --> 1.1 * [0, borders.L]
            map!(x ->  x -= borders.L * floor(x / borders.L), positions, positions)
        end

    end

    positions = get_position(sr, time)
    seriestype --> :scatter
    markersize --> 5

    for i in 1:n
        @series begin
            label --> "Body no. $i"

            if all(positions[3,:] .- positions[3,1] .< 1e-15)
                ([positions[1,i]], [positions[2,i]])
            else
                (positions[1,:], positions[2,:], positions[3,:])
            end
        end
    end
end

function distancies(result::SimulationResult, time::Real)
    n = length(result.simulation.system.bodies)
    cc = get_position(result, time)

    d = Float64[]
    for i = 1:n
        for j = 1:n
            if i != j
                push!(d, norm(cc[:,i] - cc[:,j]))
            end
        end
    end
    return d
end

@recipe function initial_distribution(sr::SimulationResult{<:WaterSPCFw}, time::Real=0.0)

    n = length(sr.simulation.system.bodies)

    borders = sr.simulation.boundary_conditions

    cc = get_position(sr, time)

    if borders isa PeriodicBoundaryConditions
        xlim --> 1.1 * [borders[1], borders[2]]
        ylim --> 1.1 * [borders[3], borders[4]]
        zlim --> 1.1 * [borders[5], borders[6]]
    elseif borders isa CubicPeriodicBoundaryConditions
        xlim --> 1.1 * [0, borders.L]
        ylim --> 1.1 * [0, borders.L]
        zlim --> 1.1 * [0, borders.L]

        map!(x ->  x -= borders.L * floor(x / borders.L), cc, cc)
    end
    seriestype --> :scatter

    @series begin
        label --> "O"
        markersize --> 8
        markercolor --> :red
        (cc[1,1:3:3 * n - 2], cc[2,1:3:3 * n - 2], cc[3,1:3:3 * n - 2])
    end

    @series begin
        label --> "H"
        markersize --> 4
        markercolor --> :green
        x = vcat(cc[1,2:3:3 * n - 1], cc[1,3:3:3 * n])
        y = vcat(cc[2,2:3:3 * n - 1], cc[2,3:3:3 * n])
        z = vcat(cc[3,2:3:3 * n - 1], cc[3,3:3:3 * n])
        (x, y, z)
    end
end

function rdf(sr::SimulationResult)
    n = length(sr.simulation.system.bodies)
    pbc = sr.simulation.boundary_conditions

    (ms, indxs) = obtain_data_for_lennard_jones_interaction(sr.simulation.system)
    indlen = length(indxs)

    maxbin = 1000
    dr = pbc.L / maxbin
    hist = zeros(maxbin)
    for t ∈ sr.solution.t
        cc = get_position(sr, t)
        for ind_i = 1:indlen
            i = indxs[ind_i]
            ri = @SVector [cc[1, i], cc[2, i], cc[3, i]]
            for ind_j = ind_i + 1:indlen
                j = indxs[ind_j]
                rj = @SVector [cc[1, j], cc[2, j], cc[3, j]]

                (rij, r, r2) = get_interparticle_distance(ri, rj, pbc)

                if r2 < (0.5 * pbc.L)^2
                    bin = ceil(Int, r / dr)
                    if bin > 1 && bin <= maxbin
                        hist[bin] += 2
                    end
                end
            end
        end
    end

    c = 4 / 3 * π * indlen / pbc.L^3

    gr = zeros(maxbin)
    rs = zeros(maxbin)
    tlen = length(sr.solution.t)
    for bin = 1:maxbin
        rlower = (bin - 1) * dr
        rupper = rlower + dr
        nideal = c * (rupper^3 - rlower^3)
        gr[bin] = (hist[bin] / (tlen * indlen)) / nideal
        rs[bin] = rlower + dr / 2
    end

    return (rs, gr)
end

function msd(sr::SimulationResult{<:PotentialNBodySystem})
    n = length(sr.simulation.system.bodies)

    (ms, indxs) = obtain_data_for_lennard_jones_interaction(sr.simulation.system)
    indlen = length(indxs)

    ts = sr.solution.t
    tlen = length(ts)
    dr2 = zeros(length(ts))

    cc0 = get_position(sr, ts[1])

    for t = 1:tlen
        cc = get_position(sr, ts[t])
        for ind_i = 1:indlen
            i = indxs[ind_i]
            dr = @SVector [cc[1, i] - cc0[1,i], cc[2, i] - cc0[2,i], cc[3, i] - cc0[3,i]]

            dr2[t] += dot(dr, dr)
        end
        dr2[t] /= indlen
    end

    (ts, dr2)
end

function msd(sr::SimulationResult{<:WaterSPCFw})
    n = length(sr.simulation.system.bodies)

    mO = sr.simulation.system.mO
    mH = sr.simulation.system.mH

    ts = sr.solution.t
    tlen = length(ts)
    dr2 = zeros(tlen)

    cc0 = get_position(sr, ts[1])

    for t = 1:tlen
        cc = get_position(sr, ts[t])
        for i = 1:n
            indO, indH1, indH2 = 3 * (i - 1) + 1, 3 * (i - 1) + 2, 3 * (i - 1) + 3
            dr = ((cc[:,indO] - cc0[:,indO]) * mO + (cc[:,indH1] - cc0[:,indH1]) * mH + (cc[:,indH2] - cc0[:,indH2]) * mH) / (2 * mH + mO)

            dr2[t] += dot(dr, dr)
        end
        dr2[t] /= n
    end

    (ts, dr2)
end

# coordinates should be in angstroms
function save(f::File{format"ProteinDataBank"}, sr::SimulationResult)
    open(f, "w") do s
        write_pdb_data(s,sr)
    end
end

function write_pdb_data(f::IO, sr::SimulationResult{<:WaterSPCFw})
    L = 10*sr.simulation.boundary_conditions.L
    strL = @sprintf("%9.3f",L)
    strA = @sprintf("%7.2f",90)
    #println(f,"CRYST1",lpad(strL,9),lpad(strL,9),lpad(strL,9),lpad(strA,7),lpad(strA,7),lpad(strA,7))
    n = length(sr.simulation.system.bodies)
    cc = zeros(3,3*n)
    count = 0
    for t in sr.solution.t
        cc0 = 10*get_position(sr, t)
        map!(x ->  x -= L * floor(x / L), cc, cc0)

        for i ∈ 1:n
            indO, indH1, indH2 = 3 * (i - 1) + 1, 3 * (i - 1) + 2, 3 * (i - 1) + 3
            @. cc[:, indH1] = cc[:,indO] + cc0[:, indH1] - cc0[:, indO]
            @. cc[:, indH2] = cc[:,indO] + cc0[:, indH2] - cc0[:, indO]
        end

        count+=1
        println(f,rpad("MODEL",10), count)
        println(f,"REMARK 250 time=$t picoseconds")
        for i ∈ 1:n
            indO, indH1, indH2 = 3 * (i - 1) + 1, 3 * (i - 1) + 2, 3 * (i - 1) + 3

            println(f,"HETATM",lpad(indO,5),"  ",rpad("O",4),"HOH",lpad(i,6),"    ", lpad(@sprintf("%8.3f",cc[1,indO]),8), lpad(@sprintf("%8.3f",cc[2,indO]),8), lpad(@sprintf("%8.3f",cc[3,indO]),8),lpad("1.00",6),lpad("0.00",6), lpad("",10), lpad("O",2))
            println(f,"HETATM",lpad(indH1,5),"  ",rpad("H1",4),"HOH",lpad(i,6),"    ", lpad(@sprintf("%8.3f",cc[1,indH1]),8), lpad(@sprintf("%8.3f",cc[2,indH1]),8), lpad(@sprintf("%8.3f",cc[3,indH1]),8),lpad("1.00",6),lpad("0.00",6), lpad("",10), lpad("H",2))
            println(f,"HETATM",lpad(indH2,5),"  ",rpad("H2",4),"HOH",lpad(i,6),"    ", lpad(@sprintf("%8.3f",cc[1,indH2]),8), lpad(@sprintf("%8.3f",cc[2,indH2]),8), lpad(@sprintf("%8.3f",cc[3,indH2]),8),lpad("1.00",6),lpad("0.00",6), lpad("",10), lpad("H",2))
            end
        println(f,"ENDMDL")
    end
end

function write_pdb_data(f::IO, sr::SimulationResult)
    n = length(sr.simulation.system.bodies)
    L = 10*sr.simulation.boundary_conditions.L
    count = 0
    for t in sr.solution.t
        cc = 10*get_position(sr, t)
        map!(x ->  x -= L * floor(x / L), cc, cc)
        count+=1
        println(f,rpad("MODEL",10), count)
        println(f,"REMARK 250 time=$t steps")
        for i ∈ 1:n

            println(f,"HETATM",lpad(i,5),"  ",rpad("Ar",4),"Ar",lpad(i,6),"    ", lpad(@sprintf("%8.3f",cc[1,i]),8), lpad(@sprintf("%8.3f",cc[2,i]),8), lpad(@sprintf("%8.3f",cc[3,i]),8),lpad("1.00",6),lpad("0.00",6), lpad("",10), lpad("Ar",2))
        end
        println(f,"ENDMDL")
    end
end

function save_to_pdb(sr::SimulationResult, f)
    open(f, "w") do s
        write_pdb_data(s,sr)
    end
end

function load(f::File{format"ProteinDataBank"})
    open(f) do s
        skipmagic(s)  # skip over the magic bytes
        extract_from_pdb(s)
    end
end

function load_water_molecules_from_pdb(path_str)
    open(path_str) do file
        extract_from_pdb(file)
    end
end

function extract_from_pdb(file)
    cc = SVector[]
    wms = WaterMolecule[]
    t_count = 0
    while !eof(file)
        ln = readline(file)
        if ln[1:5] == "MODEL"
            break
        end
    end

    ln_time = readline(file)
    parts = split(ln_time)
    t1 = parse(Float64, parts[3][6:end])
    den = 10
    while !eof(file)
        ln_cc = readline(file)
        if ln_cc[1:6]=="ENDMDL"
            break
        elseif ln_cc[14]=='O'
            ps_o = split(ln_cc)
            cc_o = SVector(parse(Float64, ps_o[6]), parse(Float64, ps_o[7]), parse(Float64, ps_o[8]))/den
            ps_h1 = split(readline(file))
            cc_h1 = SVector(parse(Float64, ps_h1[6]), parse(Float64, ps_h1[7]), parse(Float64, ps_h1[8]))/den
            ps_h2 = split(readline(file))
            cc_h2 = SVector(parse(Float64, ps_h2[6]), parse(Float64, ps_h2[7]), parse(Float64, ps_h2[8]))/den
            push!(cc, cc_o)
            push!(cc, cc_h1)
            push!(cc, cc_h2)
        end
    end

    while !eof(file)
        ln = readline(file)
        if ln[1:5] == "MODEL"
            break
        end
    end

    ln_time = readline(file)
    parts = split(ln_time)
    t2 = parse(Float64, parts[3][6:end])
    τ = t2-t1

    mc = 0
    while !eof(file)
        ln_cc = readline(file)
        if ln_cc[1:6]=="ENDMDL"
            break
        elseif  ln_cc[14]=='O'
            mc+=1
            ps_o = split(ln_cc)
            cc_o = SVector(parse(Float64, ps_o[6]), parse(Float64, ps_o[7]), parse(Float64, ps_o[8]))/den
            ps_h1 = split(readline(file))
            cc_h1 = SVector(parse(Float64, ps_h1[6]), parse(Float64, ps_h1[7]), parse(Float64, ps_h1[8]))/den
            ps_h2 = split(readline(file))
            cc_h2 = SVector(parse(Float64, ps_h2[6]), parse(Float64, ps_h2[7]), parse(Float64, ps_h2[8]))/den

            o = MassBody(cc[3*(mc-1)+1], (cc_o-cc[3*(mc-1)+1])/τ, 15.999)
            h1 = MassBody(cc[3*(mc-1)+2], (cc_h1-cc[3*(mc-1)+2])/τ, 1.00794)
            h2 = MassBody(cc[3*(mc-1)+3], (cc_h2-cc[3*(mc-1)+3])/τ, 1.00794)

            o = MassBody(cc[3*(mc-1)+1], SVector(0.0,0.0,0.0), 15.999)
            h1 = MassBody(cc[3*(mc-1)+2], SVector(0.0,0.0,0.0), 1.00794)
            h2 = MassBody(cc[3*(mc-1)+3], SVector(0.0,0.0,0.0), 1.00794)

            push!(wms, WaterMolecule(o,h1,h2))
        end
    end

    wms
end


text:
using Codecs

function main(tries)
  str_size = 10_000_000
  str = repeat("a", str_size)
  str2 = ""

  print("encode: ")
  t = time()
  s = 0
  for i in range(0, tries)
    str2 = ASCIIString(encode(Base64, str))
    s += length(str2)
  end
  print(s, ", ", time() - t, "\n")

  print("decode: ")
  t = time()
  s = 0
  for i in range(0, tries)
    s += length(ASCIIString(decode(Base64, str2)))
  end
  print(s, ", ", time() - t, "\n")
end

println(STDERR, "warming")
main(5)

println(STDERR, "bench")
x = @timed main(100)
println(STDERR, "Elapsed: $(x[2]), Allocated: $(x[3]), GC Time: $(x[4])")

ext: jl

max_stars_repo_name: trishume/benchmark-stadium

max_issues_repo_name: trishume/benchmark-stadium

max_forks_repo_name: trishume/benchmark-stadium


text:
module GroupedErrors

using Statistics
import Loess, KernelDensity, IterableTables
import DataValues: DataValue, isna
import Lazy: @>
using IndexedTables
import IndexedTables: select
using TableTraits
using MacroTools
using StatsBase
using ShiftedArrays

export @splitby, @bootstrap, @across, @x, @y, @xy, @compare, @summarize, @set_attr, @>, @plot
export @xlims, @ylims
export ProcessedTable

include("select.jl")
include("query_macro.jl")
include("plot_macro.jl")
include("pipeline.jl")
include("analysisfunctions.jl")

exampletable() = joinpath(dirname(@__DIR__), "test", "tables")
exampletable(s) = joinpath(exampletable(), s)

end


text:
using Compat.Test
using TestSetExtensions

@testset ExtendedTestSet "Nice Diffs" begin
    @test [3, 5, 6, 1, 6, 8] == [3, 5, 6, 1, 9, 8]
    @test """Lorem ipsum dolor sit amet,
             consectetur adipiscing elit, sed do
             eiusmod tempor incididunt ut
             labore et dolore magna aliqua.
             Ut enim ad minim veniam, quis nostrud
             exercitation ullamco aboris.""" ==
          """Lorem ipsum dolor sit amet,
             consectetur adipiscing elit, sed do
             eiusmod temper incididunt ut
             labore et dolore magna aliqua.
             Ut enim ad minim veniam, quis nostrud
             exercitation ullamco aboris."""
    @test Dict(:foo => "bar", :baz => [1, 4, 5], :biz => nothing) ==
          Dict(:baz => [1, 7, 5], :biz => 42)
end


text:
using MLJ
using  PrettyPrinting
import DataFrames: DataFrame
import Statistics

X, y = @load_reduced_ames
X = DataFrame(X)
@show size(X)
first(X, 3) |> pretty

@show y[1:3]
scitype(y)

creg = ConstantRegressor()

cmach = machine(creg, X, y)

train, test = partition(collect(eachindex(y)), 0.70, shuffle=true); # 70:30 split
fit!(cmach, rows=train)
ŷ = predict(cmach, rows=test)
ŷ[1:3] |> pprint

ŷ = predict_mean(cmach, rows=test)
ŷ[1:3]

rmsl(ŷ, y[test])

RidgeRegressor = @load RidgeRegressor pkg="MultivariateStats"
KNNRegressor = @load KNNRegressor

Xs = source(X)
ys = source(y)

hot = machine(OneHotEncoder(), Xs)

W = transform(hot, Xs)
z = log(ys);

knn   = machine(KNNRegressor(K=5), W, z)
ridge = machine(RidgeRegressor(lambda=2.5), W, z)

ẑ₁ = predict(ridge, W)
ẑ₂ = predict(knn, W)

ẑ = 0.3ẑ₁ + 0.7ẑ₂;

ŷ = exp(ẑ);

fit!(ŷ, rows=train)
ypreds = ŷ(rows=test)
rmsl(y[test], ypreds)

mutable struct KNNRidgeBlend <: DeterministicComposite
    knn_model::KNNRegressor
    ridge_model::RidgeRegressor
    knn_weight::Float64
end

function MLJ.fit(model::KNNRidgeBlend, verbosity::Int, X, y)
    Xs = source(X)
    ys = source(y)
    hot = machine(OneHotEncoder(), Xs)
    W = transform(hot, Xs)
    z = log(ys)
    ridge_model = model.ridge_model
    knn_model = model.knn_model
    ridge = machine(ridge_model, W, z)
    knn = machine(knn_model, W, z)

    ẑ = model.knn_weight * predict(knn, W) + (1.0 - model.knn_weight) * predict(ridge, W)
    ŷ = exp(ẑ)

    mach = machine(Deterministic(), Xs, ys; predict=ŷ)
    return!(mach, model, verbosity)
end

krb = KNNRidgeBlend(KNNRegressor(K=5), RidgeRegressor(lambda=2.5), 0.3)
mach = machine(krb, X, y)
fit!(mach, rows=train)

preds = predict(mach, rows=test)
rmsl(y[test], preds)

@show krb.knn_weight
@show krb.knn_model.K
@show krb.ridge_model.lambda

params(krb) |> pprint

k_range = range(krb, :(knn_model.K), lower=2, upper=100, scale=:log10)
l_range = range(krb, :(ridge_model.lambda), lower=1e-4, upper=10, scale=:log10)
w_range = range(krb, :(knn_weight), lower=0.1, upper=0.9)

ranges = [k_range, l_range, w_range]

tuning = Grid(resolution=3)
resampling = CV(nfolds=6)

tm = TunedModel(model=krb, tuning=tuning, resampling=resampling,
                ranges=ranges, measure=rmsl)

mtm = machine(tm, X, y)
fit!(mtm, rows=train);

krb_best = fitted_params(mtm).best_model
@show krb_best.knn_model.K
@show krb_best.ridge_model.lambda
@show krb_best.knn_weight

preds = predict(mtm, rows=test)
rmsl(y[test], preds)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@testset "Sampling Density" begin

# create a 10x10 grid of unit spaced sampling points
N = 10
g = (0:(N-1)) ./ N .- 0.5
x = vec(ones(N) * g')
y = vec(g * ones(N)')
nodes = cat(x',y', dims=1)

# approximate the density weights
p = plan_nfft(nodes, (N,N), m = 5, σ = 2.0);
weights = sdc(p, iters = 10)

# test if they approximate the true weights (1/(N*N))
@test all( (≈).(vec(weights), 1/(N*N), rtol=1e-7) )

end


text:
module AppDirs

export user_data_dir, site_data_dir,
       user_config_dir, site_config_dir,
       user_cache_dir, user_state_dir, user_log_dir

@doc raw"""
    user_data_dir(appname=nothing, appauthor=nothing, version=nothing,
                  roaming::Bool=false)

Return full path to the user-specific data dir for this application.

- `appname` is the name of application. If `nothing`, just the system directory
  is returned.

- `appauthor` (only used on Windows) is the name of the
  appauthor or distributing body for this application. Typically
  it is the owning company name. This falls back to `appname`. You may
  pass `false` to disable it.

- `version` is an optional version path element to append to the path. You
  might want to use this if you want multiple versions of your app to be able
  to run independently. If used, this would typically be `"<major>.<minor>"`.
  Only applied when `appname` is present.

- `roaming` (default `false`) can be set `true` to use the Windows roaming
  appdata directory. That means that for users on a Windows network setup for
  roaming profiles, this user data will be sync'd on login. See
  [here](http://technet.microsoft.com/en-us/library/cc766489(WS.10).aspx)
  for a discussion of issues.

Typical user data directories are:

| Setting              | Path                                                                                |
|:---------------------|:---------------------------------------------------------------------------------------------|
| macOS                | `~/Library/Application Support/<AppName>`                                                    |
| Unix                 | `~/.local/share/<AppName>` (or in `XDG_DATA_HOME`, if defined)                               |
| Win XP (not roaming) | `C:\Documents and Settings\<username>\Application Data\<AppAuthor>\<AppName>`                |
| Win XP (roaming)     | `C:\Documents and Settings\<username>\Local Settings\Application Data\<AppAuthor>\<AppName>` |
| Win 7  (not roaming) | `C:\Users\<username>\AppData\Local\<AppAuthor>\<AppName>`                                    |
| Win 7  (roaming)     | `C:\Users\<username>\AppData\Roaming\<AppAuthor>\<AppName>`                                  |
"""
function user_data_dir(appname=nothing, appauthor=nothing; version=nothing,
                       roaming::Bool=false)
  if Sys.iswindows()
    if isnothing(appauthor)
      appauthor = appname
    end
    var = roaming ? "CSIDL_APPDATA" : "CSIDL_LOCAL_APPDATA"
    path = normpath(_get_win_folder(var))
    if !isnothing(appname)
      if appauthor != false
        path = joinpath(path, appauthor, appname)
      else
        path = joinpath(path, appname)
      end
    end
  elseif Sys.isapple()
    path = expanduser("~/Library/Application Support/")
    if !isnothing(appname)
      path = joinpath(path, appname)
    end
  else
    path = get(ENV, "XDG_DATA_HOME", expanduser("~/.local/share"))
    if !isnothing(appname)
      path = joinpath(path, appname)
    end
  end
  if !isnothing(appname) && !isnothing(version)
    path = joinpath(path, version)
  end
  return path
end

@doc raw"""
    site_data_dir(appname=nothing, appauthor=nothing; version=nothing,
                  multipath::Bool=false)

Return full path to the user-shared data dir for this application.

- `appname` is the name of application. If `nothing`, just the system directory
  is returned.

- `appauthor` (only used on Windows) is the name of the appauthor or
  distributing body for this application. Typically it is the owning company
  name. This falls back to `appname`. You may pass `false` to disable it.

- `version` is an optional version path element to append to the path. You might
  want to use this if you want multiple versions of your app to be able to run
  independently. If used, this would typically be `"<major>.<minor>"`. Only
  applied when `appname` is present.

- `multipath` is an optional parameter only applicable to *nix which indicates
  that the entire list of data dirs should be returned. By default, the first
  item from `XDG_DATA_DIRS` is returned, or `/usr/local/share/<AppName>`, if
  `XDG_DATA_DIRS` is not set

Typical site data directories are:

| Setting   | Path                                                                         |
|:----------|:-----------------------------------------------------------------------------|
| macOS     | `/Library/Application Support/<AppName>`                                     |
| Unix      | `/usr/local/share/<AppName>` or `/usr/share/<AppName>`                       |
| Win XP    | `C:\Documents and Settings\All Users\Application Data\<AppAuthor>\<AppName>` |
| Win Vista | (Fail! `C:\ProgramData` is a hidden *system* directory on Vista)             |
| Win 7     | `C:\ProgramData\<AppAuthor>\<AppName>` (Hidden, but writeable on Win 7)      |

!!! warning

    Do not use this on Windows. See the Vista-Fail note above for why.
"""
function site_data_dir(appname=nothing, appauthor=nothing; version=nothing,
                       multipath::Bool=false)
  if Sys.iswindows()
    if isnothing(appauthor)
      appauthor = appname
    end
    path = normpath(_get_win_folder("CSIDL_COMMON_APPDATA"))
    if !isnothing(appname)
      if appauthor != false
        path = joinpath(path, appauthor, appname)
      else
        path = joinpath(path, appname)
      end
    end
  elseif Sys.isapple()
    path = expanduser("/Library/Application Support")
    if !isnothing(appname)
      path = joinpath(path, appname)
    end
  else
    # XDG default for $XDG_DATA_DIRS
    # only first, if multipath is False
    pathsep = Base.Filesystem.path_separator
    pathlistsep = Sys.iswindows() ? ';' : ':'
    path = get(ENV, "XDG_DATA_DIRS",
                    join(["/usr/local/share", "/usr/share"], pathlistsep))
    pathlist = [expanduser(rstrip(x, pathsep[1]))
                for x in split(path, pathlistsep)]
    if !isnothing(appname)
      if !isnothing(version)
        appname = joinpath(appname, version)
      end
      pathlist = [join([x, appname], pathsep) for x in pathlist]
    end

    if multipath
      path = join(pathlist, pathlistsep)
    else
      path = pathlist[1]
    end
    return path
  end

  if !isnothing(appname) && !isnothing(version)
    path = joinpath(path, version)
  end
  return path
end

@doc raw"""
    user_config_dir(appname=nothing, appauthor=nothing; version=nothing,
                    roaming::Bool=false)

Return full path to the user-specific config dir for this application.

- `appname` is the name of application. If `nothing`, just the system directory
  is returned.

- `appauthor` (only used on Windows) is the name of the appauthor or
  distributing body for this application. Typically it is the owning company
  name. This falls back to `appname`. You may pass `false` to disable it.

- `version` is an optional version path element to append to the path. You might
  want to use this if you want multiple versions of your app to be able to run
  independently. If used, this would typically be `"<major>.<minor>"`. Only
  applied when `appname` is present.

- `roaming` (default `false`) can be set `true` to use the Windows roaming
  appdata directory. That means that for users on a Windows network setup for
  roaming profiles, this user data will be sync'd on login. See
  [here](http://technet.microsoft.com/en-us/library/cc766489(WS.10).aspx)
  for a discussion of issues.

Typical user config directories are:

| Setting  | Path                                                        |
|:---------|:------------------------------------------------------------|
| macOS    | `~/Library/Preferences/<AppName>`                           |
| Unix     | `~/.config/<AppName>` (or in `XDG_CONFIG_HOME`, if defined) |
| Windows  | same as [`user_data_dir`](@ref)                             |
"""
function user_config_dir(appname=nothing, appauthor=nothing; version=nothing,
                         roaming::Bool=false)
  if Sys.iswindows()
    path = user_data_dir(appname, appauthor; version=nothing, roaming)
  elseif Sys.isapple()
    path = expanduser("~/Library/Preferences/")
    if !isnothing(appname)
      path = joinpath(path, appname)
    end
  else
    path = get(ENV, "XDG_CONFIG_HOME", expanduser("~/.config"))
    if !isnothing(appname)
      path = joinpath(path, appname)
    end
  end
  if !isnothing(appname) && !isnothing(version)
    path = joinpath(path, version)
  end
  return path
end

@doc raw"""
    site_config_dir(appname=nothing, appauthor=nothing; version=nothing,
                    multipath::Bool=false)

Return full path to the user-shared config dir for this application.

- `appname` is the name of application. If `nothing`, just the system directory
  is returned.

- `appauthor` (only used on Windows) is the name of the appauthor or
  distributing body for this application. Typically it is the owning company
  name. This falls back to `appname`. You may pass `false` to disable it.

- `version` is an optional version path element to append to the path. You might
  want to use this if you want multiple versions of your app to be able to run
  independently. If used, this would typically be `"<major>.<minor>"`. Only
  applied when `appname` is present.

- `multipath` is an optional parameter only applicable to *nix which indicates
  that the entire list of config dirs should be returned. By default, the first
  item from `XDG_CONFIG_DIRS` is returned, or `/etc/xdg/<AppName>`, if
  `XDG_CONFIG_DIRS` is not set

Typical site config directories are:

| Setting   | Path                                                                                |
|:----------|:-------------------------------------------------------------------------------------------|
| macOS     | `/Library/Preferences`                                                            |
| Unix      | `/etc/xdg/<AppName>` or `XDG_CONFIG_DIRS[i]/<AppName>` for each value in `XDG_CONFIG_DIRS` |
| Windows   | same as [`site_data_dir`](@ref)                                                            |
| Win Vista | (Fail! `C:\ProgramData` is a hidden *system* directory on Vista)                           |

!!! warning

    Do not use this on Windows. See the Vista-Fail note above for why.
"""
function site_config_dir(appname=nothing, appauthor=nothing; version=nothing,
                         multipath::Bool=false)
  if Sys.iswindows()
    path = site_data_dir(appname, appauthor)
    if !isnothing(appname) && !isnothing(version)
      path = joinpath(path, version)
    end
  elseif Sys.isapple()
    path = expanduser("/Library/Preferences")
    if !isnothing(appname)
      path = joinpath(path, appname)
    end
    if !isnothing(appname) && !isnothing(version)
      path = joinpath(path, version)
    end
  else
    # XDG default for $XDG_CONFIG_DIRS
    # only first, if multipath is False
    pathsep = Base.Filesystem.path_separator
    pathlistsep = Sys.iswindows() ? ';' : ':'
    path = get(ENV, "XDG_CONFIG_DIRS", "/etc/xdg")
    pathlist = [expanduser(rstrip(x, pathsep[1]))
                for x in split(path, pathlistsep)]
    if !isnothing(appname)
      if !isnothing(version)
        appname = joinpath(appname, version)
      end
      pathlist = [join([x, appname], pathsep) for x in pathlist]
    end
    if multipath
      path = join(pathlist, pathlistsep)
    else
      path = pathlist[1]
    end
  end
  return path
end

@doc raw"""
    user_cache_dir(appname=nothing, appauthor=nothing; version=nothing,
                   opinion::Bool=true)

Return full path to the user-specific cache dir for this application.

- `appname` is the name of application. If `nothing`, just the system directory
  is returned.

- `appauthor` (only used on Windows) is the name of the appauthor or
  distributing body for this application. Typically it is the owning company
  name. This falls back to `appname`. You may pass `false` to disable it.

- `version` is an optional version path element to append to the path. You might
  want to use this if you want multiple versions of your app to be able to run
  independently. If used, this would typically be `"<major>.<minor>"`. Only
  applied when `appname` is present.

- `opinion` (default `true`) can be `false` to disable the appending of `Cache`
  to the base app data dir for Windows. See discussion below.

Typical user cache directories are:

| Setting   | Path                                                                                |
|:----------|:---------------------------------------------------------------------------------------------------|
| macOS     | `~/Library/Caches/<AppName>`                                                                       |
| Unix      | `~/.cache/<AppName>`  (or in `XDG_CACHE_HOME`, if defined)                                         |
| Win XP    | `C:\Documents and Settings\<username>\Local Settings\Application Data\<AppAuthor>\<AppName>\Cache` |
| Win 7     | `C:\Users\<username>\AppData\Local\<AppAuthor>\<AppName>\Cache`                                    |

On Windows the only suggestion in the MSDN docs is that local settings go in the
`CSIDL_LOCAL_APPDATA` directory. This is identical to the non-roaming app data
dir (the default returned by [`user_data_dir`](@ref)). Apps typically put cache
data somewhere *under* the given dir here. Some examples:

- `...\Mozilla\Firefox\Profiles\<ProfileName>\Cache`
- `...\Acme\SuperApp\Cache\1.0`
"""
function user_cache_dir(appname=nothing, appauthor=nothing; version=nothing,
                        opinion::Bool=true)
  if Sys.iswindows()
    if isnothing(appauthor)
      appauthor = appname
    end
    path = normpath(_get_win_folder("CSIDL_LOCAL_APPDATA"))
    if !isnothing(appname)
      if appauthor != false
        path = joinpath(path, appauthor, appname)
      else
        path = joinpath(path, appname)
      end
      if !isnothing(version)
        path = joinpath(path, version)
      end
      if opinion
        path = joinpath(path, "Cache")
      end
    end
    return path
  elseif Sys.isapple()
    path = expanduser("~/Library/Caches")
    if !isnothing(appname)
      path = joinpath(path, appname)
    end
  else
    path = get(ENV, "XDG_CACHE_HOME", expanduser("~/.cache"))
    if !isnothing(appname)
      path = joinpath(path, appname)
    end
  end
  if !isnothing(appname) && !isnothing(version)
    path = joinpath(path, version)
  end
  return path
end

@doc raw"""
    user_state_dir(appname=nothing, appauthor=nothing; version=nothing,
                   roaming::Bool=false)

Return full path to the user-specific state dir for this application.

- `appname` is the name of application. If `nothing`, just the system directory
  is returned.

- `appauthor` (only used on Windows) is the name of the appauthor or
  distributing body for this application. Typically it is the owning company
  name. This falls back to `appname`. You may pass `false` to disable it.

- `version` is an optional version path element to append to the path. You might
  want to use this if you want multiple versions of your app to be able to run
  independently. If used, this would typically be `"<major>.<minor>"`. Only
  applied when `appname` is present.

- `roaming` (default `false`) can be set `true` to use the Windows roaming
  appdata directory. That means that for users on a Windows network setup for
  roaming profiles, this user data will be sync'd on login. See
  [here](http://technet.microsoft.com/en-us/library/cc766489(WS.10).aspx)
  for a discussion of issues.

Typical user state directories are:

| Setting   | Path                                                        |
|:----------|:------------------------------------------------------------|
| macOS     | same as [`user_data_dir`](@ref)                             |
| Unix      | `~/.local/state/<AppName>` (or `XDG_STATE_HOME` if defined) |
| Windows   | same as [`user_data_dir`](@ref)                             |
"""
function user_state_dir(appname=nothing, appauthor=nothing; version=nothing,
                        roaming::Bool=false)
  if Sys.iswindows() || Sys.isapple()
    path = user_data_dir(appname, appauthor; version=nothing, roaming)
  else
    path = get(ENV, "XDG_STATE_HOME", expanduser("~/.local/state"))
    if !isnothing(appname)
      path = joinpath(path, appname)
    end
  end
  if !isnothing(appname) && !isnothing(version)
    path = joinpath(path, version)
  end
  return path
end

@doc raw"""
    user_log_dir(appname=nothing, appauthor=nothing, version=nothing,
                 opinion::Bool=true)

Return full path to the user-specific log dir for this application.

- `appname` is the name of application. If `nothing`, just the system directory
  is returned.

- `appauthor` (only used on Windows) is the name of the appauthor or
  distributing body for this application. Typically it is the owning company
  name. This falls back to `appname`. You may pass `false` to disable it.

- `version` is an optional version path element to append to the path. You might
  want to use this if you want multiple versions of your app to be able to run
  independently. If used, this would typically be `"<major>.<minor>"`. Only
  applied when `appname` is present.

- `opinion` (default `true`) can be `false` to disable the appending of `Logs`
  to the base app data dir for Windows, and `log` to the base cache dir for
  Unix. See discussion below.

Typical user log directories are:

| Setting   | Path                                                                                |
|:----------|:---------------------------------------------------------------------------------------------------|
| macOS     | `~/Library/Logs/<AppName>`                                                                       |
| Unix      | `~/.cache/<AppName>/log`  (or under `XDG_CACHE_HOME`, if defined)                                         |
| Win XP    | `C:\Documents and Settings\<username>\Local Settings\Application Data\<AppAuthor>\<AppName>\Logs` |
| Win 7     | `C:\Users\<username>\AppData\Local\<AppAuthor>\<AppName>\Logs`                                    |

On Windows the only suggestion in the MSDN docs is that local settings go in the
`CSIDL_LOCAL_APPDATA` directory.
"""
function user_log_dir(appname=nothing, appauthor=nothing; version=nothing,
                      opinion::Bool=true)
  if Sys.isapple()
    path = joinpath(expanduser("~/Library/Logs"), appname)
  elseif Sys.iswindows()
    path = user_data_dir(appname, appauthor; version)
    version = nothing
    if opinion
      path = joinpath(path, "Logs")
    end
  else
    path = user_cache_dir(appname, appauthor; version)
    version = nothing
    if opinion
      path = joinpath(path, "log")
    end
  end
  if !isnothing(appname) && !isnothing(version)
    path = joinpath(path, version)
  end
  return path
end

"""
    AppDir(appname=nothing, appauthor=nothing, version=nothing,
           roaming::Bool=false, multipath::Bool=false)

Convenience wrapper for getting application dirs, with the following fields:

- `data`: the output of [`user_data_dir`](@ref)
- `site_data`: the output of [`site_data_dir`](@ref)
- `config`: the output of [`user_config_dir`](@ref)
- `site_config`: the output of [`site_config_dir`](@ref)
- `cache`: the output of [`user_cache_dir`](@ref)
- `state`: the output of [`user_state_dir`](@ref)
- `log`: the output of [`user_log_dir`](@ref)
"""
struct AppDir
  data::String
  site_data::String
  config::String
  site_config::String
  cache::String
  state::String
  log::String

  function AppDir(appname=nothing, appauthor=nothing; version=nothing,
                  roaming::Bool=false, multipath::Bool=false)
    new(
        user_data_dir(appname, appauthor; version, roaming),
        site_data_dir(appname, appauthor; version, multipath),
        user_config_dir(appname, appauthor; version, roaming),
        site_config_dir(appname, appauthor; version, multipath),
        user_cache_dir(appname, appauthor; version),
        user_state_dir(appname, appauthor; version, roaming),
        user_log_dir(appname, appauthor; version),
    )
  end
end

#---- internal support stuff

function _get_win_folder(csidl_name)
  # From https://github.com/JuliaLang/IJulia.jl/blob/8710cd4bc1adabfacb8d9caec56332e36daa1f90/deps/kspec.jl#L12-L17
  csidl_const = Dict(
      "CSIDL_APPDATA" =>        UInt16(26),
      "CSIDL_COMMON_APPDATA" => UInt16(35),
      "CSIDL_LOCAL_APPDATA" =>  UInt16(28),
  )[csidl_name]

  path = zeros(UInt16, 1024)
  result = ccall((:SHGetFolderPathW,:shell32), stdcall, Cint,
          (Ptr{Cvoid},Cint,Ptr{Cvoid},Cint,Ptr{UInt16}),C_NULL,csidl_const,C_NULL,0,path)

  if result == 0
    return transcode(String, resize!(path, findfirst(iszero, path) - 1))
  end

  # Fallback to ENV
  env_var_name = Dict(
      "CSIDL_APPDATA" => "APPDATA",
      "CSIDL_COMMON_APPDATA" => "ALLUSERSPROFILE",
      "CSIDL_LOCAL_APPDATA" => "LOCALAPPDATA",
  )[csidl_name]

  return get(ENV, env_var_name, "")
end

end


text:
using DiffEqProblemLibrary, Sundials, Base.Test

prob = prob_ode_linear
dt = 1//2^(4)
saveat = float(collect(0:dt:1))
sol = solve(prob,ARKODE())
@test sol.errors[:l2] < 1e-3