text:
abstract type BackgroundExtractor end

"""
    update!(be::BackgroundExtractor, img)

Update a background model using `img`.
"""
LowLevelParticleFilters.update!(be::BackgroundExtractor, img) = push!(be.fb, img)

"""
    MedianBackground{T} <: BackgroundExtractor

Models the background of a sequence of images as the median over a fixed time window. Supports the same constructors as a `FrameBuffer`

This type supports `background, foreground, update!`
"""
struct MedianBackground{T} <: BackgroundExtractor
    fb::FrameBuffer{T}
end
MedianBackground(args...) = MedianBackground(FrameBuffer(args...))

"""
    background(be::MedianBackground)

Extract the background of `img`.
"""
function background(be::MedianBackground, img=nothing)
    median(be.fb)
end

"""
    foreground(be::MedianBackground, img)

Extract the foreground of `img`.
"""
function foreground(be::MedianBackground, img)
    img - background(be, img)
end

"""
    DiffBackground{T} <: BackgroundExtractor

Models the background of a sequence of images as the diff between two consequtive images.

This type supports `background, foreground, update!`
"""
struct DiffBackground{T} <: BackgroundExtractor
    fb::FrameBuffer{T}
end
DiffBackground(img) = DiffBackground(FrameBuffer(img, 2))

"""
    background(be::DiffBackground, img)

Extract the background of `img`.
"""
function background(be::DiffBackground, img)
    dropdims(diff(be.fb).b, dims=3) .* img
end

"""
    foreground(be::DiffBackground, img)

Extract the foreground of `img`.
"""
function foreground(be::DiffBackground, img)
    dropdims(complement.(diff(be.fb).b), dims=3) .* img
end

# struct PCABackground{T} <: BackgroundExtractor
#     k::Int
#     fb::FrameBuffer{T}
# end
# PCABackground(k::Int, args...) = PCABackground(k, FrameBuffer(args...))
#
# function background(be::PCABackground, img)
#     pca = fit(PCA, Matrix(be.fb)'; maxoutdim=be.k)
#     r = reconstruct(pca, vec(img))
#     reshape(r, imgsize(be.fb))
# end
#
# function foreground(be::PCABackground, img)
#     img - background(be, img)
# end


text:
# set up McMC using 4 chains
using Distributed
nsamples, nchains, nchainsatone = 20001, 4, 1
Tmax = 2.50
addprocs(nchains)
@info "workers are $(workers())"
@everywhere using Distributed
@everywhere using transD_GP
## run McMC
@time transD_GP.main(opt, FOp, Tmax=Tmax, nsamples=nsamples, nchains=nchains, nchainsatone=nchainsatone)
rmprocs(workers())
## plot
transD_GP.getchi2forall(opt, fsize=8, alpha=0.5)
transD_GP.plot_posterior(FOp, opt,  burninfrac=0.5, figsize=(8.5,4), fsize=8, nbins=50)

ext: jl


text:
import Base: zero

abstract type MeanFunction end

"""
    ZeroMean{T<:Real} <: MeanFunction

Returns `zero(T)` everywhere.
"""
struct ZeroMean{T<:Real} <: MeanFunction end
ZeroMean() = ZeroMean{Float64}()
ew(::ZeroMean{T}, x::AV) where T = zeros(T, length(x))
zero(::MeanFunction) = ZeroMean()

"""
    OneMean{T} <: MeanFunction

Return `one(T)` everywhere.
"""
struct OneMean{T<:Real} <: MeanFunction end
OneMean() = OneMean{Float64}()
ew(::OneMean{T}, x::AV) where T = ones(T, length(x))

"""
    ConstMean{T} <: MeanFunction

Returns `c` everywhere.
"""
struct ConstMean{T<:Real} <: MeanFunction
    c::T
end
ew(m::ConstMean, x::AV) = fill(m.c, length(x))

"""
    CustomMean{Tf} <: MeanFunction

A wrapper around whatever unary function you fancy.
"""
struct CustomMean{Tf} <: MeanFunction
    f::Tf
end
@noinline ew(f::CustomMean, x::AV) = map(f.f, x)


text:
function neighbor!(bool::BoolParameter)
    perturb!(bool)
end

function neighbor!(number::NumberParameter;
                   interval::Number = typeof(number.value) <: Integer ?
                                      rand_in(1,
                                              number.max - number.min) :
                                      rand_in(1.,
                                              number.max - number.min),
                   distance::Int = 1)
    for i = 1:distance
        previous = number.value
        while number.value == previous
            perturb!(number, interval)
        end
    end
    number
end

function neighbor!(enum::EnumParameter; distance::Int = 1)
    for i = 1:distance
        previous = enum.value
        while enum.value == previous
            perturb!(enum)
        end
    end
    enum
end

function neighbor!(permutation::PermutationParameter;
                   interval = rand(1:permutation.size))
    perturb!(permutation, interval)
end

function neighbor!(configuration::Configuration,
                   intervals::Dict{String, Any}; distance::Int = 1)
    key_set = collect(keys(intervals))
    target  = key_set[rand(1:length(key_set))]
    for i = 1:distance
        neighbor!(configuration[target], distance = intervals[target])
    end
    update!(configuration)
    configuration
end

function neighbor!(configuration::Configuration; distance::Int = 1)
    key_set = collect(keys(configuration.parameters))
    target  = key_set[rand(1:length(key_set))]
    for i = 1:distance
        if !(typeof(configuration[target]) <: StringParameter)
            neighbor!(configuration[target])
        end
    end
    update!(configuration)
    configuration
end

function neighbor!(configuration::Configuration, target::String;
                   distance::Int = 1)
    for i = 1:distance
        if !(typeof(configuration[target]) <: StringParameter)
            neighbor!(configuration[target])
        end
    end
    update!(configuration)
    configuration
end


text:
export add_fixture, apply_fixtures, add_simple_fixture

SymbolOrNothing = Union(Symbol, Nothing)

immutable NamedFixture
  name::SymbolOrNothing
  fn::Function
  args::Any
  kwargs::Vector
end

fixtures = Dict{Symbol, Array{NamedFixture}}()

get_name(nfs::Array{NamedFixture}) = Set( [nf.name for nf in nfs] )

function is_name_good(name::SymbolOrNothing)
  result = false

  if name != nothing
    for fixtures in values(Fixtures.fixtures)
      if name in get_name(fixtures)
        result = true
      end
    end
  end

  return result
end

function add_fixture(scope::Symbol, name::SymbolOrNothing, fixture::Function, args...; kwargs...)

  if is_name_good(name)
    error("A fixture named $name is already defined.")
  end

  const nf = NamedFixture(name, fixture, args, kwargs)

  if haskey(fixtures, scope)
    append!(fixtures[scope], [nf])
  else
    fixtures[scope] = Array(NamedFixture, 1)
    fixtures[scope][1] = nf
  end
  return
end

add_fixture(scope::Symbol, fixture::Function, args...; kwargs...) = add_fixture(scope, nothing, fixture, args...; kwargs...)

function add_simple_fixture(scope::Symbol, name::SymbolOrNothing, before::Function, after::Function)
  add_fixture(scope, @fixture function()
                const val = before()
                yield_fixture(val)
                after()
              end)
end

add_simple_fixture(scope::Symbol, before::Function, after::Function) = add_simple_fixture(scope, nothing, before, after)

list_map(args...) = [ map(args...)... ]

@fixture function apply_fixtures(scope::Symbol; fixture_values::Bool=false)
  # Convert all the fixtures to tasks

  const fv = Dict{Symbol, Any}()

  temp = get(fixtures, scope, ())

  const tsks = list_map(temp) do nf
    if fixture_values && nf.name!=nothing
      return Task() do
        nf.fn(nf.args..., nf.kwargs...) do values::Any...
          if length(values)==0
            fv[nf.name] = nothing
          elseif length(values)==1
            fv[nf.name] = values[1]
          else
            fv[nf.name] = values
          end
          produce()
        end
      end
    else
      return Task(() -> nf.fn(produce, nf.args..., nf.kwargs...))
    end
  end

  # fixtures can be defined within in nested scopes, create a new scope
  const old = fixtures
  global fixtures = copy(fixtures)

  # Begin all the tasks
  for tsk in tsks
    consume(tsk)
  end

  yield_fixture((fixture_values ? (fv,) : ())...)

  # Complete the tasks in the opposite order we began them
  reverse!(tsks)
  for tsk in tsks
    consume(tsk)
  end

  # Reset our scope
  global fixtures = old
end


text:
#=
Should not be constructed by the user - acts like a vector (used in eigsolve)
I think it makes sense to see these things as an actual state instead of return an array of B tensors (what we used to do)
This will allow us to plot energy density (finite qp) and measure observeables.
=#

struct LeftGaugedQP{S,T1,T2}
    # !(left_gs === right_gs) => domain wall excitation
    left_gs::S
    right_gs::S

    VLs::Vector{T1} # AL' VL = 0 (and VL*X = B)
    Xs::Vector{T2} # contains variational parameters

    momentum::Float64
end

struct RightGaugedQP{S,T1,T2}
    # !(left_gs === right_gs) => domain wall excitation
    left_gs::S
    right_gs::S

    Xs::Vector{T2}
    VRs::Vector{T1}

    momentum::Float64
end

#constructors
function LeftGaugedQP(datfun,left_gs,right_gs=left_gs;sector = first(sectors(oneunit(virtualspace(left_gs,1)))),momentum=0.0)
    #find the left null spaces for the TNS
    excitation_space = ℂ[typeof(sector)](sector => 1);
    VLs = [adjoint(rightnull(adjoint(v))) for v in left_gs.AL]
    Xs = [TensorMap(datfun,eltype(left_gs.AL[1]),_lastspace(VLs[loc])',excitation_space'*virtualspace(right_gs,loc)) for loc in 1:length(left_gs)]
    left_gs isa InfiniteMPS || momentum == zero(momentum) || @warn "momentum is ignored for finite quasiparticles"
    LeftGaugedQP(left_gs,right_gs,VLs,Xs,momentum)
end

function RightGaugedQP(datfun,left_gs,right_gs=left_gs;sector = first(sectors(oneunit(virtualspace(left_gs,1)))),momentum=0.0)
    #find the left null spaces for the TNS
    excitation_space = ℂ[typeof(sector)](sector => 1);
    VRs = [adjoint(leftnull(adjoint(v))) for v in _transpose_tail.(left_gs.AR)]
    Xs = [TensorMap(datfun,eltype(left_gs.AL[1]),virtualspace(right_gs,loc-1)',excitation_space'*_firstspace(VRs[loc])) for loc in 1:length(left_gs)]
    left_gs isa InfiniteMPS || momentum == zero(momentum) || @warn "momentum is ignored for finite quasiparticles"
    RightGaugedQP(left_gs,right_gs,Xs,VRs,momentum)
end

#gauge dependent code
Base.similar(v::LeftGaugedQP,t=eltype(v)) = LeftGaugedQP(v.left_gs,v.right_gs,v.VLs,map(e->similar(e,t),v.Xs),v.momentum)
Base.similar(v::RightGaugedQP,t=eltype(v)) = RightGaugedQP(v.left_gs,v.right_gs,map(e->similar(e,t),v.Xs),v.VRs,v.momentum)

Base.getindex(v::LeftGaugedQP,i::Int) = v.VLs[mod1(i,end)]*v.Xs[mod1(i,end)];
Base.getindex(v::RightGaugedQP,i::Int) = @plansor t[-1 -2;-3 -4] := v.Xs[mod1(i,end)][-1;-3 1]*v.VRs[mod1(i,end)][1;-4 -2];

function Base.setindex!(v::LeftGaugedQP,B,i::Int)
    v.Xs[mod1(i,end)] = v.VLs[mod1(i,end)]'*B
    v
end
function Base.setindex!(v::RightGaugedQP,B,i::Int)
    @plansor v.Xs[mod1(i,end)][-1; -2 -3]:=B[-1 1;-2 2]*conj(v.VRs[mod1(i,end)][-3;2 1])
    v
end

#conversion between gauges (partially implemented)
function Base.convert(::Type{RightGaugedQP},input::LeftGaugedQP{S}) where S<:InfiniteMPS
    rg = RightGaugedQP(zeros,input.left_gs,input.right_gs,sector = first(sectors(utilleg(input))), momentum = input.momentum);
    len = length(input);

    #construct environments
    rBs = [@plansor t[-1;-2 -3] := input[len][-1 2;-2 3]*conj(input.right_gs.AR[len][-3 2;3])*exp(1im*input.momentum)]
    for i in len-1:-1:1
        t = exci_transfer_right(rBs[end],input.left_gs.AL[i],input.right_gs.AR[i]);
        @plansor t[-1;-2 -3] += input[i][-1 2;-2 3]*conj(input.right_gs.AR[i][-3 2;3])
        push!(rBs,exp(1im*input.momentum)*t);
    end
    rBs = reverse(rBs);

    (rBE,convhist) = @closure linsolve(rBs[1],rBs[1],GMRES()) do x
        y = transfer_right(x,input.left_gs.AL,input.right_gs.AR)*exp(1im*input.momentum*len)
        if input.trivial
            @plansor y[-1;-2 -3]-=y[1;-2 2]*l_LR(input.right_gs)[2;1]*r_LR(input.right_gs)[-1;-3]
        end
        x-y
    end
    convhist.converged == 0 && @warn "failed to converge $(convhist.normres)"

    rBs[1] = rBE;
    for i in len:-1:2
        rBE = transfer_right(rBE,input.left_gs.AL[i],input.right_gs.AR[i])*exp(1im*input.momentum);
        rBs[i] += rBE;
    end

    #final contraction is now easy
    for i in 1:len
        @plansor T[-1 -2;-3 -4] := input.left_gs.AL[i][-1 -2;1]*rBs[mod1(i+1,end)][1;-3 -4]
        @plansor T[-1 -2;-3 -4] += input[i][-1 -2;-3 -4]
        rg[i] = T
    end

    rg
end
function Base.convert(::Type{LeftGaugedQP},input::RightGaugedQP{S}) where S<:InfiniteMPS
    lg = LeftGaugedQP(zeros,input.left_gs,input.right_gs,sector = first(sectors(utilleg(input))), momentum = input.momentum);
    len = length(input);

    lBs = [@plansor t[-1;-2 -3] := input[1][1 2;-2 -3]*conj(input.left_gs.AL[1][1 2;-1])*exp(-1im*input.momentum)];
    for i in 2:len
        t = exci_transfer_left(lBs[end],input.right_gs.AR[i],input.left_gs.AL[i]);
        @plansor t[-1;-2 -3] += input[i][1 2;-2 -3]*conj(input.left_gs.AL[i][1 2;-1])
        push!(lBs,t*exp(-1im*input.momentum));
    end

    (lBE,convhist) = @closure linsolve(lBs[end],lBs[end],GMRES()) do x
        y = transfer_left(x,input.right_gs.AR,input.left_gs.AL)*exp(-1im*input.momentum*len)
        if input.trivial
            @plansor y[-1;-2 -3] -= y[1;-2 2]*r_RL(input.right_gs)[2;1]*l_RL(input.right_gs)[-1;-3]
        end
        x-y
    end
    convhist.converged == 0 && @warn "failed to converge $(convhist.normres)"

    lBs[end] = lBE;
    for i in 1:len-1
        lBE = transfer_left(lBE,input.right_gs.AR[i],input.left_gs.AL[i])*exp(-1im*input.momentum)
        lBs[i]+=lBE;
    end

    for i in 1:len
        @plansor T[-1 -2;-3 -4] := lBs[mod1(i-1,len)][-1;-3 1]*input.right_gs.AR[i][1 -2;-4]
        @plansor T[-1 -2;-3 -4] += input[i][-1 -2;-3 -4]
        lg[i] = T
    end

    lg
end

# gauge independent code
const QP{S,T1,T2} = Union{LeftGaugedQP{S,T1,T2},RightGaugedQP{S,T1,T2}} where {S,T1,T2};
const FiniteQP{S,T1,T2} = QP{S,T1,T2} where {S<:FiniteMPS}
const InfiniteQP{S,T1,T2} = QP{S,T1,T2} where {S<:InfiniteMPS}

utilleg(v::QP) = space(v.Xs[1],2)
Base.copy(a::QP) = copy!(similar(a),a)
Base.copyto!(a::QP,b::QP) = copy!(a,b);
function Base.copy!(a::T,b::T) where T<:QP
    for (i,j) in zip(a.Xs,b.Xs)
        copy!(i,j)
    end
    a
end
function Base.getproperty(v::QP,s::Symbol)
    if s == :trivial
        return v.left_gs === v.right_gs
    else
        return getfield(v,s)
    end
end

function Base.:-(v::T,w::T) where T<:QP
    t = similar(v)
    t.Xs[:] = (v.Xs-w.Xs)[:]
    t
end
function Base.:+(v::T,w::T) where T<:QP
    t = similar(v)
    t.Xs[:] = (v.Xs+w.Xs)[:]
    t
end
LinearAlgebra.dot(v::T, w::T)  where T<:QP = sum(dot.(v.Xs, w.Xs))
LinearAlgebra.norm(v::QP) = norm(norm.(v.Xs))
LinearAlgebra.normalize!(w::QP) = rmul!(w,1/norm(w));
Base.length(v::QP) = length(v.Xs)
Base.eltype(v::QP) = eltype(eltype(v.Xs)) # - again debateable, need scaltype

function LinearAlgebra.mul!(w::T, a, v::T)  where T<:QP
    @inbounds for (i,j) in zip(w.Xs,v.Xs)
        LinearAlgebra.mul!(i, a, j)
    end
    return w
end

function LinearAlgebra.mul!(w::T, v::T, a)  where T<:QP
    @inbounds for (i,j) in zip(w.Xs,v.Xs)
        LinearAlgebra.mul!(i, j, a)
    end
    return w
end
function LinearAlgebra.rmul!(v::QP, a)
    for x in v.Xs
        LinearAlgebra.rmul!(x, a)
    end
    return v
end

function LinearAlgebra.axpy!(a, v::T, w::T)  where T<:QP
    @inbounds for (i,j) in zip(w.Xs,v.Xs)
        LinearAlgebra.axpy!(a, j, i)
    end
    return w
end
function LinearAlgebra.axpby!(a, v::T, b, w::T)  where T<:QP
    @inbounds for (i,j) in zip(w.Xs,v.Xs)
        LinearAlgebra.axpby!(a, j, b, i)
    end
    return w
end

Base.:*(v::QP, a) = mul!(similar(v),a,v)
Base.:*(a, v::QP) = mul!(similar(v),a,v)

Base.zero(v::QP) = v*0;

function Base.convert(::Type{<:FiniteMPS},v::QP{S}) where S <: FiniteMPS
    #very slow and clunky, but shouldn't be performance critical anyway

    elt = eltype(v)

    utl = utilleg(v); ou = oneunit(utl); utsp = ou ⊕ ou;
    upper = isometry(Matrix{elt},utsp,ou); lower = leftnull(upper);
    upper_I = upper*upper'; lower_I = lower*lower'; uplow_I = upper*lower';

    Ls = v.left_gs.AL[1:end];
    Rs = v.right_gs.AR[1:end];

    #step 0 : fuse the utility leg of B with the first leg of B
    orig_Bs = map(i->v[i],1:length(v))
    Bs = @closure map(orig_Bs) do t
        frontmap = isomorphism(storagetype(t),fuse(utl*_firstspace(t)),utl*_firstspace(t));
        @plansor tt[-1 -2;-3]:=t[1 -2;2 -3]*frontmap[-1;2 1]
    end

    function simplefuse(temp)
        frontmap = isomorphism(storagetype(temp),fuse(space(temp,1)*space(temp,2)),space(temp,1)*space(temp,2))
        backmap = isomorphism(storagetype(temp),space(temp,5)'*space(temp,4)',fuse(space(temp,5)'*space(temp,4)'))

        @plansor tempp[-1 -2;-3] := frontmap[-1;1 2]*temp[1 2 -2 3;4]*backmap[4 3;-3]
    end

    #step 1 : pass utl through Ls
    passer = isomorphism(Matrix{elt},utl,utl);
    for (i,L) in enumerate(Ls)
        @plansor temp[-1 -2 -3 -4;-5]:=L[-2 -3;-4]*passer[-1;-5]
        Ls[i] = simplefuse(temp)
    end

    #step 2 : embed all Ls/Bs/Rs in the same space
    superspaces = map(zip(Ls,Rs)) do (L,R)
        supremum(space(L,1),space(R,1))
    end
    push!(superspaces,supremum(_lastspace(Ls[end])',_lastspace(Rs[end])'))

    for i in 1:(length(v)+1)
        Lf = isometry(Matrix{elt},superspaces[i],i <= length(v) ? _firstspace(Ls[i]) : _lastspace(Ls[i-1])')
        Rf = isometry(Matrix{elt},superspaces[i],i <= length(v) ? _firstspace(Rs[i]) : _lastspace(Rs[i-1])')

        if i <= length(v)
            @plansor Ls[i][-1 -2;-3] := Lf[-1;1]*Ls[i][1 -2;-3]
            @plansor Rs[i][-1 -2;-3] := Rf[-1;1]*Rs[i][1 -2;-3]
            @plansor Bs[i][-1 -2;-3] := Lf[-1;1]*Bs[i][1 -2;-3]
        end

        if i>1
            @plansor Ls[i-1][-1 -2;-3] := Ls[i-1][-1 -2;1]*conj(Lf[-3;1])
            @plansor Rs[i-1][-1 -2;-3] := Rs[i-1][-1 -2;1]*conj(Rf[-3;1])
            @plansor Bs[i-1][-1 -2;-3] := Bs[i-1][-1 -2;1]*conj(Rf[-3;1])
        end
    end

    #step 3 : fuse the correct *_I with the correct tensor (and enforce boundary conditions)
    function doboundary(temp1,pos)
        if pos == 1
            @plansor temp2[-1 -2 -3 -4;-5] := temp1[1 -2 -3 -4;-5]*conj(upper[1;-1])
        elseif pos == length(v)
            @plansor temp2[-1 -2 -3 -4;-5] := temp1[-1 -2 -3 -4;1]*lower[1;-5]
        else
            temp2 = temp1;
        end

        temp2
    end

    for i in 1:length(v)
        @plansor temp[-1 -2 -3 -4; -5] := Ls[i][-2 -3;-4]*upper_I[-1;-5]
        temp = doboundary(temp,i);
        Ls[i] = simplefuse(temp) * (i<length(v));

        @plansor temp[-1 -2 -3 -4; -5] := Rs[i][-2 -3;-4]*lower_I[-1;-5]
        temp = doboundary(temp,i);
        Rs[i] = simplefuse(temp) * (i>1);

        @plansor temp[-1 -2 -3 -4; -5] := Bs[i][-2 -3;-4]*uplow_I[-1;-5]
        temp = doboundary(temp,i);
        Bs[i] = simplefuse(temp);
    end

    return FiniteMPS(Ls+Rs+Bs,normalize=false)
end

function Base.getproperty(exci::Multiline{<:InfiniteQP},s::Symbol)
    if s == :momentum
        return first(exci.data).momentum;
    elseif s == :left_gs
        Multiline(map(x->x.left_gs,exci.data))
    elseif s == :right_gs
        Multiline(map(x->x.right_gs,exci.data))
    elseif s == :trivial
        return reduce(&,map(x->x.trivial,exci.data),init=true);
    else
        return getfield(exci,s);
    end
end


text:
using LBFGSB
using Test

@testset "driver3" begin
# translated from driver3.f
# "This time-controlled driver shows that it is possible to terminate a run by
#  elapsed CPU time, and yet be able to print all desired information. This driver
#  also illustrates the use of two stopping criteria that may be used in conjunction
#  with a limit on execution time. The sample problem used here is the same as in
#  driver1 and driver2 (the extended Rosenbrock function with bounds on the variables).

nmax = 1024   # the dimension of the largest problem to be solved
mmax = 17    # the maximum number of limited memory corrections

task = fill(Cuchar(' '), 60)    # fortran's blank padding
csave = fill(Cuchar(' '), 60)    # fortran's blank padding
lsave = zeros(Cint, 4)
isave = zeros(Cint, 44)
dsave = zeros(Cdouble, 29)

# "We specify a limite on the CPU time (in seconds)."
tlimit = 0.2

# "We suppress the default output.
#  (The user could also elect to use the default output by choosing iprint >= 0.)"
iprint = -1

# "f is a DOUBLE PRECISION variable.  If the routine setulb returns with task(1:2)= 'FG',
#  then f must be set by the user to contain the value of the function at the point x."
f = 0.0

# "g is a DOUBLE PRECISION array of length n.  If the routine setulb returns with taskb(1:2)= 'FG',
#  then g must be set by the user to contain the components of the gradient at the point x."
g = zeros(Cdouble, nmax)

wa = zeros(Cdouble, 2mmax*nmax + 5nmax + 11mmax*mmax + 8mmax)
iwa = zeros(Cint, 3*nmax)

# "We suppress both code-supplied stopping tests because we will provide our own termination conditions"
factr = 0.0
pgtol = 0.0

n = 1000    # the dimension n of the sample problem
m = 10      # the number of m of limited memory correction stored

# provide nbd which defines the bounds on the variables:
nbd = zeros(Cint, nmax)
l = zeros(Cdouble, nmax)    # the lower bounds
u = zeros(Cdouble, nmax)    # the upper bounds
# "First set bounds on the odd-numbered variables."
for i = 1:2:n
    nbd[i] = 2
    l[i] = 1e0
    u[i] = 1e2
end
# "Next set bounds on the even-numbered variables."
for i = 2:2:n
    nbd[i] = 2
    l[i] = -1e2
    u[i] = 1e2
end

# "We now define the starting point."
x = zeros(Cdouble, nmax)
for i = 1:n
    x[i] = 3e0
end

println("Solving sample problem.")
println(" (f = 0.0 at the optimal solution.)")

# "We start the iteration by initializing task."
task[1:5] = b"START"

## ------- the beginning of the loop ----------
time1Ref = Ref{Cdouble}(0)
time2Ref = Ref{Cdouble}(0)

timer(time1Ref)

let
    @label CALLLBFGSB

    # This is the call to the L-BFGS-B code.
    setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave)

    # "the minimization routine has returned to request the function f and gradient g values at the current x."
    if task[1:2] == b"FG"
        # "Before evaluating f and g we check the CPU time spent."
        timer(time2Ref)

        if time2Ref[] - time1Ref[] > tlimit
            task[1:4] = b"STOP"
            # "In this driver we have chosen to disable the printing options of the
            #  code (we set iprint=-1); instead we are using customized output: we
            #  print the latest value of x, the corresponding function value f and
            #  the norm of the projected gradient |proj g|."

            println("STOP: CPU EXCEEDING THE TIME LIMIT.")

            # "We print the latest iterate contained in wa(j+1:j+n), where j = 3*n+2*m*n+11*m**2"
            j = 3n + 2m*n + 11m^2
            println("Latest iterate X =")
            for i = j+1:j+n
                println("wa[$i] = ", wa[i])
            end
        else
            # "Compute function value f for the sample problem."
            f = 0.25 * (x[1] - 1)^2
            for i = 2:n
                f = f + (x[i] - x[i-1]^2)^2
            end

            f = 4 * f

            # "Compute gradient g for the sample problem."
            t₁ = x[2] - x[1]^2
            g[1] = 2 * (x[1] - 1) - 1.6e1 * x[1] * t₁
            for i = 2:n-1
                t₂ = t₁
                t₁ = x[i+1] - x[i]^2
                g[i] = 8 * t₂ - 1.6e1 * x[i] * t₁
            end
            g[n] = 8 * t₁
        end
            # "go back to the minimization routine."
            @goto CALLLBFGSB
    end

    # "the minimization routine has returned with a new iterate, The time limit
    #  has not been reached, and we test whether the following two stopping tests are satisfied:
    if task[1:5] == b"NEW_X"
        # "1) Terminate if the total number of f and g evaluations exceeds 900."
        isave[34] ≥ 900 && (task[1:4] = b"STOP"; terminate_info="TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT")

        # "2) Terminate if  |proj g|/(1+|f|) < 1.0d-10"
        dsave[13] ≤ 1e-10 * (1e0 + abs(f)) && (task[1:4] = b"STOP"; terminate_info="THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL")

        # "We wish to print the following information at each iteration:"
        # "1) the current iteration number, isave(30),"
        # "2) the total number of f and g evaluations, isave(34),"
        # "3) the value of the objective function f,"
        # "4) the norm of the projected gradient,  dsave(13)"

        println("Iterate ", isave[30], "  nfg = ", isave[34], "  f = ", f, "  |proj g| = ", dsave[13])

        # "If the run is to be terminated, we print also the information contained
        #  in task as well as the final value of x."
        if task[1:4] == b"STOP"
            println(terminate_info)
            for i = 1:n
                println("x($i) = ", x[i])
            end
        end

        # "go back to the minimization routine."
        @goto CALLLBFGSB
    end

    # ---------- the end of the loop -------------
    # "If task is neither FG nor NEW_X we terminate execution."

end

end # EOB


text:
import Base: view

using Base: ScalarIndex, ViewIndex, Slice, @_inline_meta, @boundscheck,
            to_indices, compute_offset1, unsafe_length, _maybe_reshape_parent, index_ndims

## construction

struct Contiguous end
struct NonContiguous end

# Detect whether the view is contiguous or not
CuIndexStyle() = Contiguous()
CuIndexStyle(I...) = NonContiguous()
CuIndexStyle(i1::Colon, ::ScalarIndex...) = Contiguous()
CuIndexStyle(i1::AbstractUnitRange, ::ScalarIndex...) = Contiguous()
CuIndexStyle(i1::Colon, I...) = CuIndexStyle(I...)

cuviewlength() = ()
cuviewlength(::Real, I...) = (@_inline_meta; cuviewlength(I...)) # skip scalars
cuviewlength(i1::AbstractUnitRange, I...) = (@_inline_meta; (unsafe_length(i1), cuviewlength(I...)...))
cuviewlength(i1::AbstractUnitRange, ::ScalarIndex...) = (@_inline_meta; (unsafe_length(i1),))

view(A::CuArray, I::Vararg{Any,N}) where {N} = (@_inline_meta; _cuview(A, I, CuIndexStyle(I...)))

function _cuview(A, I, ::Contiguous)
    @_inline_meta
    J = to_indices(A, I)
    @boundscheck checkbounds(A, J...)
    _cuview(_maybe_reshape_parent(A, index_ndims(J...)), J, cuviewlength(J...))
end

# for contiguous views just return a new CuArray
_cuview(A::CuArray{T}, I::NTuple{N,ViewIndex}, dims::NTuple{M,Integer}) where {T,N,M} =
    CuArray{T,M}(A.buf, dims; offset=A.offset + compute_offset1(A, 1, I) * sizeof(T), own=A.own)

# fallback to SubArray when the view is not contiguous
_cuview(A, I, ::NonContiguous) where {N} = invoke(view, Tuple{AbstractArray, typeof(I).parameters...}, A, I...)

## operations

# copyto! doesn't know how to deal with SubArrays, but broadcast does
# FIXME: use the rules from Adapt.jl to define copyto! methods in GPUArrays.jl
function Base.copyto!(dest::GPUArray{T,N}, src::SubArray{T,N,<:GPUArray{T,N}}) where {T,N}
    dest .= src
    dest
end

# copying to a CPU array requires an intermediate copy
# TODO: support other copyto! invocations (GPUArrays.jl copyto! defs + Adapt.jl rules)
function Base.copyto!(dest::AbstractArray{T,N}, src::SubArray{T,N,AT}) where {T,N,AT<:GPUArray{T}}
    temp = similar(AT, axes(src))
    copyto!(temp, src)
    copyto!(dest, temp)
end

# upload the SubArray indices when adapting to the GPU
# (can't do this eagerly or the view constructor wouldn't be able to boundscheck)
# FIXME: alternatively, have users do `cu(view(cu(A), inds))`, but that seems redundant
Adapt.adapt_structure(to::CUDAnative.Adaptor, A::SubArray) =
    SubArray(adapt(to, parent(A)), adapt(to, adapt(CuArray, parentindices(A))))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct SecurityProfile <: SwaggerModel
    encryptionAtHost::Any # spec type: Union{ Nothing, Bool } # spec name: encryptionAtHost

    function SecurityProfile(;encryptionAtHost=nothing)
        o = new()
        validate_property(SecurityProfile, Symbol("encryptionAtHost"), encryptionAtHost)
        setfield!(o, Symbol("encryptionAtHost"), encryptionAtHost)
        o
    end
end # type SecurityProfile

const _property_map_SecurityProfile = Dict{Symbol,Symbol}(Symbol("encryptionAtHost")=>Symbol("encryptionAtHost"))
const _property_types_SecurityProfile = Dict{Symbol,String}(Symbol("encryptionAtHost")=>"Bool")
Base.propertynames(::Type{ SecurityProfile }) = collect(keys(_property_map_SecurityProfile))
Swagger.property_type(::Type{ SecurityProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_SecurityProfile[name]))}
Swagger.field_name(::Type{ SecurityProfile }, property_name::Symbol) =  _property_map_SecurityProfile[property_name]

function check_required(o::SecurityProfile)
    true
end

function validate_property(::Type{ SecurityProfile }, name::Symbol, val)
end


text:
# Gradient boosting.
module GB

using Random
using GradientBoost.Util
using GradientBoost.LossFunctions

export GBAlgorithm,
       GBModel,
       stochastic_gradient_boost,
       fit,
       predict,
       build_base_func,
       create_sample_indices

# Gradient boost algorithm.
abstract type GBAlgorithm end

# Gradient boost model.
mutable struct GBModel{F <: AbstractFloat}
  learning_rate::F
  base_funcs::Vector{Function}
end

# Perform stochastic gradient boost.
#
# @param gb Gradient boosting algorithm.
# @param instances Instances.
# @param labels Labels.
# @return Gradient boost model.
function stochastic_gradient_boost(gb::GBAlgorithm, instances, labels)
  # Initialize base functions collection
  num_iterations = gb.num_iterations
  base_funcs = Array{Function, 1}(undef, num_iterations+1)

  # Create initial base function
  initial_val = minimizing_scalar(gb.loss_function, labels)
  initial_base_func = (instances) -> fill(initial_val, size(instances, 1))

  # Add initial base function to ensemble
  base_funcs[1] = initial_base_func

  # Build consecutive base functions
  prev_func_pred = initial_base_func(instances)
  for iter_ind = 2:num_iterations+1
    # Obtain current residuals
    psuedo = negative_gradient(
      gb.loss_function,
      labels,
      prev_func_pred
    )

    # Sample instances
    stage_sample_ind = create_sample_indices(gb, instances, labels)

    # Build current base function
    stage_base_func = build_base_func(
      gb,
      instances[stage_sample_ind, :],
      labels[stage_sample_ind],
      prev_func_pred[stage_sample_ind],
      psuedo[stage_sample_ind]
    )

    # Update previous function prediction
    prev_func_pred .+= gb.learning_rate .* stage_base_func(instances)

    # Add optimal base function to ensemble
    base_funcs[iter_ind] = stage_base_func
  end

  # Return model
  return GBModel(gb.learning_rate, base_funcs)
end

function fit(gb::GBAlgorithm, instances, labels)
  stochastic_gradient_boost(gb, instances, labels)
end
function predict(gb_model::GBModel, instances)
  outputs = gb_model.base_funcs[1](instances)
  for i = 2:length(gb_model.base_funcs)
    outputs .+= gb_model.learning_rate .* gb_model.base_funcs[i](instances)
  end
  return outputs
end

# Build base (basis) function for gradient boosting algorithm.
#
# @param gb Gradient boosting algorithm.
# @param instances Instances.
# @param labels Labels.
# @param prev_func_pred Previous base function's predictions.
# @param psuedo Psuedo-labels (psuedo-response).
# @return Function of form (instances) -> predictions.
function build_base_func(
  gb::GBAlgorithm,
  instances,
  labels,
  prev_func_pred,
  psuedo)

  err_must_be_overriden()
end

# Default sample method for gradient boosting algorithms.
# By default, it is sampling without replacement.
#
# @param gb Gradient boosting algorithm.
# @param instances Instances.
# @param labels Labels.
# @return Sample indices.
function create_sample_indices(gb::GBAlgorithm, instances, labels)
  n = size(instances, 1)
  prop = gb.sampling_rate

  ind = randperm(n)[1:Int(round(prop * n))]
end

end # module


text:
module Runtests

import ModiaMath

# Desired:
#   using Test
#
# In order that Test need not to be defined in the user environment, it is included via ModiaMath:
using ModiaMath.Test

@testset "Test ModiaMath" begin
    include(joinpath("result", "_includes.jl"))
    include(joinpath("variables", "_includes.jl"))
    include(joinpath("frames", "_includes.jl"))
    include(joinpath("nonlinearEquations", "_includes.jl"))
    include(joinpath("simulation", "_includes.jl"))

    println("\n... close all open figures.")
    ModiaMath.closeAllFigures()
end

end


text:
@testset "Transition flux report" begin

fName = "report/transitionflux"

excelFluxReport( mpSim, 12, ("", "A junior"), ("", "B junior"),
    ("", "Reserve junior"), filename = fName )
excelFluxReport( mpSim, 12, ("A junior", "A senior"), ("B junior", "B senior"),
    ("Reserve junior", "A senior"), ("Reserve junior", "B senior"),
    filename = fName, overwrite = false )
excelFluxReport( mpSim, 12, ("A junior", ""), ("B junior", ""),
    ("Reserve junior", ""), filename = fName, overwrite = false )
excelFluxReport( mpSim, 12, ("A senior", "Master"), ("B senior", "Master"),
    filename = fName, overwrite = false )
excelFluxReport( mpSim, 12, ("A senior", ""), ("B senior", ""), ("Master", ""),
    filename = fName, overwrite = false )

end  # @testset "Transition flux report"


text:
using Cairo
using Colors
using ColorSchemes
using CSV
using DataFrames
using Gadfly

include("./vb.jl")

function get_vb_df(base, comp)
    n = size(base)[1]

    ddict = Dict(
        "Real" => base,
        "Computed" => comp
    )

    source = String[]
    x = Int[]
    y = Int[]
    rate = Float64[]

    for s in ["Real", "Computed"]
        for i in 1:n
            for j in 1:n
                push!(source, s)
                push!(x, i)
                push!(y, j)
                push!(rate, ddict[s][n*(i-1)+j])
            end
        end
    end

    return DataFrame(
        :Source => source,
        :X => x,
        :Y => y,
        :Rate => rate
    )
end

function test_rates(csv_name, base_rates, write = false)
    base_rates = transpose(reduce(hcat, base_rates))

    n = size(base_rates)[1]

    base_counts = [[rand(Poisson(base_rates[i, j])) for j in 1:n] for i in 1:n]
    base_counts = transpose(reduce(hcat, base_counts))

    println("Initial rates: ")
    show(stdout, "text/plain", array_to_printable_mat(base_rates))
    println()

    est_rates, _, est_mean, est_std = variational_bayes(base_counts)
    println("Computed rates: ")
    show(stdout, "text/plain", array_to_printable_mat(est_rates))
    println()

    err = sqrt(mean([((base_rates[i]-est_rates[i])/base_rates[i])^2 for i in 1:length(est_rates)]))
    # err = sqrt(mean(([x^2 for x in vec(base_rates) - est_rates]))
    println("Log normalized RMS error: "*string(err))
    println("Estimated gaussian mean: "*string(est_mean))
    println("Estimated gaussian stdev: "*string(est_std))

    print("Exporting dataframe...")
    if write
        CSV.write(csv_name, get_vb_df(base_rates, est_rates))
    end
    println("Done!")

    return est_rates
end

function test_vb_5x5()
    base_rates = [
        [1, 1, 1, 1, 1],
        [1, 10, 50, 10, 1],
        [1, 50, 100, 50, 1],
        [1, 10, 50, 10, 1],
        [1, 1, 1, 1, 1],
    ]
    test_rates("../data/test_vb_5x5.csv", base_rates, true)
end

function test_vb_2hills()
    base_rates = [
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [0, 50, 50, 50, 0, 0, 0, 0, 0, 0],
        [0, 50, 100, 50, 0, 0, 0, 0, 0, 0],
        [0, 50, 50, 50, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 50, 50, 50, 0],
        [0, 0, 0, 0, 0, 0, 50, 100, 50, 0],
        [0, 0, 0, 0, 0, 0, 50, 50, 50, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    ]
    base_rates = [[x+1 for x in arr] for arr in base_rates]
    test_rates("../data/test_vb_2hills.csv", base_rates, true)
end

function test_vb_pacman()
    base_rates = [
        [0,  0,  0,  0,  0,  0,  0,  0, 0],
        [0, 10, 10, 10, 10, 10, 10, 10, 0],
        [0, 10, 90, 90, 90, 90, 90, 10, 0],
        [0, 10, 90, 10, 10, 10, 10, 10, 0],
        [0, 10, 90, 10,  0,  0,  0,  0, 0],
        [0, 10, 90, 10, 10, 10, 10, 10, 0],
        [0, 10, 90, 90, 90, 90, 90, 10, 0],
        [0, 10, 10, 10, 10, 10, 10, 10, 0],
        [0,  0,  0,  0,  0,  0,  0,  0, 0],
    ]
    base_rates = [[x+1 for x in arr] for arr in base_rates]
    test_rates("../data/test_vb_pacman.csv", base_rates, true)
end

function repeat_test_vb_pacman(N = 10)
    base_rates = [
        [0,  0,  0,  0,  0,  0,  0,  0, 0],
        [0, 10, 10, 10, 10, 10, 10, 10, 0],
        [0, 10, 90, 90, 90, 90, 90, 10, 0],
        [0, 10, 90, 10, 10, 10, 10, 10, 0],
        [0, 10, 90, 10,  0,  0,  0,  0, 0],
        [0, 10, 90, 10, 10, 10, 10, 10, 0],
        [0, 10, 90, 90, 90, 90, 90, 10, 0],
        [0, 10, 10, 10, 10, 10, 10, 10, 0],
        [0,  0,  0,  0,  0,  0,  0,  0, 0],
    ]
    base_rates = [[x+1 for x in arr] for arr in base_rates]
    rate_matrices = []
    for i in 1:N
        csvname = "../data/repeat_test_vb_pacman_"*string(i)*".csv"
        push!(rate_matrices, array_to_printable_mat(test_rates(csvname, base_rates)))
    end
    rate_matrix = hcat(sum((rate_matrices)/N)...)

    # println(transpose(hcat(base_rates...)))
    # println(rate_matrix)
    repeat_df = get_vb_df(transpose(hcat(base_rates...)), rate_matrix)
    CSV.write("../data/repeat_test_vb_pacman_mean.csv", repeat_df)
    return rate_matrix
end

function repeat_test_gaussian(μ = 4, σ=1, n = 10, num_iters = 3)
    base_rates = rand(LogNormal(μ, σ), n, n)
    base_rates = [base_rates[i, :] for i in 1:n]
    rate_matrices = []

    for i in 1:num_iters
        csvname = "../data/repeat_test_lgrf_"*string(i)*".csv"
        push!(rate_matrices, array_to_printable_mat(test_rates(csvname, base_rates)))
    end
    rate_matrix = hcat(sum((rate_matrices)/num_iters)...)
    repeat_df = get_vb_df(transpose(hcat(base_rates...)), rate_matrix)
    CSV.write("../data/repeat_test_lgrf_mean.csv", repeat_df)
    return rate_matrix
end


text:
using vOptGeneric
m = vModel()

@variable(m, x[1:6,1:6], Bin)
@variable(m, 0 <= y <= 10, Int)

M = 50
p1 =   [M 4 5 M M M;
        M M 2 1 2 7;
        M M M 5 2 M;
        M M 5 M M 3;
        M M M M M 4;
        M M M M M M]
p2 =   [M 3 1 M M M;
        M M 1 4 2 2;
        M M M 1 7 M;
        M M 1 M M 2;
        M M M M M 2;
        M M M M M M]

@addobjective(m, Min , y + sum(x[i,j]*p1[i,j] for i=1:6, j=1:6))
@addobjective(m, Min , y + sum(x[i,j]*p2[i,j] for i=1:6, j=1:6))

@constraint(m, sum(x[1,j] for j = 1:6) == 1)
@constraint(m, sum(x[i,6] for i = 1:6) == 1)
@constraint(m, cstr[i=2:5], sum(x[i,j] for j = 1:6) - sum(x[j,i] for j =1:6) == 0)
@constraint(m, y <= 5)

writeMOP(m, "bicriteria_shortest_path.mop")


text:
# Constans for the game
const BLOCK_SIZE = 20
const WIDTH = 640
const HEIGHT = 480

@enum Direction begin
    UP = 1
    RIGHT = 2
    DOWN = 3
    LEFT = 4
end

const MODELS_PATH = "saved_models"


text:
using Makie
using LinearAlgebra

res = 500
game = Scene(resolution = (res, res), raw = true, camera = campixel!)
cell = 10
cells = res ÷ 10
snakelen = 5
middle = cells ÷ 2
snakestart = middle - (snakelen ÷ 2)
segments = Point2f0.(middle, range(snakestart, step = cell, length = snakelen))
snake = scatter!(game, segments, markersize = cell, color = :black, marker = :rect)[end]
global last_dir = (0, 1)
dir = lift(game.events.keyboardbuttons) do but
    global last_dir
    ispressed(but, Keyboard.left) && return (last_dir = (-1, 0))
    ispressed(but, Keyboard.up) && return (last_dir = (0, 1))
    ispressed(but, Keyboard.right) && return (last_dir = (1, 0))
    ispressed(but, Keyboard.down) && return (last_dir = (0, -1))
    last_dir
end
newsnake = copy(segments)

food = scatter!(game, [Point2f0(res .÷ 2)], markersize = 2cell, marker = '☕')[end]

function run_game(game, snake, newsnake, food)
    display(game)
    spawntime = time()
    newspawn = 6
    speed = 1/10
    while isopen(game)
        curr_snake = snake[1][]
        circshift!(newsnake, curr_snake, 1)
        newsnake[1] = newsnake[2] .+ Point2f0(dir[] .* cell)
        spawndiff = time() - spawntime
        if spawndiff > 0
            if (newspawn - spawndiff) < 0.1 && (newspawn - spawndiff) > 0
                food[1][][1] = Point2f0(RNG.rand(1:cells), RNG.rand(1:cells)) .* cell

            elseif (newspawn - spawndiff) < 0.1
                food[1][][1] = Point2f0(res + cell)
                food[1][] = food[1][] # update array
                spawntime = time()
            end
        end
        if norm(newsnake[1] .- food[1][]) < 3cell
            snake[:color] = :red
            food[1][][1] = Point2f0(res + cell)
            # append segments to snake
            a, b = newsnake[end], newsnake[end - 1]
            push!(newsnake, a .+ (b .- a))
            push!(curr_snake, Point2f0(0))
            speed = speed - (speed * 0.5)
            spawntime = time()
        else
            snake[:color] = :black
        end
        food[1][] = food[1][] # update array
        snake[1] = newsnake
        newsnake = curr_snake
        sleep(speed);
        force_update!()
    end
end
run_game(game, snake, newsnake, food)


text:
baremodule JuliaProjectFormatter

export format_projects

function format_projects end

module Internal

using ..JuliaProjectFormatter: JuliaProjectFormatter

using TOML

const _project_key_order = ["name", "uuid", "keywords", "license", "desc", "deps", "compat"]
project_key_order(key::String) =
    something(findfirst(x -> x == key, _project_key_order), length(_project_key_order) + 1)

print_project(io, dict) =
    TOML.print(io, dict, sorted = true, by = key -> (project_key_order(key), key))

discover_project_toml(path) = discover_project_toml!(String[], path)

function discover_project_toml!(projects, path)
    isfile(path) && return push!(projects, path)
    for name in readdir(path)
        name == ".git" && continue
        subpath = joinpath(path, name)
        if isdir(subpath)
            discover_project_toml!(projects, subpath)
        elseif name in ("Project.toml", "JuliaProject.toml")
            push!(projects, subpath)
        end
    end
    return projects
end

function JuliaProjectFormatter.format_projects(path = pwd())
    for project in discover_project_toml(path)
        @debug "Formatting `$project`"
        dict = TOML.parsefile(project)
        open(project; write = true) do io
            print_project(io, dict)
        end
    end
end

end  # module Internal

end  # baremodule JuliaProjectFormatter


text:
@test rand() != rand()
@test 0.0 <= rand() < 1.0
@test rand(Uint32) >= 0
@test -10 <= rand(-10:-5) <= -5
@test -10 <= rand(-10:5) <= 5
@test minimum([rand(int32(1):int32(7^7)) for i = 1:100000]) > 0
@test(typeof(rand(false:true)) == Bool)

@test length(randn(4, 5)) == 20
@test length(randbool(4, 5)) == 20

@test rand(MersenneTwister()) == 0.8236475079774124
@test rand(MersenneTwister(0)) == 0.8236475079774124
@test rand(MersenneTwister(42)) == 0.5331830160438613
# Try a seed larger than 2^32
@test rand(MersenneTwister(5294967296)) == 0.3498809918210497

# randn
@test randn(MersenneTwister(42)) == -0.5560268761438381
A = zeros(2, 2)
randn!(MersenneTwister(42), A)
@test A == [-0.5560268761438381  0.02715533800914659;
            -0.4443833571072952 -0.29948409035852047]

for T in (Int8, Uint8, Int16, Uint16, Int32, Uint32, Int64, Uint64, Int128, Uint128,
          Char, Float16, Float32, Float64, Rational{Int})
    r = rand(convert(T, 97):convert(T, 122))
    @test typeof(r) == T
    @test 97 <= r <= 122
    r = rand(convert(T, 97):convert(T,2):convert(T, 122),2)[1]
    @test typeof(r) == T
    @test 97 <= r <= 122
    @test mod(r,2)==1
end

if sizeof(Int32) < sizeof(Int)
    r = rand(int32(-1):typemax(Int32))
    @test typeof(r) == Int32
    @test -1 <= r <= typemax(Int32)
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
export AbstractPolyhedraOptimizer

struct PolyhedraOptSet{T, RepT <: Rep{T}} <: MOI.AbstractVectorSet
    rep::RepT
end
MOI.dimension(set::PolyhedraOptSet) = fulldim(set.rep)

function JuMP.build_constraint(error_fun::Function, func, set::Rep)
    return JuMP.BridgeableConstraint(
        JuMP.build_constraint(error_fun, func, PolyhedraOptSet(set)),
        PolyhedraToLPBridge)
end

abstract type AbstractPolyhedraOptimizer{T} <: MOI.AbstractOptimizer end

function MOI.copy_to(dest::AbstractPolyhedraOptimizer, src::MOI.ModelLike; kws...)
    return MOI.Utilities.automatic_copy_to(dest, src; kws...)
end
MOI.Utilities.supports_default_copy_to(optimizer::AbstractPolyhedraOptimizer, copy_names::Bool) = true

function MOI.add_variable(optimizer::AbstractPolyhedraOptimizer)
    return MOI.add_variable(optimizer.lphrep.model)
end
function MOI.add_variables(optimizer::AbstractPolyhedraOptimizer, n)
    return MOI.add_variables(optimizer.lphrep.model, n)
end

function MOI.supports(::AbstractPolyhedraOptimizer{T},
                      ::Union{MOI.ObjectiveSense,
                              MOI.ObjectiveFunction{MOI.SingleVariable},
                              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}}) where T
    return true
end
function MOI.set(optimizer::AbstractPolyhedraOptimizer, ::MOI.ObjectiveSense,
                 sense::MOI.OptimizationSense)
    optimizer.objective_sense = sense
end
function MOI.set(optimizer::AbstractPolyhedraOptimizer{T}, ::MOI.ObjectiveFunction,
                 func::MOI.SingleVariable) where T
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
            convert(MOI.ScalarAffineFunction{T}, func))
end
function MOI.set(optimizer::AbstractPolyhedraOptimizer, ::MOI.ObjectiveFunction,
                 func::MOI.ScalarAffineFunction)
    indices = [term.variable_index.value for term in func.terms]
    coefs = [term.coefficient for term in func.terms]
    optimizer.objective_func = sparsevec(indices, coefs, fulldim(optimizer.lphrep))
    optimizer.objective_constant = func.constant
end

function MOI.supports_constraint(::AbstractPolyhedraOptimizer{T},
                                 ::Type{MOI.ScalarAffineFunction{T}},
                                 ::Type{<:Union{MOI.EqualTo{T}, MOI.LessThan{T}}}) where T
    return true
end
function MOI.get(optimizer::AbstractPolyhedraOptimizer{T},
                 attr::Union{MOI.ConstraintFunction,
                             MOI.ConstraintSet},
                 ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},
                                         <:Union{MOI.EqualTo{T},
                                                 MOI.LessThan{T}}}) where T
    return MOI.get(optimizer.lphrep.model, attr, ci)
end
function MOI.add_constraint(optimizer::AbstractPolyhedraOptimizer{T},
                            func::MOI.ScalarAffineFunction{T},
                            set::Union{MOI.EqualTo{T}, MOI.LessThan{T}}) where T
    optimizer.feasible_set = nothing # Invalidated by the new constraint
    return MOI.add_constraint(optimizer.lphrep, func, set)
end
function MOI.supports_constraint(::AbstractPolyhedraOptimizer,
                                 ::Type{MOI.VectorOfVariables},
                                 ::Type{<:PolyhedraOptSet})
    return true
end
function MOI.add_constraint(optimizer::AbstractPolyhedraOptimizer,
                            vov::MOI.VectorOfVariables, rep::PolyhedraOptSet)
    if vov.variables != MOI.get(optimizer.lphrep.model, MOI.ListOfVariableIndices())
        error("Cannot only add VectorOfVariables polyhedra constraints with all variables in creation order.")
    end
    if optimizer.rep !== nothing
        error("Cannot only add one polyhedra constraint.")
    end
    optimizer.feasible_set = nothing # Invalidated by the new constraint
    optimizer.rep = rep.rep
    return MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(rep)}(1)
end

coefficient_type(::MOI.ModelLike) = Float64

# Inspired by `JuMP.set_optimizer`
function layered_optimizer(factory::JuMP.OptimizerFactory)
    optimizer = factory()
    T = coefficient_type(optimizer)
    if !MOIU.supports_default_copy_to(optimizer, false)
        universal_fallback = MOIU.UniversalFallback(_MOIModel{T}())
        optimizer = MOIU.CachingOptimizer(universal_fallback, optimizer)
    end
    optimizer = MOI.Bridges.full_bridge_optimizer(optimizer, T)
    MOI.Bridges.add_bridge(optimizer, PolyhedraToLPBridge{T})
    return optimizer, T
end

"""
    isempty(p::Rep, solver::JuMP.OptimizerFactory=Polyhedra.linear_objective_solver(p))

Check whether the polyhedron `p` is empty by using the solver `solver`.
"""
function Base.isempty(p::Rep, solver::Solver=Polyhedra.linear_objective_solver(p))
    N = fulldim(p)
    if N == -1
        if p isa VRepresentation || (p isa Polyhedron && fulldim(vrep(p)) == -1)
            # Empty V-representation means empty polyhedron
            return true
        else
            # Empty H-representation means the polyhedron the full space
            return false
        end
    end
    model, T = layered_optimizer(solver)
    x, cx = MOI.add_constrained_variables(model, PolyhedraOptSet(p))
    MOI.optimize!(model)
    term = MOI.get(model, MOI.TerminationStatus())
    if term == MOI.OPTIMAL
        return false
    elseif term == MOI.INFEASIBLE || term == MOI.INFEASIBLE_OR_UNBOUNDED
        # The problem has no objective so it cannot be unbounded so
        # the `MOI.INFEASIBLE_OR_UNBOUNDED` status is accepted
        return true
    else
        error("Cannot determine whether the polyhedron is empty or not because",
              " the linear program terminated with status $term.")
    end
end


text:
RuntimeGeneratedFunctions.init(@__MODULE__)
"""
Algorithm for solving Physics-Informed Neural Networks problems.

Arguments:
* `chain`: a Flux.jl chain with a d-dimensional input and a 1-dimensional output,
* `strategy`: determines which training strategy will be used,
* `init_params`: the initial parameter of the neural network,
* `phi`: a trial solution,
* `derivative`: method that calculates the derivative.

"""

struct PhysicsInformedNN{C,P,PH,DER,T,K}
  chain::C
  init_params::P
  phi::PH
  derivative::DER
  strategy::T
  kwargs::K
end

function PhysicsInformedNN(chain,
                           strategy;
                           _init_params = nothing,
                           _phi = nothing,
                           _derivative = nothing,
                           kwargs...)
    if _init_params == nothing
        if chain isa AbstractArray
            initθ = DiffEqFlux.initial_params.(chain)
        else
            initθ = DiffEqFlux.initial_params(chain)
        end

    else
        initθ = _init_params
    end

    if _phi == nothing
        if chain isa AbstractArray
            phi = get_phi.(chain)
        else
            phi = get_phi(chain)
        end
    else
        phi = _phi
    end

    if _derivative == nothing
        derivative = get_numeric_derivative()
    else
        derivative = _derivative
    end

    PhysicsInformedNN(chain,initθ,phi,derivative, strategy, kwargs)
end

abstract type TrainingStrategies  end

"""
* `dx`: the discretization of the grid.
"""
struct GridTraining <: TrainingStrategies
    dx
end

"""
* `points`: number of points in random select training set.
"""
struct StochasticTraining <:TrainingStrategies
    points:: Int64
end

"""
* `points`:  the number of quasi-random points in minibatch,
* `sampling_alg`: the quasi-Monte Carlo sampling algorithm,
* `minibatch`: the number of subsets.

For more information look: QuasiMonteCarlo.jl https://github.com/SciML/QuasiMonteCarlo.jl
"""
struct QuasiRandomTraining <:TrainingStrategies
    points:: Int64
    sampling_alg::QuasiMonteCarlo.SamplingAlgorithm
    minibatch:: Int64
end
function QuasiRandomTraining(points;sampling_alg = UniformSample(),minibatch=500)
    QuasiRandomTraining(points,sampling_alg,minibatch)
end
"""
* `quadrature_alg`: quadrature algorithm,
* `reltol`: relative tolerance,
* `abstol`: absolute tolerance,
* `maxiters`: the maximum number of iterations in quadrature algorithm,
* `batch`: the preferred number of points to batch.

For more information look: Quadrature.jl https://github.com/SciML/Quadrature.jl
"""
struct QuadratureTraining <: TrainingStrategies
    quadrature_alg::DiffEqBase.AbstractQuadratureAlgorithm
    reltol::Float64
    abstol::Float64
    maxiters::Int64
    batch::Int64
end

function QuadratureTraining(;quadrature_alg=HCubatureJL(),reltol= 1e-6,abstol= 1e-3,maxiters=1e3,batch=0)
    QuadratureTraining(quadrature_alg,reltol,abstol,maxiters,batch)
end

"""
Create dictionary: variable => unique number for variable

# Example 1

Dict{Symbol,Int64} with 3 entries:
  :y => 2
  :t => 3
  :x => 1

# Example 2

 Dict{Symbol,Int64} with 2 entries:
  :u1 => 1
  :u2 => 2
"""
get_dict_vars(vars) = Dict( [Symbol(v) .=> i for (i,v) in enumerate(vars)])

# Wrapper for _transform_derivative
function transform_derivative(ex,dict_indvars,dict_depvars)
    if ex isa Expr
        ex.args = _transform_derivative(ex.args,dict_indvars,dict_depvars)
    end
    return ex
end

function get_ε(dim, der_num)
    epsilon = cbrt(eps(Float32))
    ε = zeros(Float32, dim)
    ε[der_num] = epsilon
    ε
end

θ = gensym("θ")

"""
Transform the derivative expression to inner representation

# Examples

1. First compute the derivative of function 'u(x,y)' with respect to x.

Take expressions in the form: `derivative(u(x,y), x)` to `derivative(phi, u, [x, y], εs, order, θ)`,
where
 phi - trial solution
 u - function
 x,y - coordinates of point
 εs - epsilon mask
 order - order of derivative
 θ - weight in neural network
"""
function _transform_derivative(_args,dict_indvars,dict_depvars)
    for (i,e) in enumerate(_args)
        if !(e isa Expr)
            if e in keys(dict_depvars)
                depvar = _args[1]
                num_depvar = dict_depvars[depvar]
                indvars = _args[2:end]
                _args = if length(dict_depvars) == 1
                    [:u, :([$(indvars...)]), :($θ), :phi]
                else
                    [:u, :([$(indvars...)]), Symbol(:($θ),num_depvar), Symbol(:phi,num_depvar)]
                end
                break
            elseif e == :derivative
                derivative_variables = Symbol[]
                order = 0
                while (_args[1] == :derivative)
                    order += 1
                    push!(derivative_variables, _args[end])
                    _args = _args[2].args
                end
                depvar = _args[1]
                num_depvar = dict_depvars[depvar]
                indvars = _args[2:end]
                dim_l = length(indvars)
                εs = [get_ε(dim_l,d) for d in 1:dim_l]
                undv = [dict_indvars[d_p] for d_p  in derivative_variables]
                εs_dnv = [εs[d] for d in undv]
                _args = if length(dict_depvars) == 1
                    [:derivative, :phi, :u, :([$(indvars...)]), εs_dnv, order, :($θ)]
                else
                    [:derivative, Symbol(:phi,num_depvar), :u, :([$(indvars...)]), εs_dnv, order, Symbol(:($θ),num_depvar)]
                end
                break
            end
        else
            _args[i].args = _transform_derivative(_args[i].args,dict_indvars,dict_depvars)
        end
    end
    return _args
end

"""
Parse ModelingToolkit equation form to the inner representation.

Example:

1)  1-D ODE: Dt(u(t)) ~ t +1

    Take expressions in the form: 'Equation(derivative(u(t), t), t + 1)' to 'derivative(phi, u_d, [t], [[ε]], 1, θ) - (t + 1)'

2)  2-D PDE: Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)

    Take expressions in the form:
     Equation(derivative(derivative(u(x, y), x), x) + derivative(derivative(u(x, y), y), y), -(sin(πx)) * sin(πy))
    to
     (derivative(phi,u, [x, y], [[ε,0],[ε,0]], 2, θ) + derivative(phi, u, [x, y], [[0,ε],[0,ε]], 2, θ)) - -(sin(πx)) * sin(πy)

3)  System of PDEs: [Dx(u1(x,y)) + 4*Dy(u2(x,y)) ~ 0,
                    Dx(u2(x,y)) + 9*Dy(u1(x,y)) ~ 0]

    Take expressions in the form:
    2-element Array{Equation,1}:
        Equation(derivative(u1(x, y), x) + 4 * derivative(u2(x, y), y), ModelingToolkit.Constant(0))
        Equation(derivative(u2(x, y), x) + 9 * derivative(u1(x, y), y), ModelingToolkit.Constant(0))
    to
      [(derivative(phi1, u1, [x, y], [[ε,0]], 1, θ1) + 4 * derivative(phi2, u, [x, y], [[0,ε]], 1, θ2)) - 0,
       (derivative(phi2, u2, [x, y], [[ε,0]], 1, θ2) + 9 * derivative(phi1, u, [x, y], [[0,ε]], 1, θ1)) - 0]
"""

function build_symbolic_equation(eq,_indvars,_depvars)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    parse_equation(eq,dict_indvars,dict_depvars)
end

function parse_equation(eq,dict_indvars,dict_depvars)
    eq_lhs = isequal(expand_derivatives(eq.lhs), 0) ? eq.lhs : expand_derivatives(eq.lhs)
    eq_rhs = isequal(expand_derivatives(eq.rhs), 0) ? eq.rhs : expand_derivatives(eq.rhs)

    left_expr = transform_derivative(toexpr(eq_lhs),
                                     dict_indvars,dict_depvars)
    right_expr = transform_derivative(toexpr(eq_rhs),
                                     dict_indvars,dict_depvars)
    loss_func = :($left_expr - $right_expr)
end

"""
Build a loss function for a PDE or a boundary condition

# Examples: System of PDEs:

Take expressions in the form:

[Dx(u1(x,y)) + 4*Dy(u2(x,y)) ~ 0,
 Dx(u2(x,y)) + 9*Dy(u1(x,y)) ~ 0]

to

:((cord, θ, phi, derivative, u)->begin
          #= ... =#
          #= ... =#
          begin
              (θ1, θ2) = (θ[1:33], θ"[34:66])
              (phi1, phi2) = (phi[1], phi[2])
              let (x, y) = (cord[1], cord[2])
                  [(+)(derivative(phi1, u, [x, y], [[ε, 0.0]], 1, θ1), (*)(4, derivative(phi2, u, [x, y], [[0.0, ε]], 1, θ2))) - 0,
                   (+)(derivative(phi2, u, [x, y], [[ε, 0.0]], 1, θ2), (*)(9, derivative(phi1, u, [x, y], [[0.0, ε]], 1, θ1))) - 0]
              end
          end
      end)
"""
function build_symbolic_loss_function(eqs,_indvars,_depvars, phi, derivative,initθ; bc_indvars=nothing)
    # dictionaries: variable -> unique number
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    bc_indvars = bc_indvars==nothing ? indvars : bc_indvars
    return build_symbolic_loss_function(eqs,indvars,depvars,
                                        dict_indvars,dict_depvars,
                                        phi, derivative,initθ,
                                        bc_indvars = bc_indvars)
end

function build_symbolic_loss_function(eqs,indvars,depvars,
                                      dict_indvars,dict_depvars,
                                      phi, derivative, initθ;
                                      bc_indvars = indvars)
    if !(eqs isa Array)
        eqs = [eqs]
    end
    loss_functions= Expr[]
    for eq in eqs
        push!(loss_functions,parse_equation(eq,dict_indvars,dict_depvars))
    end

    vars = :(cord, $θ, phi, derivative,u)
    ex = Expr(:block)
    if length(depvars) != 1
        θ_nums = Symbol[]
        phi_nums = Symbol[]
        for v in depvars
            num = dict_depvars[v]
            push!(θ_nums,:($(Symbol(:($θ),num))))
            push!(phi_nums,:($(Symbol(:phi,num))))
        end

        expr_θ = Expr[]
        expr_phi = Expr[]

        acum =  [0;accumulate(+, length.(initθ))]
        sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]

        for i in eachindex(depvars)
            push!(expr_θ, :($θ[$(sep[i])]))
            push!(expr_phi, :(phi[$i]))
        end

        vars_θ = Expr(:(=), build_expr(:tuple, θ_nums), build_expr(:tuple, expr_θ))
        push!(ex.args,  vars_θ)

        vars_phi = Expr(:(=), build_expr(:tuple, phi_nums), build_expr(:tuple, expr_phi))
        push!(ex.args,  vars_phi)
    end
    indvars_ex = [:($:cord[$i]) for (i, u) ∈ enumerate(bc_indvars)]

    left_arg_pairs, right_arg_pairs = bc_indvars,indvars_ex

    vars_eq = Expr(:(=), build_expr(:tuple, left_arg_pairs), build_expr(:tuple, right_arg_pairs))
    let_ex = Expr(:let, vars_eq, build_expr(:vect, loss_functions))
    push!(ex.args,  let_ex)

    expr_loss_function = :(($vars) -> begin $ex end)
end

function build_loss_function(eqs,_indvars,_depvars, phi, derivative,initθ;bc_indvars=nothing)
    # dictionaries: variable -> unique number
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    bc_indvars = bc_indvars==nothing ? indvars : bc_indvars
    return build_loss_function(eqs,indvars,depvars,
                               dict_indvars,dict_depvars,
                               phi, derivative,initθ,
                               bc_indvars = bc_indvars)
end

function build_loss_function(eqs,indvars,depvars,
                             dict_indvars,dict_depvars,
                             phi, derivative, initθ;
                             bc_indvars = indvars)

     expr_loss_function = build_symbolic_loss_function(eqs,indvars,depvars,
                                                       dict_indvars,dict_depvars,
                                                       phi, derivative, initθ;
                                                       bc_indvars = bc_indvars)
    u = get_u()
    _loss_function = @RuntimeGeneratedFunction(expr_loss_function)
    loss_function = (cord, θ) -> _loss_function(cord, θ, phi, derivative, u)
    return loss_function
end

function get_vars(indvars_, depvars_)
    depvars = [nameof(value(d)) for d in depvars_]
    indvars = [nameof(value(i)) for i in indvars_]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    return depvars,indvars,dict_indvars,dict_depvars
end

function get_bc_varibles(bcs,_indvars::Array,_depvars::Array)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    return get_bc_varibles(bcs,dict_indvars,dict_depvars)
end

function get_bc_varibles(bcs,dict_indvars,dict_depvars)
    bc_args = get_bc_argument(bcs,dict_indvars,dict_depvars)
    return map(barg -> filter(x -> x isa Symbol, barg), bc_args)
end

# Get arguments from boundary condition functions
function get_bc_argument(bcs,dict_indvars,dict_depvars)
    bc_args = []
    for _bc in bcs
        _bc isa Array && error("boundary conditions must be represented as a one-dimensional array")
        left_expr = transform_derivative(toexpr(_bc.lhs),
                                         dict_indvars,dict_depvars)
        bc_arg = if (left_expr.args[1] == :derivative)
            left_expr.args[4].args
        else
            left_expr.args[2].args
        end
        push!(bc_args, bc_arg)
    end

    return bc_args
end

function generate_training_sets(domains,dx,bcs,_indvars::Array,_depvars::Array)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    return generate_training_sets(domains,dx,bcs,dict_indvars,dict_depvars)
end
# Generate training set in the domain and on the boundary
function generate_training_sets(domains,dx,bcs,dict_indvars::Dict,dict_depvars::Dict)
    if dx isa Array
        dxs = dx
    else
        dxs = fill(dx,length(domains))
    end

    bound_args = get_bc_argument(bcs,dict_indvars,dict_depvars)
    bound_vars = get_bc_varibles(bcs,dict_indvars,dict_depvars)

    spans = [d.domain.lower:dx:d.domain.upper for (d,dx) in zip(domains,dxs)]
    dict_var_span = Dict([Symbol(d.variables) => d.domain.lower:dx:d.domain.upper for (d,dx) in zip(domains,dxs)])

    dif = [Float64[] for i=1:size(domains)[1]]
    for _args in bound_args
        for (i,x) in enumerate(_args)
            if x isa Number
                push!(dif[i],x)
            end
        end
    end
    collect_spans = collect.(spans)
    bc_data = map(zip(dif,collect_spans)) do (d,c)
        setdiff(c, d)
    end

    train_set = vec(map(points -> collect(points), Iterators.product(spans...)))
    pde_train_set = vec(map(points -> collect(points), Iterators.product(bc_data...)))

    bcs_train_set = map(bound_vars) do bt
        span = map(b -> dict_var_span[b], bt)
        _set = vec(map(points -> collect(points), Iterators.product(span...)))
    end

    [pde_train_set,bcs_train_set,train_set]
end

function get_bounds(domains,bcs,_indvars::Array,_depvars::Array)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    return get_bounds(domains,bcs,dict_indvars,dict_depvars)
end

function get_bounds(domains,bcs,dict_indvars,dict_depvars)
    bound_vars = get_bc_varibles(bcs,dict_indvars,dict_depvars)

    pde_lower_bounds = [d.domain.lower for d in domains]
    pde_upper_bounds = [d.domain.upper for d in domains]
    pde_bounds= [pde_lower_bounds,pde_upper_bounds]

    dict_lower_bound = Dict([Symbol(d.variables) => d.domain.lower for d in domains])
    dict_upper_bound = Dict([Symbol(d.variables) => d.domain.upper for d in domains])

    bcs_lower_bounds = map(bound_vars) do bt
        map(b -> dict_lower_bound[b], bt)
    end
    bcs_upper_bounds = map(bound_vars) do bt
        map(b -> dict_upper_bound[b], bt)
    end
    bcs_bounds= [bcs_lower_bounds,bcs_upper_bounds]

    [pde_bounds, bcs_bounds]
end

function get_phi(chain)
    # The phi trial solution
    if chain isa FastChain
        phi = (x,θ) -> chain(adapt(typeof(θ),x),θ)
    else
        _,re  = Flux.destructure(chain)
        phi = (x,θ) -> re(θ)(adapt(typeof(θ),x))
    end
    phi
end

function get_u()
	u = (cord, θ, phi)->sum(phi(cord, θ))
end
# the method to calculate the derivative
function get_numeric_derivative()
    epsilon = 2*cbrt(eps(Float32))
    derivative = (phi,u,x,εs,order,θ) ->
    begin
        ε = εs[order]
        ε = adapt(typeof(θ),ε)
        x = adapt(typeof(θ),x)
        if order > 1
            return (derivative(phi,u,x+ε,εs,order-1,θ)
                  - derivative(phi,u,x-ε,εs,order-1,θ))/epsilon
        else
            return (u(x+ε,θ,phi) - u(x-ε,θ,phi))/epsilon
        end
    end
    derivative
end

function get_loss_function(loss_functions, train_sets, strategy::GridTraining)
    # norm coefficient for loss function
    τ_ = loss_functions isa Array ? sum(length(train_set) for train_set in train_sets) : length(train_sets)
    τ = 1.0f0 / τ_

    function inner_loss(loss_function,x,θ)
        sum(abs2,loss_function(x, θ))
    end

    if !(loss_functions isa Array)
        loss_functions = [loss_functions]
        train_sets = [train_sets]
    end
    f = (loss,train_set,θ) -> sum([inner_loss(loss,x,θ) for x in train_set])
    loss = (θ) ->  τ * sum(f(loss_function,train_set,θ) for (loss_function,train_set) in zip(loss_functions,train_sets))
    return loss
end

function get_loss_function(loss_functions, bounds, strategy::StochasticTraining)
    points = strategy.points
    lbs,ubs = bounds

    function inner_loss(loss_function,x,θ)
        sum(loss_function(x, θ))
    end

    if !(loss_functions isa Array)
        loss_functions = [loss_functions]
        lbs = [lbs]
        ubs = [ubs]
    end
    τ = 1.0f0 / points

    loss = (θ) -> begin
        total = 0.
        for (lb, ub,l) in zip(lbs, ubs, loss_functions)
            len = length(lb)
            for i in 1:points
                r_point = lb .+ ub .* rand(len)
                total += inner_loss(l,r_point,θ)^2
            end
        end
        return τ * total
    end
    return loss
end

function get_loss_function(loss_functions, bounds, strategy::QuasiRandomTraining)
    sampling_alg = strategy.sampling_alg
    points = strategy.points
    minibatch = strategy.minibatch
    lbs,ubs = bounds

    function inner_loss(loss_function,x,θ)
        sum(loss_function(x, θ))
    end

    if !(loss_functions isa Array)
        loss_functions = [loss_functions]
        lbs = [lbs]
        ubs = [ubs]
    end
    τ = points
    ss =[]
    for (lb, ub) in zip(lbs, ubs)
        s = QuasiMonteCarlo.generate_design_matrices(points,lb,ub,sampling_alg,minibatch)
        push!(ss,s)
    end
    loss = (θ) -> begin
        total = 0.
        for (lb, ub,s_,l) in zip(lbs,ubs,ss,loss_functions)
            s =  s_[rand(1:minibatch)]
            step_ = size(lb)[1]
            k = size(lb)[1]-1
            for i in 1:step_:step_*points
                r_point = s[i:i+k]
                total += inner_loss(l,r_point,θ)^2
            end
        end
        return (1.0f0/τ) * total
    end
    return loss
end

function get_loss_function(loss_functions, bounds, strategy::QuadratureTraining)
    lbs,ubs = bounds

    function inner_loss(loss_function,x,θ)
        sum(loss_function(x, θ))
    end

    if !(loss_functions isa Array)
        loss_functions = [loss_functions]
        lbs = [lbs]
        ubs = [ubs]
    end

    τ = 1.0f0 / ((10)^length(ubs[1])*length(ubs))

    f = (lb,ub,loss_,θ) -> begin
        _loss = (x,θ) -> sum(abs2,inner_loss(loss_, x, θ))
        prob = QuadratureProblem(_loss,lb,ub,θ;batch = strategy.batch)
        solve(prob,
              strategy.quadrature_alg,
              reltol = strategy.reltol,
              abstol = strategy.abstol,
              maxiters = strategy.maxiters)[1]
    end
    loss = (θ) -> τ*sum(f(lb,ub,loss_,θ) for (lb,ub,loss_) in zip(lbs,ubs,loss_functions))
    return loss
end
function symbolic_discretize(pde_system::PDESystem, discretization::PhysicsInformedNN)
    eqs = pde_system.eq
    bcs = pde_system.bcs

    domains = pde_system.domain
    # dimensionality of equation
    dim = length(domains)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(pde_system.indvars, pde_system.depvars)

    chain = discretization.chain
    initθ = discretization.init_params
    phi = discretization.phi
    derivative = discretization.derivative
    strategy = discretization.strategy

    symbolic_pde_loss_function = build_symbolic_loss_function(eqs,indvars,depvars,
                                                              dict_indvars,dict_depvars,
                                                              phi, derivative,initθ)

    bc_indvars = get_bc_varibles(bcs,dict_indvars,dict_depvars)
    symbolic_bc_loss_functions = [build_symbolic_loss_function(bc,indvars,depvars,
                                                               dict_indvars,dict_depvars,
                                                               phi, derivative,initθ;
                                                               bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]
    symbolic_pde_loss_function,symbolic_bc_loss_functions
end
# Convert a PDE problem into an OptimizationProblem
function DiffEqBase.discretize(pde_system::PDESystem, discretization::PhysicsInformedNN)
    eqs = pde_system.eq
    bcs = pde_system.bcs

    domains = pde_system.domain
    # dimensionality of equation
    dim = length(domains)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(pde_system.indvars,pde_system.depvars)

    chain = discretization.chain
    initθ = discretization.init_params
    flat_initθ = if length(depvars) != 1 vcat(initθ...) else  initθ end
    phi = discretization.phi
    derivative = discretization.derivative
    strategy = discretization.strategy

    _pde_loss_function = build_loss_function(eqs,indvars,depvars,
                                                 dict_indvars,dict_depvars,phi, derivative, initθ)

    bc_indvars = get_bc_varibles(bcs,dict_indvars,dict_depvars)
    _bc_loss_functions = [build_loss_function(bc,indvars,depvars,
                                                  dict_indvars,dict_depvars,
                                                  phi, derivative, initθ;
                                                  bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]

    pde_loss_function, bc_loss_function =
    if strategy isa GridTraining
        dx = strategy.dx

        train_sets = generate_training_sets(domains,dx,bcs,
                                            dict_indvars,dict_depvars)

        # the points in the domain and on the boundary
        pde_train_set,bcs_train_set,train_set = train_sets

        pde_loss_function =get_loss_function(_pde_loss_function,
                                                        pde_train_set,
                                                        strategy)

        bc_loss_function =get_loss_function(_bc_loss_functions,
                                                       bcs_train_set,
                                                       strategy)

        (pde_loss_function, bc_loss_function)
    elseif strategy isa StochasticTraining
          bounds = get_bounds(domains,bcs,dict_indvars,dict_depvars)
          pde_bounds, bcs_bounds = bounds
          pde_loss_function = get_loss_function(_pde_loss_function,
                                                          pde_bounds,
                                                          strategy)
          lbs,ubs = bcs_bounds
          pl = length(lbs)
          bl = length(lbs[1])
          points = length(lbs[1]) == 0 ? 1 : Int(round(strategy.points^(bl/pl)))
          strategy = StochasticTraining(points)

          bc_loss_function = get_loss_function(_bc_loss_functions,
                                                         bcs_bounds,
                                                         strategy)
          (pde_loss_function, bc_loss_function)
    elseif strategy isa QuasiRandomTraining
         bounds = get_bounds(domains,bcs,dict_indvars,dict_depvars)
         pde_bounds, bcs_bounds = bounds
         pde_loss_function = get_loss_function(_pde_loss_function,
                                                        pde_bounds,
                                                        strategy)
         plbs,pubs = pde_bounds
         blbs,bubs = bcs_bounds
         pl = length(plbs)
         bl = length(blbs[1])
         points = Int(round(strategy.points^(bl/pl)))
         strategy = QuasiRandomTraining(points;
                                        sampling_alg = strategy.sampling_alg,
                                        minibatch = strategy.minibatch)

         bc_loss_function = get_loss_function(_bc_loss_functions,
                                                       bcs_bounds,
                                                       strategy)
         (pde_loss_function, bc_loss_function)
    elseif strategy isa QuadratureTraining
        dim<=1 && error("QuadratureTraining works only with dimensionality more than 1")

        (dim<=2 && strategy.quadrature_alg in [CubaCuhre(),CubaDivonne()]
        && error("$(strategy.quadrature_alg) works only with dimensionality more than 2"))

        bounds = get_bounds(domains,bcs,dict_indvars,dict_depvars)
        pde_bounds, bcs_bounds = bounds

        pde_loss_function = get_loss_function(_pde_loss_function,
                                              pde_bounds,
                                              strategy)

        bc_loss_function = get_loss_function(_bc_loss_functions,
                                             bcs_bounds,
                                             strategy)
        (pde_loss_function, bc_loss_function)
    end

    function loss_function_(θ,p)
        return pde_loss_function(θ) + bc_loss_function(θ)
    end

    f = OptimizationFunction(loss_function_, GalacticOptim.AutoZygote())
    GalacticOptim.OptimizationProblem(f, flat_initθ)
end


text:
# Julia package manager can be used to create template for our
# own package.
# We will use polgon area as an example.
#
# Start by creating a template for Package called MyPoly
#
using Pkg
Pkg.generate("MyPoly")
cd("MyPoly")
Pkg.activate(".")
import MyPoly
MyPoly.greet()


text:
using Mocking

assume_role_1_patch = @patch function assume_role(a...; b...)
    return Dict{String, Any}("PackedPolicySize" => 8, "Credentials" => Dict{String, Any}("SecretAccessKey" => "wJalrXUtnFEMI/K7MDENG/bPxRfiCYzEXAMPLEKEY", "SessionToken" => "AQoDYXdzEPT//////////wEXAMPLEtc764bNrC9SAPBSM22wDOk4x4HIZ8j4FZTwdQWLWsKWHGBuFqwAeMicRXmxfpSPfIeoIYRqTflfKD8YUuwthAx7mSEI/qkPpKPi/kMcGdQrmGdeehM4IC1NtBmUpp2wUE8phUZampKsburEDy0KPkyQDYwT7WZ0wq5VSXDvp75YU9HFvlRd8Tx6q6fE8YQcHNVXAkiY9q6d+xo0rKwT38xVqr7ZD0u0iPPkUL64lIZbqBAz+scqKmlzm8FDrypNC9Yjc8fPOLn9FX9KSYvKTr4rvx3iSIlTJabIQwj2ICCR/oLxBA==", "AccessKeyId" => "AKIAIOSFODNN7EXAMPLE", "Expiration" => "2011-07-15T23:28:33.359Z"), "AssumedRoleUser" => Dict{String, Any}("Arn" => "arn:aws:sts::123456789012:assumed-role/demo/Bob", "AssumedRoleId" => "ARO123EXAMPLE123:Bob"))
end

assume_role_with_saml_1_patch = @patch function assume_role_with_saml(a...; b...)
    return Dict{String, Any}("PackedPolicySize" => 6, "Subject" => "SamlExample", "Credentials" => Dict{String, Any}("SecretAccessKey" => "8P+SQvWIuLnKhh8d++jpw0nNmQRBZvNEXAMPLEKEY", "SessionToken" => "IQoJb3JpZ2luX2VjEOz////////////////////wEXAMPLEtMSJHMEUCIDoKK3JH9uGQE1z0sINr5M4jk+Na8KHDcCYRVjJCZEvOAiEA3OvJGtw1EcViOleS2vhs8VdCKFJQWPQrmGdeehM4IC1NtBmUpp2wUE8phUZampKsburEDy0KPkyQDYwT7WZ0wq5VSXDvp75YU9HFvlRd8Tx6q6fE8YQcHNVXAkiY9q6d+xo0rKwT38xVqr7ZD0u0iPPkUL64lIZbqBAz+scqKmlzm8FDrypNC9Yjc8fPOLn9FX9KSYvKTr4rvx3iSIlTJabIQwj2ICCR/oLxBA==", "AccessKeyId" => "ASIAV3ZUEFP6EXAMPLE", "Expiration" => "2019-11-01T20:26:47Z"), "AssumedRoleUser" => Dict{String, Any}("Arn" => "arn:aws:sts::123456789012:assumed-role/TestSaml", "AssumedRoleId" => "ARO456EXAMPLE789:TestSaml"), "SubjectType" => "transient", "Audience" => "https://signin.aws.amazon.com/saml", "Issuer" => "https://integ.example.com/idp/shibboleth", "NameQualifier" => "SbdGOnUkh1i4+EXAMPLExL/jEvs=")
end

assume_role_with_web_identity_1_patch = @patch function assume_role_with_web_identity(a...; b...)
    return Dict{String, Any}("Provider" => "www.amazon.com", "PackedPolicySize" => 123, "Credentials" => Dict{String, Any}("SecretAccessKey" => "wJalrXUtnFEMI/K7MDENG/bPxRfiCYzEXAMPLEKEY", "SessionToken" => "AQoDYXdzEE0a8ANXXXXXXXXNO1ewxE5TijQyp+IEXAMPLE", "AccessKeyId" => "AKIAIOSFODNN7EXAMPLE", "Expiration" => "2014-10-24T23:00:23Z"), "AssumedRoleUser" => Dict{String, Any}("Arn" => "arn:aws:sts::123456789012:assumed-role/FederatedWebIdentityRole/app1", "AssumedRoleId" => "AROACLKWSDQRAOEXAMPLE:app1"), "Audience" => "client.5498841531868486423.1548@apps.example.com", "SubjectFromWebIdentityToken" => "amzn1.account.AF6RHO7KZU5XRVQJGXK6HEXAMPLE")
end

decode_authorization_message_1_patch = @patch function decode_authorization_message(a...; b...)
    return Dict{String, Any}("DecodedMessage" => "{\"allowed\": \"false\",\"explicitDeny\": \"false\",\"matchedStatements\": \"\",\"failures\": \"\",\"context\": {\"principal\": {\"id\": \"AIDACKCEVSQ6C2EXAMPLE\",\"name\": \"Bob\",\"arn\": \"arn:aws:iam::123456789012:user/Bob\"},\"action\": \"ec2:StopInstances\",\"resource\": \"arn:aws:ec2:us-east-1:123456789012:instance/i-dd01c9bd\",\"conditions\": [{\"item\": {\"key\": \"ec2:Tenancy\",\"values\": [\"default\"]},{\"item\": {\"key\": \"ec2:ResourceTag/elasticbeanstalk:environment-name\",\"values\": [\"Default-Environment\"]}},(Additional items ...)]}}")
end

get_caller_identity_1_patch = @patch function get_caller_identity(a...; b...)
    return Dict{String, Any}("Arn" => "arn:aws:iam::123456789012:user/Alice", "UserId" => "AKIAI44QH8DHBEXAMPLE", "Account" => "123456789012")
end

get_caller_identity_2_patch = @patch function get_caller_identity(a...; b...)
    return Dict{String, Any}("Arn" => "arn:aws:sts::123456789012:assumed-role/my-role-name/my-role-session-name", "UserId" => "AKIAI44QH8DHBEXAMPLE:my-role-session-name", "Account" => "123456789012")
end

get_caller_identity_3_patch = @patch function get_caller_identity(a...; b...)
    return Dict{String, Any}("Arn" => "arn:aws:sts::123456789012:federated-user/my-federated-user-name", "UserId" => "123456789012:my-federated-user-name", "Account" => "123456789012")
end

get_federation_token_1_patch = @patch function get_federation_token(a...; b...)
    return Dict{String, Any}("FederatedUser" => Dict{String, Any}("Arn" => "arn:aws:sts::123456789012:federated-user/Bob", "FederatedUserId" => "123456789012:Bob"), "PackedPolicySize" => 8, "Credentials" => Dict{String, Any}("SecretAccessKey" => "wJalrXUtnFEMI/K7MDENG/bPxRfiCYzEXAMPLEKEY", "SessionToken" => "AQoDYXdzEPT//////////wEXAMPLEtc764bNrC9SAPBSM22wDOk4x4HIZ8j4FZTwdQWLWsKWHGBuFqwAeMicRXmxfpSPfIeoIYRqTflfKD8YUuwthAx7mSEI/qkPpKPi/kMcGdQrmGdeehM4IC1NtBmUpp2wUE8phUZampKsburEDy0KPkyQDYwT7WZ0wq5VSXDvp75YU9HFvlRd8Tx6q6fE8YQcHNVXAkiY9q6d+xo0rKwT38xVqr7ZD0u0iPPkUL64lIZbqBAz+scqKmlzm8FDrypNC9Yjc8fPOLn9FX9KSYvKTr4rvx3iSIlTJabIQwj2ICCR/oLxBA==", "AccessKeyId" => "AKIAIOSFODNN7EXAMPLE", "Expiration" => "2011-07-15T23:28:33.359Z"))
end

get_session_token_1_patch = @patch function get_session_token(a...; b...)
    return Dict{String, Any}("Credentials" => Dict{String, Any}("SecretAccessKey" => "wJalrXUtnFEMI/K7MDENG/bPxRfiCYzEXAMPLEKEY", "SessionToken" => "AQoEXAMPLEH4aoAH0gNCAPyJxz4BlCFFxWNE1OPTgk5TthT+FvwqnKwRcOIfrRh3c/LTo6UDdyJwOOvEVPvLXCrrrUtdnniCEXAMPLE/IvU1dYUg2RVAJBanLiHb4IgRmpRV3zrkuWJOgQs8IZZaIv2BXIa2R4OlgkBN9bkUDNCJiBeb/AXlzBBko7b15fjrBs2+cTQtpZ3CYWFXG8C5zqx37wnOE49mRl/+OtkIKGO7fAE", "AccessKeyId" => "AKIAIOSFODNN7EXAMPLE", "Expiration" => "2011-07-11T19:55:29.611Z"))
end


text:
export assemble!

"The abstract supertype of all BEM matrix operators."
abstract type BEMOperator{T} <: DictionaryOperator{T} end

kernel(A::BEMOperator) = kernel(integraloperator(A))
kernelfunction(A::BEMOperator) = kernelfunction(integraloperator(A))
measure(A::BEMOperator) = measure(integraloperator(A))
singularity(A::BEMOperator) = singularity(integraloperator(A))

"A dense BEM operator stores the dense matrix."
mutable struct DenseBEMOperator{T} <: BEMOperator{T}
    src         ::  Dictionary
    sampling    ::  SamplingOperator
    intop       ::  IntegralOperator
    quad        ::  QuadratureStrategy
    A           ::  Array{T,2}
    isassembled ::  Bool
end

initvalue(::Type{T}) where {T <: AbstractFloat} = floatmax(T)
initvalue(::Type{Complex{T}}) where {T <: AbstractFloat} = initvalue(real(T))

function DenseBEMOperator(dict::Dictionary, sampling::SamplingOperator,
        intop::IntegralOperator, quad = QuadAdaptive())
    T = promote_type(coefficienttype(dict), codomaintype(dest(sampling)))
    A = fill(initvalue(T), size(sampling)[1],length(dict))
    DenseBEMOperator{T}(dict, sampling, intop, quad, A, false)
end

samplingoperator(A::DenseBEMOperator) = A.sampling
integraloperator(A::DenseBEMOperator) = A.intop

BasisFunctions.dest(op::DenseBEMOperator) = dest(op.sampling)

isassembled(A::DenseBEMOperator{T}, i, j) where {T} =
    !(A.A[i,j] == initvalue(T))
isassembled(A::DenseBEMOperator) = A.isassembled == true

Base.:*(op::SampledIntegralOperator, dict::Dictionary) =
    DenseBEMOperator(dict, samplingoperator(op), integraloperator(op))

"Assemble the dense BEM matrix."
function assemble!(op::DenseBEMOperator, quad = op.quad, A = op.A; verbose = false, parallel = true)
    if parallel==true
		nthreads=Base.Threads.nthreads()
		verbose && println("Assembly of BEM matrix with $nthreads threads")
	    Threads.@threads for i in 1:size(A, 1)
			verbose && println("Assembly: row $i")
			for j in 1:size(A,2)
		    	A[i,j] = compute_BEM_entry(op, i, j, quad)
			end
	    end
	else
	    for i in 1:size(A, 1)
			verbose && println("Assembly: row $i")
			for j in 1:size(A,2)
		    	A[i,j] = compute_BEM_entry(op, i, j, quad)
			end
	    end
	end
	op.isassembled = true
	A
end

BasisFunctions.matrix(A::DenseBEMOperator) =
    isassembled(A) ? A.A : assemble!(A)

# We cache the computed entry, because why not, the memory is allocated anyway
Base.getindex(A::DenseBEMOperator, i::Int, j::Int) =
    isassembled(A, i, j) ? A.A[i,j] : A.A[i,j] = compute_BEM_entry(A, i, j)

function Base.getindex(A::DenseBEMOperator, i, j)
    if !isassembled(A)
        compute_BEM_range!(A.A, A, i, j)
    end
    A.A[i,j]
end


text:
using PlotlyJS

function linescatter1()
    trace1 = scatter(;x=1:4, y=[10, 15, 13, 17], mode="markers")
    trace2 = scatter(;x=2:5, y=[16, 5, 11, 9], mode="lines")
    trace3 = scatter(;x=1:4, y=[12, 9, 15, 12], mode="lines+markers")
    plot([trace1, trace2, trace3])
end

function linescatter2()
    trace1 = scatter(;x=1:5, y=[1, 6, 3, 6, 1],
                      mode="markers", name="Team A",
                      text=["A-1", "A-2", "A-3", "A-4", "A-5"],
                      marker_size=12)

    trace2 = scatter(;x=1:5+0.5, y=[4, 1, 7, 1, 4],
                      mode="markers", name= "Team B",
                      text=["B-a", "B-b", "B-c", "B-d", "B-e"])
    # setting marker.size this way is _equivalent_ to what we did for trace1
    trace2["marker"] = Dict(:size => 12)

    data = [trace1, trace2]
    layout = Layout(;title="Data Labels Hover", xaxis_range=[0.75, 5.25],
                     yaxis_range=[0, 8])
    plot(data, layout)
end

function linescatter3()
    trace1 = scatter(;x=1:5, y=[1, 6, 3, 6, 1],
                      mode="markers+text", name="Team A",
                      textposition="top center",
                      text=["A-1", "A-2", "A-3", "A-4", "A-5"],
                      marker_size=12, textfont_family="Raleway, sans-serif")

    trace2 = scatter(;x=1:5+0.5, y=[4, 1, 7, 1, 4],
                      mode="markers+text", name= "Team B",
                      textposition="bottom center",
                      text= ["B-a", "B-b", "B-c", "B-d", "B-e"],
                      marker_size=12, textfont_family="Times New Roman")

    data = [trace1, trace2]

    layout = Layout(;title="Data Labels on the Plot", xaxis_range=[0.75, 5.25],
                     yaxis_range=[0, 8], legend_y=0.5, legend_yref="paper",
                     legend=attr(family="Arial, sans-serif", size=20,
                                 color="grey"))
    plot(data, layout)
end

function linescatter4()
    trace1 = scatter(;y=fill(5, 40), mode="markers", marker_size=40,
                      marker_color=0:39)
    layout = Layout(title="Scatter Plot with a Color Dimension")
    plot(trace1, layout)
end

function linescatter5()

    country = ["Switzerland (2011)", "Chile (2013)", "Japan (2014)",
               "United States (2012)", "Slovenia (2014)", "Canada (2011)",
               "Poland (2010)", "Estonia (2015)", "Luxembourg (2013)",
               "Portugal (2011)"]

    votingPop = [40, 45.7, 52, 53.6, 54.1, 54.2, 54.5, 54.7, 55.1, 56.6]
    regVoters = [49.1, 42, 52.7, 84.3, 51.7, 61.1, 55.3, 64.2, 91.1, 58.9]

    # notice use of `attr` function to make nested attributes
    trace1 = scatter(;x=votingPop, y=country, mode="markers",
                      name="Percent of estimated voting age population",
                      marker=attr(color="rgba(156, 165, 196, 0.95)",
                                  line_color="rgba(156, 165, 196, 1.0)",
                                  line_width=1, size=16, symbol="circle"))

    trace2 = scatter(;x=regVoters, y=country, mode="markers",
                      name="Percent of estimated registered voters")
    # also could have set the marker props above by using a dict
    trace2["marker"] = Dict(:color => "rgba(204, 204, 204, 0.95)",
                           :line => Dict(:color=> "rgba(217, 217, 217, 1.0)",
                                         :width=> 1),
                           :symbol => "circle",
                           :size => 16)

    data = [trace1, trace2]
    layout = Layout(Dict{Symbol,Any}(:paper_bgcolor => "rgb(254, 247, 234)",
                                     :plot_bgcolor => "rgb(254, 247, 234)");
                    title="Votes cast for ten lowest voting age population in OECD countries",
                    width=600, height=600, hovermode="closest",
                    margin=Dict(:l => 140, :r => 40, :b => 50, :t => 80),
                    xaxis=attr(showgrid=false, showline=true,
                               linecolor="rgb(102, 102, 102)",
                               titlefont_color="rgb(204, 204, 204)",
                               tickfont_color="rgb(102, 102, 102)",
                               autotick=false, dtick=10, ticks="outside",
                               tickcolor="rgb(102, 102, 102)"),
                    legend=attr(font_size=10, yanchor="middle",
                                xanchor="right"),
                    )
    plot(data, layout)
end

function linescatter6()
    trace1 = scatter(;x=[52698, 43117], y=[53, 31],
                      mode="markers",
                      name="North America",
                      text=["United States", "Canada"],
                      marker=attr(color="rgb(164, 194, 244)", size=12,
                                  line=attr(color="white", width=0.5))
                      )

    trace2 = scatter(;x=[39317, 37236, 35650, 30066, 29570, 27159, 23557, 21046, 18007],
                      y=[33, 20, 13, 19, 27, 19, 49, 44, 38],
                      mode="markers", name="Europe",
                      marker_size=12, marker_color="rgb(255, 217, 102)",
                      text=["Germany", "Britain", "France", "Spain", "Italy",
                            "Czech Rep.", "Greece", "Poland", "Portugal"])

    trace3 = scatter(;x=[42952, 37037, 33106, 17478, 9813, 5253, 4692, 3899],
                      y=[23, 42, 54, 89, 14, 99, 93, 70],
                      mode="markers",
                      name="Asia/Pacific",
                      marker_size=12, marker_color="rgb(234, 153, 153)",
                      text=["Australia", "Japan", "South Korea", "Malaysia",
                            "China", "Indonesia", "Philippines", "India"])

    trace4 = scatter(;x=[19097, 18601, 15595, 13546, 12026, 7434, 5419],
                      y=[43, 47, 56, 80, 86, 93, 80],
                      mode="markers", name="Latin America",
                      marker_size=12, marker_color="rgb(142, 124, 195)",
                      text=["Chile", "Argentina", "Mexico", "Venezuela",
                            "Venezuela", "El Salvador", "Bolivia"])

    data = [trace1, trace2, trace3, trace4]

    layout = Layout(;title="Quarter 1 Growth",
                     xaxis=attr(title="GDP per Capital", showgrid=false, zeroline=false),
                     yaxis=attr(title="Percent", zeroline=false))

    plot(data, layout)
end

function batman()
    # reference: https://github.com/alanedelman/18.337_2015/blob/master/Lecture01_0909/The%20Bat%20Curve.ipynb
    σ(x) = √(1-x.^2)
    el(x) = 3*σ(x/7)
    s(x) = 4.2 - 0.5*x - 2.0*σ(0.5*x-0.5)
    b(x) = σ(abs(2-x)-1) - x.^2/11 + 0.5x - 3
    c(x) = [1.7, 1.7, 2.6, 0.9]

    p(i, f; kwargs...) = scatter(;x=[-i; 0.0; i], y=[f(i); NaN; f(i)],
                                  marker_color="black", showlegend=false,
                                  kwargs...)
    traces = vcat(p(3:0.1:7, el; name="wings 1"),
                  p(4:0.1:7, t->-el(t); name="wings 2"),
                  p(1:0.1:3, s; name="Shoulders"),
                  p(0:0.1:4, b; name="Bottom"),
                  p([0, 0.5, 0.8, 1], c; name="head"))

    plot(traces, Layout(title="Batman"))
end

function dumbell()
    # reference: https://plot.ly/r/dumbbell-plots/
    @eval using DataFrames

    # read Data into dataframe
    nm = tempname()
    url = "https://raw.githubusercontent.com/plotly/datasets/master/school_earnings.csv"
    download(url, nm)
    df = readtable(nm)
    rm(nm)

    # sort dataframe by male earnings
    df = sort(df, cols=[:Men], rev=false)

    men = scatter(;y=df[:School], x=df[:Men], mode="markers", name="Men",
                   marker=attr(color="blue", size=12))
    women = scatter(;y=df[:School], x=df[:Women], mode="markers", name="Women",
                     marker=attr(color="pink", size=12))

    lines = map(eachrow(df)) do r
        scatter(y=fill(r[:School], 2), x=[r[:Women], r[:Men]], mode="lines",
                name=r[:School], showlegend=false, line_color="gray")
    end

    data = Base.typed_vcat(GenericTrace, men, women, lines)
    layout = Layout(width=650, height=650, margin_l=100, yaxis_title="School",
                    xaxis_title="Annual Salary (thousands)",
                    title="Gender earnings disparity")

    plot(data, layout)
end

function errorbars1()
    trace1 = scatter(;x=vcat(1:10, 10:-1:1),
                     y=vcat(2:11, 9:-1:0),
                     fill="tozerox",
                     fillcolor="rgba(0, 100, 80, 0.2)",
                     line_color="transparent",
                     name="Fair",
                     showlegend=false)

    trace2 = scatter(;x=vcat(1:10, 10:-1:1),
                     y=[5.5, 3.0, 5.5, 8.0, 6.0, 3.0, 8.0, 5.0, 6.0, 5.5, 4.75,
                        5.0, 4.0, 7.0, 2.0, 4.0, 7.0, 4.4, 2.0, 4.5],
                     fill="tozerox",
                     fillcolor="rgba(0, 176, 246, 0.2)",
                     line_color="transparent",
                     name="Premium",
                     showlegend=false)

    trace3 = scatter(;x=vcat(1:10, 10:-1:1),
                     y=[11.0, 9.0, 7.0, 5.0, 3.0, 1.0, 3.0, 5.0, 3.0, 1.0,
                        -1.0, 1.0, 3.0, 1.0, -0.5, 1.0, 3.0, 5.0, 7.0, 9.],
                     fill="tozerox",
                     fillcolor="rgba(231, 107, 243, 0.2)",
                     line_color="transparent",
                     name="Fair",
                     showlegend=false)

    trace4 = scatter(;x=1:10, y=1:10,
                     line_color="rgb(00, 100, 80)",
                     mode="lines",
                     name="Fair")

    trace5 = scatter(;x=1:10,
                     y=[5.0, 2.5, 5.0, 7.5, 5.0, 2.5, 7.5, 4.5, 5.5, 5.],
                     line_color="rgb(0, 176, 246)",
                     mode="lines",
                     name="Premium")

    trace6 = scatter(;x=1:10, y=vcat(10:-2:0, [2, 4,2, 0]),
                     line_color="rgb(231, 107, 243)",
                     mode="lines",
                     name="Ideal")
    data = [trace1, trace2, trace3, trace4, trace5, trace6]
    layout = Layout(;paper_bgcolor="rgb(255, 255, 255)",
                    plot_bgcolor="rgb(229, 229, 229)",

                    xaxis=attr(gridcolor="rgb(255, 255, 255)",
                               range=[1, 10],
                               showgrid=true,
                               showline=false,
                               showticklabels=true,
                               tickcolor="rgb(127, 127, 127)",
                               ticks="outside",
                               zeroline=false),

                    yaxis=attr(gridcolor="rgb(255, 255, 255)",
                               showgrid=true,
                               showline=false,
                               showticklabels=true,
                               tickcolor="rgb(127, 127, 127)",
                               ticks="outside",
                               zeroline=false))

    plot(data, layout)
end

function errorbars2()
    function random_dates(d1::DateTime, d2::DateTime, n::Int)
        map(Date, sort!(rand(d1:Dates.Hour(12):d2, n)))
    end

    function _random_number(num, mul)
        value = []
        j = 0
        rand = 0
        while j <= num+1
            rand = rand() * mul
            append!(value, [rand])
            j += 1
        end
        return value
    end

    dates = random_dates(DateTime(2001, 1, 1), DateTime(2005, 12, 31), 50)

    trace1 = scatter(;x=dates,
                     y=20.0 .* rand(50),
                     line_width=0,
                     marker_color="444",
                     mode="lines",
                     name="Lower Bound")

    trace2 = scatter(;x=dates,
                     y=21.0 .* rand(50),
                     fill="tonexty",
                     fillcolor="rgba(68, 68, 68, 0.3)",
                     line_color="rgb(31, 119, 180)",
                     mode="lines",
                     name="Measurement")

    trace3 = scatter(;x=dates,
                     y=22.0 .* rand(50),
                     fill="tonexty",
                     fillcolor="rgba(68, 68, 68, 0.3)",
                     line_width=0,
                     marker_color="444",
                     mode="lines",
                     name="Upper Bound")

    data = [trace1, trace2, trace3]
    t = "Continuous, variable value error bars<br> Notice the hover text!"
    layout = Layout(;title=t, yaxis_title="Wind speed (m/s)")
    plot(data, layout)
end


text:
# include("types.jl")
include("util.jl")
# include("IntersectTest.jl")

using Sobol: SobolSeq, next!
using JSON

function slice2d(mesh::ConvexMesh; n::UInt32=UInt32(50)) :: SliceSet2D
  dims = size(mesh)[2]

  # Generate the focus points
  lbs = [b[1] for b = mesh.problemSpec.vals]
  ubs = [b[2] for b = mesh.problemSpec.vals]
  s = SobolSeq(lbs, ubs)
  fps = collect([next!(s) for _ = 1:n])

  hyperslice(mesh, fps)
end

function slice2d(mesh::ConvexMesh, fps::Array{PointND}) :: SliceSet2D
  slices = []
  for fp in fps
    ss = sliceDims(mesh, fp)
    append!(slices, ss)
  end
  SliceSet2D(mesh.problemSpec, slices)
end

# slices of all dimensions
function sliceDims(mesh::ConvexMesh, fp::PointND)
  slices = []
  for dd in dimComb(UInt64(size(mesh)[2]))
    ss = slice(mesh, fp, dd[1], dd[2])
    hss = [HypersliceSegment(fp, simp, dd[1], dd[2], s.p1d1, s.p1d2, s.p2d1, s.p2d2) for (simp, s) in ss]
    append!(slices, hss)
  end
  return slices
end

function slice(mesh::ConvexMesh, fp::PointND, d1::Dim, d2::Dim)::Vector{Tuple{Simplex,Intersect2D}}
  slices = Tuple{Simplex,Intersect2D}[]
  for s in mesh
    ts = [(s, t) for t in simplexPointIntersection(s, fp, d1, d2)]
    append!(slices, ts)
  end
  return slices
end

function hypersliceplorer(mesh::ConvexMesh; n::UInt32=UInt32(50))
  slices = slice2d(mesh; n=n);
  plot(slices)
end

function hypersliceplorer(mesh::ConvexMesh, fps::Array{PointND})
  slices = slice2d(mesh, fps);
  plot(slices)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "TDLib"
version = v"1.6.10"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/tdlib/td.git", "78e45da6b9014f8efd96cee0ad5701464d911cc7")
]

# Bash recipe for building across all platforms
script = raw"""
apk add gperf
cd td/
install_license LICENSE_1_0.txt
sed -i 's|TD_HAS_MMSG 1|TD_HAS_MMSG 0|' tdutils/td/utils/port/config.h  # otherwise not compatible with old glibc
mkdir build
cd build
cmake -DCMAKE_BUILD_TYPE=Release -DZLIB_LIBRARY=$libdir -DZLIB_INCLUDE_DIR=$includedir ..
for f in $(grep -l $(find . -name link.txt) -e 'opt/bin/c++')
do
    # without this flag it cannot find libz
    sed -i ' 1 s/$/ -lz/' $f
done
cmake --build . -j${nproc}
cp libtdjson.${dlext} $libdir/libtdjson.${dlext}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter(p -> Sys.islinux(p) && (libc(p) == "glibc" && arch(p) == "i686" || arch(p) == "x86_64"), supported_platforms())
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libtdjson", :libtdjson)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency.(["OpenSSL_jll", "Zlib_jll"])

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil