text:
isdefined(Base, :__precompile__) && __precompile__()

module NLOptControl
#TODO  enable setvalue() functionality

using JuMP
import JuMP.setRHS, JuMP.getvalue, JuMP.setvalue, JuMP.@NLexpression, JuMP.@NLobjective, JuMP.@NLparameter, JuMP.@NLconstraint, JuMP.internalmodel
using Ipopt
using KNITRO
using FastGaussQuadrature
using DataFrames # https://discourse.julialang.org/t/dataframes-0-11-released/7296
using CSV
using Interpolations

include("Base.jl")
using .Base

include("MPC_Module.jl")
using .MPC_Module

################################################################################
# Model Classfield
################################################################################
type OCP
  # general properties
  state::State              # state data
  control::Control          # control data
  tf::Any                   # final time
  t0::Any                   # initial time TODO consider getting ride of this! or replace it with n.mpc.v.t0Param
  tV::Any                   # vector for use with time varying constraints

  # boundary constraits
  X0::Array{Float64,1}      # initial state conditions
  X0_tol::Array{Float64,1}  # initial state tolerance
  x0s::Array{JuMP.Variable,1}
  XF::Array{Float64,1}      # final state conditions
  XF_tol::Array{Float64,1}  # final state tolerance
  xFs::Array{JuMP.Variable,1}

  # constant bounds on state variables
  XL::Array{Float64,1}
  XU::Array{Float64,1}

  # variables for linear bounds on state variables
  mXL::Array{Any,1}           # slope on XL -> time always starts at zero
  mXU::Array{Any,1}           # slope on XU -> time always starts at zero
  XL_var::Any                 # time varying lower bounds on states NOTE Not used currently
  XU_var::Any                 # time varying upper bounds on states NOTE Not used currently

  # constant bounds on control variables
  CL::Array{Float64,1}
  CU::Array{Float64,1}

  # ps method data
  Nck::Array{Int64,1}             # number of collocation points per interval
  Nck_cum::Array{Int64,1}         # cumulative number of points per interval
  Nck_full::Array{Int64,1}        # [0;cumsum(n.ocp.Nck+1)]
  Ni::Int64                       # number of intervals
  tau::Array{Array{Float64,1},1}  # Node points ---> Nc increasing and distinct numbers ∈ [-1,1]
  ts::Array{Array{Float64,1},1}   # time scaled based off of tau
  w::Array{Array{Float64,1},1}    # weights
  ws::Array{Array{Any,1},1}       # scaled weights
  DMatrix::Array{Array{Any,2},1}  # differention matrix
  IMatrix::Array{Array{Any,2},1}  # integration matrix

  # tm method data
  N::Int64                      # number of points in discretization
  dt::Array{Any,1}              # array of dts

  mdl::JuMP.Model               # JuMP model
  params                        # parameters for the models
  DXexpr
  NLcon# NOTE Not used currently

  # scaling factors
  XS::Array{Float64,1}           # scaling factors on states
  CS::Array{Float64,1}           # scaling factors on controls
end

# Default Constructor
function OCP()
OCP(
      State(),            # state data
      Control(),          # control data
      Any,                # final time
      0.0,                # initial time
      Any,                # optional vector for use with time varying constraints
      Float64[],          # initial state conditions
      Float64[],          # tolerances on inital state constraint
      Vector{Any}[],
      Float64[],          # final state conditions
      Vector{Any}[],
      Float64[],          # tolerances on final state constraint
      Float64[],          # XL
      Float64[],          # XU
      Any[],              # slopes on XL -> time always starts at zero
      Any[],              # slopes on XU -> time always starts at zero
      Any,                # time varying lower bounds on states
      Any,                # time varying upper bounds on states
      Float64[],          # CL
      Float64[],          # CU
      Int[],              # number of collocation points per interval
      Int[],              # Nck_cum
      Int[],              # Nck_full
      0,                  # number of intervals
      Vector{Float64}[],  # tau
      Vector{Any}[],      # ts
      Vector{Float64}[],  # weights
      Vector{Any}[],      # scaled weights
      Matrix{Any}[],      # DMatrix
      Matrix{Any}[],      # IMatrix
      0,                  # number of points in discretization
      Any[],              # array of dts
      JuMP.Model(),       # JuMP model
      Any[],
      Any[],
      Any[],
      Float64[],
      Float64[]
      )
end

type OCPFlags
  defined::Bool  # a bool to tell if define() has been called
end

function OCPFlags()
  OCPFlags(
   false
  )
end

type MPCFlags
 defined::Bool
 goalReached::Bool
 simFailed::Array{Any,1}   # a bool to indicate that the simulation failed and a symbol to indicate why
 ipDefined::Bool
 epDefined::Bool
end

function MPCFlags()
 MPCFlags(
  false,
  false,
  [false, NaN],
  false,
  false
  )
end

type Flags
  ocp::OCPFlags
  mpc::MPCFlags
end

function Flags()
  Flags(
   OCPFlags(),
   MPCFlags()
  )
end

abstract type AbstractNLOpt end
type NLOpt <: AbstractNLOpt
  # major data types
  ocp::OCP
  mpc::MPC
  s::Settings
  r::Results
  f::Flags
end

# Default Constructor
function NLOpt()
NLOpt(
  OCP(),
  MPC(),
  Settings(),
  Results(),
  Flags()
    )
end

# scripts
include("utils.jl")
include("setup.jl");
include("ps.jl");
include("diffeq.jl")

include("PrettyPlots/PrettyPlots.jl")
using .PrettyPlots

export
       # Base functions  TODO make a hamiltonian function
       evalConstraints!,
       postProcess!,
       optimize!,
       interpolateLagrange!,
       interpolateLinear!,
       interpolate_lagrange,
       opt2dfs!,

       # setup functions
       define,
       configure!,
       dynamics!,
       constraints!,
       NLExpr,

       # extra functions
       defineSolver!,

       # math functions
       integrate!,

       # optimization related functions - utils.jl
       defineTolerances!,
       linearStateTolerances!,

       # data processing  - utils.jl
       newConstraint!,
       evalMaxDualInf,
       states!,
       controls!,
       minDF,
       maxDF,
       savePlantData!,
       saveData,
       linearSpline,
       saveOptData,

       # MPC_Module.jl
       MPC,
       defineMPC!,
       initOpt!,
       defineIP!,
       mapNames!,
       simIPlant!,
       updateX0!,
       currentIPState,
       goalReached!,
       simMPC!,
       plant2dfs!,
       predictX0!,

       # Objects
       NLOpt,

       # results
       resultsDir!,  # function to make a results folder

       #JuMP macros and functions
       @NLexpression,
       @NLobjective,
       @NLparameter,
       @NLconstraint,
       setvalue,
       getvalue,
       setRHS,

       # PrettyPlots
        minDF,
        maxDF,
        plotSettings,
        _pretty_defaults,
        currentSettings,

        # NLOptControl plots
        statePlot,
        controlPlot,
        costatesPlot,
        costatesPlots,
        allPlots,
        adjust_axis,

        # MPC plots
        mpcPlot,
        tPlot,
        optPlot,

        # Plots.jl exported functions
        xlims!,
        ylims!,
        plot
end # module


text:
abstract type AbstractConstraintInfo end

get_component_name(d::AbstractConstraintInfo) = d.component_name

abstract type AbstractRangeConstraintInfo <: AbstractConstraintInfo end
abstract type AbstractRampConstraintInfo <: AbstractConstraintInfo end
abstract type AbstractStartConstraintInfo <: AbstractConstraintInfo end

""" Data Container to construct range constraints"""
struct DeviceRangeConstraintInfo <: AbstractRangeConstraintInfo
    component_name::String
    limits::MinMax
    additional_terms_ub::Vector{Symbol}
    additional_terms_lb::Vector{Symbol}
end

function DeviceRangeConstraintInfo(name::String, limits::MinMax)
    return DeviceRangeConstraintInfo(name, limits, Vector{Symbol}(), Vector{Symbol}())
end

function DeviceRangeConstraintInfo(name::String)
    return DeviceRangeConstraintInfo(
        name,
        (min = -Inf, max = Inf),
        Vector{Symbol}(),
        Vector{Symbol}(),
    )
end

get_component_name(d::DeviceRangeConstraintInfo) = d.component_name

struct DeviceTimeSeriesConstraintInfo
    bus_number::Int
    multiplier::Float64
    timeseries::Vector{Float64}
    range::DeviceRangeConstraintInfo
end

function DeviceTimeSeriesConstraintInfo(
    bus_number::Int,
    multiplier::Float64,
    timeseries,
    range::DeviceRangeConstraintInfo,
)
    ts::Vector{Float64} = timeseries
    return DeviceTimeSeriesConstraintInfo(bus_number, multiplier, ts, range)
end

get_component_name(d::DeviceTimeSeriesConstraintInfo) = get_component_name(d.range)
get_limits(d::DeviceTimeSeriesConstraintInfo) = d.range.limits
get_timeseries(d::DeviceTimeSeriesConstraintInfo) = d.timeseries

function DeviceTimeSeriesConstraintInfo(
    device::PSY.Device,
    multiplier_function::Function,
    ts_vector,
    get_constraint_values::Union{Function, Nothing} = nothing,
)
    name = PSY.get_name(device)
    bus_number = PSY.get_number(PSY.get_bus(device))
    multiplier = multiplier_function(device)
    if get_constraint_values === nothing
        range_constraint_info = DeviceRangeConstraintInfo(name)
    else
        range_constraint_info =
            DeviceRangeConstraintInfo(name, get_constraint_values(device))
    end
    return DeviceTimeSeriesConstraintInfo(
        bus_number,
        multiplier,
        ts_vector,
        range_constraint_info,
    )
end

struct DeviceMultiStartRangeConstraintsInfo <: AbstractRangeConstraintInfo
    component_name::String
    limits::PSY.Min_Max
    lag_ramp_limits::NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}}
    additional_terms_ub::Vector{Symbol}# [:R1, :R2]
    additional_terms_lb::Vector{Symbol}
end

function DeviceMultiStartRangeConstraintsInfo(
    component_name::String,
    limits::PSY.Min_Max,
    lag_ramp_limits::NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}},
)
    return DeviceMultiStartRangeConstraintsInfo(
        component_name,
        limits,
        lag_ramp_limits,
        Vector{Symbol}(),
        Vector{Symbol}(),
    )
end

function DeviceMultiStartRangeConstraintsInfo(component_name::String)
    return DeviceMultiStartRangeConstraintsInfo(
        component_name,
        (min = -Inf, max = Inf),
        (startup = Inf, shutdown = Inf),
        Vector{Symbol}(),
        Vector{Symbol}(),
    )
end

struct DeviceRampConstraintInfo <: AbstractRampConstraintInfo
    component_name::String
    limits::PSI.MinMax
    ic_power::InitialCondition
    ramp_limits::PSI.UpDown
    additional_terms_ub::Vector{Symbol}
    additional_terms_lb::Vector{Symbol}
end

function DeviceRampConstraintInfo(
    component_name::String,
    limits::PSY.Min_Max,
    ic_power::InitialCondition,
    ramp_limits::PSI.UpDown,
)
    return DeviceRampConstraintInfo(
        component_name,
        limits,
        ic_power,
        ramp_limits,
        Vector{Symbol}(),
        Vector{Symbol}(),
    )
end

function DeviceRampConstraintInfo(component_name::String, ic_power::InitialCondition)
    return DeviceRampConstraintInfo(
        component_name,
        (min = -Inf, max = Inf),
        ic_power,
        (up = Inf, down = Inf),
    )
end

get_limits(d::DeviceRampConstraintInfo) = d.limits
get_ic_power(d::DeviceRampConstraintInfo) = d.ic_power
get_ramp_limits(d::DeviceRampConstraintInfo) = d.ramp_limits
get_additional_terms_ub(d::DeviceRampConstraintInfo) = d.additional_terms_ub
get_additional_terms_lb(d::DeviceRampConstraintInfo) = d.additional_terms_lb

struct DeviceStartUpConstraintInfo <: AbstractStartConstraintInfo
    component_name::String
    time_limits::StartUpStages
    startup_types::Int
end

struct DeviceStartTypesConstraintInfo <: AbstractStartConstraintInfo
    component_name::String
    startup_types::Int
end

struct DeviceEnergyTargetConstraintInfo <: AbstractStartConstraintInfo
    component_name::String
    multiplier::Float64
    storage_target::Float64
end

struct EnergyBalanceConstraintInfo <: AbstractStartConstraintInfo
    component_name::String
    efficiency_data::InOut
    ic_energy::InitialCondition
    multiplier::Union{Nothing, Float64}
    timeseries::Union{Nothing, Vector{Float64}}
end

function EnergyBalanceConstraintInfo(; component_name, efficiency_data, ic_energy)
    return EnergyBalanceConstraintInfo(
        component_name,
        efficiency_data,
        ic_energy,
        nothing,
        nothing,
    )
end

struct ServiceRampConstraintInfo <: AbstractRampConstraintInfo
    component_name::String
    ramp_limits::UpDown
end

struct ReserveRangeConstraintInfo
    component_name::String
    limits::MinMax
    efficiency::InOut
    time_frames::Dict{Symbol, Float64}
    additional_terms_up::Vector{Symbol}
    additional_terms_dn::Vector{Symbol}
end

function ReserveRangeConstraintInfo(name::String, limits::MinMax, efficiency::InOut)
    return ReserveRangeConstraintInfo(
        name,
        limits,
        efficiency,
        Dict{Symbol, Float64}(),
        Vector{Symbol}(),
        Vector{Symbol}(),
    )
end

get_component_name(d::ReserveRangeConstraintInfo) = d.component_name
get_time_frames(v::ReserveRangeConstraintInfo) = v.time_frames
get_time_frame(v::ReserveRangeConstraintInfo, name::Symbol) = v.time_frames[name]
set_time_frame!(v::ReserveRangeConstraintInfo, value::Pair{Symbol, Float64}) =
    push!(v.time_frames, value)

ext: jl


text:
import Base: convert

struct SDL2GUIBackend <: AbstractGUIBackend
    function SDL2GUIBackend()
        SDL2.init()
        new()
    end
end

struct SDL2Window <: AbstractGUIWindow
    win::SDL2.SDLWindow
    surfaces::SDL2.CairoSDLSurface

    function SDL2Window(w,h,name)
        win = SDL2.SDLWindow(w,h,name)
        surfaces = SDL2.CairoSDLSurface(w,h)
        new(win,surfaces)
    end
end

width(win::SDL2Window)  = win.surfaces.cairo_surface.width
height(win::SDL2Window) = win.surfaces.cairo_surface.height

window(b::SDL2GUIBackend,w,h,name) = SDL2Window(w,h,name)
window(b::SDL2GUIBackend,w,h) =  window(b,w,h,"")

destroy(win::SDL2Window) = SDL2.DestroyWindow(win.win.win)#nice

function update!(win::SDL2Window)

    w, h = width(win), height(win)

    renderer = win.win.renderer
    tex = SDL2.CreateTextureFromSurface(renderer,win.surfaces.SDL_surface)
    SDL2.RenderCopy(renderer, tex, C_NULL, pointer_from_objref(SDL2.Rect(0,0,w,h)))
    SDL2.RenderPresent(renderer)
end

function rectangle(win::SDL2Window,x,y,w,h,color=RGB(0.9,0.9,0.9))
    cr = win.surfaces.cairo_context
    sdl_surf = win.surfaces.SDL_surface

    save(cr)
    set_source_rgb(cr,color.r,color.g,color.b)
    rectangle(cr,x,y,w,h)
    fill(cr)
    restore(cr)
end

function clear!(win::SDL2Window,color=RGB(1.0,1.0,1.0))
    w, h = width(win), height(win)
    rectangle(win,0,0,w,h,color)
end

function text(win::SDL2Window,x,y,str::String,color=RGB(0.0,0.0,0.0))

    cr = win.surfaces.cairo_context

    set_source_rgb(cr,color.r,color.g,color.b)
    select_font_face(cr, "monospace", Cairo.FONT_SLANT_NORMAL,Cairo.FONT_WEIGHT_NORMAL)
    set_font_size(cr, 12.0)

    move_to(cr, x, y)
    show_text(cr,str)
end

mouse_position(win::SDL2Window) = SDL2.mouse_position()

function events(win::SDL2Window)
    out = []
    ev = SDL2.event()
    while ev != nothing
        push!(out,ev)
        ev = SDL2.event()
    end

    out
end

function convert(KeyEvent,ev::SDL2.KeyboardEvent)

    action = ev._type == SDL2.KEYDOWN ? RELEASE : PRESS
    KeyEvent(ev.keysym.sym, action, ev.keysym.mod )
end

##


text:
using DataFrames
using Distributions
using StatsBase

const workertypedict = Dict(:shinsotsu => 1,
                            :tyuto => 2,
                            :chuto => 2,
                            :yakuin => 3)

function generateold(;N::Int, workertype::Int, λ::Union{Float64, Nothing}=nothing)::Vector{Int}
  if workertype == 1
    data = rand(Poisson(isnothing(λ) ? 2.5 : λ), N) .+ 22
  elseif workertype == 2
    data = rand(Poisson(isnothing(λ) ? 8. : λ), N) .+ 22
  elseif workertype == 3
    data = rand(Poisson(isnothing(λ) ? 8. : λ), N) .+ 35
  else
    error("定義されていないworkertype: $workertype")
  end
  data
end
# generateold(;N::Int, workertype::Symbol, λ::Union{Float64, Nothing}=nothing)::Vector{Int} = generateold(;N=N, workertype=get(workertypedict, workertype, -1), λ=λ)

function generategender(;N::Int, workertype::Int, p::Union{Float64, Nothing}=nothing)::Vector{Int}
  if workertype == 1
    data = rand(Bernoulli(isnothing(p) ? 0.65 : p), N)
  elseif workertype == 2
    data = rand(Bernoulli(isnothing(p) ? 0.7 : p), N)
  elseif workertype == 3
    data = rand(Bernoulli(isnothing(p) ? 0.9 : p), N)
  else
    error("定義されていないworkertype: $workertype")
  end
  data
end
# generategender(;N::Int, workertype::Symbol, p::Union{Float64, Nothing}=nothing)::Vector{Int} = generategender(;N=N, workertype=get(workertypedict, workertype, -1), p=p)

function generateyearofjoin(;df::DataFrame, stats::Union{Vector, Nothing}=nothing)::Vector{Int}
  if isnothing(stats)
    stats = [
      Dict(:min_year => 0, :max_year => 15),

      Dict(:min_year => 0, :max_year => 15),

      Dict(:min_year => 10, :max_year => 15)
    ]
  end

  function judgeyearofjoin(record, stats)
    workertype = record[:workertype]
    rand(stats[workertype][:min_year]:stats[workertype][:max_year])
  end

  [judgeyearofjoin(record, stats) for record in eachrow(df)]
end

function generategakui(;df::DataFrame, stats::Union{Vector, Nothing}=nothing)::Vector{Int}
  if isnothing(stats)
    stats = [
      Dict(:young => [0., 1.0, 0.0, 0.0],
           :middle => [0., 0.15, 0.85, 0.0],
           :senior => [0., 0.1, 0.4, 0.5]),

      Dict(:young => [0.4, 0.6, 0.0, 0.0],
           :middle => [0.4, 0.5, 0.1, 0.0],
           :senior => [0.3, 0.3, 0.3, 0.1]),

      Dict(:young => [0.05, 0.7, 0.2, 0.05],
           :middle => [0.05, 0.7, 0.2, 0.05],
           :senior => [0.05, 0.7, 0.2, 0.05])
    ]
  end

  function judgegakui(record, stats)
    old, workertype = record[:old], record[:workertype]
    if 22 ≤ old ≤ 23
      output = sample([0, 1, 2, 3], Weights(stats[workertype][:young]))
    elseif 24 ≤ old ≤ 26
      output = sample([0, 1, 2, 3], Weights(stats[workertype][:middle]))
    elseif 27 ≤ old
      output = sample([0, 1, 2, 3], Weights(stats[workertype][:senior]))
    end
    output
  end

  [judgegakui(record, stats) for record in eachrow(df)]
end

"""
scoreを算出し、scoreが高いほど勤続年数を伸ばす設定とする
昔ブラックだったが、近年ホワイトになる会社を想定
学位を積めば積むほど勤続年数が伸びる設定
入社時の年齢が若いほど勤続年数が短くなる設定
性別は勤続年数に依らない設定
勤続年数は「役員->新卒->中途」の順で短くなる設定
"""
function generateseniority(;df::DataFrame)::Vector{Int}
  function judgeseniority(record)
    yearofjoin, gakui, old, gender, workertype = record[:yearofjoin], record[:gakui], record[:old], record[:gender], record[:workertype]
    score = - yearofjoin / 10.0 + gakui / 3. + old / 20.
    if workertype == 3
      score += 5.
    elseif workertype == 1
      score += 2.
    else
      score += 1.
    end
    score
  end

  N = nrow(df)
  X = [judgeseniority(record) for record in eachrow(df)]

  dt = fit(UnitRangeTransform, X, dims=1)
  normalizedX = StatsBase.transform(dt, X)

  noisedX = normalizedX .+ rand(Uniform(-0.1, 0.1), N)
  dt = fit(UnitRangeTransform, noisedX, dims=1)
  normalizedX = StatsBase.transform(dt, noisedX)

  seniority = round.(Int, normalizedX .* 20)
  [s > y ? y : s for (s, y) in zip(seniority, df.yearofjoin)]
end

function generatecensored(;df::DataFrame)::Vector{Int}
  function judgecensored(record)
    seniority, yearofjoin = record[:seniority], record[:yearofjoin]
    seniority == yearofjoin ? 1 : 0
  end

  [judgecensored(record) for record in eachrow(df)]
end


text:
using Air
using Test
using Random.Random
import Base.IdSet
import Base.delete!
import Base.isready

@testset "Air.jl" begin
    include("api.jl")
    include("parray.jl")
    include("pset.jl")
    include("pdict.jl")
    include("lazydict.jl")
    include("pwset.jl")
    include("pwdict.jl")
    include("TX.jl")
    include("countdown.jl")
end


text:
# This file was generated, do not modify it. # hide
pyplot()
Plots.plot(sin)
Plots.savefig(joinpath(@OUTPUT, "pyplotbkend.svg")) # hide


text:
using BSON
using Test

roundtrip_equal(x) = BSON.roundtrip(x) == x

mutable struct Foo
  x
end

struct T{A}
  x::A
end

struct S end

@testset "BSON" begin

@testset "Primitive Types" begin
  @test roundtrip_equal(nothing)
  @test roundtrip_equal(1)
  @test roundtrip_equal(Dict(:a => 1,:b => 2))
  @test roundtrip_equal(UInt8[1,2,3])
  @test roundtrip_equal("b")
  @test roundtrip_equal([1,"b"])
  @test roundtrip_equal(Tuple)
end

@testset "Complex Types" begin
  @test roundtrip_equal(:foo)
  @test roundtrip_equal(Int64)
  @test roundtrip_equal(Complex{Float32})
  @test roundtrip_equal(Complex)
  @test roundtrip_equal(Array)
  @test roundtrip_equal([1,2,3])
  @test roundtrip_equal(rand(2,3))
  @test roundtrip_equal(Array{Real}(rand(2,3)))
  @test roundtrip_equal(1+2im)
  @test roundtrip_equal(Nothing[])
  @test roundtrip_equal(S[])
  @test roundtrip_equal(fill(nothing, (3,2)))
  @test roundtrip_equal(fill(S(), (1,3)))
  @test roundtrip_equal(Set([1,2,3]))
  @test roundtrip_equal(Dict("a"=>1))
  @test roundtrip_equal(T(()))
end

@testset "Circular References" begin
  x = [1,2,3]
  (x1, x2) = BSON.roundtrip((x,x))
  @test x1 == x
  @test x1 === x2

  d = Dict{Symbol,Any}(:a=>1)
  d[:d] = d
  d = BSON.roundtrip(d)
  @test d[:a] == 1
  @test d[:d] === d

  x = Foo(1)
  x.x = x
  x = BSON.roundtrip(x)
  @test x.x === x
end

@testset "Anonymous Functions" begin
  f = x -> x+1
  f2 = BSON.roundtrip(f)
  @test f2(5) == f(5)
  @test typeof(f2) !== typeof(f)

  chicken_tikka_masala(y) = x -> x+y
  f = chicken_tikka_masala(5)
  f2 = BSON.roundtrip(f)
  @test f2(6) == f(6)
  @test typeof(f2) !== typeof(f)
end

@testset "Int Literals in Type Params #41" begin
  @test BSON.constructtype(Array, (Any, Int32(1))) === Vector{Any}
  @test BSON.constructtype(Array, (Any, Int64(1))) === Vector{Any}

  @test BSON.load(joinpath(@__DIR__, "test_41_from_32bit.bson")) == Dict(:obj => Dict("name"=>[0x01, 0x02]))
  @test BSON.load(joinpath(@__DIR__, "test_41_from_64bit.bson")) == Dict(:obj => Dict("name"=>[0x01, 0x02]))
end

end


text:
# This function rotates igs GPS ephemeris from the ITRF14 reference frame to the GCRF reference frame
# using igs provided Earth rotation parameters

function rotateData(rawData::DataFrame, erps::Matrix)

    # Instantiate data matrix
    data = zeros(length(rawData.MJD), 10)

    # Read in nutation data
    #Ax = CSV.File(datadir("tab5")*"\\Ax.txt"; delim = " ", ignorerepeated = true, header = 0, types = Float64) |> Tables.matrix
    #Ay = CSV.File(datadir("tab5")*"\\Ay.txt"; delim = " ", ignorerepeated = true, header = 0, types = Float64) |> Tables.matrix
    #As = CSV.File(datadir("tab5")*"\\As.txt"; delim = " ", ignorerepeated = true, header = 0, types = Float64) |> Tables.matrix
    Ax = CSV.File(datadir("tab5", "Ax.txt"); delim = " ", ignorerepeated = true, header = 0, types = Float64) |> Tables.matrix
    Ay = CSV.File(datadir("tab5", "Ay.txt"); delim = " ", ignorerepeated = true, header = 0, types = Float64) |> Tables.matrix
    As = CSV.File(datadir("tab5", "As.txt"); delim = " ", ignorerepeated = true, header = 0, types = Float64) |> Tables.matrix

    # Loop through data and rotate
    @showprogress "Rotating ephemeris to J2000..." for i in 1:length(rawData.MJD)
        # ===== Get epoch in TAI
        tai     = TAIEpoch(rawData.MJD[i]*days, origin=:modified_julian)

        # ===== First push data that does not require rotation to data matrix
        data[i,1]   = tai.second
        data[i,2]   = parse(Float64, replace(string(rawData.Sat[i]), "G" => ""))
        data[i,6]   = rawData.Clk[i]
        data[i,10]  = rawData.Clkstd[i]

        # ===== Compute polar motion rotation matrix (W)
        ac = 0.26 # Chandler Wobble (arcseconds)
        aa = 0.12 # Annual Wobble (arcseconds)
        arcSec2Rad = π/(180.0*3600.0)

        # Get pole locations
        refMJD  = floor(data[i,1])
        dayFrac = data[i,1] % refMJD
        dayFrac < 0.5 ? (refMJD -= 0.5) : (refMJD += 0.5)
        refIdx  = findmin(abs.(view(erps,:,1) .- refMJD))[2]

        xp = (erps[refIdx,2] + dayFrac*erps[refIdx, 13]) * 10^-6 # (arcseconds)
        yp = (erps[refIdx,3] + dayFrac*erps[refIdx, 14]) * 10^-6 # (arcseconds)

        # Compute T_TT
        tt      = TTEpoch(tai)
        jdtt    = julian(tt)
        T_TT    = (jdtt.second / 86400.0 - 2451545.0) / 36525.0

        # Compute s'
        sp =  -0.0015*(ac^2/1.2 + aa^2)*T_TT

        # Compute W
        cxp = cos(xp*arcSec2Rad)
        sxp = sin(xp*arcSec2Rad)
        cyp = cos(yp*arcSec2Rad)
        syp = sin(yp*arcSec2Rad)
        csp = cos(sp*arcSec2Rad)
        ssp = sin(sp*arcSec2Rad)

        W11 = cxp*csp
        W12 = -cyp*ssp + syp*sxp*csp
        W13 = -syp*ssp - cyp*sxp*csp
        W21 = cxp*ssp
        W22 = cyp*csp + syp*sxp*ssp
        W23 = syp*csp - cyp*sxp*ssp
        W31 = sxp
        W32 = -syp*cxp
        W33 = cyp*cxp
        W   = @SMatrix [W11 W12 W13; W21 W22 W23; W31 W32 W33]

        # ===== Compute Earth rotation matrix (R)
        ut1 = UT1Epoch(tai)

        # Get Julian date
        JDut1 = julian(ut1)

        θ = 2*π*(0.7790572732640 + 1.00273781191135448*(JDut1.second/86400.0 - 2451545.0))
        R = @SMatrix [ cos(-θ) sin(-θ) 0;
                      -sin(-θ) cos(-θ) 0;
                       0 0 1]

        # ===== Compute Precession-Nutation matrix (PN)
        PN = computePrecessionNutation(T_TT, Ax, Ay, As)

        # ===== Rotate Vectors
        T       = SMatrix{3,3}(PN*R*W) # Could improve this with BLAS
        r       = @SVector [rawData.X[i], rawData.Y[i], rawData.Z[i]]
        rstd    = @SVector [rawData.Xstd[i], rawData.Ystd[i], rawData.Zstd[i]]
        data[i, 3:5] .= T*r
        data[i, 7:9] .= T*rstd
    end

    return data
end

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
import Pkg
Pkg.add("CSV")
Pkg.add("StatsBase")
Pkg.add("DataFrames")
Pkg.add("DocOpt")
Pkg.add("UMAP")
Pkg.add("RCall")
Pkg.add("MultipleTesting")
Pkg.add("MultivariateStats")
Pkg.add("Distributions")
Pkg.add("CategoricalArrays")
Pkg.add("BioProfiling")
Pkg.add("HTTP")
Pkg.add("JSON")
Pkg.add("JLD")
Pkg.add("FreqTables")
Pkg.add("DataStructures")
Pkg.add("BenchmarkTools")
Pkg.add("ParallelDataTransfer")
Pkg.add("Distances")
Pkg.add("LightGraphs")
Pkg.add("SimpleWeightedGraphs")
Pkg.add("GraphIO")
Pkg.add("ParserCombinator")

max_stars_repo_name: menchelab/BioProfilingNotebooks

max_issues_repo_name: menchelab/BioProfilingNotebooks

max_forks_repo_name: menchelab/BioProfilingNotebooks


text:
# functions for construction of Hamiltonians

## Initialize Hamiltonian
"""
    TBHamiltonian(A::Vector...; hoppingtype=Float64)

Given a tuple of vectors a1, a2, ..., initializes a tight-binding system of dimension `D`, and hoppings of type `hoppingtype`, without any orbitals or hoppings added.
"""
function TBHamiltonian(A::AbstractVector...; hoppingtype=Float64)

    dim = length(A)

    latbasis = SMatrix{dim, dim, Float64}(hcat(A...))

    return TBHamiltonian(latticebasis(A...),
                        Orbital[],
                        Hopping{length(A), hoppingtype}[],
                        [false]
                        )

end

## Add orbitals
"""
    add_orbital!(sys::TBHamiltonian, pos::AbstractVector, lm=(0, 0))

Given a label, position (pos), angular momentum (l, m) adds an Orbital(pos, (l, m)) to the tight-binding TBHamiltonian
"""
function add_orbital!(sys::TBHamiltonian, pos::AbstractVector, lm=(0, 0))

    push!(sys.orbitals, Orbital(pos, lm))
end

"""
    add_orbital!(sys::TBHamiltonian, pos::AbstractVector, wf::Symbol)

Given a label, position (pos), atomic wavefunction type adds an Orbital(pos, wf) to the tight-binding TBHamiltonian
"""
function add_orbital!(sys::TBHamiltonian, pos::AbstractVector, wf::Symbol)

    push!(sys.orbitals, Orbital(pos, wf))
end

## Add hopping
"""
    add_hopping!(sys::TBHamiltonian, from::Int, to::Int, dir::AbstractVector{Int}, val)

Given a initial orbital index (from), final orbital (to), direction (dir) and a value (val), adds a Hopping(dir, from, to, val) to the tight-binding system.
"""
function add_hopping!(sys::TBHamiltonian, from::Int, to::Int, dir::AbstractVector{Int}, val)

    push!(sys.hoppings,
        Hopping(
            from, to, toSVector(dir, dimension(sys)), convert(hopping_type(sys), val))
        )
end


text:
"""
A sub-module for a specific idea for how the basics of Breadth works.
"""


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_antdpie

"""
    ''_antdpie(;kwargs...)

An AntdPie component.

Keyword arguments:
- `id` (String; optional)
- `angleField` (String; optional)
- `annotations` (optional)
- `appendPadding` (Real | Array of Reals | String; optional)
- `autoFit` (Bool; optional)
- `className` (String; optional)
- `color` (optional): . color has the following type: String | Array of Strings | lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `colorField` (String; optional)
- `data` (Array of Dicts; required)
- `endAngle` (Real; optional)
- `height` (Real; optional)
- `innerRadius` (Real; optional)
- `label` (optional)
- `legend` (optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (String; optional)
- `meta` (optional)
- `padding` (Real | Array of Reals | String; optional)
- `pieStyle` (optional): . pieStyle has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `radius` (Real; optional)
- `renderer` (String; optional)
- `startAngle` (Real; optional)
- `statistic` (optional): . statistic has the following type: lists containing elements 'title', 'content'.
Those elements have the following types:
  - `title` (optional): . title has the following type: Bool | lists containing elements 'style', 'content', 'formatter', 'rotate', 'offsetX', 'offsetY'.
Those elements have the following types:
  - `style` (Dict; optional)
  - `content` (String; optional)
  - `formatter` (optional): . formatter has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
  - `rotate` (Real; optional)
  - `offsetX` (Real; optional)
  - `offsetY` (Real; optional)
  - `content` (optional): . content has the following type: Bool | lists containing elements 'style', 'content', 'formatter', 'rotate', 'offsetX', 'offsetY'.
Those elements have the following types:
  - `style` (Dict; optional)
  - `content` (String; optional)
  - `formatter` (optional): . formatter has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
  - `rotate` (Real; optional)
  - `offsetX` (Real; optional)
  - `offsetY` (Real; optional)
- `style` (Dict; optional)
- `tooltip` (optional)
- `width` (Real; optional)
"""
function ''_antdpie(; kwargs...)
        available_props = Symbol[:id, :angleField, :annotations, :appendPadding, :autoFit, :className, :color, :colorField, :data, :endAngle, :height, :innerRadius, :label, :legend, :loading_state, :locale, :meta, :padding, :pieStyle, :radius, :renderer, :startAngle, :statistic, :style, :tooltip, :width]
        wild_props = Symbol[]
        return Component("''_antdpie", "AntdPie", "feffery_antd_charts", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-antd-charts

max_issues_repo_name: CNFeffery/feffery-antd-charts

max_forks_repo_name: CNFeffery/feffery-antd-charts


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Orbit representation as Keplerian elements.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export kepler_to_sv

################################################################################
#                                     API
################################################################################

# Direct getters.
@inline get_epoch(k::KeplerianElements) = k.t
@inline get_a(k::KeplerianElements) = k.a
@inline get_e(k::KeplerianElements) = k.e
@inline get_i(k::KeplerianElements) = k.i
@inline get_Ω(k::KeplerianElements) = k.Ω
@inline get_ω(k::KeplerianElements) = k.ω
@inline get_f(k::KeplerianElements) = k.f

# Getters that require conversions.
@inline get_M(k::KeplerianElements) = f_to_M(k.e, k.f)
@inline get_r(k::KeplerianElements) = kepler_to_rv(k)[1]
@inline get_v(k::KeplerianElements) = kepler_to_rv(k)[2]
@inline get_rv(k::KeplerianElements) = kepler_to_rv(k)

# Conversions.
"""
    kepler_to_sv(k::KeplerianElements{Tepoch, T}) where {Tepoch, T}

Convert the Keplerian elements `k` to a state vector.
"""
function kepler_to_sv(k::KeplerianElements{Tepoch, T}) where {Tepoch, T}
    r_i, v_i = kepler_to_rv(k)
    return OrbitStateVector{Tepoch, T}(t = k.t, r = r_i, v = v_i)
end

################################################################################
#                                  Overloads
################################################################################

function getindex(k::KeplerianElements{Tepoch, T}, ::Colon) where {Tepoch, T}
    return SVector{7, T}(k.t, k.a, k.e, k.i, k.Ω, k.ω, k.f)
end

################################################################################
#                                      IO
################################################################################

function show(io::IO, k::KeplerianElements{Tepoch, T}) where {Tepoch, T}
    compact = get(io, :compact, true)
    epoch_str = sprint(print, k.t, context = :compact => compact)
    jd_str = sprint(print, jd_to_date(DateTime, k.t))
    print(io, "KeplerianElements{", Tepoch, ", ", T, "}: Epoch = $epoch_str ($jd_str)")
end

function show(
    io::IO,
    mime::MIME"text/plain",
    k::KeplerianElements{Tepoch, T}
) where {Tepoch, T}
    d2r = 180/π

    # Check if the IO supports color.
    color = get(io, :color, false)

    # Compact printing.
    compact = get(io, :compact, true)

    # Definition of colors that will be used for printing.
    b = color ? _b : ""
    d = color ? _d : ""
    g = color ? _g : ""
    y = color ? _y : ""

    # Convert the data to string.
    date_str  = sprint(print, jd_to_date(DateTime, k.t))
    epoch_str = sprint(print, k.t, context = :compact => compact)
    a_str     = sprint(print, k.a/1000, context = :compact => compact)
    e_str     = sprint(print, k.e, context = :compact => compact)
    i_str     = sprint(print, rad2deg(k.i), context = :compact => compact)
    Ω_str     = sprint(print, rad2deg(k.Ω), context = :compact => compact)
    ω_str     = sprint(print, rad2deg(k.ω), context = :compact => compact)
    f_str     = sprint(print, rad2deg(k.f), context = :compact => compact)

    # Padding to align in the floating point.
    Δepoch = findfirst('.', epoch_str)
    Δepoch === nothing && (Δepoch = length(epoch_str) + 1)

    Δa = findfirst('.', a_str)
    Δa === nothing && (Δa = length(a_str) + 1)

    Δe = findfirst('.', e_str)
    Δe === nothing && (Δe = length(e_str) + 1)

    Δi = findfirst('.', i_str)
    Δi === nothing && (Δi = length(i_str) + 1)

    ΔΩ = findfirst('.', Ω_str)
    ΔΩ === nothing && (ΔΩ = length(Ω_str) + 1)

    Δω = findfirst('.', ω_str)
    Δω === nothing && (Δω = length(ω_str) + 1)

    Δf = findfirst('.', f_str)
    Δf === nothing && (Δf = length(f_str) + 1)

    dp_pos = max(Δepoch, Δa, Δe, Δi, ΔΩ, Δω, Δf)

    epoch_str = " "^(dp_pos - Δepoch) * epoch_str
    a_str     = " "^(dp_pos - Δa) * a_str
    e_str     = " "^(dp_pos - Δe) * e_str
    i_str     = " "^(dp_pos - Δi) * i_str
    Ω_str     = " "^(dp_pos - ΔΩ) * Ω_str
    ω_str     = " "^(dp_pos - Δω) * ω_str
    f_str     = " "^(dp_pos - Δf) * f_str

    max_length = max(length(a_str),
                     length(e_str),
                     length(i_str),
                     length(Ω_str),
                     length(ω_str),
                     length(f_str))

    # Add the units.
    a_str *= " "^(max_length - length(a_str)) * " km"
    i_str *= " "^(max_length - length(i_str)) * " °"
    Ω_str *= " "^(max_length - length(Ω_str)) * " °"
    ω_str *= " "^(max_length - length(ω_str)) * " °"
    f_str *= " "^(max_length - length(f_str)) * " °"

    # Print the Keplerian elements.
    println(io, "KeplerianElements{", Tepoch, ", ", T, "}:")
    println(io, "$(b)           Epoch : $(d)" * epoch_str * " (" * date_str * ")");
    println(io, "$(b) Semi-major axis : $(d)" * a_str)
    println(io, "$(b)    Eccentricity : $(d)" * e_str)
    println(io, "$(b)     Inclination : $(d)" * i_str)
    println(io, "$(b)            RAAN : $(d)" * Ω_str)
    println(io, "$(b) Arg. of Perigee : $(d)" * ω_str)
    print(io,   "$(b)    True Anomaly : $(d)" * f_str)
end


text:
# Test Scheduling

using SimLynx
SimLynx.greet()

@event test_process(i) begin
    println("$(current_time()): Event $i occurred")
end

function main()
    println("--- Test Scheduling ---")
    println("The events should execute in order.")
    @simulation begin
        # Future events
        @schedule at 0.0 test_process(4)
        @schedule at 0.0 test_process(5)
        @schedule at 0.5 test_process(7)
        @schedule at 2.0 test_process(9)
        @schedule at 1.0 test_process(8)
        @schedule in 0.0 test_process(6)
        @schedule in 2.0 test_process(10)
        # Now events
        @schedule now test_process(2)
        @schedule now test_process(3)
        @schedule immediate test_process(1)
        start_simulation()
    end
end

main()


text:
module JuliaMaps
    include("geocode.jl")
    include("geom.jl")
    export geocode
    export geometry
    export address
    export geomgoogle
    export geomOSM
    export longitude
    export latitude
end


text:
using REPL
using Test

using Terming
const T = Terming

@testset "manual" begin

    manual_debug_mode = false
    (manual_debug_mode) && include("manual.jl")

end

@testset "Terming.jl" begin

    T.set_term!(T.PseudoTerminal(
        Base.BufferStream(), Base.BufferStream(), Base.BufferStream()
    ))

    include("terminal.jl")
    include("event.jl")
    include("parser.jl")

end


text:
using QuBase
using QuDynamics
using Base.Test
using Compat

include("propagatortest.jl")
include("qutipinterfacetest.jl")


text:
using SearchSortedNearest
using Test

@testset "SearchSortedNearest.jl" begin
    @testset "Range" begin
        for rng in [1:10, (1:10) .+ 0.4, (1:10) .+ 0.5, 1.0:10.0], item in rng
            searchsortednearest(rng, item) == floor(Int, item)
        end
        searchsortednearest(1:10, 0) == 1
        searchsortednearest(1:10, 11) == 10
        searchsortednearest(1:10, 0, rev=true) == 10
        searchsortednearest(1:10, 11, rev=true) == 1
    end
    @testset "By" begin
        struct StructWithoutOrder end
        rng = map(x -> x => StructWithoutOrder(), 1:10)
        for item in rng
            searchsortednearest(rng, item, by=first) == first(item)
            searchsortednearest(rng, item, by=first, rev=true) == 11 - first(item)
        end
    end
end


text:
let
    roadway = gen_straight_roadway(1, 1000.0)
    @test length(roadway.segments) == 1

    seg = roadway[1]
    @test length(seg.lanes) == 1

    lane = seg.lanes[1]
    @test isapprox(lane.width, DEFAULT_LANE_WIDTH)
    @test lane.curve[1] == CurvePt(VecSE2(0.0,0.0,0.0), 0.0)
    @test lane.curve[2] == CurvePt(VecSE2(1000.0,0.0,0.0), 1000.0)
end

let
    roadway = gen_straight_roadway(3, 500.0, lane_width=1.0)
    @test length(roadway.segments) == 1

    seg = roadway[1]
    @test length(seg.lanes) == 3

    lane = seg.lanes[1]
    @test isapprox(lane.width, 1.0)
    @test lane.curve[1] == CurvePt(VecSE2(0.0,0.0,0.0), 0.0)
    @test lane.curve[2] == CurvePt(VecSE2(500.0,0.0,0.0), 500.0)

    lane = seg.lanes[2]
    @test isapprox(lane.width, 1.0)
    @test lane.curve[1] == CurvePt(VecSE2(0.0,1.0,0.0), 0.0)
    @test lane.curve[2] == CurvePt(VecSE2(500.0,1.0,0.0), 500.0)

    lane = seg.lanes[3]
    @test isapprox(lane.width, 1.0)
    @test lane.curve[1] == CurvePt(VecSE2(0.0,2.0,0.0), 0.0)
    @test lane.curve[2] == CurvePt(VecSE2(500.0,2.0,0.0), 500.0)

    origin = VecSE2(1.0,2.0,deg2rad(90))
    roadway = gen_straight_roadway(1, 100.0, lane_width=1.0, origin=origin)
    @test length(roadway.segments) == 1
    lane = roadway[1].lanes[1]
    @test lane.curve[1].pos == origin
    @test isapprox(lane.curve[2].pos.x,         1.0, atol=1e-6)
    @test isapprox(lane.curve[2].pos.y,       102.0, atol=1e-6)
    @test isapprox(lane.curve[2].pos.θ, deg2rad(90), atol=1e-6)

    roadway = gen_straight_roadway(2, 100.0, lane_widths=[1.0, 2.0])
    @test length(roadway.segments) == 1
    @test isapprox(roadway[1].lanes[1].curve[1].pos.y, 0.0, atol=1e-6)
    @test isapprox(roadway[1].lanes[2].curve[1].pos.y, 0.5 + 1.0, atol=1e-6)
end

let
    roadway = gen_stadium_roadway(1, length=100.0, width=10.0,
                                  radius=25.0, ncurvepts_per_turn=2)
    @test length(roadway.segments) == 4

    seg = roadway[1]
    @test length(seg.lanes) == 1

    lane = seg.lanes[1]
    @test lane.curve[1] == CurvePt(VecSE2(100.0,0.0,0.0), 0.0)
    @test lane.curve[2] == CurvePt(VecSE2(125.0,25.0,π/2), 25*π/2)

    @test roadway[LaneTag(2,1)].curve[1] == CurvePt(VecSE2(125.0,35.0,π/2), 0.0)

    roadway = gen_stadium_roadway(4)
end


text:
using SparseArrays
using LinearAlgebra

mutable struct GPI{F<:Number, LU<:LearningUpdate} <: Learner
    ψ::Matrix{F}
    r_w::Matrix{F}
    update::LU

    num_demons::Int
    num_actions::Int

    feature_size::Int

    num_tasks::Int

    feature_projector::FeatureCreator
end

function get_weights(learner::GPI)
    return vcat(learner.ψ, learner.r_w)
end

function GPI{F}(lu, feature_size, num_demons, num_actions, num_tasks, fp, w_init) where {F<:Number}
    ψ_init = if w_init == 0 0 else 1 end
    GPI(ones(F, num_demons-num_tasks, feature_size * num_actions) * (ψ_init / feature_size),
              # ones(F, num_tasks, length(fp) * num_actions) * w_init,
              ones(F, num_tasks, length(fp)) * w_init,

              lu,
              num_demons,
              num_actions,
              feature_size,
              num_tasks,
              fp)
end

GPI(lu, feature_size, num_demons, num_actions, num_tasks, fp, w_init) =
    GPI{Float64}(lu, feature_size, num_demons, num_actions, num_tasks, fp, w_init)

update(l::GPI) = l.update

Base.size(learner::GPI) = learner.num_demons

function predict_SF(learner::GPI, ϕ::SparseVector, action)
    active_state_action = get_active_action_state_vector(ϕ, action, length(ϕ), learner.num_actions)
    learner.ψ[:, active_state_action.nzind] * active_state_action.nzval
end

function predict(learner::GPI, ϕ::SparseVector, action)

    SF = predict_SF(learner, ϕ, action)
    #Column is SF per task
    # reshaped_SF = reshape(SF, :, learner.num_tasks)
    reshaped_SF = reshape(SF, length(learner.r_w), :)

    Q = learner.r_w * reshaped_SF
    return maximum(Q)
end

(l::GPI)(ϕ, a) = predict(l, ϕ, a)
(l::GPI)(ϕ) = [predict(l, ϕ, a) for a ∈ 1:l.num_actions]

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using Test, Printf

using NCDatasets
using SingleColumnModels.DomainDecompositions
using SingleColumnModels.FiniteDifferenceGrids
using SingleColumnModels.StateVecs

n_elems_real = 10 # number of elements

grid = UniformGrid(0.0, 1.0, n_elems_real)
FT = eltype(grid)
vars = (
    (:ρ_0, DomainSubSet(gm = true)),
    (:a, DomainSubSet(gm = true, en = true, ud = true)),
    (:w, DomainSubSet(gm = true, en = true, ud = true)),
)
domain_set = DomainSet(gm = 1, en = 1, ud = 2)
q = StateVec(vars, grid, domain_set)
q_compare = deepcopy(q)
gm, en, ud, sd, al = allcombinations(q)

@testset "State Vectors" begin
    for k in over_elems(grid)
        q[:ρ_0, k] = 2.0
    end
    @test all([q[:ρ_0, k] for k in over_elems(grid)] .== 2.0)

    for k in over_elems(grid)
        q[:w, k, gm] = 3.0
    end
    @test all([q[:w, k, gm] for k in over_elems(grid)] .== 3.0)

    @test_throws BoundsError q[:w, 1, 1000] = 3.0
    @test_throws BoundsError q[:ρ_0, 1, en] = 3.0

    for k in over_elems(grid)
        ρ_0_e = q[:ρ_0, k]
        for i in alldomains(q)
            w_0_e_i = q[:w, k, i]
        end
    end
    sprint(show, q)

    assign!(q, grid, 0.0)
    k = 1
    apply_Dirichlet!(q, :ρ_0, grid, FT(2), Zmin())
    @test q[:ρ_0, k] ≈ 4
    apply_Neumann!(q, :ρ_0, grid, -2 / grid.Δz, Zmin())
    @test q[:ρ_0, k] ≈ 2

    k = grid.n_elem
    apply_Dirichlet!(q, :ρ_0, grid, FT(2), Zmax())
    @test q[:ρ_0, k] ≈ 4
    apply_Neumann!(q, :ρ_0, grid, 2 / grid.Δz, Zmax())
    @test q[:ρ_0, k] ≈ 2

    q[:ρ_0, 1] = 1.0
    q[:ρ_0, 2] = 2.0
    q[:ρ_0, 3] = 3.0
    q[:a, 1] = 1.0
    q[:a, 2] = 2.0
    q[:a, 3] = 3.0
    ρα_0_cut = q[:ρ_0, Cut(2)] .* q[:a, Cut(2)]
    @test all(ρα_0_cut .== ρα_0_cut)
    ρα_0_dual = q[:ρ_0, Dual(2)] .* q[:a, Dual(2)]
    @test all(ρα_0_dual .== ρα_0_dual)

    q[:a, 3] = 1
    @test !isnan(q)
    @test !isinf(q)
    q[:a, 3] = 1 / 0
    @test isinf(q)
    q[:a, 3] = NaN
    @test isnan(q)

    assign!(q, grid, 2.0)
    vals = [
        q[ϕ, k, i] for ϕ in var_names(q) for i in over_sub_domains(q, ϕ)
        for k in over_elems(grid)
    ]
    @test all(
        [
            q[ϕ, k, i] for ϕ in var_names(q) for i in over_sub_domains(q, ϕ)
            for k in over_elems(grid)
        ] .≈ 2.0,
    )

    assign!(q, grid, 0.0)

    assign!(q, (:w), grid, 2.0)
    @test all(
        [
            q[:w, k, i] for i in over_sub_domains(q, :w)
            for k in over_elems(grid)
        ] .≈ 2.0,
    )

    assign!(q, grid, 0.0)

    assign_real!(q, :w, grid, [2.0 for x in over_elems_real(grid)], gm)
    @test all(
        [
            q[:w, k, gm] for i in over_sub_domains(q, :w)
            for k in over_elems_real(grid)
        ] .≈ 2.0,
    )
    @test all(
        [
            q[:w, k, gm] for i in over_sub_domains(q, :w)
            for k in over_elems_ghost(grid)
        ] .≈ 0.0,
    )

    assign!(q, grid, 2.0)
    assign_ghost!(q, :w, grid, 0.0, gm)
    @test all(
        [
            q[:w, k, gm] for i in over_sub_domains(q, :w)
            for k in over_elems_real(grid)
        ] .≈ 2.0,
    )
    @test all(
        [
            q[:w, k, gm] for i in over_sub_domains(q, :w)
            for k in over_elems_ghost(grid)
        ] .≈ 0.0,
    )

    assign!(q_compare, var_names(q_compare), grid, FT(2.0))
    assign!(q, var_names(q_compare), grid, FT(1.0))

    assign!(q, grid, 0.0)
    assign_real!(q, :w, grid, [2.0 for x in over_elems_real(grid)], gm)
    extrap_0th_order!(q, :w, grid, gm)
    @test all([q[:w, k, gm] for k in over_elems(grid)] .≈ 2.0)

    for k in over_elems(grid)
        q[:a, k] = 2
    end

    assign_ghost!(q, :a, grid, 0.0)
    @test all(q[:a, k] ≈ 0.0 for k in over_elems_ghost(grid))

    extrap!(q, :a, grid)
    k = 1 + grid.n_ghost
    @test q[:a, k - 1] ≈ 2 * q[:a, k] - q[:a, k + 1]
    k = grid.n_elem - grid.n_ghost
    @test q[:a, k + 1] ≈ 2 * q[:a, k] - q[:a, k - 1]

    assign!(q, grid, 0.0)
    ki_1 = first_interior(grid, Zmin())
    ki_2 = first_interior(grid, Zmax())
    kg_1 = first_ghost(grid, Zmin())
    kg_2 = first_ghost(grid, Zmax())

    q[:w, ki_1, gm] = 1.0
    q[:w, ki_2, gm] = 1.0
    extrap!(q, :w, grid, gm)
    @test q[:w, kg_1, gm] ≈ 2.0
    @test q[:w, kg_2, gm] ≈ 2.0

    # Extrapolate field
    q[:w, ki_1, gm] = 2.0
    q[:w, kg_2, gm] = 2.0
    extrap!(q, :w, grid, Zmin(), gm)
    q[:w, ki_1, gm] ≈ 4.0
    q[:w, kg_2, gm] ≈ 2.0

    q[:w, kg_1, gm] = 2.0
    q[:w, ki_2, gm] = 2.0
    extrap!(q, :w, grid, Zmax(), gm)
    q[:w, kg_1, gm] ≈ 2.0
    q[:w, kg_2, gm] ≈ 4.0

    # Extrapolate field with specified boundary value
    q[:w, ki_1, gm] = 0.0
    q[:w, kg_2, gm] = 0.0
    extrap!(q, :w, grid, 1.0, Zmin(), gm)
    q[:w, ki_1, gm] ≈ 2.0
    q[:w, kg_2, gm] ≈ 0.0

    q[:w, kg_1, gm] = 0.0
    q[:w, ki_2, gm] = 0.0
    extrap!(q, :w, grid, 1.0, Zmax(), gm)
    q[:w, kg_1, gm] ≈ 0.0
    q[:w, kg_2, gm] ≈ 2.0

    @test var_string(q, :a, gm) == "a_gm"
    @test var_string(q, :a, en) == "a_en"
    @test var_string(q, :a, ud[1]) == "a_ud_1"
    @test var_string(q, :a, ud[2]) == "a_ud_2"
    @test var_suffix(q, :a, gm) == "_gm"
    @test var_suffix(q, :a, en) == "_en"
    @test var_suffix(q, :a, ud[1]) == "_ud_1"
    @test var_suffix(q, :a, ud[2]) == "_ud_2"
    @test_throws BoundsError var_suffix(q, :a, 1000)
end

output_root = joinpath("..", "output", "tests", "StateVecFuncs")

n_subdomains = 3 # number of sub-domains
n_elems_real = 10 # number of elements

domain_set = DomainSet(gm = 1, en = 1, ud = 1)
grid = UniformGrid(0.0, 1.0, n_elems_real)
vars = (
    (:ρ_0, DomainSubSet(gm = true)),
    (:a, DomainSubSet(gm = true, en = true, ud = true)),
    (:w, DomainSubSet(gm = true, en = true, ud = true)),
    (:ϕ, DomainSubSet(gm = true, en = true, ud = true)),
    (:ψ, DomainSubSet(gm = true, en = true, ud = true)),
)
q = StateVec(vars, grid, domain_set)
vars = (
    (:cv_ϕ_ψ, DomainSubSet(gm = true, en = true, ud = true)),
    (:TCV_ϕ_ψ, DomainSubSet(gm = true)),
)
aux = StateVec(vars, grid, domain_set)

gm, en, ud, sd, al = allcombinations(q)

@testset "IO" begin
    for k in over_elems(grid)
        q[:ρ_0, k] = 1.0 * k
        for i in al
            q[:w, k, i] = 2.0 * k
            q[:a, k, i] = 3.0 * k
        end
    end

    if !Sys.iswindows()
        nc = init_data(NetCDFWriter("q"), grid, q)
        @test isfile("q.nc")

        append_data(nc, q)
        ds = Dataset("q.nc", "r")
        @test all(ds["w_gm"][:] .≈ [q[:w, k, gm] for k in over_elems(grid)])
        @test all(ds["a_gm"][:] .≈ [q[:a, k, gm] for k in over_elems(grid)])
        rm("q.nc")

        export_state(NetCDFWriter("q"), grid, q)
        ds = Dataset("q.nc", "r")
        @test all(ds["w_gm"][:] .≈ [q[:w, k, gm] for k in over_elems(grid)])
        @test all(ds["a_gm"][:] .≈ [q[:a, k, gm] for k in over_elems(grid)])
        rm("q.nc")
    end
end

@testset "Domain average, distribute, diagnose environment, total covariance" begin
    q[:a, 2, en] = 0.25
    q[:a, 2, ud[1]] = 0.75
    q[:w, 2, en] = 2
    q[:w, 2, ud[1]] = 2
    grid_mean!(q, q, :a, :w, grid)
    @test q[:w, 2, gm] ≈ 2

    q[:w, 2, gm] = 2
    distribute!(q, grid, :w)
    @test q[:w, 2, en] ≈ 2
    for i in ud
        @test q[:w, 2, i] ≈ 2
    end

    q[:a, 2, ud[1]] = 0.1
    q[:a, 2, en] = 0.9
    q[:a, 2, gm] = 1.0

    q[:w, 2, gm] = 2
    for i in ud
        q[:w, 2, i] = 2
    end

    diagnose_environment!(q, grid, :a, (:w))
    @test q[:w, 2, en] ≈ 2

    q[:a, 2, 1] = 0.1
    q[:a, 2, 2] = 0.2
    q[:ϕ, 2, 1] = 1
    q[:ϕ, 2, 2] = 2
    q[:ψ, 2, 1] = 2
    q[:ψ, 2, 2] = 3
    aux[:cv_ϕ_ψ, 2, 1] = 1.0
    aux[:cv_ϕ_ψ, 2, 2] = 1.0
    decompose_ϕ_ψ(tcv) = tcv == :TCV_ϕ_ψ ? (:ϕ, :ψ) : error("Bad init")
    total_covariance!(aux, q, aux, :TCV_ϕ_ψ, :cv_ϕ_ψ, :a, grid, decompose_ϕ_ψ)
    @test aux[:TCV_ϕ_ψ, 2] ≈ 0.32
end

@testset "Boundary conditions" begin

    n_subdomains = 3 # number of sub-domains
    n_elems_real = 10 # number of elements

    domain_set = DomainSet(gm = 1, en = 1, ud = 1)
    grid = UniformGrid(0.0, 1.0, n_elems_real)
    FT = eltype(grid)
    vars = (
        (:a, DomainSubSet(gm = true, en = true, ud = true)),
        (:w, DomainSubSet(gm = true, en = true, ud = true)),
    )
    q = StateVec(vars, grid, domain_set)
    vars = (
        (:ρ_0, DomainSubSet(gm = true, en = true, ud = true)),
        (:K, DomainSubSet(gm = true, en = true, ud = true)),
    )
    aux = StateVec(vars, grid, domain_set)

    gm, en, ud, sd, al = allcombinations(q)
    assign!(aux, :K, grid, FT(3))
    assign!(aux, :ρ_0, grid, FT(2))
    assign!(q, :a, grid, FT(0.1))
    @test bc_source(
        q,
        grid,
        aux,
        :w,
        :ρ_0,
        :a,
        :K,
        Zmin(),
        FT(2),
        Dirichlet(),
        DiffusionAbsorbed(),
    ) ≈ FT(240)
    @test bc_source(
        q,
        grid,
        aux,
        :w,
        :ρ_0,
        :a,
        :K,
        Zmin(),
        FT(2),
        Neumann(),
        DiffusionAbsorbed(),
    ) ≈ FT(-4)

    # TODO: Validate bc_source for `AdvectionAbsorbed`
    # @test bc_source(q, grid, aux, :w, :ρ_0, :a, :K,
    #                 Zmin(), FT(2), Dirichlet(), AdvectionAbsorbed()) ≈
    # @test bc_source(q, grid, aux,
    #                 :w, :ρ_0, :a, :K,
    #                 Zmin(), FT(2), Neumann(), AdvectionAbsorbed()) ≈

end


text:
"""
    module FMM3D

Wrappers for the Flatiron Institute's FMM3D
library.

# wrappers

All N-body codes return output in an
`FMMVals` structure. See documentation
of N-body codes for details.

N-body interactions with the Helmholtz kernel
- [`hfmm3d`](@ref): ``O(N)`` fast mutlipole code
- [`h3ddir`](@ref): ``O(N^2)`` direct code

N-body interactions with the Laplace kernel
- [`lfmm3d`](@ref): ``O(N)`` fast mutlipole code
- [`l3ddir`](@ref): ``O(N^2)`` direct code

# lower level routines

For a list of lower level routines see the
documentation for [`lower_level_routs`](@ref)
"""
module FMM3D

using FMM3D_jll

export FMMVals, hfmm3d, lfmm3d, h3ddir, l3ddir, lowerlevel_routs
export besseljs3d

"""

# Available lower level routines

* [`besseljs3d`](@ref) spherical Bessel function eval
"""
function lower_level_routs() end

# fortran input/return types

Fd = Ref{Float64}
Fi = Ref{Int32}
Fc = Ref{ComplexF64}

# common input types

TFN = Union{Array{Float64},Nothing}
TCN = Union{Array{ComplexF64},Nothing}

"""
    FMMVals()

Return type for FMM and direct computation function
calls. Fields are `nothing` on return unless requested.
See individual FMM/direct computation function
documentation for specifics.
"""
mutable struct FMMVals
    pot
    grad
    hess

    pottarg
    gradtarg
    hesstarg

    # stokes-specific return values

    pre
    pretarg

    # maxwell-specific return values
    e
    h
    grade
    gradh

    etarg
    htarg
    gradetarg
    gradhtarg
end

function FMMVals()
    return FMMVals(nothing,nothing,nothing,
                   nothing,nothing,nothing,
                   nothing,nothing,
                   nothing,nothing,nothing,nothing,
                   nothing,nothing,nothing,nothing)
end

"""
```julia
    vals = hfmm3d(eps,zk,sources;charges=nothing,dipvecs=nothing,
                  targets=nothing,pg=0,pgt=0,nd=1)
```

This function computes the N-body Helmholtz interactions
in three dimensions where the interaction kernel is given
by ``e^{ikr}/r`` and its gradients. This is the
``O(N)`` fast multipole code which computes the interactions
to the requested precision.

```math
 u(x) = \\sum_{j=1}^{N} c_{j} \\frac{e^{ik \\|x-x_{j}\\|}}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{e^{ik \\|x-x_{j}\\|}}{\\|x-x_{j}\\|}
\\right)  \\, ,
```

where ``c_{j}`` are the charge densities,
      ``v_{j}`` are the dipole orientation vectors, and
      ``x_{j}`` are the source locations.
      When ``x=x_{m}``, the term corresponding to
      ``x_{m}`` is dropped from the sum

# Input

* `eps::Float64` precision requested
* `zk::ComplexF64` Helmholtz parameter
* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)
* `charges::Array{ComplexF64}` size (nd,n) or (n) charge densities (c_{j})
* `dipvecs::Array{ComplexF64}` size (nd,3,n) or (3,n)) dipole orientation vectors (``v_{j}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `pg::Integer` source eval flag.
    + Potential (``u``) at sources evaluated if `pg == 1`.
    + Potential and gradient (``\\nabla u``) at sources evaluated if `pg == 2`
* `pgt::Integer` target eval flag.
    + Potential at targets evaluated if `pgt == 1`.
    + Potential and gradient at targets evaluated if `pgt == 2`
* `nd::Integer` number of densities

# Output

`vals::FMMVals` with the fields

* `vals.pot::Array{ComplexF64}` size (nd,n) or (n) potential at source locations if requested
* `vals.grad::Array{ComplexF64}` size (nd,3,n) or (3,n) gradient at source locations if requested
* `vals.pottarg::Array{ComplexF64}` size (nd,nt) or (nt) potential at target locations if requested
* `vals.gradtarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) gradient at target locations if requested
"""
function hfmm3d(eps::Float64,zk::Union{Float64,ComplexF64},
                sources::Array{Float64};
                charges::TCN=nothing,dipvecs::TCN=nothing,
                targets::TFN=nothing,pg::Integer=0,pgt::Integer=0,
                nd::Integer=1)

    zk = complex(zk)

    # default values

    vals = FMMVals()

    ifcharge = 0
    ifdipole = 0

    zero = ComplexF64(0)

    pot = zero
    grad = zero
    hess = zero
    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    # check inputs

    @assert size(sources,1) == 3
    @assert nd >= 0

    if (pg > 2 || pg < 0); @warn "flag pg not in expected range" end
    if (pgt > 2 || pgt < 0); @warn "flag pgt not in expected range" end

    n = div(length(sources),3)

    nt = 0

    if targets != nothing
        @assert size(targets,1) == 3
        nt = div(length(targets),3)
    else
        targets = 0.0
    end

    if charges != nothing
        @assert div(length(charges),nd) == n
        ifcharge = 1
    else
        charges = zero
        ifcharge = 0
    end

    if dipvecs != nothing
        @assert div(length(dipvecs),nd) == n*3
        ifdipole = 1
    else
        dipvecs = zero
        ifdipole = 0
    end

    if ifcharge == 0 && ifdipole == 0
        @warn "no charges or dipoles provided, doing nothing"
        return vals
    end

    if (pg != 1 && pg != 2) && (pgt != 1 && pgt != 2)
        @warn "no output requested, doing nothing"
        return vals
    end

    if (pgt == 1 || pgt == 2) && targets == nothing
        @warn "target values requested but no targets provided"
    end

    # allocate memory for return values

    if pg == 1 || pg == 2
        if nd > 1
            pot = zeros(ComplexF64,nd,n)
        else
            pot = zeros(ComplexF64,n)
        end
    end

    if pg == 2
        if nd > 1
            grad = zeros(ComplexF64,nd,3,n)
        else
            grad = zeros(ComplexF64,3,n)
        end
    end

    if pgt == 1 || pgt == 2
        if nd > 1
            pottarg = zeros(ComplexF64,nd,nt)
        else
            pottarg = zeros(ComplexF64,nt)
        end
    end

    if pgt == 2
        if nd > 1
            gradtarg = zeros(ComplexF64,nd,3,nt)
        else
            gradtarg = zeros(ComplexF64,3,nt)
        end
    end

    # actually call the function
    #
    # fortran calling sequence:
    # subroutine hfmm3d(nd,eps,zk,nsource,source,ifcharge,
    #     $    charge,ifdipole,dipvec,ifpgh,pot,grad,hess,ntarg,
    #     $    targ,ifpghtarg,pottarg,gradtarg,hesstarg)
    #

    ccall((:hfmm3d_,libfmm3d),Cvoid,(Fi,Fd,Fc,Fi,Fd,Fi,Fc,
                                    Fi,Fc,Fi,Fc,Fc,Fc,Fi,
                                    Fd,Fi,Fc,Fc,Fc),
          nd,eps,zk,n,sources,ifcharge,charges,ifdipole,
          dipvecs,pg,pot,grad,hess,nt,targets,pgt,
          pottarg,gradtarg,hesstarg)

    # load requested values

    if pg == 1 || pg == 2; vals.pot = pot end
    if pg == 2; vals.grad = grad end
    if pgt == 1 || pgt == 2; vals.pottarg = pottarg end
    if pgt == 2; vals.gradtarg = gradtarg end

    return vals

end

"""
```julia
    vals = h3ddir(zk,sources,targets;charges=nothing,
                    dipvecs=nothing,pgt=0,nd=1,
                    thresh=1e-16)
```

This function computes the N-body Helmholtz interactions
in three dimensions where the interaction kernel is given
by ``e^{ikr}/r`` and its gradients. This is the
``O(N^2)`` direct evaluation code.

```math
 u(x) = \\sum_{j=1}^{N} c_{j} \\frac{e^{ik \\|x-x_{j}\\|}}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{e^{ik \\|x-x_{j}\\|}}{\\|x-x_{j}\\|}
\\right)  \\, ,
```

where ``c_{j}`` are the charge densities,
      ``v_{j}`` are the dipole orientation vectors, and
      ``x_{j}`` are the source locations.
      When ``x=x_{m}``, the term corresponding to
      ``x_{m}`` is dropped from the sum

# Input

* `zk::ComplexF64` Helmholtz parameter
* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `charges::Array{ComplexF64}` size (nd,n) or (n) charge densities (c_{j})
* `dipvecs::Array{ComplexF64}` size (nd,3,n) or (3,n)) dipole orientation vectors (``v_{j}``)
* `pgt::Integer` target eval flag.
    + Potential at targets evaluated if `pgt == 1`.
    + Potential and gradient at targets evaluated if `pgt == 2`
* `nd::Integer` number of densities
* `thresh::Float64` threshold for ignoring interactions when ``\\|x-x_{j}\\| \\leq thresh``

# Output

`vals::FMMVals` with the fields

* `vals.pottarg::Array{ComplexF64}` size (nd,nt) or (nt) potential at target locations if requested
* `vals.gradtarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) gradient at target locations if requested
"""
function h3ddir(zk::Union{ComplexF64,Float64},sources::Array{Float64},
                targets::Array{Float64};
                charges::TCN=nothing,dipvecs::TCN=nothing,
                pgt::Integer=0,nd::Integer=1,
                thresh::Float64=1e-16)

    zk = complex(zk)

    # default values

    vals = FMMVals()

    ifcharge = 0
    ifdipole = 0

    zero = ComplexF64(0)

    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    # check inputs

    @assert size(sources,1) == 3
    @assert size(targets,1) == 3
    @assert nd >= 0

    if (pgt > 2 || pgt < 0); @warn "flag pgt not in expected range" end

    n = div(length(sources),3)
    nt = div(length(targets),3)

    if charges != nothing
        @assert div(length(charges),nd) == n
        ifcharge = 1
    else
        charges = zero
        ifcharge = 0
    end

    if dipvecs != nothing
        @assert div(length(dipvecs),nd) == n*3
        ifdipole = 1
    else
        dipvecs = zero
        ifdipole = 0
    end

    if ifcharge == 0 && ifdipole == 0
        @warn "no charges or dipoles provided, doing nothing"
        return vals
    end

    if (pgt != 1 && pgt != 2)
        @warn "no output requested, doing nothing"
        return vals
    end

    # allocate memory for return values

    if pgt == 1 || pgt == 2
        if nd > 1
            pottarg = zeros(ComplexF64,nd,nt)
        else
            pottarg = zeros(ComplexF64,nt)
        end
    end

    if pgt == 2
        if nd > 1
            gradtarg = zeros(ComplexF64,nd,3,nt)
        else
            gradtarg = zeros(ComplexF64,3,nt)
        end
    end

    # dispatch to appropriate wrapper

    if pgt == 1
        if ifcharge == 1
            if ifdipole == 1
                h3ddirectcdp!(nd,zk,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,thresh)
                vals.pottarg = pottarg
            else
                h3ddirectcp!(nd,zk,sources,charges,
                             n,targets,nt,pottarg,
                             thresh)
                vals.pottarg = pottarg
            end
        else
            if ifdipole == 1
                h3ddirectdp!(nd,zk,sources,
                             dipvecs,n,targets,nt,
                             pottarg,thresh)
                vals.pottarg = pottarg
            end
        end
    elseif pgt == 2
        if ifcharge == 1
            if ifdipole == 1
                h3ddirectcdg!(nd,zk,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            else
                h3ddirectcg!(nd,zk,sources,charges,
                             n,targets,nt,pottarg,
                             gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            end
        else
            if ifdipole == 1
                h3ddirectdg!(nd,zk,sources,
                              dipvecs,n,targets,nt,
                             pottarg,gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            end
        end
    end

    return vals

end

function h3ddirectcp!(nd,zk,sources,charges,n,targets,nt,
                      pot,thresh)
    ccall((:h3ddirectcp_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,
                                          Fd,Fi,Fc,Fd),
          nd,zk,sources,charges,n,targets,nt,pot,thresh)
    return
end

function h3ddirectdp!(nd,zk,sources,dipvecs,n,targets,nt,
                      pot,thresh)
    ccall((:h3ddirectdp_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,
                                          Fd,Fi,Fc,Fd),
          nd,zk,sources,dipvecs,n,targets,nt,pot,thresh)
    return
end

function h3ddirectcdp!(nd,zk,sources,charges,dipvecs,n,targets,
                       nt,pot,thresh)
    ccall((:h3ddirectcdp_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fc,Fi,
                                          Fd,Fi,Fc,Fd),
          nd,zk,sources,charges,dipvecs,n,targets,nt,pot,thresh)
    return
end

function h3ddirectcg!(nd,zk,sources,charges,n,targets,nt,
                      pot,grad,thresh)
    ccall((:h3ddirectcg_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,
                                          Fd,Fi,Fc,Fc,Fd),
          nd,zk,sources,charges,n,targets,nt,pot,grad,thresh)
    return
end

function h3ddirectdg!(nd,zk,sources,dipvecs,n,targets,nt,
                      pot,grad,thresh)
    ccall((:h3ddirectdg_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,
                                          Fd,Fi,Fc,Fc,Fd),
          nd,zk,sources,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end

function h3ddirectcdg!(nd,zk,sources,charges,dipvecs,n,targets,
                       nt,pot,grad,thresh)
    ccall((:h3ddirectcdg_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fc,Fi,
                                          Fd,Fi,Fc,Fc,Fd),
          nd,zk,sources,charges,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end

"""
```julia
    vals = lfmm3d(eps,sources;charges=nothing,dipvecs=nothing,
                  targets=nothing,pg=0,pgt=0,nd=1)
```

This function computes the N-body Laplace interactions
in three dimensions where the interaction kernel is given
by ``\$1/r\$`` and its gradients.  This is the
``O(N)`` fast multipole code which computes the interactions
to the requested precision.

```math
u(x) = \\sum_{j=1}^{N} c_{j} / \\|x-x_{j}\\| + v_{j}
\\cdot \\nabla( 1/\\|x-x_{j}\\|)  \\, ,
u(x) = \\sum_{j=1}^{N} c_{j} \\frac{1}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{1}{\\|x-x_{j}\\|}
\\right)  \\, ,
```

where ``c_{j}`` are the charge densities,
      ``v_{j}`` are the dipole orientation vectors, and
      ``x_{j}`` are the source locations.
      When ``x=x_{m}``, the term corresponding to
      ``x_{m}`` is dropped from the sum

# Input

# Input

* `eps::Float64` precision requested
* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)
* `charges::Array{Float64}` size (nd,n) or (n) charge densities (c_{j})
* `dipvecs::Array{Float64}` size (nd,3,n) or (3,n)) dipole orientation vectors (``v_{j}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `pg::Integer` source eval flag.
    + Potential (``u``) at sources evaluated if `pg == 1`.
    + Potential and gradient (``\\nabla u``) at sources evaluated if `pg == 2`
* `pgt::Integer` target eval flag.
    + Potential at targets evaluated if `pgt == 1`.
    + Potential and gradient at targets evaluated if `pgt == 2`
* `nd::Integer` number of densities

# Output

`vals::FMMVals` with the fields

Note that the Hessian is returned as a length
6 vector at each point with the second derivatives
in the order: ``\\partial_{xx}``, ``\\partial_{yy}``,
``\\partial_{zz}``, ``\\partial_{xy}``,
``\\partial_{xz}``, ``\\partial_{yz}``.

* `vals.pot::Array{Float64}` size (nd,n) potential at source locations if requested
* `vals.grad::Array{Float64}` size (nd,3,n) gradient at source locations if requested
* `vals.hess::Array{Float64}` size (nd,6,n) Hessian at source locations if requested
* `vals.pottarg::Array{Float64}` size (nd,nt) potential at target locations if requested
* `vals.gradtarg::Array{Float64}` size (nd,3,nt) gradient at target locations if requested
* `vals.hesstarg::Array{Float64}` size (nd,6,nt) Hessian at target locations if requested
"""
function lfmm3d(eps::Float64,sources::Array{Float64};charges::TFN=nothing,
                dipvecs::TFN=nothing,targets::TFN=nothing,
                pg::Integer=0,pgt::Integer=0,nd::Integer=1)

    # default values

    vals = FMMVals()

    ifcharge = 0
    ifdipole = 0

    zero = Float64(0)

    pot = zero
    grad = zero
    hess = zero
    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    # check inputs

    @assert size(sources,1) == 3
    @assert nd >= 0

    if (pg > 3 || pg < 0); @warn "flag pg not in expected range" end
    if (pgt > 3 || pgt < 0); @warn "flag pgt not in expected range" end

    n = div(length(sources),3)

    nt = 0

    if targets != nothing
        @assert size(targets,1) == 3
        nt = div(length(targets),3)
    else
        targets = 0.0
    end

    if charges != nothing
        @assert div(length(charges),nd) == n
        ifcharge = 1
    else
        charges = zero
        ifcharge = 0
    end

    if dipvecs != nothing
        @assert div(length(dipvecs),nd) == n*3
        ifdipole = 1
    else
        dipvecs = zero
        ifdipole = 0
    end

    if ifcharge == 0 && ifdipole == 0
        @warn "no charges or dipoles provided, doing nothing"
        return vals
    end

    if ((pg != 1 && pg != 2 && pg !=3) &&
        (pgt != 1 && pgt != 2 && pgt != 3))
        @warn "no output requested, doing nothing"
        return vals
    end

    if (pgt == 1 || pgt == 2 || pgt == 3) && targets == nothing
        @warn "target values requested but no targets provided"
    end

    # allocate memory for return values

    if pg == 1 || pg == 2 || pg == 3
        if nd > 1
            pot = zeros(Float64,nd,n)
        else
            pot = zeros(Float64,n)
        end
    end

    if pg == 2 || pg == 3
        if nd > 1
            grad = zeros(Float64,nd,3,n)
        else
            grad = zeros(Float64,3,n)
        end
    end

    if pg == 3
        if nd > 1
            hess = zeros(Float64,nd,6,n)
        else
            hess = zeros(Float64,6,n)
        end
    end

    if pgt == 1 || pgt == 2 || pgt == 3
        if nd > 1
            pottarg = zeros(Float64,nd,nt)
        else
            pottarg = zeros(Float64,nt)
        end
    end

    if pgt == 2 || pgt == 3
        if nd > 1
            gradtarg = zeros(Float64,nd,3,nt)
        else
            gradtarg = zeros(Float64,3,nt)
        end
    end

    if pgt == 3
        if nd > 1
            hesstarg = zeros(Float64,nd,6,nt)
        else
            hesstarg = zeros(Float64,6,nt)
        end
    end

    # actually call the function
    #
    # fortran calling sequence:
    # subroutine lfmm3d(nd,eps,nsource,source,ifcharge,
    #     $    charge,ifdipole,dipvec,ifpgh,pot,grad,hess,ntarg,
    #     $    targ,ifpghtarg,pottarg,gradtarg,hesstarg)
    #

    ccall((:lfmm3d_,libfmm3d),Cvoid,(Fi,Fd,Fi,Fd,Fi,Fd,
                                    Fi,Fd,Fi,Fd,Fd,Fd,Fi,
                                    Fd,Fi,Fd,Fd,Fd),
          nd,eps,n,sources,ifcharge,charges,ifdipole,
          dipvecs,pg,pot,grad,hess,nt,targets,pgt,
          pottarg,gradtarg,hesstarg)

    # load requested values

    if pg == 1 || pg == 2 || pg == 3; vals.pot = pot end
    if pg == 2 || pg == 3; vals.grad = grad end
    if pg == 3; vals.hess = hess end
    if pgt == 1 || pgt == 2 || pgt == 3; vals.pottarg = pottarg end
    if pgt == 2 || pgt == 3; vals.gradtarg = gradtarg end
    if pgt == 3; vals.hesstarg = hesstarg end

    return vals

end

"""
```julia
    vals = l3ddir(sources,targets;charges=nothing,
                dipvecs=nothing,pgt=0,nd=1,
                thresh=1e-16)
```

This function computes the N-body Laplace interactions
in three dimensions where the interaction kernel is given
by ``1/r`` and its gradients.  This is the
``O(N^2)`` direct evaluation code.

```math
 u(x) = \\sum_{j=1}^{N} c_{j} \\frac{1}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{1}{\\|x-x_{j}\\|}
\\right)  \\, ,
```

where ``c_{j}`` are the charge densities,
      ``v_{j}`` are the dipole orientation vectors, and
      ``x_{j}`` are the source locations.
      When ``x=x_{m}``, the term corresponding to
      ``x_{m}`` is dropped from the sum

# Input

* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `charges::Array{ComplexF64}` size (nd,n) or (n) charge densities (c_{j})
* `dipvecs::Array{ComplexF64}` size (nd,3,n) or (3,n)) dipole orientation vectors (``v_{j}``)
* `pgt::Integer` target eval flag.
    + Potential at targets evaluated if `pgt == 1`.
    + Potential and gradient at targets evaluated if `pgt == 2`
* `nd::Integer` number of densities
* `thresh::Float64` threshold for ignoring interactions when ``\\|x-x_{j}\\| \\leq thresh``

# Output

`vals::FMMVals` with the fields

* `vals.pottarg::Array{ComplexF64}` size (nd,nt) or (nt) potential at target locations if requested
* `vals.gradtarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) gradient at target locations if requested
"""
function l3ddir(sources::Array{Float64},targets::Array{Float64};
                charges::TFN=nothing,dipvecs::TFN=nothing,
                pgt::Integer=0,nd::Integer=1,
                thresh::Float64=1e-16)

    # default values

    vals = FMMVals()

    ifcharge = 0
    ifdipole = 0

    zero = Float64(0)

    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    # check inputs

    @assert size(sources,1) == 3
    @assert size(targets,1) == 3
    @assert nd >= 0

    if (pgt > 2 || pgt < 0); @warn "flag pgt not in expected range" end

    n = div(length(sources),3)
    nt = div(length(targets),3)

    if charges != nothing
        @assert div(length(charges),nd) == n
        ifcharge = 1
    else
        charges = zero
        ifcharge = 0
    end

    if dipvecs != nothing
        @assert div(length(dipvecs),nd) == n*3
        ifdipole = 1
    else
        dipvecs = zero
        ifdipole = 0
    end

    if ifcharge == 0 && ifdipole == 0
        @warn "no charges or dipoles provided, doing nothing"
        return vals
    end

    if (pgt != 1 && pgt != 2)
        @warn "no output requested, doing nothing"
        return vals
    end

    # allocate memory for return values

    if pgt == 1 || pgt == 2
        if nd > 1
            pottarg = zeros(Float64,nd,nt)
        else
            pottarg = zeros(Float64,nt)
        end
    end

    if pgt == 2
        if nd > 1
            gradtarg = zeros(Float64,nd,3,nt)
        else
            gradtarg = zeros(Float64,3,nt)
        end
    end

    # dispatch to appropriate wrapper

    if pgt == 1
        if ifcharge == 1
            if ifdipole == 1
                l3ddirectcdp!(nd,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,thresh)
                vals.pottarg = pottarg
            else
                l3ddirectcp!(nd,sources,charges,
                             n,targets,nt,pottarg,
                             thresh)
                vals.pottarg = pottarg
            end
        else
            if ifdipole == 1
                l3ddirectdp!(nd,sources,
                             dipvecs,n,targets,nt,
                             pottarg,thresh)
                vals.pottarg = pottarg
            end
        end
    elseif pgt == 2
        if ifcharge == 1
            if ifdipole == 1
                l3ddirectcdg!(nd,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            else
                l3ddirectcg!(nd,sources,charges,
                             n,targets,nt,pottarg,
                             gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            end
        else
            if ifdipole == 1
                l3ddirectdg!(nd,sources,
                              dipvecs,n,targets,nt,
                             pottarg,gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            end
        end
    end

    return vals

end

function l3ddirectcp!(nd,sources,charges,n,targets,nt,
                      pot,thresh)
    ccall((:l3ddirectcp_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd),
          nd,sources,charges,n,targets,nt,pot,thresh)
    return
end

function l3ddirectdp!(nd,sources,dipvecs,n,targets,nt,
                      pot,thresh)
    ccall((:l3ddirectdp_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd),
          nd,sources,dipvecs,n,targets,nt,pot,thresh)
    return
end

function l3ddirectcdp!(nd,sources,charges,dipvecs,n,targets,
                       nt,pot,thresh)
    ccall((:l3ddirectcdp_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd),
          nd,sources,charges,dipvecs,n,targets,nt,pot,thresh)
    return
end

function l3ddirectcg!(nd,sources,charges,n,targets,nt,
                      pot,grad,thresh)
    ccall((:l3ddirectcg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd,Fd),
          nd,sources,charges,n,targets,nt,pot,grad,thresh)
    return
end

function l3ddirectdg!(nd,sources,dipvecs,n,targets,nt,
                      pot,grad,thresh)
    ccall((:l3ddirectdg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd,Fd),
          nd,sources,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end

function l3ddirectcdg!(nd,sources,charges,dipvecs,n,targets,
                       nt,pot,grad,thresh)
    ccall((:l3ddirectcdg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd,Fd),
          nd,sources,charges,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end

"""
    function fj, fjder = besseljs3d(nterms,z;scale=1.0,ifder=0)

This subroutine evaluates the first `nterms` spherical Bessel
functions, and if requested, their derivatives.
It incorporates a scaling parameter `scale` so that

      	fjs_n(z)=j_n(z)/SCALE^n
      	fjder_n(z)=\\frac{\\partial fjs_n(z)}{\\partial z}

# Input

* `nterms::Integer` order of expansion of output array `fjs`
* `z::ComplexF64` argument of the spherical Bessel functions
* `scale::Float64` scaling factor
* `ifder::Integer1` flag indicating whether to calculate `fjder`
      	          0	NO
      	          1	YES
OUTPUT:

* `fjs::Array{ComplexF64}` array of length `nterms+1` of scaled Bessel functions.
* `fjder::Array{ComplexF64}` array of derivatives of scaled Bessel functions, if requested.
"""
function besseljs3d(nterms::Integer,z::ComplexF64;scale::Float64=1.0,
                    ifder::Integer=0)

    @assert (nterms >= 0)
    if (ifder !=0 && ifder != 1)
        @warn "unexpected value in ifder, no ders computed"
    end

    fjs = Array{ComplexF64}(undef,nterms+1)
    fjder = ComplexF64(0)

    if ifder == 1
        fjder = Array{ComplexF64}(undef,nterms+1)
    end

    # fortran interface
    # subroutine besseljs3d(nterms,z,scale,fjs,ifder,fjder)

    ccall((:besseljs3d_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,Fc),
          nterms,z,scale,fjs,ifder,fjder)

    if ifder != 1
        fjder = nothing
    end

    return fjs, fjder

end

end # module


text:
## PARTE 1
cd(@__DIR__)


text:
function initializeWorld()
    Δs_lim = 1.0 # Distance limit per time unit

    s_t_min = s_0
    function execute(u_hat_t::Float64)
        Δs_t = clamp(u_hat_t, -Δs_lim, Δs_lim) # Allowed displacement

        # Compute next state
        s_t = s_t_min + Δs_t

        # Reset state for next step
        s_t_min = s_t
    end

    s_t = s_0
    observe() = s_t

    return (execute, observe)
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
# # Exercise: compressing an image using the SVD

using Pkg
Pkg.add("Images")
Pkg.add("ImageMagick")

using Images, LinearAlgebra

banana = load("Courses/Introduction to Julia/images/banana.jpg")

# Images work just like arrays, but display specially

size(banana)

banana[60,50]

gb = Gray.(banana)

channelview(gb)

# The SVD factorizes a matrix A such that
#     A == U * S * transpose(V)
# where U and V are unitary, and S is diagonal.
#
# You can think of this as a sum of outer products (column of U times row
# of V), where each one is scaled by a value in diag(S).
# Each outer product is a full matrix (image).
# So, the ones scaled by larger singular values are more important,
# and we can throw away ones scaled by smaller values.

# Here's what it looks like keeping 30, 10, 5, and 3 values:

load("Courses/Introduction to Julia/images/banana_30svals.png")

load("Courses/Introduction to Julia/images/banana_10svals.png")

load("Courses/Introduction to Julia/images/banana_5svals.png")

load("Courses/Introduction to Julia/images/banana_3svals.png")

# ## Problem statement
#
# Write a function `compress_image`. Its arguments should be an image and the
# factor by which you want to compress the image. A compressed grayscale image
# should display when compress_image is called.
# For example,
#     `compress_image("images/104_100.jpg", 33)`
# will return a compressed image of a grayscale banana built using 3 singular
# values. (This image has 100 singular values, so use fld(100, 33) to determine
# how many singular values to keep. fld performs "floor" division.)
#
# Hints:
# - Perform the SVD on the `channelview` of the image.
# - Read the documentation for `svd`.

# Remember: you can index arrays with `:` and ranges `a:b`:

A = [ i+j for i = 1:10, j = 1:10 ]

A[2:4, :]
A[1:2, 1:2]
A[:, 4:5]

max_stars_repo_path: Courses/Introduction to Julia/05. Image compression exercise.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
module Skew

mvnrnd(n,mu,C) = sqrtm(C)*randn(size(C,1),n) .+ mu

# Generate n random samples from SkewNormal(Omega,alpha).
#    Omega = d-by-d "correlation" matrix (a covariance matrix with ones on the diagonal).
#    alpha = length d vector of skew parameters.
# (This implements Proposition 1 of Azzalini & Capitanio (1999).)
function skewrndNd(n,Omega,alpha)
    delta = (Omega*alpha)/sqrt(1+dot(alpha,Omega*alpha))
    x = mvnrnd(n,0,[1 delta'; delta Omega])
    return sign(x[1,:]).*x[2:end,:]
end

# Generate n random samples from SkewNormal(xi,C,alpha).
#    xi = length d vector of location parameters (it is not necessarily the mean of the resulting samples).
#    C = d-by-d posdef matrix (it is not necessarily the covariance matrix of the resulting samples).
#    alpha = length d vector of skew parameters.
# (This implements the location-scale extension in Section 5.1 of Azzalini & Capitanio (1999).)
function skewrndNd(n,xi,C,alpha)
    s = sqrt(diag(C))
    Omega = (1.0./s) .* C .* (1.0./s)'
    return xi .+ s.*skewrndNd(n,Omega,alpha)
end

# Compute the mean of Z ~ SkewNormal(Omega,alpha).
# (See Equation (4) of Azzalini & Capitanio (1999).)
mu(Omega,alpha) = (delta = (Omega*alpha)/sqrt(1+dot(alpha,Omega*alpha)); sqrt(2/pi)*delta)

# Compute the covariance matrix of Z ~ SkewNormal(Omega,alpha).
# (See Equation (4) of Azzalini & Capitanio (1999).)
Cov(Omega,alpha) = (m = mu(Omega,alpha); Omega-m*m')

# Generate n random samples from SkewNormal(Omega,alpha) and transform them so that each entry has zero mean and unit variance.
skewrndNormalized(n,Omega,alpha) = (skewrndNd(n,Omega,alpha) .- mu(Omega,alpha))./sqrt(diag(Cov(Omega,alpha)))

end

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
module TextUnidecode

const table = Vector{Vector{String}}(undef, 256)

"""
Convert non-ascii characters to "good enough" ascii.
```jldoctest
julia> unidecode("南无阿弥陀佛")
"Nan Wu A Mi Tuo Fo"

julia> unidecode("あみだにょらい")
amidaniyorai
```
"""
function unidecode(str::AbstractString)::AbstractString
    buf = IOBuffer()
    ascii_char = Char(0)
    @inbounds for c in str
        code_point = codepoint(c)
        if code_point < 0x80
            ascii_char = c
            print(buf, ascii_char)
            continue
        elseif code_point > 0xffff
            continue
        else
            # Chop off last two hex digits
            section = code_point >> 8
            # Last two hex digits
            pos = code_point % 256
            cache = get_cache(section)
            print(buf, cache[pos + 1])
        end
    end
    # rstrip if from unicode table
    if buf.size > 1 && buf.data[buf.size] == 0x20 && ascii_char != ' '
        buf.size -= 1
    end
    String(take!(buf))
end

function get_cache(section::UInt32)::Vector{String}
    if isassigned(table, section+1)
        table[section + 1]
    else
        path = joinpath(@__DIR__, "..", "resources", "X$(string(section, base = 16, pad = 3))" )
        if isfile(path)
            table[section + 1] = readlines(path)
        else
            # No match, remove all
            table[section + 1] = fill("", 256)
        end
    end
end

export unidecode
end