text:
function __solve end
function __init end
function solve! end

NO_TSPAN_PROBS = Union{AbstractLinearProblem, AbstractNonlinearProblem,
                       AbstractQuadratureProblem,
                       AbstractSteadyStateProblem,AbstractJumpProblem}

function init_call(_prob,args...;kwargs...)
  if :kwargs ∈ propertynames(_prob)
    __init(_prob,args...;_prob.kwargs...,kwargs...)
  else
    __init(_prob,args...;kwargs...)
  end
end

function init(prob::DEProblem,args...;kwargs...)
  _prob = get_concrete_problem(prob,kwargs)
  if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    init_call(_prob,alg,args...;kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    init_call(_prob,args...;kwargs...)
  else
    init_call(_prob,args...;kwargs...)
  end
end

function solve_call(_prob,args...;kwargs...)
  if :kwargs ∈ propertynames(_prob)
    __solve(_prob,args...;_prob.kwargs...,kwargs...)
  else
    __solve(_prob,args...;kwargs...)
  end
end

function solve(prob::DEProblem,args...;kwargs...)
  _prob = get_concrete_problem(prob,kwargs)
  if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    solve_call(_prob,alg,args...;kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    solve_call(_prob,args...;kwargs...)
  elseif isempty(args) # Default algorithm handling
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    solve_call(_prob,args...;kwargs...)
  else
    solve_call(_prob,args...;kwargs...)
  end
end

function solve(prob::EnsembleProblem,args...;kwargs...)
  __solve(prob,args...;kwargs...)
end

function get_concrete_problem(prob::AbstractJumpProblem,kwargs)
  prob
end

function get_concrete_problem(prob::AbstractSteadyStateProblem, kwargs)
  u0 = get_concrete_u0(prob, Inf, kwargs)
  remake(prob; u0 = u0)
end

function get_concrete_problem(prob::AbstractEnsembleProblem, kwargs)
  prob
end

function DiffEqBase.solve(prob::PDEProblem,alg::DiffEqBase.DEAlgorithm,args...;
                                          kwargs...)
    solve(prob.prob,alg,args...;kwargs...)
end

function discretize end

function get_concrete_problem(prob, kwargs)
  tspan = get_concrete_tspan(prob, kwargs)
  u0 = get_concrete_u0(prob, tspan[1], kwargs)
  remake(prob; u0 = u0, tspan = tspan)
end

function get_concrete_problem(prob::DDEProblem, kwargs)
  tspan = get_concrete_tspan(prob, kwargs)

  u0 = get_concrete_u0(prob, tspan[1], kwargs)

  if prob.constant_lags isa Function
    constant_lags = prob.constant_lags(prob.p)
  else
    constant_lags = prob.constant_lags
  end

  remake(prob; u0 = u0, tspan = tspan, constant_lags = constant_lags)
end

function get_concrete_tspan(prob, kwargs)
  if prob.tspan isa Function
    tspan = prob.tspan(prob.p)
  elseif prob.tspan === (nothing, nothing)
    if haskey(kwargs, :tspan)
      tspan = kwargs[:tspan]
    else
      error("No tspan is set in the problem or chosen in the init/solve call")
    end
  else
    tspan = prob.tspan
  end

  tspan
end

function get_concrete_u0(prob, t0, kwargs)
  if eval_u0(prob.u0)
    u0 = prob.u0(prob.p, t0)
  elseif prob.u0 === nothing
    u0 = kwargs[:u0]
  else
    u0 = prob.u0
  end

  handle_distribution_u0(u0)
end

handle_distribution_u0(_u0) = _u0
eval_u0(u0::Function) = true
eval_u0(u0) = false

"""
$(SIGNATURES)

Check whether the values of `u0` and `tspan` are appropriate for use with
adaptive integrators and emit specific warnings if they are not.
"""
function adaptive_warn(u0,tspan)
  adaptive_integer_warn(tspan)
end

"""
$(SIGNATURES)

Emit a warning about incompatibility with adaptive integers if `tspan` contains
integers.
"""
function adaptive_integer_warn(tspan)
  if eltype(tspan) <: Integer
    @warn("Integer time values are incompatible with adaptive integrators. Utilize floating point numbers instead of integers in this case, i.e. (0.0,1.0) instead of (0,1).")
  end
end


text:
using Documenter
using Layered

makedocs(
    sitename = "Layered.jl",
    pages = [
        "Start" => "index.md",
        "Tutorial" => "tutorial.md",
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/jkrumbiegel/Layered.jl.git",
)

ext: jl


text:
"""
    module Utility

Some useful functions, mostly for converting stuff and loading files.
"""
module Utility

using Printf
"""
Some useful functions
"""

function cutoff_fn(num, min_c, max_c)
    if num < min_c
        return 1.0
    elseif num > max_c
        return 0.0
    else
        t = (num - min_c)/(max_c-min_c)
        return 1.0 - 10.0 * t^3 + 15.0 *  t^4  - 6.0 * t^5
    end
end

function arr2str(a::Array{Int,1})
    st=""
    for i = 1:size(a,1)
        t= @sprintf(" % 7s", a[i])
        st=st*t
    end
    return st

end

function arr2str(a::Array{Float64,1})
    st=""
    for i = 1:size(a,1)
        t= @sprintf(" % 2.16e", a[i])
        st=st*t
    end
    return st

end

#now matches all number arrays!!!
function arr2str(a::AbstractArray{<:Number,2})
    st=""
    for i = 1:size(a,1)
        t= arr2str(a[i,:])
        st=st*t*"\n"
    end
    return st

end

#now matches all number arrays!!!
function arr2str(a::AbstractArray{Any,2})
    st=""
    for i = 1:size(a,1)
        t = str_w_spaces(a[i,:])
        st=st*t*"\n"
    end
    return st

end

function str_w_spaces(a)
    st=""
    for i = 1:(size(a,1)-1)
        st= st*string(a[i])* "  "
    end
    st= st*string(a[size(a,1)])

    return st

end

function str_w_spaces(a::Set)
    st=""
    for s in a
        st= st*string(s)* "  "
    end
#    st= st*string(a[size(a,1)])

    return st

end

function parse_str_ARR_float(sp)

    return map(x->parse(Float64,x),split(sp))

end

function parse_str_ARR_complex(sp)

    t = map(x->parse(Float64,x),split(replace(sp, "," => " ")))
    return t[1:2:end] + im*t[2:2:end]

end

function dict2str(a::Dict)
    st=""
    for k in keys(a)
        st *= string(k) *" => "* string(a[k]) *"\n"
    end
    return st
end

function str2tuplesdict(st::String, d=missing)
    if ismissing(d)
        d = Dict()
    end

    for line in split(st, "\n")
        sp = split(line, " => ")
        if length(sp) == 2
            key = eval(Meta.parse(sp[1]))
            val = eval(Meta.parse(sp[2]))
            d[key] = val
        end
    end
    return d
end

function write_to_file(str, filename, directory="./")

    f=open("$directory/$filename", "w")
    write(f, str)
    close(f)

end

end


text:
module MyProject
export hello, addnums, mulnums
"""
    hello(who::String)

Return "Hello, `who`".
"""
hello(who::String) = "Hello, $who"

"""
    addnums(x::Number, y::Number)
Return `x + y`
"""
addnums(x::Number, y::Number) = x + y

"""
    mulnums(x::Number, y::Number)
Return `z`
"""
function mulnums(x::Number, y::Number)
  z = x * y
  return z
end

end


text:
using CFITSIO
using Test
using Aqua

function tempfitsfile(fn)
    mktempdir() do dir
        filename = joinpath(dir, "temp.fits")
        fitsfile = fits_clobber_file(filename)
        fn(fitsfile)

        if fitsfile.ptr != C_NULL
            # write some data to file to avoid errors on closing
            data = ones(1)
            fits_create_img(fitsfile, data)
            fits_write_pix(fitsfile, data)
            fits_delete_file(fitsfile)
        end
    end
end

@testset "project quality" begin
    Aqua.test_all(CFITSIO)
end

# `create_test_file` : Create a simple FITS file for testing, with the
# given header string added after the required keywords. The length of
# `header` must be a multiple of 80.  The purpose of creating such
# files is to test the parsing of non-standard FITS keyword records
# (non-standard files can't be created with cfitsio).

function create_test_file(fname::AbstractString, header::String)
    if length(header) % 80 != 0
        error("length of header must be multiple of 80")
    end

    f = open(fname, "w")

    stdhdr = "SIMPLE  =                    T / file does conform to FITS standard             BITPIX  =                  -64 / number of bits per data pixel                  NAXIS   =                    2 / number of data axes                            NAXIS1  =                   10 / length of data axis 1                          NAXIS2  =                   10 / length of data axis 2                          EXTEND  =                    T / FITS dataset may contain extensions            "
    endline = "END                                                                             "
    data = fill(0., (10, 10))  # 10x10 array of big-endian Float64 zeros

    # write header
    write(f, stdhdr)
    write(f, header)
    write(f, endline)

    # add padding
    block_position = (length(stdhdr) + length(header) + length(endline)) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, " "^padding)

    # write data
    write(f, data)

    # add padding
    block_position = sizeof(data) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, fill(0x00, (padding,)))

    close(f)
end

function writehealpix(filename, pixels, nside, ordering, coordsys)
    if eltype(pixels) == Float32
        tform = "1E"
    elseif eltype(pixels) == Float64
        tform = "1D"
    end

    file = fits_clobber_file(filename)
    try
        fits_create_img(file, Int16, Int[])
        fits_write_date(file)
        fits_movabs_hdu(file, 1)
        fits_create_binary_tbl(file, length(pixels), [("SIGNAL", tform, "")], "BINTABLE")
        fits_write_key(file, "PIXTYPE",  "HEALPIX", "HEALPIX pixelization")
        fits_write_key(file, "ORDERING", ordering,  "Pixel ordering scheme (either RING or NESTED)")
        fits_write_key(file, "NSIDE",    nside,     "Resolution parameter for HEALPIX")
        fits_write_key(file, "COORDSYS", coordsys,  "Pixelization coordinate system")
        fits_write_comment(file, "G = galactic, E = ecliptic, C = celestial = equatorial")
        fits_write_col(file, 1, 1, 1, pixels)
    finally
        fits_close_file(file)
    end
end

function readhealpix(filename)
    file = fits_open_file(filename)
    try
        hdutype = fits_movabs_hdu(file, 2)
        tform, tform_comment = fits_read_key_str(file, "TFORM1")
        if tform == "1E"
            T = Float32
        elseif tform == "1D"
            T = Float64
        end

        naxes, naxis_comment = fits_read_key_lng(file, "NAXIS")
        naxis, nfound = fits_read_keys_lng(file, "NAXIS", 1, naxes)
        nside, nside_comment = fits_read_key_lng(file, "NSIDE")
        npix = 12*nside*nside

        ordering, ordering_comment = fits_read_key_str(file, "ORDERING")
        coordsys, coordsys_comment = fits_read_key_str(file, "COORDSYS")

        pixels = zeros(T, npix)
        fits_read_col(file, 1, 1, 1, pixels)

        return pixels, nside, ordering, coordsys

    finally
        fits_close_file(file)
    end
end

@testset "CFITSIO.jl" begin

    @testset "file name and mode" begin
        mktempdir() do dir
            filename = joinpath(dir, "temp.fits")
            f = fits_clobber_file(filename)
            @test fits_file_name(f) == filename
            @test fits_file_mode(f) == 1

            # Write some data to the file
            a = ones(2,2)
            fits_create_img(f, eltype(a), [size(a)...])
            fits_write_pix(f, a)
            close(f)

            f = fits_open_file(filename, 0)
            @test fits_file_mode(f) == 0 == Int(CFITSIO.R)
            close(f)

            f = fits_open_file(filename, CFITSIO.R)
            @test fits_file_mode(f) == 0 == Int(CFITSIO.R)
            close(f)

            f = fits_open_file(filename, 1)
            @test fits_file_mode(f) == 1 == Int(CFITSIO.RW)
            close(f)

            f = fits_open_file(filename, CFITSIO.RW)
            @test fits_file_mode(f) == 1 == Int(CFITSIO.RW)
            close(f)
        end
    end

    @testset "types" begin
        for (T, code) in (
            (UInt8, 11),
            (Int8, 12),
            (Bool, 14),
            (String, 16),
            (Cushort, 20),
            (Cshort, 21),
            (Cuint, 30),
            (Cint, 31),
            (Int64, 81),
            (Float32, 42),
            (Float64, 82),
            (ComplexF32, 83),
            (ComplexF64, 163),
        )
            @test cfitsio_typecode(T) == Cint(code)
        end

        for (T, code) in ((UInt8,     8), # BYTE_IMG
                        (Int16,    16), # SHORT_IMG
                        (Int32,    32), # LONG_IMG
                        (Int64,    64), # LONGLONG_IMG
                        (Float32, -32), # FLOAT_IMG
                        (Float64, -64), # DOUBLE_IMG
                        (Int8,     10), # SBYTE_IMG
                        (UInt16,   20), # USHORT_IMG
                        (UInt32,   40)) # ULONG_IMG
            @test bitpix_from_type(T) == code
            @test type_from_bitpix(Cint(code)) == type_from_bitpix(Val(Cint(code))) == T
            @test type_from_bitpix(Int16(code)) == T
        end
        @test_throws MethodError type_from_bitpix(7)
        @test_throws MethodError type_from_bitpix("BITPIX")
    end

    # test reading/writing Healpix maps as FITS binary tables using the Libcfitsio interface
    mktempdir() do dir
        filename = joinpath(dir, "temp.fits")
        for T in (Float32, Float64)
            nside = 4
            npix = 12*nside*nside
            pixels = rand(T, npix)
            ordering = "NESTED"
            coordsys = "G"

            writehealpix(filename, pixels, nside, ordering, coordsys)
            @test readhealpix(filename) == (pixels, nside, ordering, coordsys)
        end
    end

    @testset "Miscellaneous" begin
        # test that this function works and returns the right type.
        @test typeof(libcfitsio_version()) === VersionNumber
        # test it parses a number as intended.
        @test libcfitsio_version(3.341)  === VersionNumber(3, 34, 1)
        @test libcfitsio_version(3.41f0) === VersionNumber(3, 41, 0)
    end

    @testset "hdu operations" begin
        tempfitsfile() do f

            # Create a few HDUs
            a = ones(2,2)
            fits_create_img(f, eltype(a), [size(a)...])
            fits_write_pix(f, a)

            a = ones(3,3)
            fits_create_img(f, eltype(a), [size(a)...])
            fits_write_pix(f, a)

            a = ones(4,4)
            fits_create_img(f, eltype(a), [size(a)...])
            fits_write_pix(f, a)

            @test fits_get_num_hdus(f) == 3

            for i in 1:3
                fits_movabs_hdu(f, i)
                @test fits_get_hdu_num(f) == i
                @test fits_get_hdu_type(f) == :image_hdu
            end

            for i = 1:2
                fits_movrel_hdu(f, -1)
                @test fits_get_hdu_num(f) == 3 - i
            end

            fits_movabs_hdu(f, 2)
            fits_delete_hdu(f)
            @test fits_get_num_hdus(f) == 2
            @test fits_get_hdu_num(f) == 2
            @test fits_get_img_size(f) == [4,4]

            fits_movabs_hdu(f, 1)
            fits_delete_hdu(f)
            @test fits_get_num_hdus(f) == 2
            @test fits_get_hdu_num(f) == 1
            @test fits_get_img_dim(f) == 0
        end

        @testset "insert image" begin
            tempfitsfile() do f
                a = ones(2,2); b = similar(a)
                fits_insert_img(f, a)
                fits_write_pix(f, a)
                fits_read_pix(f, b)
                @test b == a
                @test fits_get_num_hdus(f) == 1
                a .*= 2
                fits_insert_img(f, eltype(a), [size(a)...])
                fits_write_pix(f, a)
                fits_read_pix(f, b)
                @test b == a
                @test fits_get_num_hdus(f) == 2
                fits_movabs_hdu(f, 1)
                a .*= 2
                fits_insert_img(f, a)
                fits_write_pix(f, a)
                fits_read_pix(f, b)
                @test b == a
                @test fits_get_num_hdus(f) == 3
                # test that the HDU is added in the middle
                fits_movabs_hdu(f, 1)
                fits_read_pix(f, b)
                @test b == ones(2,2)
                fits_movabs_hdu(f, 3)
                fits_read_pix(f, b)
                @test b == ones(2,2) .* 2
                fits_movabs_hdu(f, 2)
                fits_read_pix(f, b)
                @test b == ones(2,2) .* 4
            end
        end
    end

    @testset "image type/size" begin
        tempfitsfile() do f
            a = ones(2,2)
            fits_create_img(f, eltype(a), [size(a)...])
            fits_write_pix(f, a)
            @test fits_get_img_dim(f) == 2
            @test fits_get_img_size(f) == [2, 2]
            @test fits_get_img_type(f) == -64
            @test fits_get_img_equivtype(f) == -64

            a = ones(Int64, 2,2)
            fits_create_img(f, eltype(a), [size(a)...])
            fits_write_pix(f, a)
            @test fits_get_img_dim(f) == 2
            @test fits_get_img_size(f) == [2, 2]
            @test fits_get_img_type(f) == 64
            @test fits_get_img_equivtype(f) == 64
        end
    end

    @testset "read/write subset" begin
        tempfitsfile() do f
            a = rand(10,10)
            b = similar(a)
            fits_create_img(f, a)
            fits_write_pix(f, a)

            fits_read_subset(f, [1,1], [size(a)...], [1,1], b)
            @test a == b

            fits_write_subset(f, [1,1], [size(a,1),4], a)
            c = zeros(eltype(a), size(a,1), 4)
            fits_read_subset(f, [1,1], [size(a,1),4], [1,1], c)
            @test a[:, 1:4] == c

            a[:,1] .= NaN
            fits_write_subset(f, [1,1], [size(a,1),4], a)
            c .= 0
            fits_read_subset(f, [1,1], [size(a,1),4], [1,1], 100.0, c)
            @test all(==(100), c[:, 1])
            @test c[:, 2:4] == a[:,2:4]

            fits_read_subset(f, [1,1], [size(a,1),4], [1,1], 0.0, c)
            @test all(isnan, c[:, 1])
            @test c[:, 2:4] == a[:,2:4]

            tempfitsfile() do f2
                a = rand(20,20)
                fits_create_img(f, a)
                fits_write_pix(f, a)
                fits_copy_image_section(f, f2, "1:20,1:10")
                b = similar(a, 20, 10)
                fits_read_pix(f2, b)
                close(f2)
                @test a[1:20, 1:10] == b
            end
        end
    end

    @testset "error message" begin
        mktempdir() do dir
            filename = joinpath(dir, "temp.fits")
            try
                f = fits_clobber_file(filename)
                fits_close_file(f)
            catch e
                @test e isa CFITSIO.CFITSIOError
                @test e isa Exception # bugfix test as CFITSIOError didn't subtype Exception in #3
                io = IOBuffer()
                Base.showerror(io, e)
                errstr = String(take!(io))
                @test occursin(r"Error code"i, errstr)
                @test occursin(r"Error message"i, errstr)
            finally
                rm(filename, force=true)
            end
        end
    end

    @testset "closed file errors" begin
        tempfitsfile() do f
            # write arbitrary data to the file
            a = ones(2,2)
            fits_create_img(f, eltype(a), [size(a)...])
            fits_write_pix(f, a)
            close(f)

            # check that closed files throw a julia error and don't segfault
            for fn in [fits_file_name, fits_file_mode, fits_get_hdrspace, fits_write_date,
                fits_hdr2str, fits_get_hdu_num, fits_get_hdu_type, fits_get_img_type,
                fits_get_img_equivtype, fits_get_img_dim, fits_get_num_cols, fits_get_num_hdus,
                fits_get_rowsize, fits_get_img_size, fits_get_img_size, fits_get_num_rows,
                fits_delete_hdu,
                ]

                @test_throws Exception fn(f)
            end

            for fn in [fits_read_key_str, fits_read_key_lng, fits_read_keyword, fits_write_comment,
                fits_write_history, fits_write_record, fits_delete_key, fits_movnam_hdu, fits_get_colnum,
                ]

                @test_throws Exception fn(f, "abc")
            end

            for fn in [fits_read_record, fits_read_keyn, fits_delete_record,
                fits_movabs_hdu, fits_movrel_hdu, fits_get_coltype,
                fits_get_eqcoltype, fits_read_tdim, ]

                @test_throws Exception fn(f, 1)
            end

            for fn in [fits_insert_rows, fits_delete_rows, fits_read_descript, CFITSIO.fits_write_null_img]
                @test_throws Exception fn(f, 1, 2)
            end

            for fn in [fits_write_pix, fits_read_pix]
                @test_throws Exception fn(f, a)
            end
            for fn in [fits_write_pix, fits_read_pix]
                @test_throws Exception fn(f, [1,1], length(a), a)
            end

            @test_throws Exception fits_read_pix(f, ones(Int, ndims(a)), length(a), zero(eltype(a)),  a)

            @test_throws Exception fits_read_keys_lng(f, "a", 1, 2)

            for fn in [fits_write_key, fits_update_key]
                @test_throws Exception fn(f, "a", 1, "b")
            end

            for fn in [fits_read_col, fits_write_col]
                @test_throws Exception fn(f, 1, 1, 1, ["abc"])
                @test_throws Exception fn(f, 1, 1, 1, ["abc", 1])
            end

            for fn in [fits_create_binary_tbl, fits_create_ascii_tbl]
                @test_throws Exception fn(f, 1, [("name", "3D", "c")], "extname")
            end

            @test_throws Exception fits_update_key(f, "a", 1.0, "b")
            @test_throws Exception fits_update_key(f, "a", nothing, "b")
            @test_throws Exception fits_write_tdim(f, 1, [1, 2])

            @test_throws Exception fits_read_subset(f, [1,1], [2,2], [1,1], a)
            @test_throws Exception fits_write_subset(f, [1,1], [2,2], a)

            @test_throws Exception fits_create_img(f, Int64, [2,3])

            tempfitsfile() do f2
                fits_create_img(f2, eltype(a), [size(a)...])
                fits_write_pix(f2, a)
                close(f2)

                @test_throws Exception fits_copy_image_section(f, f2, "1:2")
                @test_throws Exception fits_copy_image_section(f2, f, "1:2")
            end
        end
    end

    @testset "null values" begin
        tempfitsfile() do f

            # all values are nullified
            a = ones(2,2)
            nullarray = similar(a, UInt8)
            fits_create_img(f, eltype(a), [size(a)...])
            fits_write_pixnull(f, ones(Int, ndims(a)), length(a), a, 1.0)
            fits_read_pix(f, a)
            @test all(isnan, a)

            # one values is nullified
            a .= [1 2; 3 4]
            fits_write_pixnull(f, a, 3.0)
            fits_read_pix(f, a)
            @test isnan(a[2,1])
            @test !isnan(a[2,2]) && all(!isnan, a[1,:])

            # one values is nullified, but replaced while being read back in
            a .= [1 2; 3 4]
            fits_write_pixnull(f, a, 3.0)
            fits_read_pix(f, a, 3.0)
            @test !any(isnan, a)
            @test a[2,1] == 3.0

            fits_write_pixnull(f, a, 3.0)
            fits_read_pix(f, a, 0.0)
            @test isnan(a[2,1])
            fits_read_pix(f, a, C_NULL)
            @test isnan(a[2,1])

            # get the indices of null values
            nullarray .= 0
            fits_read_pixnull(f, a, nullarray)
            @test nullarray[2,1] == 1
            @test iszero(nullarray[2,2]) && all(iszero, nullarray[1,:])

            nullarray .= 0
            fits_read_pixnull(f, a, vec(nullarray))
            @test nullarray[2,1] == 1
            @test iszero(nullarray[2,2]) && all(iszero, nullarray[1,:])

            @test_throws Exception fits_read_pixnull(f, a, similar(nullarray, 1))

            # don't treat null values as special while writing out
            a .= [1 2; 3 4]
            fits_write_pixnull(f, a, C_NULL)
            fits_read_pix(f, a)
            @test !any(isnan, a)
            a .= [1 2; NaN 4]
            fits_write_pixnull(f, a, C_NULL)
            fits_read_pix(f, a)
            @test isnan(a[2,1])
            @test !isnan(a[2,2]) && all(!isnan, a[1,:])

            # test that tuples and vectors of pixels behave identically
            a .= [1 2; NaN 4]
            nullarray .= 0
            nullarray2 = similar(nullarray)
            fits_write_pixnull(f, a, C_NULL)
            b = similar(a)
            fits_read_pixnull(f, [1,1], length(b), b, nullarray)
            fits_read_pixnull(f, (1,1), length(b), b, nullarray2)
            @test nullarray == nullarray2

            # Read in data by replacing null values with the specified value
            for nullval in Any[7, 7.0], fpixel in Any[[1,1], (1,1)]
                b .= 0
                fits_read_pix(f, fpixel, length(b), nullval, b)
                @test b[2,1] == nullval
            end

            for fpixel in Any[[1,1], (1,1)]
                # Write data first by treating the value 2 as null
                fits_write_pixnull(f, fpixel, length(b), a, 2)
                b .= 0
                fits_read_pix(f, fpixel, length(b), b)
                @test isnan(b[1,2])
                @test isnan(b[2,1])
                # replace the null value by a specified one
                for nullval in Any[7, 7.0]
                    b .= 0
                    fits_read_pix(f, fpixel, length(b), nullval, b)
                    @test b[2,1] == nullval
                    @test b[1,2] == nullval
                end
            end

            # set a stretch of values to null (NaN in this case)
            # for the test we set all values to null
            fits_write_null_img(f, 1, length(a))
            fits_read_pix(f, a)
            @test all(isnan, a)
        end
    end

    @testset "empty file" begin
        tempfitsfile() do f
            a = zeros(2,2)
            @test_throws ErrorException fits_read_pix(f, a)
            @test_throws ErrorException fits_read_pix(f, a, 1)
            @test_throws ErrorException fits_read_pixnull(f, a, similar(a, UInt8))
        end
    end

    @testset "non-Int64 types" begin
        tempfitsfile() do f
            # write arbitrary data to the file
            a = ones(2,2)
            fits_create_img(f, eltype(a), [size(a)...])
            fits_write_pix(f, Int32[1,1], big(length(a)),  a)
            b = similar(a)
            fits_read_pix(f, Int32[1,1], Int128(length(b)), b)
            @test b == a
        end
    end

    @testset "resize image" begin
        tempfitsfile() do f
            a = [1 2; 3 4];
            fits_create_img(f, a);
            fits_write_pix(f, a);
            fits_resize_img(f, [3,3]);
            @test fits_get_img_size(f) == [3,3]
            fits_resize_img(f, (4,4));
            @test fits_get_img_size(f) == [4,4]
            fits_resize_img(f, (7,));
            @test fits_get_img_size(f) == [7]
            fits_resize_img(f, Float64);
            @test type_from_bitpix(fits_get_img_type(f)) == Float64
            fits_write_pix(f, Float64.(1:7))
            b = zeros(Float64, 7)
            fits_read_pix(f, b)
            @test b == 1:7
        end
    end

    @testset "tuples vs vectors" begin
        tempfitsfile() do f
            a = Float64[1 3; 2 4]
            b = similar(a); c = similar(a);

            @testset "create" begin
                fits_create_img(f, eltype(a), size(a))
                fits_write_pix(f, a)
                fits_read_pix(f, b)
                fits_create_img(f, eltype(a), [size(a)...])
                fits_write_pix(f, a)
                fits_read_pix(f, c)
                @test b == c
            end
            @testset "write" begin
                fits_write_pix(f, [1,1], length(a), a)
                fits_read_pix(f, b)
                fits_write_pix(f, (1,1), length(a), a)
                fits_read_pix(f, c)
                @test b == c
            end

            @testset "size" begin
                @test fits_get_img_size(f, Val(2)) == (2,2)
            end
            @testset "read" begin
                @testset "full image" begin
                    fits_read_pix(f, b)
                    @test b == a
                    # test that vectors and tuples of pixels behave identically
                    fits_read_pix(f, [1,1], length(b), b)
                    fits_read_pix(f, (1,1), length(c), c)
                    @test b == c
                end
                @testset "subset" begin
                    b .= 0
                    # test that vectors and tuples of pixels behave identically
                    fits_read_subset(f, (1,1), (2,1), (1,1), @view b[:,1])
                    fits_read_subset(f, [1,1], [2,1], [1,1], @view b[:,2])
                    @test @views b[:,1] == b[:,2]
                end
            end
            @testset "subset" begin
                fits_create_img(f, eltype(a), (size(a,1),))
                fits_write_subset(f, [1,1], [2,1], a)
                fits_read_pix(f, @view b[:,1])
                fits_write_subset(f, (1,1), (2,1), a)
                fits_read_pix(f, @view b[:, 2])
                @test @views b[:,1] == b[:,2]

                b .= 0
                fits_write_subset(f, [1,1], [2,1], a)
                fits_read_subset(f, [1,1], [2,1], [1,1], @view b[:,1])
                fits_read_subset(f, (1,1), (2,1), (1,1), @view b[:,2])
                @test @views b[:,1] == b[:,2]
                b .= 0
                fits_write_subset(f, [1,1], [2,1], replace(a, 2 => NaN))
                fits_read_subset(f, [1,1], [2,1], [1,1], 4, @view b[:,1])
                fits_read_subset(f, (1,1), (2,1), (1,1), 4, @view b[:,2])
                @test @views b[:,1] == b[:,2]
                @test b[2,1] == b[2,2] == 4
            end
        end
    end

    @testset "extended filename parser" begin
        filename = tempname()
        a = [1 3; 2 4]
        b = similar(a)
        try
            # for filenames that don't contain extended format specifications,
            # the diskfile functions are equivalent to the file ones
            for createfn in [fits_create_file, fits_create_diskfile]
                f = createfn(filename)
                fits_create_img(f, a)
                fits_write_pix(f, a)
                close(f)
                # Extended format: flipping the image
                f = fits_open_file(filename*"[*,-*]")
                fits_read_pix(f, b)
                @test a[:, [2,1]] == b
                close(f)
                f = fits_open_file(filename*"[-*,*]")
                fits_read_pix(f, b)
                @test a[[2,1], :] == b
                close(f)
                # without extended format
                f = fits_open_diskfile(filename)
                fits_read_pix(f, b)
                @test a == b
                close(f)
                rm(filename)
            end
        finally
            rm(filename, force = true)
        end

        # the diskfile functions may include [] in the filenames
        filename2 = filename * "[abc].fits"
        @test_throws CFITSIO.CFITSIOError fits_create_file(filename2)
        try
            f = fits_create_diskfile(filename2)
            fits_create_img(f, a)
            fits_write_pix(f, a)
            fits_read_pix(f, b)
            @test a == b
            close(f)
            f = fits_open_diskfile(filename2)
            fits_read_pix(f, b)
            @test a == b
            close(f)
            @test_throws CFITSIO.CFITSIOError fits_open_file(filename2)
        finally
            rm(filename2, force = true)
        end
    end

    @testset "stdout/stdin streams" begin
        # We redirect the output streams to a temp file to avoid cluttering the output
        # At present this doesn't work completely, as there is some output from fits_create_img
        # that is not captured
        mktemp() do _, io
            redirect_stdout(io) do
                for fname in ["-", "stdout.gz"]
                    f = fits_create_file(fname);
                    for a in Any[[1 2; 3 4], Float64[1 2; 3 4]]
                        b = similar(a)
                        fits_create_img(f, a)
                        fits_write_pix(f, a)
                        fits_read_pix(f, b)
                        @test a == b
                    end
                    close(f)
                end
            end
        end
    end
end


text:
export Custom

module Custom

using ..Nodes
using ..Ranger
using ..Geometry
using ..Rendering
using ..Math
using ..Nodes.Filters
using ..Events
using ..Particles

include("cross_node.jl")
include("aa_rectangle.jl")
include("vector_text_node.jl")
include("outlined_rectangle.jl")
include("outlined_triangle.jl")
include("outlined_circle.jl")
include("line_node.jl")
include("particle_node.jl")

include("anchor_node.jl")

include("zoom_node.jl")

end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
#=

Author: Benjamin Hemingway

This code adjusts the square matrix in the case it is not positive semi definite.

This code is based on the Matlab package nearestSPD version 1.1.0.0 by John D'Errico
=#
function spdMatrix(M::Array{T,2}) where T<:Real
  # symmetrize A into B
  B = (M + M')/2;
  # Compute the symmetric polar factor of B. Call it H.
  # Clearly H is itself SPD.
  out = LinearAlgebra.svd(B)
  H = out.Vt'*Diagonal(out.S)*out.Vt
  # get Ahat in the above formula
  Ahat = (B+H)/2
  # ensure symmetry
  Ahat = (Ahat + Ahat')/2;
  # test that Ahat is in fact PD. if it is not so, then tweak it just a bit.
  p = isposdef(Ahat)
  k = 0;
  while p==false
    k+=1
    # Ahat failed the chol test. It must have been just a hair off,
    # due to floating point trash, so it is simplest now just to
    # tweak by adding a tiny multiple of an identity matrix.
    mineig = minimum(eigvals(Ahat))
    Ahat += (-mineig*k.^2 + eps(mineig))*Matrix{Float64}(I,size(Ahat))
    p = isposdef(Ahat)
    if k>99
      error("Could not create Positive Definite Matrix")
    end
  end
  return Ahat
end


text:
using Snopt
using SparseArrays

model, obj = Dynamics.dubinscar_parallelpark

solver = Solver(model, obj, N=101)
n,m,N = get_sizes(solver)
X0 = line_trajectory(solver)
U0 = ones(m,N-1)

solver.opts.cost_tolerance = 1e-6
solver.opts.cost_tolerance_intermediate = 1e-4
solver.opts.constraint_tolerance = 1e-4
solver.opts.square_root = true
solver.opts.penalty_initial = 0.1
solver.opts.penalty_scaling = 20
res, stats = solve(solver,U0)
stats["runtime"]
_cost(solver, res)

method = :hermite_simpson
N,N_ = get_N(solver,method)
n,m = get_sizes(solver)

# Create results structure
results = DircolResults(n,m,solver.N,method)
var0 = DircolVars(X0,U0)
Z0 = var0.Z

# Generate the objective/constraint function and its gradients
usrfun = gen_usrfun(solver, results, method, grads=:quad)
J, c, ceq, grad_J, jacob_c, jacob_ceq, fail = usrfun(Z0)

NN = (n+m)N
nG, = get_nG(solver,method)
eval_f, eval_g, eval_grad_f, eval_jac_g = gen_usrfun_ipopt(solver,method)
x_L, x_U, g_L, g_U = get_bounds(solver,method)
P = length(g_L)  # Total number of constraints
g = zeros(P)
grad_f = zeros(NN)
rows = zeros(nG)
cols = zeros(nG)
vals = zeros(nG)

eval_f(Z0) ≈ J
eval_g(Z0,g)
g == ceq
eval_grad_f(Z0, grad_f)
grad_f == grad_J
eval_jac_g(Z0,:Structure,rows,cols,vals)
eval_jac_g(Z0,:vals,rows,cols,vals)
cJ = sparse(rows,cols,vals)
jacob_ceq == cJ

function eval_all(Z0)
    eval_f(Z0)
    eval_g(Z0,g)
    eval_grad_f(Z0,grad_f)
    eval_jac_g(Z0,:vals,rows,cols,vals)
end

@btime usrfun($Z0)
@btime eval_all($Z0)

# Set up the problem
lb,ub = get_bounds(solver,method)

options = Dict{String, Any}()
options["Derivative option"] = 0
options["Verify level"] = 1
options["Minor feasibility tol"] = solver.opts.constraint_tolerance
# options["Minor optimality  tol"] = solver.opts.eps_intermediate
options["Major optimality  tol"] = solver.opts.cost_tolerance

row,col = constraint_jacobian_sparsity(solver,method)
# prob = Snopt.createProblem(usrfun, Z0, lb, ub, iE=row, jE=col)
# prob = Snopt.createProblem(usrfun, Z0, lb, ub)
# prob.x = Z0
# t_eval = @elapsed z_opt, fopt, info = snopt(prob, options, start=start)

usrfun = gen_usrfun(solver, results, method, grads=:quad)
@time z_opt, fopt, info = snopt(usrfun, Z0, lb, ub, options)
var_opt = DircolVars(z_opt,n,m,N)

cost(solver, var_opt)
_cost(solver, res)

res_s, stats_s = solve_dircol(solver, X0, U0, nlp=:snopt)
stats_s["runtime"]
res_p, stats_p = solve_dircol(solver, X0, U0, nlp=:ipopt)
stats_p["runtime"]
stats["runtime"]
parse_snopt_summary()

p = plot()
plot_trajectory!(res_s.X,label="snopt")
plot_trajectory!(res_p.X,label="ipopt")
plot_trajectory!(res,label="iLQR")


text:
export VectorialCompactSARTSATrajectory

const VectorialCompactSARTSATrajectory = Trajectory{
    SARTSA,
    types,
    NamedTuple{RTSA,trace_types},
} where {types,trace_types<:Tuple{Vararg{<:Vector}}}

"""
    VectorialCompactSARTSATrajectory(; state_type = Int, action_type = Int, reward_type = Float32, terminal_type = Bool)

This function creates a [`VectorialTrajectory`](@ref) of [`RTSA`](@ref) fields. Here the **Compact** in the function name means that, `state` and `next_state`, `action` and `next_action` reuse a same vector underlying.

# Example

```julia-repl
julia> t = VectorialCompactSARTSATrajectory()
0-element Trajectory{(:state, :action, :reward, :terminal, :next_state, :next_action),Tuple{Int64,Int64,Float32,Bool,Int64,Int64},NamedTuple{(:reward, :terminal, :state, :action),Tuple{Array{Float32,1},Array{Bool,1},Array{Int64,1},Array{Int64,1}}}}

julia> push!(t, state=0, action=0)

julia> push!(t, reward=0.f0, terminal=false, state=1, action=1)

julia> t
1-element Trajectory{(:state, :action, :reward, :terminal, :next_state, :next_action),Tuple{Int64,Int64,Float32,Bool,Int64,Int64},NamedTuple{(:reward, :terminal, :state, :action),Tuple{Array{Float32,1},Array{Bool,1},Array{Int64,1},Array{Int64,1}}}}:
 (state = 0, action = 0, reward = 0.0, terminal = 0, next_state = 1, next_action = 1)

julia> push!(t, reward=1.f0, terminal=true, state=2, action=2)

julia> t
2-element Trajectory{(:state, :action, :reward, :terminal, :next_state, :next_action),Tuple{Int64,Int64,Float32,Bool,Int64,Int64},NamedTuple{(:reward, :terminal, :state, :action),Tuple{Array{Float32,1},Array{Bool,1},Array{Int64,1},Array{Int64,1}}}}:
 (state = 0, action = 0, reward = 0.0, terminal = 0, next_state = 1, next_action = 1)
 (state = 1, action = 1, reward = 1.0, terminal = 1, next_state = 2, next_action = 2)

julia> get_trace(t, :state, :action)
(state = [0, 1], action = [0, 1])

julia> get_trace(t, :next_state, :next_action)
(next_state = [1, 2], next_action = [1, 2])

julia> pop!(t)
1-element Trajectory{(:state, :action, :reward, :terminal, :next_state, :next_action),Tuple{Int64,Int64,Float32,Bool,Int64,Int64},NamedTuple{(:reward, :terminal, :state, :action),Tuple{Array{Float32,1},Array{Bool,1},Array{Int64,1},Array{Int64,1}}}}:
 (state = 0, action = 0, reward = 0.0, terminal = 0, next_state = 1, next_action = 1)
```
"""
function VectorialCompactSARTSATrajectory(;
    state_type = Int,
    action_type = Int,
    reward_type = Float32,
    terminal_type = Bool,
)
    VectorialCompactSARTSATrajectory{
        Tuple{state_type,action_type,reward_type,terminal_type,state_type,action_type},
        Tuple{
            Vector{reward_type},
            Vector{terminal_type},
            Vector{state_type},
            Vector{action_type},
        },
    }((
        reward = Vector{reward_type}(),
        terminal = Vector{terminal_type}(),
        state = Vector{state_type}(),
        action = Vector{action_type}(),
    ))
end


text:
using StatsBase, Plots, FreqTables

# Univariate
# h = fit(Histogram, rand(100))
# h = fit(Histogram, rand(100), 0:0.1:1.0)
h = fit(Histogram, rand(100), nbins=10)
# h = fit(Histogram, rand(100), weights(rand(100)), 0:0.1:1.0)
# h = fit(Histogram, [20], 0:20:100)
# h = fit(Histogram, [20], 0:20:100, closed=:left)

# Multivariate
# h = fit(Histogram, (rand(100),rand(100)))
# h = fit(Histogram, (rand(100),rand(100)),nbins=10)

display(h)
# plot(h)


text:
type DropoutLayer <: Layer
    inputActivation::Vol
    outputActivation::Vol
    keepProbability::Float64
    keptWeights::BitArray{3}

    function DropoutLayer(keepProbability::Float64, dimensions)
        new(Vol(dimensions), Vol(dimensions), keepProbability, BitArray(dimensions))
    end
end

function forward(self::DropoutLayer, inputActivation::Vol)
    self.inputActivation = inputActivation
    self.keptWeights = rand(size(self.inputActivation)) .< self.keepProbability
    self.outputActivation = self.inputActivation .* self.keptWeights
    return self.outputActivation
end

function backward(self::DropoutLayer)
    fill!(self.inputActivation.dw, 0.0)
    self.inputActivation.dw = self.outputActivation.dw .* self.keptWeights
end

gradientStates(self::DropoutLayer) = []


text:
using SVDD, OneClassActiveLearning, JSON, Logging, JuMP, Ipopt, Random

function run_experiment(experiment::Dict)
    Random.seed!(0)
    @info "host: $(gethostname()) worker: $(getpid()) exp_hash: $(experiment[:hash])"
    if isfile(experiment[:output_file])
        @warn "Aborting experiment because the output file already exists. Filename: $(experiment[:output_file])"
        return nothing
    end

    res = Result(experiment)
    errorfile = "$(experiment[:log_dir])worker/$(gethostname())_$(getpid())"
    try
        time_exp = @elapsed res = OneClassActiveLearning.active_learn(experiment)
        res.al_summary[:runtime] = Dict(:time_exp => time_exp)
    catch e
        res.status[:exit_code] = Symbol(typeof(e))
        @warn "Experiment $(experiment[:hash]) finished with unkown error."
        @warn e
    finally
        if res.status[:exit_code] != :success
            @info "Writing error hash to $errorfile.error."
            open("$errorfile.error", "a") do f
                print(f, "$(experiment[:hash])\n")
            end
        end
        @info "Writing result to $(experiment[:output_file])."
        OneClassActiveLearning.write_result_to_file(experiment[:output_file], res)
    end
end

const NUM_OBSERVATIONS = 157

experiment = Dict{Symbol, Any}(
    :hash => 1,
    :data_file => "$(@__DIR__)/dummy.csv",
    :data_set_name => "Dummy",
    :output_file => "$(@__DIR__)/dummy.json",
    :log_dir => "$(@__DIR__)/dummy.log",
    :split_strategy_name => "Sf",
    :initial_pool_strategy_name => "Pu",
    :model => Dict(:type => :SVDDneg,
                   :param => Dict{Symbol, Any}(),
                   :init_strategy => SimpleCombinedStrategy(RuleOfThumbScott(), BoundedTaxErrorEstimate(0.05, 0.02, 0.98))),
    :query_strategy => Dict(:type => :RandomOutlierQs,
                            :param => Dict{Symbol, Any}()),
    :split_strategy => OneClassActiveLearning.DataSplits(trues(NUM_OBSERVATIONS)),
    :oracle => :PoolOracle,
    :param => Dict(:num_al_iterations => 10,
                   :solver => with_optimizer(Ipopt.Optimizer; print_level=0),
                   :initial_pools => fill(:U, NUM_OBSERVATIONS),
                   :adjust_K => true,
                   :initial_pool_resample_version => 1))

run_experiment(experiment)


text:
## extensions/HBond.jl --- models for hydrogen bonding

export HBondDict, hbondenergy, hbonded
export donors, donorindices, acceptors, acceptorindices

export nturn, alphaturn, shortturn, piturn
export parallelbridge, antiparallelbridge

export recorddonor!, recordacceptor!

"""
    HBondDict{N,T<:AbstractFloat}

Hydrogen bonding dictionary for a specific residue. Type parameter `N` describes at
most how many hydrogen bonds can be donated/accepted by the amide part of the residue.
Index of the residue (in a given topology) is `this.handle`.

**Note**: Values of `N` that make sense are ``1`` and ``2``; any greater would not
quite make sense. DSSP uses two to store the H-bonds with lowest energy (most stable).
"""
struct HBondDict{N,T<:AbstractFloat} <: AbstractVector{Tuple{UInt64,T}}
    interacting_residues::MVector{N,UInt64}
    energies::MVector{N,T}
    base_size::Int

    function HBondDict(::Val{M}, ::Type{Tp}) where {M,Tp}
        inter = @MVector zeros(UInt64, M)
        energies = @MVector zeros(Tp, M)
        return new{M,Tp}(inter, energies, M >> 1)
    end
end
function HBondDict(base::Integer)
    !(0 < base) && error("try using at least 1 hydrogen bond")
    return HBondDict(Val(2*base), float(Int))
end
Base.length(::HBondDict{N}) where N = N
Base.size(::HBondDict{N}) where N = (N,)

function Base.iterate(dict::HBondDict{N}, state=1) where N
    state > N && return nothing
    (; interacting_residues, energies) = dict
    ## donors are stored first relative to acceptors
    tpl = @inbounds interacting_residues[state], energies[state]
    return (tpl, state+1)
end

Base.@propagate_inbounds function Base.getindex(dict::HBondDict, i::Integer)
    return (dict.interacting_residues[i], dict.energies[i])
end

Base.@propagate_inbounds function Base.setindex!(dict::HBondDict, val, i::Integer)
    res, energy = val
    dict.interacting_residues[i] = res
    dict.energies[i] = energy
    nothing
end

donors(dict::HBondDict{N}) where N = Iterators.take(dict, N>>1)
acceptors(dict::HBondDict{N}) where N = Iterators.drop(dict, N>>1)
donorindices(dict) = (i for (i, _) in donors(dict))
acceptorindices(dict) = (i for (i, _) in acceptors(dict))

## TODO: CREATE A VECTOR OF INITIALIZED HBONDDICT FROM FRAME INPUT

"""
    hbondenergy(frame, donor, acceptor; method=:dist)

Calculate the (amide) hydrogen bond energy between a donor residue and an
acceptor residue in a frame. Both nonframe arguments are _indices_ in the
`StructureFrame`. Defaults to a mimicked method of that _currently_ used
in DSSP (using distances, although there are plans to use angles).
An _approximate_ potential using dot products `method=:dot` is provided
in this module.

# Keyword arguments
+ `method`: Defaults to `:dist`. Method used for energy calculation.
"""
function hbondenergy(frame::StructureFrame, donor_i, acceptor_i; method=:dist)
    (; at_pos, res_list) = frame

    calcfunc = method == :dot ? (_hbenergy_dot) : (_hbenergy_dist)
    donor = res_list[donor_i]; acceptor = res_list[acceptor_i]
    ## get N and H in donor; C and O in acceptor
    Nid = donor[:N]; Hid = donor[:H]
    Cid = acceptor[:C]; Oid = acceptor[:O]

    Npos = at_pos[Nid]; Hpos = at_pos[Hid]
    Cpos = at_pos[Cid]; Opos = at_pos[Oid]

    return calcfunct(Npos, Hpos, Cpos, Opos)
end

function _hbenergy_dist(Npos, Hpos, Cpos, Opos)
    ## DSSP default implementation
    KK = DDSPConstants[:coupleConstant]     # coupling constant from dict
    HO = distance(Hpos, Opos); HC = distance(Hpos, Cpos)
    NO = distance(Npos, Opos); NC = distance(Npos, Cpos)

    ## coulomb interaction from two dipoles
    return KK * (1/HO - 1/HC + 1/NO - 1/NC)
end

function _hbenergy_dot(Npos, Hpos, Cpos, Opos)
    ## method using dot
    KK = DSSPConstants[:coupleConstant]     # coupling constant from dict
    ## direction of dipole moment
    NHp = Hpos - Npos; COp = Cpos - Opos
    ## centers of dipoles
    NHc = (Npos + Hpos) / 2; COc = (Cpos + Opos) / 2
    R = COc - NHc; r = norm(disp)

    ## "keesom" interaction of two dipoles
    pre = KK * (NHp⋅COp) / r^3 - 3 * (NHp⋅R) * (COp⋅R) / r^5
    return KK * pre
end

"""
    hbonded(hbdicts, res_i, res_j) -> Bool

Check if the residue at _index_ `res_i` acts as a hydrogen bond acceptor to
the residue at `res_j` based on the list of dictionaries `hbdicts`.
Equivalently, if residue at `res_j` acts as a hydrogen bond donor to
the residue at `res_i`.
"""
function hbonded(hbdicts::Vector{<:HBondDict}, i, j)
    ## annotate with @boundscheck??
    checkbounds(hbdicts, i); checkbounds(hbdicts, j)
    @inbounds i in acceptorindices(hbdicts[j])
end

## TODO: check for chain breaks between indices i and j
### best solution i can think of is a new data struct that is produced linearly
### relative to number of residues; and can be searched logarithmically
### relative to number of (separate) chains (tree-like)

"""
    nturn(hbdicts, i, n) -> Bool
    alphaturn(hbdicts, i) = nturn(hbdicts, i, 4)
    shortturn(hbdicts, i) = nturn(hbdicts, i, 3)
    piturn(hbdicts, i) = nturn(hbdicts, i, 5)

Check if there is an ``n``-turn starting from residue with index `i`. Uses
the vector `hbdicts` as reference. DSSP defines turns with specific `n` such
as ``3,4,5``-turns. These are `shortturn`, `alphaturn`, and `piturn`,
respectively, in this module.
"""
nturn(hbdicts, i, n) = hbonded(hbdicts, i, i+n)
@doc (@doc nturn) alphaturn(hbdicts, i) = nturn(hbdicts, i, 4)
@doc (@doc nturn) shortturn(hbdicts, i) = nturn(hbdicts, i, 3)
@doc (@doc nturn) piturn(hbdicts, i) = nturn(hbdicts, i, 5)

"""
    parallelbridge(hbdicts, i, j)
    antiparallelbridge(hbdicts, i, j)

Check if there is an elementary bridge between the residue with indices `i`
and `j`. Check the original DSSP reference (link in README) or the Julia
source codes for the definitions of these bridges.
"""
parallelbridge(hbdicts, i, j) = begin
    t1 = hbonded(hbdicts, i, j-1) && hbonded(hbdicts, j, i+1)
    t1 || hbonded(hbdicts, j-1, i) && hbonded(hbdicts, i, j+1)
end
@doc (@doc parallelbridge) antiparallelbridge(hbdicts, i, j) = begin
    t1 = hbonded(hbdicts, i, j) && hbonded(hbdicts, j, i)
    t1 || hbonded(hbdicts, i-1, j+1) && hbonded(hbdicts, j-1, i+1)
end

### ADDING DONORS/ACCEPTORS TO DICTIONARY

"""
    recorddonor!(hbdict, i, bondenergy)
    recordacceptor!(hbdict, i, bondenergy)

Record a donor/acceptor to the hydrogen bonding dictionary. The bond is with
the residue at index `i`. The corresponding bond energy is `bondenergy`.

Using these functions records entries with increasing energy (decreasing bond
strength).
"""
function recorddonor! end
@doc (@doc recorddonor!) function recordacceptor! end

tpes = map((:donor, :acceptor)) do nm
    Symbol(:record, nm, :!), Symbol(nm, :s)
end

for (fnm, called) in tpes
    offexp = fnm === :recorddonor! ? :(0) : :(hbdict.base_size)
    @eval begin
        function $fnm(hbdict::HBondDict, index, energy)
            offset = $offexp
            tr_ix, tr_energy = index, energy
            for (key, tpl) in $called(hbdict) |> enumerate
                ix, E = tpl
                tr_energy < E || continue
                # reaching here means energy < E; and partial roll starts
                @inbounds hbdict[key+offset] = tr_ix, tr_energy
                tr_ix, tr_energy = ix, E
            end
        end
    end
end


text:
"""
 Holds the tableau of a stochastic implicit split partitioned Runge-Kutta method.
 qdrift, pdrift1, pdrift2 hold the RK coefficients for the drift parts,
 and qdiff, pdiff1, pdiff2 hold the RK coefficients for the diffusion part of the SDE.
"""
struct TableauSISPRK{T} <: AbstractTableauIRK{T}
    name::Symbol
    s::Int
    qdrift ::CoefficientsRK{T}
    qdiff  ::CoefficientsRK{T}
    pdrift1::CoefficientsRK{T}
    pdrift2::CoefficientsRK{T}
    pdiff1 ::CoefficientsRK{T}
    pdiff2 ::CoefficientsRK{T}

    # Order of the tableau is not included, because unlike in the deterministic
    # setting, it depends on the properties of the noise (e.g., the dimension of
    # the Wiener process and the commutativity properties of the diffusion matrix)
    #
    # Orders stored in qdrift and qdiff are understood as the classical orders of these methods.

    function TableauSISPRK{T}(name, s, qdrift, qdiff, pdrift1, pdrift2, pdiff1, pdiff2) where {T}
        # THE COMMENTED OUT PART WAS FOR TableauFIRK. MAY IMPLEMENT SOMETHING
        # SIMILAR FOR TableauSFIRK LATER.

        # if (q.s > 1 && istrilstrict(q.a)) || (q.s==1 && q.a[1,1] == 0)
        #     warn("Initializing TableauFIRK with explicit tableau ", q.name, ".\n",
        #          "You might want to use TableauERK instead.")
        # elseif q.s > 1 && istril(q.a)
        #     warn("Initializing TableauFIRK with diagonally implicit tableau ", q.name, ".\n",
        #          "You might want to use TableauDIRK instead.")
        # end

        @assert s == qdrift.s == qdiff.s == pdrift1.s == pdrift2.s == pdiff1.s == pdiff2.s

        new(name, s, qdrift, qdiff, pdrift1, pdrift2, pdiff1, pdiff2)
    end
end

function TableauSISPRK(name::Symbol, qdrift::CoefficientsRK{T},
                                     qdiff::CoefficientsRK{T},
                                     pdrift1::CoefficientsRK{T}, pdrift2::CoefficientsRK{T},
                                     pdiff1::CoefficientsRK{T}, pdiff2::CoefficientsRK{T}) where {T}
    TableauSISPRK{T}(name, qdrift.s, qdrift, qdiff, pdrift1, pdrift2, pdiff1, pdiff2)
end

function TableauSISPRK(name::Symbol, qorder_drift::Int, qa_drift::Matrix{T}, qb_drift::Vector{T}, qc_drift::Vector{T},
                                      qorder_diff::Int, qa_diff::Matrix{T}, qb_diff::Vector{T}, qc_diff::Vector{T},
                                      porder1_drift::Int, pa1_drift::Matrix{T}, pb1_drift::Vector{T}, pc1_drift::Vector{T},
                                      porder2_drift::Int, pa2_drift::Matrix{T}, pb2_drift::Vector{T}, pc2_drift::Vector{T},
                                      porder1_diff::Int, pa1_diff::Matrix{T}, pb1_diff::Vector{T}, pc1_diff::Vector{T},
                                      porder2_diff::Int, pa2_diff::Matrix{T}, pb2_diff::Vector{T}, pc2_diff::Vector{T}) where {T}
    TableauSISPRK{T}(name, length(qc_drift), CoefficientsRK(name, qorder_drift, qa_drift, qb_drift, qc_drift),
                                              CoefficientsRK(name, qorder_diff, qa_diff, qb_diff, qc_diff),
                                              CoefficientsRK(name, porder1_drift, pa1_drift, pb1_drift, pc1_drift),
                                              CoefficientsRK(name, porder2_drift, pa2_drift, pb2_drift, pc2_drift),
                                              CoefficientsRK(name, porder1_diff, pa1_diff, pb1_diff, pc1_diff),
                                              CoefficientsRK(name, porder2_diff, pa2_diff, pb2_diff, pc2_diff))
end

# TODO function readTableauSFIRKFromFile(dir::AbstractString, name::AbstractString)

"""
Parameters for right-hand side function of implicit Runge-Kutta methods.
  A - if positive, the upper bound of the Wiener process increments; if A=0.0, no truncation
"""
mutable struct ParametersSISPRK{DT, TT, ET <: SPSDE{DT,TT}, D, M, S} <: Parameters{DT,TT}
    equ::ET
    tab::TableauSISPRK{TT}
    Δt::TT
    ΔW::Vector{DT}
    ΔZ::Vector{DT}
    A::DT

    t::TT
    q::Vector{DT}
    p::Vector{DT}
end

function ParametersSISPRK(equ::ET, tab::TableauSISPRK{TT}, Δt::TT, ΔW::Vector{DT}, ΔZ::Vector{DT}, A::DT) where {DT, TT, ET <: SPSDE{DT,TT}}
    @assert equ.m == length(ΔW) == length(ΔZ)
    ParametersSISPRK{DT, TT, ET, equ.d, equ.m, tab.s}(equ, tab, Δt, ΔW, ΔZ, A, 0, zeros(DT, equ.d), zeros(DT, equ.d))
end

"""
Structure for holding the internal stages Q, the values of the drift vector
and the diffusion matrix evaluated at the internal stages V=v(Q), B=B(Q),
and the increments Y = Δt*a_drift*v(Q) + a_diff*B(Q)*ΔW
"""
struct NonlinearFunctionCacheSISPRK{DT}
    Q ::Vector{Vector{DT}}
    P ::Vector{Vector{DT}}
    V ::Vector{Vector{DT}}
    F1::Vector{Vector{DT}}
    F2::Vector{Vector{DT}}
    B ::Vector{Matrix{DT}}
    G1::Vector{Matrix{DT}}
    G2::Vector{Matrix{DT}}
    Y ::Vector{Vector{DT}}
    Z ::Vector{Vector{DT}}

    v ::Vector{DT}
    f1::Vector{DT}
    f2::Vector{DT}
    b ::Matrix{DT}
    g1::Matrix{DT}
    g2::Matrix{DT}
    y ::Vector{DT}
    z ::Vector{DT}

    function NonlinearFunctionCacheSISPRK{DT}(D, M, S) where {DT}

        # create internal stage vectors
        Q  = create_internal_stage_vector(DT, D, S)
        P  = create_internal_stage_vector(DT, D, S)
        V  = create_internal_stage_vector(DT, D, S)
        F1 = create_internal_stage_vector(DT, D, S)
        F2 = create_internal_stage_vector(DT, D, S)
        B  = create_internal_stage_vector(DT, D, M, S)
        G1 = create_internal_stage_vector(DT, D, M, S)
        G2 = create_internal_stage_vector(DT, D, M, S)
        Y  = create_internal_stage_vector(DT, D, S)
        Z  = create_internal_stage_vector(DT, D, S)

        # create velocity and update vector
        v  = zeros(DT,D)
        f1 = zeros(DT,D)
        f2 = zeros(DT,D)
        b  = zeros(DT,D,M)
        g1 = zeros(DT,D,M)
        g2 = zeros(DT,D,M)
        y  = zeros(DT,D)
        z  = zeros(DT,D)

        new(Q, P, V, F1, F2, B, G1, G2, Y, Z, v, f1, f2, b, g1, g2, y, z)
    end
end

"""
Unpacks the data stored in x = (Y[1][1], Y[1][2], ... Y[1][D], Y[2][1], ..., Z[1][1], Z[1][2], ... Z[1][D], Z[2][1], ...)
into Y, Z::Vector{Vector}, calculates the internal stages Q, P, the values of the RHS
of the SDE ( vi(Q,P), fi(Q,P), Bi(Q,P) and Gi(Q,P) ), and assigns them to V[i], F[i], B[i] and G[i].
Unlike for FIRK, here
Y = Δt a_drift v(Q,P) + a_diff B(Q,P) ΔW,
Z = Δt ̃a1_drift f1(Q,P) + Δt ̃a2_drift f2(Q,P) + ̃a1_diff G1(Q,P) ΔW + ̃a2_diff G2(Q,P) ΔW.
"""
function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, P::Vector{Vector{ST}},
                                        V::Vector{Vector{ST}}, F1::Vector{Vector{ST}}, F2::Vector{Vector{ST}},
                                        B::Vector{Matrix{ST}}, G1::Vector{Matrix{ST}}, G2::Vector{Matrix{ST}},
                                        Y::Vector{Vector{ST}}, Z::Vector{Vector{ST}},
                                        params::ParametersSISPRK{DT,TT,ET,D,M,S}) where {ST,DT,TT,ET,D,M,S}

    local tqᵢ::TT       #times for the q internal stages
    local tpᵢ::TT       #times for the p internal stages

    @assert S == length(Q) == length(V) == length(B) == length(P) == length(F1) == length(F2) == length(G1) == length(G2)

    # copy x to Y, Z and calculate Q, P
    for i in 1:S
        @assert D == length(Q[i]) == length(V[i]) == size(B[i],1) == length(P[i]) == length(F1[i]) == length(F2[i]) == size(G1[i],1) == size(G2[i],1)
        @assert M == size(B[i],2) == size(G1[i],2) == size(G2[i],2)

        for k in 1:D
            Y[i][k] = x[D*(  i-1)+k]
            Z[i][k] = x[D*(S+i-1)+k]
        end

        Q[i] .= params.q .+ Y[i]
        P[i] .= params.p .+ Z[i]
    end

    # compute Vi = vi(Q) and Bi=Bi(Q)
    for i in 1:S
        tqᵢ = params.t + params.Δt * params.tab.qdrift.c[i]
        # Not sure what about pdrift2.c[i] --- should there be another time point?
        tpᵢ = params.t + params.Δt * params.tab.pdrift1.c[i]
        # calculates v(t,Q,P), f1(t,Q,P), f2(t,Q,P) and assigns to V, F1, F2
        params.equ.v(tqᵢ, Q[i], P[i], V[i])
        params.equ.f1(tpᵢ, Q[i], P[i], F1[i])
        params.equ.f2(tpᵢ, Q[i], P[i], F2[i])
        # calculates B(t,Q,P), G1(t,Q,P), G2(t,Q,P) and assigns to the matrices B[i], G1[i], G2[i]
        params.equ.B(tqᵢ, Q[i], P[i], B[i])
        params.equ.G1(tpᵢ, Q[i], P[i], G1[i])
        params.equ.G2(tpᵢ, Q[i], P[i], G2[i])
    end
end

"Compute stages of stochastic implicit split partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersSISPRK{DT,TT,ET,D,M,S}) where {ST,DT,TT,ET,D,M,S}

    cache = NonlinearFunctionCacheSISPRK{ST}(D, M, S)

    quote
        compute_stages!(x, $cache.Q, $cache.P, $cache.V, $cache.F1, $cache.F2, $cache.B, $cache.G1, $cache.G2, $cache.Y, $cache.Z, params)

        local y1::ST
        local y2::ST
        local z1::ST
        local z2::ST

        # compute b = - (Y-AV)
        for i in 1:S
            for k in 1:D
                y1 = 0
                y2 = 0
                z1 = 0
                z2 = 0
                for j in 1:S
                    y1 += params.tab.qdrift.a[i,j]  * $cache.V[j][k]  * params.Δt + params.tab.qdiff.a[i,j]  * dot($cache.B[j][k,:], params.ΔW)
                    y2 += params.tab.qdrift.â[i,j]  * $cache.V[j][k]  * params.Δt + params.tab.qdiff.â[i,j]  * dot($cache.B[j][k,:], params.ΔW)
                    z1 += params.tab.pdrift1.a[i,j] * $cache.F1[j][k] * params.Δt + params.tab.pdrift2.a[i,j] * $cache.F2[j][k] * params.Δt
                        + params.tab.pdiff1.a[i,j]  * dot($cache.G1[j][k,:], params.ΔW) + params.tab.pdiff2.a[i,j] * dot($cache.G2[j][k,:], params.ΔW)
                    z2 += params.tab.pdrift1.â[i,j] * $cache.F1[j][k] * params.Δt + params.tab.pdrift2.â[i,j] * $cache.F2[j][k] * params.Δt
                        + params.tab.pdiff1.â[i,j]  * dot($cache.G1[j][k,:], params.ΔW) + params.tab.pdiff2.â[i,j] * dot($cache.G2[j][k,:], params.ΔW)
                end
                b[D*(  i-1)+k] = - $cache.Y[i][k] + (y1 + y2)
                b[D*(S+i-1)+k] = - $cache.Z[i][k] + (z1 + z2)
            end
        end
    end
end

"Stochastic implicit partitioned Runge-Kutta integrator."
struct IntegratorSISPRK{DT, TT, PT <: ParametersSISPRK{DT,TT},
                                ST <: NonlinearSolver{DT}, N} <: StochasticIntegrator{DT,TT}
    params::PT
    solver::ST
    # InitialGuessPSDE not implemented for SFIPRK
    #iguess::IT
    fcache::NonlinearFunctionCacheSISPRK{DT}

    q::Matrix{Vector{TwicePrecision{DT}}}
    p::Matrix{Vector{TwicePrecision{DT}}}
end

# K - the integer in the bound A = √(2 K Δt |log Δt|) due to Milstein & Tretyakov; K=0 no truncation
function IntegratorSISPRK(equation::SPSDE{DT,TT,VT,F1T,F2T,BT,G1T,G2T,N}, tableau::TableauSISPRK{TT}, Δt::TT; K::Int=0) where {DT,TT,VT,F1T,F2T,BT,G1T,G2T,N}
    D = equation.d
    M = equation.m
    NS= equation.ns
    NI= equation.n
    S = tableau.s

    # create params
    K==0 ? A = 0.0 : A = sqrt( 2*K*Δt*abs(log(Δt)) )
    params = ParametersSISPRK(equation, tableau, Δt, zeros(DT,M), zeros(DT,M), A)

    # create solver
    solver = create_nonlinear_solver(DT, 2*D*S, params)

    # Not implementing InitialGuessSDE
    # create initial guess
    #iguess = InitialGuessODE(get_config(:ig_interpolation), equation, Δt)

    # create cache for internal stage vectors and update vectors
    fcache = NonlinearFunctionCacheSISPRK{DT}(D, M, S)

    # create solution vectors
    q = create_solution_vector(DT, D, NS, NI)
    p = create_solution_vector(DT, D, NS, NI)

    # create integrator
    IntegratorSISPRK{DT, TT, typeof(params), typeof(solver), N}(params, solver, fcache, q, p)
end

equation(integrator::IntegratorSISPRK) = integrator.params.equ
timestep(integrator::IntegratorSISPRK) = integrator.params.Δt
tableau(integrator::IntegratorSISPRK) = integrator.params.tab
Base.eltype(integrator::IntegratorSISPRK{DT, TT, PT, ST, N}) where {DT, TT, PT, ST, N} = DT

function initialize!(int::IntegratorSISPRK{DT,TT}, sol::SolutionPSDE, k::Int, m::Int) where {DT,TT}
    check_solution_dimension_asserts(sol, k, m)

    # copy the m-th initial condition for the k-th sample path
    get_initial_conditions!(sol, int.q[k,m], int.p[k,m], k, m)

    # Not implementing InitialGuessSDE
    # # initialise initial guess
    # initialize!(int.iguess, m, sol.t[0], int.q[m])
end

"""
This function computes initial guesses for Y, Z and assigns them to int.solver.x
The prediction is calculated using an explicit integrator.
"""
function initial_guess!(int::IntegratorSISPRK{DT,TT}) where {DT,TT}

    # NOT IMPLEMENTING InitialGuessSDE

    # SIMPLE SOLUTION
    # The simplest initial guess for Y, Z is 0
    int.solver.x .= zeros(eltype(int), 2*int.params.tab.s*ndims(int))

    # FIX NEEDED !!!
    # APPARENTLY PREDICTION USING AN EXPLICIT METHOD DOESNT WORK WELL HERE
    # PROBABLY BECAUSE THE pdrift2.c[i] POINTS ARE NOT INCLUDED

    # # USING AN EXPLICIT INTEGRATOR TO COMPUTE AN INITIAL GUESS
    # # Below we use the R2 method of Burrage & Burrage to calculate
    # # the internal stages at the times c[1]...c[s].
    # # This approach seems to give very good approximations if the time step
    # # and magnitude of noise are not too large. If the noise intensity is too big,
    # # one may have to perform a few iterations of the explicit method with a smaller
    # # time step, use a higher-order explicit method (e.g. CL or G5), or use
    # # the simple solution above.
    #
    # local tV1  = zeros(DT,int.params.equ.d)
    # local tV2  = zeros(DT,int.params.equ.d)
    # local tF1  = zeros(DT,int.params.equ.d)
    # local tF2  = zeros(DT,int.params.equ.d)
    # local tmpF = zeros(DT,int.params.equ.d)
    # local tB1  = zeros(DT,int.params.equ.d, int.params.equ.m)
    # local tB2  = zeros(DT,int.params.equ.d, int.params.equ.m)
    # local tG1  = zeros(DT,int.params.equ.d, int.params.equ.m)
    # local tG2  = zeros(DT,int.params.equ.d, int.params.equ.m)
    # local tmpG = zeros(DT,int.params.equ.d, int.params.equ.m)
    # local Q    = zeros(DT,int.params.equ.d)
    # local P    = zeros(DT,int.params.equ.d)
    # local t2::TT
    # local Δt_local::TT
    # local ΔW_local = zeros(DT,int.params.equ.m)
    #
    # # When calling this function, int.params should contain the data:
    # # int.params.q - the q solution at the previous time step
    # # int.params.p - the p solution at the previous time step
    # # int.params.t - the time of the previous step
    # # int.params.ΔW- the increment of the Brownian motion for the current step
    #
    # #Evaluating the functions v and B at t,q - same for all stages
    # int.params.equ.v(int.params.t, int.params.q, int.params.p, tV1)
    # int.params.equ.B(int.params.t, int.params.q, int.params.p, tB1)
    # int.params.equ.f1(int.params.t, int.params.q, int.params.p, tF1)
    # int.params.equ.f2(int.params.t, int.params.q, int.params.p, tmpF)
    # #tF1 = tF1 + tmpF
    # int.params.equ.G1(int.params.t, int.params.q, int.params.p, tG1)
    # int.params.equ.G2(int.params.t, int.params.q, int.params.p, tmpG)
    # #tG1 = tG1 + tmpG
    #
    #
    # # Calculating the positions q at the points qdrift.c[i]
    # # if qdrift.c==pdrift.c, then also calculating the momenta
    # for i in 1:int.params.tab.s
    #
    #     # Taking the c[i] from the qdrift tableau.
    #     Δt_local  = int.params.tab.qdrift.c[i]*int.params.Δt
    #     ΔW_local .= int.params.tab.qdrift.c[i]*int.params.ΔW
    #
    #     Q = int.params.q + 2./3. * Δt_local * tV1 + 2./3. * tB1 * ΔW_local
    #     P = int.params.p + 2./3. * Δt_local * tF1 + 2./3. * tG1 * ΔW_local
    #
    #     t2 = int.params.t + 2./3.*Δt_local
    #
    #     int.params.equ.v(t2, Q, P, tV2)
    #     int.params.equ.B(t2, Q, P, tB2)
    #     int.params.equ.f1(t2, Q, P, tF2)
    #     int.params.equ.f2(t2, Q, P, tmpF)
    #     #tF2 = tF2 + tmpF
    #     int.params.equ.G1(t2, Q, P, tG2)
    #     int.params.equ.G2(t2, Q, P, tmpG)
    #     #tG2 = tG2 + tmpG
    #
    #     #Calculating the Y's and assigning them to the array int.solver.x as initial guesses
    #     for j in 1:int.params.equ.d
    #         int.solver.x[(i-1)*int.params.equ.d+j] =  Δt_local*(1./4.*tV1[j] + 3./4.*tV2[j]) + dot( (1./4.*tB1[j,:] + 3./4.*tB2[j,:]), ΔW_local )
    #     end
    #
    #     # if the collocation points are the same for both q and p parts
    #     # not sure what about pdrift2.c
    #     if int.params.tab.qdrift.c==int.params.tab.pdrift1.c
    #         for j in 1:int.params.equ.d
    #             int.solver.x[(int.params.tab.s+i-1)*int.params.equ.d+j] =  Δt_local*(1./4.*tF1[j] + 3./4.*tF2[j]) + dot( (1./4.*tG1[j,:] + 3./4.*tG2[j,:]), ΔW_local )
    #         end
    #     end
    # end
    #
    #
    # # If qdrift.c != pdrift1.c, then calculating the momenta p at the points pdrift1.c[i]
    # # not sure what about pdrift2.c
    # if int.params.tab.qdrift.c != int.params.tab.pdrift1.c
    #
    #     for i in 1:int.params.tab.s
    #
    #         # Taking the c[i] from the pdrift1 tableau.
    #         Δt_local  = int.params.tab.pdrift1.c[i]*int.params.Δt
    #         ΔW_local .= int.params.tab.pdrift1.c[i]*int.params.ΔW
    #
    #         Q = int.params.q + 2./3. * Δt_local * tV1 + 2./3. * tB1 * ΔW_local
    #         P = int.params.p + 2./3. * Δt_local * tF1 + 2./3. * tG1 * ΔW_local
    #
    #         t2 = int.params.t + 2./3.*Δt_local
    #
    #         int.params.equ.v(t2, Q, P, tV2)
    #         int.params.equ.B(t2, Q, P, tB2)
    #         int.params.equ.f1(t2, Q, P, tF2)
    #         int.params.equ.f2(t2, Q, P, tmpF)
    #         #tF2 = tF2 + tmpF
    #         int.params.equ.G1(t2, Q, P, tG2)
    #         int.params.equ.G2(t2, Q, P, tmpG)
    #         #tG2 = tG2 + tmpG
    #
    #         # Calculating the Z's and assigning them to the array int.solver.x as initial guesses
    #         # The guesses for the Y's have already been written to x above
    #         for j in 1:int.params.equ.d
    #             int.solver.x[(int.params.tab.s+i-1)*int.params.equ.d+j] =  Δt_local*(1./4.*tF1[j] + 3./4.*tF2[j]) + dot( (1./4.*tG1[j,:] + 3./4.*tG2[j,:]), ΔW_local )
    #         end
    #     end
    # end

end

"""
Integrate PSDE with a stochastic implicit partitioned Runge-Kutta integrator.
 Integrating the k-th sample path for the m-th initial condition
"""
function integrate_step!(int::IntegratorSISPRK{DT,TT}, sol::SolutionPSDE{DT,TT,NQ,NW}, k::Int, m::Int, n::Int) where {DT,TT,NQ,NW}
    check_solution_dimension_asserts(sol, k, m, n)

    # set time for nonlinear solver
    int.params.t  = sol.t[0] + (n-1)*int.params.Δt
    int.params.q .= int.q[k, m]
    int.params.p .= int.p[k, m]

    # copy the increments of the Brownian Process
    if NW==1
        #1D Brownian motion, 1 sample path
        int.params.ΔW[1] = sol.W.ΔW[n-1]
        int.params.ΔZ[1] = sol.W.ΔZ[n-1]
    elseif NW==2
        #Multidimensional Brownian motion, 1 sample path
        int.params.ΔW .= sol.W.ΔW[n-1]
        int.params.ΔZ .= sol.W.ΔZ[n-1]
    elseif NW==3
        #1D or Multidimensional Brownian motion, k-th sample path
        int.params.ΔW .= sol.W.ΔW[n-1,k]
        int.params.ΔZ .= sol.W.ΔZ[n-1,k]
    end

    # truncate the increments ΔW with A
    if int.params.A>0
        for i in eachindex(int.params.ΔW)
            if int.params.ΔW[i]<-int.params.A
                int.params.ΔW[i] = -int.params.A
            elseif int.params.ΔW[i]>int.params.A
                int.params.ΔW[i] = int.params.A
            end
        end
    end

    # compute initial guess and assign to int.solver.x
    initial_guess!(int)

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute the drift vector field and the diffusion matrix at internal stages
    compute_stages!(int.solver.x, int.fcache.Q, int.fcache.P, int.fcache.V, int.fcache.F1, int.fcache.F2, int.fcache.B, int.fcache.G1, int.fcache.G2, int.fcache.Y, int.fcache.Z, int.params)

    # compute final update
    update_solution!(int.q[k,m], int.p[k,m], int.fcache.V, int.fcache.F1, int.fcache.F2, int.fcache.B, int.fcache.G1, int.fcache.G2,
                    int.params.tab.qdrift.b, int.params.tab.qdrift.b̂, int.params.tab.qdiff.b, int.params.tab.qdiff.b̂,
                    int.params.tab.pdrift1.b, int.params.tab.pdrift1.b̂, int.params.tab.pdrift2.b, int.params.tab.pdrift2.b̂,
                    int.params.tab.pdiff1.b, int.params.tab.pdiff1.b̂, int.params.tab.pdiff2.b, int.params.tab.pdiff2.b̂,
                    int.params.Δt, int.params.ΔW)

    # # NOT IMPLEMENTING InitialGuessSDE
    # # # copy solution to initial guess
    # # update!(int.iguess, m, sol.t[0] + n*int.params.Δt, int.q[m])

    # take care of periodic solutions
    cut_periodic_solution!(int.q[k,m], int.params.equ.periodicity)

    # # copy to solution
    set_solution!(sol, int.q[k,m], int.p[k,m], n, k, m)
end


text:
# specific stuff for quad trees

mutable struct QuadTreeNode{T<:AbstractPoint2D} <: SpatialTreeNode
    id::Int64
    r::Float64
    midx::Float64
    midy::Float64
    is_empty::Bool
    is_divided::Bool
    point::T
    lxly::QuadTreeNode{T}
    lxhy::QuadTreeNode{T}
    hxly::QuadTreeNode{T}
    hxhy::QuadTreeNode{T}
    function QuadTreeNode{T}(r::Number, midx::Number, midy::Number) where T <:AbstractPoint2D
        n = new(0, r, midx, midy, true, false, T())
        n.lxly = n
        n.lxhy = n
        n.hxly = n
        n.hxhy = n
        n
    end
end

QuadTreeNode(r::Number, midx::Number, midy::Number, ::Type{T}) where T<:AbstractPoint2D =  QuadTreeNode{T}(r, midx, midy)
QuadTreeNode(::Type{T}) where T<:AbstractPoint2D = QuadTreeNode(0.5, 0.5, 0.5, T)
QuadTreeNode() = QuadTreeNode(Point2D);

mutable struct QuadTree{T<:AbstractPoint2D} <: SpatialTree
	head::QuadTreeNode{T}
	number_of_nodes_used::Int64
	nodes::Array{QuadTreeNode, 1}
    faststack::Array{QuadTreeNode{T}, 1}
    function QuadTree{T}(r::Number, midx::Number, midy::Number, n::Int64=100000) where T <:AbstractPoint2D
    	nodes = QuadTreeNode[QuadTreeNode(r, midx, midy, T) for i in 1:n]
        new(nodes[1], 1, nodes, [QuadTreeNode(T) for i in 1:10000])
    end
end

QuadTree(r::Number, midx::Number, midy::Number, ::Type{T};n=1000) where T<:AbstractPoint2D = QuadTree{T}(r, midx, midy, n)
QuadTree(::Type{T};n=10000) where T <: AbstractPoint2D = QuadTree(0.5, 0.5, 0.5, T; n=n)
QuadTree(n::Int64) = QuadTree(Point2D;n=n)
QuadTree() = QuadTree(Point2D)
Base.eltype(::QuadTree{T}) where T <: AbstractPoint2D = T

function initnode!(q::QuadTreeNode{T}, r::Number, midx::Number, midy::Number) where T <: AbstractPoint2D
    q.r = r
    q.midx = midx
    q.midy = midy
    q.is_empty = true
    q.is_divided = false
end

function divide!(h::QuadTree{T}, q::QuadTreeNode{T}) where T <: AbstractPoint2D
    # make sure we have enough nodes
    if length(h.nodes) - h.number_of_nodes_used < 4
    	new_size = length(h.nodes)+(length(h.nodes) >>> 1)
    	sizehint!(h.nodes, new_size)
    	for i in 1:(new_size-length(h.nodes))
    		push!(h.nodes, QuadTreeNode(T))
    	end
    end

    # populate new nodes
    q.lxly = h.nodes[h.number_of_nodes_used+1]
    q.lxhy = h.nodes[h.number_of_nodes_used+2]
    q.hxly = h.nodes[h.number_of_nodes_used+3]
    q.hxhy = h.nodes[h.number_of_nodes_used+4]

    # set new nodes properties (dimensions etc.)
    r2 = q.r/2
    initnode!(q.lxly, r2, q.midx-r2, q.midy-r2)
    initnode!(q.lxhy, r2, q.midx-r2, q.midy+r2)
    initnode!(q.hxly, r2, q.midx+r2, q.midy-r2)
    initnode!(q.hxhy, r2, q.midx+r2, q.midy+r2)

    # update tree and parent node
    h.number_of_nodes_used += 4
    q.is_divided = true
    if !q.is_empty
        # move point in parent node to child node
        sq = _getsubnode(q, q.point)
        sq.is_empty = false
        q.is_empty = true
        sq.point, q.point = q.point, sq.point
    end
    q
end

function _getsubnode(q::QuadTreeNode{T}, point::T) where T <: AbstractPoint2D
    x=getx(point)
    y=gety(point)
    if x<q.midx
        y<q.midy && return q.lxly
        return q.lxhy
    end
    y<q.midy && return q.hxly
    return q.hxhy
end

function map(t::QuadTree{T}, cond_data) where T <: AbstractPoint2D
    curr_stack_ix = 1
    t.faststack[1] = t.head
    @inbounds while curr_stack_ix > 0
        q = t.faststack[curr_stack_ix]
        curr_stack_ix -= 1
        if !stop_cond(q, cond_data) && q.is_divided
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.lxly
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.lxly
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.lxhy
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.lxhy
        end
    end
end


text:
module StippleCharts

import Genie
import Stipple

#===#

const COMPONENTS = [:apexchart => :VueApexCharts]

#===#

function deps() :: String
  Genie.Router.route("/js/stipple/apexcharts.min.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "js", "apexcharts.min.js"), String),
      :javascript) |> Genie.Renderer.respond
  end

  Genie.Router.route("/js/stipple/vue-apexcharts.min.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "js", "vue-apexcharts.min.js"), String),
      :javascript) |> Genie.Renderer.respond
  end

  string(
    Genie.Renderer.Html.script(src="/js/stipple/apexcharts.min.js"),
    Genie.Renderer.Html.script(src="/js/stipple/vue-apexcharts.min.js")
  )
end

#===#

include("Charts.jl")

function __init__()
  push!(Stipple.DEPS, deps)
end

end # module


text:
using KhepriFreeCAD
using Test

@testset "KhepriFreeCAD.jl" begin
    # Write your tests here.
end


text:
export RandomWalk1D

"""
    RandomWalk1D(;rewards=-1. => 1.0, N=7, start_pos=(N+1) ÷ 2, actions=[-1,1])

An agent is placed at the `start_pos` and can move left or right (stride is
defined in actions). The game terminates when the agent reaches either end and
receives a reward correspondingly.

Compared to the [`MultiArmBanditsEnv`](@ref):

1. The state space is more complicated (well, not that complicated though).
1. It's a sequential game of multiple action steps.
1. It's a deterministic game instead of stochastic game.
"""
Base.@kwdef mutable struct RandomWalk1D <: AbstractEnv
    rewards::Pair{Float64,Float64} = -1.0 => 1.0
    N::Int = 7
    actions::Vector{Int} = [-1, 1]
    start_pos::Int = (N + 1) ÷ 2
    pos::Int = start_pos
end

RLBase.action_space(env::RandomWalk1D) = Space(OneTo(length(env.actions)))

function (env::RandomWalk1D)(action)
    env.pos = max(min(env.pos + env.actions[action], env.N), 1)
end

RLBase.state(env::RandomWalk1D) = env.pos
RLBase.state_space(env::RandomWalk1D) = Space(OneTo(env.N))
RLBase.is_terminated(env::RandomWalk1D) = env.pos == 1 || env.pos == env.N
RLBase.reset!(env::RandomWalk1D) = env.pos = env.start_pos

function RLBase.reward(env::RandomWalk1D)
    if env.pos == 1
        first(env.rewards)
    elseif env.pos == env.N
        last(env.rewards)
    else
        0.0
    end
end

RLBase.NumAgentStyle(::RandomWalk1D) = SINGLE_AGENT
RLBase.DynamicStyle(::RandomWalk1D) = SEQUENTIAL
RLBase.ActionStyle(::RandomWalk1D) = MINIMAL_ACTION_SET
RLBase.InformationStyle(::RandomWalk1D) = PERFECT_INFORMATION
RLBase.StateStyle(::RandomWalk1D) = Observation{Int}()
RLBase.RewardStyle(::RandomWalk1D) = TERMINAL_REWARD
RLBase.UtilityStyle(::RandomWalk1D) = GENERAL_SUM
RLBase.ChanceStyle(::RandomWalk1D) = DETERMINISTIC


text:
module Day5
export is_reacting
export apply_reaction
export part_1
export part_2

function is_reacting(input::AbstractString)
	for i in 1:length(input) - 1
		a, b = input[i], input[i + 1]
		if a != b
			if islowercase(a) & (a == lowercase(b))
				return true
			elseif isuppercase(a) & (a == uppercase(b))
				return true
			end
		end
	end
	return false
end

function apply_reaction(input::AbstractString; count=1)
	for i in 1:length(input) - 1
		a, b = input[i], input[i + 1]
		if a == b
			continue
		end

		if (islowercase(a) & (a == lowercase(b))) | (isuppercase(a) & (a == uppercase(b)))
			return replace(input, (a * b) => ""; count=count)
		end
	end
	return input
end

function part_1(input::AbstractString)
	last = ""
	reacted = apply_reaction(input)
	while last != reacted
		last = reacted
		reacted = apply_reaction(reacted, count=length(input))
	end
	return length(reacted)
end

function part_2(input::AbstractString)
	polymers = collect(Set(lowercase(input)))

	return minimum(map(p -> begin
		without_p = replace(input, p => "")
		without_p = replace(without_p, uppercase(p) => "")
		return part_1(without_p)
	end, polymers))
end

end


text:
using MPI, MPIArrays

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
n_proc = MPI.Comm_size(comm)

# 1
# default blocksize = (1, 1)
A = CyclicMPIArray(Float64, 7, 5)
rma!(A) do
    forlocalpart!(x -> fill!(x, rank), A)
    sync(A)
    if rank == 0
        println("size $(size(A))")
        show(stdout, "text/plain", A)
        println()
    end
end

# 2
A = CyclicMPIArray(Float64, 9, 7, blocksizes=(2, 2))
rma!(A) do
    forlocalpart!(x -> fill!(x, rank), A)
    sync(A)
    if rank == 0
        println("size $(size(A))")
        show(stdout, "text/plain", A)
        println()
    end
end

# 3
A = CyclicMPIArray(ComplexF64, 7, 5)
rma!(A) do
    forlocalpart!(x -> fill!(x, rank), A)
    sync(A)
    if rank == 0
        println("size $(size(A))")
        show(stdout, "text/plain", A)
        println()
    end
end

# 4
if n_proc % 2 == 0
    A = CyclicMPIArray(Int64, 11, 11, proc_grids=(Int(n_proc/2), 2), blocksizes=(2, 2))
    rma!(A) do
        forlocalpart!(x -> fill!(x, rank), A)
        sync(A)
        if rank == 0
            println("size $(size(A))")
            show(stdout, "text/plain", A)
            println()
        end
    end
end

# 5
if n_proc % 2 == 0
    A = CyclicMPIArray(Int64, 11, 9, proc_grids=(Int(n_proc/2), 2), blocksizes=(3, 3))

    print("rank: $rank, $(size(A.localarray))\n")
    println("$(collect(localindices(A)[1]))")

    rma!(A) do
        forlocalpart!(x -> fill!(x, rank), A)
        sync(A)
        if rank == 0

            println("size $(size(A))")
            show(stdout, "text/plain", A)
            println()
        end
    end
end

MPI.Finalize()


text:
using Base.Test
#using Grisu
cd("C:\\Users\\karbarcca\\.julia\\v0.3\\Grisu\\test")
include("C:/Users/karbarcca/.julia/v0.3/Grisu/src/Grisu.jl")

function trimrep(buffer)
  len = length(bytestring(pointer(buffer)))
  ind = len
  for i = len:-1:1
    buffer[i] != 0x30 && break
    ind -= 1
  end
  buffer[ind+1] = 0
  return bytestring(pointer(buffer))
end

const bufsize = 500
buffer = Array(Uint8,bufsize);
fill!(buffer,0);

include("float.jl")
include("fastdtoa.jl")
include("fastfixedtoa.jl")
include("bignum-dtoa.jl")

# Huge precomputed tests; 100,000 each
if !isdefined(:ShortestTestNumbers)
  include("gay-shortest.jl")
end
if !isdefined(:ShortestSingleTestNumbers)
  include("gay-shortest-single.jl")
end
if !isdefined(:PrecisionTestNumbers)
  include("gay-precision.jl")
end
if !isdefined(:FixedTestNumbers)
  include("gay-fixed.jl")
end

function test(t::PrecomputedShortest,buffer,f)
    status,len,point,buffer = f(t.v,Grisu.SHORTEST,0,buffer)
    if status
      @test t.representation == bytestring(pointer(buffer))
      @test t.decimal_point == point
    end
end
function test(t::PrecomputedShortestSingle,buffer,f)
    status,len,point,buffer = f(t.v,Grisu.SHORTEST,0,buffer)
    if status
      @test t.representation == bytestring(pointer(buffer))
      @test t.decimal_point == point
    end
end
function test(t::PrecomputedFixed,buffer,f)
    n = t.number_digits
    status,len,point,buffer = f(t.v,Grisu.FIXED,n,buffer)
    if status
      @test n >= len - point - 1
      @test t.representation == trimrep(buffer)
      @test t.decimal_point == point
    end
end
function test(t::PP,buffer,f)
    n = t.number_digits
    status,len,point,buffer = f(t.v,Grisu.PRECISION,n,buffer)
    if status
      @test t.number_digits >= len-1
      @test t.representation == trimrep(buffer)
      @test t.decimal_point == point
    end
end

function test(tests,f)
  const bufsize = 500
  buffer = Array(Uint8,bufsize);
  fill!(buffer,0);
  iter = 1
  for t in tests
    test(t,buffer,f)
    fill!(buffer,0);
    iter += 1
  end
  return iter
end

@time test(ShortestTestNumbers,Grisu.fastdtoa)
@time test(ShortestSingleTestNumbers,Grisu.fastdtoa)
@time test(PrecisionTestNumbers,Grisu.fastdtoa)
@time test(FixedTestNumbers,Grisu.fastfixedtoa)

@time test(ShortestTestNumbers,Grisu.bignumdtoa)
@time test(ShortestSingleTestNumbers,Grisu.bignumdtoa)
@time test(PrecisionTestNumbers,Grisu.bignumdtoa)
@time test(FixedTestNumbers,Grisu.bignumdtoa)

include("dtoa.jl")


text:
using ParallelUtils
using Test

@testset "ParallelUtils.jl" begin
    # Write your tests here.
end


text:
@testset "Noise Interpolation Test" begin

using DiffEqNoiseProcess

W = WienerProcess(0.0,0.0,0.0)

dt = 0.1
calculate_step!(W,dt)

for i in 1:10
  accept_step!(W,dt)
end

dt = dt/100
for t in dt:dt:1-dt
  W(t)
end

dt = dt/100
for t in dt:dt:1-dt
  W(t)
end

end


text:
module BasicTraits
export ismutable, isimmutable, isiterable, iscallable,
  isbitstype, isconcretetype

"""
    WhereTraits.BasicTraits.@overwrite_Base

Overloading functions already defined in base, so that they work similar to `Base.eltype`
in that they are best defined on Types, with a fallback from values which just grabs the type of the value
and tries it again.

The following functions are currently over-defined:
- isimmutable
- ismutable
- isbitstype
- isconcretetype

In addition, also `using WhereTraits.BasicTraits` is run.
"""
macro overwrite_Base()
  esc(quote
    using WhereTraits.BasicTraits
    const isimmutable = WhereTraits.BasicTraits.isimmutable
    const ismutable = WhereTraits.BasicTraits.ismutable
    const isbitstype = WhereTraits.BasicTraits.isbitstype
    const isconcretetype = WhereTraits.BasicTraits.isconcretetype
    nothing  # empty return
  end)
end

# mimicking typeclasses from SimpleTraits.BaseTraits
# --------------------------------------------------

"""
    ismutable(v) -> Bool

Return `true` iff value `v` is mutable.  Note that this function works on Types instead of values
(oppositely to `Base.isimmutable`). When applied to a value, a default clause will match on its type instead.

# Examples
```jldoctest
julia> using WhereTraits; WhereTraits.BasicTraits.@overwrite_Base

julia> ismutable(1)
false

julia> ismutable([1,2])
true
```
"""
ismutable(T::DataType) = T.mutable
ismutable(::Type) = false
ismutable(value) = ismutable(typeof(value))

# isimmutable is already part of Base export, hence we better use it,
# however it says "Note that this function works on values"
# i.e. we could overwrite the Base.isimmutable for Type, however this is Type piracy and potentially dangerous
# hence we define our own isimmutable
"""
    isimmutable(v) -> Bool

Return `true` iff value `v` is immutable.  Note that this function works on Types instead of values
(oppositely to `Base.isimmutable`). When applied to a value, a default clause will match on its type instead.

# Examples
```jldoctest
julia> using WhereTraits; WhereTraits.BasicTraits.@overwrite_Base

julia> isimmutable(1)
true

julia> isimmutable([1,2])
false
```
"""
isimmutable(any) = !ismutable(any)

"""
    isiterable(T) -> Bool

Test if type `T` is an iterable collection type or not,
that is whether it has an `iterate` method or not.

When given a value instead of a Type, it fallbacks to use `typeof(value)`.
"""
isiterable(T::Type) = Base.isiterable(T)
isiterable(any) = isiterable(typeof(any))

"""
    callable(T) -> Bool

Checks whether the call-syntax is defined for the given Type.
For convenience `iscallable(value) = iscallable(typeof(value))`

# Examples
```jldoctest
julia> using WhereTraits; WhereTraits.BasicTraits.@overwrite_Base

julia> iscallable(typeof(+))
true

julia> iscallable(+)
true

julia> iscallable(Some)
false

julia> iscallable(typeof(Some))
true
```
"""
@generated function iscallable(::Type{T}) where T
  # without `@generated` it fails to type infer properly
  lim = -1
  world = typemax(UInt)
  methods = Base._methods_by_ftype(Tuple{T, Vararg}, lim, world)
  !isempty(methods)
end
iscallable(value) = iscallable(typeof(value))

"""
    isbitstype(T)

Return `true` if type `T` is a "plain data" type,
meaning it is immutable and contains no references to other values,
only `primitive` types and other `isbitstype` types.
Typical examples are numeric types such as `UInt8`,
`Float64`, and `Complex{Float64}`.
This category of types is significant since they are valid as type parameters,
may not track `isdefined` / `isassigned` status,
and have a defined layout that is compatible with C.

For convenience `isbitstype(value) = isbitstype(typeof(value))`

# Examples
```jldoctest
julia> using WhereTraits; WhereTraits.BasicTraits.@overwrite_Base

julia> isbitstype(Complex{Float64})
true

julia> isbitstype(Complex)
false

julia> isbitstype(1)
true
```
"""
isbitstype(T::Type) = Base.isbitstype(T)
isbitstype(value) = isbitstype(typeof(value))

"""
    isconcretetype(T)

Determine whether type `T` is a concrete type, meaning it could have direct instances
(values `x` such that `typeof(x) === T`).

For convenience `isconcretetype(value) = isconcretetype(typeof(value))`

# Examples
```jldoctest
julia> using WhereTraits; WhereTraits.BasicTraits.@overwrite_Base

julia> isconcretetype(Complex)
false

julia> isconcretetype(Complex{Float32})
true

julia> isconcretetype(Vector{Complex})
true

julia> isconcretetype(Vector{Complex{Float32}})
true

julia> isconcretetype(Union{})
false

julia> isconcretetype(Union{Int,String})
false

julia> isconcretetype([1,2,3])
true

julia> isconcretetype("hi")
true
```
"""
isconcretetype(T::Type) = Base.isconcretetype(T)
isconcretetype(value) = isconcretetype(typeof(value))

end  # module