text:
module KernelDensityEstimatePlotting

using LinearAlgebra, Statistics
using KernelDensityEstimate
using Gadfly, Colors, Cairo, Fontconfig, Compose

import Gadfly: plot
# KDE imports will be remove when plotting is permanently removed from the KDE package
# import KernelDensityEstimate: plotKDE,plotKDEContour,drawAllPairs,stackMarginals,vstackedPlots,drawHorDens,toggleYTicks

export
  # Gadfly plotting functions
  plot,
  getColorsByLength, # duplicated in RoMEPlotting
  plotKDEContour,
  plotKDE,
  drawPair,
  drawAllPairs,
  stackMarginals,
  vstackedPlots,
  drawHorDens,
  toggleYTicks,
  DOYTICKS

NothingUnion{T}   = Union{Nothing, T}

# Guide.manual_color_key("Legend", ["Points", "Line"], ["green", "deepskyblue"])

DOYTICKS=true

function toggleYTicks()
  global DOYTICKS
  DOYTICKS = DOYTICKS ? false : true
  return DOYTICKS
end

# function plot(p::BallTreeDensity)
#    plotKDE(p)
# end

function plot(darr::Union{BallTreeDensity, Vector{BallTreeDensity}};
      c::NothingUnion{Vector}=nothing,
      N::Int=200,
      rmax=-Inf,rmin=Inf,  # should be deprecated
      axis::NothingUnion{Array{Float64,2}}=nothing,
      dims::NothingUnion{VectorRange{Int}}=nothing,
      xlbl::T="X", # to be deprecated
      title::NothingUnion{T}=nothing,
      legend::NothingUnion{Vector{T}}=nothing,
      dimLbls::NothingUnion{Vector{T}}=nothing,
      levels::NothingUnion{Int}=nothing,
      fill=false, layers::Bool=false ) where {T <: AbstractString}
  #
  plotKDE(darr, c=c, N=N, rmin=rmin, rmax=rmax, axis=axis, dims=dims, xlbl=xlbl, title=title, legend=legend, dimLbls=dimLbls, levels=levels, fill=fill, layers=layers )
end

function draw1D!(bd::BallTreeDensity,
                 bins::AbstractArray{Float64,1},
                 e,
                 c::T="deepskyblue",
                 myStyle::T="";
                 xlbl="X",
                 points::Bool=true,
                 legend=nothing,
                 title::NothingUnion{T}=nothing,
                 fill=false, layers::Bool=false ) where {T <: AbstractString}
  #
  global DOYTICKS

  yV = evaluateDualTree(bd,bins)
  # clamp max y values
  yV[3.0 .< yV] .= 3.0
  if e == nothing
    ptArr = Any[]

    if points
      tps = getPoints(bd)
      lrug = Gadfly.layer(x=tps, y=-0.01.*ones(size(tps,2)), Geom.point)
      push!(ptArr, lrug)
    end
    l1 = Gadfly.layer(x=bins,y=yV,Geom.line, Gadfly.Theme(default_color=parse(Colorant,c),line_width=2pt, point_size=1pt))
    push!(ptArr, l1)
    push!(ptArr, Guide.xlabel(xlbl))
    if !DOYTICKS
      # e=Gadfly.plot(x=bins,y=yV,Geom.line, Gadfly.Theme(default_color=parse(Colorant,c)),Guide.xlabel(xlbl))
    # else
      push!(ptArr, Guide.ylabel(""))
      push!(ptArr,Guide.yticks(ticks=nothing))
      # e=Gadfly.plot(x=bins,y=yV,Geom.line, Gadfly.Theme(default_color=parse(Colorant,c)),Guide.xlabel(xlbl),Guide.ylabel(""),Guide.yticks(ticks=nothing))
    end
    if legend != nothing
      push!(ptArr, legend)
    end
    if title != nothing
      push!(ptArr, Guide.title(title))
    end

    e = !layers ? Gadfly.plot(ptArr...) : ptArr
  else
    push!(e.layers, layer(x=bins, y=yV, Geom.line, Gadfly.Theme(default_color=parse(Colorant,c),line_width=2pt))[1])
  end

  #mx = maximum(y)
  #wts = getWeights(bd)
  return e
end

function plotKDEContour(pp::Vector{BallTreeDensity};
    xmin=-Inf,xmax=Inf,ymin=-Inf,ymax=Inf,
    xlbl::T="x", ylbl::T="y",
    N::Int=200,
    c::NothingUnion{Vector}=nothing,
    legend=nothing,
    title::NothingUnion{T}=nothing,
    levels::NothingUnion{Int}=nothing,
    fill=false, layers::Bool=false ) where {T <: AbstractString}

  rangeV = getKDERange(pp[1])
  size(rangeV,1) == 2 ? nothing : error("plotKDEContour must receive two dimensional kde, you gave $(Ndim(x))")
  xmin = xmin != -Inf ? xmin : rangeV[1,1]
  xmax = xmax != Inf ? xmax : rangeV[1,2]
  ymin = ymin != -Inf ? ymin : rangeV[2,1]
  ymax = ymax != Inf ? ymax : rangeV[2,2]
  for i in 2:length(pp)
    rangeV = getKDERange(pp[i])
    xmin = xmin <= rangeV[1,1] ? xmin : rangeV[1,1]
    xmax = xmax >= rangeV[1,2] ? xmax : rangeV[1,2]
    ymin = ymin <= rangeV[2,1] ? ymin : rangeV[2,1]
    ymax = ymax >= rangeV[2,2] ? ymax : rangeV[2,2]
  end

  PL = []

  # default options
  CO = levels == nothing ? Geom.contour : Geom.contour(levels=levels)
  if c == nothing
    c = ["deepskyblue" for i in 1:length(pp)]
  else
    push!(PL, Gadfly.Scale.color_none)
  end
  # Gadfly.plot(

  i = 0
  for p in pp
    i+=1
    push!(PL, layer(z=(x,y)->evaluateDualTree(p,([[x]';[y]']))[1],
    x=range(xmin,stop=xmax,length=N),
    y=range(ymin,stop=ymax,length=N),
    CO,
    Theme(default_color=parse(Colorant,c[i]),line_width=2.25pt))[1] )
    ## trying to get rug to work
    # rugpl = plot(x=getPoints(p)[1,:], y=getPoints(p)[2,:], Guide.xrug, Guide.yrug)
    # push!( PL, rugpl.layers[1] )
  end

  push!(PL,Coord.Cartesian(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax))
  push!(PL,Guide.xlabel(xlbl), Guide.ylabel(ylbl))

  # Might be a per layer theme
  !fill ? nothing : push!(PL, Theme(background_color=colorant"white"))

  if legend != nothing
    push!(PL, legend)
  end

  if title != nothing
    push!(PL, Guide.title(title))
  end
  if !layers
    return Gadfly.plot(PL...)
  else
    return PL
  end
end
# p1 = plot(
#        layer(z=ff,x=linspace(-5,5,100),y=linspace(-5,5,100), Geom.contour(levels=5),Theme(default_color=parse(Colorant,"blue"))),
#        layer(z=ff2,x=linspace(-5,5,100),y=linspace(-5,5,100), Geom.contour(levels=5),Theme(default_color=parse(Colorant,"red"))),
#        Scale.color_none)

function plotKDEContour(p::BallTreeDensity;
    xmin=-Inf,xmax=Inf,ymin=-Inf,ymax=Inf,
    xlbl::T="x", ylbl::T="y",
    N::Int=200,
    c::NothingUnion{Vector}=nothing,
    legend=nothing,
    title::NothingUnion{T}=nothing,
    levels::NothingUnion{Int}=nothing,
    fill=false, layers::Bool=false ) where {T <: AbstractString}

    plotKDEContour([p],
      xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,
      xlbl=xlbl, ylbl=ylbl,
      N=N,
      c=c,
      legend=legend,
      title=title,
      levels=levels,
      fill=fill, layers=layers )
end

function drawPair(xx::Vector{BallTreeDensity},
                  dims::Vector{Int};
                  axis::NothingUnion{Array{Float64,2}}=nothing,
                  dimLbls::NothingUnion{Vector{T}}=nothing,
                  legend=nothing,
                  title::NothingUnion{T}=nothing,
                  levels::NothingUnion{Int}=nothing,
                  c::NothingUnion{Vector}=nothing,
                  fill=false,
                  layers::Bool=false,
                  overlay=nothing ) where {T <: AbstractString}
  #
  # pts = getPoints(x);
  xmin, xmax, ymin, ymax = -Inf,Inf,-Inf,Inf
  if axis != nothing
    xmin, xmax, ymin, ymax = axis[dims[1],1], axis[dims[1],2], axis[dims[2],1], axis[dims[2],2]
  end

  xlbl, ylbl = nothing, nothing
  if dimLbls!=nothing
    xlbl, ylbl = dimLbls[dims[1]], dimLbls[dims[2]]
  end

  X = BallTreeDensity[]
  for x in xx
    push!(X, marginal(x,dims))
  end
  plr = plotKDEContour(X,
    xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,
    xlbl=xlbl,ylbl=ylbl,
    legend=legend,
    title=title,
    levels=levels,c=c,
    fill=fill, layers=layers  )

  # add overlay
  if overlay != nothing
    union!(plr.layers, overlay.layers)
  end
  plr
end

function stacking(spp::Vector{<:Compose.Context})
  hstack(spp...)
end

# function to draw all pairs of mulitdimensional kernel density estimate
# axis is matrix with rows as dimensions and two columns for min and max axis cutoffs
function drawAllPairs(xx::Vector{BallTreeDensity};
      dims::NothingUnion{VectorRange{Int}}=nothing,
      axis::NothingUnion{Array{Float64,2}}=nothing,
      dimLbls::NothingUnion{Vector{T}}=nothing,
      legend=nothing,
      title::NothingUnion{T}=nothing,
      levels::NothingUnion{Int}=nothing,
      c::NothingUnion{Vector}=nothing,
      fill=false,
      layers::Bool=false,
      overlay=nothing  ) where {T <: AbstractString}

  # pts = getPoints(xx[1]);
  # e = [];
  dims = dims != nothing ? collect(dims) : collect(1:Ndim(xx[1]))
  Nout = length(dims);
  PlotI2 = triu(repeat( (dims)' ,Nout, 1), 1);
  PlotI1 = triu(repeat((dims)' , Nout, 1)', 1);
  PlotI1, PlotI2 = PlotI1[(LinearIndices(PlotI1))[findall(x->x!=0,PlotI1)]],  PlotI2[(LinearIndices(PlotI2))[findall(x->x!=0,PlotI2)]];
  # PlotI1, PlotI2 = PlotI1[find(PlotI1)],  PlotI2[find(PlotI2)];
  Ncol = round(Int, sqrt(length(PlotI2)));
  Nrow = ceil(Int, length(PlotI2)/Ncol);

  subplots = Array{Gadfly.Plot,2}(undef, Nrow,Ncol)
  for iT=1:length(PlotI2)
    # only returns layers for first pair
    ovl = Nout==2 ? overlay : nothing
    if !layers
      subplots[iT] = drawPair(xx,[PlotI1[iT];PlotI2[iT]], axis=axis, dimLbls=dimLbls, legend=legend, title=title, levels=levels, c=c, fill=fill, layers=layers, overlay=ovl );
    else
      return drawPair(xx,[PlotI1[iT];PlotI2[iT]], axis=axis, dimLbls=dimLbls, legend=legend, title=title, levels=levels, c=c, fill=fill, layers=layers, overlay=ovl )
    end
  end;

  Nrow==1 && Ncol==1 ? nothing : println("Multiple planes stacked into Compose.Context, use Gadfly.draw(PNG(file.png,10cm,10cm),plothdl). Or PDF.")
  hh = Vector{Gadfly.Context}(undef, Nrow)
  for i in 1:Nrow
    sp = Compose.Context[]
    for j in 1:Ncol
      if Nrow == 1 && Ncol==1

      else
        try
          push!(sp,hstack(subplots[i,j])) # very hacky line, but Compose.Context and Gadfly.Plot not playing nice together in hstack
        catch e
          print("KDEPlotting01.jl/drawAllPairs -- supressing all exceptions for stacking empty contour plots")
          println(e)
          push!(sp,Gadfly.context())
        end
      end
    end
    hh[i] = stacking(sp) #hstack(sp) #subplots[i,:])
  end

  return Nrow==1 && Ncol==1 ? subplots[1,1] : vstack(hh...)
end

# """
#     $(SIGNATURES)
#
# Standardize the length colors used by RoMEPlotting.
# """
function getColorsByLength(len::Int=7)::Vector{String}
  len > 99 ? error("Don't have enough colors, 100 is the max.") : nothing
  COLORS = String["red";"green";"blue";"magenta";"yellow";"deepskyblue"]
  if len > 6
    scale = len -7 + 2 # + 2 is artificial and avoids gray100==white
    scale = 100/scale
    for i in 7:len
      push!(COLORS, "gray$(floor(Int,(i-7)*scale))")
    end
  end
  return COLORS[1:len]
end

# function to draw all pairs of mulitdimensional kernel density estimate
# axis is matrix with rows as dimensions and two columns for min and max axis cutoffs
function plotKDE(darr::Array{BallTreeDensity,1};
                 c::NothingUnion{Vector}=getColorsByLength(length(darr)), # nothing
                 N::Int=200,
                 rmax=-Inf,rmin=Inf,  # should be deprecated
                 axis::NothingUnion{Array{Float64,2}}=nothing,
                 dims::NothingUnion{VectorRange{Int}}=nothing,
                 xlbl::T="X", # to be deprecated
                 title::NothingUnion{T}=nothing,
                 legend::NothingUnion{Vector{T}}=nothing,
                 dimLbls::NothingUnion{Vector{T}}=nothing,
                 levels::NothingUnion{Int}=nothing,
                 fill=false,
                 points::Bool=true,
                 layers::Bool=false,
                 overlay=nothing ) where {T <: AbstractString}
    #
    # defaults
    defaultcolor = false
    if c==nothing || length(c) != length(darr)
      c = getColorsByLength(length(darr))
      defaultcolor = true
    end
    # c = (length(c)>=2) ? c : repeat(c,length(darr))
    lg = (legend == nothing) ? nothing : Guide.manual_color_key("Legend", legend, c)

    H = nothing
    i = 0

    Ndims = Ndim(darr[1])
    dim = dims!=nothing ? dims : 1:Ndims #.bt.dims
    dimLbls = dimLbls!=nothing ? dimLbls : String["$(i)" for i in 1:Ndims]
    dim = collect(dim)
    if length(dim) == 1
      for bd in darr
          i+=1
          mbd = marginal(bd,dim)
          rangeV = getKDERange(mbd)
          if (length(dim) == 1) # Ndim(bd)
            if rangeV[1] > rmin  rangeV[1] = rmin end
            if rmax > rangeV[2]  rangeV[2] = rmax end
            if axis!=nothing
              di = dim[1]
              if rangeV[1] > axis[di,1]  rangeV[1] = axis[di,1] end
              if axis[di,2] > rangeV[2]  rangeV[2] = axis[di,2] end
            end
            H=draw1D!(mbd,range(rangeV[1],stop=rangeV[2],length=N), H, c[i],xlbl=xlbl,legend=lg, title=title,            points=points, fill=fill) #,argsPlot,argsKDE
          else
            #
          end
      end

    else
      color = defaultcolor ? nothing : c
      H = drawAllPairs(darr, axis=axis, dims=dim, dimLbls=dimLbls, legend=lg, title=title, levels=levels, c=color, fill=fill, layers=layers)
    end

    if overlay != nothing && length(dim) <= 2
      union!(H.layers, )
    end

    return H
end

function plotKDE(bd::BallTreeDensity;
                 c::NothingUnion{Vector}=nothing,
                 N::Int=200,
                 rmax=-Inf,rmin=Inf,  # should be deprecated
                 axis::NothingUnion{Array{Float64,2}}=nothing,
                 dims::NothingUnion{VectorRange{Int}}=nothing,
                 xlbl::T="X",
                 legend::NothingUnion{Vector{T}}=nothing,
                 title::NothingUnion{T}=nothing,
                 dimLbls::NothingUnion{Vector{T}}=nothing,
                 levels::NothingUnion{Int}=nothing,
                 fill=false,
                 points::Bool=true,
                 layers::Bool=false,
                 overlay=nothing  ) where {T <: AbstractString}

  plotKDE([bd],N=N,c=c,rmax=rmax,rmin=rmin,xlbl=xlbl,legend=legend, dims=dims, axis=axis, dimLbls=dimLbls, levels=levels, title=title, fill=fill, layers=layers, overlay=overlay, points=points )
end

macro ifund(exp)
    local e = :($exp)
    isdefined(e.args[1]) ? :($(e.args[1])) : :($(esc(exp)))
end

function drawHorDens(pDens::Array{BallTreeDensity,1}; N::Int=200,
                    gt=Union{}, lbls=String[], extend::Float64=0.1)
    len = length(pDens)
    h = Array{Gadfly.Plot,1}(undef, len)
    # r = Array{Future,1}(len) #Array{RemoteRef,1}(len)

    if gt!=Union{}
      for i in 1:length(pDens)
        g = kde!(gt[i][1,:],[gt[i][2,1]])
        v = getKDERange(pDens[i],extend=extend)
        p = plotKDE([pDens[i];g],c=["black";"red"],rmin=v[1],rmax=v[2],xlbl=lbls[i])
        h[i] =p;
      end
    else
      for i in 1:length(pDens)
         h[i]=plotKDE([pDens[i]],xlbl=lbls[i])
      end
    end
    #println("draw $(len) plots")
    # [r[i] = @spawn plotKDE(pDens[i], N=N) for i in 1:len]
    # [h[i] = fetch(r[i]) for i in 1:len]

    hstack(h)
end

function stackMarginals(P::Array{BallTreeDensity,1}, m::Int64)
  ret = Array{BallTreeDensity,1}(length(P))
  for i in 1:length(P)
    ret[i] =  marginal(P[i],[m])#push!()
  end
  return ret
end

function vstackedPlots(plots::Vector{<: Gadfly.Compose.Context})
    evalstr = ""
    for i in 1:length(plots)
        evalstr = string(evalstr, ",plots[$(i)]")
    end
    evalstr = string("vstack(",evalstr[2:end],")")
    eval(parse(evalstr))
end

end  # module


text:
using Invest
using Test

@testset "Invest" begin
    @testset "Option" begin
        include("option_test.jl")
    end
end

ext: jl


text:
const threshold = 1e-8

myeqzero{T<:Real}(x::T) = x == zero(T)
myeqzero{T<:AbstractFloat}(x::T) = -threshold < x < threshold
myeqzero{T<:Real}(x::Vector{T}) = myeqzero(sum(abs(x)))

myeq{T<:Real}(x::T, y::T) = x == y
# I check with zero because isapprox(0, 1e-100) is false...
# but isapprox(1e-100, 2e-100) should be false
myeq{T<:AbstractFloat}(x::T, y::T) = (x == zero(T) ? myeqzero(y) : (y == zero(T) ? myeqzero(x) : isapprox(x, y)))
myeq{T<:Real, S<:Real}(x::T, y::S) = myeq(promote(x, y)...)

mylt{T<:Real}(x::T, y::T) = x < y
mylt{T<:AbstractFloat}(x::T, y::T) = x < y && !myeq(x, y)
mylt{S<:Real,T<:Real}(x::S, y::T) = mylt(promote(x, y)...)
mygt{S<:Real, T<:Real}(x::S, y::T) = mylt(y, x)
myleq{T<:Real}(x::T, y::T) = x <= y
myleq{T<:AbstractFloat}(x::T, y::T) = x <= y || myeq(x, y)
myleq{S<:Real,T<:Real}(x::S, y::T) = myleq(promote(x, y)...)
mygeq{T<:Real}(x::T, y::T) = myleq(y, x)
mypos{T<:Real}(x::T) = mygt(x, zero(T))
myneg{T<:Real}(x::T) = mylt(x, zero(T))
mynonneg{T<:Real}(x::T) = mygeq(x, zero(T))
mynonpos{T<:Real}(x::T) = myleq(x, zero(T))

function mypromote_type{T1,T2}(::Type{T1}, ::Type{T2})
  if T1 <: AbstractFloat || T2 <: AbstractFloat
    Float64
  else
    Rational{BigInt}
  end
end


text:
using Weave
using Test

function mmtest(source, resfile, doctype)
  weave("documents/$source", out_path = "documents/multimedia/$resfile", mod=:sandbox,
    doctype=doctype, template = "templates/mini.tpl")
  result =  read("documents/multimedia/$resfile", String)
  ref =  read("documents/multimedia/$resfile.ref", String)
  @test result == ref
  rm("documents/multimedia/$resfile")
end

mmtest("rich_output.jmd", "rich_output.html", "md2html")
mmtest("rich_output.jmd", "rich_output.md", "pandoc")
mmtest("rich_output.jmd", "rich_output.tex", "tex")
mmtest("rich_output.jmd", "rich_output.github", "github")


text:
using MonteCarloMarkovKernels
using StaticArrays
import Statistics.mean
using Test

zero2 = [0.0, 0.0]
μ1 = [-1.0, 0.0]
Σ1 = [1.0 0.25 ; 0.25 3.0]

Szero2 = SVector{2, Float64}(zero2)

function makelogMVN(μ::SVector{d, Float64}, Σ::SMatrix{d, d, Float64}) where d
  invΣ = inv(Σ)
  lognc = - 0.5 * d * log(2 * π) - 0.5 * logdet(Σ)
  function lpi(x::SVector{d, Float64})
    v = x - μ
    return lognc - 0.5*dot(v, invΣ * v)
  end
  return lpi
end

logtarget = makelogMVN(SVector{2, Float64}(μ1), SMatrix{2, 2, Float64}(Σ1))

seed!(12345)

P_AM = makeAMKernel(logtarget, 2)
chain = simulateChain(P_AM, Szero2, 2^20)

@test mean(chain) ≈ μ1 atol = 0.01
@test MonteCarloMarkovKernels.cov(chain) ≈ Σ1 atol = 0.05

vs = (i->(x->x[i]).(chain)).(1:2)

d1(x) = 1/sqrt(2*π*Σ1[1,1])*exp(-1/2/Σ1[1,1]*(x-μ1[1])^2)
d2(x) = 1/sqrt(2*π*Σ1[2,2])*exp(-1/2/Σ1[2,2]*(x-μ1[2])^2)

xs, ys = kde(vs[1], P_AM(:acceptanceRate))
@test ys ≈ d1.(xs) atol=0.05

xs, ys = kde(vs[2], P_AM(:acceptanceRate))
@test ys ≈ d2.(xs) atol=0.05

xs, ys, f1 = kde(vs[1], vs[2], P_AM(:acceptanceRate))
M1 = broadcast(f1, xs, ys')
M2 = broadcast((x,y)->exp(logtarget((SVector{2,Float64}(x,y)))), xs, ys')
@test M1 ≈ M2 atol=0.05


text:
function idx_con()
    CON_LABEL = 1
    CON_TYPE = 2
    CON_FR = 3
    CON_TO = 4
    CON_ID = 5

    return CON_LABEL, CON_TYPE, CON_FR, CON_TO, CON_ID
end

function parse_con(filename::AbstractString; delim=" \t")
    f = open(filename, "r")
    lines = readlines(f)
    close(f)

    pos = 1
    ndata = 0
    while pos <= length(lines)
        line = strip(lines[pos])

        if line[1] == 'R' || line[1] == 'O'
            ndata += 1
        end

        pos += 1
    end

    data = Array{Any}(undef, ndata, 5)

    label = ""
    ndata = 0
    pos = 1
    while pos <= length(lines)
        line = strip(lines[pos])
        fields = split(line)

        if line[1] == 'C'      # CONTINGENCY
            label = fields[2]
        elseif line[1] == 'R'  # REMOVE
            ndata += 1
            data[ndata,1] = label
            data[ndata,2] = 2                     # ctg type: generator
            data[ndata,3] = parse(Int, fields[6]) # I
            data[ndata,4] = 0                     # No J for generator
            data[ndata,5] = fields[3]             # ID
        elseif line[1] == 'O'  # OPEN
            ndata += 1
            data[ndata,1] = label
            data[ndata,2] = 1                     # ctg type: branch
            data[ndata,3] = parse(Int, fields[5]) # I
            data[ndata,4] = parse(Int, fields[8]) # J
            data[ndata,5] = fields[10]            # CKT
        end

        pos += 1
    end

    condata = Dict{String,Array}()
    condata["CONTINGENCY"] = data

    return condata
end


text:
function nlp_cvx_108_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - convex objective
    # - intersection of constraints
    # - convex constraint, only in positive domain
    # Variants
    #   010 - no binding constraints
    #   011 - intersection point
    #   012 - intersection point
    #   013 - one binding constraint

    model = Model(optimizer)

    @variable(model, x >= 0)
    @variable(model, y >= 0)

    @objective(model, Min, (x - 1.0)^2 + (y - 0.75)^2)
    @NLconstraint(model, 2 * x^2 - 4x * y - 4 * x + 4 <= y)
    @NLconstraint(model, y^2 <= -x + 2)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 0, tol = objective_tol)
    return check_solution([x, y], [1, 0.75], tol = primal_tol)
end


text:
import Base.indexin

function indexin(dd::Dict{String,Any},species_symbol::String;
    key="total_species_list")::Union{Nothing,Int}

    # get the total species list -
    if (haskey(dd,key) == false)
        return nothing
    end
    total_species_list = dd[key]

    # check -
    return findfirst(x->x==species_symbol,total_species_list)
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
#using ScenTreesMakie
#using LinearAlgebra

# Benchmarks different versions of the tree structure Tree() with type annotations or without
# and some constructor optimization

# Test branching structure
#bs = [1,2,2,2,2,2,2,2,2,2,2,2,2,2,2];  n = length(bs)

# Current ScenTreesMakie version, with parameterized types
#    trr = ScenTreesMakie.Tree(bs);
#    trr2 = trr;
#    path3() = ScenTreesMakie.gaussian_path(n)

    # benchmarking approxiation for ScenTreesMakie version with type annotation

 #   @time ScenTreesMakie.tree_approximation!(trr,path3,10,2,2);

    # Alois Version of tree_approximation
#    @time tree_approximation2!(trr2,path3,10,2,2);


text:
export SHA1, PartialSHA1, Oid, raw, iszero, @sha1_str

const OID_RAWSZ = 20
const OID_HEXSZ = OID_RAWSZ * 2
const OID_MINPREFIXLEN = 4

# immutable Oid
#   id1::UInt8
#   id2::UInt8
#   ...
#   id20::UInt8
# end

@eval begin
    $(Expr(:type, false, :Oid,
        Expr(:block,
            [Expr(:(::), symbol("id$i"), :UInt8) for i=1:OID_RAWSZ]...)))
end

# default Oid constructor (all zeros)
@eval begin
    Oid() = $(Expr(:call, :Oid, [:(0x00) for _=1:OID_RAWSZ]...))
end

Oid(ptr::Ptr{Oid}) = unsafe_load(ptr)::Oid

Oid(ptr::Ptr{UInt8}) = begin
    if ptr == C_NULL
        throw(ArgumentError("NULL pointer passed to Oid() constructor"))
    end
    oid_ptr = Oid[Oid()]
    ccall((:git_oid_fromraw, libgit2), Void, (Ptr{Oid}, Ptr{UInt8}), oid_ptr, ptr)
    return oid_ptr[1]
end

Oid(id::Array{UInt8,1}) = begin
    if length(id) != OID_RAWSZ
        throw(ArgumentError("invalid raw buffer size"))
    end
    return Oid(pointer(id))
end

Oid(id::AbstractString) = begin
    bstr = bytestring(id)
    if sizeof(bstr) != OID_HEXSZ
        throw(ArgumentError("invalid hex size"))
    end
    oid_ptr = Oid[Oid()]
    @check ccall((:git_oid_fromstrp, libgit2), Cint,
                 (Ptr{Oid}, Ptr{Cchar}), oid_ptr, bstr)
    return oid_ptr[1]
end

Oid(id::Oid) = id

Base.copy!(arr::Array{UInt8, 1}, id::Oid) = begin
    unsafe_copy!(convert(Ptr{UInt8}, arr), convert(Ptr{UInt8},  hex2bytes(hex(id))), OID_RAWSZ)
    return arr
end

raw(id::Oid) = copy!(Array(UInt8, OID_RAWSZ), id)

const _hexstr = Array(UInt8, OID_HEXSZ)
Base.hex(id::Oid) = begin
    ccall((:git_oid_nfmt, libgit2), Void,
          (Ptr{UInt8}, Csize_t, Ptr{Oid}), _hexstr, OID_HEXSZ, &id)
    return bytestring(convert(Ptr{UInt8}, _hexstr), OID_HEXSZ)
end

Base.string(id::Oid) = hex(id)

Base.show(io::IO, id::Oid) = print(io, "Oid($(string(id)))")

Base.hash(id::Oid) = hash(hex(id))

Base.cmp(id1::Oid, id2::Oid) = int(ccall((:git_oid_cmp, libgit2), Cint,
                                         (Ptr{Oid}, Ptr{Oid}), &id1, &id2))

Base.(:(==))(id1::Oid, id2::Oid) = cmp(id1, id2) == 0
Base.isequal(id1::Oid, id2::Oid) = cmp(id1, id2) == 0
Base.isless(id1::Oid, id2::Oid)  = cmp(id1, id2) < 0

iszero(id::Oid) = begin
    bytes = raw(id)
    for i=1:OID_RAWSZ
        if bytes[i] != 0x00
            return false
        end
    end
    return true
end

#throw(ArgumentError("not a hexadecimal string: $(repr(s))"))

function allhexchar(s::ByteString)
    for c in s
        if !('0' <= c <= '9' || 'a' <= c <= 'f' || 'A' <= c <= 'F')
            return false
        end
    end
    return true
end

immutable PartialSHA1
    data::ASCIIString

    PartialSHA1(s::AbstractString) = begin
        len = length(s)
        if len == 0
            throw(ArgumentError("PartialSHA1 string cannot be empty"))
        elseif len >= OID_HEXSZ
            throw(ArgumentError("PartialSHA1 string length must be <= $(OID_HEXSZ)"))
        end
        str = ASCIIString(s)
        if !allhexchar(str)
            throw(ArgumentError("SHA1 string is not a hexadecimal string: $(repr(str))"))
        end
        return new(str)
    end
end

immutable SHA1
    data::ASCIIString

    SHA1(s::AbstractString) = begin
        len = length(s)
        if len == 0
            throw(ArgumentError("SHA1 string cannot be empty"))
        elseif len < OID_HEXSZ
            throw(ArgumentError("SHA1 string length < $(OID_HEXSZ)"))
        elseif len > OID_HEXSZ
            throw(ArgumentError("SHA1 string length > $(OID_HEXSZ)"))
        end
        str = ASCIIString(s)
        if !allhexchar(str)
            throw(ArgumentError("SHA1 string is not a hexadecimal string: $(repr(str))"))
        end
        return new(str)
    end
end

macro sha1_str(s)
    len = length(s)
    if len < OID_HEXSZ
        :(PartialSHA1($s))
    else
        :(SHA1($s))
    end
end

Oid(sha::SHA1) = Oid(hex2bytes(sha.data))


text:
function partitionnementMonoObjectif(instance)
    println("\n Choix du partitionnement exact :")
    println(" --------------------------------")
    println("  1) Agatz Exact Partitionning (AEP = A1)")
    println("  2) Poikonen Exact Partitionning (EP = A2)")
    choix = choixListe(2)
    if choix == "1"
        AEP(instance)
    elseif choix == "2"
        EP(instance)
    end
end

function tempsOperation(instance, i,j,k)
    solutionTSP = instance.solution.ordreVisiteCamion # Solution TSP : ordre de visite  du camion
    # Debug
    #println(" Calcul du temps de l'opération : (i=$i, j=$j, k=$k)")
    #println(" Solution TSP : R_tsp = $solutionTSP")
    #print(" r_$i = $(solutionTSP[i]);")
    #print(" r_$j = $(solutionTSP[j]);")
    #if k != -1
    #    println(" r_$k = $(solutionTSP[k])")
    #end
    # Calcul distance parcourue par le drone
    if k == -1 # Pas de noeud drone
        distanceDrone = 0
        ri = solutionTSP[i]
        rj = solutionTSP[j]
        distanceCamion = instance.D[ri,rj]
        tempsAttente = 0.0
    else
        ri = solutionTSP[i]
        rj = solutionTSP[j]
        rk = solutionTSP[k]
        distanceDrone = instance.D[ri,rk] + instance.D[rk,rj]
        distanceCamion = 0
        for l in i:(k-2)
            r1 = solutionTSP[l]
            r2 = solutionTSP[l+1]
            distanceCamion += instance.D[r1,r2]
        end
        rkAvant = solutionTSP[k-1]
        rkApres = solutionTSP[k+1]
        distanceCamion += instance.D[rkAvant,rkApres]
        for l in (k+1):(j-1)
            r1 = solutionTSP[l]
            r2 = solutionTSP[l+1]
            distanceCamion += instance.D[r1,r2]
        end
        tempsAttente = abs((distanceCamion/instance.vitesseCamion) - (distanceDrone/instance.vitesseDrone))
    end
    tempsParcours = max(distanceCamion/instance.vitesseCamion, distanceDrone/instance.vitesseDrone)
    return tempsParcours, tempsAttente
    #println(" Distance parcourue par le drone : $distanceDrone")
    #println(" Distance parcourue par le camion : $distanceCamion")
end

function calculMeilleuresOperations(instance)
    solutionTSP = instance.solution.ordreVisiteCamion
    # Permet de calculer la meilleure opération pour aller de ri à rj pour toutes les sous-séquences (i,...,j) de solutionTSP
    tempsMinIJ = Array{Float64,2}(undef, length(solutionTSP), length(solutionTSP)) # Temps de la meilleure opération pour aller de i à j
    attenteIJ = Array{Float64,2}(undef, length(solutionTSP), length(solutionTSP)) # Attente liée à la meilleure opération pour aller de i à j
    operationMinIJ = Array{Int64,2}(undef, length(solutionTSP), length(solutionTSP)) # Noeud drone correspondant à la meilleure opération pour aller de i à j (-1 si pas de noeud drone)
    for i in 1:length(solutionTSP)
        for j in 1:length(solutionTSP)
            tempsMinIJ[i,j] = -Inf32
            operationMinIJ[i,j] = -100
        end
    end
    for i in 1:length(solutionTSP)-1
        for j in i+1:length(solutionTSP)
            if j == i+1 # k = -1
                tempsMinIJ[i, j], attenteIJ[i, j] = tempsOperation(instance, i, j, -1)
                operationMinIJ[i, j] = -1
            else
                minTij = Inf32
                attente = Inf32
                operation = -100
                # k > 0
                for k in (i+1):(j-1)
                    Tij, Aij = tempsOperation(instance, i, j, k)
                    if Tij < minTij
                        minTij = Tij
                        attente = Aij
                        operation = k
                    end
                end
                tempsMinIJ[i, j] = minTij
                attenteIJ[i, j] = attente
                operationMinIJ[i, j] = operation
            end
        end
    end
    return tempsMinIJ, operationMinIJ, attenteIJ
end

function AEP(instance)
    tempsMinIJ, operationMinIJ, attenteIJ = calculMeilleuresOperations(instance)
    V = Vector{Float64}(undef, length(instance.solution.ordreVisiteCamion))
    A = Vector{Float64}(undef, length(instance.solution.ordreVisiteCamion))
    P = Vector{Int64}(undef, length(instance.solution.ordreVisiteCamion))
    V[1] = 0
    A[1] = 0
    P[1] = -1

    for i in 2:length(V)
        minTemps = Inf32
        attenteAssociee = Inf32
        pred = -1
        for k in 1:(i-1)
            temps = V[k] + tempsMinIJ[k,i]
            attente = A[k] + attenteIJ[k,i]
            if temps < minTemps
                minTemps = temps
                attenteAssociee = attente
                pred = k
            end
        end
        V[i] = minTemps
        A[i] = attenteAssociee
        P[i] = pred
    end
    #println(" V = $V")
    #println(" A = $A")
    #println(" P = $P")

    temp = []
    pred = P[end]
    push!(temp, length(instance.solution.ordreVisiteCamion))
    push!(temp, pred)
    while pred != 1
        pred = P[pred]
        push!(temp, pred)
    end
    temp = reverse(temp)

    ordreVisiteCamion = []
    operationsDrone = []
    push!(ordreVisiteCamion, 1)
    for i in 1:length(temp)-1
        operation = operationMinIJ[temp[i],temp[i+1]]
        if operation == -1 # Trajet direct
            noeudCombine = instance.solution.ordreVisiteCamion[temp[i+1]]
            push!(ordreVisiteCamion, noeudCombine)
        else # Opération drone
            noeudDrone = instance.solution.ordreVisiteCamion[operation]
            noeudDepart = instance.solution.ordreVisiteCamion[temp[i]]
            noeudFin = instance.solution.ordreVisiteCamion[temp[i+1]]
            push!(operationsDrone, [noeudDepart, noeudDrone, noeudFin])
            for l in temp[i]+1:temp[i+1]-1
                noeudActu = instance.solution.ordreVisiteCamion[l]
                if noeudActu != noeudDrone # Noeud camion
                    push!(ordreVisiteCamion, noeudActu)
                end
            end
            push!(ordreVisiteCamion, noeudFin)
        end
    end

    instance.solution.tempsParcours = V[end]
    instance.solution.tempsAttente = A[end]
    instance.solution.ordreVisiteCamion = ordreVisiteCamion
    instance.solution.operationsDrone = operationsDrone
    instance.codeResolution = 3
end

#################### Poikonen (EP) ####################

function tempsOperationEP(instance, i,j,k)
    solutionTSP = instance.solution.ordreVisiteCamion # Solution TSP : ordre de visite  du camion
    # Calcul distance parcourue par le drone
    bestY = -100 # Pour le cas k == -3
    bestZ = -100
    if k == -1 # Pas de noeud drone
        ri = solutionTSP[i]
        rj = solutionTSP[j]
        distanceCamion = instance.D[ri,rj]
        tempsParcours = distanceCamion/instance.vitesseCamion
        tempsAttente = 0.0
    elseif k == -2
        ri = solutionTSP[i]
        rj = solutionTSP[j]
        # Le camion va de ri à rj
        distanceCamion = instance.D[ri,rj]
        distanceDrone = 0
        for l in (i+1):(j-1)
            rl = solutionTSP[l]
            distanceDrone += 2 * instance.D[ri,rl]
        end
        tempsParcours = (distanceCamion/instance.vitesseCamion) + (distanceDrone/instance.vitesseDrone)
        tempsAttente = (distanceDrone/instance.vitesseDrone) # Le camion doit attendre que le drone termine la livraison
    elseif k == -3
        minTempsParcours = Inf32
        minTempsAttente = Inf32
        tempsP = 0.0
        tempsA = 0.0
        for z in (i+1):(j-2)
            for y in (z+1):(j-1)
                tempsP, tempsA = T(instance,i,j,z,y)
                if tempsP < minTempsParcours
                    minTempsParcours = tempsP
                    minTempsAttente = tempsA
                    bestY = y
                    bestZ = z
                end
            end
        end
        tempsParcours = tempsP
        tempsAttente = tempsA
    else
        ri = solutionTSP[i]
        rj = solutionTSP[j]
        rk = solutionTSP[k]
        distanceDrone = instance.D[ri,rk] + instance.D[rk,rj]
        distanceCamion = 0
        for l in i:(k-2)
            r1 = solutionTSP[l]
            r2 = solutionTSP[l+1]
            distanceCamion += instance.D[r1,r2]
        end
        rkAvant = solutionTSP[k-1]
        rkApres = solutionTSP[k+1]
        distanceCamion += instance.D[rkAvant,rkApres]
        for l in (k+1):(j-1)
            r1 = solutionTSP[l]
            r2 = solutionTSP[l+1]
            distanceCamion += instance.D[r1,r2]
        end
        tempsParcours = max(distanceCamion/instance.vitesseCamion, distanceDrone/instance.vitesseDrone)
        tempsAttente = abs((distanceCamion/instance.vitesseCamion) - (distanceDrone/instance.vitesseDrone))
    end
    return tempsParcours, tempsAttente, bestY, bestZ
    #println(" Distance parcourue par le drone : $distanceDrone")
    #println(" Distance parcourue par le camion : $distanceCamion")
end

function T(instance, i,j,z,y) # k = -3
    solutionTSP = instance.solution.ordreVisiteCamion # Solution TSP : ordre de visite  du camion
    ri = solutionTSP[i]
    rj = solutionTSP[j]
    # Allers-retours drone
    distanceDroneEtape1 = 0
    for l in (i+1):z
        rl = solutionTSP[l]
        distanceDroneEtape1 += 2 * instance.D[ri,rl]
    end
    # Le reste
    ry = solutionTSP[y]
    distanceDroneEtape2 = instance.D[ri,ry] + instance.D[ry,rj]
    distanceCamion = 0.0
    if y == z+1 && j != y+1 ## Cas w1 vide (i -> i+1 -> ... -> z -> y -> w_2_1 -> ... -> j)
        rw21 = solutionTSP[y+1]
        distanceCamion = instance.D[ri,rw21]
        for l in (y+1):(j-1)
            rw2_l1 = solutionTSP[l]
            rw2_l2 = solutionTSP[l+1]
            distanceCamion += instance.D[rw2_l1,rw2_l2]
        end
    elseif y != z+1 && j == y+1 ## Cas w2 vide (i -> i+1 -> ... -> z -> w_1_1 > ... -> y -> j)
        rw11 = solutionTSP[z+1]
        distanceCamion = instance.D[ri,rw11]
        for l in (z+1):(y-2) # On ne livre pas y
            rw1_l1 = solutionTSP[l]
            rw1_l2 = solutionTSP[l+1]
            distanceCamion += instance.D[rw1_l1,rw1_l2]
        end
        rw1 = solutionTSP[y-1]
        distanceCamion += instance.D[rw1,rj]
    elseif y == z+1 && j == y+1 ## Cas w1 et w2 vides (i -> i+1 -> ... -> z -> y -> j)
        distanceCamion = instance.D[ri,rj]
    elseif y != z+1 && j != y+1 ## Cas w1 et w2 avec au moins 1 noeud (i -> i+1 -> ... -> z -> w_1_1 -> ... -> y -> w_2_1 -> ... -> j)
        rw11 = solutionTSP[z+1]
        distanceCamion = instance.D[ri,rw11]
        for l in (z+1):(y-2) # On ne livre pas y
            rw1_l1 = solutionTSP[l]
            rw1_l2 = solutionTSP[l+1]
            distanceCamion += instance.D[rw1_l1,rw1_l2]
        end
        rw1 = solutionTSP[y-1]
        rw21 = solutionTSP[y+1]
        distanceCamion = instance.D[rw1,rw21]
        for l in (y+1):(j-1)
            rw2_l1 = solutionTSP[l]
            rw2_l2 = solutionTSP[l+1]
            distanceCamion += instance.D[rw2_l1,rw2_l2]
        end
    end
    tempsParcours = distanceDroneEtape1/instance.vitesseDrone + max(distanceDroneEtape2/instance.vitesseDrone, distanceCamion/instance.vitesseCamion)
    tempsAttente = (distanceDroneEtape1/instance.vitesseDrone) + abs(distanceDroneEtape2/instance.vitesseDrone - distanceCamion/instance.vitesseCamion)
    return tempsParcours, tempsAttente
end

function calculMeilleuresOperationsEP(instance)
    solutionTSP = instance.solution.ordreVisiteCamion
    # Permet de calculer la meilleure opération pour aller de ri à rj pour toutes les sous-séquences (i,...,j) de solutionTSP
    tempsMinIJ = Array{Float64,2}(undef, length(solutionTSP), length(solutionTSP)) # Temps de la meilleure opération pour aller de i à j
    attenteIJ = Array{Float64,2}(undef, length(solutionTSP), length(solutionTSP)) # Attente liée à la meilleure opération pour aller de i à j
    operationMinIJ = Array{Vector{Int64},2}(undef, length(solutionTSP), length(solutionTSP)) # Noeud drone correspondant à la meilleure opération pour aller de i à j (-1 si pas de noeud drone)
    bestY = -100
    bestZ = -100
    for i in 1:length(solutionTSP)
        for j in 1:length(solutionTSP)
            tempsMinIJ[i,j] = -Inf32
            operationMinIJ[i,j] = [-100]
        end
    end
    for i in 1:length(solutionTSP)-1
        for j in i+1:length(solutionTSP)
            if j == i+1 # k = -1
                tempsMinIJ[i, j], attenteIJ[i, j] = tempsOperationEP(instance, i, j, -1)
                operationMinIJ[i, j] = [-1]
            else
                minTij = Inf32
                attente = Inf32
                operation = -100
                # k > 0
                for k in (i+1):(j-1)
                    Tij, Aij = tempsOperationEP(instance, i, j, k)
                    if Tij < minTij
                        minTij = Tij
                        attente = Aij
                        operation = k
                    end
                end
                # k = -2
                Tij, Aij = tempsOperationEP(instance, i, j, -2)
                if Tij < minTij
                    minTij = Tij
                    attente = Aij
                    operation = -2
                end
                # k = -3
                if j >= i+3 # Minimum pour placer z et y
                    Tij, Aij, bestY, bestZ = tempsOperationEP(instance, i, j, -3)
                    if Tij < minTij
                        minTij = Tij
                        attente = Aij
                        operation = -3
                    end
                end
                tempsMinIJ[i, j] = minTij
                attenteIJ[i, j] = attente
                operationMinIJ[i, j] = [operation, bestY, bestZ] # Pour gérer le cas k = -3
            end
        end
    end
    return tempsMinIJ, operationMinIJ, attenteIJ
end

function EP(instance)
    tempsMinIJ, operationMinIJ, attenteIJ = calculMeilleuresOperationsEP(instance)
    V = Vector{Float64}(undef, length(instance.solution.ordreVisiteCamion))
    A = Vector{Float64}(undef, length(instance.solution.ordreVisiteCamion))
    P = Vector{Int64}(undef, length(instance.solution.ordreVisiteCamion))
    V[1] = 0
    A[1] = 0
    P[1] = -1
    for i in 2:length(V)
        minTemps = Inf32
        attenteAssociee = Inf32
        pred = -1
        for k in 1:(i-1)
            temps = V[k] + tempsMinIJ[k,i]
            attente = A[k] + attenteIJ[k,i]
            if temps < minTemps
                minTemps = temps
                attenteAssociee = attente
                pred = k
            end
        end
        V[i] = minTemps
        A[i] = attenteAssociee
        P[i] = pred
    end
    #println(" V = $V")
    #println(" A = $A")
    #println(" P = $P")

    temp = []
    pred = P[end]
    push!(temp, length(instance.solution.ordreVisiteCamion))
    push!(temp, pred)
    while pred != 1
        pred = P[pred]
        push!(temp, pred)
    end
    temp = reverse(temp)

    ordreVisiteCamion = []
    operationsDrone = []
    push!(ordreVisiteCamion, 1)
    for i in 1:length(temp)-1
        operation = operationMinIJ[temp[i],temp[i+1]][1] # (operation, bestY, bestZ)
        if operation == -1 # Trajet direct
            noeudCombine = instance.solution.ordreVisiteCamion[temp[i+1]]
            push!(ordreVisiteCamion, noeudCombine)
        elseif operation == -2
            noeudDepart = instance.solution.ordreVisiteCamion[temp[i]]
            noeudFin = instance.solution.ordreVisiteCamion[temp[i+1]]
            push!(ordreVisiteCamion, noeudFin)
            for l in temp[i]+1:temp[i+1]-1
                noeudDrone = instance.solution.ordreVisiteCamion[l]
                push!(operationsDrone, [noeudDepart, noeudDrone, noeudDepart]) # [i,k,j]
            end
        elseif operation == -3
            ri = instance.solution.ordreVisiteCamion[temp[i]]
            rj = instance.solution.ordreVisiteCamion[temp[i+1]]
            ry = instance.solution.ordreVisiteCamion[operationMinIJ[temp[i],temp[i+1]][2]]
            rz = instance.solution.ordreVisiteCamion[operationMinIJ[temp[i],temp[i+1]][3]]
            #println(" Ordre visite camion TSP : ", instance.solution.ordreVisiteCamion)
            #println(" ri = $ri")
            #println(" rj = $rj")
            #println(" ry = $ry")
            #println(" rz = $rz")
            #i = temp[i]
            #j = temp[i+1]
            y = operationMinIJ[temp[i],temp[i+1]][2]
            z = operationMinIJ[temp[i],temp[i+1]][3]
            # Etape 1 : Allers-retours du drone
            for l in (temp[i]+1):z
                noeudDrone = instance.solution.ordreVisiteCamion[l]
                push!(operationsDrone, [ri, noeudDrone, ri]) # [i,k,j]
            end
            # Etape 2
            if y == z+1 && temp[i+1] == y+1 ## Cas w1 et w2 vides (i -> i+1 -> ... -> z -> y -> j)
                #println(" --> Cas w1 et w2 vides")
                push!(ordreVisiteCamion, rj) # Le camion va en j directement
                push!(operationsDrone, [ri, ry, rj]) # Le drone livre y puis rejoint le camion en j
            elseif y != z+1 && temp[i+1] != y+1 ## Cas w1 et w2 avec au moins 1 noeud (i -> i+1 -> ... -> z -> w_1_1 -> ... -> y -> w_2_1 -> ... -> j)
                push!(operationsDrone, [ri, ry, rj]) # Le drone livre y puis rejoint le camion en j
                for l in (z+1):(y-1) # On ne livre pas y
                    rw1_l1 = instance.solution.ordreVisiteCamion[l]
                    push!(ordreVisiteCamion, rw1_l1) # Le camion va en j directement
                end
                for l in (y+1):temp[i+1]
                    rw2_l1 = instance.solution.ordreVisiteCamion[l]
                    push!(ordreVisiteCamion, rw2_l1)
                end
                #println(" --> Cas w1 et w2 avec au moins 1 noeud")
            elseif y == z+1 && temp[i+1] != y+1 ## Cas w1 vide (i -> i+1 -> ... -> z -> y -> w_2_1 -> ... -> j)
                #println(" --> Cas w1 vide")
                push!(operationsDrone, [ri, ry, rj]) # Le drone livre y puis rejoint le camion en j
                for l in (y+1):temp[i+1]
                    rw2_l1 = instance.solution.ordreVisiteCamion[l]
                    push!(ordreVisiteCamion, rw2_l1)
                end
            elseif y != z+1 && temp[i+1] == y+1 ## Cas w2 vide (i -> i+1 -> ... -> z -> w_1_1 > ... -> y -> j)
                #println(" --> Cas w2 vide")
                push!(operationsDrone, [ri, ry, rj]) # Le drone livre y puis rejoint le camion en j
                for l in (z+1):(y-1) # On ne livre pas y
                    rw1_l1 = instance.solution.ordreVisiteCamion[l]
                    push!(ordreVisiteCamion, rw1_l1) # Le camion va en j directement
                end
            end
        else # Opération drone
            noeudDrone = instance.solution.ordreVisiteCamion[operation]
            noeudDepart = instance.solution.ordreVisiteCamion[temp[i]]
            noeudFin = instance.solution.ordreVisiteCamion[temp[i+1]]
            push!(operationsDrone, [noeudDepart, noeudDrone, noeudFin]) # [i,k,j]
            for l in temp[i]+1:temp[i+1]-1
                noeudActu = instance.solution.ordreVisiteCamion[l]
                if noeudActu != noeudDrone # Noeud camion
                    push!(ordreVisiteCamion, noeudActu)
                end
            end
            push!(ordreVisiteCamion, noeudFin)
        end
    end

    instance.solution.tempsParcours = V[end]
    instance.solution.tempsAttente = A[end]
    instance.solution.ordreVisiteCamion = ordreVisiteCamion
    instance.solution.operationsDrone = operationsDrone
    instance.codeResolution = 4
end

max_stars_repo_name: xgandibleux/tspd-solver

max_issues_repo_name: xgandibleux/tspd-solver

max_forks_repo_name: xgandibleux/tspd-solver


text:
# Simulate markov decision process
# given a defined state graph

module MDPSimulator
export MDPGraph, createMDPGraph, addMDPGraphTransition!, copyMDPGraph, nextMDP!, runMDP!, displayMDP

"""
Graph structure containing relations between states for MDP\\
`terminal` defines the terminal state name\\
`states` contains all possible states\\
`actions` contains all possible actions\\
`relations` contains each trasition (S,A,S')\\
`rewards` defines reward for each transition\\
`probs` defines probability of each transition (S,A,S')\\
`values` contains value for each state\\
`decided` contains decided action for each state\\
`gamma` is the discount of next state's value\\
`iternum` is the number of iterations performed
"""
mutable struct MDPGraph
    terminal::String
    states::Array{String}
    actions::Array{String}
    relations::Array{Tuple{String,String,String}}
    rewards::Dict{Tuple{String,String,String},AbstractFloat}
    probs::Dict{Tuple{String,String,String},AbstractFloat}
    values::Dict{String,AbstractFloat}
    decided::Dict{String,Union{String,Nothing}}
    gamma::AbstractFloat
    iternum::Integer
end

"""
MDPGraph creation function
"""
function createMDPGraph(states::Array{String}, actions::Array{String}, terminals::Array{String}, terminal_vals::Array{T},
        gamma::T)::MDPGraph where T<:AbstractFloat
    @assert length(terminals) == length(terminal_vals)
    @assert length(states) > 0
    @assert length(actions) > 0
    # create a terminal state name
    terminal = "T"
    while terminal in states
        terminal = terminal * terminal
    end
    states = copy(states)
    push!(states, terminal)
    # prepare data
    relations::Array{Tuple{String,String,String}} = []
    rewards::Dict{Tuple{String,String,String},AbstractFloat} = Dict()
    probs::Dict{Tuple{String,String,String},AbstractFloat} = Dict()
    values::Dict{String,AbstractFloat} = Dict()
    decided::Dict{String,Union{String,Nothing}} = Dict()
    # initialize values and decided actions
    for state in states
        values[state] = 0.0
        decided[state] = nothing
    end
    for i in 1:length(terminals)
        push!(relations, (terminals[i], actions[1], terminal))
        probs[(terminals[i], actions[1], terminal)] = 1.0
        rewards[(terminals[i], actions[1], terminal)] = terminal_vals[i]
    end
    graph = MDPGraph(terminal, states, copy(actions), relations, rewards, probs, values, decided, gamma, 0)
    return graph
end

"""
MDPGraph copy function
"""
function copyMDPGraph(graph::MDPGraph)::MDPGraph
    newGraph = MDPGraph(graph.terminal, copy(graph.states), copy(graph.actions), copy(graph.relations),
        copy(graph.rewards), copy(graph.probs), copy(graph.values), copy(graph.decided), graph.gamma, graph.iternum)
    return newGraph
end

"""
MDPGraph add transition
"""
function addMDPGraphTransition!(graph::MDPGraph, state_0::String, action::String, state_1::String;
        reward::T=0.0, prob::T=0.5) where T<:AbstractFloat
    @assert state_0 != graph.terminal
    @assert state_0 in graph.states
    @assert state_1 in graph.states
    @assert action in graph.actions
    @assert 0.0 <= prob <= 1.0
    push!(graph.relations, (state_0, action, state_1))
    graph.rewards[(state_0, action, state_1)] = reward
    graph.probs[(state_0, action, state_1)] = prob
    return nothing
end

"""
Get possible actions given a state
"""
function getPossibleActions(graph::MDPGraph, state::String)::Array
    @assert state in graph.states
    actions = []
    for (s0, a, s1) in graph.relations
        if s0 == state
            push!(actions, a)
        end
    end
    return actions
end

"""
Get possible states given a state and action
"""
function getPossibleStates(graph::MDPGraph, state::String, action::String)::Array
    @assert state in graph.states
    @assert action in graph.actions
    states = []
    for (s0, a, s1) in graph.relations
        if s0 == state && a == action
            push!(states, s1)
        end
    end
    return states
end

"""
Compute Q value given a state and action
"""
function computeQValue(graph::MDPGraph, state::String, action::String)::AbstractFloat
    @assert state in graph.states
    @assert action in graph.actions
    result = 0.0
    for nextState in getPossibleStates(graph, state, action)
        @assert (state, action, nextState) in graph.relations
        @assert (state, action, nextState) in keys(graph.probs)
        @assert (state, action, nextState) in keys(graph.rewards)
        p = graph.probs[(state, action, nextState)]
        r = graph.rewards[(state, action, nextState)]
        result += p * (r + graph.gamma * graph.values[nextState])
    end
    return result
end

"""
Get best action with corresponding Q value
"""
function getBestActionWithQValue(graph::MDPGraph, state::String)::Pair{String,AbstractFloat}
    @assert state in graph.states
    @assert state != graph.terminal
    actions = getPossibleActions(graph, state)
    @assert length(actions) > 0
    result::Dict{String,AbstractFloat} = Dict()
    for action in actions
        result[action] = computeQValue(graph, state, action)
    end
    return sort(collect(result), by=m->m[2], rev=true)[1]
end

"""
Run 1 MDP iteration on MDPGraph, updating values in batch
"""
function nextMDP!(graph::MDPGraph)::Bool
    new_values = copy(graph.values)
    for state in graph.states
        # if terminal state, skip
        if state == graph.terminal
            continue
        end
        action, qval = getBestActionWithQValue(graph, state)
        graph.decided[state] = action
        new_values[state] = qval
    end
    converged = (graph.values == new_values)
    graph.values = new_values
    graph.iternum += 1
    return converged
end

"""
Run MDP by multiple iterations
"""
function runMDP!(graph::MDPGraph, iterations::Integer)
    @assert iterations >= 0
    for i in 1:iterations
        if nextMDP!(graph)
            println("MDP Converged at iter = $(i-1)")
            break
        end
    end
end

"""
Display computed values and decisions for each state
"""
function displayMDP(graph::MDPGraph)
    iternum = graph.iternum
    println("MDP Display Start (Iter = $iternum)")
    for state in graph.states
        if state == graph.terminal
            continue
        end
        print("State = $state")
        val = graph.values[state]
        action = graph.decided[state]
        if (state, graph.actions[1], graph.terminal) in graph.relations
            print(" (terminal)")
        end
        println("\n\tValue  = $val")
        println("\tAction = $action")
    end
    println("MDP Display End")
end

end

max_stars_repo_name: teamclouday/MDPSimulator

max_issues_repo_name: teamclouday/MDPSimulator

max_forks_repo_name: teamclouday/MDPSimulator


text:
# Jonas August, 2/2022

#=
This problem was asked by Google.

Given the head of a singly linked list, reverse it in-place.
=#

mutable struct Node
    value
    next::Union{Nothing, Node}
end

"
Reverse singly-linked list in place.

           l
nothing    n1 -> n2 -> n3 -> ... -> nmax -> nothing
rev        temp

                 l
nothing    n1 -> n2 -> n3 -> ... -> nmax -> nothing
rev        temp

                 l
nothing <- n1    n2 -> n3 -> ... -> nmax -> nothing
rev        temp

                 l
nothing <- n1    n2 -> n3 -> ... -> nmax -> nothing
           rev   temp
"
function reverselist(l::Union{Nothing, Node})
    rev = nothing
    while !isnothing(l)
        temp = l
        l = l.next
        temp.next = rev
        rev = temp
    end
    rev
end

function printlist(l::Union{Nothing, Node})
    while !isnothing(l)
        println(l.value)
        l = l.next
    end
end

lists = [nothing, Node(1, nothing), Node(1, Node(3, Node(5, nothing)))]
for l in lists
    @show l
    printlist(l)
    println()
    rev = reverselist(l)
    @show rev
    printlist(rev)
    println()
    println()
end

max_stars_repo_path: Reverse singly linked list in place - Problem 73.jl


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 45beb2f8-bc6f-11eb-3483-71eebecb5732
using Luxor

# ╔═╡ bfe33804-3d25-4820-b316-aa98231575a6
md"Based on [inconvergent/hyphae](https://inconvergent.net/generative/hyphae/)"

# ╔═╡ dfee5344-7984-4a0d-a299-9c937c1fddf8
@draw begin
	line(Point(0, 0), Point(100, 0), :stroke)
end

# ╔═╡ 507ca324-a3bd-48aa-b78a-b090e3227b56

# ╔═╡ Cell order:
# ╠═bfe33804-3d25-4820-b316-aa98231575a6
# ╠═45beb2f8-bc6f-11eb-3483-71eebecb5732
# ╠═dfee5344-7984-4a0d-a299-9c937c1fddf8
# ╠═507ca324-a3bd-48aa-b78a-b090e3227b56


text:
function DLDR(ALPHA,BETA)
      A=[0.0050    0.0070    0.0130    0.0180    0.0150    0.0210    0.0230
         0.0170    0.0160    0.0130    0.0150    0.0140    0.0110    0.0100
         0.0140    0.0140    0.0110    0.0150    0.0130    0.0100    0.0110
         0.0100    0.0140    0.0120    0.0140    0.0130    0.0110    0.0110
        -0.0050    0.0130    0.0110    0.0140    0.0120    0.0100    0.0110
         0.0090    0.0090    0.0090    0.0140    0.0110    0.0090    0.0100
         0.0190    0.0120    0.0080    0.0140    0.0110    0.0080    0.0080
         0.0050    0.0050    0.0050    0.0150    0.0100    0.0100    0.0100
              0         0   -0.0020    0.0130    0.0080    0.0060    0.0060
        -0.0050    0.0040    0.0050    0.0110    0.0080    0.0050    0.0140
        -0.0110    0.0090    0.0030    0.0060    0.0070         0    0.0200
         0.0080    0.0070    0.0050    0.0010    0.0030    0.0010         0]
      #
      S= 0.2 * ALPHA
      K= floor(S)+3
      if(K <= 1)
            K= 2
      end
      if(K >=  12)
            K=  11
      end
      DA= S - K +3
      L = K + fix( 1.1*sign(DA) )
      #
      S= 0.1 * BETA
      M= floor(S)+4
      if(M <= 1)
            M= 2
      end
      if(M >=  7)
            M=  6
      end
      DB= S - M + 4
      N= M + fix( 1.1*sign(DB) )
      #
      K=convert(Int64, K)
      L=convert(Int64, L)
      M=convert(Int64, M)
      N=convert(Int64, N)
      T= A[K,M]
      U= A[K,N]
      V= T + abs(DA) * (A[L,M] - T)
      W= U + abs(DA) * (A[L,N] - U)
      return dldrv = V + (W-V)  * abs(DB)
end


text:
###### type definition ##########

using Base: @propagate_inbounds

import Base: convert, copy, length, show, getindex, iterate,
             lastindex, size, eachindex, ==, isequal, hash, ndims,
             getproperty, propertynames, values

abstract type AbstractTimeSeries{T,N,D} end

"""
    TimeArray{T,N,D<:TimeType,A<:AbstractArray{T,N}} <: AbstractTimeSeries{T,N,D}

# Constructors

    TimeArray(timestamp, values[, colnames, meta = nothing])
    TimeArray(ta::TimeArray; timestamp, values, colnames, meta)
    TimeArray(data::NamedTuple, timestamp = :datetime, meta)
    TimeArray(table; timestamp::Symbol, timeparser::Callable = identity)

The second constructor will yields a new TimeArray with the new given fields.
Note that the unchanged fields will be shared, there aren't any copy for the
underlying arrays.

The third constructor builds a `TimeArray` from a `NamedTuple`.

# Arguments

- `timestamp::AbstractVector{<:TimeType}`: a vector of sorted timestamps,
  Each element in this vector should be unique.

- `timestamp::Symbol`: the column name of the time index from the source table.
  The constructor is used for the Tables.jl package integration.

- `values::AbstractArray`: a data vector or matrix. Its number of rows
  should match the length of `timestamp`.

- `colnames::Vector{Symbol}`: the column names. Its length should match
  the column of `values`.

- `meta::Any`: a user-defined metadata.

- `timeparser::Callable`: a mapping function for converting the source time index.
  For instance, `Dates.unix2datetime` is a common case.

# Examples

    data = (datetime = [DateTime(2018, 11, 21, 12, 0), DateTime(2018, 11, 21, 13, 0)],
            col1 = [10.2, 11.2],
            col2 = [20.2, 21.2],
            col3 = [30.2, 31.2])
    ta = TimeArray(data; timestamp = :datetime, meta = "Example")

"""
struct TimeArray{T,N,D<:TimeType,A<:AbstractArray{T,N}} <: AbstractTimeSeries{T,N,D}

    timestamp::Vector{D}
    values::A
    colnames::Vector{Symbol}
    meta::Any

    function TimeArray{T,N,D,A}(
            timestamp::AbstractVector{D},
            values::A,
            colnames::Vector{Symbol},
            meta::Any;
            unchecked = false) where {T,N,D<:TimeType,A<:AbstractArray{T,N}}
        nrow = size(values, 1)
        ncol = size(values, 2)
        colnames = copy(colnames)

        unchecked && return new(timestamp, values, replace_dupes!(colnames), meta)

        nrow != length(timestamp) && throw(DimensionMismatch("values must match length of timestamp"))
        ncol != length(colnames) && throw(DimensionMismatch("column names must match width of array"))

        _issorted_and_unique(timestamp) && return new(
            timestamp, values, replace_dupes!(colnames), meta)

        timestamp_r = reverse(timestamp)
        _issorted_and_unique(timestamp_r) && return new(
            timestamp_r, reverse(values, dims = 1), replace_dupes!(colnames), meta)

        throw(ArgumentError("timestamps must be strictly monotonic"))
    end
end

###### outer constructor ########

TimeArray(d::AbstractVector{D}, v::AbstractArray{T,N},
          c::Vector{Symbol} = gen_colnames(size(v, 2)),
          m::Any = nothing; args...) where {T,N,D<:TimeType} =
    TimeArray{T,N,D,typeof(v)}(d, v, c, m; args...)

TimeArray(d::D, v::AbstractArray{T,N},
          c::Vector{Symbol} = gen_colnames(size(v, 2)),
          m::Any = nothing; args...) where {T,N,D<:TimeType} =
    TimeArray{T,N,D,typeof(v)}([d], v, c, m; args...)

TimeArray(ta::TimeArray;
          timestamp = _timestamp(ta), values = _values(ta),
          colnames = _colnames(ta), meta = _meta(ta), args...) =
    TimeArray(timestamp, values, colnames, meta; args...)

function TimeArray(data::NamedTuple; timestamp::Symbol, meta = nothing, args...)
    columns = (key for key in keys(data) if key != timestamp)
    dat = hcat((data[key] for key in columns)...)
    TimeArray(data[timestamp], dat, collect(columns), meta; args...)
end

###### conversion ###############

convert(::Type{TimeArray{Float64,N}}, x::TimeArray{Bool,N}) where N =
    TimeArray(timestamp(x), Float64.(values(x)), colnames(x), meta(x); unchecked = true)

convert(x::TimeArray{Bool,N}) where N =
    convert(TimeArray{Float64,N}, x::TimeArray{Bool,N})

###### copy ###############

copy(ta::TimeArray) =
    TimeArray(timestamp(ta), values(ta), colnames(ta), meta(ta); unchecked = true)

###### length ###################

length(ata::AbstractTimeSeries) = length(timestamp(ata))

###### size #####################

size(ta::TimeArray) = size(values(ta))
size(ta::TimeArray, dim) = size(values(ta), dim)

###### ndims #####################

ndims(ta::AbstractTimeSeries{T,N}) where {T,N} = N

###### iteration protocol ########

@generated function iterate(ta::AbstractTimeSeries{T,N}, i = 1) where {T,N}
    val = (N == 1) ? :(values(ta)[i]) : :(values(ta)[i, :])

    quote
        i > length(ta) && return nothing
        ((timestamp(ta)[i], $val), i + 1)
    end
end

###### equal ####################

"""
    ==(x::TimeArray, y::TimeArray)

If `true`, all fields of `x` and `y` should be equal,
meaning that the two `TimeArray`s have the same values at the same points in time,
the same colnames and the same metadata.

Implies

```julia
x.timestamp == y.timestamp &&
x.values    == y.values    &&
x.colnames  == y.colnames  &&
x.meta      == y.meta
```
"""
==

# Other type info is not helpful for assertion.
# e.g.
#      1.0 == 1
#      Date(2111, 1, 1) == DateTime(2111, 1, 1)
==(x::TimeArray{T,N}, y::TimeArray{S,M}) where {T,S,N,M} = false
==(x::TimeArray{T,N}, y::TimeArray{S,N}) where {T,S,N} =
    all(f -> getfield(x, f) == getfield(y, f), fieldnames(TimeArray))

isequal(x::TimeArray{T,N}, y::TimeArray{S,M}) where {T,S,N,M} = false
isequal(x::TimeArray{T,N}, y::TimeArray{S,N}) where {T,S,N} =
    all(f -> isequal(getfield(x, f), getfield(y, f)), fieldnames(TimeArray))

# support for Dict
hash(x::TimeArray, h::UInt) =
    sum(f -> hash(getfield(x, f), h), fieldnames(TimeArray))

###### eltype #####################

Base.eltype(::AbstractTimeSeries{T,1,D}) where {T,D} = Tuple{D,T}
Base.eltype(::AbstractTimeSeries{T,2,D}) where {T,D} = Tuple{D,Vector{T}}

###### show #####################

@inline _showval(v::Any) = repr(v)
@inline _showval(v::Number) = string(v)
@inline _showval(v::AbstractFloat) =
    ifelse(isnan(v), MISSING, string(round(v, digits=DECIMALS)))

"""
calculate the paging

```
> using MarketData
> AAPL  # this function will return `UnitRange{Int64}[1:9, 10:12]`
```
"""
@inline function _showpages(dcol::Int, timewidth::Int, colwidth::Array{Int})
    ret = UnitRange{Int}[]
    c = dcol - timewidth - 4
    last_i = 1
    for i in eachindex(colwidth)
        w = colwidth[i] + 3
        if c - w < 0
            push!(ret, last_i:i-1)
            # next page
            c = dcol - timewidth - 4 - w
            last_i = i
        elseif i == length(colwidth)
            push!(ret, last_i:i)
        else
            c -= w
        end
    end
    ret
end

function print_time_array(io::IO, ta::TimeArray{T}, short=false) where T
    # summary line
    nrow = size(values(ta), 1)
    ncol = size(values(ta), 2)

    print(io, "$(nrow)×$(ncol) $(typeof(ta))")
    if nrow != 0
        print(io, " $(timestamp(ta)[1]) to $(timestamp(ta)[end])")
    else  # e.g. TimeArray(Date[], [])
        return
    end

    short && return
    println(io)

    # calculate column withs
    drow, dcol = displaysize(io)
    res_row    = 7  # number of reserved rows: summary line, lable line ... etc
    half_row   = floor(Int, (drow - res_row) / 2)
    add_row    = (drow - res_row) % 2

    if nrow > (drow - res_row)
        tophalf = 1:(half_row + add_row)
        bothalf = (nrow - half_row + 1):nrow
        strs = _showval.(@view values(ta)[[tophalf; bothalf], :])
        ts   = @view timestamp(ta)[[tophalf; bothalf]]
    else
        strs = _showval.(values(ta))
        ts   = timestamp(ta)
    end

    colwidth = maximum(
        [textwidth.(string.(colnames(ta)))'; textwidth.(strs); fill(5, ncol)'],
        dims = 1)

    # paging
    spacetime = textwidth(string(ts[1]))
    pages = _showpages(dcol, spacetime, colwidth)

    for p ∈ pages
        # row label line
        ## e.g. | Open  | High  | Low   | Close  |
        print(io, "│", " "^(spacetime + 2))
        for (name, w) in zip(colnames(ta)[p], colwidth[p])
            print(io, "│ ", _rpad(name, w + 1))
        end
        println(io, "│")
        ## e.g. ├───────┼───────┼───────┼────────┤
        print(io, "├", "─"^(spacetime + 2))
        for w in colwidth[p]
            print(io, "┼", "─"^(w + 2))
        end
        print(io, "┤")

        # timestamp and values line
        if nrow > (drow - res_row)
            for i in tophalf
                println(io)
                print(io, "│ ", ts[i], " ")
                for j in p
                    print(io, "│ ", _rpad(strs[i, j], colwidth[j] + 1))
                end
                print(io, "│")
            end

            print(io, "\n   \u22EE")

            for i in (length(bothalf) - 1):-1:0
                i = size(strs, 1) - i
                println(io)
                print(io, "│ ", ts[i], " ")
                for j in p
                    print(io, "│ ", _rpad(strs[i, j], colwidth[j] + 1))
                end
                print(io, "│")
            end

        else
            for i in 1:nrow
                println(io)
                print(io, "│ ", ts[i], " ")
                for j in p
                    print(io, "│ ", _rpad(strs[i, j], colwidth[j] + 1))
                end
                print(io, "│")
            end
        end

        if length(pages) > 1 && p != pages[end]
            print(io, "\n\n")
        end
    end  # for p ∈ pages
end
Base.show(io::IO, ta::TimeArray) = print_time_array(io, ta, true)
Base.show(io::IO, ::MIME"text/plain", ta::TimeArray) =
    print_time_array(io, ta, false)

###### getindex #################

# the getindex function should return a new TimeArray, and copy data from
# the source, includes `timestamp`, `values` and `colnames`.

getindex(ta::TimeArray) = throw(BoundsError(typeof(ta), []))

# single row
@propagate_inbounds getindex(ta::TimeArray, n::Integer) =
    # avoid conversion to column vector
    TimeArray(timestamp(ta)[n], values(ta)[n:n, :], colnames(ta), meta(ta))

# single row 1d
@propagate_inbounds getindex(ta::TimeArray{T,1}, n::Integer) where {T} =
    TimeArray(timestamp(ta)[n], values(ta)[[n]], colnames(ta), meta(ta))

# range of rows
@propagate_inbounds getindex(ta::TimeArray, r::UnitRange{<:Integer}) =
    TimeArray(timestamp(ta)[r], values(ta)[r, :], colnames(ta), meta(ta))

# range of 1d rows
@propagate_inbounds getindex(ta::TimeArray{T,1}, r::UnitRange{<:Integer}) where T =
    TimeArray(timestamp(ta)[r], values(ta)[r], colnames(ta), meta(ta))

# array of rows
@propagate_inbounds getindex(ta::TimeArray, a::AbstractVector{<:Integer}) =
    TimeArray(timestamp(ta)[a], values(ta)[a, :], colnames(ta), meta(ta))

# array of 1d rows
@propagate_inbounds getindex(ta::TimeArray{T,1}, a::AbstractVector{<:Integer}) where T =
    TimeArray(timestamp(ta)[a], values(ta)[a], colnames(ta), meta(ta))

# single column by name
@propagate_inbounds function getindex(ta::TimeArray, s::Symbol)
    n = findcol(ta, s)
    TimeArray(timestamp(ta), values(ta)[:, n], Symbol[s], meta(ta), unchecked = true)
end

# array of columns by name
@propagate_inbounds getindex(ta::TimeArray, ss::Symbol...) = getindex(ta, collect(ss))
@propagate_inbounds getindex(ta::TimeArray, ss::Vector{Symbol}) =
    TimeArray(ta; values = values(ta)[:, map(s -> findcol(ta, s), ss)], colnames = ss)

# ta[rows, cols]
@propagate_inbounds getindex(ta::TimeArray,
         rows::Union{AbstractVector{<:Integer},Colon},
         cols::AbstractVector{Symbol}) =
    TimeArray(
        ta;
        timestamp = timestamp(ta)[rows],
        values = values(ta)[rows, map(s -> findcol(ta, s), cols)],
        colnames = cols,
        unchecked = true)

# ta[n, cols]
@propagate_inbounds getindex(ta::TimeArray, n::Integer, cols) =
    getindex(ta, [n], cols)

# ta[rows, col]
@propagate_inbounds getindex(ta::TimeArray, rows, col::Symbol) =
    getindex(ta, rows, [col])

# ta[n, col]
@propagate_inbounds getindex(ta::TimeArray, n::Integer, col::Symbol) =
    getindex(ta, [n], [col])

# single date
@propagate_inbounds function getindex(ta::TimeArray{T,N,D}, d::D) where {T,N,D}
    idxs = searchsorted(timestamp(ta), d)
    length(idxs) == 1 ? ta[idxs[1]] : nothing
end

# multiple dates
@propagate_inbounds function getindex(ta::TimeArray{T,N,D}, dates::Vector{D}) where {T,N,D}
    dates = sort(dates)
    idxs, _ = overlap(timestamp(ta), dates)
    ta[idxs]
end

# StepRange{Date,...}
@propagate_inbounds getindex(ta::TimeArray{T,N,D}, r::StepRange{D}) where {T,N,D} = ta[collect(r)]

@propagate_inbounds getindex(ta::TimeArray, k::TimeArray{Bool,1}) = ta[findwhen(k)]

# day of week
# getindex{T,N}(ta::TimeArray{T,N}, d::DAYOFWEEK) = ta[dayofweek(timestamp(ta)) .== d]

# Define end keyword
lastindex(ta::TimeArray, d::Integer = 1) =
    (d == 1) ? length(timestamp(ta)) :
    (d == 2) ? length(colnames(ta)) :
    1

eachindex(ta::TimeArray) = Base.OneTo(length(timestamp(ta)))

###### getproperty/propertynames #################

getproperty(ta::AbstractTimeSeries, c::Symbol) = ta[c]

propertynames(ta::TimeArray) = colnames(ta)

###### element wrapers ###########

"""
    timestamp(ta::TimeArray)

Get the time index of a `TimeArray`.
"""
timestamp(ta::TimeArray) = getfield(ta, :timestamp)

"""
    values(ta::TimeArray)

Get the underlying value table of a `TimeArray`.
"""
values(ta::TimeArray) = getfield(ta, :values)

"""
    colnames(ta::TimeArray)

Get the column names of a `TimeArray`.

# Examples

```julia-repl
julia> colnames(ohlc)
4-element Array{Symbol,1}:
 :Open
 :High
 :Low
 :Close
```
"""
colnames(ta::TimeArray) = getfield(ta, :colnames)

"""
    meta(ta::TimeArray)

Get the user-defined metadata of a `TimeArray`.
"""
meta(ta::TimeArray) = getfield(ta, :meta)

# internal use, to avoid name collision
_timestamp(ta::TimeArray) = getfield(ta, :timestamp)
_values(ta::TimeArray)    = getfield(ta, :values)
_colnames(ta::TimeArray)  = getfield(ta, :colnames)
_meta(ta::TimeArray)      = getfield(ta, :meta)


text:
@testset "Inverter Components" begin
    converter = AverageConverter(690.0, 2750000.0) #S_rated goes in Watts
    @test converter isa PowerSystems.DynamicComponent
    dc_source = FixedDCSource(600.0) #Not in the original data, guessed.
    @test dc_source isa PowerSystems.DynamicComponent
    filter = LCLFilter(0.08, 0.003, 0.074, 0.2, 0.01)
    @test filter isa PowerSystems.DynamicComponent
    pll = KauraPLL(500.0, 0.084, 4.69)
    @test pll isa PowerSystems.DynamicComponent
    virtual_H = VirtualInertia(2.0, 400.0, 20.0, 2 * pi * 50.0)
    @test virtual_H isa PowerSystems.DeviceParameter
    Q_control = ReactivePowerDroop(0.2, 1000.0)
    @test Q_control isa PowerSystems.DeviceParameter
    outer_control = OuterControl(virtual_H, Q_control)
    @test outer_control isa PowerSystems.DynamicComponent
    vsc = CurrentControl(0.59, 736.0, 0.0, 0.0, 0.2, 1.27, 14.3, 0.0, 50.0, 0.2)
    @test vsc isa PowerSystems.DynamicComponent
    BESS_source = ZeroOrderBESS(
        (sqrt(8) / sqrt(3)) * 690.0,
        (sqrt(3) / sqrt(8)) * 2750000.0,
        370.0,
        0.001,
        4.63,
        3200.0,
        0.6,
        4.0,
        0.39,
        10.34,
        1.08,
    )
    @test BESS_source isa PowerSystems.DynamicComponent
end

@testset "Dynamic Inverter" begin
    sys = create_system_with_dynamic_inverter()
    inverters = collect(get_components(DynamicInverter, sys))
    @test length(inverters) == 1
    test_inverter = inverters[1]
    @test test_inverter isa PowerSystems.Component
end


text:
module TestNA

using FactCheck
using MultidimensionalTables
using MultidimensionalTables: AbstractArrayWrapper, FloatNAArray, simplify_array, wrap_array, naop_plus

facts("NA tests") do
  @fact AbstractArrayWrapper(nalift([1,2,3])) --> AbstractArrayWrapper([Nullable(1), Nullable(2), Nullable(3)])
  @fact map(AbstractArrayWrapper,Any[nalift([1 2 3]),nalift([4,5])]) --> map(AbstractArrayWrapper,Any[[Nullable(1) Nullable(2) Nullable(3)],[Nullable(4), Nullable(5)]])
  @fact AbstractArrayWrapper(simplify_array(nalift(Any[[1 2 3],[4,5]]))) --> AbstractArrayWrapper(simplify_array(Any[Nullable([1 2 3]), Nullable([4,5])]))
  @fact @nalift([1 2 3;4 NA 5]) --> AbstractArrayWrapper([Nullable(1) Nullable(2) Nullable(3);Nullable(4) Nullable{Int}() Nullable(5)])
  @fact @nalift(Any[[1 2 NA],[3,NA]]) --> nalift(map(AbstractArrayWrapper,Any[[Nullable(1) Nullable(2) Nullable{Int}()], [Nullable(3), Nullable{Int}()]]))
  @fact AbstractArrayWrapper([x for x in FloatNAArray([1.0,2.0,NaN])]) --> AbstractArrayWrapper([Nullable(1.0), Nullable(2.0), Nullable{Float64}()])
  @fact AbstractArrayWrapper(FloatNAArray(reverse([1.0,2.0,3.0]))) --> AbstractArrayWrapper(FloatNAArray([3.0,2.0,1.0]))
  @fact simplify_array(nalift([1,2,3])) --> AbstractArrayWrapper([Nullable(1),Nullable(2),Nullable(3)])
  @fact simplify_array(Any[Nullable(3.5), Nullable(2.0), Nullable{Float64}()]) --> AbstractArrayWrapper(FloatNAArray([3.5, 2.0, NaN]))
  @fact (t=@larr(a=[1 2 3;4 5 6], b=[4 5 6;7 8 9]);tbltool.setna!(t.data.data[:a],2,2);tbltool.setna!(t,2,1);t) --> @larr(a=[1 2 3;NA NA 6], b=[4 5 6;NA 8 9])
  @fact igna(Nullable{Int}(), 3) --> 3
  @fact igna(Nullable(1), 3) --> 1
  @fact ignabool(nalift([true, false])) --> [true, false]
  @fact ignabool(Nullable{Bool}()) --> false
  @fact ignabool(Nullable(true)) --> true
  @fact ignabool(Nullable(false)) --> false
  @fact ignabool([Nullable(true), Nullable{Bool}()]) --> [true, false]
  @fact ignabool([Nullable(true)]) --> [true]
  @fact ignabool([Nullable(false)]) --> [false]
  @fact @nalift([NA NA]) --> wrap_array([Nullable{Any}() Nullable{Any}()])
  @fact @nalift([NA,NA]) --> wrap_array([Nullable{Any}(),Nullable{Any}()])
  @fact @nalift([NA NA;NA NA]) --> wrap_array(reshape(fill(Nullable{Any}(),4),2,2))
  @fact @nalift([1,NA]) --> wrap_array([Nullable(1), Nullable{Int}()])
  @fact @nalift([1 NA]) --> wrap_array([Nullable(1)  Nullable{Int}()])
  @fact @nalift([1 NA;NA NA]) --> wrap_array([Nullable(1)  Nullable{Int}();Nullable{Int}() Nullable{Int}()])
  @fact map(x->x, nalift([1.0,2.0,3.0])) --> nalift([1.0,2.0,3.0])
  @fact map((x,y)->naop_plus(x, y), nalift([1.0,2.0,3.0]), nalift([2.0,3.0,4.0])) --> nalift([3.0,5.0,7.0])
end

end


text:
module Purses

using Base: @_inline_meta

export Purse

include("purse.jl")
include("register.jl")
include("operators.jl")

end # module


text:
function mutualinfo(Z::Array{T,2};
                    remove_dups::Bool=true,
                    min_separation::Int=1,
                    pseudocount::Real=0.2,
                    theta=:auto,
                    output=:matrix # either :matrix, :score (i,j,val) or :apcscore
                    ) where T <: Integer

    q = Int(maximum(Z))
    Pi_true, Pij_true, Meff, _ = compute_weighted_frequencies(Z, :auto)
    N = div(length(Pi_true),q-1)
    if 0.0 < pseudocount <= 1.
        Pi, Pij = add_pseudocount(Pi_true, Pij_true, pseudocount)
        mi = mut_inf(Pij,Pi,q)
    elseif pseudocount == 0
        mi = mut_inf(Pij_true,Pi_true,q)
    else
        error("pseudocount = $pseudocount should be in [0,1]")
    end
    if output == :matrix
        return mi
    elseif output == :score
        Neff = N - min_separation
        score= Vector{Tuple{Int,Int,Float64}}(undef, (Neff*(Neff+1))>>1)

        ctr = 0
        for i=1:N-min_separation, j=i+min_separation:N
            ctr += 1
            score[ctr]=(i,j,mi[i,j])
        end
        sort!(score,by=x->x[3],rev=true)
        return score
    elseif output == :apcscore
        mi = correct_APC(mi)
        score = compute_ranking(mi,min_separation)
        return score
    else
        error("$output: output format not supported. Only :matrix, :score, :apcscore")
    end
end

function mutualinfo(filename::String;
                    max_gap_fraction::Real=0.9,
                    kwds...)
    Z = read_fasta_alignment(filename, max_gap_fraction)
    mutualinfo(Z; kwds...)
end

# This computation is complicated by the fact that in Pij we
# store only probabilities 1,...,q-1 × 1,...,q-1.
# We need therefore to recontruct the frame values (bottom, and right)
# A special case is for the element q,q

function mut_inf(Pij::Matrix{Float64},Pi::Vector{Float64},Q::Int)
    q=Q-1

    N1,N2 = size(Pij)
    N1 == N2 || error("matrix should be square. size(matrix) = $N1 != $N2")
    N3 = length((Pi))
    N3 == N2 || error("vector dim=$N3 size(matrix) = ($N1,$N1)")

    rem(N3,q) == 0 || error("dimension dim=$N3 should be multiple of $q")
    N = div(N1,q)

    scoreM = zeros(Float64,N,N)
    piq = zeros(Float64, N)    # vector of values for Pi[Q]
    pijq = zeros(Float64, N,q) # vector of values for Pij[i,j,Q,a] a in 1,...,q

    vec_row = zeros(Float64,q)
    vec_col = zeros(Float64,q)

    for i=1:N
        blk = (i-1)*q+1:i*q
        for j in blk
            piq[i] -= Pi[j]
        end
        piq[i] += 1.
        if 0.0 > piq[i] > 1.0
            error("probability is not normalized???")
        end
    end

    ctr = 0
    for i=1:N-1
        row0 = (i-1)*q
        for j=i+1:N
            compute_frame_col!(vec_col, i, j, Pij,Pi,q)
            compute_frame_row!(vec_row, i, j, Pij,Pi,q)
            col0 = (j-1)*q
            mi = 0.0
            for b=1:q
                for a=1:q
                    pij = Pij[row0 + a, col0 + b]
                    pi  = Pi[row0 + a]
                    pj  = Pi[col0 + b]
                    if pij > 0
                        mi += pij * log(pij / (pi*pj))

                    end
                end
            end
            for a=1:q
                pij = vec_row[a]
                pi  = piq[i]
                pj  = Pi[col0 + a]
                pipj = pi * pj
                if pij > 0 && pipj >0
                    mi += pij * log(pij / (pi*pj)) # contribution bottom frame
                end
                pij = vec_col[a]
                pi  = Pi[row0 + a]
                pj  = piq[j]
                pipj = pj * pj
                if pij > 0 && pipj > 0
                    mi += pij * log(pij / (pi*pj)) # contribution rigth frame
                end
            end

            _srow = sum(vec_row)  # contribution from Pij[q,q]
            _scol = sum(vec_col)
            pqqrow = piq[i] - _srow
            pqqcol = piq[j] - _scol
            pipj = piq[i]*piq[j]
            if pqqrow > 0 && pipj > 0
                mi += pqqrow * log(pqqrow/pipj)
            end
            scoreM[i,j] = mi
            scoreM[j,i] = mi
        end
    end
    return scoreM
end

function compute_frame_row!(vec_row::Vector{Float64},i::Int, j::Int, Pij::Matrix{Float64}, Pi::Vector{Float64}, q::Int)

    blki = (i-1)*q+1:i*q
    blkj = (j-1)*q+1:j*q
    ctr  = 0
    for b in blkj
        s = 0.0
        for a in blki
            s += Pij[a,b]
        end
        ctr += 1
        vec_row[ctr] = Pi[(j-1)*q + ctr] - s
        if 0.0 >= vec_row[ctr] >= 1.0
            println("element $ctr = ", vec_row[ctr])
            error("")
        end
    end
end

function compute_frame_col!(vec_col::Vector{Float64},i::Int, j::Int, Pij::Matrix{Float64}, Pi::Vector{Float64}, q::Int)

    blki = (i-1)*q+1:i*q
    blkj = (j-1)*q+1:j*q
    ctr  = 0
    for a in blki
        s = 0.0
        for b in blkj
            s += Pij[a,b]
        end
        ctr += 1
        vec_col[ctr] = Pi[(i-1)*q + ctr] - s
        if 0.0 >= vec_col[ctr] >= 1.0
            println("element $ctr = ", vec_col[ctr])
            error("");
        end
    end
end

max_stars_repo_name: giovact/PlmDCA

max_issues_repo_name: giovact/PlmDCA

max_forks_repo_name: giovact/PlmDCA


text:
using DashBootstrapComponents, DashHtmlComponents

cards = html_div([
    dbc_card(
        dbc_cardbody([
            html_h5("75% width card", className = "card-title"),
            html_p(
                [
                    "This card uses the ",
                    html_code("w-75"),
                    " class to set the width to 75%",
                ],
                className = "card-text",
            ),
        ]),
        className = "w-75 mb-3",
    ),
    dbc_card(
        dbc_cardbody([
            html_h5("50% width card", className = "card-title"),
            html_p(
                [
                    "This card uses the ",
                    html_code("w-50"),
                    " class to set the width to 50%",
                ],
                className = "card-text",
            ),
        ]),
        className = "w-50",
    ),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
using PerformanceIndex
using Test

error = [1., -2., 3]
data = [1., 1., 1.]
time = [1., 2., 3.]
function desired(t)
    t
end

@test iae(error) == 6.
@test prod(iae([error, error]) .== [6., 6.])
@test iae(data, time, desired) == 3.

@test ise(error) == 14.
@test prod(ise([error, error]) .== [14., 14.])
@test ise(data, time, desired) == 5.

@test itae(error, time) == 14.
@test prod(itae([error, error], time) .== [14., 14.])
@test itae(data, time, desired) == 8.

@test itse(error, time) == 36.
@test prod(itse([error, error], time) .== [36., 36.])
@test itse(data, time, desired) == 14.


text:
# # This is a title

# This is a description
# that spreads along
# multiple lines
# # [This is another title](@id custom_id)

x = 1+1


text:
if basename(pwd()) == "aoc"
    cd("2018/8")
end

abstract type Graph end

struct NodeTree <: Graph
    children::Vector{<:Graph}
    metadata::Vector{Int}
end

function NodeTree(xs::AbstractVector{<:Integer})
    childcount = popfirst!(xs)
    datacount = popfirst!(xs)
    children = [NodeTree(xs) for _ in 1:childcount]
    metadata = [popfirst!(xs) for _ in 1:datacount]
    NodeTree(children, metadata)
end

function Base.sum(tree::NodeTree; init::Real = 0)
    init + sum(tree.metadata, init = 0) + sum(sum(child) for child in tree.children; init = 0)
end

part1(tree::NodeTree) = sum(tree)

loadnumbers(filename::AbstractString) = parse.(Int, split(readline(filename)))
part1(filename::AbstractString) = part1(NodeTree(loadnumbers(filename)))

@assert part1("example.txt") == 138
part1("input.txt")

function value(tree::NodeTree)
    sum(
        if isempty(tree.children)
            tree.metadata
        else
            sum(
                if checkbounds(Bool, tree.children, index)
                    value(tree.children[index])
                else
                    0
                end
                for index in tree.metadata;
                init = 0
            )
        end,
        init = 0
    )
end

part2(tree::NodeTree) = value(tree)
part2(filename::AbstractString) = part2(NodeTree(loadnumbers(filename)))

@assert part2("example.txt") == 66

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using XGPaint
using Cosmology
using Test
using HDF5

println("")  # useful in Atom console

# relative background evolutions differ by 1e-3 between Julia and Python 2
rtol = 1e-3
cosmo = XGPaint.get_cosmology(h=0.7f0, OmegaM=0.25f0)

@testset "cib" begin
    # These tests come from computing examples in the Python xgpaint
    model = XGPaint.CIB_Planck2013{Float32}()
    r2z = XGPaint.build_r2z_interpolator(0.0f0, 4.5f0, cosmo)
    hod_shang = XGPaint.build_shang_interpolator(log(1.0f13), log(1.0f15), model)
    clnm2r = XGPaint.build_c_lnm2r_interpolator(nbin=30)
    sigma_sat = XGPaint.build_sigma_sat_ln_interpolator(log(4f15), model)
    muofn = XGPaint.build_muofn_interpolator(model)

    @test model.shang_Mpeak ≈ 10^12.3  # check default is viero as docs say

    # python: h2fm.utils.r2z(3407.0f0)
    @test isapprox(
        r2z(3407.0f0), 0.9996, rtol=rtol)

    # python: h2fm.utils.jiang_shmf(0.5e12, 1e13)
    @test XGPaint.jiang_shmf(0.5f12, 1.0f13, model) ≈ 2.5335820678141165

    # python: h2fm.hod.hod_shang(np.array([1e13, 1e14, 1e15]))
    @test hod_shang(log(1.0f13)) ≈ 9.78325852
    @test hod_shang(log(1.0f14)) ≈ 54.21666623
    @test hod_shang(log(1.0f15)) ≈ 336.55409152

    # python: h2fm.utils.mz2c(1e13, 4.0)
    @test isapprox(
        XGPaint.mz2c(1.0f13, 4.0f0, cosmo),
        2.262202312880881, rtol=rtol)

    # python: h2fm.utils.m2r(1e13)
    @test isapprox(
        XGPaint.m2r(1.0f13, cosmo),
        4.12328640425709, rtol=rtol)

    # python: h2fm.hod.clnm2r(np.array([1e-3, 25.0]), np.log(0.2))
    @test isapprox(
        clnm2r(1f-3 / 0.999f0, log(0.2f0)),
        0.44787404, rtol=rtol)
    @test isapprox(
        clnm2r(25.0f0 / 1.001f0, log(0.2f0)),
        0.085837, rtol=rtol)

    # python: h2fm.fluxmodel.sigma_cen(1e13)
    @test isapprox(XGPaint.sigma_cen(1.0f13, model),
        3218663770378.3066, rtol=rtol)

    # python: h2fm.fluxmodel.nu2theta(150e9, 0.5)
    @test isapprox(XGPaint.nu2theta(150f9, 0.5f0, model),
        6.5705729824589e-16, rtol=rtol)
    # python: h2fm.fluxmodel.nu2theta(1.5e6, 1.0)
    @test isapprox(XGPaint.nu2theta(1.5e6, 1.0, model),
        5.945023942373837e-34, rtol=rtol)

    # python: h2fm.fluxmodel.integrand_L(30, 32)
    @test isapprox(
        XGPaint.integrand_L(30.0f0, 32.0f0, model),
        3929486334377.691, rtol=rtol)

    # python: h2fm.fluxmodel.l2f(1.0, 1.0e-3, 1.0e-3, 1.0e-3)
    @test isapprox(
        XGPaint.l2f( 1.0f0, sqrt(3.0f-6), r2z(sqrt(3.0f-6))) * 4π,
        333198.42738065525, rtol=rtol)

    # python: h2fm.fluxmodel.sigma_sat(np.array([1e13, 4e15]))
    @test sigma_sat(log(4f15)) ≈ 3.50530158e+14
    @test sigma_sat(log(1f13)) ≈ 2.48648819e+12

    @test muofn(500.0f0) ≈ 6.14975653e-05
    @test muofn(1.0f0) ≈ 0.11765558

    @test XGPaint.shang_z_evo(0.0f0, model) ≈ 1.0f0
end

radio_model = Radio_Sehgal2009{Float32}()
@testset "radio" begin
    @test XGPaint.FR_I_redshift_evolution(0.0f0, radio_model) ≈ 1.0f0
    @test XGPaint.FR_II_redshift_evolution(0.0f0, radio_model) ≈ 1.0f0
    @test (XGPaint.FR_I_redshift_evolution(radio_model.I_z_p, radio_model) ≈
        XGPaint.FR_I_redshift_evolution(radio_model.I_z_p+1.0f0, radio_model))
    @test XGPaint.FR_II_redshift_evolution(1.3f0, radio_model) > 190

end

@testset "sanity" begin
    @test XGPaint.l2f(1.0, 1.0, 0.0) ≈ 1/4π
    @test XGPaint.random_phi(Float32) <= 2π
    @test XGPaint.random_theta(Float32) <= π
    @test XGPaint.random_phi(Float64) <= 2π
    @test XGPaint.random_theta(Float64) <= π

    # file reading test
    testData = rand( Float32, 4, 100 )
    h5open("test_file_writing.h5", "w") do file
        write(file, "halos", testData)  # alternatively, say "@write file A"
    end
    read_pos, read_mass = XGPaint.read_halo_catalog_hdf5("test_file_writing.h5")
    @test all( read_pos .≈ testData[1:3,:] )
    @test all( read_mass .≈ testData[4,:] )
    rm("test_file_writing.h5")

    @test all(XGPaint.chunk(10, 3) == [(1,3), (4,6), (7,9), (10,10)])
    @test all(XGPaint.chunk(10, 4) == [(1,4), (5,8), (9,10)])
    @test all(XGPaint.chunk(10, 10) == [(1,10)])

end


text:
function coupled_henon1D(M,N, u0=rand(M), v0=rand(M))
    function henon(U,V)
        Un = copy(U)
        Vn = copy(V)
        Un[1] = Un[M] = 0.5
        Vn[1] = Vn[M] = 0
        for m=2:M-1
            Un[m] = 1 - 1.45*(.5*U[m]+ .25*U[m-1] + .25U[m+1])^2 + 0.3*V[m]
            Vn[m] = U[m]
        end
        return Un, Vn
    end

    U = Vector{Vector{Float64}}(undef,N)
    V = Vector{Vector{Float64}}(undef,N)
    U[1] = u0
    V[1] = v0
    for n = 2:N
        U[n],V[n] = henon(U[n-1],V[n-1])
    end
    return U,V
end

function coupled_henon2D(X,Y,N,u0=rand(X,Y), v0=rand(X,Y))
    function henon(U,V)
        Un = copy(U)
        Vn = copy(V)
        Un[1,:] = Un[:,1] = Un[X, :] = Un[:, Y] .= 0.5
        Vn[1,:] = Vn[:,1] = Vn[X, :] = Vn[:, Y] .= 0
        for mx=2:X-1 , my=2:Y-1
            Un[mx,my] = 1 - 1.45*(.5*U[mx,my]+ .125*U[mx-1,my] +
             .125U[mx+1,my]+ .125U[mx, my-1]+ .125U[mx,my+1])^2 + 0.3*V[mx,my]
            Vn[mx,my] = U[mx,my]
        end
        return Un, Vn
    end
    U = Vector{Matrix{Float64}}(undef,N)
    V = Vector{Matrix{Float64}}(undef,N)
    U[1] = u0
    V[1] = v0
    for n = 2:N
        U[n],V[n] = henon(U[n-1],V[n-1])
    end
    return U,V

end

"""
```
barkley(T;
        tskip=0,
        periodic=true,
        ssize=(50,50),
        a=0.75, b=0.06, ε=0.08, D=1/50, h=0.1, Δt=0.1)
```
Simulate the Barkley model (nonlinear `u^3` term).
"""
function barkley(T;
                tskip=0,
                periodic=true,
                ssize=(50,50),
                a=0.75, b=0.06, ε=0.08, D=1/50, h=0.1, Δt=0.1,
                init = rand(10,10))

    Nx, Ny = ssize
    @assert Nx ≥ 40
    @assert Ny ≥ 40
    U = Vector{Array{Float64,2}}()
    V = Vector{Array{Float64,2}}()

    u = zeros(Nx,Ny)
    v = zeros(Nx,Ny)

    αα, ββ = size(init)

    v .= u .= repeat(init, inner=(Nx÷αα,Ny÷ββ))
    v .= v.>0.2
    Σ = zeros(Nx, Ny, 2)
    r,s = 1,2

    function F(u, uth)
        if u < uth
            u/(1-(Δt/ε)*(1-u)*(u-uth))
        else
            (u + (Δt/ε)*u*(u-uth))/(1+(Δt/ε)*u*(u-uth))
        end
    end

    function periodic_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,Δt,h)
        for j=1:Ny, i=1:Nx
            uth = (v[i,j] + b)/a
            v[i,j] = v[i,j] + Δt*(u[i,j]^3 - v[i,j])
            u[i,j] = F(u[i,j], uth) + D * Δt/h^2 *Σ[i,j,r]
            Σ[i,j,s] -= 4u[i,j]
            Σ[  mod(i-1-1,Nx)+1,j,s] += u[i,j]
            Σ[  mod(i+1-1,Nx)+1,j,s] += u[i,j]
            Σ[i,mod(j-1-1,Ny)+1,  s] += u[i,j]
            Σ[i,mod(j+1-1,Ny)+1,  s] += u[i,j]
            Σ[i,j,r] = 0
        end
    end
    function constant_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,Δt,h)
        for i=1:Nx, j=1:Ny
            uth = (v[i,j] + b)/a
            v[i,j] = v[i,j] + Δt*(u[i,j] - v[i,j])
            u[i,j] = F(u[i,j], uth) + Δt/h^2 *Σ[i,j,r]
            Σ[i,j,s] -= 4u[i,j]
            i > 1  && (Σ[i-1,j,s] += u[i,j])
            i < Nx && (Σ[i+1,j,s] += u[i,j])
            j > 1  && (Σ[i,j-1,s] += u[i,j])
            j < Ny && (Σ[i,j+1,s] += u[i,j])
            Σ[i,j,r] = 0
        end
    end

    for m=1:T+tskip
        periodic ? periodic_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,Δt,h) :
                   constant_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,Δt,h)
        r,s = s,r
        if m > tskip
            push!(U,copy(u))
            push!(V,copy(v))
        end
    end
    return U,V
end


text:
include("euler/euler.jl")
using .Primes: prime_sieve

function compute(n::Int)::Int
    sieve = prime_sieve(n)
    candidates = [i - 1 for (i, is_prime) ∈ enumerate(sieve) if is_prime && (i - 1) % 4 != 0]
    result = 0
    for candidate ∈ candidates
        if !sieve[candidate ÷ 2 + 2] continue end
        is_valid = true
        for i ∈ 3:trunc(Int, √n + 1)
            if candidate % i ≠ 0
                continue
            elseif !sieve[candidate ÷ i + i]
                is_valid = false
                break
            end
        end
        if is_valid result += candidate end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler