text:
_hold = false
hold() = (global _hold = !_hold)
hold(h::Bool) = (global _hold = h)

function ghf()
    if !_hold
        global _pwinston = FramedPlot()
    end
    _pwinston
end
ghf(p) = (global _pwinston = p)

savefig(fname::AbstractString, args...; kvs...) = savefig(_pwinston, fname, args...; kvs...)
@deprecate file savefig

for f in (:xlabel,:ylabel,:title)
    @eval $f(s::AbstractString) = (setattr(_pwinston, $f=s); _pwinston)
end
for (f,k) in ((:xlim,:xrange),(:ylim,:yrange))
    @eval $f(a, b) = (setattr(_pwinston, $k=(a,b)); _pwinston)
    @eval $f(a) = (setattr(_pwinston, $k=(a[1],a[2])); _pwinston)
    @eval $f() = $k(limits(_pwinston))
end

const chartokens = @Dict(
    '-' => (:linekind, "solid"),
    ':' => (:linekind, "dotted"),
    ';' => (:linekind, "dotdashed"),
    '+' => (:symbolkind, "plus"),
    'o' => (:symbolkind, "circle"),
    '*' => (:symbolkind, "asterisk"),
    '.' => (:symbolkind, "dot"),
    'x' => (:symbolkind, "cross"),
    's' => (:symbolkind, "square"),
    'd' => (:symbolkind, "diamond"),
    '^' => (:symbolkind, "triangle"),
    'v' => (:symbolkind, "down-triangle"),
    '>' => (:symbolkind, "right-triangle"),
    '<' => (:symbolkind, "left-triangle"),
    'y' => (:color, colorant"yellow"),
    'm' => (:color, colorant"magenta"),
    'c' => (:color, colorant"cyan"),
    'r' => (:color, colorant"red"),
    'g' => (:color, colorant"green"),
    'b' => (:color, colorant"blue"),
    'w' => (:color, colorant"white"),
    'k' => (:color, colorant"black"),
)

function _parse_spec(spec::AbstractString)
    style = Dict()

    try
        style[:color] = parse(Colors.Colorant, spec)
        return style
    end

    for (k,v) in (("--","dashed"), ("-.","dotdashed"))
        splitspec = split(spec, k)
        if length(splitspec) > 1
            style[:linekind] = v
            spec = join(splitspec)
        end
    end

    for char in spec
        if haskey(chartokens, char)
            (k,v) = chartokens[char]
            style[k] = v
        else
            warn("unrecognized style '$char'")
        end
    end

    style
end

function default_color(i::Int)
    cs = [0x000000, 0xED2C30, 0x008C46, 0x1859A9,
          0xF37C21, 0x652B91, 0xA11C20, 0xB33794]
    cs[mod1(i,length(cs))]
end

function _process_keywords(kvs, p, components...)
    for (k,v) in kvs
        if k in [:angle,
                 :color,
                 :face,
                 :halign,
                 :linecolor,
                 :linekind,
                 :linewidth,
                 :size,
                 :symbolkind,
                 :symbolsize,
                 :valign]
            for c in components
                style(c, k, v)
            end
        else
            setattr(p, k, v)
        end
    end
end

@compat typealias PlotArg Union{AbstractString,AbstractVector,AbstractMatrix,Function,Real}

isrowvec(x::AbstractArray) = ndims(x) == 2 && size(x,1) == 1 && size(x,2) > 1

isvector(x::AbstractVector) = true
isvector(x::AbstractMatrix) = size(x,1) == 1

function plot(p::FramedPlot, args::PlotArg...; kvs...)
    args = Any[args...]
    components = Any[]
    color_idx = 0

    default_style = Dict()
    attr = Any[]
    xrange = nothing
    for (k,v) in kvs
        if k in (:linestyle, :linetype)
            default_style[:linekind] = v
        elseif k in (:marker, :symboltype)
            default_style[:symbolkind] = v
        elseif k in (:markersize,)
            default_style[:symbolsize] = v
        elseif k in (:color, :linekind, :linewidth, :symbolkind, :symbolsize)
            default_style[k] = v
        else
            k == :xrange && (xrange = v)
            push!(attr, (k,v))
        end
    end

    # parse the args into tuples of the form (x, y, spec) or (func, lims, spec)
    parsed_args = Any[]

    i = 0
    need_xrange = false
    while length(args) > 0
        local x, y
        a = shift!(args); i += 1
        if isa(a, Function)
            x = a
            if length(args) > 1 && isa(args[1],Real) && isa(args[2],Real)
                y = (shift!(args),shift!(args)); i += 2
            else
                y = ()
                need_xrange = true
            end
        elseif isa(a, AbstractVecOrMat)
            elt = eltype(a)
            if elt <: Complex
                x = real(a)
                y = imag(a)
            elseif length(args) > 0 && isa(args[1], AbstractVecOrMat) &&
               elt <: Real && eltype(args[1]) <: Real
                x = a
                y = shift!(args); i += 1
            elseif elt <: Real
                y = a
                x = 1:(isrowvec(y) ? size(y,2) : size(y,1))
            else
                error("eltype of argument #$i is not Real or Complex")
            end
        else
            error("expected array or function for argument #$i; got $(typeof(a))")
        end
        spec = ""
        if length(args) > 0 && isa(args[1], AbstractString)
            spec = shift!(args); i += 1
        end
        push!(parsed_args, (x,y,spec))
    end

    need_xrange && xrange === nothing && error("need to specify xrange")

    for (a,b,spec) in parsed_args
        local x, y
        if isa(a, Function)
            xlim = b == () ? xrange : b
            x, y = fplot_points(a, xlim[1], xlim[2])
        else
            x, y = a, b
        end

        sopts = copy(default_style)
        spec != "" && merge!(sopts, _parse_spec(spec))

        no_color = !haskey(sopts, :color)
        add_curve = haskey(sopts, :linekind) || !haskey(sopts, :symbolkind)
        add_points = haskey(sopts, :symbolkind)

        isvector(x) && (x = vec(x))
        isvector(y) && (y = vec(y))

        local xys
        if isa(x, AbstractVector) && isa(y, AbstractVector)
            xys = [ (x,y) ]
        elseif isa(x, AbstractVector)
            xys = length(x) == size(y,1) ?
                  [ (x, sub(y,:,j)) for j = 1:size(y,2) ] :
                  [ (x, sub(y,i,:)) for i = 1:size(y,1) ]
        elseif isa(y, AbstractVector)
            xys = size(x,1) == length(y) ?
                  [ (sub(x,:,j), y) for j = 1:size(x,2) ] :
                  [ (sub(x,i,:), y) for i = 1:size(x,1) ]
        else
            @assert size(x) == size(y)
            xys = [ (sub(x,:,j), sub(y,:,j)) for j = 1:size(y,2) ]
        end

        for (x,y) in xys
            if no_color
                color_idx += 1
                sopts[:color] = default_color(color_idx)
            end
            if add_curve
                push!(components, Curve(x, y, sopts))
            end
            if add_points
                push!(components, Points(x, y, sopts))
            end
        end
    end

    for (k,v) in attr
        setattr(p, k, v)
    end

    for c in components
        add(p, c)
    end

    global _pwinston = p
    p
end
plot(args::PlotArg...; kvs...) = plot(ghf(), args...; kvs...)

# shortcut for overplotting
oplot(args::PlotArg...; kvs...) = plot(_pwinston, args...; kvs...)

# shortcuts for creating log plots
semilogx(args::PlotArg...; kvs...) = plot(args...; xlog=true, kvs...)
semilogy(args::PlotArg...; kvs...) = plot(args...; ylog=true, kvs...)
loglog(args::PlotArg...; kvs...) = plot(args...; xlog=true, ylog=true, kvs...)

typealias Interval @compat(Tuple{Real,Real})

function data2rgb{T<:Real}(data::AbstractArray{T}, limits::Interval, colormap::Array{UInt32,1})
    img = similar(data, UInt32)
    ncolors = length(colormap)
    limlower = limits[1]
    limscale = ncolors/(limits[2]-limits[1])
    for i = 1:length(data)
        datai = data[i]
        if isfinite(datai)
            idxr = limscale*(datai - limlower)
            idx = trunc(Int, idxr)
            idx += idxr > convert(T, idx)
            idx = clamp(idx, 1, ncolors)
            img[i] = colormap[idx]
        else
            img[i] = 0x00000000
        end
    end
    img
end

# from http://www.metastine.com/?p=7
function jetrgb(x)
    fourValue = 4x
    r = min(fourValue - 1.5, -fourValue + 4.5)
    g = min(fourValue - 0.5, -fourValue + 3.5)
    b = min(fourValue + 0.5, -fourValue + 2.5)
    RGB(clamp(r,0.,1.), clamp(g,0.,1.), clamp(b,0.,1.))
end

colormap() = (global _current_colormap; _current_colormap)
colormap(c::Array{UInt32,1}) = (global _current_colormap = c; nothing)
colormap{C<:Color}(cs::Array{C,1}) =
    colormap(UInt32[convert(RGB24,c) for c in cs])
function colormap(name::AbstractString, n::Int=256)
    if name == "jet"
        colormap([jetrgb(x) for x in linspace(0.,1.,n)])
    else
        colormap(Colors.colormap(name, n))
    end
end
colormap("jet")

function imagesc{T<:Real}(xrange::Interval, yrange::Interval, data::AbstractArray{T,2}, clims::Interval)
    p = ghf()
    if !_hold
        setattr(p, :xrange, xrange)
        setattr(p, :yrange, reverse(yrange))
    end
    img = data2rgb(data, clims, _current_colormap)
    xrange[1] > xrange[2] && (img = flipdim(img,2))
    yrange[1] < yrange[2] && (img = flipdim(img,1))
    add(p, Image(xrange, reverse(yrange), img))
    ghf(p)
end

imagesc(xrange, yrange, data) = imagesc(xrange, yrange, data, (minimum(data),maximum(data)+1))
imagesc(data) = ((h, w) = size(data); imagesc((0,w), (0,h), data))
imagesc{T}(data::AbstractArray{T,2}, clims::Interval) = ((h, w) = size(data); imagesc((0,w), (0,h), data, clims))

function spy(S::SparseMatrixCSC, nrS::Integer, ncS::Integer)
    m, n = size(S)
    colptr = S.colptr
    rowval = S.rowval
    nzval  = S.nzval

    if nrS > m; nrS = m; end
    if ncS > n; ncS = n; end

    target = zeros(nrS, ncS)
    x = nrS / m
    y = ncS / n

    for col = 1:n
        for k = colptr[col]:colptr[col+1]-1
            row = rowval[k]
            target[ceil(row * x), ceil(col * y)] += 1
        end
    end

    imagesc((1,m), (1,n), target)
end

scatter(x::AbstractVecOrMat, y::AbstractVecOrMat, spec::String="o"; kvs...) = scatter(x, y, 1., spec; kvs...)
scatter{C<:Complex}(z::AbstractVecOrMat{C}, spec::String="o"; kvs...) = scatter(real(z), imag(z), 1., spec; kvs...)
function scatter(x::AbstractVecOrMat, y::AbstractVecOrMat,
                 s::Real, spec::String="o"; kvs...)
    sopts = _parse_spec(spec)
    p = ghf()
    c = Points(x, y, sopts, symbolsize=s)
    add(p, c)
    for (k,v) in kvs
        if k in [:linekind,:symbolkind,:color,:linecolor,:linewidth,:symbolsize]
            style(c, k, v)
        else
            setattr(p, k, v)
        end
    end
    ghf(p)
end
function scatter(x::AbstractVecOrMat, y::AbstractVecOrMat,
                 s::AbstractVecOrMat, spec::String="o"; kvs...)
    c = convert(RGB24, color(get(_parse_spec(spec), :color, RGB(0,0,0))))
    scatter(x, y, s, fill(c,size(x)...), spec; kvs...)
end
function scatter(x::AbstractVecOrMat, y::AbstractVecOrMat,
                 s::(@compat Union{Real,AbstractVecOrMat}), c::AbstractVecOrMat,
                 spec::String="o"; kvs...)
    if typeof(s) <: Real
        s = fill(s, size(x)...)
    end
    if eltype(c) <: Real
        c = data2rgb(c, extrema(c), _current_colormap)
    elseif !(eltype(c) <: Color)
        error("bad color array")
    end
    sopts = _parse_spec(spec)
    p = ghf()
    c = ColoredPoints(x, y, s, c, sopts)
    add(p, c)
    for (k,v) in kvs
        if k in [:linekind,:symbolkind,:color,:linecolor,:linewidth,:symbolsize]
            style(c, k, v)
        else
            setattr(p, k, v)
        end
    end
    ghf(p)
end

## stem ##

stem(y::AbstractVecOrMat, spec::String="o"; kvs...) = stem(1:length(y), y, spec; kvs...)
function stem(x::AbstractVecOrMat, y::AbstractVecOrMat, spec::String="o"; kvs...)
    p = ghf()
    sopts = _parse_spec(spec)
    s = Stems(x, y, sopts)
    haskey(sopts,:symbolkind) || (sopts[:symbolkind] = "circle")
    o = Points(x, y, sopts)
    _process_keywords(kvs, p, s, o)
    add(p, s, o)
    ghf(p)
end

