text:
using GUI: googlefont, googlefontlink

@testset "Google Fonts" begin
  @test googlefont("Roboto") === "Roboto"
  @test googlefont("IBM Flex Mono") === "IBM+Flex+Mono"
  @test googlefont(["Roboto", "IBM Flex Mono"]) === "Roboto|IBM+Flex+Mono"
  @test googlefont("Roboto", [400]) === "Roboto:400"
  @test googlefont("Roboto", [400,700]) === "Roboto:400,700"
  @test googlefont("Roboto", [400], ["","i"]) === "Roboto:400,400i"
  @test googlefont("Roboto", [400,700], ["","i"]) == "Roboto:400,400i,700,700i"
end

@testset "Google Font Links" begin
  @test googlefontlink("Roboto") === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto\", rel=\"stylesheet\">" ||
        googlefontlink("Roboto") === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto\", rel=\"stylesheet\"/>" ||
        googlefontlink("Roboto") === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto\", rel=\"stylesheet\"></link>"
  @test googlefontlink("IBM Flex Mono") === "<link ref=\"https://fonts.googleapis.com/css?family=IBM+Flex+Mono\", rel=\"stylesheet\">" ||
        googlefontlink("IBM Flex Mono") === "<link ref=\"https://fonts.googleapis.com/css?family=IBM+Flex+Mono\", rel=\"stylesheet\"/>" ||
        googlefontlink("IBM Flex Mono") === "<link ref=\"https://fonts.googleapis.com/css?family=IBM+Flex+Mono\", rel=\"stylesheet\"></link>"
  @test googlefontlink(["Roboto", "IBM Flex Mono"]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto|IBM+Flex+Mono\", rel=\"stylesheet\">" ||
        googlefontlink(["Roboto", "IBM Flex Mono"]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto|IBM+Flex+Mono\", rel=\"stylesheet\"/>" ||
        googlefontlink(["Roboto", "IBM Flex Mono"]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto|IBM+Flex+Mono\", rel=\"stylesheet\"></link>"
  @test googlefontlink("Roboto", weight=[400]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400\", rel=\"stylesheet\">" ||
        googlefontlink("Roboto", weight=[400]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400\", rel=\"stylesheet\"/>" ||
        googlefontlink("Roboto", weight=[400]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400\", rel=\"stylesheet\"></link>"
  @test googlefontlink("Roboto", weight=[400,700]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400,700\", rel=\"stylesheet\">" ||
        googlefontlink("Roboto", weight=[400,700]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400,700\", rel=\"stylesheet\"/>" ||
        googlefontlink("Roboto", weight=[400,700]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400,700\", rel=\"stylesheet\"></link>"
  @test googlefontlink("Roboto", weight=[400], style=["","i"]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400,400i\", rel=\"stylesheet\">" ||
        googlefontlink("Roboto", weight=[400], style=["","i"]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400,400i\", rel=\"stylesheet\"/>" ||
        googlefontlink("Roboto", weight=[400], style=["","i"]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400,400i\", rel=\"stylesheet\"></link>"
  @test googlefontlink("Roboto", weight=[400,700], style=["","i"]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400,400i,700,700i\", rel=\"stylesheet\">" ||
        googlefontlink("Roboto", weight=[400,700], style=["","i"]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400,400i,700,700i\", rel=\"stylesheet\"/>" ||
        googlefontlink("Roboto", weight=[400,700], style=["","i"]) === "<link ref=\"https://fonts.googleapis.com/css?family=Roboto:400,400i,700,700i\", rel=\"stylesheet\"></link>"

end


text:
# Define MyType
struct MyType  # Don't inherit from anything
    value
end
value(x::MyType) = x.value

# Define its methods
Base.one(::Type{MyType}) = MyType(1)
Base.zero(::Type{MyType}) = MyType(0)

Base.:+(x::MyType, y::MyType) = MyType(value(x)+value(y))
Base.:*(x::MyType, y::MyType) = MyType(value(x)*value(y))
Base.:-(x::MyType, y::MyType) = MyType(value(x)-value(y))
Base.:/(x::MyType, y::MyType) = MyType(value(x)/value(y))

Base.adjoint(x::MyType) = MyType(adjoint(value(x)))

function Base.isapprox(x::MyType, y::MyType; atol, rtol)
    return isapprox(value(x), value(y), atol=atol, rtol=rtol)
end
function Base.rtoldefault(::Type{MyType}, ::Type{MyType}, z)
    # I don't really want to be restricted to Float64 tolerance.
    return Base.rtoldefault(Float64, Float64, z)
end

# Shouldn't have to define these
Base.one(::MyType) = one(MyType)
Base.zero(::MyType) = zero(MyType)

Base.:+(x::MyType, y) = MyType(value(x)+y)
Base.:*(x::MyType, y) = MyType(value(x)*y)

Base.abs(x::MyType) = MyType(abs(value(x)))
Base.:<(x::MyType, y::MyType) = value(x) < value(y)
Base.inv(x::MyType) = MyType(inv(value(x)))

Base.promote_rule(::Type{Any}, ::Type{MyType}) = MyType
Base.convert(::Type{MyType}, x::MyType) = x
Base.convert(::Type{MyType}, x) = MyType(x)

Base.iterate(x::MyType) = (value(x), nothing)
Base.iterate(::MyType, ::Any) = nothing
Base.length(x::MyType) = 1

# Begin check for ((X-I)')^-1 ≈ Y
X = [
    0.4 1;
    -0.2 0.8;
]
X = map(x -> MyType(x), X)

Y = [
    -0.625 0.625;
    -3.125 -1.875;
]
Y = map(x -> MyType(x), Y)

using LinearAlgebra
println((X-I)')
println(inv((X-I)'))
println(inv((X-I)') ≈ Y)

dataframe = [
    1 2;
    1 3;
    2 3;
    2 4;
    2 5;
    3 2;
    4 2;
    5 3;
]

max_trial = maximum(dataframe[:, 1])
max_c1 = Array{Any}(undef, max_trial)
for c in 1:max_trial
    c1 = findall(dataframe[:, 1] .== c)
    max_c1[c] = maximum(dataframe[c1, 2])
end

max_c1 = [maximum(dataframe[c1, 2]) for c in unique()]

# Name of trials
naming = unique(dataframe[!,:9])
dataframe = dataframe[!,:1:8]
# Number of trials
max_trial = maximum(dataframe[!,:1])
# Number of variables on each trial
max_c1 = []
for c in 1:max_trial
    c1 = findall(dataframe[:,1] .== c)
    max_c1[c] = maximum(dataframe[c1,2])
end

ext: jl


text:
struct Solver{T}
    p::Problem{T}
    nlp_bounds::Vector{MOI.NLPBoundsPair}
    block_data::MOI.NLPBlockData
    solver::Ipopt.Optimizer
    z::Vector{MOI.VariableIndex}
end

function Solver(trajopt::TrajectoryOptimizationProblem; eval_hess=false, options=nothing)
    p = Problem(trajopt, eval_hess=eval_hess)

    nlp_bounds = MOI.NLPBoundsPair.(p.con_bnds...)
    block_data = MOI.NLPBlockData(nlp_bounds, p, true)

    solver = Ipopt.Optimizer()

    # TODO: options
    solver.options["max_iter"] = 1000
    solver.options["tol"] = 1.0e-3
    solver.options["constr_viol_tol"] = 1.0e-3

    z = MOI.add_variables(solver, p.num_var)

    for i = 1:p.num_var
        MOI.add_constraint(solver, z[i], MOI.LessThan(p.var_bnds[2][i]))
        MOI.add_constraint(solver, z[i], MOI.GreaterThan(p.var_bnds[1][i]))
    end

    MOI.set(solver, MOI.NLPBlock(), block_data)
    MOI.set(solver, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return Solver(p, nlp_bounds, block_data, solver, z)
end

function initialize!(s::Solver, z)
    for i = 1:s.p.num_var
        MOI.set(s.solver, MOI.VariablePrimalStart(), s.z[i], z[i])
    end
end

function solve!(s::Solver)
    MOI.optimize!(s.solver)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((3, 4, true, true, true, true),),
    ((3, 4, false, true, true, true),),
    ((3, 4, true, false, true, true),),
    ((3, 4, false, false, true, true),),
    ((3, 4, true, true, false, true),),
    ((3, 4, false, false, false, true),),
    ((3, 4, true, true, false, false),),
    ((3, 4, false, false, false, false),),
    ]
insts["fast"] = [
    ((12, 24, true, true, true, true),),
    ((12, 24, false, true, true, true),),
    ((12, 24, true, false, true, true),),
    ((12, 24, false, false, true, true),),
    ((12, 24, true, true, false, true),),
    ((12, 24, false, false, false, true),),
    ((12, 24, true, true, false, false),),
    ((12, 24, false, false, false, false),),
    ]
insts["various"] = [
    ((5, 100, true, true, true, true),),
    ((5, 100, true, true, false, true),),
    ((10, 200, true, true, true, true),),
    ((10, 200, true, true, false, true),),
    ]
return (MatrixCompletionNative, insts)


text:
using BDF, Base.Test

#test split BDF at trigger

origFilePath = joinpath(dirname(@__FILE__), "Newtest17-256.bdf")
copyFilePath = joinpath(dirname(@__FILE__), "Newtest17-256_midtrig.bdf")
copyFilePath1 = joinpath(dirname(@__FILE__), "Newtest17-256_midtrig_1.bdf")
copyFilePath2 = joinpath(dirname(@__FILE__), "Newtest17-256_midtrig_2.bdf")
copyFilePathComb = joinpath(dirname(@__FILE__), "Newtest17-256_comb.bdf")

bdfHeader = readBDFHeader(origFilePath)
dats, evtTab, trigs, statusChan = readBDF(origFilePath)

trigs[bdfHeader["sampRate"][1]*30] = 200

writeBDF(copyFilePath, dats, trigs, statusChan, bdfHeader["sampRate"][1],
         subjID=bdfHeader["subjID"], recID=bdfHeader["recID"], startDate=bdfHeader["startDate"],
         startTime=bdfHeader["startTime"], versionDataFormat=bdfHeader["versionDataFormat"],
                  chanLabels=bdfHeader["chanLabels"][1:end-1], transducer=bdfHeader["transducer"][1:end-1],
                  physDim=bdfHeader["physDim"][1:end-1], physMin=bdfHeader["physMin"][1:end-1], physMax=bdfHeader["physMax"][1:end-1],
                  prefilt=bdfHeader["prefilt"][1:end-1])

splitBDFAtTrigger(copyFilePath, 200)
dats1, evtTab1, trigs1, statusChan1 = readBDF(copyFilePath1)
dats2, evtTab2, trigs2, statusChan2 = readBDF(copyFilePath2)

writeBDF(copyFilePathComb, [dats1 dats2], [trigs1; trigs2], [statusChan1; statusChan2], bdfHeader["sampRate"][1])
datsCopy, evtTabCopy, trigsCopy, statusChanCopy = readBDF(copyFilePath)

datsComb, evtTabComb, trigsComb, statusChanComb = readBDF(copyFilePathComb)

rm(copyFilePath)
rm(copyFilePath1)
rm(copyFilePath2)
rm(copyFilePathComb)

@test isequal(datsCopy, datsComb)
@test isequal(evtTabCopy, evtTabComb)
@test isequal(trigsCopy, trigsComb)
@test isequal(statusChanCopy, statusChanComb)


text:
typealias PaTime Cdouble
typealias PaError Cint
typealias PaSampleFormat Culong
# PaStream is always used as an opaque type, so we're always dealing with the
# pointer
typealias PaStream Ptr{Void}
typealias PaDeviceIndex Cint
typealias PaHostApiIndex Cint
typealias PaTime Cdouble
typealias PaHostApiTypeId Cint
typealias PaStreamCallback Void

const PA_NO_ERROR = 0
const PA_INPUT_OVERFLOWED = -10000 + 19
const PA_OUTPUT_UNDERFLOWED = -10000 + 20

const paFloat32 = convert(PaSampleFormat, 0x01)
const paInt32   = convert(PaSampleFormat, 0x02)
const paInt24   = convert(PaSampleFormat, 0x04)
const paInt16   = convert(PaSampleFormat, 0x08)
const paInt8    = convert(PaSampleFormat, 0x10)
const paUInt8   = convert(PaSampleFormat, 0x20)

# PaHostApiTypeId values
const pa_host_api_names = {
    0 => "In Development", # use while developing support for a new host API
    1 => "Direct Sound",
    2 => "MME",
    3 => "ASIO",
    4 => "Sound Manager",
    5 => "Core Audio",
    7 => "OSS",
    8 => "ALSA",
    9 => "AL",
    10 => "BeOS",
    11 => "WDMKS",
    12 => "Jack",
    13 => "WASAPI",
    14 => "AudioScience HPI"
}

# track whether we've already inited PortAudio
portaudio_inited = false

################## Types ####################

type PortAudioStream <: AudioStream
    root::AudioMixer
    info::DeviceInfo
    show_warnings::Bool
    stream::PaStream

    function PortAudioStream(sample_rate::Integer=44100,
                             buf_size::Integer=1024,
                             show_warnings::Bool=false)
        require_portaudio_init()
        stream = Pa_OpenDefaultStream(1, 1, paFloat32, sample_rate, buf_size)
        Pa_StartStream(stream)
        root = AudioMixer()
        this = new(root, DeviceInfo(sample_rate, buf_size),
                   show_warnings, stream)
        info("Scheduling PortAudio Render Task...")
        # the task will actually start running the next time the current task yields
        @schedule(portaudio_task(this))
        finalizer(this, destroy)

        this
    end
end

function destroy(stream::PortAudioStream)
    # in 0.3 we can't print from a finalizer, as STDOUT may have been GC'ed
    # already and we get a segfault. See
    # https://github.com/JuliaLang/julia/issues/6075
    #info("Cleaning up stream")
    Pa_StopStream(stream.stream)
    Pa_CloseStream(stream.stream)
    # we only have 1 stream at a time, so if we're closing out we can just
    # terminate PortAudio.
    Pa_Terminate()
    portaudio_inited = false
end

############ Internal Functions ############

function portaudio_task(stream::PortAudioStream)
    info("PortAudio Render Task Running...")
    n = bufsize(stream)
    buffer = zeros(AudioSample, n)
    try
        while true
            while Pa_GetStreamReadAvailable(stream.stream) < n
                sleep(0.005)
            end
            Pa_ReadStream(stream.stream, buffer, n, stream.show_warnings)
            # assume the root is always active
            rendered = render(stream.root.renderer, buffer, stream.info)::AudioBuf
            for i in 1:length(rendered)
                buffer[i] = rendered[i]
            end
            for i in (length(rendered)+1):n
                buffer[i] = 0.0
            end
            while Pa_GetStreamWriteAvailable(stream.stream) < n
                sleep(0.005)
            end
            Pa_WriteStream(stream.stream, buffer, n, stream.show_warnings)
        end
    catch ex
        warn("Audio Task died with exception: $ex")
        Base.show_backtrace(STDOUT, catch_backtrace())
    end
end

type PaDeviceInfo
    struct_version::Cint
    name::Ptr{Cchar}
    host_api::PaHostApiIndex
    max_input_channels::Cint
    max_output_channels::Cint
    default_low_input_latency::PaTime
    default_low_output_latency::PaTime
    default_high_input_latency::PaTime
    default_high_output_latency::PaTime
    default_sample_rate::Cdouble
end

type PaHostApiInfo
    struct_version::Cint
    api_type::PaHostApiTypeId
    name::Ptr{Cchar}
    deviceCount::Cint
    defaultInputDevice::PaDeviceIndex
    defaultOutputDevice::PaDeviceIndex
end

type PortAudioInterface <: AudioInterface
    name::String
    host_api::String
    max_input_channels::Int
    max_output_channels::Int
end

function get_portaudio_devices()
    require_portaudio_init()
    device_count = ccall((:Pa_GetDeviceCount, libportaudio), PaDeviceIndex, ())
    pa_devices = [Pa_GetDeviceInfo(i) for i in 0:(device_count - 1)]
    [PortAudioInterface(bytestring(d.name),
                        bytestring(Pa_GetHostApiInfo(d.host_api).name),
                        d.max_input_channels,
                        d.max_output_channels)
     for d in pa_devices]
end

function require_portaudio_init()
    # can be called multiple times with no effect
    global portaudio_inited
    if !portaudio_inited
        info("Initializing PortAudio. Expect errors as we scan devices")
        Pa_Initialize()
        portaudio_inited = true
    end
end

# Low-level wrappers for Portaudio calls
Pa_GetDeviceInfo(i) = unsafe_load(ccall((:Pa_GetDeviceInfo, libportaudio),
                                 Ptr{PaDeviceInfo}, (PaDeviceIndex,), i))
Pa_GetHostApiInfo(i) = unsafe_load(ccall((:Pa_GetHostApiInfo, libportaudio),
                                   Ptr{PaHostApiInfo}, (PaHostApiIndex,), i))

function Pa_Initialize()
    err = ccall((:Pa_Initialize, libportaudio), PaError, ())
    handle_status(err)
end

function Pa_Terminate()
    err = ccall((:Pa_Terminate, libportaudio), PaError, ())
    handle_status(err)
end

function Pa_StartStream(stream::PaStream)
    err = ccall((:Pa_StartStream, libportaudio), PaError,
                (PaStream,), stream)
    handle_status(err)
end

function Pa_StopStream(stream::PaStream)
    err = ccall((:Pa_StopStream, libportaudio), PaError,
                (PaStream,), stream)
    handle_status(err)
end

function Pa_CloseStream(stream::PaStream)
    err = ccall((:Pa_CloseStream, libportaudio), PaError,
                (PaStream,), stream)
    handle_status(err)
end

function Pa_GetStreamReadAvailable(stream::PaStream)
    avail = ccall((:Pa_GetStreamReadAvailable, libportaudio), Clong,
                (PaStream,), stream)
    avail >= 0 || handle_status(avail)
    avail
end

function Pa_GetStreamWriteAvailable(stream::PaStream)
    avail = ccall((:Pa_GetStreamWriteAvailable, libportaudio), Clong,
                (PaStream,), stream)
    avail >= 0 || handle_status(avail)
    avail
end

function Pa_ReadStream(stream::PaStream, buf::Array, frames::Integer=length(buf),
                       show_warnings::Bool=true)
    frames <= length(buf) || error("Need a buffer at least $frames long")
    err = ccall((:Pa_ReadStream, libportaudio), PaError,
                (PaStream, Ptr{Void}, Culong),
                stream, buf, frames)
    handle_status(err, show_warnings)
    buf
end

function Pa_WriteStream(stream::PaStream, buf::Array, frames::Integer=length(buf),
                        show_warnings::Bool=true)
    frames <= length(buf) || error("Need a buffer at least $frames long")
    err = ccall((:Pa_WriteStream, libportaudio), PaError,
                (PaStream, Ptr{Void}, Culong),
                stream, buf, frames)
    handle_status(err, show_warnings)
    nothing
end

Pa_GetVersion() = ccall((:Pa_GetVersion, libportaudio), Cint, ())

function Pa_GetVersionText()
    versionPtr = ccall((:Pa_GetVersionText, libportaudio), Ptr{Cchar}, ())
    bytestring(versionPtr)
end

function Pa_OpenDefaultStream(inChannels::Integer, outChannels::Integer,
                              sampleFormat::PaSampleFormat,
                              sampleRate::Real, framesPerBuffer::Integer)
    streamPtr::Array{PaStream} = PaStream[0]
    err = ccall((:Pa_OpenDefaultStream, libportaudio),
                PaError, (Ptr{PaStream}, Cint, Cint,
                          PaSampleFormat, Cdouble, Culong,
                          Ptr{PaStreamCallback}, Ptr{Void}),
                streamPtr, inChannels, outChannels, sampleFormat, sampleRate,
                framesPerBuffer, 0, 0)
    handle_status(err)

    streamPtr[1]
end

function handle_status(err::PaError, show_warnings::Bool=true)
    if err == PA_OUTPUT_UNDERFLOWED || err == PA_INPUT_OVERFLOWED
        if show_warnings
            msg = ccall((:Pa_GetErrorText, libportaudio),
                        Ptr{Cchar}, (PaError,), err)
            warn("libportaudio: " * bytestring(msg))
        end
    elseif err != PA_NO_ERROR
        msg = ccall((:Pa_GetErrorText, libportaudio),
                    Ptr{Cchar}, (PaError,), err)
        error("libportaudio: " * bytestring(msg))
    end
end


text:
using SimpleLockFiles

## ----------------------------------------
function _read(valfn)
    try; return parse(Int, read(valfn, String))
    catch ignored end
    return -1
end

function _write(valfn, val)
    try; write(valfn, string(val))
    catch ignored end
    return nothing
end

function _append(valfn, val)
    try; open((io) -> println(io, val), valfn, "a")
    catch ignored end
    return nothing
end

## ----------------------------------------
# OS
@async let
    while true
        println("pid: ", getpid(), ", running!")
        sleep(1.0)
    end
end

## ----------------------------------------
try
    println("Hi from ", getpid())

    # must match the values on tace_test.jl
    N = 10
    lkfn = joinpath(@__DIR__, "lock")
    valfn = joinpath(@__DIR__, "sum.txt")
    logfn = joinpath(@__DIR__, "log.txt")
    slf = SimpleLockFile(lkfn)

    val = 0
    it = 0
    t0 = time()
    frec = 0.0

    lock_kwargs = (;tout = 10.0, vtime = 1.0, wtime = 0.1, ctime = 0.5, force = false)
    while true
        lkid = string("PROC-", getpid(), "-", it)
        println("Iter init", lkid)
        ok_flag = lock(slf, lkid; lock_kwargs...) do

            touch(valfn)
            val = _read(valfn)
            val += 1
            _write(valfn, val)

            # info
            frec = time() - t0
            msg = string("pid: ", getpid(), ", it: ", it, ", val: ", val, ", frec [s]: ", frec)
            _append(logfn, msg); println(msg); flush.([stdout, stderr])
            t0 = time()
            it += 1
        end

        println("Iter finished, ok_flag: ", ok_flag)
        it == N && break
        sleep(0.1 * rand())
    end # while true

    msg = string("pid: ", getpid(), ", it: ", it, ", val: ", val, " GOOD BYE")
    _append(logfn, msg); println(msg); flush.([stdout, stderr])

catch err
    println()
    showerror(stdout, err, catch_backtrace())
    println()
finally
    flush.([stdout, stderr])
    sleep(15.0)
    exit()
end


text:
module Graphs

using Reexport

include("BasicGraphs.jl")
include("BipartiteGraphs.jl")
include("PropertyGraphs.jl")
include("GraphAlgorithms.jl")
include("GraphGenerators.jl")

@reexport using .BasicGraphs
@reexport using .BipartiteGraphs
@reexport using .PropertyGraphs
@reexport using .GraphAlgorithms
@reexport using. GraphGenerators

end


text:
# Squared Exponential Covariance Function

include("se_iso.jl")
include("se_ard.jl")

@doc """
# Description
Constructor the Squared Exponential kernel

# See also
SEIso, SEArd
""" ->
SE(ll::Float64, lσ::Float64) = SEIso(ll, lσ)
SE(ll::Vector{Float64}, lσ::Float64) = SEArd(ll, lσ)


text:
using EAGOSmoothMcCormickGrad
using Base.Test

# write your own tests here
println("Testing Utility Functions...")
t = @elapsed include("Utilities.jl")
println("done (took $t seconds).")

println("Testing Sparse Preconditioner...")
t = @elapsed include("SparseCntr.jl")
println("done (took $t seconds).")

println("Implicit Bounding Utilities...")
t = @elapsed include("Imp_Util_Test.jl")
println("done (took $t seconds).")

println("Testing McCormick Multiplication Operator...")
t = @elapsed include("Mult_Test.jl")
println("done (took $t seconds).")

println("Testing McCormick Extrema Operator...")
t = @elapsed include("Extrema_Test.jl")
println("done (took $t seconds).")

println("Testing McCormick Operators...")
t = @elapsed include("Operators.jl")
println("done (took $t seconds).")

println("Implicit Bounding Subroutines...")
t = @elapsed include("ImplicitBnd.jl")
println("done (took $t seconds).")

#=
println("Implicit Bounding Utilities...")
t = @elapsed include("D1_Interval_Test.jl")
println("done (took $t seconds).")

println("Implicit Function Routines...")
t = @elapsed include("D1_Interval_Test.jl")
println("done (took $t seconds).")
=#


text:
using Gadfly, DataArrays, RDatasets

 plot(dataset("ggplot2", "diamonds"), x=:Price, Geom.density)


text:
############################################################################################
"Default Configuration for Custom sampler"
struct ConfigCustom{F} <: AbstractConfiguration
    "Target Acceptance Rate - ignored for custom sampler."
    δ::Float64
    "Default size for tuning iterations in each cycle - ignored for custom sampler."
    window::ConfigTuningWindow
    "A function of a parameter vector θ that returns a valid Distribution to sample and evaluate from"
    proposal::F
    function ConfigCustom(
        δ::Float64,
        window::ConfigTuningWindow,
        proposal::F
    ) where {F}
        @argcheck 0.0 < δ <= 1.0 "Acceptance rate not bounded between 0 and 1"
        return new{F}(
            δ,
            window,
            proposal
        )
    end
end

############################################################################################
"""
$(SIGNATURES)
A placeholder function to configure a custom proposal distribution.
Can also be a closure that captures, i.e., Σ., or a a callable struct:

# Examples
```julia

#Mixture via callable struct
struct MyMixture{W, F, T}
    weights::W
    Σ1::F
    Σ2::T
end
function (mixture::MyMixture)(θ)
    @unpack weights, Σ1, Σ2 = mixture
    return MixtureModel([MvNormal(θ, Σ1), MvNormal(θ, Σ2)], weights)
end

#Independence sampler via callable struct without extra allocation
struct MyIndependenceSampler{W}
    dist::W
end
function (sampler::MyIndependenceSampler)(θ)
    return sampler.dist
end
```

"""
function customdefaultkernel(θ::AbstractVector{T}) where {T<:Real}
    Nparams = length(θ)
    return MvNormal(θ, T(1/Nparams^2) * LinearAlgebra.I(Nparams))
end

"""
$(SIGNATURES)
Initialize Custom custom configurations.

# Examples
```julia
```

"""
function init(
    ::Type{AbstractConfiguration},
    mcmc::Type{Custom},
    objective::Objective,
    proposalconfig::ConfigProposal;
    ## Target Acceptance Rate
    δ=0.234,
    ## MCMC Phase tune variables based on standard HMC best practice target acceptance rate of 0.80
    window= ConfigTuningWindow(
        [1, 5, 1],
        Int.(floor.([75, 50, 25] * (0.80 / 0.234))),
        [Warmup(), Adaptionˢˡᵒʷ(), Adaptionᶠᵃˢᵗ(), Exploration()],
        ),
    ## Custom distribution
    proposal = customdefaultkernel
)
    return ConfigCustom(
        δ,
        window,
        proposal
    )
end

############################################################################################
# Export
export ConfigCustom, customdefaultkernel


text:
module UiPathOrchestratorJobSchedulingPlanCreate
using LinearAlgebra,Dates,StatsBase
using JuMP,Ipopt,Cbc
using DataFrames,XLSX
ENV["GKS_ENCODING"]="utf8"
using Plots,GR

include("input.jl")
include("commonfunc.jl")
include("engine1.jl")
include("output.jl")

"""
    UiPathOrchestratorJobSchedulingPlanCreate

# 処理概要
UiPathOrchestratorのジョブスケジューリング案の作成

# 使い方
ExcelFilePath=""
uipathorchestratorschedulrecreate(ExcelFilePath,"parameters","schedule",plotengine="GR")
"""

"""
    uipathorchestratorschedulrecreate(ExcelFilePath::String,parameters::String,schedule::String;planexport::Bool=false,ExportExcelFilePath::String="",plotengine="PlotlyJS",schedulcolumn::Int=6,fontfamily="serif-roman")

# 処理概要
スケジュール作成

# 引数
* `ExcelFilePath`:スケジュール作成に必要な情報が記載されたExcelファイルのフルパスを指定。
* `parameters`:パラメーターが記載されたシート名を指定
* `schedule`:ジョブ作成に必要な情報が記載されたシート名を指定
* ` planexport`:(オプション)スケジュール作成結果を出力するかどうか。出力する場合、trueを指定。
* `ExportExcelFilePath`:(オプション)スケジュール作成先を指定。指定しない場合は、カレントディレクトリ（フォルダ）に"UiPathOrchestratorJobSchedulingPlan.xlsx"名で出力
* `plotengine`:(オプション)スケジュール調整結果を出力。指定値は、"GR","PlotlyJS","off"の３種類。offを選択するとグラフ出力しない
* `schedulcolumn`:スケジュールON,OFF開始列を指定（デフォルト：６）
* `fontfamily`:フォントを指定

# 結果（戻り値）
スケジュール調整結果を出力
"""
function uipathorchestratorschedulrecreate(ExcelFilePath::String,parameters::String,schedule::String;planexport::Bool=false,ExportExcelFilePath::String="",plotengine="PlotlyJS",schedulcolumn::Int=6,checkreturn::Bool=false,fontfamily::String="serif-roman")
  scheduleplan,robotn,run_unit_time,jobn,timen,schedule,blocktime_start,blocktime_end,blocktime_dow=readprerequisite(ExcelFilePath,parameters,schedule,schedulcolumn=schedulcolumn)
  plan,runtime=uipathorchestratorschedulreadjustment(scheduleplan,robotn,run_unit_time,jobn,timen,schedule,blocktime_start,blocktime_end,blocktime_dow,schedulcolumn=schedulcolumn)
  if(checkreturn)
    plan,adjustedresultcheckflag=adjustedresultcheck(plan,runtime,scheduleplan,robotn,jobn,timen,schedule,blocktime_start,blocktime_end,blocktime_dow,schedulcolumn=schedulcolumn,checkreturn=checkreturn)
  else
    plan=adjustedresultcheck(plan,runtime,scheduleplan,robotn,jobn,timen,schedule,blocktime_start,blocktime_end,blocktime_dow,schedulcolumn=schedulcolumn,checkreturn=checkreturn)
  end

  if(plotengine=="PlotlyJS")
    plotplanplotlyjs(plan,schedulcolumn=schedulcolumn,fontfamily="serif-roman")
  elseif(plotengine=="GR")
    plotplangr(plan,schedulcolumn=schedulcolumn,fontfamily="serif-roman")
  elseif(plotengine=="off"||plotengine=="")
  else
    println("plotengineはPlotlyJS,GR,offのどれかを選択してください")
  end

  if(planexport)
    exportplan(plan,ExcelFilePath=ExportExcelFilePath)
  end

  if(checkreturn)
    return plan,adjustedresultcheckflag
  else
    return plan
  end

end

export uipathorchestratorschedulrecreate,readprerequisite,uipathorchestratorschedulreadjustment,adjustedresultcheck,ContinuousOperation
export exportplan,plotplanmaster,plotplangr,plotplanplotlyjs
export createtimeset
end # module


text:
# # Regular Expressions [`CombinedParser`](@ref)
# Parsing is reading by a computer program,
# formally transforming character sequences
# into defined and structured representations.
#
# -- name parser slide
#
# Representation and Transformation
#
# -- number list slide
#
# So parsing is a fundamental
# obstacle before you can begin working with data.
# Writing parsers is ubiquitous.
#
# ## What is a regular expression?
# Regular expressions and state machines are the most common tools of the parsing trade.
# Most programmers and data scientists are familiar with regular expressions.
#
# A regular expression is a textual representation of a parser program.
# The syntax of regular expressions is condensed and arcane.
#
# Parsing has been a painpoint in all of my work as a developer and a researcher with a focus on text
# while I accumulated plenty of regular expressions.
# I still use these regexps at many places and do not want to rewrite all regular expressions I still need.
# So I wrote a CombinedParser that parses regex, and transforms it into an equivalent julia parser program,
using CombinedParsers

# This guide demonstrates constructing a recursive `CombinedParser` with the **[`@syntax`](@ref) or [`push!`](@ref) to [`Either`](@ref)** technique.
# This page outlines essentials of the [`@re_str`](@ref) macro
using CombinedParsers.Regexp

# ## Characters and Escaped Characters
# The simplest regular expression pattern is a character matcher.
import CombinedParsers.Regexp: meta_chars
char_matcher =  Either(
    CharNotIn(meta_chars), # character
    Sequence( # escaped_meta
	2,    # transform: emit unescaped at index
        '\\', CharIn(meta_chars)
    )
)

# The `char_matcher` unescapes regex `meta_chars`:
parse(Repeat(char_matcher),raw"a\[b"; log=true)

# ## Repeated patterns
# ### Number of Repetitions: Transforming a Parsing
# A basic feature of regular expression is whether a repeatable pattern (e.g. a char) is optional `?` (repeated 0 or 1 times), or should be repeated at least once `+`, any times `*`, or `{min,max}` times ([pcre spec](https://www.pcre.org/original/doc/html/pcrepattern.html#SEC17)).
# The syntax of regex repetitions require mapping such `SubString`s to representions as `UnitRange{Int}`.
#
#
# Above, the `:escaped_meta` mapped the parsing to emit the unescaped character at the second position in its sequence.
# Any Julia function can be used for transformations with the `do` syntax after a `CombinedParser` constructor to implicitly call [`map`](@ref) creating a [`Transformation`](@ref) parser.

import CombinedParsers.Regexp: integer, Repeat_max
@with_names repetitions = Either(
    '?' => 0:1,
    '+' => 1:Repeat_max,
    '*' => 0:Repeat_max,
    Sequence(
        '{',
        :min => integer,
        :max => Optional(
	    Sequence(2,
		     ',', integer | Repeat_max)),
        '}'
    ) do v
    if v.max isa Missing
    v.min:v.min
    else
    v.min:v.max
    end::UnitRange{Int}
    end
)

# `repetitions` parses PCRE repetitions syntax into a `(min, max)`:
parse(repetitions, "?")
parse(repetitions, "*")
parse(repetitions, "+")
parse(repetitions, "{2}")
parse(repetitions, "{5,}")

# ## [`Either`](@ref)
# Regular expressions can repeat character matchers and other sub-patterns when appending the `repetitions` suffix.
repeatable = Either{CombinedParser}(
    map(CharIn,char_matcher),
    '.' => AnyChar()
)

# We will add capture groups and other sub-patterns to `repeatable` later.

# ## Repeatable patterns, [`Optional`](@ref), [`Lazy`](@ref) and [`Atomic`](@ref)
# The CombinedParser to parse a repeated pattern from PCRE syntax
@with_names quantified=Sequence(
        repeatable,
        Optional(repetitions, default=(1,1)),
        Optional(map(v->convert(Char,v),with_name("possessive_quantifier",CharIn('+','?'))))
    ) do v
        pat = v[1]
        result = if v[2]==(1,1)
            pat
        elseif v[2]==(0,1)
            Optional(pat)
        else
            Repeat(v[2],pat)
        end
        if v[3] === missing
            result
        elseif v[3]=='+'
            Atomic(result)
        elseif v[3]=='?'
            Lazy(result)
        else
            result
        end
    end

# PCRE repetitions syntax can be `quantified`:
p = quantified(raw"\++?", log=true)
parse(p, "+++")

# ## [`Sequence`](@ref)s and Alternations
# Regular expression patterns can be written in sequence, delimited by `|` for matching either one of several alternatives.
subpattern = Either{CombinedParser}(
    '^' => AtStart(),
    '$' => AtEnd(),
    quantified);

@with_names sequence =
    map( v->sSequence(v...),
         Repeat( subpattern ));
p = sequence(raw"\++?$")
parse(p, "+++")

# The alternation parser transforms a regex string into a parser
@with_names pattern = join(sequence, '|') do v
    sEither(v...)::CombinedParser
end;

result_type(pattern)
p = parse(pattern, "ab|c", log=(:sequence, :pattern))
parse(p, "c")

# ## [`Capture`](@ref)s
# Parentheses allow groupings that are repeatable.
import CombinedParsers.Regexp: name

# The [`push!`](@ref) to [`Either`](@ref) technique allows for the construction of recursive parsers.
@with_names capture_group =
    map(Capture, Sequence(2, '(', pattern, ')'));
push!(repeatable,capture_group);

# Conveniently, the [`@syntax`](@ref) with a `for <name> in <either>` calls `push!`.
@syntax for noncapture_group in repeatable
    Sequence(2, "(?:", pattern, ')')
end;

import CombinedParsers.Regexp: ParserWithCaptures
regcomb = map(ParserWithCaptures,pattern);

match(regcomb("(1a(2b?)*)*0"),"1a1a21a2b22b0")
match(CombinedParsers.Regexp.Regcomb("(1a(2b?)*)*0"),"1a1a21a2b22b0")
match(Regex("(1a(2b?)*)*0"),"1a1a21a2b22b0")
pattern("(1a(2b?)*)*0")("1a1a21a2b22b0")

# ## [`Lookahead`](@ref) and [`Lookbehind`](@ref)

@syntax for lookaround in repeatable
    Sequence(
        "(?",
        :direction => Either(
            '<' => Lookbehind,
            Always() => Lookahead),
        :doesmatch => Either(
            '=' => true,
            '!' => false),
        :pattern => pattern,
        ')') do v
            v.direction(v.doesmatch,v.pattern)::CombinedParser
            ## annotation needed for type inference
            ## obfuscated by :direction
        end
end;

match(pattern("a(?=c|d)."),"abad")
match(Regex("a(?=c|d)."),"abad")

# ## Regular Expression Brackets
# The regular expression parser packaged with CombinedParser is widely compliant with the PCRE C library used in Julia `Base.Regex` -- and in more than half the tests actually faster than C
# [see testing and benchmarking against PCRE tests](pcre-compliance.md).
import CombinedParsers.Regexp: bracket, ParserWithCaptures
push!(repeatable,bracket);

# [`@re_str`](@ref) is tested and benchmarked against the PCRE C library testset (see [compliance report](pcre-compliance.md)).

# Syntax is essential, semantics is high-brow.
# What is the semantics of a regular expression?
# The semantics of regular expression is a parser program that reads a context free regular language.
#


text:
"""
`module  JAC.ElectronCapture`
    ... a submodel of JAC that contains all methods for computing electron-capture properties between some initial and final-state
        multiplets.
"""
module ElectronCapture

    using  Printf, ..AngularMomentum, ..Basics, ..Continuum, ..Defaults, ..InteractionStrength, ..ManyElectron, ..Nuclear,
                   ..Radial, ..TableStrings

    """
    `struct  ElectronCapture.Settings`  ... defines a type for the details and parameters of computing electron-capture lines.

        + printBefore             ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + lineSelection           ::LineSelection                ... Specifies the selected levels, if any.
        + minCaptureEnergy        ::Float64                      ... Minimum energy of free (Auger) electrons to be captured.
        + maxCaptureEnergy        ::Float64                      ... Maximum energy of free (Auger) electrons to be captured.
        + maxKappa                ::Int64                        ... Maximum kappa value of partial waves to be included.
        + operator                ::AbstractEeInteraction        ... Auger/capture operator that is to be used for evaluating the capture amplitudes:
                                                                     allowed values are: CoulombInteraction(), BreitInteraction(), ...
    """
    struct Settings
        printBefore               ::Bool
        lineSelection             ::LineSelection
        minCaptureEnergy          ::Float64
        maxCaptureEnergy          ::Float64
        maxKappa                  ::Int64
        operator                  ::AbstractEeInteraction
    end

    """
    `ElectronCapture.Settings()`  ... constructor for the default values of ElectronCapture line computations
    """
    function Settings()
        Settings(false, LineSelection(), 0., 1.0e5, 2, CoulombInteraction())
    end

    """
    `ElectronCapture.Settings(set::ElectronCapture.Settings;`

            printBefore=..,       minCaptureEnergy=..,       maxCaptureEnergy=..,       maxKappa=..,       operator=..)

        ... constructor for modifying the given ElectronCapture.Settings by 'overwriting' the previously selected parameters.
    """
    function Settings(set::ElectronCapture.Settings;
        printBefore::Union{Nothing,Bool}=nothing,               lineSelection::Union{Nothing,LineSelection}=nothing,
        minCaptureEnergy::Union{Nothing,Float64}=nothing,       maxCaptureEnergy::Union{Nothing,Float64}=nothing,
        maxKappa::Union{Nothing,Int64}=nothing,                 operator::Union{Nothing,String}=nothing)

        if  printBefore      == nothing   printBeforex      = set.printBefore       else  printBeforex      = printBefore        end
        if  lineSelection    == nothing   lineSelectionx    = set.lineSelection     else  lineSelectionx    = lineSelection      end
        if  minCaptureEnergy == nothing   minCaptureEnergyx = set.minCaptureEnergy  else  minCaptureEnergyx = minCaptureEnergy   end
        if  maxCaptureEnergy == nothing   maxCaptureEnergyx = set.maxCaptureEnergy  else  maxCaptureEnergyx = maxCaptureEnergy   end
        if  maxKappa         == nothing   maxKappax         = set.maxKappa          else  maxKappax         = maxKappa           end
        if  operator         == nothing   operatorx         = set.operator          else  operatorx         = operator           end

        Settings( printBeforex, lineSelectionx, minCaptureEnergyx, maxCaptureEnergyx, maxKappax, operatorx)
    end

    # `Base.show(io::IO, settings::ElectronCapture.Settings)`  ... prepares a proper printout of the variable settings::ElectronCapture.Settings.
    function Base.show(io::IO, settings::ElectronCapture.Settings)
        println(io, "printBefore:                   $(settings.printBefore)  ")
        println(io, "lineSelection:                 $(settings.lineSelection)  ")
        println(io, "minCaptureEnergy:              $(settings.minCaptureEnergy)  ")
        println(io, "maxCaptureEnergy:              $(settings.maxCaptureEnergy)  ")
        println(io, "maxKappa:                      $(settings.maxKappa)  ")
        println(io, "operator:                      $(settings.operator)  ")
    end

    """
    `struct  Channel`
        ... defines a type for a ElectronCapture channel to help characterize a scattering (continuum) state of many
            electron-states with a single free electron.

        + kappa          ::Int64                ... partial-wave of the free electron
        + symmetry       ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + phase          ::Float64              ... phase of the partial wave
        + amplitude      ::Complex{Float64}     ... electron-capture amplitude associated with the given channel.
    """
    struct  Channel
        kappa            ::Int64
        symmetry         ::LevelSymmetry
        phase            ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  Line`
        ... defines a type for a ElectronCapture line that may include the definition of sublines and their
            corresponding amplitudes.

        + initialLevel   ::Level           ... initial-(state) level
        + finalLevel     ::Level           ... final-(state) level
        + electronEnergy ::Float64         ... Energy of the (incoming free) electron.
        + totalRate      ::Float64         ... Total rate of this line.
        + hasChannels    ::Bool
            ... Determines whether the individual scattering (sub-) channels are defined in terms of their free-electron energy, kappa
                and the total angular momentum/parity as well as the amplitude, or not.
        + channels       ::Array{ElectronCapture.Channel,1}  ... List of ElectronCapture channels of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        electronEnergy   ::Float64
        totalRate        ::Float64
        hasChannels      ::Bool
        channels         ::Array{ElectronCapture.Channel,1}
    end

    """
    `ElectronCapture.Line(initialLevel::Level, finalLevel::Level, totalRate::Float64)`
        ... constructor for an ElectronCapture line between a specified initial and final level.
    """
    function Line(initialLevel::Level, finalLevel::Level, totalRate::Float64)
        Line(initialLevel, finalLevel, 0., totalRate, false, ElectronCapture.Channel[])
    end

    # `Base.show(io::IO, line::ElectronCapture.Line)`  ... prepares a proper printout of the variable line::ElectronCapture.Line.
    function Base.show(io::IO, line::ElectronCapture.Line)
        println(io, "initialLevel:           $(line.initialLevel)  ")
        println(io, "finalLevel:             $(line.finalLevel)  ")
        println(io, "electronEnergy:         $(line.electronEnergy)  ")
        println(io, "totalRate:              $(line.totalRate)  ")
        println(io, "hasChannels:            $(line.hasChannels)  ")
        println(io, "channels:               $(line.channels)  ")
    end

    """
    `ElectronCapture.amplitude(kind::String, channel::ElectronCapture.Channel, finalLevel::Level, continuumLevel::Level,
                               grid::Radial.Grid; printout::Bool=true)`
        ... to compute the kind in  CoulombInteraction(), BreitInteraction(), CoulombBreit()   ElectronCapture amplitude
            <alpha_f J_f || O^(capture, kind) || (alpha_i J_i, kappa) J_f> = <(alpha_i J_i, kappa) J_f || O^(Auger, kind) || alpha_f J_f>^*
            due to the interelectronic interaction for the given final and initial level. A value::ComplexF64 is returned.
    """
    function amplitude(kind::String, channel::ElectronCapture.Channel, finalLevel::Level, continuumLevel::Level, grid::Radial.Grid; printout::Bool=true)
        aChannel  = AutoIonization.Channel(channel.kappa, channel.symmetry, channel.phase, 0.)
        amplitude = conj(AutoIonization.amplitude(kind, aChannel, continuumLevel::Level, finalLevel::Level, grid; printout=printout))

        return( amplitude )
    end

    """
    `ElectronCapture.computeAmplitudesProperties(line::ElectronCapture.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                                 settings::ElectronCapture.Settings; printout::Bool=true)`
        ... to compute all amplitudes and properties of the given line; a line::ElectronCapture.Line is returned for which the amplitudes
            and properties are now evaluated.
    """
    function computeAmplitudesProperties(line::ElectronCapture.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                         settings::ElectronCapture.Settings; printout::Bool=true)
        newChannels = ElectronCapture.Channel[];   contSettings = Continuum.Settings(false, nrContinuum);   rate = 0.
        for channel in line.channels
            newfLevel = Basics.generateLevelWithSymmetryReducedBasis(line.finalLevel, line.finalLevel.basis.subshells)
            newiLevel = Basics.generateLevelWithSymmetryReducedBasis(line.initialLevel, newfLevel.basis.subshells)
            newfLevel = Basics.generateLevelWithExtraSubshell(Subshell(101, channel.kappa), newfLevel)
            cOrbital, phase  = Continuum.generateOrbitalForLevel(line.electronEnergy, Subshell(101, channel.kappa), newiLevel, nm, grid, contSettings)
            newcLevel  = Basics.generateLevelWithExtraElectron(cOrbital, channel.symmetry, newiLevel)
            newChannel = ElectronCapture.Channel(channel.kappa, channel.symmetry, phase, 0.)
            amplitude  = ElectronCapture.amplitude(settings.operator, newChannel, newcLevel, newfLevel, grid, printout=printout)
            rate       = rate + conj(amplitude) * amplitude
            push!( newChannels, ElectronCapture.Channel(newChannel.kappa, newChannel.symmetry, newChannel.phase, amplitude) )
        end
        totalRate = 2pi* rate
        newLine   = ElectronCapture.Line(line.initialLevel, line.finalLevel, line.electronEnergy, totalRate, true, newChannels)

        return( newLine )
    end

    """
    `ElectronCapture.displayRates(stream::IO, lines::Array{ElectronCapture.Line,1}, settings::ElectronCapture.Settings)`
        ... to list all results, energies, rates, etc. of the selected lines. A neat table is printed but nothing is returned
            otherwise.
    """
    function  displayRates(stream::IO, lines::Array{ElectronCapture.Line,1}, settings::ElectronCapture.Settings)
        nx = 96
        println(stream, " ")
        println(stream, "  Auger rates: \n")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(12, "Energy"   ; na=2);
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(14, "Electron energy"   ; na=2);
        sb = sb * TableStrings.center(14,TableStrings.inUnits("energy"); na=2)
        sa = sa * TableStrings.center(14, "Capture rate"; na=2);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("rate"); na=2)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.initialLevel.energy))  * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.electronEnergy))       * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("rate: from atomic", line.totalRate))              * "    "
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module


text:
# Helper functions for nMDS based on # https://link.springer.com/content/pdf/10.1007/BF02289694.pdf
# Varibales are generally named simillaryly to thos in the original paper. However
# it felt simples to have a a vector of Is and one of Js, rather than an IJ array.
# The complete algorithm, using these functions can be found in main.jl

#### Gradient calculations (for the Euclidean case) pp.126 of Kruskal(1964)

# There are 3 terms within the sum on the right hand side of the equation
# The first term simply denotes which values affect which elements of X2 and
# whether to add or subtract them. The second and third terms can be calculated
# separately to form a length M vector and an M by L arrau respectively. Then we
# loop over L and M, using the indeces in Is and Js to add t appropriate values to X2.
# We multiply by minus S, since we are interested in the negative gradient.

# Second term
function grhs2!(d, d̂, S_star,T_star,M,res)
    @turbo for m in 1:M
        res[m] = (d[m]-d̂[m])/S_star - d[m]/T_star
    end
end
# Third term
function grhs3!(x,d, M,L,res,Is,Js)
    Imat = view(x,Is,:)
    Jmat = view(x,Js,:)
    @turbo for m in 1:M
        for l in 1:L
            res[m,l] = (Imat[m,l] - Jmat[m,l])/d[m]
        end
    end
end
# put them together
function g(grhs2,grhs3,L,M,S,Is,Js,X2,grhs)
    @inbounds for l in 1:L
       grhs .= grhs2 .* view(grhs3,:,l)
       for m in 1:M
           X2[Is[m],l] += grhs[m]
           X2[Js[m],l] -= grhs[m]
       end
   end
   return X2 *-S
end

function checkup(start,block,partition,d,M)
    println("up")

    ∑ = 0
    i = start
    while i < M
    #while partition[i] == block
        println(i,"","hi")
        if partition[i] !== block
            break
        end
        ∑+=d[i]
        i+=1
    end
    return ∑/(i -start), i-1
end
function checkdown(start,block,partition,d)
    println("down")
    ∑ = 0
    i = start
    while  i>0
        println(i,"","hip")
        if partition[i] !== block
            break
        end
        ∑+=d[i]
        i-=1

    end
    return ∑/(start-i), i+1
end

#### Monotone rgression pp.126-128 Kruscal(1964)
#The R implementation monMDS also explains this procedure https://github.com/vegandevs/vegan/blob/master/src/monoMDS.f

# The aim is to sort the distances (in DIST) into blocks, such that the mean (d̂) of each block increases
# from first to last i.e. for any given block - the block before it will have a lower d̂ and the one after it will
# have a greater d̂. In Kruskals terms, each block must be downsatisfied (block before it has lower d̂) and up
# satisfied (block after it has higher d̂).
# To achieve this we start with each element of d having it's own singleton block.
# Starting with the first block, we alternate between checking whether the block is upsatisfied and downsatisfied, whenever the block
# is not satisfied for a given direction we merge it with the adjacent block in that direction and begin
# checking that this block is up and down satified. If the current block is both up and down satisfied,
# we move on to the next block and do the same. This procedure continues until all blocks are up and down
# satisfied.
# The fitted values returned to DHAT are the d̂ values associated with the block each DIST value belongs to.

function monotone!(d,d̂,M, partition = false)
    loc = 1
    # if a starting partition is supplied e.g. because of ties, use it
    partition = partition == false ? collect(1:M) : partition
    block = view(d,partition .==loc)
    localmean = d[1]
    bounds = (1,1)
    #start at position 1, up-active i.e. check to see that the block to the right has a higher d̂
    while true#(loc < maximum(partition)) | (localmean < checkdown(bounds[1],loc-1,partition,d)[1])
    println(loc)
    println(bounds)
    println(partition[1:10])
    if loc < maximum(partition)
        upmean, upbound = checkup(bounds[2]+1,loc+1,partition,d,M)
        println(localmean,"",upmean)
        if upmean > localmean # if up satissfied
            if loc > 1 # and not first block
                downmean, downbound = checkdown(bounds[1]-1,loc-1,partition,d)#check down
                if downmean < localmean # if down satisfied also
                    loc += 1 # move to next block
                    localmean = upmean
                    bounds = (bounds[2]+1,upbound)

                else # if not down satisfied
                    partition[bounds[1]:end] .-=1
                    localmean = (localmean*(bounds[2] -bounds[1]+1)+downmean*(bounds[1]-downbound))/(bounds[2]-downbound+1)
                    bounds = (downbound,bounds[2]) # merge down
                    loc -= 1
                end
            else #if first block and up satisfied
                loc += 1 # move to next block
                localmean = upmean
                bounds = (bounds[2]+1,upbound)

            end
        else # if not upsatisfied
            partition[bounds[2]+1:end] .-=1
            localmean = (localmean*(bounds[2] -bounds[1]+1)+upmean*(upbound-bounds[2]))/(upbound-bounds[1]+1)
            bounds = (bounds[1],upbound) # merge up
            #loc += 1
            if loc > 1 # and not first block
                downmean, downbound = checkdown(bounds[1]-1,loc-1,partition,d)#check down
              #  if downmean < localmean # if down satisfied also
                 #   loc += 1 # move to next block
                if   downmean > localmean
               # else # if not down satisfied
                    partition[bounds[1]:end] .-=1
                    localmean = (localmean*(bounds[2] -bounds[1]+1)+downmean*(bounds[1]-downbound))/(bounds[2]-downbound+1)
                    bounds = (downbound,bounds[2]) # merge down
                    loc -= 1
                 end
            end
        end
    else # if last block
        println("last block")
        downmean, downbound = checkdown(bounds[1]-1,loc-1,partition,d)#check down
        if downmean < localmean # if down satisfied
            break # exit loop
        else # if not down satisfied
            partition[bounds[1]:end] .-=1
            localmean = (localmean*(bounds[2] -bounds[1]+1)+downmean*(bounds[1]-downbound))/(bounds[2]-downbound+1)
            bounds = (downbound,bounds[2]) # merge down
            loc -= 1
        end
    end

    end

    blockmax = partition[end]
    blocks = 1:blockmax

    for block in blocks
        boolmask = partition .== block
        d̂[boolmask] .= mean(d[boolmask])
    end
end

# normalize point coordinates
function normalize(x,dim)
    x .-= mean(x,dims = dim)
    return x ./std(x,dims = dim)
end

#### Step size related function pp.120-122 of Kruskal(1964)
cosθ(g,g″) = sum(g .* g″)/(sqrt(sum(g .^2))*sqrt(sum(g″ .^2)))
angle_factor(g,g″) = 4.0^(cosθ(g,g″)^3.0)
relaxation_factor(x) = 1.3/(1+(five_step_ratio(x)^5))
five_step_ratio(x) = min(1,x[end]/x[1])
good_luck_factor(x) = min(1,x[end]/x[end-1])
mag(g,K) = sqrt(sum(g .^2)/K)
x′(x,g,α,K) = x .+ g*α/mag(g,K)

# Calculate the stress (S), S⋆ an T⋆ according to Kruskals stress formula pp.115 & 125 Kruskal(1964)
function stress(d,d̂,M)
    S_star = 0
    T_star = 0
    @turbo for m in 1:M
        S_star += (d[m] - d̂[m])^2
        T_star += d[m]^2
    end
    return S_star, T_star, sqrt(S_star/T_star)
end
# calculate Euclidean distances between points
function get_dist!(X1,DIST,Is,Js,M)
    @inbounds for m in 1:M
        DIST[m] = euclidean(X1[Is[m],:],X1[Js[m],:])
    end
end


text:
export project_l1_Duchi!
function project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF) where {TF<:Real}
# % PROJECTONTOL1BALL Projects point onto L1 ball of specified radius.
# %
# % w = ProjectOntoL1Ball(v, b) returns the vector w which is the solution
# %   to the following constrained minimization problem:
# %
# %    min   ||w - v||_2
# %    s.t.  ||w||_1 <= b.
# %
# %   That is, performs Euclidean projection of v to the 1-norm ball of radius
# %   b.
# %
# % Author: John Duchi (jduchi@cs.berkeley.edu)
# Translated (with some modification) to Julia 1.1 by Bas Peters
if (b <= TF(0))
  error("Radius of L1 ball is negative");
end
if norm(v, 1) <= b
  return v
end
lv=length(v)
u = Vector{TF}(undef,lv)
sv= Vector{TF}(undef,lv)

#u = sort(abs(v),'descend');
#u = sort(abs.(v), rev=true) #faster than the line below
#u = sort(v, by=abs , rev=true)

#use RadixSort for Float32 (short keywords)
u=copy(v)
if TF==Float32
 u = sort!(abs.(u), rev=true,alg=RadixSort)
else
 u = sort!(abs.(u), rev=true,alg=QuickSort)
end

#sv = cumsum(u);
sv = cumsum(u)

#rho = find(u > (sv - b) ./ (1:length(u))', 1, 'last');
#tmp = (u .> ((sv.-b)./ LinSpace(1,lv,lv)))#::BitVector

rho = max(1,min(lv,findlast(u .> ((sv.-b)./ (1.0:1.0:lv)))))
#convert(TF,rho) why was this here...
#rho = findlast(tmp)::Int64

#theta = max(0, (sv(rho) - b) / rho);
theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF

#w = sign(v) .* max(abs(v) - theta, 0);
v .= sign.(v) .* max.(abs.(v) .- theta, TF(0))

return v
end #end project_l1_Duchi
#
# function project_l1_Duchi!{TF<:Real}(v::Union{Array{TF,3},Array{Complex{TF},3}}, b::TF,mode::String)
#
# if (b < TF(0))
#   error("Radius of L1 ball is negative");
# end
# (n1,n2,n3)=size(v)
#
# ##Slice based projection
# if mode == "x_slice"
#   u = Vector{TF}(n2*n3)
#   sv= Vector{TF}(n2*n3)
#   v = reshape(v,n1,n2*n3)
#   Threads.@threads for i=1:n1
#     u = sort(abs.(v[i,:]), rev=true)
#     sv = cumsum(u)
#     rho = findlast(u .> ((sv.-b)./ (1.0:1.0:(n2*n3))))
#     theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF
#     @inbounds v[i,:] .= sign.(v[i,:]) .* max.(abs.(v[:,i]).-theta, TF(0))
#   end
#   v = reshape(v,n1,n2,n3)
# elseif mode == "z_slice"
#   u = Vector{TF}(n1*n2)
#   sv= Vector{TF}(n1*n2)
#   v = reshape(v,n1*n2,n3)
#   Threads.@threads for i=1:n3
#     u = sort(abs.(v[:,i]), rev=true)
#     sv = cumsum(u)
#     rho = findlast(u .> ((sv.-b)./ (1.0:1.0:(n1*n2))))
#     theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF
#     @inbounds v[:,i] .= sign.(v[:,i]) .* max.(abs.(v[:,i]).-theta, TF(0))
#   end
#   v = reshape(v,n1,n2,n3)
# elseif mode == "y_slice"
#   permutedims(v,[2;1;3]);
#   v = reshape(v,n2,n1*n3)
#   u = Vector{TF}(n1*n3)
#   sv= Vector{TF}(n1*n3)
#   Threads.@threads for i=1:size(v,1)
#     u = sort(abs.(v[i,:]), rev=true)
#     sv = cumsum(u)
#     rho = findlast(u .> ((sv.-b)./ (1.0:1.0:(n1*n3))))
#     theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF
#     @inbounds v[i,:] .= sign.(v[i,:]) .* max.(abs.(v[:,i]).-theta, TF(0))
#   end
#   v = reshape(v,n2,n1,n3);
#   permutedims(v,[2;1;3]);
# else
#   error("provided incorred mode for l1_projection")
# end #end if
# v=vec(v)
#
# return v
# end #end project_l1_Duchi


text:
using FFTW

"""
  ksbenchmark(Nx, ksintegrator,Nruns=4,printnroms=false): benchmark a kuramoto-sivashinksy integration algorithm
    on Nx gridpoints. Usage example: ksbenchmark(512, ksintegrateNaive)
"""
function ksbenchmark(Nx, ksintegrator, Nruns=4, printnorms=false)

    Lx = Nx/16*pi             # spatial domain [0, L] periodic
    dt = 1/16                 # discrete time step
    T  = 200                  # integrate from t=0 to t=T
    nplot = round(Int,1/dt)   # save every nploth time step
    Nt = round(Int, T/dt)     # total number of timesteps

    x = Lx*(0:Nx-1)/Nx
    u0 = cos.(x) + 0.1*sin.(x/8) + 0.01*cos.((2*pi/Lx)*x) .+ 0im
    u = copy(u0)

    skip = 1
    avgtime = 0

    for run=1:Nruns
        cputime = @elapsed ksintegrator(u0, Lx, dt, Nt)
        if run > skip
            avgtime = avgtime + cputime
        end
    end

    if printnorms
        @show ksnorm(u0)
        @show ksnorm(u)
    end

    avgtime = avgtime/(Nruns-skip)
    @show avgtime

end

function ksnorm(u)
    s = 0.0
    for n = 1:length(u)
        s += (u[n])^2
    end
    sqrt(s/length(u))
end

"""
ksintegrateNaive: integrate kuramoto-sivashinsky equation (Julia)
       u_t = -u*u_x - u_xx - u_xxxx, domain x in [0,Lx], periodic BCs

 inputs
          u = initial condition (vector of u(x) values on uniform gridpoints))
         Lx = domain length
         dt = time step
         Nt = number of integration timesteps
      nsave = save every nsave-th time step

 outputs
          u = final state, vector of u(x, Nt*dt) at uniform x gridpoints

This a line-by-line translation of a Matlab code into Julia. It uses out-of-place
FFTs and doesn't pay any attention to the allocation of temporary vectors within
the time-stepping loop. Hence the name "naive".
"""
function ksintegrateNaive(u, Lx, dt, Nt)
    Nx = length(u)                  # number of gridpoints
    x = collect(0:(Nx-1)/Nx)*Lx
    kx = vcat(0:Nx/2-1, 0, -Nx/2+1:-1)  # integer wavenumbers: exp(2*pi*kx*x/L)
    alpha = 2*pi*kx/Lx              # real wavenumbers:    exp(alpha*x)
    D = 1im*alpha                  # D = d/dx operator in Fourier space
    L = alpha.^2 - alpha.^4         # linear operator -D^2 - D^4 in Fourier space
    G = -0.5*D                      # -1/2 D operator in Fourier space

    # Express PDE as u_t = Lu + N(u), L is linear part, N nonlinear part.
    # Then Crank-Nicolson Adams-Bashforth discretization is
    #
    # (I - dt/2 L) u^{n+1} = (I + dt/2 L) u^n + 3dt/2 N^n - dt/2 N^{n-1}
    #
    # let A = (I - dt/2 L)
    #     B = (I + dt/2 L), then the CNAB timestep formula
    #
    # u^{n+1} = A^{-1} (B u^n + 3dt/2 N^n - dt/2 N^{n-1})

    # convenience variables
    dt2  = dt/2
    dt32 = 3*dt/2
    A =  ones(Nx) + dt2*L
    B = (ones(Nx) - dt2*L).^(-1)

    Nn  = G.*fft(u.*u) # -u u_x (spectral), notation Nn = N^n     = N(u(n dt))
    Nn1 = copy(Nn)     #                   notation Nn1 = N^{n-1} = N(u((n-1) dt))
    u  = fft(u)        # transform u to spectral

    # timestepping loop
    for n = 1:Nt
        Nn1 = copy(Nn)                 # shift nonlinear term in time: N^{n-1} <- N^n
        Nn  = G.*fft(real(ifft(u)).^2) # compute Nn = -u u_x

        u = B .* (A .* u + dt32*Nn - dt2*Nn1)
    end

    real(ifft(u))
end

"""
ksintegrate: integrate kuramoto-sivashinsky equation (Julia)
       u_t = -u*u_x - u_xx - u_xxxx, domain x in [0,Lx], periodic BCs

 inputs
          u = initial condition (vector of u(x) values on uniform gridpoints))
         Lx = domain length
         dt = time step
         Nt = number of integration timesteps
      nsave = save every nsave-th time step

 outputs

          u = final state, vector of u(x, Nt*dt) at uniform x gridpoints

This implementation has two improvements over ksintegrateNaive.jl. It uses
   (1) in-place FFTs
   (2) loop fusion: Julia can translate arithmetic vector expressions in dot syntax
       to single for loop over the components, which should be much faster than
       constructing a temporary vector for each operation in the vector expression.
"""
function ksintegrateInplace(u, Lx, dt, Nt)
    u = (1+0im)*u                       # force u to be complex
    Nx = length(u)                      # number of gridpoints
    kx = vcat(0:Nx/2-1, 0:0, -Nx/2+1:-1)# integer wave #s,  exp(2*pi*i*kx*x/L)
    alpha = 2*pi*kx/Lx                  # real wavenumbers, exp(i*alpha*x)
    D = 1im*alpha                       # spectral D = d/dx operator
    L = alpha.^2 - alpha.^4             # spectral L = -D^2 - D^4 operator
    G = -0.5*D                          # spectral -1/2 D operator, to eval -u u_x = 1/2 d/dx u^2

    # convenience variables
    dt2  = dt/2
    dt32 = 3*dt/2
    A_inv = (ones(Nx) - dt2*L).^(-1)
    B     =  ones(Nx) + dt2*L

    # compute FFTW plans
    FFT! = plan_fft!(u, flags=FFTW.ESTIMATE)
    IFFT! = plan_ifft!(u, flags=FFTW.ESTIMATE)

    # compute nonlinear term Nu == -u u_x and Nuprev (Nu at prev timestep)
    Nu     = G.*fft(u.^2) # Nu == -1/2 d/dx (u^2) = -u u_x
    Nuprev = copy(Nu)      # use Nuprev = Nu at first time step
    FFT!*u

    # timestepping loop
    for n = 1:Nt

        Nuprev .= Nu   # shift nonlinear term in time
        Nu .= u         # put u into N in prep for comp of nonlineat

        IFFT!*Nu       # transform Nu to gridpt values, in place
        Nu .= Nu.*Nu   # collocation calculation of u^2
        FFT!*Nu        # transform Nu back to spectral coeffs, in place

        Nu .= G.*Nu

        # loop fusion! Julia translates the folling line of code to a single for loop.
        u .= A_inv .*(B .* u .+ dt32.*Nu .- dt2.*Nuprev)
    end

    IFFT!*u
    real(u)
end

"""
ksintegrateUnrolled: integrate kuramoto-sivashinsky equation (Julia)

      u_t = -u*u_x - u_xx - u_xxxx, domain x in [0,Lx], periodic BCs

 inputs
          u = initial condition (vector of u(x) values on uniform gridpoints))
         Lx = domain length
         dt = time step
         Nt = number of integration timesteps
      nsave = save every nsave-th time step

 outputs
          u = final state, vector of u(x, Nt*dt) at uniform x gridpoints

This implementation has one improvements over ksintegrateInplace.jl. It explicitly
unrolls all the vector equations into for loops. For some reason this works
better (runs faster) than the unrolled loops produced by the Julia compiler.
"""
function ksintegrateUnrolled(u, Lx, dt, Nt)
    u = (1+0im)*u                       # force u to be complex
    Nx = length(u)                      # number of gridpoints
    kx = vcat(0:Nx/2-1, 0:0, -Nx/2+1:-1)# integer wavenumbers: exp(2*pi*kx*x/L)
    alpha = 2*pi*kx/Lx                  # real wavenumbers:    exp(alpha*x)
    D = 1im*alpha                       # spectral D = d/dx operator
    L = alpha.^2 - alpha.^4             # spectral L = -D^2 - D^4 operator
    G = -0.5*D                          # spectral -1/2 D operator, to eval -u u_x = 1/2 d/dx u^2

    # convenience variables
    dt2  = dt/2
    dt32 = 3*dt/2
    A =  ones(Nx) + dt2*L
    B = (ones(Nx) - dt2*L).^(-1)

    # compute FFTW plans
    FFT! = plan_fft!(u, flags=FFTW.ESTIMATE)
    IFFT! = plan_ifft!(u, flags=FFTW.ESTIMATE)

    # compute nonlinear term Nu == -u u_x and Nuprev (Nu at prev timestep)
    Nn  = G.*fft(u.^2) # Nnf == -1/2 d/dx (u^2) = -u u_x, spectral
    Nn1 = copy(Nn)     # use Nnf1 = Nnf at first time step
    FFT!*u

    # timestepping loop
    for n = 1:Nt

        Nn1 .= Nn
        Nn .= u

        IFFT!*Nn # in-place FFT

        @inbounds for i = 1:length(Nn)
            @fastmath Nn[i] = Nn[i]*Nn[i]
        end

        FFT!*Nn

        @inbounds for i = 1:length(Nn)
            @fastmath Nn[i] = G[i]*Nn[i]
        end

        @inbounds for i = 1:length(u)
            @fastmath u[i] = B[i]* (A[i] * u[i] + dt32*Nn[i] - dt2*Nn1[i])
        end

    end

    IFFT!*u
    real(u)
end

"""
   Construct initial conditions for benchmarking ksintegrate* algorithms
   Useful for using Julia's benchmark utilities.
"""
function ksinitconds(Nx)
    Lx = Nx/16*pi             # spatial domain [0, L] periodic
    dt = 1/16                 # discrete time step
    T  = 200                  # integrate from t=0 to t=T
    nplot = round(Int,1/dt)   # save every nploth time step
    Nt = round(Int, T/dt)     # total number of timesteps

    x = Lx*(0:Nx-1)/Nx
    u0 = cos.(x) + 0.1*sin.(x/8) + 0.01*cos.(x/16)

    u0, Lx, dt, Nt
end

max_stars_repo_name: KZiemian/julia-pde-benchmark

max_issues_repo_name: KZiemian/julia-pde-benchmark

max_forks_repo_name: KZiemian/julia-pde-benchmark


text:
"""
    belief_ellipse(b::GaussianBelief, P::Float=0.95; δ::Number=5)

Construct and return the x and y points of a 2D gaussian belief,
with P being the total probability captured by the ellipse (P ∈ (0,1)),
and δ the degree increment between points.
"""
function belief_ellipse(b::GaussianBelief, P::Number=0.95; δ::Number=5)
    @assert 0<P<1
    @assert 0<δ<360

    θ = (pi/180)*collect(0:δ:360)
    rad_w = (-2*log(1-P))^0.5

    w1 = rad_w*cos.(θ)
    w2 = rad_w*sin.(θ)

    Sig12 = cholesky(b.Σ)
    x=Sig12.L*[w1 w2]' .+ b.μ
    return x[1,:], x[2,:]
end

"""
    belief_ellipse(μ::AbstractVector, Σ::AbstractMatrix, P::Float=0.95; δ::Number=5)

Construct and return the x and y points of a 2D gaussian belief,
with P being the total probability captured by the ellipse (P ∈ (0,1)),
and δ the degree increment between points.
"""

function belief_ellipse(μ::AbstractVector, Σ::AbstractMatrix, P::Number=0.95; δ::Number=5)
    @assert 0<P<1
    @assert 0<δ<360

    θ = (pi/180)*collect(0:δ:360)
    rad_w = (-2*log(1-P))^0.5

    w1 = rad_w*cos.(θ)
    w2 = rad_w*sin.(θ)

    Sig12 = cholesky(Σ)
    x = Sig12.L*[w1 w2]' .+ μ
    return x[1,:], x[2,:]
end


text:
using Flux: @functor
using MacroTools: @forward

using ..Basic
using ..Basic: AbstractTransformer
using ..Stacks

struct Bert{T<:Stack, DP<:Dropout} <: AbstractTransformer
  ts::T
  drop::DP
end

@functor Bert

@forward Bert.ts Base.getindex, Base.length

"""
    Bert(size::Int, head::Int, ps::Int, layer::Int;
        act = gelu, pdrop = 0.1, attn_pdrop = 0.1)
    Bert(size::Int, head::Int, hs::Int, ps::Int, layer::Int;
        act = gelu, pdrop = 0.1, attn_pdrop = 0.1)

the Bidirectional Encoder Representations from Transformer(BERT) model.

    (bert::Bert)(x, mask=nothing; all::Bool=false)

eval the bert layer on input `x`. If length `mask` is given (in shape (1, seq_len, batch_size)), mask the attention with `getmask(mask, mask)`. Moreover, set `all` to `true` to get all
outputs of each transformer layer.
"""
function Bert(size::Int, head::Int, ps::Int, layer::Int;
              act = gelu, pdrop = 0.1, attn_pdrop = 0.1)
  rem(size,  head) != 0 && error("size not divisible by head")
  Bert(size, head, div(size, head), ps, layer; act=act, pdrop=pdrop, attn_pdrop=attn_pdrop)
end

function Bert(size::Int, head::Int, hs::Int, ps::Int, layer::Int; act = gelu, pdrop = 0.1, attn_pdrop = 0.1)
  Bert(
    Stack(
      @nntopo_str("((x, m) => x':(x, m)) => $layer"),
      [
        Transformer(size, head, hs, ps; future=true, act=act, pdrop=attn_pdrop)
        for i = 1:layer
      ]...
    ),
    Dropout(pdrop))
end

function (bert::Bert)(x::T, mask=nothing; all::Bool=false) where T
  e = bert.drop(x)

  if mask === nothing
    t, ts = bert.ts(e, nothing)
  else
    t, ts = bert.ts(e, getmask(mask, mask))
  end

  if all
    if mask !== nothing
      ts = map(ts) do ti
        ti .* mask
      end
    end
    ts[end], ts
  else
    t = mask === nothing ? t : t .* mask
    t
  end
end

"""
    masklmloss(embed::Embed{T}, transform,
               t::AbstractArray{T, N}, posis::AbstractArray{Tuple{Int,Int}}, labels) where {T,N}
    masklmloss(embed::Embed{T}, transform, output_bias,
               t::AbstractArray{T, N}, posis::AbstractArray{Tuple{Int,Int}}, labels) where {T,N}

helper function for computing the maks language modeling loss.
Performance `transform(x) .+ output_bias` where `x` is the mask specified by
`posis`, then compute the similarity with `embed.embedding` and crossentropy between true `labels`.
"""
function masklmloss(embed::Embed{T}, transform, t::AbstractArray{T, N}, posis::AbstractArray{Tuple{Int,Int}}, labels) where {T,N}
  masktok = gather(t, posis)
  sim = logsoftmax(transpose(embed.embedding) * transform(masktok))
  return logcrossentropy(labels, sim)
end

function masklmloss(embed::Embed{T}, transform, output_bias, t::AbstractArray{T, N}, posis::AbstractArray{Tuple{Int,Int}}, labels) where {T,N}
  masktok = gather(t, posis)
  sim = logsoftmax(transpose(embed.embedding) * transform(masktok) .+ output_bias)
  return logcrossentropy(labels, sim)
end

function Base.show(io::IO, bert::Bert)
  hs = div(size(bert.ts[1].mh.iqproj.W)[1], bert.ts[1].mh.head)
  h, ps = size(bert.ts[1].pw.dout.W)

  print(io, "Bert(")
  print(io, "layers=$(length(bert.ts)), ")
  print(io, "head=$(bert.ts[1].mh.head), ")
  print(io, "head_size=$(hs), ")
  print(io, "pwffn_size=$(ps), ")
  print(io, "size=$(h))")
end


text:
#
# Maneuver calculations for orbits within the Twobody Problem
#

"""
Provides the radius of escape.
"""
escape_radius(r₀, v₀, aₜ) = r₀ * v₀ / (20 * aₜ^2 * r₀^2)^(1/4)
escape_radius(orbit::T, m::ConstantManeuver) where T <: RestrictedTwoBodySystem = escape_radius(radius(orbit), velocity(orbit), m.aₜ)

"""
Provides the velocity at escape.
"""
escape_velocity(r₀, v₀, aₜ, μ) = √(2 * μ / escape_radius(r₀, v₀, aₜ))
escape_velocity(orbit::T, m::ConstantManeuver) where T <: RestrictedTwoBodySystem = escape_velocity(radius(orbit), velocity(orbit), m.aₜ, orbit.body.μ)

"""
Provides time delta from the provided initial orbit to escape.
"""
escape_time(r₀, v₀, aₜ) = upreferred(v₀ / aₜ) * (1 - (upreferred(20aₜ^2 * r₀^2) / upreferred(v₀^4))^(1/8))
escape_time(orbit::T, m::ConstantManeuver) where T <: RestrictedTwoBodySystem = escape_time(radius(orbit), velocity(orbit), m.aₜ)

"""
Provides the path length from the initial condition to escape.
"""
escape_path_length(r₀, v₀, aₜ) = upreferred(v₀^2 / 2aₜ) * (1 - upreferred((1/v₀) * (20aₜ^2 * r₀^2)^(1/4)))
escape_path_length(orbit::T, m::ConstantManeuver) where T <: RestrictedTwoBodySystem = escape_path_length(radius(orbit), velocity(orbit), m.aₜ)


text:
#########################################################################
# Spatial functions
#########################################################################
function spatialidxs(::UnstructuredGrid, A)
    return ((Coord(i),) for i in 1:size(A, Coord))
end

function zonalmean(::UnstructuredGrid, A::AbDimArray, ::Nothing)
    idxs, lats = uniquelats(A)
    other = otherdims(A, Coord())
    r = zeros(eltype(A), (length(lats), size.(Ref(A), other)...))
    R = ClimArray(r, (Lat(lats), other...); name=A.name, attrib=A.attrib)
    for (i, r) in enumerate(idxs)
        for j in otheridxs(A, Coord())
            R[Lat(i), j...] = mean(view(A, Coord(r), j...))
        end
    end
    return R
end
function zonalmean(::UnstructuredGrid, A::AbDimArray{T, 1}, ::Nothing) where {T}
    idxs, lats = uniquelats(A)
    res = zeros(T, length(lats))
    for (i, r) in enumerate(idxs)
        res[i] = mean(view(A.data, r))
    end
    return ClimArray(res, (Lat(lats),); name=A.name, attrib=A.attrib)
end

# zonal mean with weights
function zonalmean(::UnstructuredGrid, A::ClimArray, W)
	@assert size(A) == size(W)
    idxs, lats = uniquelats(A)
    other = otherdims(A, Coord())
    r0 = zeros(eltype(A), (length(lats), size.(Ref(A), other)...))
    R = ClimArray(r0, (Lat(lats), other...); name=A.name, attrib=A.attrib)
    for (i, r) in enumerate(idxs)
        for j in otheridxs(A, Coord())
            R[Lat(i), j...] = mean(view(A, Coord(r), j...), weights(view(W, Coord(r), j...)))
        end
    end
    return R
end
function zonalmean(::UnstructuredGrid, A::ClimArray{T, 1}, W) where {T}
	@assert size(A) == size(W)
    idxs, lats = uniquelats(A)
    res = zeros(T, length(lats))
    for (i, r) in enumerate(idxs)
        res[i] = mean(view(A.data, r), weights(view(W.data, r)))
    end
    return ClimArray(res, (Lat(lats),); name=A.name, attrib=A.attrib)
end

"""
    uniquelats(A::AbDimArray) → idxs, lats
    uniquelats(c::Vector{<:AbstractVector}) → idxs, lats
Find the unique latitudes of `A`. Return the indices (vector of ranges) that each latitude
in `lats` covers, as well as the latitudes themselves.
"""
uniquelats(A::AbDimArray) = uniquelats(dims(A, Coord).val)
function uniquelats(c)
    @assert issorted(c; by = x -> x[2])
    idxs = Vector{UnitRange{Int}}()
    lats = eltype(eltype(c))[]
    iprev = 1
    for i in 2:length(c)
        if c[i][2] != c[i-1][2]
            push!(lats, c[i-1][2])
            push!(idxs, iprev:(i-1))
            iprev = i
        end
    end
    push!(lats, c[end][2])
    push!(idxs, iprev:length(c))
    return idxs, lats
end

spaceagg(::UnstructuredGrid, f, A, ::Nothing) = dropagg(f, A, Coord)
# I think the best scenario is to modify `dropagg` to take in weights.
function spaceagg(::UnstructuredGrid, f, A, exw)
    error("TODO")
    w = pweights(Array(exw))
    dropagg(f, A, Coord)
end

function hemispheric_functions(::UnstructuredGrid, A)
    c = dims(A, Coord).val
    nhi, shi = hemisphere_indices(c)
    nh = A[Coord(nhi)]
    sh = A[Coord(shi)]
    oldc = dims(sh, Coord).val
    si = sortperm(oldc, by = x -> x[2], rev = true)
    newc = [SVector(x[1], abs(x[2])) for x in oldc[si]]
    di = dimindex(sh, Coord)
    newdims = Any[dims(sh)...]
    newdims[di] = Coord(newc)
    data = reverse(Array(sh); dims = di)
    sh = ClimArray(data, (newdims...,))
    return nh, sh
end

function hemispheric_means(::UnstructuredGrid, A::AbDimArray)
    nhi, shi = hemisphere_indices(A)
    nh = dropagg(mean, A[Coord(nhi)], Coord)
    sh = dropagg(mean, A[Coord(shi)], Coord)
    return nh, sh
end

function hemispheric_means(::UnstructuredGrid, A::AbDimArray, W)
	@assert size(A) == size(W)
    nhi, shi = hemisphere_indices(A)
    nh = dropagg(mean, A[Coord(nhi)], Coord, W[Coord(nhi)])
    sh = dropagg(mean, A[Coord(shi)], Coord, W[Coord(shi)])
    return nh, sh
end

"""
	hemisphere_indices(coords) → nhi, shi
Return the indices of coordinates belonging to the north and south hemispheres.
"""
function hemisphere_indices(c)
    idxs, lats = uniquelats(c)
    i = findfirst(x -> x > 0, lats)
    shi = idxs[1][1]:idxs[i][end]
    nhi = idxs[i+1][1]:idxs[end][end]
    return nhi, shi
end

latitudes(::UnstructuredGrid, A) = unique!([x[2] for x in dims(A, Coord)])

#########################################################################
# Extention of convenience indexing of `Coord`
#########################################################################
# The code here is exclusively a performance optimization that relies
# on the fact that we sort coordinates by latitude
function coord_latitudes_between(c, l1, l2)
    idxs, lats = uniquelats(c)
    i1 = searchsortedfirst(lats, l1)
    i2 = searchsortedlast(lats, l2)
    i1 > i2 && ((i1, i2) = (i2, i1)) # in case bounds are given in reverse order
    return idxs[i1][1]:idxs[i2][end]
end

# This modifies what happens on A[Coord(Lat(Between(x,y)))]
function DimensionalData.sel2indices(c::Coord, sel::Tuple{<:Lat{ <: Between}})
    l1, l2 = sel[1].val.val
    return coord_latitudes_between(c.val, l1, l2) # this is Vector{Int}
end


text:
##########
## instructions for global `.juliarc.jl`
##########

# it is necessary to include the following in the global `.juliarc.jl`
# it makes julia search the current directory for a local `juliarc.jl` file

#  if pwd() != ENV["HOME"]
    #  local_rc = joinpath(pwd(),"juliarc.jl")
    #  if isfile(local_rc)
        #  println("Loading Local juliarc.jl file....")
        #  include(local_rc)
    #  end
#  end

##########
## local `juliarc.jl`
##########

# adds current folder to the julia path so we can load the local modules
push!(LOAD_PATH,pwd())

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
function init_state_prognostic_arr!(
    balance_law,
    state::AbstractArray,
    aux::AbstractArray,
    local_geom,
    args...,
)
    FT = eltype(state)
    init_state_prognostic!(
        balance_law,
        Vars{vars_state(balance_law, Prognostic(), FT)}(state),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(aux),
        local_geom,
        args...,
    )

end

function source_arr!(
    balance_law,
    source::AbstractArray,
    state::AbstractArray,
    diffusive::AbstractArray,
    aux::AbstractArray,
    t::Real,
    direction,
)
    FT = eltype(state)
    source!(
        balance_law,
        Vars{vars_state(balance_law, Prognostic(), FT)}(source),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state),
        Vars{vars_state(balance_law, GradientFlux(), FT)}(diffusive),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(aux),
        t,
        direction,
    )
end

function flux_first_order_arr!(
    balance_law,
    flux::AbstractArray,
    state::AbstractArray,
    aux::AbstractArray,
    t::Real,
    direction,
)
    FT = eltype(state)
    flux_first_order!(
        balance_law,
        Grad{vars_state(balance_law, Prognostic(), FT)}(flux),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(aux),
        t,
        direction,
    )
end

function flux_second_order_arr!(
    balance_law,
    flux::AbstractArray,
    state::AbstractArray,
    diffusive::AbstractArray,
    hyperdiffusive::AbstractArray,
    aux::AbstractArray,
    t::Real,
)
    FT = eltype(state)
    flux_second_order!(
        balance_law,
        Grad{vars_state(balance_law, Prognostic(), FT)}(flux),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state),
        Vars{vars_state(balance_law, GradientFlux(), FT)}(diffusive),
        Vars{vars_state(balance_law, Hyperdiffusive(), FT)}(hyperdiffusive),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(aux),
        t,
    )
end

function compute_gradient_argument_arr!(
    balance_law,
    transform::AbstractArray,
    state::AbstractArray,
    aux::AbstractArray,
    t::Real,
)
    FT = eltype(state)
    compute_gradient_argument!(
        balance_law,
        Vars{vars_state(balance_law, Gradient(), FT)}(transform),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(aux),
        t,
    )
end

function compute_gradient_flux_arr!(
    balance_law,
    diffusive::AbstractArray,
    ∇transform::AbstractArray,
    state::AbstractArray,
    aux::AbstractArray,
    t::Real,
)

    FT = eltype(state)
    compute_gradient_flux!(
        balance_law,
        Vars{vars_state(balance_law, GradientFlux(), FT)}(diffusive),
        Grad{vars_state(balance_law, Gradient(), FT)}(∇transform),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(aux),
        t,
    )
end

function integral_load_auxiliary_state_arr!(
    balance_law,
    local_kernel::AbstractArray,
    state_prognostic::AbstractArray,
    state_auxiliary::AbstractArray,
)
    FT = eltype(state_auxiliary)
    integral_load_auxiliary_state!(
        balance_law,
        Vars{vars_state(balance_law, UpwardIntegrals(), FT)}(local_kernel),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary),
    )
end

function integral_set_auxiliary_state_arr!(
    balance_law,
    state_auxiliary::AbstractArray,
    local_kernel::AbstractArray,
)

    FT = eltype(state_auxiliary)
    integral_set_auxiliary_state!(
        balance_law,
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary),
        Vars{vars_state(balance_law, UpwardIntegrals(), FT)}(local_kernel),
    )
end

function reverse_integral_load_auxiliary_state_arr!(
    balance_law,
    l_T::AbstractArray,
    state::AbstractArray,
    state_auxiliary::AbstractArray,
)

    FT = eltype(state_auxiliary)
    reverse_integral_load_auxiliary_state!(
        balance_law,
        Vars{vars_state(balance_law, DownwardIntegrals(), FT)}(l_T),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary),
    )

end

function reverse_integral_set_auxiliary_state_arr!(
    balance_law,
    state_auxiliary::AbstractArray,
    l_V::AbstractArray,
)

    FT = eltype(state_auxiliary)
    reverse_integral_set_auxiliary_state!(
        balance_law,
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary),
        Vars{vars_state(balance_law, DownwardIntegrals(), FT)}(l_V),
    )
end

function transform_post_gradient_laplacian_arr!(
    balance_law,
    hyperdiffusion::AbstractArray,
    l_grad_lap::AbstractArray,
    prognostic::AbstractArray,
    auxiliary::AbstractArray,
    t,
)

    FT = eltype(prognostic)
    transform_post_gradient_laplacian!(
        balance_law,
        Vars{vars_state(balance_law, Hyperdiffusive(), FT)}(hyperdiffusion),
        Grad{vars_state(balance_law, GradientLaplacian(), FT)}(l_grad_lap),
        Vars{vars_state(balance_law, Prognostic(), FT)}(prognostic),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(auxiliary),
        t,
    )
end


text:
@testset "303.range-sum-query-immutable.jl" begin
    @test sum_range(NumArray([-2, 0, 3, -5, 2, -1]), 1, 3) == 1
    @test sum_range(NumArray([-2, 0, 3, -5, 2, -1]), 3, 6) == -1
    @test sum_range(NumArray([-2, 0, 3, -5, 2, -1]), 1, 6) == -3
end


text:
push!(LOAD_PATH, "./gridmodule/")

using BaseGrid

g1Da = grid1D("./mesh1d/mesh1")
g1Db = grid1D("./mesh1d/mesh2")
g1Dc = grid1D("./mesh1d/mesh3")

g2D = grid2D("./mesh2d/test")
g3D = grid3D("./mesh3d/ball")

gNEW = grid2D_avs("./mesh2d/len_10_s_1")

print("Interval mesh output:\n")
show(g1Da)
show(g1Db)
show(g1Dc)

print("Triangle mesh output:\n")
show(g2D)
show(gNEW)

print("Tetrahedron mesh output:\n")
show(g3D)

outvtk_mesh("output/test1Da.vtk",g1Da.geometry,g1Da.topology)
outvtk_mesh("output/test1Db.vtk",g1Db.geometry,g1Db.topology)
outvtk_mesh("output/test1Dc.vtk",g1Dc.geometry,g1Dc.topology)
outvtk_mesh("output/test2D.vtk" ,g2D.geometry ,g2D.topology)
outvtk_mesh("output/test3D.vtk" ,g3D.geometry ,g3D.topology)
outvtk_mesh("output/testNEW.vtk" ,gNEW.geometry ,gNEW.topology)

max_stars_repo_name: dpeschka/jPDE

max_issues_repo_name: dpeschka/jPDE

max_forks_repo_name: dpeschka/jPDE


text:
function show_table(app, prop)
    callback!(
        app,
        Output("show_$prop", "children"),
        Output("store_$prop", "data"),
        Input("upload_$prop", "contents"),
        State("upload_$prop", "filename"),
        State("upload_$prop", "last_modified"),
        prevent_initial_call = true
    ) do contents, filename, last_modified
        df = parse_contents(contents, filename)
        df_show = show_table(df, filename, last_modified)
        df_json = JSON.json(df)
        return df_show, df_json
    end
end

function run_simulation(app)
    callback!(
        app,
        Output("sim_complete", "displayed"),
        Output("sim_complete", "message"),
        Output("build_and_run_msg", "children"),
        Output("store_inv_on_hand", "data"),
        Output("store_inv_level", "data"),
        Output("store_inv_position", "data"),
        Output("store_ech_position", "data"),
        Output("store_inv_pipeline", "data"),
        Output("store_market_demand", "data"),
        Output("store_replenishments", "data"),
        Input("build_and_run", "n_clicks"),
        State("store_bill_of_materials", "data"),
        State("store_lead_times", "data"),
        State("store_demand", "data"),
        State("store_policy", "data"),
        State("policy_variable","value"),
        State("policy_type","value"),
        State("mode","value"),
        State("num_periods","value"),
        prevent_initial_call = true
    ) do n_clicks, bom_json, lt_json, demand_json, policy_json, policy_variable, policy_type, op_mode, num_periods
        #initialize outputs
        msg = html_div("")
        msg_txt = ""
        show_msg = false
        inv_on_hand = DataFrame()
        inv_level = DataFrame()
        inv_position = DataFrame()
        ech_position = DataFrame()
        inv_pipeline = DataFrame()
        demand = DataFrame()
        replenishments = DataFrame()

        #trigger function call
        if num_periods > 0
            try
                #parse inputs
                bom_df = DataFrame(JSON.parse(bom_json))
                lt_df = DataFrame(JSON.parse(lt_json))
                demand_df = DataFrame(JSON.parse(demand_json))
                policy_df = DataFrame(JSON.parse(policy_json))
                backlog = op_mode == "backlog"
                policy_variable = Symbol(policy_variable)
                policy_type = Symbol(policy_type)

                #build and run model
                net = build_network(lt_df)
                build_bom!(net, bom_df)
                build_demand!(net, demand_df)
                env = run_policy!(net, policy_df, policy_variable, policy_type, num_periods, backlog)

                #competion message
                msg_txt = "Simulation Complete!"
                current_time = Libc.strftime(Libc.TmStruct(Libc.TimeVal().sec))
                msg = html_div("Simulation completed at $(current_time)", style = (color = "green",))
                show_msg = true

                #prepare results
                node_dict = get_prop(net, :node_dictionary)
                node_names = Dict(val => key for (key,val) in node_dict)
                inv_on_hand = env.inv_on_hand
                inv_level = env.inv_level
                inv_position = env.inv_position
                ech_position = env.ech_position
                demand = env.demand
                inv_pipeline = env.inv_pipeline
                replenishments = env.replenishments
                for df in [inv_on_hand, inv_level, inv_position, ech_position, demand]
                    transform!(df,
                        :node => ByRow(i -> node_names[i]) => :node
                    )
                end
                for df in [inv_pipeline, replenishments]
                    transform!(df,
                        :arc => ByRow(i -> "$(node_names[i[1]]) => $(node_names[i[2]])") => :arc
                    )
                end
            catch e
                current_time = Libc.strftime(Libc.TmStruct(Libc.TimeVal().sec))
                msg = html_div("Simulation failed at $(current_time) due to an error: $e.", style = (color = "red",))
                throw(e)
            end
        else
            msg = html_div("Enter a positive integer for the number of periods to run.", style = (color = "red",))
        end

        return show_msg, msg_txt, msg,
            JSON.json(inv_on_hand), JSON.json(inv_level),
            JSON.json(inv_position), JSON.json(ech_position),
            JSON.json(inv_pipeline), JSON.json(demand),
            JSON.json(replenishments)
    end
end

function download_results(app)
    callback!(
        app,
        Output("inv_onhand_result", "data"),
        Output("inv_level_result", "data"),
        Output("inv_position_result", "data"),
        Output("ech_position_result", "data"),
        Output("inv_pipeline_result", "data"),
        Output("market_demand_result", "data"),
        Output("replenishments_result", "data"),
        Input("download", "n_clicks"),
        State("store_inv_on_hand", "data"),
        State("store_inv_level", "data"),
        State("store_inv_position", "data"),
        State("store_ech_position", "data"),
        State("store_inv_pipeline", "data"),
        State("store_market_demand", "data"),
        State("store_replenishments", "data"),
        prevent_initial_call = true
    ) do n_clicks, inv_onhand_json, inv_level_json, inv_position_json, ech_position_json, inv_pipeline_json, market_demand_json, replenishments_json
        #parse results data
        inv_onhand = DataFrame(JSON.parse(inv_onhand_json))
        inv_level = DataFrame(JSON.parse(inv_level_json))
        inv_position = DataFrame(JSON.parse(inv_position_json))
        ech_position = DataFrame(JSON.parse(ech_position_json))
        inv_pipeline = DataFrame(JSON.parse(inv_pipeline_json))
        market_demand = DataFrame(JSON.parse(market_demand_json))
        replenishments = DataFrame(JSON.parse(replenishments_json))
        replace!(replenishments.reallocated, nothing => missing)
        replace!(inv_onhand.level, nothing => Inf)
        replace!(inv_level.level, nothing => Inf)
        replace!(inv_position.level, nothing => Inf)
        replace!(ech_position.level, nothing => Inf)

        return dcc_send_string(CSV.write, inv_onhand, "inv_onhand.csv"),
            dcc_send_string(CSV.write, inv_level, "inv_level.csv"),
            dcc_send_string(CSV.write, inv_position, "inv_position.csv"),
            dcc_send_string(CSV.write, ech_position, "ech_position.csv"),
            dcc_send_string(CSV.write, inv_pipeline, "inv_pipeline.csv"),
            dcc_send_string(CSV.write, market_demand, "market_demand.csv"),
            dcc_send_string(CSV.write, replenishments, "replenishments.csv")
    end
end

function build_callbacks(app)
    show_table(app, "bill_of_materials")
    show_table(app, "lead_times")
    show_table(app, "demand")
    show_table(app, "policy")
    run_simulation(app)
    download_results(app)
end


text:
module DivPopGa

include("clusteredga.jl")
include("summary.jl")

end

max_stars_repo_name: akadal/divpopga

max_issues_repo_name: akadal/divpopga

max_forks_repo_name: akadal/divpopga