function text(x::Real, y::Real, s::AbstractString; kvs...)
    p = _pwinston
    c = DataLabel(x, y, s, halign="left")
    _process_keywords(kvs, p, c)
    add(p, c)
end

spy(S::SparseMatrixCSC) = spy(S, 100, 100)
spy(A::AbstractMatrix, nrS, ncS) = spy(sparse(A), nrS, ncS)
spy(A::AbstractMatrix) = spy(sparse(A))

function plothist(p::FramedPlot, h::@compat(Tuple{Range,Vector}); kvs...)
    c = Histogram(h...)
    add(p, c)

    for (k,v) in kvs
        if k in [:color,:linecolor,:linekind,:linetype,:linewidth]
            style(c, k, v)
        else
            setattr(p, k, v)
        end
    end

    global _pwinston = p
    p
end
plothist(p::FramedPlot, args...; kvs...) = plothist(p::FramedPlot, hist(args...); kvs...)
plothist(args...; kvs...) = plothist(ghf(), args...; kvs...)

# 3x3 gaussian
#_default_kernel2d=[.05 .1 .05; .1 .4 .1; .05 .1 .05]

# 5x5 gaussian
_default_kernel2d=(1.0/273.)*[1.0 4.0 7.0 4.0 1.0;
                             4.0 16. 26. 16. 4.0;
                             7.0 26. 41. 26. 7.0;
                             1.0 4.0 7.0 4.0 1.0;
                             4.0 16. 26. 16. 4.0]

#hist2d
function plothist2d(p::FramedPlot, h::@compat(Tuple{(@compat Union{Range,Vector}),(@compat Union{Range,Vector}),Array{Int,2}}); colormap=_current_colormap, smooth=0, kernel=_default_kernel2d, kvs...)
    xr, yr, hdata = h

    for i in 1:smooth
        hdata = conv2(hdata*1.0, kernel)
    end

    clims = (minimum(hdata), maximum(hdata)+1)
    img = data2rgb(hdata, clims, colormap)'
    add(p, Image((xr[1], xr[end]), (yr[1], yr[end]), img;))

    #XXX: check if there is any Image-related named arguments
    setattr(p; kvs...)

    global _pwinston = p
    p
end
plothist2d(p::FramedPlot, args...; kvs...) = plothist2d(p::FramedPlot, hist2d(args...); kvs...)
plothist2d(args...; kvs...) = plothist2d(ghf(), args...; kvs...)

#errorbar
errorbar(args...; kvs...) = errorbar(ghf(), args...; kvs...)
function errorbar(p::FramedPlot, x::AbstractVector, y::AbstractVector; xerr=nothing, yerr=nothing, kvs...)

    xn=length(x)
    yn=length(y)

    if xerr != nothing
        xen = length(xerr)
        if xen == xn
            cx = SymmetricErrorBarsX(x, y, xerr)
        elseif xen == 2xn
            cx = ErrorBarsX(y, x.-xerr[1:xn], x.+xerr[xn+1:xen])
        else
            warn("Dimensions of x and xerr do not match!")
        end
        style(cx; kvs...)
        add(p,cx)
    end

    if yerr != nothing
        yen=length(yerr)
        if yen == yn
            cy = SymmetricErrorBarsY(x, y, yerr)
        elseif yen == 2yn
            cy = ErrorBarsY(x, y.-yerr[1:yn], y.+yerr[yn+1:yen])
        else
            warn("Dimensions of y and yerr do not match!")
        end
        style(cy; kvs...)
        add(p,cy)
    end

    global _pwinston = p
    p
end

function fplot_points(f::Function, xmin::Real, xmax::Real;
        max_recursion::Int=6, min_points::Int=10, tol::Float64=0.01)
    @assert xmin < xmax
    @assert min_points > 1
    @assert max_recursion >= 0

    xs = Float64[]
    ys = Float64[]
    cs = Array(Float64, max_recursion)
    fcs = Array(Float64, max_recursion)
    ls = Array(Int, max_recursion)

    local c::Float64
    local fc::Float64

    function good(a, b, c, fa, fb, fc, tol)
        u = b - a
        fu = fb - fa
        v = c - b
        fv = fc - fb
        n = u*v + fu*fv
        d2 = (u*u + fu*fu)*(v*v + fv*fv)
        n*n > d2*(1. - tol)
    end

    p = linspace(xmin, xmax, min_points)
    q = [f(x) for x in p]

    if max_recursion == 0
        return p, q
    end

    for i = 1:length(p)-1
        a::Float64  = p[i]
        fa::Float64 = q[i]
        c  = p[i+1]
        fc = q[i+1]
        level::Int = 0
        n::Int = 0
        while true
            b::Float64 = 0.5(a + c)
            fb::Float64 = f(b)
            g1::Bool = good(a,b,c,fa,fb,fc,tol)
            g2::Bool = length(xs) > 0 ? good(xs[end],a,b,ys[end],fa,fb,tol) : true
            if (g1 && g2) || level == max_recursion
                push!(xs, a, b)
                push!(ys, fa, fb)
                a = c
                fa = fc
                n == 0 && break
                c = cs[n]
                fc = fcs[n]
                level = ls[n]
                n -= 1
            else
                level += 1
                n += 1
                ls[n] = level
                cs[n] = c
                fcs[n] = fc
                c = b
                fc = fb
            end
        end
    end

    push!(xs, c)
    push!(ys, fc)

    xs, ys
end

function fplot(f::Function, limits, args...; kvs...)
    pargs = []
    fopts = Dict()
    for arg in args
        if typeof(arg) <: AbstractString
            pargs = [arg]
        elseif typeof(arg) <: Integer
            fopts[:min_points] = arg
        elseif typeof(arg) <: @compat AbstractFloat
            fopts[:tol] = arg
        else
            error("unrecognized argument ", arg)
        end
    end
    xmin = limits[1]
    xmax = limits[2]
    x,y = fplot_points(f, xmin, xmax; fopts...)
    plot(x, y, pargs...; kvs...)
end

# bar, barh
ax = @compat Dict{Any,Any}(:bar => :x, :barh => :y)
ax1 = @compat Dict{Any,Any}(:bar => :x1, :barh => :y1)
vert = @compat Dict{Any,Any}(:bar => true, :barh => false)
for fn in (:bar, :barh)
    eval(quote
          function $fn(p::FramedPlot, b::FramedBar, args...; kvs...)
              setattr(b, vertical=$(vert[fn]))
              setattr(p.$(ax[fn]), draw_subticks=false)
              setattr(p.$(ax[fn]), ticks=collect(1.:length(b.h)))
              setattr(p.$(ax1[fn]), ticklabels=b.g)
              add(p, b)
              global _pwinston = p
              p
          end
          function $fn(p::FramedPlot, g::AbstractVector, h::AbstractVector, args...; kvs...)
              b = FramedBar(g, h[:,end], args...; kvs...)
              $fn(p, b, args...; kvs...)
          end
          function $fn(p::FramedPlot, g::AbstractVector, h::AbstractMatrix, args...; kvs...)
              nc = size(h,2)
              barwidth = config_value("FramedBar", "barwidth")/nc
              offsets = barwidth * (nc - 1) * linspace(-.5, .5, nc)
              for c = 1:nc-1
                  b = FramedBar(g, h[:,c], args...; kvs...)
                  setattr(b, offset=offsets[c])
                  setattr(b, barwidth=barwidth)
                  setattr(b, vertical=$(vert[fn]))
                  style(b, fillcolor=default_color(c))
                  style(b, draw_baseline=false)
                  add(p, b)
              end
              b = FramedBar(g, h[:,nc], args...; kvs...)
              setattr(b, offset=offsets[nc])
              setattr(b, barwidth=barwidth)
              style(b, fillcolor=default_color(nc))
              $fn(p, b, args...; kvs...)
          end
          $fn(p::FramedPlot, h::AbstractVecOrMat, args...; kvs...) =
              $fn(p, [1:size(h,1)], h, args...; kvs...)
          $fn(args...; kvs...) = $fn(ghf(), args...; kvs...)
    end )
end

grid(p::FramedPlot, tf::Bool) = (setattr(p.frame, draw_grid=tf); p)
grid(p::FramedPlot) = grid(p, !any(map(x->getattr(x, "draw_grid"), p.frame.objs)))
grid(args...) = grid(_pwinston, args...)

function legend(p::FramedPlot, lab::AbstractVector, args...; kvs...)
    if length(args) > 0 && length(args[1]) == 2 && eltype(args[1]) <: Real
        position = args[1]
        args = args[2:end]
    elseif length(args) > 1 && eltype(args[1:2]) <: Real
        position = args[1:2]
        args = args[3:end]
    else
        position = [0.1, 0.9]
    end
    # TODO: define other legend positions
    plotcomp = getcomponents(p)
    nitems = min(length(lab), length(plotcomp))
    for c in 1:nitems
        setattr(plotcomp[c], label=lab[c])
    end
    add(p, Legend(position..., plotcomp[1:nitems], args...; kvs...))
end
legend(lab::AbstractVector, args...; kvs...) = legend(_pwinston, lab, args...; kvs...)

function timeplot(p::FramedPlot, x::Vector{DateTime}, y::AbstractArray, args...; kvs...)
    limits = datetime2unix([minimum(x), maximum(x)])

    ticks = collect(0.0:0.2:1.0)
    ticklabels = x[round(Int64, ticks * (length(x) - 1) + 1)]
    normalized_x = (datetime2unix(x) - limits[1]) / (limits[2] - limits[1])

    span = @compat Int(x[end] - x[1]) / 1000
    kvs = Dict(kvs)

    if :format in keys(kvs)
        format = kvs[:format]
        delete!(kvs, :format)
    else
        if span > 365 * 24 * 60 * 60 # 1 year
            format = "%Y-%m"
        elseif 365 * 24 * 60 * 60 > span > 30 * 24 * 60 * 60 # 1 month
            format = "%Y-%m-%d"
        elseif 30 * 24 * 60 * 60 > span > 24 * 60 * 60 # 1 day
            format = "%Y-%m-%d\n%H:%M"
        elseif 24 * 60 * 60 > span > 60 * 60 # 1 hour
            format = "%H:%M"
        elseif 60 * 60 > 60 # 1 minute
            format = "%H:%M:%S"
        else
            format = "%H:%M:%S"
        end
    end

    ticklabels = map(d -> strftime(format, datetime2unix(d)), ticklabels)

    setattr(p.x1, :ticklabels, ticklabels)
    setattr(p.x1, :ticks, ticks)
    setattr(p.x1, :ticklabels_style, @compat Dict(:fontsize=>1.5))

    plot(p, normalized_x, y, args...; kvs...)
end

timeplot(x::Vector{DateTime}, y::AbstractArray, args...; kvs...) = timeplot(ghf(), x, y, args...; kvs...)
timeplot(x::Vector{Date}, y::AbstractArray, arg...; kvs...) = timeplot(ghf(), DateTime(x), y, arg...; kvs...)
timeplot(p::FramedPlot, x::Vector{Date}, y::AbstractArray, arg...; kvs...) =
    timeplot(p, DateTime(x), y, arg...; kvs...)


text:
using Molly
using Test
using SingleCrystal

function forces_are_zero(forces; dims=3)
    return all([isapprox(f, zero(rand(3)), atol=1e-4) for f in forces])
end

function glue_remains_reasonable(glue_init,glue_end)
    return all(isapprox.(glue_init,glue_end,rtol=.1))
end

function groundstate_energy_as_expected(s, ref)
    return isapprox(s.loggers["energy"].energies[1]/length(s.coords), -ref["u"], atol=1e-2)
end

function vacancy_formation_energy_as_expected(s, s_vac, ref)
    u_gs = s.loggers["energy"].energies[1]
    n_atoms = length(s.coords)
    u_vac = s_vac.loggers["energy"].energies[1]
    n_atoms_vac = length(s_vac.coords)
    u_vac_form = u_vac - n_atoms_vac/n_atoms * u_gs
    return isapprox(u_vac_form, ref["u_vac"], atol=7e-2)
end

make_atom(name,mass) = Atom(name=name,mass=mass)

function run_bcc(element::String, fs_inter, a::Real; T::Real=.01,
        nx::Integer=3, ny::Integer=3, nz::Integer=3, vac::Bool=false)

    masses = Dict("V" => 50.9415, "Nb" => 92.9064, "Ta" => 180.9479,
        "Cr" => 51.996, "Mo" => 95.94, "W" => 183.85,
        "Fe" => 55.847)

    unitcell = SingleCrystal.make_bcc_unitcell(element, a, make_atom)
    supercell = SingleCrystal.make_supercell(unitcell, nx=nx, ny=ny, nz=nz)

    if vac
        # introducing a vacancy
        supercell = SingleCrystal.add_vacancies(supercell, ixs=[1])
    end
    n_atoms = length(supercell.positions)

    velocities = [velocity(atom.mass, T, dims=3) for atom in supercell.atoms]
    nb_matrix = trues(n_atoms,n_atoms)
    dist_cutoff = 2. * a

    loggers = Dict(
        "glue" => GlueDensityLogger(1),
        "forces" => ForceLogger(5),
        "energy" => EnergyLogger(1)
    )

    nf = DistanceNeighborFinder(nb_matrix, 1, dist_cutoff)

    s = Simulation(
        simulator=VelocityVerlet(),
        atoms=supercell.atoms,
        general_inters=(fs_inter,),
        coords=[SVector{3}(v) for v in supercell.positions],
        velocities=velocities,
        temperature=T,
        box_size=supercell.edge_lengths[1],
        timestep=.002,
        n_steps=100,
        neighbor_finder=nf,
        loggers=loggers,
    )
    simulate!(s)
    return s
end

@testset "Glue: Finnis-Sinclair" begin
    T = .01 # K
    fs_inter, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true)

    T *= fs_inter.kb
    for element in elements
        @testset "$element" begin
            a = bcc_lattice_constants[element]

            s = run_bcc(element, fs_inter, a, T=T, vac=false)
            @testset "forces zero" begin
                @test forces_are_zero(s.loggers["forces"].forces[1])
            end
            @testset "reasonable glue" begin
                @test glue_remains_reasonable(s.loggers["glue"].glue_densities[1],
                    s.loggers["glue"].glue_densities[end])
         end
            @testset "checking groundstate" begin
                @test groundstate_energy_as_expected(s, reference_energies[element])
            end

            s_vac = run_bcc(element, fs_inter, a, T=T, vac=true)
            @testset "vacancy forces nonzero" begin
                @test !forces_are_zero(s_vac.loggers["forces"].forces[1])
            end
            @testset "vacancy formation energy" begin
                @test vacancy_formation_energy_as_expected(s, s_vac, reference_energies[element])
            end
        end
    end

end

ext: jl


text:
mutable struct Circuit
    karva::String
    parameters::Array{Any,1}
    fitness::Union{Nothing,AbstractFloat}
end

struct TreeNode{I<:Integer,P<:Real}
    Index::I
    ParentIndex::Union{Nothing,I}
    Child1Index::Union{Nothing,I}
    Child2Index::Union{Nothing,I}
    Type::Char
    Parameter::Union{P,Array{P,1}}
end

isless(a::TreeNode, b::TreeNode) = isless(a.Index, b.Index)
isless(a::Circuit, b::Circuit) = isless(a.fitness, b.fitness)

function number_circuit(circuit)
    l_positions = findall(isterminal ,circuit)
    numbers = [string(i) for i in 1:length(l_positions)]
    offsets = cumsum([length(n) for n in numbers])
    l_positions[2:end] .+= offsets[1:end-1]
    for (p,n) in zip(l_positions,numbers)
        circuit = circuit[1:p]*n*circuit[p+1:end]
    end
    return circuit
end

isoperation(chr) = occursin(chr,"-+")
isterminal(chr) = occursin(chr,"RCLP")

function generatekarva(head,terminals="RCLP")
    operators =  "+-+-"
    all_elements = operators*terminals
    karva = rand(operators)*String(rand(all_elements,head-1))*String(rand(terminals,head+1))
    return karva
end

function karva_parameters(karva)
    parameters = Array{Any}(undef, length(karva))
    ranges = Dict('R'=>4000,'C'=> 0.0001,'L'=> 1,'+'=> 0,'-'=> 0,'P'=>[4000,1], 'W'=>4000)
    for (e,i) in enumerate(karva)
            if e == 'P'
                Element_values[n] = [ranges[e][1]*rand(),ranges[e][2]*rand()]
            else
            parameters[e] = rand()*ranges[i]
        end
    end
    return parameters
end

function karva_to_tree(karva,parameters)
    tree = []
    Root_node = TreeNode(1,nothing,2,3,karva[1],parameters[1])
    push!(tree,Root_node)
    child1_karva_index = 2 ; child2_karva_index = 3
    ChildNode = 4
    free_spaces = 2
    for i in 1:length(karva)
        free_spaces -= 1
        K = karva[i]
        if isoperation(K)
            free_spaces += 2
            if isoperation(karva[child1_karva_index])
                push!(tree,TreeNode(child1_karva_index,i,ChildNode,ChildNode+1,karva[child1_karva_index],parameters[child1_karva_index]))
                child1_karva_index += 2
                ChildNode += 2
            else
                push!(tree,TreeNode(child1_karva_index,i,nothing,nothing,karva[child1_karva_index],parameters[child1_karva_index]))
                child1_karva_index += 2
            end

            if isoperation(karva[child2_karva_index])
                push!(tree,TreeNode(child2_karva_index,i,ChildNode,ChildNode+1,karva[child2_karva_index],parameters[child2_karva_index]))
                child2_karva_index += 2
                ChildNode += 2
            else
                push!(tree,TreeNode(child2_karva_index,i,nothing,nothing,karva[child2_karva_index],parameters[child2_karva_index]))
                child2_karva_index += 2
            end
        end
        if free_spaces == 1
            break
        end
    end
    return tree
end

function karva_to_tree(karva)
    parameters = karva_parameters(karva)
    tree = []
    Root_node = TreeNode(1,nothing,2,3,karva[1],parameters[1])
    push!(tree,Root_node)
    child1_karva_index = 2 ; child2_karva_index = 3
    ChildNode = 4
    free_spaces = 2
    for i in 1:length(karva)
        free_spaces -= 1
        K = karva[i]
        if isoperation(K)
            free_spaces += 2
            if isoperation(karva[child1_karva_index])
                push!(tree,TreeNode(child1_karva_index,i,ChildNode,ChildNode+1,karva[child1_karva_index],parameters[child1_karva_index]))
                child1_karva_index += 2
                ChildNode += 2
            else
                push!(tree,TreeNode(child1_karva_index,i,nothing,nothing,karva[child1_karva_index],parameters[child1_karva_index]))
                child1_karva_index += 2
            end

            if isoperation(karva[child2_karva_index])
                push!(tree,TreeNode(child2_karva_index,i,ChildNode,ChildNode+1,karva[child2_karva_index],parameters[child2_karva_index]))
                child2_karva_index += 2
                ChildNode += 2
            else
                push!(tree,TreeNode(child2_karva_index,i,nothing,nothing,karva[child2_karva_index],parameters[child2_karva_index]))
                child2_karva_index += 2
            end
        end
        if free_spaces == 1
            break
        end
    end
    return tree
end

get_circuit_tree(circuit::Circuit) = karva_to_tree(circuit.karva,circuit.parameters)

listcomponents(tree) = tree[findall(x->isterminal(x.Type),tree)]
listoperations(tree) = tree[findall(x->isoperation(x.Type),tree)]

function removecomponent(tree,node_to_remove)
    Parent = tree[node_to_remove.ParentIndex]
    new_tree_legnth = (length(tree)-2)
    sibling = Parent.Child1Index == node_to_remove.Index ? tree[Parent.Child2Index] : tree[Parent.Child1Index]
    all_indices = collect(1:new_tree_legnth)
    old_to_new_location_dict = Dict(Parent.Index=>Parent.ParentIndex,sibling.Index => Parent.Index)
    children_dict = Dict()
    new_tree = Array{TreeNode}(undef, new_tree_legnth)
    new_tree[Parent.Index] = tree[sibling.Index]
    edited_old_tree = reverse(tree[findall(x-> !(x.Index in [Parent.Index,sibling.Index,node_to_remove.Index]),tree)])
    splice!(all_indices, Parent.Index)
    rev_inds = reverse(all_indices)
    for (i,j) in zip(rev_inds,edited_old_tree)
        new_tree[i] = j
        old_to_new_location_dict[j.Index] = i
    end
    correct_final_tree = Array{TreeNode}(undef, new_tree_legnth)
    for (e,node) in enumerate(new_tree)
        if isterminal(node.Type)
            correct_final_tree[e] = TreeNode(e,old_to_new_location_dict[node.ParentIndex],nothing,nothing,node.Type,node.Parameter)
            haskey(children_dict,old_to_new_location_dict[new_tree[e].ParentIndex]) ? push!(children_dict[old_to_new_location_dict[node.ParentIndex]],e) : children_dict[old_to_new_location_dict[node.ParentIndex]] = [e]
        else
            if node.Index != 1
                haskey(children_dict,old_to_new_location_dict[node.ParentIndex]) ? push!(children_dict[old_to_new_location_dict[node.ParentIndex]],e) : children_dict[old_to_new_location_dict[node.ParentIndex]] = [e]
            end
        end
    end
    for e in 1:new_tree_legnth
        if !isterminal(new_tree[e].Type)
            if e == 1
                correct_final_tree[e] = TreeNode(e,nothing,2,3,new_tree[e].Type,new_tree[e].Parameter)
            else
                correct_final_tree[e] = TreeNode(e,old_to_new_location_dict[new_tree[e].ParentIndex],children_dict[e][1],children_dict[e][2],new_tree[e].Type,new_tree[e].Parameter)
            end
        end
    end
    return correct_final_tree
end

function descendents(Tree,Node)
    if isterminal(Node.Type)
        return [Node]
    else
        descendent1 = Tree[Node.Child1Index]
        descendent2 = Tree[Node.Child2Index]
        if all(isterminal.([descendent1.Type,descendent2.Type]))
            return [Node,descendent1,descendent2]
        else
            return sort(vcat(Node,descendents(Tree,descendent1),descendents(Tree,descendent2)))
        end
    end
end

function readablecircuit(circuit::Circuit)
    tree = karva_to_tree(circuit.karva,circuit.parameters)
    return tree_to_circuit(tree)[1]
end

function tree_to_karva(tree,head=8,terminals = "RCPL")
    coding_karva = join([node.Type for node in tree])
    remaining_length = 2*head+1 - length(coding_karva)
    karva = coding_karva*join(rand(terminals,remaining_length))
    return karva
end

function tree_to_circuit(tree)
    nodecount = length(tree)
    essential_info = [[node.ParentIndex,node.Type,[node.Parameter]] for node in tree]
    for i in 1:((length(essential_info)-1)/2)
        parent1,type1,params1 = pop!(essential_info)
        parent2,type2,params2 = pop!(essential_info)
        essential_info[parent1][2] = essential_info[parent1][2] == '+' ? type1*'-'*type2 : '['*type1*','*type2*']'
        essential_info[parent1][3] = vcat(params1,params2)
    end
return  number_circuit(essential_info[1][2]),  essential_info[1][3]
end

denumber_circuit(Circuit) = replace(Circuit, r"[0-9]" => "")

head_length(karva) = findlast(x->isoperation(x),karva)

function get_karva_elements_and_parameters(circuit,circuit_parameters)
    circuit_bare = foldl(replace,["["=>"","]"=>"","-"=>"+",","=>"-"],init = denumber_circuit(circuit))
    karva_params = Array{Any}(undef,length(circuit_bare))
    counter = 1
    for (e,i) in enumerate(circuit_bare)
        if isterminal(i)
            karva_params[e] = circuit_parameters[counter]
                counter +=1
        else
            karva_params[e] = 0
        end
    end
    operation_indexes = findall(karva_params .== 0)
    terminal_indexes = findall(karva_params .!= 0)
    return circuit_bare, karva_params, operation_indexes,terminal_indexes
end

function exhaustive_search(circuit,circuit_parameters)
    # Get simulated measurement from the circuit
    circuit_bare, karva_params, operation_indexes,terminal_indexes = get_karva_elements_and_parameters(circuit,circuit_parameters)
    simulated_measurement = get_target_impedance(circuit,circuit_parameters)

    # Get all possible array permutations (the next few parts can be swapped for a heuristic if there are too many permutations).
    index_permutations = collect(permutations(1:length(circuit_bare)))
    # Remove permutations that lead to synctactically incorrect circuits.
    permitted_index_permutations = index_permutations[findall(x-> (x[1] in operation_indexes)&&(x[end] in terminal_indexes)&&(x[end-1] in terminal_indexes),index_permutations)]
    # Generate the arrays of potential solution.
    potential_karvas = [circuit_bare[inds] for inds in permitted_index_permutations]
    potential_paramsets = [karva_params[inds] for inds in permitted_index_permutations]
    potential_trees = [karva_to_tree(k,p) for (k,p) in zip(potential_karvas,potential_paramsets)]
    potential_circuits = [tree_to_circuit(t)[1] for t in potential_trees]

    # Limit the potential solutions to those with a correct circuit configuration.
    correct_circuits = findall(potential_circuits .== circuit)
    potential_karvas = potential_karvas[correct_circuits]
    potential_paramsets = potential_paramsets[correct_circuits]
    potential_trees = potential_trees[correct_circuits]

    # Limit the potential solutions to those with a correct function evaluation.
    potential_functions = [tree_to_function(T) for T in potential_trees];
    function_evaluations = [f(params,1000) for (f,params) in zip(potential_functions,potential_paramsets)]
    correct_evaluations = findall(function_evaluations .== simulated_measurement)
    potential_karvas = potential_karvas[correct_evaluations]
    potential_paramsets = potential_paramsets[correct_evaluations]
    potential_trees = potential_trees[correct_evaluations]

    # Obtain solution(s) with the sortest head. Maybe randomly select one of them.
    potential_heads = head_length.(potential_karvas)
    potential_head_minimum = min(potential_heads...)
    minimal_heads = findall(potential_heads .== potential_head_minimum)
    karva_solutions = potential_karvas[minimal_heads]
    parameter_solutions = potential_paramsets[minimal_heads]

    return karva_solutions[1] , parameter_solutions[1]
end

function genetic_algorithm(circuit,circuit_parameters,pop_size = 50) #Maybe try to reimplement this with the Evolutionary package.
    # Preliminaries
    circuit_bare, karva_params, operation_indexes,terminal_indexes = get_karva_elements_and_parameters(circuit,circuit_parameters)
    target_impedance = get_target_impedance(circuit,circuit_parameters)
    # Generate initial population.
    Candidate_index_sets = generate_candidate_index_sets(pop_size,length(karva_params),operation_indexes,terminal_indexes)
    # Evaluate initial population's fitness.
    fitnesses = [indexset_fitness(Candidate_index_sets[i],circuit_bare,karva_params,target_impedance) for i in 1:pop_size]
    # Get the sorted fitnesses and candidate indexsets sorted by fitness.
    order = sortperm(fitnesses)
    Candidate_index_sets = Candidate_index_sets[order]
    fitnesses = fitnesses[order]
    best_possible_fitness = count(z->z==0,karva_params)
    if fitnesses[1] == best_possible_fitness
        return circuit_bare[Candidate_index_sets[1]], karva_params[Candidate_index_sets[1]]
    end
    max_fitness = fitnesses[1]
    unchanged_generation_counter = 0
    while fitnesses[1] != best_possible_fitness && !(unchanged_generation_counter >= 20 && max_fitness < 10)
        Candidate_index_sets,fitnesses = next_generation(Candidate_index_sets,fitnesses,circuit_bare,karva_params,operation_indexes,terminal_indexes,target_impedance,0.05)
        println(max_fitness)
        if fitnesses[1] >= max_fitness
            unchanged_generation_counter += 1
        else
            max_fitness = fitnesses[1]
            unchanged_generation_counter = 0
        end
        if unchanged_generation_counter >= 30 #restart if you get stuck in a local minimum. as long as fitness>10, the algorithm doesnt really have any direction, maybe an additional distance measure should be considered.
            fitnesses = [indexset_fitness(Candidate_index_sets[i],circuit_bare,karva_params,target_impedance) for i in 1:pop_size]
            order = sortperm(fitnesses)
            Candidate_index_sets = Candidate_index_sets[order]
            fitnesses = fitnesses[order]
            best_possible_fitness = count(z->z==0,karva_params)
        end
    end
    return circuit_bare[Candidate_index_sets[1]], karva_params[Candidate_index_sets[1]]
end

function generate_indexset(size,operation_indexes,terminal_indexes)
    initial_operation_constraint = false
    terminal_tail_constraint = false
    indexset = randperm(size)
    while !initial_operation_constraint || !terminal_tail_constraint
        initial_operation_constraint = false
        terminal_tail_constraint = false
        if indexset[1] in operation_indexes
            initial_operation_constraint = true
            if indexset[end-1] in terminal_indexes && indexset[end] in terminal_indexes
                terminal_tail_constraint = true
            else indexset = randperm(size)
            end
        else indexset = randperm(size)
        end

    end
    return indexset
end

function generate_candidate_index_sets(population_size,size,operation_indexes,terminal_indexes)
    return [generate_indexset(size,operation_indexes,terminal_indexes) for i in 1:population_size]
end

function indexset_fitness(indexset,circuit_bare,karva_params,target_impedance)
    candidate_solution = circuit_bare[indexset]
    candidate_parameters = karva_params[indexset]
    candidate_function = tree_to_function(karva_to_tree(candidate_solution,candidate_parameters))
    candidate_impedance = candidate_function(candidate_parameters,1000)

    fitness = candidate_impedance == target_impedance ? head_length(candidate_solution) : 10

    return fitness
end

function next_generation(Candidate_index_sets,fitnesses,circuit_bare,karva_params,operation_indexes,terminal_indexes,target_impedance,mutationRate=0.05)
    #SELECTION
    individual_size = length(karva_params)
    pop_size = length(fitnesses)
    elite_size = Int(ceil(pop_size*0.1))
    children = Candidate_index_sets[1:elite_size] #The elites are automatically accepted into the new generation (should they be exempt from mutation as well?)
    Mating_pool = get_parent_inds_from_tournamant(pop_size,fitnesses,elite_size)
    Mating_pool_size = length(Mating_pool)
    #BREEDING
    for i in 1:Mating_pool_size
        parent1,parent2 = Candidate_index_sets[[Mating_pool[i],Mating_pool[Mating_pool_size-i+1]]]
        push!(children,ordered_restricted_crossover(individual_size,parent1,parent2,operation_indexes,terminal_indexes))
    end
    #MUTATION
    children = [restricted_swap_mutation(child,mutationRate,operation_indexes,terminal_indexes) for child in children]
    #FITNESS EVAULATION
    child_fitnesses = [indexset_fitness(children[i],circuit_bare,karva_params,target_impedance) for i in 1:pop_size]
    order = sortperm(child_fitnesses)

    return children[order], child_fitnesses[order]
end

function get_parent_inds_from_tournamant(pop_size,fitnesses,elite_size)
    parents = []
    mating_pool_size = pop_size-elite_size
    tournament_size = Int(ceil(0.1*pop_size))
    for i in 1:mating_pool_size
    tournament = rand(1:pop_size,tournament_size)
    max_,ind = findmin(fitnesses[tournament])
    push!(parents,tournament[ind])
    end
    return parents
end

function restricted_swap_mutation(Candidate,mutationrate,operation_indicess,terminal_indicess)
    individual = copy(Candidate)
    individual_length = length(individual)
    swapWith = 0
    swapped = 0
    for swapped in 1:individual_length
        if rand() < mutationrate
            if swapped == 1
                swapWith = rand(findall(x->x in operation_indicess,individual))
            elseif swapped in [individual_length-1,individual_length]
                swapWith = rand(findall(x->x in terminal_indicess,individual))
            else
                if swapped in findall(x->x in operation_indicess,individual)
                    swapWith = rand(1:individual_length-2)
                else
                    swapWith = rand(2:individual_length)
                end

            end
            gene1 = individual[swapped]
            gene2 = individual[swapWith]
            individual[swapped] = gene2
            individual[swapWith] = gene1
        end
    end
    return individual
end

function ordered_crossover(size,parent1,parent2)
    geneA = Int(ceil(rand()*size))
    geneB = Int(ceil(rand()*size))
    StartGene = min(geneA,geneB)
    EndGene = max(geneA,geneB)
    child = Int[]
    parent1_range = StartGene:EndGene
    parent1 = parent1[StartGene:EndGene]
    parent2 = parent2[findall(x->!(x in parent1),parent2)]
    for i in 1:size
        if i in parent1_range
            push!(child,popfirst!(parent1))
        else
            push!(child,popfirst!(parent2))
        end
    end
    return child
end

function ordered_restricted_crossover(size,parent1,parent2,operation_indexes,terminal_indexes)
    valid_crossover = false
    child = ordered_crossover(size,parent1,parent2)
    while !valid_crossover
        if child[1] in operation_indexes && child[end-1] in terminal_indexes && child[end] in terminal_indexes
            valid_crossover = true
        else
            child = ordered_crossover(size,parent1,parent2)
        end
    end
    return child
end

function circuit_to_karva(circuit,circuit_parameters,head=8)
    karva = ""
    params = []
    if length(foldl(replace,["["=>"","]"=>"","-"=>"+",","=>"-"], init = denumber_circuit(circuit))) <= 9
        karva,params = exhaustive_search(circuit,circuit_parameters)
    else
        karva,params = genetic_algorithm(circuit,circuit_parameters)
    end
    amount_to_extend = (2*head+1) - length(karva)
    extension_karva = String(rand("RCLP",amount_to_extend))
    extension_params = karva_parameters(extension_karva)
    return karva*extension_karva , vcat(params,extension_params)
end

function get_tree_parameters(tree)
    return [node.Parameter for node in tree]
end

function tree_to_circuit_with_inds(tree)
    nodecount = length(tree)
    essential_info = [[node.ParentIndex,node.Type,node.Index,[node.Parameter]] for node in tree]
    for i in 1:((length(essential_info)-1)/2)
        parent1,type1,index1,params1 = pop!(essential_info)
        parent2,type2,index2,params2 = pop!(essential_info)
        essential_info[parent1][2] = essential_info[parent1][2] == '+' ? type1*'-'*type2 : '['*type1*','*type2*']'
        essential_info[parent1][4] = vcat(params1,params2)
        essential_info[parent1][3] = vcat(index1,index2)
    end
return  number_circuit(essential_info[1][2]),  essential_info[1][4] , essential_info[1][3]
end


text:
using DSUtils
using Documenter

makedocs(;
    modules=[DSUtils],
    authors="Daymond Ling",
    repo="https://github.com/DaymondLing/DSUtils.jl/blob/{commit}{path}#L{line}",
    sitename="DSUtils.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DaymondLing.github.io/DSUtils.jl/stable",
        assets=String[],
    ),
    pages = Any[
        "Home" => "index.md",
        "Manual" => Any[
            "Binary Classifier" => "man/bcdiag.md",
            "Data preparation" => "man/dataprep.md",
            "Binning methods" => "man/binning.md",
            "Categorical encoding" => "man/categorical.md",
        ],
        "Function Reference" => "Reference.md",
    ],
)

deploydocs(;
   repo="github.com/DaymondLing/DSUtils.jl.git",
   devbranch = "main",
)


text:
"Calling Dex from Julia"
module DexCall
    using CombinedParsers
    using CombinedParsers.Regexp

    export evaluate, DexError, DexModule, juliaize, NativeFunction, @dex_func_str

    include("api_types.jl")
    include("api.jl")
    include("evaluate.jl")
    include("native_function.jl")

    # use this to disable free'ing haskell objects after we have closed the RTS
    const NO_FREE = Ref(false)

    function __init__()
        init()
        @eval const JIT = create_JIT()
        atexit() do
            destroy_JIT(JIT)
            NO_FREE[] = true
            fini()
        end

        @eval const PRELUDE = create_context()
    end
end

max_stars_repo_name: AlexDikelsky/dex-lang

max_issues_repo_name: AlexDikelsky/dex-lang

max_forks_repo_name: AlexDikelsky/dex-lang


text:
"""
    evaulate(model::PSGeometricBrownianMotionModelParameters)

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function evaluate(model::PSGeometricBrownianMotionModelParameters, initial_condition::Float64, tspan::Tuple{Float64,Float64}, timeStep::Float64;
    number_of_trials::Int64=10000, return_time_step::Float64=1.0)::PSResult

    try

        # get parameters from the model -
        μ = model.μ
        σ = model.σ

        # setup functions for the model equations -
        f(S, p, t) = μ * S
        g(S, p, t) = σ * S

        # setup the problem -
        problem = SDEProblem(f, g, initial_condition, tspan)

        # how many time steps do we have?
        T = collect(tspan[1]:return_time_step:tspan[2])
        number_of_times_steps = length(T)

        # initilize -
        X = zeros(number_of_times_steps, number_of_trials)

        # solve the model number_of_trials times, and then return the mean and stdev -
        for trial_index = 1:number_of_trials

            # solve the problem -
            solution = solve(problem, EM(), dt=timeStep, saveat=return_time_step)
            for step_index = 1:number_of_times_steps

                # grab the solution 0
                soln_array = solution.u[step_index]
                X[step_index, trial_index] = soln_array[1]
            end
        end

        # compute the mean, and std -
        μ = mean(X, dims=2)
        σ = std(X, dims=2)

        # create a named tuple to return -
        return_tuple = (T = T, X = X, μ = μ, σ = σ)

        # return -
        return PSResult(return_tuple)
    catch error
        return PSResult(error)
    end
end

"""
    evaulate(model::PSHestonAssetPricingModelParameters)

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function evaluate(model::PSHestonAssetPricingModelParameters, initialCondition::Array{Float64,1}, tspan::Tuple{Float64,Float64}, timeStep::Float64;
    number_of_trials::Int64=10000, return_time_step::Float64=1.0)::PSResult

    try

        # setup the problem -
        problem = HestonProblem(model.μ, model.𝝹, model.ϴ, model.σ, model.𝜌, initialCondition, tspan)

        # solve the problem -
        solution = solve(problem, EM(), dt=timeStep, saveat=return_time_step)

        # how many time steps do we have?
        T = collect(tspan[1]:return_time_step:tspan[2])
        number_of_times_steps = length(T)

        # initilize -
        X = zeros(number_of_times_steps, 1)
        for step_index = 1:number_of_times_steps

            # grab the solution 0
            soln_array = solution.u[step_index]
            X[step_index] = soln_array[1]
        end

        # compute the mean, and std -
        μ = mean(X, dims=2)
        σ = std(X, dims=2)

        # create a named tuple to return -
        return_tuple = (T = T, X = X, μ = μ, σ = σ)

        # return -
        return PSResult(return_tuple)
    catch error
        return PSResult(error)
    end
end

function evaluate(model::PSSingleIndexModelParameters, factorArray::Array{Float64,1};
    number_of_samples::Int64=100)::PSResult

    try

        # get parameters from model -
        ⍺ = model.⍺
        β = model.β
        𝝐 = model.𝝐
        riskFreeRate = model.riskFreeRate
        pV = [⍺ ; β]

        # initialize -
        number_of_steps = length(factorArray)
        X = zeros(number_of_steps, 2)
        Y = zeros(number_of_steps, number_of_samples)

        # formulate X -
        for step_index = 1:number_of_steps
            X[step_index,1] = 1.0
            X[step_index,2] = (factorArray[step_index,1] - riskFreeRate)
        end

        # formulate Y -
        for sample_index = 1:number_of_samples

            # compute delta -
            d = rand(𝝐, number_of_steps)

            # compute tmp -
            Ytmp = X * pV .+ d

            # cache -
            for step_index = 1:number_of_steps
                Y[step_index,sample_index] = Ytmp[step_index]
            end
        end

        # # compute stats from the samples -
        μ = mean(Y, dims=2)
        σ = std(Y, dims=2)

        # create a named tuple to return -
        return_tuple = (Y = Y, μ = μ, σ = σ)

        # return -
        return PSResult(return_tuple)
    catch error
        return PSResult(error)
    end
end

function price(model::PSSingleIndexModelParameters, factorArray::Array{Float64,1}, initialPrice::Float64;
    number_of_samples::Int64=100)::PSResult

    try

        # compute the returns -
        return_tuple = evaluate(model, factorArray; number_of_samples=number_of_samples) |> checkresult
        R = return_tuple.Y
        μ = return_tuple.μ
        σ = return_tuple.σ

        # compute return actual (not excess return) -
        R_actual = R .+ model.riskFreeRate

        # get the size -
        (number_of_rows, number_of_cols) = size(R_actual)

        # initialize -
        price_array = zeros(number_of_rows, number_of_samples)

        # samples -
        for sample_index = 1:number_of_samples

            # add the initial value -
            price_array[1,sample_index] = initialPrice

            # main loop -
            price_value = initialPrice
            for step_index = 2:number_of_rows

                # compute new price -
                theta_value = R_actual[step_index - 1, sample_index]
                tmp = price_value * (1 + theta_value)

                # grab -
                price_array[step_index, sample_index] = tmp

                # update -
                price_value = tmp
            end
        end

        # return -
        return PSResult(price_array)
    catch error
        return PSResult(error)
    end
end

function price(model::PSRandomWalkModelParameters, initialPrice::Float64, number_of_steps::Int64; number_of_samples::Int64=100)::PSResult

    # initialize -
    computed_price_array = zeros(number_of_steps, number_of_samples)

    try

        # draw samples -
        d_matrix = rand(model.𝝐, number_of_steps, number_of_samples)

        # main loop -
        for sample_index = 1:number_of_samples

            # set the initial price -
            computed_price_array[1,sample_index] = initialPrice

            # walk through time -
            for step_index = 2:number_of_steps
                value = computed_price_array[step_index - 1,sample_index] + d_matrix[step_index,sample_index]
                computed_price_array[step_index,sample_index] = value
            end
        end

        # return -
        return PSResult(computed_price_array)
    catch error
        return PSResult(error)
    end
end

function estimate_single_index_model(assetReturnArray::Array{Float64,1}, factorArray::Array{Float64,1};
    riskFreeRate::Float64=0.00169)::PSResult

    try

        # create X array -
        number_of_time_steps = length(assetReturnArray)
        X = zeros(number_of_time_steps, 2)
        for time_index = 1:number_of_time_steps
            X[time_index,1] = 1.0
            X[time_index,2] = (factorArray[time_index,1] - riskFreeRate)
        end

        # Rename the assetReturn array -
        Y = assetReturnArray .- riskFreeRate

        # compute the parameters theta -
        theta = (inv(transpose(X) * X)) * transpose(X) * Y

        # ok, so now that we have the parameters, let's compute the residual distribution -
        # TODO: Pass in the distribution type?
        Ymodel = X * theta
        residual_array = (Y - Ymodel)
        D = fit(Laplace, residual_array)

        # create a model wrapper -
        ⍺ = theta[1]
        β = theta[2]
        model_wrapper = PSSingleIndexModelParameters(⍺, β, riskFreeRate, D)

        return PSResult(model_wrapper)
    catch error
        return PSResult(error)
    end
end

function estimate_random_walk_model(assetPriceArray::Array{Float64,1})::PSResult

    # initialize -
    price_delta_array = Array{Float64,1}()
    number_of_time_steps = length(assetPriceArray)

    try

        # compute the price difference array -
        for time_index = 2:number_of_time_steps
            value = assetPriceArray[time_index] - assetPriceArray[time_index - 1]
            push!(price_delta_array, value)
        end

        # build the model -
        model = kde_lscv(price_delta_array)

        # build wrapper and return -
        return PSResult(PSRandomWalkModelParameters(model))
    catch error
        return PSResult(error)
    end
end

"""
    sample_underlying_gbm_model

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function sample_underlying_gbm_model(model::PSGeometricBrownianMotionModelParameters, initialUnderlyingPrice::Float64, tspan::Tuple{Float64,Float64}, timeStep::Float64;
    number_of_trials::Int64=10000, return_time_step::Float64=1.0)

    # initialize -
    Random.seed!(123) # Setting the seed
    time_sample_array = range(tspan[1], stop=tspan[2], step=timeStep)
    number_of_time_steps = length(time_sample_array)
    sample_price_array = zeros(2 * number_of_trials, (number_of_time_steps + 1))

    try

        # setup the random number -
        d = Normal()

        # generate samples -
        positive_array = rand(d, number_of_time_steps, number_of_trials)
        negative_array = -1 * positive_array
        Z_sample_array = [positive_array negative_array]

        # grab stuff from the model -
        μ = model.μ
        σ = model.σ
        Δt = timeStep

        # add the initial condition -
        for sample_index = 1:(2 * number_of_trials)
            sample_price_array[sample_index,1] = initialUnderlyingPrice
        end

        # setup the time array -
        for sample_index = 1:(2 * number_of_trials)
            for time_index = 1:number_of_time_steps

                # grab Z value -
                Z_value = Z_sample_array[time_index,sample_index]

                # compute new state -
                soln_term = (μ - 0.5 * σ^2) * Δt + σ * (sqrt(Δt)) * Z_value
                old_state = sample_price_array[sample_index,time_index]
                new_state = old_state * exp(soln_term)

                # cache -
                sample_price_array[sample_index,time_index + 1] = new_state
            end
        end

        # return -
        return PSResult(sample_price_array)
    catch error
        return PSResult(error)
    end
end


text:
using SimJulia

type StoreObject
  i :: Int
end

function consumer(sim::Simulation, sto::Store)
  for i = 1:10
    yield(sim, timeout(sim, rand()))
    println("$(now(sim)), consumer is demanding object")
    obj = yield(sim, get(sim, sto))
    println("$(now(sim)), consumer is being served with object $(obj.i)")
  end
end

function producer(sim::Simulation, sto::Store)
  for i = 1:10
    println("$(now(sim)), producer is offering object $i")
    yield(sim, put(sim, sto, StoreObject(i)))
    println("$(now(sim)), producer is being served")
    yield(sim, timeout(sim, 2*rand()))
  end
end

sim = Simulation()
sto = Store(StoreObject)
Process(sim, consumer, sto)
Process(sim, producer, sto)
run(sim)


text:
using LightGraphs
using DataFrames
using Random
using Statistics
using JLD
using JLD2
using StatsBase
using CSV
using GraphIO

for script in readdir("src")
    if script != "simulation.jl"
        include(script)
    end
end

mutable struct Simulation

    config::Config
    init_state::Any
    final_state::Any
    agent_log::DataFrame
    post_log::Any
    graph_list::Array{AbstractGraph}

    function Simulation(config=Config())
        new(
            config,
            (nothing, nothing),
            (nothing, nothing),
            DataFrame(),
            DataFrame(),
            Array{AbstractGraph, 1}(undef, 0)
        )
    end

end

function tick!(
    state::Tuple{AbstractGraph, AbstractArray}, post_list::AbstractArray,
    tick_nr::Int64, config::Config
)
    agent_list = state[2]
    for agent_idx in shuffle(1:length(agent_list))
        this_agent = agent_list[agent_idx]
        update_feed!(state, agent_idx, config)
        update_perceiv_publ_opinion!(state, agent_idx)
        update_opinion!(state, agent_idx, config)

        drop_friends!(state, agent_idx, config)
        if indegree(state[1], agent_idx) < config.network.m0
            if config.simulation.addfriends == "neighborsofneighbors"
                add_friends_neighbors_of_neighbors!(
                    state, agent_idx, post_list,
                    config, config.network.new_follows
                )
            elseif config.simulation.addfriends == "random"
                add_friends_random!(
                    state, agent_idx, post_list,
                    config, config.network.new_follows
                )
            else
                add_friends_neighbors_of_neighbors!(
                    state, agent_idx, post_list,
                    config, floor(Int,config.network.new_follows/2)
                )
                add_friends_random!(
                    state, agent_idx, post_list,
                    config, floor(Int,config.network.new_follows/2)
                )
            end
        end

        inclin_interact = deepcopy(this_agent.inclin_interact)
        while inclin_interact > 0
            if rand() < inclin_interact
                publish_post!(state, post_list, agent_idx, tick_nr)
            end
            inclin_interact -= 1.0
        end
    end
    return log_network(state, tick_nr)
end

function run!(
    simulation::Simulation=Simulation();
    name::String = "result"
    )

    config = simulation.config

    graph = DiGraph(barabasi_albert(
                    config.network.agent_count,
                    config.network.m0))
    simulation.init_state = (graph, create_agents(graph))
    state = deepcopy(simulation.init_state)
    post_log = Array{Post, 1}(undef, 0)
    simulation.graph_list = Array{AbstractGraph, 1}([simulation.init_state[1]])

    agent_log = DataFrame(
        TickNr = Int64[],
        AgentID = Int64[],
        Opinion = Float64[],
        PerceivPublOpinion = Float64[],
        Indegree = Int64[],
        Outdegree = Int64[],
        Centrality = Float64[],
        CC = Float64[],
        Component = Int64[]
    )

    if !in("tmp", readdir())
        mkdir("tmp")
    end
    if name == "result"
        print("Current Tick: 0")
    end

    for i in 1:config.simulation.ticks

        if name == "result"
            print('\r')
            print("Current Tick: $i, current AVG agents connection count::" * string(round(mean(degree(state[1])))) * ", max outdegree: " * string(maximum(outdegree(state[1]))) * ", current Posts: " * string(length([post for post in post_log if length(post.seen_by) > 0])))
        end
        append!(agent_log, tick!(state, post_log, i, config))
        if i % ceil(config.simulation.ticks / 10) == 0
            if name != "result"
                print(".")
            end

            push!(simulation.graph_list, deepcopy(state[1]))

            simulation.final_state = state
            simulation.agent_log = agent_log
            simulation.post_log = post_log

            save(joinpath("tmp", name * ".jld2"), string(i), simulation)
        end
    end

    simulation.final_state = state
    simulation.agent_log = agent_log
    simulation.post_log = DataFrame(
        Opinion = [p.opinion for p in post_log],
        Weight = [p.weight for p in post_log],
        Source_Agent = [p.source_agent for p in post_log],
        Published_At = [p.published_at for p in post_log],
        Seen = [p.seen_by for p in post_log]
    )

    if !in("results", readdir())
        mkdir("results")
    end
    save(joinpath("results", name * ".jld2"), name, simulation)
    rm(joinpath("tmp", name * ".jld2"))

    print("\n---\nFinished simulation run with the following specifications:\n $config\n---\n")

    return simulation
end

function run_batch(
    configlist::Array{Config, 1};
    resume_at::Int64=1,
    stop_at::Int64=length(configlist),
    batch_name::String = ""
    )

    for i in resume_at:stop_at
        run_nr = lpad(string(i),length(string(length(configlist))),"0")
        run!(
            Simulation(configlist[i]),
            name = (batch_name * "_run$run_nr")
        )
    end
end

function run_resume!(
    path::String = "result"
)

    if !("tmp" in readdir())
        mkdir("tmp")
    elseif "tmp" in readdir() && path == "result"
        path = joinpath("tmp", readdir("tmp")[1])
    end

    raw_data = load(path)

    name = (
        path[first(findlast("\\", path))+1:first(findfirst(".jld2", path))-1]
    )

    tick_nr = parse(Int, first(keys(raw_data))) + 1
    simulation = collect(values(raw_data))[1]
    config = simulation.config

    println("Resumed Run: $name \n $config) \n ---")

    for i in 1:Int((ticknr - 1) / config.simulation.ticks * 10)
        print(".")
    end

    state = simulation.final_state
    agent_log = simulation.agent_log
    post_log = simulation.post_log

    for i in tick_nr:config.simulation.ticks

        if name == "result"
            print('\r')
            print(
                "Current Tick: $i, current AVG agents connection count::"
                * string(round(mean(degree(state[1])))) * ", max outdegree: "
                * string(maximum(outdegree(state[1]))) * ", current Posts: "
                * string(length(
                    [post for post in post_log if length(post.seen_by) > 0])
                )
            )
        end

        append!(agent_log, tick!(state, post_log, i, config))

        if i % ceil(config.simulation.ticks / 10) == 0
            if name != "result"
                print(".")
            end

            push!(simulation.graph_list, deepcopy(state[1]))

            simulation.final_state = state
            simulation.agent_log = agent_log
            simulation.post_log = post_log

            save(joinpath("tmp", name * ".jld2"), string(i), simulation)
        end
    end

    simulation.final_state = state
    simulation.agent_log = agent_log
    simulation.post_log = DataFrame(
        Opinion = [p.opinion for p in post_log],
        Weight = [p.weight for p in post_log],
        Source_Agent = [p.source_agent for p in post_log],
        Published_At = [p.published_at for p in post_log],
        Seen = [p.seen_by for p in post_log]
    )

    if !in("results", readdir())
        mkdir("results")
    end
    save(joinpath("results", name * ".jld2"), name, simulation)
    rm(joinpath("tmp", name * ".jld2"))

    if length(readdir("tmp")) == 0
        rm("tmp")
    end

    print(
        "\n---\nFinished simulation run with the following specifications:\n
        $(simulation.config)\n---\n"
    )

    return simulation

end

# suppress output of include()
;


text:
using Test
using SafeTestsets

@safetestset "FEM Matrix assembly" begin include("fem_assembly.jl") end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# # Triangulation wombling

using SpatialBoundaries
using StatsPlots

# Plot defaults

default(; dpi=500, size=(600, 600), aspectratio=1, c=:davos, frame=:box)

# Get some points at random

n = 500
x = rand(n)
y = rand(n)
z = [(x[i]<=0.5)&(y[i]<=0.5) ? rand() : rand().+1.2 for i in eachindex(x)]

scatter(x, y, marker_z = z, lab="")

# Womble

W = wombling(x, y, z)

# Get the rate of change

scatter(x, y, c=:lightgrey, msw=0.0, lab="", m=:square, ms=3)
scatter!(W.x, W.y, marker_z = log1p.(W.m), lab="")

# Angle histogram

stephist(
    deg2rad.(sort(vec(W.θ)));
    proj=:polar,
    lab="",
    c=:teal,
    fill=(0, 0.2, :teal),
    nbins=100,
)

# Show the rotation with a color

scatter(W.x, W.y, marker_z = W.θ, c=:vik, clim=(0, 360))


text:
using BlackBoxOptim

p = BlackBoxOptim.example_problems["Ackley"]

BlackBoxOptim.repeated_bboptimize(5, p, 100, [
  :generating_set_search,
  :probabilistic_descent,
  :adaptive_de_rand_1_bin_radiuslimited,
  :random_search,
  :separable_nes],
  20.0)


text:
module QuWalkEC

using LightGraphs
using IterTools

export have_stationary_state

export
  have_ss_connected,
  have_ss,
  search_subgraph_ss,
  search_subgraph_noss

include("ec_search.jl")
include("ec_check.jl")

end

max_stars_repo_name: QuantumWalks/ExceptionalConfiguration

max_issues_repo_name: QuantumWalks/ExceptionalConfiguration

max_forks_repo_name: QuantumWalks/ExceptionalConfiguration


text:
@time using DamagedShearBand
const DSB = DamagedShearBand
@time using Plots
using Printf
using LaTeXStrings
##
# imposed values
ϵ̇₁₁ = -1e-5
ϵ̇ⁱξη = 1e-5
σ₃ = -1e6 # negative in compression

# initial higher compressive stress multiplier (S = σ₁₁/σ₂₂)
S_i = 1.0
θ = 60.0

# Create a Rheology type instance containing elastic moduli and damage parameters. Change default values by  supplying keywords arguments.
r = DSB.Rheology(D₀=0.2,n=5.0)

# initial damage
#D_i = r.D₀
Dᵒ, _ = DSB.get_KI_mininizer_D_on_S_range(r,S_i,100,σ₃)
#Dⁱ = min(Dᵒ+0.15, 0.9)
D_i = Dᵒ

# time integration parameters :
ϵⁱξη_goal = 0.01
tspan = (0.0,ϵⁱξη_goal/ϵ̇ⁱξη) # initial and final simulation time
tspan = (0.0,141.26)
Δt_i = 0.01 # initial timestep

# Build stress and strain tensors
# σᵢⱼ_i = DSB.build_principal_stress_tensor(r,S_i,σ₃,D_i ; abstol=1e-15) # takes care of the plane strain constraint by solving non linear out of plane strain wrt σ₃₃ using Newton algorithm
# ϵᵢⱼ_i = DSB.compute_ϵij(r,D_i,σᵢⱼ_i)

n_iter_saved = 20000
n_iter_printed = 10
op = DSB.OutputParams(save_frequency = n_iter_saved/(tspan[2]-tspan[1]),
                      save_period = 1000,
                      print_frequency = n_iter_printed/(tspan[2]-tspan[1]),
                      print_period = 1000)

sp = DSB.SolverParams(newton_abstol = 1e-12,
                      newton_maxiter = 100,
                      time_maxiter = nothing,
                      e₀ = (D=1e-2, S=1e-2, σ=100.0, ϵ=1e-2))
# sp = DSB.SolverParams(newton_abstol = 1e-12,
#                       newton_maxiter = 100,
#                       time_maxiter = nothing,
#                       e₀ = (D=1e-3, S=1e-3, σ=1.0, ϵ=1e-5))

p = DSB.Params(op,sp)
# Integrate over time, by solving for the unknowns [σ₁₁next, σ₃₃next, ϵ₂₂next] at each timestep
# -- e₀ is the absolute tolerance used to adapt time stepping, may be a scalar or a NamedTuple with keys (D,σ,ϵ)
# -- abstol is the absolute tolerance on the unknowns in the Newton solver

#t_vec, S_vec, σⁱᵢⱼ_vec, σᵒᵢⱼ_vec, ϵⁱᵢⱼ_vec, Dⁱ_vec, Dᵒ_vec = DSB.adaptative_time_integration_2_points(r,p,σ₃,ϵ̇ⁱξη,Δt_i,θ,tspan);
t_vec, S_vec, σⁱᵢⱼ_vec, σᵒᵢⱼ_vec, ϵⁱᵢⱼ_vec, Dⁱ_vec, Dᵒ_vec = DSB.adaptative_time_integration_2_points(r,p,S_i,σ₃,Dᵒ,ϵ̇₁₁,ϵ̇ⁱξη,Δt_i,θ,tspan ; damage_growth_out=false, bifurcate_on=:KI)

#
σᵒ_p_vec = DSB.principal_coords.(σᵒᵢⱼ_vec,Ref(θ))
σⁱ_p_vec = DSB.principal_coords.(σⁱᵢⱼ_vec,Ref(θ))
σᵒ₁₁_p_vec = [σᵢⱼ[1,1] for σᵢⱼ in σᵒ_p_vec]
σᵒ₂₂_p_vec = [σᵢⱼ[2,2] for σᵢⱼ in σᵒ_p_vec]
σᵒ₃₃_p_vec = [σᵢⱼ[3,3] for σᵢⱼ in σᵒ_p_vec]

σᵒ₁₁_vec = [σᵢⱼ[1,1] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₁₁_vec = [σᵢⱼ[1,1] for σᵢⱼ in σⁱᵢⱼ_vec]
σᵒ₁₂_vec = [σᵢⱼ[1,2] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₁₂_vec = [σᵢⱼ[1,2] for σᵢⱼ in σⁱᵢⱼ_vec]
σᵒ₂₁_vec = [σᵢⱼ[2,1] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₂₁_vec = [σᵢⱼ[2,1] for σᵢⱼ in σⁱᵢⱼ_vec]
σᵒ₂₂_vec = [σᵢⱼ[2,2] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₂₂_vec = [σᵢⱼ[2,2] for σᵢⱼ in σⁱᵢⱼ_vec]
σᵒ₃₃_vec = [σᵢⱼ[3,3] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₃₃_vec = [σᵢⱼ[3,3] for σᵢⱼ in σⁱᵢⱼ_vec]

ϵⁱ₁₁_vec = [ϵᵢⱼ[1,1] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱ₁₂_vec = [ϵᵢⱼ[1,2] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱ₂₁_vec = [ϵᵢⱼ[2,1] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱ₂₂_vec = [ϵᵢⱼ[2,2] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱ₃₃_vec = [ϵᵢⱼ[3,3] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]

τⁱ_vec =  DSB.get_τ.(σⁱᵢⱼ_vec; input=:σ)
τᵒ_vec =  DSB.get_τ.(σᵒᵢⱼ_vec; input=:σ)
pⁱ_vec = -1/3 * tr.(σⁱᵢⱼ_vec)
pᵒ_vec = -1/3 * tr.(σⁱᵢⱼ_vec)
KIⁱ_vec = [DSB.compute_KI(r,σⁱᵢⱼ_vec[i], Dⁱ_vec[i]) for i in eachindex(σⁱᵢⱼ_vec)]
KIᵒ_vec = [DSB.compute_KI(r,σᵒᵢⱼ_vec[i], Dᵒ_vec[i]) for i in eachindex(σᵒᵢⱼ_vec)]
γⁱ_vec = [sqrt(2*ϵᵢⱼ⊡ϵᵢⱼ) for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱkk_vec = [tr(ϵᵢⱼ) for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
deviation_out_vec = DSB.get_stress_deviation_from_far_field.(σᵒ_p_vec ; offdiagtol=1e-5)
deviation_in_vec = DSB.get_stress_deviation_from_far_field.(σⁱ_p_vec ; offdiagtol=1e-5)

σᵒᵢⱼ_i_p, σⁱᵢⱼ_i_p, ϵᵒᵢⱼ_i_p, ϵⁱᵢⱼ_i_p = DSB.initialize_state_var_D(r,p,S_i,σ₃,Dⁱ,Dᵒ,θ ; coords=:principal)
σᵒᵢⱼ_i, σⁱᵢⱼ_i, ϵᵒᵢⱼ_i, ϵⁱᵢⱼ_i = DSB.initialize_state_var_D(r,p,S_i,σ₃,Dⁱ,Dᵒ,θ ; coords=:band)
deviation_out_check = DSB.get_stress_deviation_from_far_field(σᵒᵢⱼ_i_p ; offdiagtol=1e-12)
deviation_in_check = DSB.get_stress_deviation_from_far_field(σⁱᵢⱼ_i_p ; offdiagtol=1e-12)
deviation_out_eps = DSB.get_stress_deviation_from_far_field(ϵᵒᵢⱼ_i_p ; offdiagtol=1e-12)
deviation_in_eps = DSB.get_stress_deviation_from_far_field(ϵⁱᵢⱼ_i_p ; offdiagtol=1e-12)

deviation_out_check = DSB.get_stress_deviation_from_far_field(DSB.principal_coords(σᵒᵢⱼ_i,θ))
deviation_in_check = DSB.get_stress_deviation_from_far_field(DSB.principal_coords(σⁱᵢⱼ_i,θ))
deviation_out_eps = DSB.get_stress_deviation_from_far_field(DSB.principal_coords(ϵᵒᵢⱼ_i,θ))
deviation_in_eps = DSB.get_stress_deviation_from_far_field(DSB.principal_coords(ϵⁱᵢⱼ_i,θ))
## Plots ##
# figure size and backend
plotlyjs(size=(1000,300))

# S vs time
plot(t_vec,γⁱ_vec)
plot(t_vec,ϵⁱ₁₂_vec)
plot!(t_vec,S_vec./100)
plot(ϵⁱ₁₂_vec,S_vec)
# principal stress orientation deviation from far field vs time
plot(t_vec,deviation_in_vec)
plot!(t_vec,deviation_out_vec)
#plot!(t_vec,-(Dⁱ_vec .- Dⁱ).*60)

plot(deviation_in_vec)
# Damage vs time
plot(t_vec,Dⁱ_vec)
plot!(t_vec,Dᵒ_vec)
ylims!(0,1)
# KI/KIC vs time
plot(t_vec,KIⁱ_vec./r.K₁c)
plot!(t_vec,KIᵒ_vec./r.K₁c)

plot!(ϵⁱ₁₂_vec,KIⁱ_vec./r.K₁c)
plot!(ϵⁱ₁₂_vec,KIᵒ_vec./r.K₁c)
# shear stress on band vs time
plot(t_vec,σⁱ₁₂_vec)
plot!(t_vec,σᵒ₁₂_vec)
# vol/shear strain invariants vs time
plot(t_vec, ϵⁱkk_vec./γⁱ_vec)
plot(t_vec[1:end-1], diff(ϵⁱkk_vec)./diff(γⁱ_vec))

# plot τ in and out
plot(t_vec, τⁱ_vec)
plot!(t_vec, τᵒ_vec)
## composite plot :
ϵⁱ₁₂_vec_prev = ϵⁱ₁₂_vec
S_vec_prev = S_vec
Dⁱ_vec_prev = Dⁱ_vec
Dᵒ_vec_prev = Dᵒ_vec
deviation_in_vec_prev = deviation_in_vec
deviation_out_vec_prev = deviation_out_vec
KIⁱ_vec_prev = KIⁱ_vec
KIᵒ_vec_prev = KIᵒ_vec
##
#plotlyjs(size=(1000,600))
pyplot(size=(1000,600))
layout = @layout [a ; b ; c ; d]

p1 = plot(ϵⁱ₁₂_vec.*100,S_vec,c=:blue,label="KI controled")
plot!(p1,ϵⁱ₁₂_vec_prev.*100,S_vec_prev,c=:blue,ls=:dash,label="S controled")
#title!(p1,"D₀ = $(r.D₀) & D_in/out = ("*@sprintf("%.2f",Dⁱ)*", "*@sprintf("%.2f",Dᵒ)*") & gammadot = $(2*ϵ̇ⁱξη) & p=$(-σ₃)")
title!(p1,"D₀=$(r.D₀)  &  D_initial="*@sprintf("%.2f",Dⁱ)*"  &  gammadot=$(2*ϵ̇ⁱξη)  &  p=$(-σ₃)")
ylabel!("σᵒᵘᵗ₁ multiplier (S)")

p2 = plot(ϵⁱ₁₂_vec.*100,Dⁱ_vec,c=:blue,label="inside")
plot!(p2,ϵⁱ₁₂_vec.*100,Dᵒ_vec,c=:red,label="outside")
plot!(p2,ϵⁱ₁₂_vec_prev.*100,Dⁱ_vec_prev,c=:blue,ls=:dash,label=nothing)
plot!(p2,ϵⁱ₁₂_vec_prev.*100,Dᵒ_vec_prev,c=:red,ls=:dash,label=nothing)
ylabel!("Damage")

p3 = plot(ϵⁱ₁₂_vec.*100,deviation_in_vec,c=:blue,label="inside")
plot!(p3,ϵⁱ₁₂_vec.*100,deviation_out_vec,c=:red,label="outside")
plot!(p3,ϵⁱ₁₂_vec_prev.*100,deviation_in_vec_prev,c=:blue,ls=:dash,label=nothing)
plot!(p3,ϵⁱ₁₂_vec_prev.*100,deviation_out_vec_prev,c=:red,ls=:dash,label=nothing)
ylabel!("stress orientation in band")

p4 = plot(ϵⁱ₁₂_vec,KIⁱ_vec./r.K₁c,c=:blue,label="inside")
plot!(p4,ϵⁱ₁₂_vec,KIᵒ_vec./r.K₁c,c=:red,label="outside")
plot!(p4,ϵⁱ₁₂_vec_prev,KIⁱ_vec_prev./r.K₁c,c=:blue,ls=:dash,label=nothing)
plot!(p4,ϵⁱ₁₂_vec_prev,KIᵒ_vec_prev./r.K₁c,c=:red,ls=:dash,label=nothing)
ylabel!(L"K_{I}/K_{Ic}")
xlabel!("ϵⁱ₁₂ (%)")

pl = plot(p1,p2,p3,p4,layout=layout,lw=0.)

#savefig(pl,"./examples/figures/KI&S_controled_bifurcation_D0=$(r.D₀)_p=$(-σ₃)_ϵ̇ⁱξη=$(ϵ̇ⁱξη).pdf")


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

# Test string/show representation of Date
@test string(Dates.Date(1, 1, 1)) == "0001-01-01" # January 1st, 1 AD/CE
@test sprint(show, Dates.Date(1, 1, 1)) == "0001-01-01"
@test string(Dates.Date(0, 12, 31)) == "0000-12-31" # December 31, 1 BC/BCE
@test Dates.Date(1, 1, 1) - Dates.Date(0, 12, 31) == Dates.Day(1)
@test Dates.Date(Dates.UTD(-306)) == Dates.Date(0, 2, 29)
@test string(Dates.Date(0, 1, 1)) == "0000-01-01" # January 1st, 1 BC/BCE
@test string(Dates.Date(-1, 1, 1)) == "-0001-01-01" # January 1st, 2 BC/BCE
@test string(Dates.Date(-1000000, 1, 1)) == "-1000000-01-01"
@test string(Dates.Date(1000000, 1, 1)) == "1000000-01-01"
@test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 1)) == "2000-01-01T00:00:00.001"
@test sprint(show, Dates.DateTime(2000, 1, 1, 0, 0, 0, 1)) == "2000-01-01T00:00:00.001"
@test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 2)) == "2000-01-01T00:00:00.002"
@test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 500)) == "2000-01-01T00:00:00.5"
@test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 998)) == "2000-01-01T00:00:00.998"
@test string(Dates.DateTime(2000, 1, 1, 0, 0, 0, 999)) == "2000-01-01T00:00:00.999"
@test string(Dates.Time(0)) == "00:00:00"
@test string(Dates.Time(0, 1)) == "00:01:00"
@test string(Dates.Time(0, 1, 2)) == "00:01:02"
@test string(Dates.Time(0, 1, 2, 3)) == "00:01:02.003"
@test string(Dates.Time(0, 1, 2, 3, 4)) == "00:01:02.003004"
@test string(Dates.Time(0, 1, 2, 3, 4, 5)) == "00:01:02.003004005"
@test string(Dates.Time(0, 0, 0, 0, 1)) == "00:00:00.000001"
@test string(Dates.Time(0, 0, 0, 0, 0, 1)) == "00:00:00.000000001"
@test string(Dates.Time(0, 0, 0, 1)) == "00:00:00.001"

# DateTime parsing
# Useful reference for different locales: http://library.princeton.edu/departments/tsd/katmandu/reference/months.html

# Allow parsing of strings which are not representable as a TimeType
let str = "02/15/1996 24:00", df = Dates.DateFormat("mm/dd/yyyy HH:MM")
    parsed = Any[
        Dates.Month(2), Dates.Day(15), Dates.Year(1996), Dates.Hour(24), Dates.Minute(0)
    ]
    @test Dates.parse_components(str, df) == parsed
    @test_throws ArgumentError Dates.parse(DateTime, str, df)
end

# DateFormat printing

@test sprint(show, DateFormat("yyyzzxmmdd\\MHH:MM:SS\\P")) == "dateformat\"yyyzzxmmdd\\MHH:MM:SSP\""
@test sprint(show, DateFormat("yyy").tokens[1]) == "DatePart(yyy)"
@test sprint(show, DateFormat("mmzzdd").tokens[2]) == "Delim(zz)"
@test sprint(show, DateFormat("ddxmm").tokens[2]) == "Delim(x)"
@test sprint(show, DateFormat("xxmmxx").tokens[2]) == "DatePart(mm)"

# Common Parsing Patterns
#'1996-January-15'
dt = Dates.DateTime(1996, 1, 15)
f = "yy-mm-dd"
a = "96-01-15"
@test DateTime(a, f) + Dates.Year(1900) == dt
@test Dates.format(dt, f) == a
a1 = "96-1-15"
@test Dates.DateTime(a1, f) + Dates.Year(1900) == dt
@test Dates.format(dt, "yy-m-dd") == a1
a2 = "96-1-1"
@test Dates.DateTime(a2, f) + Dates.Year(1900) + Dates.Day(14) == dt
@test Dates.format(dt - Dates.Day(14), "yy-m-d") == a2
a3 = "1996-1-15"
@test Dates.DateTime(a3, f) == dt
@test Dates.format(dt, "yyyy-m-d") == a3
a4 = "1996-Jan-15"
@test_throws ArgumentError Dates.DateTime(a4, f) # Trying to use month name, but specified only "mm"

f = "yy/uuu/dd"
b = "96/Feb/15"
@test Dates.DateTime(b, f) + Dates.Year(1900) == dt + Dates.Month(1)
@test Dates.format(dt + Dates.Month(1), f) == b
b1 = "1996/Feb/15"
@test Dates.DateTime(b1, f) == dt + Dates.Month(1)
@test Dates.format(dt + Dates.Month(1), "yyyy/uuu/dd") == b1
b2 = "96/Feb/1"
@test Dates.DateTime(b2, f) + Dates.Year(1900) + Dates.Day(14) == dt + Dates.Month(1)
@test Dates.format(dt + Dates.Month(1) - Dates.Day(14), "yy/uuu/d") == b2
# Here we've specifed a text month name, but given a number
b3 = "96/2/15"
@test_throws ArgumentError Dates.DateTime(b3, f)
try
    Dates.parse(DateTime, "2012/2/20T9:9:31.25i90", dateformat"yyyy/mm/ddTHH:MM:SS.s")
    @test false
catch err
    @test isa(err, ArgumentError)
    @test err.msg == "Found extra characters at the end of date time string"
end
try
    Dates.parse(DateTime, "2012/2/20T9:9:3i90", dateformat"yyyy/mm/ddTHH:MM:SS.s")
    @test false
catch err
    @test isa(err, ArgumentError)
    @test err.msg == "Unable to parse date time. Expected directive Delim(.) at char 16"
end

f = "yy:dd:mm"
c = "96:15:01"
@test Dates.DateTime(c, f) + Dates.Year(1900) == dt
@test Dates.format(dt, f) == c
c1 = "1996:15:01"
@test Dates.DateTime(c1, f) == dt
@test Dates.format(dt, "yyyy:dd:mm") == c1
c2 = "96:15:1"
@test Dates.DateTime(c2, f) + Dates.Year(1900) == dt
@test Dates.format(dt, "yy:dd:m") == c2
c3 = "96:1:01"
@test Dates.DateTime(c3, f) + Dates.Year(1900) + Dates.Day(14) == dt
@test Dates.format(dt - Dates.Day(14), "yy:m:dd") == c3
c4 = "1996:15:01 # random comment"
@test_throws ArgumentError Dates.DateTime(c4, f)

f = "yyyy, uuu, dd"
d = "1996, Jan, 15"
@test Dates.DateTime(d, f) == dt
@test Dates.format(dt, f) == d
d1 = "96, Jan, 15"
@test Dates.DateTime(d1, f) + Dates.Year(1900) == dt
@test Dates.format(dt, "yy, uuu, dd") == d1
d2 = "1996, Jan, 1"
@test Dates.DateTime(d2, f) + Dates.Day(14) == dt
@test Dates.format(dt - Dates.Day(14), "yyyy, uuu, d") == d2
d3 = "1996, 2, 15"
@test_throws ArgumentError Dates.DateTime(d3, f)

f = "yyyy.U.dd"
e = "1996.January.15"
@test Dates.DateTime(e, f) == dt
@test Dates.format(dt, f) == e
e1 = "96.January.15"
@test Dates.DateTime(e1, f) + Dates.Year(1900) == dt
@test Dates.format(dt, "yy.U.dd") == e1

fo = "yyyy m dd"
f = "1996 1 15"
@test Dates.DateTime(f, fo) == dt
@test Dates.format(dt, fo) == f
f1 = "1996 01 15"
@test Dates.DateTime(f1, fo) == dt
@test Dates.format(dt, "yyyy mm dd") == f1
f2 = "1996 1 1"
@test Dates.DateTime(f2, fo) + Dates.Day(14) == dt
@test Dates.format(dt - Dates.Day(14), "yyyy m d") == f2

j = "1996-01-15"
f = "yyyy-mm-dd zzz"
@test Dates.DateTime(j, f) == dt
@test Dates.format(dt, f) == j * " zzz"
k = "1996-01-15 10:00:00"
f = "yyyy-mm-dd HH:MM:SS zzz"
@test Dates.DateTime(k, f) == dt + Dates.Hour(10)
@test Dates.format(dt + Dates.Hour(10), f) == k * " zzz"
l = "1996-01-15 10:10:10.25"
f = "yyyy-mm-dd HH:MM:SS.ss zzz"
@test Dates.DateTime(l, f) == dt + Dates.Hour(10) + Dates.Minute(10) + Dates.Second(10) + Dates.Millisecond(250)
@test Dates.format(dt + Dates.Hour(10) + Dates.Minute(10) + Dates.Second(10) + Dates.Millisecond(250), f) == l * " zzz"

r = "1/15/1996" # Excel
f = "m/dd/yyyy"
@test Dates.DateTime(r, f) == dt
@test Dates.format(dt, f) == r
s = "19960115"
f = "yyyymmdd"
@test Dates.DateTime(s, f) == dt
@test Dates.format(dt, f) == s
v = "1996-01-15 10:00:00"
f = "yyyy-mm-dd HH:MM:SS"
@test Dates.DateTime(v, f) == dt + Dates.Hour(10)
@test Dates.format(dt + Dates.Hour(10), f) == v
w = "1996-01-15T10:00:00"
f = "yyyy-mm-ddTHH:MM:SS zzz"
@test Dates.DateTime(w, f) == dt + Dates.Hour(10)
@test Dates.format(dt + Dates.Hour(10), f) == w * " zzz"

f = "yyyy/m"
y = "1996/1"
@test Dates.DateTime(y, f) == dt - Dates.Day(14)
@test Dates.format(dt, f) == y
y1 = "1996/1/15"
@test_throws ArgumentError Dates.DateTime(y1, f)
y2 = "96/1"
@test Dates.DateTime(y2, f) + Dates.Year(1900) == dt - Dates.Day(14)
@test Dates.format(dt, "yy/m") == y2

f = "yyyy"
z = "1996"
@test Dates.DateTime(z, f) == dt - Dates.Day(14)
@test Dates.format(dt, f) == z
z1 = "1996-3"
@test_throws ArgumentError Dates.DateTime(z1, f)
z2 = "1996-3-1"
@test_throws ArgumentError Dates.DateTime(z2, f)

aa = "1/5/1996"
f = "m/d/yyyy"
@test Dates.DateTime(aa, f) == dt - Dates.Day(10)
@test Dates.format(dt - Dates.Day(10), f) == aa
bb = "5/1/1996"
f = "d/m/yyyy"
@test Dates.DateTime(bb, f) == dt - Dates.Day(10)
@test Dates.format(dt - Dates.Day(10), f) == bb
cc = "01151996"
f = "mmddyyyy"
@test Dates.DateTime(cc, f) == dt
@test Dates.format(dt, f) == cc
dd = "15011996"
f = "ddmmyyyy"
@test Dates.DateTime(dd, f) == dt
@test Dates.format(dt, f) == dd
ee = "01199615"
f = "mmyyyydd"
@test Dates.DateTime(ee, f) == dt
@test Dates.format(dt, f) == ee
ff = "1996-15-Jan"
f = "yyyy-dd-uuu"
@test Dates.DateTime(ff, f) == dt
@test Dates.format(dt, f) == ff
gg = "Jan-1996-15"
f = "uuu-yyyy-dd"
@test Dates.DateTime(gg, f) == dt
@test Dates.format(dt, f) == gg
hh = "1996#1#15"
f = "yyyy#m#d"
@test Dates.DateTime(hh, f) == dt
@test Dates.format(dt, f) == hh

# test prefix.
s = "/1996/1/15"
f = "/yyyy/m/d"
@test Dates.DateTime(s, f) == dt
@test Dates.format(dt, f) == s
@test_throws ArgumentError Dates.DateTime("1996/1/15", f)

# from Jiahao
@test Dates.Date("2009年12月01日", "yyyy年mm月dd日") == Dates.Date(2009, 12, 1)
@test Dates.format(Dates.Date(2009, 12, 1), "yyyy年mm月dd日") == "2009年12月01日"
@test Dates.Date("2009-12-01", "yyyy-mm-dd") == Dates.Date(2009, 12, 1)

# French: from Milan
f = "dd/mm/yyyy"
f2 = "dd/mm/yy"
@test Dates.Date("28/05/2014", f) == Dates.Date(2014, 5, 28)
@test Dates.Date("28/05/14", f2) + Dates.Year(2000) == Dates.Date(2014, 5, 28)

# Customizing locale
Dates.LOCALES["french"] = Dates.DateLocale(
    ["janvier", "février", "mars", "avril", "mai", "juin",
     "juillet", "août", "septembre", "octobre", "novembre", "décembre"],
    ["janv", "févr", "mars", "avril", "mai", "juin",
     "juil", "août", "sept", "oct", "nov", "déc"],
    ["lundi", "mardi", "mercredi", "jeudi", "vendredi", "samedi", "dimanche"],
    [""],
)

f = "dd uuuuu yyyy"
@test Dates.Date("28 mai 2014", f; locale="french") == Dates.Date(2014, 5, 28)
@test Dates.format(Dates.Date(2014, 5, 28), f; locale="french") == "28 mai 2014"
@test Dates.Date("28 févr 2014", f; locale="french") == Dates.Date(2014, 2, 28)
@test Dates.format(Dates.Date(2014, 2, 28), f; locale="french") == "28 févr 2014"
@test Dates.Date("28 août 2014", f; locale="french") == Dates.Date(2014, 8, 28)
@test Dates.format(Dates.Date(2014, 8, 28), f; locale="french") == "28 août 2014"
@test Dates.Date("28 avril 2014", f; locale="french") == Dates.Date(2014, 4, 28)
@test Dates.format(Dates.Date(2014, 4, 28), f; locale="french") == "28 avril 2014"

f = "dd u yyyy"
@test Dates.Date("28 avril 2014", f; locale="french") == Dates.Date(2014, 4, 28)
f = "dduuuuyyyy"
# parses 3 and 4 character month names
@test Dates.Date("28mai2014", f; locale="french") == Dates.Date(2014, 5, 28)
@test Dates.Date("28août2014", f; locale="french") == Dates.Date(2014, 8, 28)
# doesn't parse month name greater than 4 chars
@test_throws ArgumentError Dates.Date("28avril2014", f; locale="french")

# From Tony Fong
f = "dduuuyy"
@test Dates.Date("01Dec09", f) + Dates.Year(2000) == Dates.Date(2009, 12, 1)
@test Dates.format(Dates.Date(2009, 12, 1), f) == "01Dec09"
f = "dduuuyyyy"
@test Dates.Date("01Dec2009", f) == Dates.Date(2009, 12, 1)
@test Dates.format(Dates.Date(2009, 12, 1), f) == "01Dec2009"
f = "duy"
const globex = ["f", "g", "h", "j", "k", "m", "n", "q", "u", "v", "x", "z"]
locale = Dates.DateLocale(globex, map(uppercase, globex), globex[1:7], globex[1:7])
@test Dates.Date("1F4", f; locale=locale) + Dates.Year(2010) == Dates.Date(2014, 1, 1)
@test Dates.format(Dates.Date(2014, 1, 1), f; locale=locale) == "1F4"

# From Matt Bauman
f = "yyyy-mm-ddTHH:MM:SS"
@test Dates.DateTime("2014-05-28T16:46:04", f) == Dates.DateTime(2014, 5, 28, 16, 46, 04)

# Try to break stuff

# Specified mm/dd, but date string has day/mm
@test_throws ArgumentError Dates.DateTime("18/05/2009", "mm/dd/yyyy")
@test_throws ArgumentError Dates.DateTime("18/05/2009 16", "mm/dd/yyyy hh")
# Used "mm" for months AND minutes
@test_throws ArgumentError Dates.DateTime("18/05/2009 16:12", "mm/dd/yyyy hh:mm")
# Date string has different delimiters than format string
@test_throws ArgumentError Dates.DateTime("18:05:2009", "mm/dd/yyyy")

f = "y m d"
@test Dates.Date("1 1 1", f) == Dates.Date(1)
@test Dates.Date("10000000000 1 1", f) == Dates.Date(10000000000)
@test_throws ArgumentError Dates.Date("1 13 1", f)
@test_throws ArgumentError Dates.Date("1 1 32", f)
@test_throws ArgumentError Dates.Date(" 1 1 32", f)
@test_throws ArgumentError Dates.Date("# 1 1 32", f)
@test Dates.Date("1", f) == Dates.Date(1)
@test Dates.Date("1 ", f) == Dates.Date(1)
@test Dates.Date("1 2", f) == Dates.Date(1, 2)
# can't find space delimiter (finds '/'), so fails
@test_throws ArgumentError Dates.Date("2000/1", f)

f = "ymd"
@test Dates.Date("111", f) == Dates.Date(1)
@test Dates.Date("1", f) == Dates.Date(1)

@test Dates.DateTime("20140529 120000", "yyyymmdd HHMMSS") == Dates.DateTime(2014, 5, 29, 12)

@test Dates.Date(string(Dates.Date(dt))) == Dates.Date(dt)
@test Dates.DateTime(string(dt)) == dt

# Vectorized
dr = ["2000-01-01", "2000-01-02", "2000-01-03", "2000-01-04", "2000-01-05",
      "2000-01-06", "2000-01-07", "2000-01-08", "2000-01-09", "2000-01-10"]
dr2 = [Dates.Date(2000) : Dates.Date(2000, 1, 10);]
@test Dates.Date(dr) == dr2
@test Dates.Date(dr, "yyyy-mm-dd") == dr2
@test Dates.DateTime.(dr) == Dates.DateTime.(dr2)
@test Dates.DateTime(dr, "yyyy-mm-dd") == Dates.DateTime.(dr2)

@test Dates.format(dr2) == dr
@test Dates.format(dr2, "yyyy-mm-dd") == dr

@test typeof(Dates.Date(dr)) == Array{Date, 1}

# Issue 13
t = Dates.DateTime(1, 1, 1, 14, 51, 0, 118)
@test Dates.DateTime("[14:51:00.118]", "[HH:MM:SS.sss]") == t
@test Dates.DateTime("14:51:00.118", "HH:MM:SS.sss") == t
@test Dates.DateTime("[14:51:00.118?", "[HH:MM:SS.sss?") == t
@test Dates.DateTime("?14:51:00.118?", "?HH:MM:SS.sss?") == t
@test Dates.DateTime("x14:51:00.118", "xHH:MM:SS.sss") == t
@test Dates.DateTime("14:51:00.118]", "HH:MM:SS.sss]") == t

# RFC1123Format
dt = Dates.DateTime(2014, 8, 23, 17, 22, 15)
@test Dates.format(dt, Dates.RFC1123Format) == "Sat, 23 Aug 2014 17:22:15"
@test Dates.DateTime(Dates.format(dt, Dates.RFC1123Format), Dates.RFC1123Format) == dt
@test Dates.format(dt, "yyyy-mm-ddTHH:MM:SS E") == "2014-08-23T17:22:15 Saturday"
@test Dates.format(dt, "yyyy-mm-ddTHH:MM:SS e") == "2014-08-23T17:22:15 Sat"
@test Dates.format(dt, "yyyy-mm-dd E") == "2014-08-23 Saturday"
@test Dates.format(dt, "yyyy-mm-dd e") == "2014-08-23 Sat"
@test Dates.format(dt, "yyyy-e-mm-dd") == "2014-Sat-08-23"

@test Dates.format(Dates.DateTime(2014, 1, 2, 0, 0, 0, 999), Dates.RFC1123Format) == "Thu, 02 Jan 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 2, 18, 0, 0, 0, 9), Dates.RFC1123Format) == "Tue, 18 Feb 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 3, 8, 0, 0, 0, 9), Dates.RFC1123Format) == "Sat, 08 Mar 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 4, 28, 0, 0, 0, 9), Dates.RFC1123Format) == "Mon, 28 Apr 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 5, 10, 0, 0, 0, 9), Dates.RFC1123Format) == "Sat, 10 May 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 6, 4, 0, 0, 0, 9), Dates.RFC1123Format) == "Wed, 04 Jun 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 7, 13, 0, 0, 0, 9), Dates.RFC1123Format) == "Sun, 13 Jul 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 8, 17, 0, 0, 0, 9), Dates.RFC1123Format) == "Sun, 17 Aug 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 9, 20, 0, 0, 0, 9), Dates.RFC1123Format) == "Sat, 20 Sep 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 10, 31, 0, 0, 0, 9), Dates.RFC1123Format) == "Fri, 31 Oct 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 11, 2, 0, 0, 0, 9), Dates.RFC1123Format) == "Sun, 02 Nov 2014 00:00:00"
@test Dates.format(Dates.DateTime(2014, 12, 5, 0, 0, 0, 9), Dates.RFC1123Format) == "Fri, 05 Dec 2014 00:00:00"

dt = Dates.DateTime(2016, 11, 12, 7, 45, 36)
@test parse(Dates.DateTime, "Sat, 12 Nov 2016 07:45:36", Dates.RFC1123Format) == dt
@test parse(Dates.DateTime, "Mon, 12 Nov 2016 07:45:36", Dates.RFC1123Format) == dt  # Wrong day of week
@test_throws ArgumentError parse(Date, "Foo, 12 Nov 2016 07:45:36", Dates.RFC1123Format)

# Issue 15195
let f = "YY"
    @test Dates.format(Dates.Date(1999), f) == "1999"
    @test Dates.format(Dates.Date(9), f) == "09"
    @test Dates.format(typemax(Dates.Date), f) == "252522163911149"
end

# Issue: https://github.com/quinnj/TimeZones.jl/issues/19
let
    const Zulu = String

    function Dates.tryparsenext(d::Dates.DatePart{'Z'}, str, i, len)
        Dates.tryparsenext_word(str, i, len, Dates.min_width(d), Dates.max_width(d))
    end

    str = "2015-07-24T05:38:19.591Z"
    dt = Dates.DateTime(2015, 7, 24, 5, 38, 19, 591)
    parsed = Any[
        Dates.Year(2015), Dates.Month(7), Dates.Day(24),
        Dates.Hour(5), Dates.Minute(38), Dates.Second(19), Dates.Millisecond(591)
    ]

    format = "yyyy-mm-ddTHH:MM:SS.sssZ"
    escaped_format = "yyyy-mm-dd\\THH:MM:SS.sss\\Z"

    # Typically 'Z' isn't treated as a specifier so it doesn't have to be escaped
    @test Dates.parse_components(str, Dates.DateFormat(format)) == parsed
    @test Dates.parse_components(str, Dates.DateFormat(escaped_format)) == parsed

    try
        # Make 'Z' into a specifier
        Dates.CONVERSION_SPECIFIERS['Z'] = Zulu
        Dates.CONVERSION_DEFAULTS[Zulu] = ""

        @test Dates.parse_components(str, Dates.DateFormat(format)) == [parsed; Zulu("Z")]
        @test Dates.parse_components(str, Dates.DateFormat(escaped_format)) == parsed
    finally
        delete!(Dates.CONVERSION_SPECIFIERS, 'Z')
        delete!(Dates.CONVERSION_DEFAULTS, Zulu)
    end

    # Ensure that the default behaviour has been restored
    @test Dates.parse_components(str, Dates.DateFormat(format)) == parsed
    @test Dates.parse_components(str, Dates.DateFormat(escaped_format)) == parsed
end

# Issue 10817
@test Dates.Date("Apr 01 2014", "uuu dd yyyy") == Dates.Date(2014, 4, 1)
@test_throws ArgumentError Dates.Date("Apr 01 xx 2014", "uuu dd zz yyyy")
@test_throws ArgumentError Dates.Date("Apr 01 xx 2014", "uuu dd    yyyy")

# Issue 21001
for (ms, str) in zip([0, 1, 20, 300, 450, 678], ["0", "001", "02", "3", "45", "678"])
    dt = DateTime(2000, 1, 1, 0, 0, 0, ms)
    @test Dates.format(dt, "s") == str
    @test Dates.format(dt, "ss") == rpad(str, 2, '0')
    @test Dates.format(dt, "sss") == rpad(str, 3, '0')
    @test Dates.format(dt, "ssss") == rpad(str, 4, '0')
end

# Issue #21504
@test isnull(tryparse(Dates.Date, "0-1000"))

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia