text:
using DynamicTimeWarp
using Base.Test

#############################################
# Test the Distance functions
#############################################

Distance = DynamicTimeWarp.Distance
@test Distance.square(0,0) == 0
@test Distance.square(0,1) == 1
@test Distance.square(9,10) == 1
@test Distance.square(0,5) == 25
@test Distance.square(0.0,5.0) == 25.0
@test Distance.square(0.0,5) == 25.0

@test Distance.absval(0,0) == 0
@test Distance.absval(0,3) == 3
@test Distance.absval(0,5) == 5
@test Distance.absval(-2,1) == 3
@test Distance.absval(-2,3) == 5
@test Distance.absval(-2,5) == 7
@test Distance.absval(-2,-5) == 3

s1, s2 = collect(1:5), collect(2:2:10)
pp1 = Distance.poissonpenalty(s1, s2)
pp2 = Distance.poissonpenalty(15, 30)
@test pp1(0,0) == 0
@test pp2(0,0) == 0
@test pp1(1,2) == 0
@test pp1(5,10) == 0
@test pp2(5,10) == 0
@test pp1(9,0) == 8.0
@test pp2(9,0) == 8.0
@test pp1(0,9) == 2.0
@test pp2(0,9) == 2.0
@test pp1(9,9) == 1.0
@test pp2(9,9) == 1.0

#############################################
# Test dtw itself
#############################################

a=[1,1,1,2,4,6,5,5,5,4,4,3,1,1,1]
b=[1,1,2,4,6,6,6,5,4,4,4,3,3,3,1]
cost, match1, match2 = dtw(a,b)
@test cost==0
@test match1==[1,2,3,4,5,6,6,6,7,8,9,10,10,11,12,12,12,13,14,15]
@test match2==[1,1,2,3,4,5,6,7,8,8,8, 9,10,11,12,13,14,15,15,15]

a[end] += 2
cost, match1, match2 = dtw(a,b)
@test cost==4

a=collect(1:10)
b=a+1
cost, match1, match2 = dtw(a,b)
@test cost==2

a=zeros(Int,6)
b=1+a
cost, match1, match2 = dtw(a,b)
@test cost==length(a)

# Verify that a tie prefers diagonal moves
a=[1,1,1]
b=[1,1,1]
cost, pa, pb = dtw(a,b)
@test cost==0
@test pa==[1,2,3]
@test pb==[1,2,3]

# Verify that trackback ends properly if it reaches an edge before reaching [1,1]
# Also check that trackback prefers diagonal moves
a=[0,1,1,1]
b=[0,0,1,1]
cost, pa, pb = dtw(a,b)
@test cost==0
@test pa==[1,1,2,3,4]
@test pb==[1,2,3,3,4]

#############################################
# Test DTW with windows
#############################################

# Verify that a tie prefers diagonal moves
a=[1,1,1]
b=[1,1,1]
cost, pa, pb = dtwwindowed(a,b,[1,1,1],[3,3,3])
@test cost==0
@test pa==[1,2,3]
@test pb==[1,2,3]

# Verify that trackback ends properly if it reaches an edge before reaching [1,1]
# Also check that trackback prefers diagonal moves
a=[0,1,1,1]
b=[0,0,1,1]
cost, pa, pb = dtwwindowed(a,b,[1,1,1,1],[4,4,4,4])
@test cost==0
@test pa==[1,1,2,3,4]
@test pb==[1,2,3,3,4]

# First do the windowed test w/o windows
a=[0,1,2,3,4,4,4,4]
b=[0,0,1,2,2,2,3,4]
best_pa = [1,1,2,3,3,3,4,5,6,7,8]
best_pb = [1,2,3,4,5,6,7,8,8,8,8]
cost, pa, pb = dtw(a,b)
@test cost == 0
@test pa == best_pa
@test pb == best_pb

# Wide window, not touching optimal path
rmin = [1,1,1,2,3,4,5,6]
rmax = [4,6,7,8,8,8,8,8]
cost, pa, pb = dtwwindowed(a,b,rmin,rmax)
@test cost == 0
@test pa == best_pa
@test pb == best_pb

# Bottom of window is optimal path
rmin = [1,3,4,7,8,8,8,8]
rmax = [4,6,7,8,8,8,8,8]
cost, pa, pb = dtwwindowed(a,b,rmin,rmax)
@test cost == 0
@test pa == best_pa
@test pb == best_pb

# Top of window is optimal path
rmin = [1,1,1,2,3,4,5,6]
rmax = [2,3,6,7,8,8,8,8]
cost, pa, pb = dtwwindowed(a,b,rmin,rmax)
@test cost == 0
@test pa == best_pa
@test pb == best_pb

# Top and bottom of window are optimal path
rmin = [1,3,4,7,8,8,8,8]
rmax = [2,3,6,7,8,8,8,8]
cost, pa, pb = dtwwindowed(a,b,rmin,rmax)
@test cost == 0
@test pa == best_pa
@test pb == best_pb

# Now top of window cuts into optimal path
rmin = [1,1,1,2,3,4,5,6]
rmax = [4,4,5,6,7,8,8,8]
cost, pa, pb = dtwwindowed(a,b,rmin,rmax)
@test cost == 2
@test pa == [1,1,2,3,3,4,5,6,7,8]
@test pb == [1,2,3,4,5,6,7,8,8,8]

# Compare windowed and regular on non-square geometry:
seq1=collect(1:16)
seq2=seq1[1:2:end]
cost,pa,pb = dtw(seq1,seq2)
n1,n2 = length(seq1), length(seq2)
cost2,qa,qb = dtwwindowed(seq1, seq2, fill(1,n1), fill(n2,n1))
@test cost==cost2
@test pa==qa
@test pb==qb

seq1=[1.2051,1.7890,2.1314,2.2423,2.2101,2.0234,1.6898,1.3422,1.1530,1.2197,1.5435,2.07147,2.7359,3.469,4.2033,4.8704,5.4056,5.75243,5.8693,5.7362,5.35974,4.7761,4.0482,3.2548,2.4764]
seq2=[1.4927,2.1859,2.1182,1.5186, 1.1858,1.8034,3.0968,4.5317,5.5763,5.8038,5.0724,3.6576,2.4823]

cost,pa,pb = dtw(seq1,seq2)
n1,n2 = length(seq1), length(seq2)
cost2,qa,qb = dtwwindowed(seq1, seq2, fill(1,n1), fill(n2,n1))
@test cost==cost2
@test pa==qa
@test pb==qb

#############################################
# Test sequence compressions
#############################################

compress = DynamicTimeWarp.compress
s=collect(0:2:98)
s1 = compress(s)
s2 = compress(s1)
@test s1==float(collect(1:4:97))
@test s2==float(vcat(collect(3:8:91),[97]))

s = [1]
s1 = compress(s)
@test s1==[1.0]

#############################################
# Test window computation
#############################################
computewindow = DynamicTimeWarp.computewindow

# Simplest path (along the diagonal)
p=collect(1:8)
rmin,rmax = computewindow(p,p,1)
@test rmin==[1,1,1,2,3,4,5,6]
@test rmax==[3,4,5,6,7,8,8,8]

rmin,rmax = computewindow(p,p,2)
@test rmin==[1,1,1,1,1,2,3,4]
@test rmax==[5,6,7,8,8,8,8,8]

# A warpy path
pa=[1,1,2,3,4,5,6,7,8,8,8]
pb=[1,2,3,3,3,4,4,5,6,7,8]
rmin,rmax = computewindow(pa, pb, 1)
@test pa[end]==length(rmin)
@test pa[end]==length(rmax)
@test rmin==[1,1,2,2,2,3,3,4]
@test rmax==[4,4,4,5,5,6,8,8]

rmin,rmax = computewindow(pa, pb, 2)
@test pa[end]==length(rmin)
@test pa[end]==length(rmax)
@test rmin==[1,1,1,1,1,1,2,2]
@test rmax==[5,5,6,6,7,8,8,8]

rmin,rmax = computewindow(pa, pb, 20)
@test pa[end]==length(rmin)
@test pa[end]==length(rmax)
@test rmin==fill(1,8)
@test rmax==fill(8,8)

# Extreme path: follows left then upper edge
pa=[1,1,1,1,1,1,1,1,2,3,4,5,6,7,8]
pb=[1,2,3,4,5,6,7,8,8,8,8,8,8,8,8]
rmin,rmax = computewindow(pa, pb, 1)
@test pa[end]==length(rmin)
@test pa[end]==length(rmax)
@test rmin==[1,1,7,7,7,7,7,7]
@test rmax==[8,8,8,8,8,8,8,8]

rmin,rmax = computewindow(pa, pb, 2)
@test pa[end]==length(rmin)
@test pa[end]==length(rmax)
@test rmin==[1,1,1,6,6,6,6,6]
@test rmax==[8,8,8,8,8,8,8,8]

# More columns than rows
pa=[1,2,3,4,5,6,7,8]
pb=[1,2,3,4,4,4,4,4]
rmin,rmax = computewindow(pa, pb, 1)
@test rmin==[1,1,1,2,3,3,3,3]
@test rmax==[3,4,4,4,4,4,4,4]

rmin,rmax = computewindow(pa, pb, 2)
@test rmin==[1,1,1,1,1,2,2,2]
@test rmax==fill(4,8)

rmin,rmax = computewindow(pa, pb, 3)
@test rmin==fill(1,8)
@test rmax==fill(4,8)

rmin,rmax = computewindow(pa, pb, 4)
@test rmin==fill(1,8)
@test rmax==fill(4,8)

rmin,rmax = computewindow(pa, pb, 47)
@test rmin==fill(1,8)
@test rmax==fill(4,8)

# More rows than columns
pa,pb = pb,pa
rmin,rmax = computewindow(pa, pb, 1)
@test rmin==[1,1,1,2]
@test rmax==[3,4,8,8]

rmin,rmax = computewindow(pa, pb, 2)
@test rmin==fill(1,4)
@test rmax==[5,8,8,8]

rmin,rmax = computewindow(pa, pb, 3)
@test rmin==fill(1,4)
@test rmax==fill(8,4)

rmin,rmax = computewindow(pa, pb, 4)
@test rmin==fill(1,4)
@test rmax==fill(8,4)

rmin,rmax = computewindow(pa, pb, 47)
@test rmin==fill(1,4)
@test rmax==fill(8,4)

#############################################
# Test DTW and FastDTW on simple, large cases
#############################################
t=collect(1:1600)
pktimes=[100,300,1000,1300]
x = 1*exp(-0.5*((t-pktimes[1])/100).^2);
x+= 2*exp(-0.5*((t-pktimes[2])/150).^2);
x+= 3*exp(-0.5*((t-pktimes[3])/250).^2);
x+= 4*exp(-0.5*((t-pktimes[4])/250).^2);
y = x[1:2:end];
cost,px,py = dtw(x,y)
cost1,qx,qy = DynamicTimeWarp.fastdtw(x, y, 15)

## Need to actually TEST something here!


text:
# This example came from the R package ggm.

using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir) #do

derived_df = DataFrame!(CSV.File(scm_path("..", "data", "derived.csv")));

df = DataFrame(
  :y => log.(derived_df.Sys ./ derived_df.Dia),
  :x => log.(derived_df.Dia),
  :z => derived_df.Wei ./ (derived_df.Hei / 100).^2,
  :w => derived_df.Age
  )

R_cov = "

           Sys       Dia        Age       Hei       Wei
Sys 176.691332 96.300211  71.738901 -2.415433  27.49471
Dia  96.300211 92.758985  54.571882 -4.371036  33.11839
Age  71.738901 54.571882 109.697146 -8.497357  60.60994
Hei  -2.415433 -4.371036  -8.497357 32.335624  24.83404
Wei  27.494715 33.118393  60.609937 24.834038 113.66808
";

# A DAG model with a latent variable U
dag = DAG("ggm_derived", "DAG(y ~ z + u, x ~ u + w, z ~ w)")

to_dagitty(dag.d) |> display

fname = ProjDir * "/derived.dot"
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

set_dag_df!(dag, df; force=true)

println("Definition of DAG:")
display(dag)

bs = basis_set(dag; debug=true)
display(bs)

#t = shipley_test(dag)
#display("shipley_test = $t"); println()

pt = pcor_test(dag, [:w, :y], 1, 44)
display("pcor_test = $pt"); println()

as = adjustment_sets(dag, :w, :z)
as |> display

println("\nAncestral graph:")
ag = ancestral_graph(dag; m=[:u])
display(ag)

ext: jl


text:
# An implementation of a FISTA-like method, where the smooth part of the objective function can be strongly convex.

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

"""
    SFISTAIteration(; <keyword-arguments>)

Instantiate the FISTA-like method in [3] for solving strongly-convex composite optimization problems of the form

    minimize f(x) + h(x),

where h is proper closed convex and f is a continuously differentiable function that is μ-strongly convex and whose gradient is
Lf-Lipschitz continuous.

The scheme is based on Nesterov's accelerated gradient method [1, Eq. (4.9)] and Beck's method for the convex case [2]. Its full
definition is given in [3, Algorithm 2.2.2.], and some analyses of this method are given in [3, 4, 5]. Another perspective is that
it is a special instance of [4, Algorithm 1] in which μh=0.

# Arguments
- `y0`: initial point; must be in the domain of h.
- `f=Zero()`: smooth objective term.
- `h=Zero()`: proximable objective term.
- `μf=0` : strong convexity constant of f (see above).
- `Lf` : Lipschitz constant of ∇f (see above).
- `adaptive=false` : enables the use of adaptive stepsize selection.

# References
1. Nesterov, Y. (2013). Gradient methods for minimizing composite functions. Mathematical Programming, 140(1), 125-161.
2. Beck, A., & Teboulle, M. (2009). A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1), 183-202.
3. Kong, W. (2021). Accelerated Inexact First-Order Methods for Solving Nonconvex Composite Optimization Problems. arXiv preprint arXiv:2104.09685.
4. Kong, W., Melo, J. G., & Monteiro, R. D. (2021). FISTA and Extensions - Review and New Insights. arXiv preprint arXiv:2107.01267.
5. Florea, M. I. (2018). Constructing Accelerated Algorithms for Large-scale Optimization-Framework, Algorithms, and Applications.
"""
Base.@kwdef struct SFISTAIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},Tf,Th}
    y0::Tx
    f::Tf = Zero()
    h::Th = Zero()
    Lf::R
    μf::R = real(eltype(Lf))(0.0)
end

Base.IteratorSize(::Type{<:SFISTAIteration}) = Base.IsInfinite()

Base.@kwdef mutable struct SFISTAState{R, Tx}
    λ::R                                # stepsize (mutable if iter.adaptive == true).
    yPrev::Tx                           # previous main iterate.
    y:: Tx = zero(yPrev)                # main iterate.
    xPrev::Tx = copy(yPrev)             # previous auxiliary iterate.
    x::Tx  = zero(yPrev)                # auxiliary iterate (see [3]).
    xt::Tx = zero(yPrev)                # prox center used to generate main iterate y.
    τ::R = real(eltype(yPrev))(1.0)     # helper variable (see [3]).
    a::R = real(eltype(yPrev))(0.0)     # helper variable (see [3]).
    APrev::R = real(eltype(yPrev))(1.0) # previous A (helper variable).
    A::R = real(eltype(yPrev))(0.0)     # helper variable (see [3]).
    gradf_xt::Tx = zero(yPrev)          # array containing ∇f(xt).
end

function Base.iterate(
    iter::SFISTAIteration,
    state::SFISTAState = SFISTAState(λ=1/iter.Lf, yPrev=copy(iter.y0))
)
    # Set up helper variables.
    state.τ = state.λ * (1 + iter.μf * state.APrev)
    state.a = (state.τ + sqrt(state.τ ^ 2 + 4 * state.τ * state.APrev)) / 2
    state.A = state.APrev + state.a
    state.xt .= (state.APrev / state.A) .* state.yPrev + (state.a / state.A) .* state.xPrev
    gradient!(state.gradf_xt, iter.f, state.xt)
    λ2 = state.λ / (1 + state.λ * iter.μf)
    # FISTA acceleration steps.
    prox!(state.y, iter.h, state.xt - λ2 * state.gradf_xt, λ2)
    state.x .= state.xPrev .+ (state.a / (1 + state.A * iter.μf)) .* ((state.y .- state.xt) ./ state.λ .+ iter.μf .* (state.y .- state.xPrev))
    # Update state variables.
    state.yPrev .= state.y
    state.xPrev .= state.x
    state.APrev = state.A
    return state, state
end

## Solver.

struct SFISTA{R, K}
    maxit::Int
    tol::R
    termination_type::String
    verbose::Bool
    freq::Int
    kwargs::K
end

# Different stopping conditions (sc). Returns the current residual value and whether or not a stopping condition holds.
function check_sc(state::SFISTAState, iter::SFISTAIteration, tol, termination_type)
    if termination_type == "AIPP"
        # AIPP-style termination [4]. The main inclusion is: r ∈ ∂_η(f + h)(y).
        r = (iter.y0 - state.x) / state.A
        η = (norm(iter.y0 - state.y) ^ 2 - norm(state.x - state.y) ^ 2) / (2 * state.A)
        res = (norm(r) ^ 2 + max(η, 0.0)) / max(norm(iter.y0 - state.y + r) ^ 2, 1e-16)
    else
        # Classic (approximate) first-order stationary point [4]. The main inclusion is: r ∈ ∇f(y) + ∂h(y).
        λ2 = state.λ / (1 + state.λ * iter.μf)
        gradf_y, = gradient(iter.f, state.y)
        r = gradf_y - state.gradf_xt + (state.xt - state.y) / λ2
        res = norm(r)
    end
    return res, (res <= tol || res ≈ tol)
end

# Functor ('function-like object') for the above type.
function (solver::SFISTA)(y0; kwargs...)
    raw_iter = SFISTAIteration(; y0=y0, solver.kwargs..., kwargs...)
    stop(state::SFISTAState) = check_sc(state, raw_iter, solver.tol, solver.termination_type)[2]
    disp((it, state)) = @printf("%5d | %.3e\n", it, check_sc(state, raw_iter, solver.tol, solver.termination_type)[1])
    iter = take(halt(raw_iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.y, num_iters
end

SFISTA(; maxit=10_000, tol=1e-6, termination_type="", verbose=false, freq=(maxit < Inf ? Int(maxit/100) : 100), kwargs...) =
    SFISTA(maxit, tol, termination_type, verbose, freq, kwargs)


text:
using Conda, Compat, VersionParsing
using Compat
using Compat: @info
using Compat.Test

exe = Compat.Sys.iswindows() ? ".exe" : ""

Conda.update()

env = :test_conda_jl
rm(Conda.prefix(env); force=true, recursive=true)

@test Conda.exists("curl", env)
Conda.add("curl", env)

@testset "Install Python package" begin
    Conda.add("python=3.6", env)  # 3.7 doesn't work on Windows at the moment
    pythonpath = joinpath(Conda.python_dir(env), "python" * exe)
    @test isfile(pythonpath)

    cmd = Conda._set_conda_env(`$pythonpath -c "import zmq"`, env)
    @test_throws Exception run(cmd)
    Conda.add("pyzmq", env)
    run(cmd)
end

curlvers = Conda.version("curl",env)
@test curlvers >= v"5.0"
@test Conda.exists("curl==$curlvers", env)

curl_path = joinpath(Conda.bin_dir(env), "curl" * exe)
@test isfile(curl_path)

@test "curl" in Conda.search("cu*", env)

Conda.rm("curl", env)
@test !isfile(curl_path)

pythonpath = joinpath(Conda.PYTHONDIR, "python" * exe)
@test isfile(pythonpath)
pyversion = read(`$pythonpath -c "import sys; print(sys.version)"`, String)
@test pyversion[1:1] == Conda.MINICONDA_VERSION

Conda.add_channel("foo", env)
@test Conda.channels(env) == ["foo", "defaults"]
# Testing that calling the function twice do not fail
Conda.add_channel("foo", env)

Conda.rm_channel("foo", env)
@test Conda.channels(env) == ["defaults"]

@testset "Batch install and uninstall" begin
    Conda.add(["affine", "ansi2html"], env)
    installed = Conda._installed_packages(env)
    @test "affine" ∈ installed
    @test "ansi2html" ∈ installed

    Conda.rm(["affine", "ansi2html"], env)
    installed = Conda._installed_packages(env)
    @test "affine" ∉ installed
    @test "ansi2html" ∉ installed
end


text:
module TestFromDict

using Test
using Configurations
using Configurations: from_dict_union_type_dynamic

@option "option_a" struct OptionAliasA
    x::Int = 1
    y::Int = 2
end

@option struct IgnoreExtra
    x::Int = 1
    y::Int = 2
end

Configurations.ignore_extra(::Type{IgnoreExtra}) = true

@option "option_b" struct OptionAliasB
    x::Int = 1
end

@option "option_c" struct OptionAliasC
    y::Int = 1
end

@option struct OptionReflectA
    type::Reflect
    x::Int = 1
    y::Int = 2
end

@option struct OptionReflectB
    type::Reflect
    x::Int = 1
end

@option struct OptionReflectC
    type::Reflect
    y::Int = 1
end

@option "option_type" struct ReflectWithAlias
    type::Reflect
    y::Int = 1
end

@option struct TestUnionType
    maybe_reflect::Maybe{OptionReflectA}
    maybe_alias::Maybe{OptionAliasA}
    union_of_alias::Union{OptionAliasA,OptionAliasB,OptionAliasC}
    union_of_reflects::Union{OptionReflectA,OptionReflectB,OptionReflectC}
    maybe_of_alias::Maybe{Union{OptionAliasA,OptionAliasB,OptionAliasC}}
    maybe_of_reflects::Maybe{Union{OptionReflectA,OptionReflectB,OptionReflectC}}
    totally_mixed::Union{OptionAliasA,OptionAliasB,OptionReflectA,OptionReflectB}
    maybe_totally_mixed::Maybe{
        Union{OptionAliasA,OptionAliasB,OptionReflectA,OptionReflectB}
    }
end

@testset "from_dict" begin
    @testset "ignore_extra" begin
        d = Dict{String, Any}(
            "x" => 1,
            "y" => 2,
            "z" => 3,
        )
        @test_throws InvalidKeyError from_dict(OptionAliasA, d)
        @test from_dict(IgnoreExtra, d) == IgnoreExtra()
    end

    @testset "test maybe union type" begin
        d = Dict{String,Any}(
            "maybe_reflect" => Dict{String,Any}(),
            "maybe_alias" => Dict{String,Any}(),
            "union_of_alias" => Dict{String,Any}("option_b" => Dict{String,Any}("x" => 1)),
            "union_of_reflects" =>
                Dict{String,Any}("type" => "Main.TestFromDict.OptionReflectC", "y" => 1),
            "maybe_of_alias" => nothing,
            "maybe_of_reflects" => nothing,
            "totally_mixed" =>
                Dict{String,Any}("type" => "Main.TestFromDict.OptionReflectB", "x" => 1),
            "maybe_totally_mixed" =>
                Dict{String,Any}("option_a" => Dict{String,Any}("x" => 1, "y" => 2)),
        )
        option = TestUnionType(;
            union_of_alias=OptionAliasB(),
            union_of_reflects=OptionReflectC(),
            totally_mixed=OptionReflectB(),
            maybe_totally_mixed=OptionAliasA(),
        )
        @test from_dict(TestUnionType, d) == option
    end

    @testset "test dynamic union type" for option in [
        TestUnionType(;
            maybe_reflect=OptionReflectA(),
            maybe_alias=OptionAliasA(),
            union_of_alias=OptionAliasB(),
            union_of_reflects=OptionReflectC(),
            maybe_of_alias=OptionAliasB(),
            maybe_of_reflects=OptionReflectC(),
            totally_mixed=OptionReflectB(),
            maybe_totally_mixed=OptionAliasA(),
        ),
        TestUnionType(;
            maybe_reflect=OptionReflectA(),
            maybe_alias=OptionAliasA(),
            union_of_alias=OptionAliasA(),
            union_of_reflects=OptionReflectB(),
            maybe_of_alias=OptionAliasA(),
            maybe_of_reflects=OptionReflectB(),
            totally_mixed=OptionAliasA(),
            maybe_totally_mixed=OptionReflectB(),
        ),
    ]
        d = to_dict(option)
        @test from_dict(TestUnionType, d) == option

        for idx in 1:fieldcount(TestUnionType)
            field = fieldname(TestUnionType, idx)
            type = fieldtype(TestUnionType, idx)
            types = Base.uniontypes(type)
            if Nothing in types
                type = Union{filter(x -> x !== Nothing, types)...}
            end
            @test from_dict_union_type_dynamic(
                TestUnionType, OptionField(field), type, d[string(field)]
            ) == getfield(option, field)
        end
    end
end

end # TestFromDict


text:
function test7()
  filename = "test7.asc"
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)

  @test sim["A"] == 10e-12
  @test ~haskey(sim,"b") # keys only for numbers
  @test ~haskey(sim,"C") # keys only for numbers
  @test ~haskey(sim,"d") # keys only for numbers
  @test ~haskey(sim,"a") # keys are case sensitive
  @test ~haskey(sim,"B") # both
  @test ~haskey(sim,"D") # both

  @test sim["Ee"] == 1e3
  @test sim["fF"] == 1e3
  @test sim["GG"] == 1e6
  @test sim["h"] == 1.123e6
  @test sim["i"] == 1e9
  @test sim["j"] == 1e9
  @test sim["k"] == 1e12
  @test sim["l"] == 1e12
  @test sim["m"] == 1e-3
  @test sim["n"] == 1e-3
  @test sim["o"] == 1e-6
  @test sim["p"] == 1e-6
  @test sim["q"] == 1e-9
  @test sim["r"] == 1e-9
  @test sim["s"] == 1e-12
  @test sim["t"] == 1e-12
  @test sim["u"] == 1e-15
  @test sim["v"] == 1e-15
  show(IOBuffer(),sim)
end
test7()


text:
# This file is part of TaylorSeries.jl, MIT licensed
#
# Tests for TaylorSeries

testfiles = (
    "onevariable.jl",
    "manyvariables.jl",
    "mixtures.jl",
    "mutatingfuncts.jl",
    "identities_Euler.jl",
    "fateman40.jl",
    "intervals.jl"
    )

for file in testfiles
    include(file)
end


text:
"""
    add_variables_from_data!(gm::Union{JuMP.Model, GlobalModel},
                            X::DataFrame)

Adds/finds variables depending on the columns of X.
"""
function add_variables_from_data!(m::JuMP.Model,
                                 X::DataFrame)
    varkeys = names(X)
    vars = []
    data_vars = []
    for key in varkeys
        try
            push!(vars, m[Symbol(key)])
        catch KeyError
            nv = @variable(m, base_name = string(key))
            m[Symbol(key)] = nv
            push!(vars, nv)
            push!(data_vars, Symbol(key))
        end
    end
    length(data_vars) != 0 && @info("Added $(length(data_vars)) variables from data: " * string(data_vars))
    return vars
end

add_variables_from_data!(gm::GlobalModel, X::DataFrame) = add_variables_from_data!(gm.model, X)

"""
    bound_to_data!(gm::Union{JuMP.Model, GlobalModel},
                  X::DataFrame)

Constrains the domain of relevant variables to the box interval defined by X.
"""
function bound_to_data!(m::JuMP.Model,
                       X::DataFrame)
    vars = []
    try
        vars = [m[Symbol(key)] for key in names(X)]
    catch UndefVarError
        throw(OCTHaGOnException("Please make sure the keys of X match with variables in model $(gm)."))
    end
    ranges = Dict(m[Symbol(key)] => [minimum((X[!, Symbol(key)])), maximum((X[!, Symbol(key)]))] for key in names(X))
    bound!(m, ranges);
    return
end

bound_to_data!(gm::GlobalModel, X::DataFrame) = bound_to_data!(gm.model, X)

"""
    add_datadriven_constraint(gm::GlobalModel,
                     X::DataFrame, Y::Array;
                     constraint::Union{Nothing, JuMP.ConstraintRef, Expr} = nothing,
                     vars::Union{Nothing, Array{JuMP.VariableRef, 1}} = nothing,
                     dependent_var::Union{Nothing, JuMP.VariableRef} = nothing,
                     name::String = "bbl" * string(length(gm.bbls) + 1),
                     equality::Bool = false)

 Adds a data-driven constraint to GlobalModel. Data driven BBLs
do not allow for resampling.
"""
function add_datadriven_constraint(gm::GlobalModel,
                     X::DataFrame, Y::Array;
                     constraint::Union{Nothing, JuMP.ConstraintRef, Expr} = nothing,
                     vars::Union{Nothing, Array{JuMP.VariableRef, 1}} = nothing,
                     dependent_var::Union{Nothing, JuMP.VariableRef} = nothing,
                     name::String = "bbl" * string(length(gm.bbls) + 1),
                     equality::Bool = false)
    if isnothing(vars)
        vars = add_variables_from_data!(gm, X)
    else
        @assert size(X, 2) == length(vars)
    end
    expr_vars = vars
    if isnothing(dependent_var)
        new_bbl = BlackBoxClassifier(constraint = constraint, vars = vars, expr_vars = expr_vars,
                                        equality = equality, name = name)
        add_data!(new_bbl, X, Y)
        push!(gm.bbls, new_bbl)
        return
    else
        new_bbl = BlackBoxRegressor(constraint = constraint, vars = vars, expr_vars = expr_vars,
                                    dependent_var = dependent_var, equality = equality, name = name)
        add_data!(new_bbl, X, Y)
        push!(gm.bbls, new_bbl)
        return
    end
end


text:
using consonance

push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")
push!(LOAD_PATH, "/pool/scout0/DiskJockey/")

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--run_index", "-r"
    help = "Output run index"
    arg_type = Int
end

parsed_args = parse_args(ARGS, s)
run_index = parsed_args["run_index"]

using constants
using image

import PyPlot.plt
using LaTeXStrings

global norm = plt.Normalize(vmin=log10(vmax_d) - 10, vmax=log10(vmax_d)-2, clip=false)

# Plot the channel maps using sky convention
function plot_dust(img::image.SkyImage, id::Int)

    fig, ax = plt.subplots(ncols=1, figsize=(2.0, 2.0))
    # Image needs to be flipped along RA dimension
    ext = (img.ra[end], img.ra[1], img.dec[1], img.dec[end])

    ax[:set_xlabel](L"$\Delta \alpha$ ('')", size=8)
    ax[:set_ylabel](L"$\Delta \delta$ ('')", size=8)
    ax[:tick_params](axis="both", which="major", labelsize=8)

    #Flip the frame for Sky convention
    frame = fliplr(img.data[:,:,1])
    frame += 1e-99 #Add a tiny bit so that we don't have log10(0)
    lframe = log10(frame)
    vax = maximum(lframe)

    ax[:imshow](lframe, extent=ext, interpolation="none", origin="lower", norm=norm, cmap=plt.get_cmap("Greys"))

    # Annotate the parameters
    p = pars[id]
    mass = p.M_star
    r_c = p.r_c
    incl = p.incl
    label = L"$M_\ast$: " * @sprintf("%.2f", mass) * L" $M_\odot$   $r_c$: " * @sprintf("%.0f", r_c) * L" AU   $i$: " * @sprintf("%.0f", incl) * L"${}^\circ$"
    fig[:text](0.1, 0.07, label, size=8)

    fig[:subplots_adjust](top=0.95, bottom=0.26, left=0.25, right=0.75)

    plt.savefig(outdir * @sprintf("d%04d.png", id), dpi=300)
    plt.close("all")

end

# How many frames per process?
start = run_index * nframes_per_proc + 1

ids = Int[i for i=start:(start + nframes_per_proc)]

for id in ids
    if id <= nframes
        fname = outdir * @sprintf("dimage%04d.out", id)
        im = imread(fname)
        skim = imToSky(im, 73.)
        plot_dust(skim, id)
    else
        break
    end
end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
"""
ITensors is a library for rapidly creating correct and efficient tensor network algorithms.

An ITensor is a tensor whose interface is independent of its memory layout. ITensor indices are objects which carry extra information and which 'recognize' each other (compare equal to each other).

The ITensor library also includes composable and extensible algorithms for optimizing and transforming tensor networks, such as matrix product state and matrix product operators, such as the DMRG algorithm.
"""
module ITensors

#####################################
# External packages
#
using Compat
using HDF5
using KrylovKit
using LinearAlgebra
using NDTensors
using PackageCompiler
using Pkg
using Printf
using Random
using StaticArrays
using TimerOutputs

#####################################
# NDTensors (definitions that will be moved to NDTensors
# module)
#
include("NDTensors/NDTensors.jl")

#####################################
# ContractionSequenceOptimization
#
include("ContractionSequenceOptimization/ContractionSequenceOptimization.jl")
using .ContractionSequenceOptimization

#####################################
# Directory helper functions (useful for
# running examples)
#
src_dir() = dirname(pathof(@__MODULE__))
pkg_dir() = joinpath(src_dir(), "..")
examples_dir() = joinpath(pkg_dir(), "examples")

#####################################
# Determine version and uuid of the package
#
function _parse_project_toml(field::String)
  return Pkg.TOML.parsefile(joinpath(pkg_dir(), "Project.toml"))[field]
end
version() = VersionNumber(_parse_project_toml("version"))
uuid() = Base.UUID(_parse_project_toml("uuid"))

#####################################
# Exports
#
include("exports.jl")

#####################################
# Imports
#
include("imports.jl")

#####################################
# Global Variables
#
include("global_variables.jl")

#####################################
# Index and IndexSet
#
include("smallstring.jl") # Not currently using in TagSet
include("readwrite.jl")
include("not.jl")
include("tagset.jl")
include("arrow.jl")
include("symmetrystyle.jl")
include("index.jl")
include("set_operations.jl")
include("indexset.jl")

#####################################
# ITensor
#
include("itensor.jl")
include("broadcast.jl")
include("decomp.jl")
include("iterativesolvers.jl")

#####################################
# QNs
#
include("qn/qn.jl")
include("qn/qnindex.jl")
include("qn/qnindexset.jl")
include("qn/qnitensor.jl")

#####################################
# MPS/MPO
#
include("mps/abstractmps.jl")
include("mps/deprecated.jl")
include("mps/mps.jl")
include("mps/mpo.jl")
include("mps/sweeps.jl")
include("mps/projmpo.jl")
include("mps/projmposum.jl")
include("mps/projmps.jl")
include("mps/projmpo_mps.jl")
include("mps/observer.jl")
include("mps/dmrg.jl")

#####################################
# Physics
#
include("physics/sitetype.jl")
include("physics/lattices.jl")
include("physics/site_types/generic_sites.jl")
include("physics/site_types/qubit.jl")
include("physics/site_types/spinhalf.jl")
include("physics/site_types/spinone.jl")
include("physics/site_types/fermion.jl")
include("physics/site_types/electron.jl")
include("physics/site_types/tj.jl")
include("physics/fermions.jl")
include("physics/autompo.jl")

#####################################
# Deprecations
#
include("deprecated.jl")

#####################################
# Argument parsing
#
include("argsdict/argsdict.jl")

#####################################
# Package compilation
#
include("packagecompile/compile.jl")

#####################################
# Developer tools, for internal
# use only
#
include("developer_tools.jl")

function __init__()
  return resize!(empty!(INDEX_ID_RNGs), Threads.nthreads()) # ensures that we didn't save a bad object
end

#####################################
# Precompile certain functions
# (generated from precompile/make_precompile.jl
# using SnoopCompile.jl)
#
#include("../precompile/precompile.jl")
#_precompile_()

end # module ITensors


text:
"""
Authors: Yaroslav Salii, 2020+
         Kara Ignatenko, 2021

Submodule responsible for handling input and output, i.e.
reading tract and travel data from disk and outputting travel matrices.

data-io.jl
2021-08-13 v.0.1: First modular version
2021-08-24 v.0.2: Made the I/O directory path independent of the working dir;
                  added select_mkid to export
2021-08-24 v.0.3: renamed IO to DataIO to avoid confusion with any builtin
2021-08-24 v.0.4: added heuristic to determine project root
2021-09-08 v.0.5: add infect2! to try to infect tract 41~001~950100 or thereabouts
2021-10-13 v.0.6: automatically normalize single-digit FIPS codes without a leading zero

"""
module DataIO

export fn,
       writeMe,
       checkIfFilesExist,
       rdAPs,
       rdBTS,
       rdTidyWfsByFIPS,
       normalizeOneDigitFIPS,
       rdFluteTract,
       rdWholeUS,
       censorFluteTractByFIPS,
       ns2iv,
       ns2iv_sterile,
       infect2!,
       select_mkid

using DataFrames
#CSV: IO FluTE & my, DelimitedFiles: writing the matrices (.dat)
using CSV,DelimitedFiles

#find the project root by searching for `epi-net-m` in this file's path
#then obtain the absolute path of the `data` subdirectory
thisPath = splitpath(@__DIR__)
projRoot = thisPath[1:findfirst(isequal("epi-net-m"), thisPath)]
const dataDir = joinpath(projRoot..., "data")

#===TYPES AND NAMING CONVENTIONS===#
"""Data type for holding paths and naming conventions for I/O"""
struct NamingSpec #all fields are String, don't say I didn't warn you
    # e.g. sep="-", sSep="_"; $name_$size-$suff
    sep::String #to the right of last $sep is filename suffix, to the left is the instance name
    sSep::String
    # read from ifDir, write to ofDir
    ifDir::String
 #  ifDirAir::String
    ofDir::String
    # what's after instance's name in its Initial Values file name
    fltInitSuff::String
    myInitSuff::String
end

#the naming conventions I am going to use
#as well as input and output directories
const fn=NamingSpec("-","_"
    ,joinpath(dataDir,"by-tract","flute")
    ,joinpath(dataDir,"by-tract")
    ,"tracts.dat","init.csv")

#===OUTPUT===#
"""Return `true` if there are files whose names start with `iname` in `odfir`."""
function checkIfFilesExist(iname::String; ofdir=fn.ofDir)
    dir = readdir(ofdir)
    any(filename -> startswith(filename, iname), dir)
end

"""Save the initial values table `ivs` and the travel matrix `A` to `ofdir`."""
function writeMe(iname::String,ivs::DataFrame,A::Array{Float64,2};ofdir=fn.ofDir)
    ofivs= join([iname,nrow(ivs)],"_") * "-init.csv"
    oftrv= join([iname,nrow(ivs)],"_") * "-trav.dat"
    CSV.write(joinpath(ofdir,ofivs),ivs)
    writedlm(joinpath(ofdir,oftrv),A)
end

#===AIRPORTS===#
#locating BTS and OpenFlights input files
const APdir= joinpath(dataDir,"by-tract","air")::String
#raw BTS data, with separate per-carrier flights
const ifBTS=joinpath(APdir,"2019 BTS domestic.csv")::String
#raw OpenFlights AP data
const ifAPs=joinpath(APdir,"Openflights airports.dat")::String
#TODO: consider wgetting from the original https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat
const apAllColNames = [:ID,:Name,:City,:Country,:IATA_Code,:ICAO_Code,:LAT,:LNG,:Altitude,:Timezone,:Daylight_Savings,:TZ,:Type,:Source]
#reordered in the order of necessity; 3-letter IATA code as ID
const apRetainedColNames=[:IATA_Code,:LAT,:LNG,:Name,:City,:Country]

"""
Read the OpenFlights.org's airports.dat, and retain only the columns
that will be used; put :IATA_Code,:LAT,:LNG first, these are definitely useful.
"""
#CAVEAT: some missing values in :City, some weird values in :Country
#TODO: resolve this caveat
function rdAPs(ifName=ifAPs::String)
    out=CSV.File(ifName,header=false) |> DataFrame
    rename!(out, apAllColNames)
    #retain only the useful columns
    select!(out,apRetainedColNames)
end

#---READ---BTS--FLIGHTS-------------------------#
"""
Read the BTS file, and retain only :1 Passengers, :5 ORIGIN, and :7 DEST;
set the columns to [:ORG,:DST,:PSG] for uniformity.
"""
function rdBTS(ifName=ifBTS::String)
    rawBTS = select(CSV.read(ifName,DataFrame), :1,:5,:7)
    if map( x -> floor(x)==x, rawBTS[:,1]) |> all #if :PSG are Integer
        rawBTS.PASSENGERS=convert(Array{Int64,1},rawBTS.PASSENGERS)
    else #∃ non-integer passengers-per-year entry
        error("CAVEAT: *non-integer* PASSENGERS in Flights input.")
    end
    rename!(rawBTS,[:PSG,:ORG,:DST])
    return select(rawBTS,[:ORG,:DST,:PSG]) #
end

#===COMMUTER FLOW (*-wf-*.dat files)===#
"""List the `-wf-` filenames found in `ins.ifDir`, default to `fn`"""
function lsWf(ins::NamingSpec = fn)
    filter(s::String -> startswith(s,"usa-wf-"),readdir(ins.ifDir))
end

ls_fipsAll = map(s -> split(split(s,".")[1],"-")[3] |> string,lsWf())

"""
Read and tidy all "usa-wf-\$fips.dat" for \$fips in fipss; default to NW: Oregon + Washington.
If the tracts argument is provided, commutes referencing tracts that aren't in it
will be discarded.
"""
function rdTidyWfsByFIPS(fipss::Array{String,1}=["41","53"], tracts::DataFrame=DataFrame(),
     ins::NamingSpec=fn)
    #generate FluTE filename by State's FIPS
    wf_by_FIPS(fips::String) = "usa-wf-$fips.dat"

    wfs = map(wf_by_FIPS, fipss)
    wfs2 = [CSV.File(ipath,
        header = false,
        types =[String,String,String,String,String,String,Int64]) |>
    DataFrame for ipath in map(f -> joinpath(ins.ifDir,f), wfs)]

    to_keep(r) = r[4] ∈ fipss && r[3]≠r[6]
    #TIDY: retain only the commute between states in `fipss`; this chucks the REFLEXIVE commute
    wfs3 = map(df -> filter(to_keep,eachrow(df) ) |> DataFrame,wfs2)
    wfs4 = (length(wfs3) > 1) ? reduce(vcat,wfs3) : wfs3[1] #add them all together
    #combine the State,County,Tract triples into single columns
    insertcols!(wfs4, 1,
        :ORG => map((s,z,w)-> join([s,z,w],"~"),wfs4.Column1, wfs4.Column2,wfs4.Column3))
    insertcols!(wfs4, 2,
        :DST => map((s,z,w)-> join([s,z,w],"~"),wfs4.Column4, wfs4.Column5,wfs4.Column6))
    rename!(wfs4, :Column7 => :CMT) #these are daily commuters between :ORG and :DST
    select!(wfs4,:ORG,:DST,:CMT) #chuck the unnecessary
    #Chuck commute pairs that have tract IDs not present in tracts
    wfs5 = isempty(tracts) ? # check if tracts has been passed
        wfs4 :
        scrubWfs(tracts, wfs4) |> DataFrame
    sort!(wfs5,[:ORG,:DST])

    return wfs5
end

"""
Ensure that the FIPS code `fips` has the canonical form by prefixing it with a 0
if it is one digit long.
"""
function normalizeOneDigitFIPS(fips::String)
    if length(fips) == 1
        return "0" * fips
    end
    fips
end

"""
Return `wfs` excluding the commutes for tracts not present in `ns`.
`ns` MUST have [:Name]; `wfs` MUST have [:ORG,:DST,:CMT]
"""
function scrubWfs(ns::DataFrame,wfs::DataFrame)
    allCmtTracts = (wfs.ORG |> unique) ∪ (wfs.DST |> unique)
     #these tracts are not present in FluTE's `us-tracts.dat`
    notPresent = setdiff(allCmtTracts,ns.Name)
    retwfs = filter(r -> r.ORG ∉ notPresent && r.DST ∉ notPresent,eachrow(wfs))
    return retwfs
end

#===FluTE TRACTS (*-tracts.dat files)===#
# sample usage: Oboe.lsTracts()[4] |> Oboe.rdFluteTract |> Oboe.aggBySte
# a *node* (tract) must have the fields :Pop,:LAT,:LNG

"""List the FluTE's tract filenames found in `ins.ifDir`, default to `fn`"""
function lsTracts(ins::NamingSpec = fn)
    filter(s::String -> endswith(s,ins.fltInitSuff),readdir(ins.ifDir))
end

"""
Load a FluTE census tract info into a `DataFrame`, setting the types and colnames.
Default to NW tracts (Oregon + Washington)
"""
function rdFluteTract(ifName::String=lsTracts()[findfirst(s -> startswith(s,"a~NW"),lsTracts())],
            ins::NamingSpec=fn)
    idf=CSV.File(joinpath(ins.ifDir,ifName)
    ,header=false
    ,types=[String,String,String,Int64,Float64,Float64]) |> DataFrame
    rename!(idf, [:Ste,:Cty,:Tra,:Pop,:LAT,:LNG])
    idf.Name = map((s,z,w)-> join([s,z,w],"~"),idf.Ste,idf.Cty,idf.Tra)
    return idf
end

"""Read the FluTE's all-US-by-tract instance (`usa-tracts.dat`)"""
function rdWholeUS(ins::NamingSpec=fn)
    ifName = filter(s -> split(s,"-")[1]=="usa",lsTracts(ins) )[1]
    rdFluteTract(ifName,ins)
end

"""Given a `DataFrame` of tracts, remove those that are not from the given states."""
function censorFluteTractByFIPS(tracts::DataFrame, fips::Vector{String})
    filter(row -> row.Ste ∈ fips, tracts)
end

#===EPIDEMIC SIMULATION AND INITIAL VALUES===#
"""
input MUST have [:Pop,:Name,:LAT,:LNG]
input SHOULD have [:IATA_Code], if absent, inserted as "dummy"
Initialize a `DataFrame` for SIR with the columns
`[:id,:IATA_Code,:N_i,:S_i,:I_i,:R_i,:Name,:LAT,:LNG]`
"""
function ns2iv_sterile(ns::DataFrame)
    out = copy(ns)
    rename!(out, :Pop => :N_i)
    insertcols!(out,1,:id => map(select_mkid(names(out)), eachrow(out)))
    insertcols!(out,4, :S_i => out.N_i)
    insertcols!(out,5, :I_i => zeros(Int64,nrow(out)))
    insertcols!(out,6, :R_i => zeros(Int64,nrow(out)))
    if "IATA_Code" ∉ names(out)
        insertcols!(out,2,:IATA_Code => repeat(["dummy"],nrow(out)))
    end
    select(out,[:id,:IATA_Code,:N_i,:S_i,:I_i,:R_i,:Name,:LAT,:LNG])
end

#quick delegate for the most common use case
"""Initialize a `DataFrame` for SIR and infect the first node."""
ns2iv(ns::DataFrame) = ns2iv_sterile(ns) |> infect!

"""
If the first node is sterile (`_I_i[1]==0`), seed it with 1 infected (Hi, I'm _idempotent_!)
input MUST have `[:I_i,:S_i]`
"""
function infect!(ns::DataFrame)
    if ns.I_i[1] == 0
        ns[1,:S_i]-=1
        ns[1,:I_i]+=1
    end
    return ns
end

"""
Try to infect tract 41~001~950100 or what's nearest to it in view of aggregation:
if aggregated to airports, try ALW
if not, try state 41 or county 41001, or tract "41~001~950001"
"""
function infect2!(ns::DataFrame;aggType::String)
    tryTargetIx = findfirst(r -> r.id == dispatchTargetNodeID(aggType),eachrow(ns))
    targetIx = isnothing(tryTargetIx) ? 1 : tryTargetIx #infect node 1 if target not found
    if ns.I_i[targetIx] == 0
        ns.S_i[targetIx] -= 1
        ns.I_i[targetIx] += 1
    end
    return ns
end

#---AUX---INITIAL--VALUES--GENERATION-------------#

"""
Select an `id`-generating function
0. if it's a census tract, just write \$Ste~\$Cty~\$Tra
1. if it's a county or a state, write as Integer to match `us-10m.json`
2. if it's a Voronoi cell of an AP, write the AP's IATA_Code
"""
function select_mkid(nms::Array{String})
    if ["Ste","Cty","Tra"] ⊆ nms #node's a census tract
         return r-> join([r.Ste,r.Cty,r.Tra],"~")
    elseif ["Ste","Cty"] ⊆ nms #node's a county, must accomodate `us-10m.json`
         return r-> parse(Int,r.Ste*"000") + parse(Int,r.Cty)
    elseif  ["Ste"] ⊆ nms #node's a state, must accomodate `us-10m.json`
         return r-> parse(Int,r.Ste)
    elseif ["IATA_Code"] ⊆ nms #node is Voronoi cell around an AP
         return r -> r.IATA_Code
    else
        error("Nodes header not recognized. Can't generate an :id without FIPS or IATA_Code.")
    end
end

"Match target node ID by aggregation type"
function dispatchTargetNodeID(agg::String)
    if agg == "tra"
        return "41~001~950100"
    elseif agg == "cty"
        return 41001
    elseif agg == "ap"
        return "ALW"
    elseif agg == "ste"
        return 41
    else
        error("Unknown aggregation requested. Terminating.\n")
    end
end

end #end module Oboe.DataIO


text:
function vnormal(v::Int,points::Array{Float64,2},faces::Array{Int,2},Iterator)

    #N_MWA = zeros(3)
    N_MWA = fill(0,3)

    for ti in Iterator
        face = faces[:,ti]

        w = face.==v
        cw = w[[3,1,2]]
        ccw = w[[2,3,1]]

        p1 = points[:,face[w]...]
        p2 = points[:,face[cw]...]
        p3 = points[:,face[ccw]...]

        a1 = p2 - p1
        a2 = p3 - p1

        alpha = asin(norm(cross(a1,a2))/norm(a1)/norm(a2))

        N_MWA += alpha*cross(a1,a2)

    end

    N_MWA = N_MWA/norm(N_MWA)

    return N_MWA
end

function NormalVectors!(n,points,faces,FaceVRing)
    for i in 1:size(points,2)
        n[:,i] = vnormal(i,points,faces,FaceVRing(i))
    end
end

function ZinchenkoDifertential!(vects,vnormal,R)

    Lx = [0 0 0; 0 0 -1; 0 1 0]
    Ly = [0 0 1; 0 0 0; -1 0 0]
    Lz = [0 -1 0; 1 0 0; 0 0 0]

    d = [0,0,1] + vnormal
    d /= norm(d)

    Ln = d[1]*Lx + d[2]*Ly + d[3]*Lz
    #R[:,:] = expm(pi*Ln)
    R[:,:] = exp(pi*Ln)

    # println("Rotation matrix is")
    # dump(R)

    for vj in 1:size(vects)[end]
        vects[:,vj] = R*vects[:,vj]
    end

    ### Construction of the system
    A = Array{Float64}(undef,3,3)
    B = Array{Float64}(undef,3)

    vects_norm2 = Array{Float64}(undef,size(vects)[end])
    for vj in 1:size(vects)[end]
       vects_norm2[vj] = norm(vects[:,vj])^2
    end

    A[1,1] = sum(vects[1,:].^4 ./vects_norm2)
    A[1,2] = sum(vects[1,:].^3 .*vects[2,:]./vects_norm2)
    A[1,3] = sum(vects[2,:].^2 .*vects[1,:].^2 ./vects_norm2)
    A[2,1] = A[1,2]
    A[2,2] = A[1,3]
    A[2,3] = sum(vects[2,:].^3 .*vects[1,:]./vects_norm2)
    A[3,1] = A[1,3]
    A[3,2] = A[2,3]
    A[3,3] = sum(vects[2,:].^4 ./vects_norm2)

    B[1] = sum(vects[3,:].*vects[1,:].^2 ./vects_norm2)
    B[2] = sum(vects[1,:].*vects[2,:].*vects[3,:]./vects_norm2)
    B[3] = sum(vects[2,:].^2 .*vects[3,:]./vects_norm2)

    C,D,E = A\B
    return C,D,E
end

### Orginals ieks utils.jl
function vcurvature(v::Int,points::Array{Float64,2},faces::Array{Int,2},normal,Iterator)

    #points = msh.points
    vring = Array{Int}(undef,0)
    for ti in Iterator #TriRing(v,msh)
        face = faces[:,ti]
        cw = (face.==v)[[3,1,2]]
        vi, = face[cw]
        push!(vring,vi)
    end

    p0 = points[:,v]
    vects = Array{Float64}(undef,3,length(vring))
    for (j,vi) in enumerate(vring)
        vects[:,j] = points[:,vi] - p0
    end

    #vnormall = vnormal(v,points,faces,Iterator) # If I would like to implement more precise curvature calculation
    #R = eye(3)
    R = Matrix{Float64}(I,3,3)
    C,D,E = ZinchenkoDifertential!(vects,normal,R)

    A = [C D/2;D/2 E]
    k1,k2 = eigvals(-A)
    H = (k1 + k2)/2

    return H*2
end

function MeanCurvatures!(curvatures,points,faces,n,VertexVRing)
    for i in 1:size(points,2)
        curvatures[i] = vcurvature(i,points,faces,n[:,i],VertexVRing(i))
    end
end

function FitEllipsoid(points)

    x = points[1,:][:]
    y = points[2,:][:]
    z = points[3,:][:]

    D = [x.*x+y.*y - 2*z.*z x.*x+z.*z-2*y.*y 2*x.*y 2*x.*z 2*y.*z 2*x 2*y 2*z 1 .+ 0*x]

    d2 = x .* x + y .* y + z .* z; #% the RHS of the llsq problem (y's)
    u = ( D' * D ) \ ( D' * d2 );  #% solution to the normal equations

    v = Array{Float64}(undef,10)
    v[1] = u[1] +     u[2] - 1;
    v[2] = u[1] - 2 * u[2] - 1;
    v[3] = u[2] - 2 * u[1] - 1;
    v[4:10] = u[3:9];

    # v the 10 parameters describing the ellipsoid / conic algebraically:
    #               Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz + J = 0
    ### Asuming that ellipsoid is placed at the middle and

    A = v[1]
    B = v[2]
    C = v[3]
    J = v[10]

    a = sqrt(-J/A)
    b = sqrt(-J/B)
    c = sqrt(-J/C)

    return a,b,c
end

function ellipsoid_parameters(points)

    warn("Deprecated syntax use FitEllipsoid instead")
    FitEllipsoid(points)
end


text:
# Singleton implementation that is thread-safe unique per async tasks

const context_store = Dict{UInt, Context}()
const context_store_lock = Base.Threads.SpinLock()

"""
    context([id::UInt], [container])

Get a new context object.  If `id` is not passed, then return a global
context for the current thread/task.  The `container` must be an object that
supports `push!`, `getindex`, `empty!`, and `length` functions.
The default is `Dict{Any,Any}`.
"""
function context(; id = default_context_id(), container = Dict())
    try
        lock(context_store_lock)
        return get!(context_store, id) do
            Context(id, container)
        end
    finally
        unlock(context_store_lock)
    end
end

"""
    default_context_id()

Return a context name for the current thread.
"""
function default_context_id()
    # This id must be task/thread-unique.
    # Current method reuse the most significant byte from the current task's pointer address.
    # Limitations are:
    # 1. Up to 256 threads are supported
    # 2. Current task's pointer address may spill over to the MSB.
    return convert(UInt, Base.Threads.threadid()) << 56 +
           convert(UInt, pointer_from_objref(current_task()))
end


text:
using Test
using Distributions

include("../src/utils.jl")
include("../src/sample_sobol.jl")
include("../src/test_functions/ishigami.jl")
include("../src/analyze_sobol.jl")

##
## 1. utils
##

# SobolData
N = 100
calc_second_order = false
parameters = OrderedDict(
    :param1 => Normal(0,1),
    :param2 => LogNormal(5, 20),
    :param3 => TriangularDist(0, 4, 1)
)

data1 = SobolData(params = parameters, calc_second_order = calc_second_order, N = N)
@test data1.params[:param1] == parameters[:param1]
@test data1.params[:param2] == parameters[:param2]
@test data1.params[:param2] == parameters[:param2]
@test data1.params[:param3] == parameters[:param3]
@test data1.N == N
@test data1.calc_second_order == calc_second_order

data2 = SobolData()
@test data2.params == nothing
@test data2.N == 1000
@test data2.calc_second_order == true

data3 = SobolData(params = parameters)
data4 = SobolData(calc_second_order = false)
data5 = SobolData(N = 100)

# scale_sobol_seq!
seq = rand(100, 8)
original_seq = copy(seq)
dists = [Normal(1, 0.2), Uniform(0.75, 1.25), LogNormal(20,4), TriangularDist(0, 4, 1)]
scale_sobol_seq!(seq, dists)
@test size(seq) == size(original_seq)
@test seq != original_seq

##
## 2. Sample Sobol
##

samples = sample(data1)
@test size(samples, 2) == length(data1.params)
if data1.calc_second_order
    @test size(samples,1) == data1.N * (2 * length(data1.params) + 2)
else
    @test size(samples,1) == data1.N * (length(data1.params) + 2)
end

samples3 = sample(data3)
@test size(samples3, 2) == length(data3.params)
if data3.calc_second_order
    @test size(samples3,1) == data3.N * (2 * length(data3.params) + 2)
else
    @test size(samples3,1) == data3.N * (length(data3.params) + 2)
end

@test_throws ErrorException sample(data2) # params are nothing
@test_throws ErrorException sample(data4) # params are nothing
@test_throws ErrorException sample(data5) # params are nothing

##
## 4. Analyze Sobol
##

Y1 = ishigami(samples)
results1 = analyze(data1, Y1)
for Si in results1[:firstorder]
    @test Si <= 1
end
for CI in results1[:firstorder_conf]
    @test CI > 0
end
for St in results1[:totalorder]
    @test St <= 1
end
for CI in results1[:totalorder_conf]
    @test CI > 0
end
@test sum(results1[:totalorder]) > sum(results1[:firstorder])

Y3 = ishigami(samples3)
results3 = analyze(data3, Y3)
for Si in results3[:firstorder]
    @test Si <= 1
end
for CI in results3[:firstorder_conf]
    @test ismissing(CI) || CI > 0
end
for S2i in results3[:secondorder]
    @test ismissing(S2i) || S2i <= 1
end
for CI in results3[:secondorder_conf]
    @test ismissing(CI) || CI > 0
end
for St in results3[:totalorder]
    @test St <= 1
end
for CI in results3[:totalorder_conf]
    @test ismissing(CI) || CI > 0
end
@test sum(results3[:totalorder]) > sum(results3[:firstorder])


text:
"""
    Kernel

Abstract type for representing kernel functions.
"""
abstract type Kernel end

"""
    ExactKernel <: Kernel

Abstract type for representing kernels that are evaluated exactly.
`ExactKernel`s must be callable with two inputs.
"""
abstract type ExactKernel <: Kernel end

"""
    RFFKernel <: Kernel

Abstract type for representing kernels that are approximated via random Fourier
features. `RFFKernel`s must be callable with one or three inputs.
"""
abstract type RFFKernel <: Kernel end

"""
    EuclideanKernel() <: ExactKernel

Create a kernel function where the kernel is just the Euclidean inner product
(i.e., ridge regression).
"""
struct EuclideanKernel <: ExactKernel end

"""
    (k::EuclideanKernel)(p, q)

Evaluate the Euclidean inner product between `p` and `q`.

# Arguments
- `p::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: First
  kernel input [Q,M] or \\[M\\] (if Q = 1) or scalar (if Q = M = 1)
- `q::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Second
  kernel input [Q,N] or \\[N\\] (if Q = 1) or scalar (if Q = N = 1)

# Note
- Q is the number of features
- M is the number of feature vectors in the first input
- N is the number of feature vectors in the second input

# Return
- `K::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Kernel
  output [M,N] or \\[M\\] (if N = 1) or \\[N\\] (if M = 1) or scalar
  (if M = N = 1)
"""
function (k::EuclideanKernel)(
    p::AbstractMatrix{<:Real},
    q::AbstractMatrix{<:Real}
)

    # Get the sizes of p and q
    M = size(p, 2)
    N = size(q, 2)

    # Compute the kernel function
    return [p[:,m]' * q[:,n] for m = 1:M, n = 1:N]

end

function (k::EuclideanKernel)(
    p::AbstractVector{<:Real}, # [M]
    q::AbstractMatrix{<:Real} # [Q,N]
)

    size(q, 1) == 1 ||
        throw(DimensionMismatch("p has feature dimension equal to 1, but q " *
                                "has a larger feature dimension"))

    return k(p, vec(q))

end

function (k::EuclideanKernel)(
    p::AbstractVector{<:Real},
    q::AbstractVector{<:Real}
)

    return p * q'

end

function (k::EuclideanKernel)(
    p::AbstractVector{<:Real},
    q::Real
)

    return p * q'

end

function (k::EuclideanKernel)(
    p::Real,
    q::AbstractVector{<:Real}
)

    return p * conj(q)

end

function (k::EuclideanKernel)(
    p::Real,
    q::Real
)

    return p * q'

end

"""
    GaussianKernel(Λ) <: ExactKernel

Create a Gaussian kernel function.

# Properties
- `Λ::Union{<:Real,AbstractVector{<:Real}}`: Length scales \\[Q\\] or scalar
  (if Q = 1)

# Note
- Q is the number of features
"""
struct GaussianKernel{T<:Union{<:Real,<:AbstractVector{<:Real}}} <: ExactKernel
    Λ::T

    GaussianKernel(Λ::Union{<:Real,<:AbstractVector{<:Real}}) = begin
        length(Λ) == 1 && (Λ = Λ[])
        new{typeof(Λ)}(Λ)
    end
end

"""
    GuassianKernel(Λy, [Λν])

Create a Guassian kernel function.

# Arguments
- `Λy::Union{<:Number,<:AbstractVector{<:Number}}`: Length scales for features
  \\[Q\\] or scalar (if Q = 1)
- `Λν::Union{<:Real,<:AbstractVector{<:Real}}`: Length scales for known
  parameters \\[K\\] or scalar (if K = 1)
"""
function GaussianKernel(
    Λ::Union{<:Complex,<:AbstractVector{<:Complex}}
)

    Λ = reduce(vcat, ([real(Λ[i]), imag(Λ[i])] for i = 1:length(Λ)))

    return GaussianKernel(Λ)

end

function GaussianKernel(
    Λy::Union{<:Number,<:AbstractVector{<:Number}},
    Λν::Union{<:Real,<:AbstractVector{<:Real}}
)

    if eltype(Λy) <: Complex
        Λy = reduce(vcat, ([real(Λy[i]), imag(Λy[i])] for i = 1:length(Λy)))
    end

    return GaussianKernel([Λy; Λν])

end

"""
    (k::GaussianKernel)(p, q)

Evaluate the Gaussian kernel.

# Arguments
- `p::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: First
  kernel input [Q,M] or \\[M\\] (if Q = 1) or scalar (if Q = M = 1)
- `q::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Second
  kernel input [Q,N] or \\[N\\] (if Q = 1) or scalar (if Q = N = 1)

## Note
- Q is the number of features
- M is the number of feature vectors in the first input
- N is the number of feature vectors in the second input

# Return
- `K::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Kernel
  output [M,N] or \\[M\\] (if N = 1) or \\[N\\] (if M = 1) or scalar
  (if M = N = 1)
"""
function (k::GaussianKernel)(
    p::AbstractMatrix{<:Real},
    q::AbstractMatrix{<:Real}
)

    # Get the sizes of p and q
    M = size(p, 2)
    N = size(q, 2)

    # Construct the square root of the Gaussian covariance matrix and invert it
    sqrtΣ = Diagonal(1 ./ k.Λ)

    # Compute the kernel function
    return [gaussiankernel(p[:,m], q[:,n], sqrtΣ) for m = 1:M, n = 1:N]

end

function (k::GaussianKernel)(
    p::AbstractVector{<:Real}, # [M]
    q::AbstractMatrix{<:Real} # [Q,N]
)

    size(q, 1) == 1 ||
        throw(DimensionMismatch("p has feature dimension equal to 1, but q " *
                                "has a larger feature dimension"))

    return k(p, vec(q))

end

function (k::GaussianKernel)(
    p::AbstractVector{<:Real},
    q::AbstractVector{<:Real}
)

    # Get the lengths of p and q
    M = length(p)
    N = length(q)

    # Construct the square root of the Gaussian covariance matrix and invert it
    sqrtΣ = 1 / k.Λ

    # Compute the kernel function
    return [gaussiankernel(p[m], q[n], sqrtΣ) for m = 1:M, n = 1:N]

end

function (k::GaussianKernel)(
    p::AbstractVector{<:Real},
    q::Real
)

    # Get the length of p
    M = length(p)

    # Construct the square root of the Gaussian covariance matrix and invert it
    sqrtΣ = 1 / k.Λ

    # Compute the kernel function
    return [gaussiankernel(p[m], q, sqrtΣ) for m = 1:M]

end

function (k::GaussianKernel)(
    p::Real,
    q::AbstractVector{<:Real}
)

    # Gaussian kernel is symmetric, so just swap the order of inputs
    return k(q, p)

end

function (k::GaussianKernel)(
    p::Real,
    q::Real
)

    # Construct the square root of the Gaussian covariance matrix and invert it
    sqrtΣ = 1 / k.Λ

    # Compute the kernel function
    return gaussiankernel(p, q, sqrtΣ)

end

"""
    gaussiankernel(p, q, sqrtΣ)

Compute the Gaussian kernel with covariance matrix Σ evaluated at `p` and `q`.
Note that function input is `sqrtΣ`, i.e., the square root of Σ.
"""
function gaussiankernel(p, q, sqrtΣ)

    return exp(-0.5norm(sqrtΣ * (p - q))^2)

end

"""
    GaussianRFF(H, Λ) <: RFFKernel

Create an approximate (via random Fourier features) Gaussian kernel function.

# Properties
- `H::Integer`: Approximation order
- `Λ::Union{<:Real,AbstractVector{<:Real}}`: Length scales \\[Q\\] or scalar (if
  Q = 1)

## Note
- Q is the number of features
"""
struct GaussianRFF{T1<:Integer,T2<:Union{<:Real,<:AbstractVector{<:Real}}} <: RFFKernel
    H::T1
    Λ::T2

    GaussianRFF(H::Integer, Λ::Union{<:Real,<:AbstractVector{<:Real}}) = begin
        length(Λ) == 1 && (Λ = Λ[])
        new{typeof(H),typeof(Λ)}(H, Λ)
    end
end

"""
    GuassianRFF(H, Λy, [Λν])

Create an approximate (via random Fourier features) Gaussian kernel function.

# Arguments
- `H::Integer`: Approximation order
- `Λy::Union{<:Number,<:AbstractVector{<:Number}}`: Length scales for features
  \\[Q\\] or scalar (if Q = 1)
- `Λν::Union{<:Real,<:AbstractVector{<:Real}}`: Length scales for known
  parameters \\[K\\] or scalar (if K = 1)
"""
function GaussianRFF(
    H::Integer,
    Λ::Union{<:Complex,<:AbstractVector{<:Complex}}
)

    Λ = reduce(vcat, ([real(Λ[i]), imag(Λ[i])] for i = 1:length(Λ)))

    return GaussianRFF(H, Λ)

end

function GaussianRFF(
    H::Integer,
    Λy::Union{<:Number,<:AbstractVector{<:Number}},
    Λν::Union{<:Real,<:AbstractVector{<:Real}}
)

    if eltype(Λy) <: Complex
        Λy = reduce(vcat, ([real(Λy[i]), imag(Λy[i])] for i = 1:length(Λy)))
    end

    return GaussianRFF(H, [Λy; Λν])

end

"""
    (k::GaussianRFF)([rng], q)
    (k::GaussianRFF)(q, f, phase)

Evaluate the approximate Gaussian kernel.

# Arguments
- `rng::AbstractRNG = Random.GLOBAL_RNG`: Random number generator to use to
  generate `f` and `phase`
- `q::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Kernel
  input [Q,N] or \\[N\\] (if Q = 1) or scalar (if Q = N = 1)
- `f::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}} = randn(k.H, Q)`:
  Unscaled random frequency values [H,Q] or \\[H\\] (if Q = 1)
- `phase::AbstractVector{<:Real} = rand(k.H)`: Random phase values \\[H\\]

## Note
- Q is the number of features
- N is the number of feature vectors in the input
- H is the approximation order for the random Fourier features

# Return
- `z::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`:
  Higher-dimensional features [H,N] or \\[H\\] (if N = 1)
- `freq::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Random
  frequency values [H,Q] or \\[H\\] (if Q = 1)
- `phase::AbstractVector{<:Real}`: Random phase values \\[H\\]
"""
function (k::GaussianRFF)(
    rng::AbstractRNG,
    q::AbstractMatrix{<:Real}
)

    # Get the first dimension of q
    Q = size(q, 1)

    # Generate random phase and unscaled frequency values
    f = randn(rng, k.H, Q)
    phase = rand(rng, k.H)

    return k(q, f, phase)

end

function (k::GaussianRFF)(
    rng::AbstractRNG,
    q::Union{<:Real,<:AbstractVector{<:Real}}
)

    # Generate random phase and unscaled frequency values
    f = randn(rng, k.H)
    phase = rand(rng, k.H)

    return k(q, f, phase)

end

function (k::GaussianRFF)(
    q
)

    return k(Random.GLOBAL_RNG, q)

end

function (k::GaussianRFF)(
    q::AbstractMatrix{<:Real},
    f::AbstractMatrix{<:Real},
    phase::AbstractVector{<:Real}
)

    # Construct the covariance matrix from which to draw the Gaussian samples
    # and take the square root
    sqrtΣ = Diagonal(div0.(1, 2π .* k.Λ))

    # Scale the random frequency values
    freq = f * sqrtΣ

    # Map the features to a higher dimensional space via random Fourier features
    # Also return freq and phase for use later
    return (rffmap(q, freq, phase), freq, phase)

end

function (k::GaussianRFF)(
    p::Union{<:Real,<:AbstractVector{<:Real}},
    f::AbstractMatrix{<:Real},
    phase::AbstractVector{<:Real}
)

    size(f, 2) == 1 ||
        throw(DimensionMismatch("p has feature dimension equal to 1, but f " *
                                "has a larger feature dimension"))

    return k(p, vec(f), phase)

end

function (k::GaussianRFF)(
    q::Union{<:Real,<:AbstractVector{<:Real}},
    f::AbstractVector{<:Real},
    phase::AbstractVector{<:Real}
)

    # Construct the covariance matrix from which to draw the Gaussian samples
    # and take the square root
    sqrtΣ = div0(1, 2π * k.Λ)

    # Scale the random frequency values
    freq = f * sqrtΣ

    # Map the features to a higher dimensional space via random Fourier features
    # Also return freq and phase for use later
    return (rffmap(q, freq, phase), freq, phase)

end

"""
    rffmap(q, freq, phase)

Map features to a higher dimensional space via random Fourier features.

# Arguments
- `q::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`:
  Lower-dimensional features [Q,N] or \\[N\\] (if Q = 1) or scalar
  (if Q = N = 1)
- `freq::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}`: Random
  frequency values [H,Q] or \\[H\\] (if Q = 1)
- `phase::AbstractVector{<:Real}`: Random phase values \\[H\\]

## Note
- Q is the number of features
- N is the number of feature vectors in the input
- H is the approximation order for the random Fourier features

# Return
- `z::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`:
  Higher-dimensional features [H,N] or \\[H\\] (if N = 1)
"""
function rffmap(
    q::AbstractMatrix{<:Real},
    freq::AbstractMatrix{<:Real},
    phase::AbstractVector{<:Real}
)

    return sqrt(div0(2, length(phase))) .* cos.(2π .* (freq * q .+ phase))

end

function rffmap(
    q::AbstractMatrix{<:Real},
    freq::AbstractVector{<:Real},
    phase::AbstractVector{<:Real}
)

    size(q, 1) == 1 ||
        throw(DimensionMismatch("freq has feature dimension equal to 1, but " *
                                "q has a larger feature dimension"))

    return rffmap(vec(q), freq, phase)

end

function rffmap(
    q::AbstractVector{<:Real},
    freq::AbstractVector{<:Real},
    phase::AbstractVector{<:Real}
)

    return sqrt(div0(2, length(phase))) .* cos.(2π .* (freq * transpose(q) .+ phase))

end

function rffmap(
    q::Real,
    freq::AbstractVector{<:Real},
    phase::AbstractVector{<:Real}
)

    return sqrt(div0(2, length(phase))) .* cos.(2π .* (freq .* q .+ phase))

end


text:
# This file is part of JuliaFMI.
# Licensed under MIT: https://github.com/AnHeuermann/JuliaFMI/blob/master/LICENSE.txt

"""
    fmi2SetCallbackFunction!(fmu::FMU;stepFinished::Function=nothing,componentEnvironment::Function=nothing)

Helper function to add a function declared in julia to the CallbackFunctions struct on the fmu.

...
# Arguments
- `fmu::FMU`:
Optional args:
- `stepFinished::Function=nothing`:
- `componentEnvironment::fmi2ComponentEnvironment=nothing`:
...

# Example
'''
'''
"""
function fmi2SetCallbackFunctions!(fmu::FMU; stepFinished::T=nothing,
   componentEnvironment::T=nothing) where {T<:Union{Nothing,Function}, U<:Union{Nothing,AbstractFMI2ComponentEnvironment}}
   fmi2SetCallbackFunctions!(fmu.fmi2CallbackFunctions, stepFinished=stepFinished, componentEnvironment=componentEnvironment)
end

"""
    fmi2SetCallbackFunction!(functions::CallbackFunctions;stepFinished::Function=nothing,componentEnvironment::ComponentEnvironment=nothing,)where{ComponentEnvironment<:Function}

Helper function to add a function declared in julia to the CallbackFunctions struct.

...
# Arguments
- `functions::CallbackFunctions`:
Optional args:
- `stepFinished::Function=nothing`:
- `componentEnvironment::AbstractFMI2ComponentEnvironment=nothing`:
- `where{ComponentEnvironment<:Function}`:
...

# Example
'''
'''
"""
function fmi2SetCallbackFunctions!(functions::CallbackFunctions; stepFinished::T=nothing, componentEnvironment::U=nothing
    ) where {T<:Union{Nothing,Function}, U<:Union{Nothing,AbstractFMI2ComponentEnvironment}}
    if !isnothing(stepFinished)
        GC.@preserve cfunc = @cfunction($stepFinished, Cvoid, (Ref(U), Cint))
        functions.stepFinished = cfunc.ptr
    end
    if !isnothing(componentEnvironment)
        GC.@preserve cfunc = @cfunction($componentEnvironment, Cvoid, ())
        functions.componentEnvironment = cfunc.ptr
    end
end

max_stars_repo_name: AnHeuermann/JuliaFMU

max_issues_repo_name: AnHeuermann/JuliaFMU

max_forks_repo_name: jwscook/JuliaFMI


text:
using Makie
using MakieLayout
using DifferentialEquations

f(u,p,t) = 1.01*u
u0 = 1/2
tspan=(0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

scene, layout = layoutscene()

ax = layout[1,1] = LAxis(scene)

ls = [lines!(ax, sol.t, sol, color=:blue, linewidth=5),
      lines!(ax, sol.t, t->0.5*exp(1.01*t), color=:red, linewidth=3, linestyle=:dash)]

leg = layout[1, 1] = LLegend(scene, ls, ["solve", "true"];
                             width = Auto(), height = Auto(), halign = :left, valign = :top,
                             tellwidth = false, tellheight = false,
                             margin = (20, 20, 15, 15))
display(scene)
readline()

save("plot.png", scene)

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
using Heuristics
using Test, SafeTestsets

@time begin
#@time @safetestset "Particle tests..." begin include("particleTests.jl") end
#@time @safetestset "Swarm tests..." begin include("swarmTests.jl") end
@time @safetestset "PSO optimization test..." begin include("psoOptTests.jl") end
@time @safetestset "MS_PSO optimization test..." begin include("mspsoOptTests.jl") end
end


text:
using Test
using ModelingToolkit

using Distributed

# Test that the precompiled ODE system works
@everywhere push!(LOAD_PATH, joinpath(@__DIR__, "precompile_test"))

using ODEPrecompileTest

u  = collect(1:3)
p  = collect(4:6)

# These cases do not work, because they get defined in the ModelingToolkit's RGF cache.
@test parentmodule(typeof(ODEPrecompileTest.f_bad.f.f_iip).parameters[2]) == ModelingToolkit
@test parentmodule(typeof(ODEPrecompileTest.f_bad.f.f_oop).parameters[2]) == ModelingToolkit
@test parentmodule(typeof(ODEPrecompileTest.f_noeval_bad.f.f_iip).parameters[2]) == ModelingToolkit
@test parentmodule(typeof(ODEPrecompileTest.f_noeval_bad.f.f_oop).parameters[2]) == ModelingToolkit
@test_throws KeyError ODEPrecompileTest.f_bad(u, p, 0.1)
@test_throws KeyError ODEPrecompileTest.f_noeval_bad(u, p, 0.1)

# This case works, because it gets defined with the appropriate cache and context tags.
@test parentmodule(typeof(ODEPrecompileTest.f_noeval_good.f.f_iip).parameters[2]) == ODEPrecompileTest
@test parentmodule(typeof(ODEPrecompileTest.f_noeval_good.f.f_oop).parameters[2]) == ODEPrecompileTest
@test ODEPrecompileTest.f_noeval_good(u, p, 0.1) == [4, 0, -16]


text:
using DelimitedFiles
using LOWDER

include("generators.jl")

function runtest_hs(
                    filename::String
                    )

    file = open( filename, "w" )

    for i = 1:100

        try

            # Generates the problem
            ( x, l, u, fmin ) = problem_generator_hs( i )

            n = length(x)
            p = length(fmin)

            n_points = n + 1

            # Solves the problem using 'lowder'
            sol = LOWDER.lowder( fmin, x, l, u; m = n_points )

            # Saves info
            text = "$( n ) $( p ) $( sol.status ) $( sol.true_val ) $( sol.iter ) $( sol.nf ) $( sol.nf / p ) $( sol.stationarity ) $( sol.sample_radius ) $( sol.tr_radius ) $( sol.index ) $( sol.f )"
            print( file, text )
            for j = 1 : ( n - 1 )

                print( file, " $(sol.solution[ j ])" )

            end
            println( file, " $( sol.solution[ n ] )" )

            println("succes!")

        catch

            println( file, "execution_fail" )

            println("fail!")

        end

    end

    close( file )

end

# -----------------------------------------------
# Path to file
# -----------------------------------------------

filename = "../data_files/hs_cons.dat"

# -----------------------------------------------
# Funtion call
# -----------------------------------------------

runtest_hs( filename )


text:
include("../../src/model_selection/train_test_split.jl")
using Test, DataFrames, Random

test = DataFrame(
    x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
    y = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
)

@testset "train test split" begin
    @test isequal(
        partition_train_test(test, 0.5, 42)[1],
        DataFrame(x = [7, 3, 5, 6, 9], y = [7, 3, 5, 6, 9]),
    )
    @test isequal(
        train_test_split(test, "y", seed = 42)[1],
        DataFrame(x = [7, 3, 5, 6, 9, 8]),
    )
    @test isequal(
        train_test_split(test, "y", seed = 42)[2],
        DataFrame(x = [2, 1, 10, 4]),
    )
    @test isequal(train_test_split(test, "y", seed = 42)[3], [7, 3, 5, 6, 9, 8])
    @test isequal(train_test_split(test, "y", seed = 42)[4], [2, 1, 10, 4])
    @test isequal(
        train_test_split(test, "y", train_share = 0.5, seed = 42)[1],
        DataFrame(x = [7, 3, 5, 6, 9]),
    )
end


text:
module ModelAbstraction

using CompEcon

# ------------------------------------------------------------------- #
# Economic Entities
# ------------------------------------------------------------------- #
abstract AbstractEntity
abstract AbstractAgent <: AbstractEntity
abstract AbstractConsumer <: AbstractAgent
abstract AbstractProducer <: AbstractAgent
abstract AbstractGovernment <: AbstractAgent  # is the government an entity or an agent?

immutable NotImplementedError <: Exception
    msg::String
end

# AbstractConsumer Interface
for f in [:is_recursive, :utility]
    @eval $(f){T<:AbstractConsumer}(::T, args...) =
        throw(NotImplementedError("$f not implemented for type $T"))
end

include("MA_Consumers.jl")
include("MA_Producers.jl")

# ------------------------------------------------------------------- #
# State Variables
# ------------------------------------------------------------------- #
abstract AbstractVariable
abstract AbstractStateVariable <: AbstractVariable
abstract AbstractJumpVariable <: AbstractVariable

immutable EndogenousState{T,BF,BP,TM} <: AbstractStateVariable
    grid::Vector{T}
    basis::Basis{1,BF,BP}
    Φ::TM
    # name::Symbol  # could default to `gensym()` (which is julia for I don't care)
end

EndogenousState(grid::Vector) =
    EndogenousState(grid, Basis(LinParams(grid, 0)))

# Note depending on the type of Basis you pass, the grid inside the state
# instance might be different than the grid you passed in.
function EndogenousState(grid::Vector, b::Basis{1})
    real_grid = nodes(b)[1]
    Φ = CompEcon.evalbase(b.params[1], real_grid, 0)[1]
    EndogenousState(real_grid, b, Φ)
end

immutable ExogenousState{T,BF,BP,TM} <: AbstractStateVariable
    grid::Vector{T}   # grid for realizations in the current period
    gridp::Matrix{T}  # grid moving each i ∈ grid  with all possible shocks
    basis::Basis{1,BF,BP}
    Φ::TM
    # name::Symbol
end

ExogenousState{T}(grid::Vector{T}, gridp::Matrix{T}) =
    ExogenousState(grid, gridp, Basis(LinParams(grid, 0)))

function ExogenousState(grid::Vector, gridp, b::Basis{1})
    real_grid = nodes(b)[1]
    if real_grid != grid
        msg = string("You must use `grid = nodes(b)` and compute `gridp`",
                     " based on that grid or things won't line up")
        throw(ArgumentError(msg))
    end
    Φ = CompEcon.evalbase(b.params[1], real_grid, 0)[1]
    ExogenousState(real_grid, gridp, b, Φ)
end

# ------------------------------------------------------------------- #
# Exogenous Process
# ------------------------------------------------------------------- #
abstract AbstractExogenousProcess{N}

include("MA_Exogenous.jl")

# ------------------------------------------------------------------- #
# State Spaces
# ------------------------------------------------------------------- #
abstract AbstractStateSpace{Nendog,Nexog}

immutable DefaultStateSpace{Nendog,Nexog,T,N} <: AbstractStateSpace{N}
    endo::NTuple{Nendog,EndogenousState{T}}
    exog::NTuple{Nexog,ExogenousState{T}}
    grid::Matrix{T}
    grid_transpose::Matrix{T}
    basis::Basis{N}
end

function DefaultStateSpace{N1,N2,T,N}(endog::NTuple{N1,EndogenousState{T}},
                                      exog::NTuple{N2,ExogenousState{T}},
                                      grid::Matrix{T},
                                      b::Basis{N})
    DefaultStateSpace(endog, exog, grid, grid', b)
end

function DefaultStateSpace{N1,N2,T}(endog::NTuple{N1,EndogenousState{T}},
                                    exog::NTuple{N2,ExogenousState{T}})
    grids = vcat(Vector[i.grid for i in endog], Vector[i.grid for i in exog])
    grid = gridmake(grids...)
    basis = Basis([i.basis for i in endog]..., [i.basis for i in exog]...)
    DefaultStateSpace{N1,N2,T}(endog, exog, grid, basis)
end

# function DefaultStateSpace(b::Basis, nendog::Int)
#     grid, grids = nodes(b)
#     N = ndims(B)
#     nexog = N - nendog
#     endog = tuple([EndogenousState(grids[i], b[i]) for i=1:nendog]...)
#     exog = tuple([(grids[i], b[i]) for i=nendog+1:N]...)

# function DefaultStateSpace{N1,N2,T}(endog::NTuple{N1,Vector{T}},
#                                     exog::NTuple{N2,Vector{T}})

# ------------------------------------------------------------------- #
# Frictions
# ------------------------------------------------------------------- #
include("MA_Frictions.jl")

export AbstractEntity,
       AbstractConsumer, EpsteinZinAgent, utility, is_recursive,
       AbstractProducer, CESProducer, f_k,
       AbstractVariable,
       AbstractStateVariable, EndogenousState, ExogenousState,
       AbstractExogenousProcess,
       ConstantVolatility1, ConstantVolatility, StochasticVolatility, simulate,
       AbstractStateSpace, DefaultStateSpace,
       AbstractAdjustmentCost, AdjCost, _ac, _dIac, _dkac,
       _unpack

end  # module

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
function retrieve(c::Client, ::Type{User}, user::Integer)
    return if me(c) !== nothing && me(c).id == user
        get_current_user(c)
    else
        get_user(c, user)
    end
end
retrieve(c::Client, ::Type{User}) = get_current_user(c)

update(c::Client; kwargs...) = modify_current_user(c; kwargs...)


text:
module TestMathleticsFiles

# using DataFrames
using MathleticsFiles: dataset
using Test

@testset "Test of single dataset" begin
    df = dataset("rushing")
    @test size(df) == (18, 4)
    df = dataset("nfl_team_totals")
    @test size(df) == (128, 13)
end

end  # module TestMathleticsFiles


text:
module NonlinearEigenproblems

using Reexport

"Include the specified module, then `@reexport` it."
macro include_export(file_name, module_name)
    :( include($file_name); @reexport using .$module_name )
end

# Re-export core functionality, to remove the need for "using" statements.
# Non-core functionality not re-exported below can be accessed by "using"
# those modules, for example "using NonlinearEigenproblems.RKHelper".
# Note: Order matters below, due to inter-module dependencies.
include(joinpath("utils", "Serialization.jl"))
@include_export("NEPCore.jl", NEPCore)
@include_export("NEPTypes.jl", NEPTypes)
@include_export("LinSolvers.jl", LinSolvers)
include(joinpath("rk_helper", "RKHelper.jl"))
@include_export("NEPTransformations.jl", NEPTransformations)
@include_export("NEPSolver.jl", NEPSolver)
@include_export("Gallery.jl", Gallery)

end


text:
using Test
using LogicCircuits
using DataFrames: DataFrame, DataFrameRow, nrow, ncol
using CUDA: CUDA

@testset "Data utils" begin

    m = [1.1 2.1; 3.1 4.1; 5.1 6.1]
    df = DataFrame(m)
    dfb = DataFrame(BitMatrix([true false; true true; false true]))

    batched_df = batch(df, 1)
    batched_dfb = batch(dfb, 1)

    @test num_examples(df) == 3
    @test num_examples(dfb) == 3
    @test num_examples(batched_df) == 3
    @test num_examples(batched_dfb) == 3

    @test num_features(df) == 2
    @test num_features(dfb) == 2
    @test num_features(batched_df) == 2
    @test num_features(batched_dfb) == 2

    @test example(df,2) isa Vector
    @test example(df,2)[1] == 3.1
    @test example(df,2)[2] == 4.1

    @test feature_values(df,2) == [2.1,4.1,6.1]
    @test feature_values(dfb,2) isa BitVector
    @test feature_values(dfb,2) == BitVector([false,true,true])

    @test isfpdata(df)
    @test !isfpdata(dfb)
    @test isfpdata(batched_df)
    @test !isfpdata(batched_dfb)
    @test !isfpdata(DataFrame([1 "2"; 3 "4"]))

    @test !isbinarydata(df)
    @test isbinarydata(dfb)
    @test !isbinarydata(batched_df)
    @test isbinarydata(batched_dfb)

    @test num_examples(shuffle_examples(df)) == 3
    @test 1.1 in feature_values(shuffle_examples(df), 1)

    dft, _, _ = threshold(df, nothing, nothing)

    @test feature_values(dft,1) == [false, false, true]

    @test LogicCircuits.Utils.fully_factorized_log_likelihood(dfb; pseudocount=1) ≈ -1.280557674335465 #not verified
    @test LogicCircuits.Utils.fully_factorized_log_likelihood(dfb) ≈ -1.2730283365896256 #not verified

    @test ll_per_example(-12.3, dfb) ≈ -4.1 #not verified

    @test bits_per_pixel(-12.3, dfb) ≈ 2.9575248338223754 #not verified

    dfb = DataFrame(BitMatrix([true false; true true; false true]))

    @test !isweighted(dfb)

    weights1 = DataFrame(weight = [0.6, 0.6, 0.6])
    weights2 = [0.6, 0.6, 0.6]
    wdfb1 = add_sample_weights(dfb, weights1)
    wdfb2 = add_sample_weights(dfb, weights2)

    dfb_split1 = split_sample_weights(wdfb1)[1]

    @test get_weights(wdfb1)[1] ≈ 0.6

    @test isweighted(wdfb1)
    @test isweighted(wdfb2)
    @test !isweighted(dfb_split1)

    dataset = DataFrame([true false false; true true true; false false true; false true false])
    weights = DataFrame(weight = [0.6, 0.6, 0.6, 0.8])
    dataset = add_sample_weights(dataset, weights)
    bag_datasets = bagging_dataset(dataset; num_bags = 5, frac_examples = 1.0);
    @test bag_datasets isa Array{DataFrame}
    @test nrow(bag_datasets[1]) == 4
    @test ncol(bag_datasets[1]) == ncol(dataset) - 1

    bag_datasets = bagging_dataset(dataset; num_bags = 5, frac_examples = 1.0, batch_size = 1);
    @test length(bag_datasets[1]) == 4
    @test ncol(bag_datasets[1][1]) == ncol(dataset) - 1

    df = DataFrame(BitMatrix([true true; false true; false false]))
    mar = marginal_prob(df)
    @test mar[1] ≈ 0.333333333 atol = 1e-6
    @test mar[2] ≈ 0.666666666 atol = 1e-6
    sdf = soften(df, 0.001; scale_by_marginal = false)
    @test sdf[1,1] ≈ 0.999 atol = 1e-6
    @test sdf[2,1] ≈ 0.001 atol = 1e-6
    @test sdf[3,1] ≈ 0.001 atol = 1e-6
    @test sdf[1,2] ≈ 0.999 atol = 1e-6
    @test sdf[2,2] ≈ 0.999 atol = 1e-6
    @test sdf[3,2] ≈ 0.001 atol = 1e-6
    sdf = soften(df, 0.001; scale_by_marginal = true)
    @test sdf[1,1] ≈ 0.999 + 0.001 * 0.3333333 atol = 1e-6
    @test sdf[2,1] ≈ 0.001 * 0.3333333 atol = 1e-6
    @test sdf[3,1] ≈ 0.001 * 0.3333333 atol = 1e-6
    @test sdf[1,2] ≈ 0.999 + 0.001 * 0.6666666 atol = 1e-6
    @test sdf[2,2] ≈ 0.999 + 0.001 * 0.6666666 atol = 1e-6
    @test sdf[3,2] ≈ 0.001 * 0.6666666 atol = 1e-6

    if CUDA.functional()
        wdfb1_gpu = to_gpu(wdfb1)
        wdfb2_gpu = to_gpu(wdfb2)
        dfb_split1_gpu = to_gpu(dfb_split1)

        dfb_gpu_split1 = split_sample_weights(wdfb1_gpu)[1]

        @test isgpu(wdfb1_gpu)
        @test isgpu(wdfb2_gpu)
        @test isgpu(dfb_split1_gpu)
        @test isgpu(dfb_gpu_split1)
        @test isgpu(to_gpu(wdfb1_gpu))

        dfb = DataFrame(BitMatrix([true false; true true; false true]))
        weights1 = DataFrame(weight = [0.6, 0.6, 0.6])
        weights2 = [0.6, 0.6, 0.6]
        batched_dfb = batch(dfb, 1)
        batched_weights1 = batch(weights1, 1)
        batched_weights2 = batch(weights2, 1)

        @test isbatched(batched_dfb)
        @test !isbatched(dfb)

        @test isweighted(add_sample_weights(batched_dfb, weights1))
        @test isweighted(add_sample_weights(batched_dfb, weights2))
        @test isweighted(add_sample_weights(batched_dfb, batched_weights1))
        @test isweighted(add_sample_weights(batched_dfb, batched_weights2))

        batched_wdfb = add_sample_weights(batched_dfb, weights1)

        @test !isweighted(split_sample_weights(batched_wdfb)[1])

        @test isgpu(to_gpu(batched_wdfb))

        @test get_weights(batched_wdfb)[1][1] ≈ 0.6
    end
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    SDFileReader,
    sdfilereader,
    nohaltsupplier,
    sdftomol,
    RDFileReader,
    rdfilereader,
    nohalt_rxn_supplier,
    rxntoreaction

const SDF_CHARGE_TABLE = Dict(
    0 => 0, 1 => 3, 2 => 2, 3 => 1, 4 => 0, 5 => -1, 6 => -2, 7 => -3
)

function sdfatom(line)
    xpos = parse(Float64, line[1:10])
    ypos = parse(Float64, line[11:20])
    zpos = parse(Float64, line[21:30])
    coords = Float64[xpos, ypos, zpos]
    sym = Symbol(rstrip(line[32:34]))
    # atom.mass_diff = parse(Int, line[35:36]) use ISO property
    old_sdf_charge = parse(Int, line[37:39])
    charge = SDF_CHARGE_TABLE[old_sdf_charge]
    multi = old_sdf_charge == 4 ? 2 : 1
    # atom.stereo_flag = parse(Int, line[40:42])
    # valence = parse(Int, line[46:48])
    return (sym, charge, multi, coords)
end

function sdfbond(line)
    u = parse(Int, line[1:3])
    v = parse(Int, line[4:6])
    order = parse(Int, line[7:9])
    notation = parse(Int, line[10:12])
    return (u, v, order, notation)
end

function sdfprops(lines)
    props = Dict{Int,Dict{Symbol,Real}}() # atomindex, {type => value}
    for line in lines
        startswith(line, "M  ") || continue
        proptype = line[4:6]
        if !(proptype in ("CHG", "RAD", "ISO"))
            continue # Other properties are not supported yet
        end
        count = parse(Int, line[7:9])
        for c in 1:count
            i = c - 1
            idx = parse(Int, line[8i + 11 : 8i + 13])
            val = line[8i + 15 : 8i + 17]
            haskey(props, idx) || (props[idx] = Dict{Symbol,Real}())
            if proptype == "CHG"
                props[idx][:CHG] = parse(Int, val)
            elseif proptype == "RAD"
                props[idx][:RAD] = parse(Int, val)
            elseif proptype == "ISO"
                props[idx][:ISO] = parse(Float64, val)
            end
        end
    end
    return props
end

function sdfoptions(lines)
    function add_attribute!(d, key, list)
        if key != :NOTANOPTION
            while !isempty(list) && isempty(list[end])
                pop!(list)
            end
            d[key] = length(list) == 1 ? list[1] : copy(list)
            empty!(list)
        end
        return d
    end

    data = Dict{Symbol,Union{String,Vector{String}}}()
    lastoption = Ref(:NOTANOPTION)
    prevlines = String[]
    for (i, line) in enumerate(lines)
        # Some inappropriate signs are accepted for practical use
        m = match(r">.*?<([\w -.%=/]+)>", line)
        if m !== nothing
            add_attribute!(data, lastoption[], prevlines)
            lastoption[] = Symbol(m[1])
        elseif line == "\$\$\$\$"  # $$$$ indicates the next compound
            add_attribute!(data, lastoption[], prevlines)
            lastoption[] = :NOTANOPTION
            break
        else
            push!(prevlines, line)
        end
    end
    # Ensure that an unterminated file includes the final option
    add_attribute!(data, lastoption[], prevlines)

    return data
end

"""
    parse(::Type{SDFile}, lines)

Parse lines of a SDFile mol block data into a molecule object.
"""
function Base.parse(::Type{SDFile}, sdflines)
    sdflines = collect(sdflines)
    molend = findnext(x -> x == "M  END", sdflines, 1)
    lines = @view sdflines[1:molend-1]
    optlines = @view sdflines[molend+1:end]

    # Get element blocks
    countline = lines[4]

    occursin("V3000", countline) && return sdftomol3000(join(sdflines, "\n"))

    atomcount = parse(UInt16, countline[1:3])
    bondcount = parse(UInt16, countline[4:6])
    # chiralflag = countline[12:15] Not used
    # propcount = countline[30:33] No longer supported
    atomoffset = 5
    bondoffset = atomcount + atomoffset
    propoffset = bondoffset + bondcount
    atomblock = @view lines[atomoffset:bondoffset-1]
    bondblock = @view lines[bondoffset:propoffset-1]
    propblock = @view lines[propoffset:end]

    # Parse atoms
    nodeattrs = SDFileAtom[]
    props = sdfprops(propblock)
    for (i, line) in enumerate(atomblock)
        (sym, charge, multi, coords) = sdfatom(line)
        mass = nothing
        if !isempty(props)
            # If prop block exists, any annotations in atom blocks are ignored
            charge = 0
            multi = 1
            if haskey(props, i)
                prop = props[i]
                haskey(prop, :CHG) && (charge = prop[:CHG])
                haskey(prop, :RAD) && (multi = prop[:RAD])
                haskey(prop, :ISO) && (mass = prop[:ISO])
            end
        end
        push!(nodeattrs, SDFileAtom(sym, charge, multi, mass, coords))
    end

    # Parse bonds
    edges = Tuple{Int,Int}[]
    edgeattrs = SDFileBond[]
    for line in bondblock
        (u, v, order, notation) = sdfbond(line)
        push!(edges, (u, v))
        push!(edgeattrs, SDFileBond(order, notation))
    end

    molobj = graphmol(edges, nodeattrs, edgeattrs)
    merge!(molobj.attributes, sdfoptions(optlines))
    return molobj
end

"""
    parse(::Type{GraphReaction}, lines)

Parse lines of a rxn file into a reaction object.
"""
function Base.parse(::Type{GraphReaction}, rxn::AbstractString)
    reactants, products = if startswith(rxn, raw"$RXN V3000")
        reactants = match(blockregex("REACTANT"), rxn).match
        products = match(blockregex("PRODUCT"), rxn).match
        rr = String[m.captures[1] for m in eachmatch(blockregex("CTAB"), reactants)]
        pp = String[m.captures[1] for m in eachmatch(blockregex("CTAB"), products)]
        sdftomol3000.(rr), sdftomol3000.(pp)
    else
        parts = split(rxn, r"\$MOL\r?\n")
        (ne, np) = parse.(Int, split(split(first(parts), r"\r?\n")[end-1]))
        rr = parts[2:1 + ne]
        pp = parts[(2 + ne):(1 + ne + np)]
        sdftomol.(IOBuffer.(rr)), sdftomol.(IOBuffer.(pp))
    end
    GraphReaction(reactants, products)
end

Base.parse(::Type{GraphReaction}, rxnlines) = parse(GraphReaction, join(rxnlines, "\n"))

function nohaltsupplier(block)
    return try
        mol = parse(SDFile, block)
        setdiastereo!(mol)
        setstereocenter!(mol)
        return mol
    catch e
        if e isa ErrorException
            println("$(e.msg) (#$(i) in sdfilereader)")
            return graphmol(SDFileAtom, SDFileBond)
        else
            throw(e)
        end
    end
end

function nohalt_rxn_supplier(block)
    return try
        reaction = parse(GraphReaction, block)
        setdiastereo!(reaction)
        setstereocenter!(reaction)
        return reaction
    catch e
        if e isa ErrorException
            println("$(e.msg) (#$(i) in rxnfilereader)")
            return Reaction()
        else
            throw(e)
        end
    end
end

struct SDFileReader
    lines::Base.EachLine
    parser::Function
end

struct RDFileReader
    io::IO
    parser::Function
end

function Base.iterate(reader::SDFileReader, state=nothing)
    block = String[]
    next = iterate(reader.lines)
    while next !== nothing
        (line, state) = next
        if startswith(line, raw"$$$$")
            return (reader.parser(block), state)
        end
        push!(block, rstrip(line))
        next = iterate(reader.lines, state)
    end
    if !isempty(block)
        return (reader.parser(block), state)
    end
    return
end

function Base.iterate(reader::RDFileReader, state=nothing)
    block = readuntil(reader.io, "\$RXN")

    # if the file is a normal rxn file, block is empty and one needs to read one more time
    isempty(block) && (block = readuntil(reader.io, "\$RXN"))
    while startswith(block, "\$RDFILE") && block != ""
        block = readuntil(reader.io, "\$RXN")
    end
    isempty(block) ? nothing : (reader.parser("\$RXN" * block), state)
end

Base.IteratorSize(::Type{SDFileReader}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{SDFileReader}) = Base.EltypeUnknown()

Base.IteratorSize(::Type{RDFileReader}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{RDFileReader}) = Base.EltypeUnknown()

"""
    sdfilereader(file::IO)
    sdfilereader(path::AbstractString)

Read SDFile data from input stream (or a file path as a string) and return a
lazy iterator that yields molecule objects.

`sdfilereader` does not stop and raise errors when an erroneous or incompatible
SDFile block is read but produces an error message and yields an empty molecule.
If this behavior is not desirable, you can use the customized supplier function
instead of default supplier `nohaltsupplier`

```
function customsupplier()
    mol = try
        parse(SDFile, block)
    catch e
        throw(ErrorException("incompatible molecule found, aborting..."))
    end
end

function sdfilereader(file::IO)
    return SDFileReader(eachline(file), customsupplier)
end
```
"""
sdfilereader(file::IO) = SDFileReader(eachline(file), nohaltsupplier)
sdfilereader(path::AbstractString) = sdfilereader(open(path))

rdfilereader(file::IO) = RDFileReader(file, nohalt_rxn_supplier)
rdfilereader(path::AbstractString) = rdfilereader(open(path))

"""
    sdftomol(lines) -> GraphMol{SDFileAtom,SDFileBond}
    sdftomol(file::IO) -> GraphMol{SDFileAtom,SDFileBond}
    sdftomol(path::AbstractString) -> GraphMol{SDFileAtom,SDFileBond}

Read a SDFile(.sdf or .mol) and parse it into a molecule object. The given
argument should be a file input stream, a file path as a string or an iterator
that yields each sdfile text lines.
"""
function sdftomol(lines)
    mol = parse(SDFile, lines)
    setdiastereo!(mol)
    setstereocenter!(mol)
    return mol
end
sdftomol(file::IO) = sdftomol(eachline(file))
sdftomol(path::AbstractString) = sdftomol(open(path))

"""
    rxntoreaction(lines) -> GraphReaction
    rxntoreaction(file::IO) -> GraphReaction
    sdftomol(path::AbstractString) -> GraphReaction

Read a RXN file and parse it into a reaction object. The given
argument should be a file input stream, a file path as a string or an iterator
that yields each sdfile text lines.
"""
function rxntoreaction(lines)
    reaction = parse(GraphReaction, lines)
    setdiastereo!(reaction)
    setstereocenter!(reaction)
    return reaction
end
rxntoreaction(file::IO) = rxntoreaction(eachline(file))
rxntoreaction(path::AbstractString) = rxntoreaction(open(path))

# support of extended mol file format (V3000)

function sympair(s)::Pair{Symbol, Union{String, Int}}
    p = split(s,"=")
    length(p) != 2 && return :nothing => ""
    int = tryparse(Int, p[2])
    Symbol(p[1]) => isnothing(int) ? p[2] : int
end

blockregex(s::AbstractString) = Regex("M  V30 BEGIN $(s)\r?\n(.*?)\r?\nM  V30 END $(s)", "s")

function parseatomblock3000(atomblock)
    nodeattrs = SDFileAtom[]
    for line in eachline(IOBuffer(atomblock))
        ss = split(line)
        coords = parse.(Float64, ss[5:7])
        sym = Symbol(ss[4])

        props = Dict(sympair.(ss[9:end])...)

        charge = get(props, :CHG, 0)
        mass   = get(props, :MASS, nothing)
        multi  = get(props, :RAD, 1)

        push!(nodeattrs, SDFileAtom(sym, charge, multi, mass, coords))
    end
    nodeattrs
end

function parsebondblock3000(bondblock)
    edges = Tuple{Int,Int}[]
    edgeattrs = SDFileBond[]

    for line in eachline(IOBuffer(bondblock))
        ss = split(line)
        order, u, v = parse.(Int, ss[4:6])
        props = Dict(sympair.(ss[7:end])...)
        notation = get(props, :CFG, 0)
        push!(edges, (u, v))
        push!(edgeattrs, SDFileBond(order, notation))
    end
    edges, edgeattrs
end

function sdftomol3000(s::AbstractString)
    atomblock = match(blockregex("ATOM"), s).captures[1]
    bondblock = match(blockregex("BOND"), s).captures[1]

    nodeattrs = parseatomblock3000(atomblock)
    edges, edgeattrs = parsebondblock3000(bondblock)

    mol = graphmol(edges, nodeattrs, edgeattrs)
    # setdiastereo!(mol)
    # setstereocenter!(mol)
end


text:
# Copyright 2020 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

using Test
using Blades
using LinearAlgebra
using Combinatorics

module g3
  using Blades
  @generate_basis("+++")
end

module pg2
  using Blades
  @generate_basis("0++")
end

module g4
  using Blades
  @generate_basis("++++", false, true)
end

module pg3
  using Blades
  @generate_basis("0+++")
end

module g2
  using Blades
  @generate_basis("++")
end

module G5
  using Blades
  @generate_basis("+++++",false,true,true)
end

using .pg3
using .g4
using .g3
using .pg2
using .g2
using .G5

partial(f,x...) = (y...)->f(x...,y...)

"""
We can complete the definition of the geometric product.
All Blades are simple k-vectors.  So only one of a∧b or a⋅b is ever nonzero.
By defining addition of a Blade with 0 we can complete the Geometric Algebra of Blades.
u*v == u∧v + u⋅v
Not defined in the Blades package since it is defined in Multivectors
"""
Base.:(+)(a::B, s::Z) where {B<:Blade,Z<:Val{0.0}} = a
Base.:(+)(a::B, s::T) where {B<:Blade,T<:Real} = a + Val(Float64(s))
Base.:(+)(s::T, a::B) where {B<:Blade,T<:Real} = a + s

@testset "PGA3D" begin
  e₁, e₂, e₃, e₄, e₁₂, e₁₃, e₁₄, e₂₃, e₂₄, e₃₄, e₁₂₃, e₁₂₄, e₁₃₄, e₂₃₄, e₁₂₃₄ = alle(pg3,4)

  a = e₁(1.0); b = e₂(2.0); c = e₃(3.0); d = e₄(4.0)

  @test e₁(1.0)*e₁(1.0) == 0.0
  @test e₂(1)*e₂(1) == 1
  @test a*b*a == 0.0
  @test a*b == -b*a
  @test a*b*c == b*c*a
  @test b*c*b == -b.x*b.x*c
  @test a*b == e₁₂(a.x*b.x)
  @test reverse(a*b*c*d) == d*c*b*a == reverse((a*b)*(c*d)) == (d*c*b)*a
  @test e₃(2.0)*e₂₃₄(3.0) == -e₂₄(6.0) == e₂₄(-6.0)
  @test grade(1) == 0
  @test grade(a*b*a) == 0
  @test grade(a*b*c) == 3

  B = b*c

  @test a*B == a∧B
  @test b*B == b⋅B
  @test d*B == d∧B
  @test a*b == a∧b + a⋅b
  @test c*B == c∧B + c⋅B
  @test B*a == B∧a + B⋅a
  @test c*c == c∧c + c⋅c

  @test 1.0∧a == a
  @test 1.0⋅a == a

  @test norm(b) == 2.0 == det(b)
  @test norm(normalize(b)) == 1.0
  @test norm(normalize(-b)) == 1.0
  @test magnitude(normalize(-b)) == -1.0
  @test Blades.norm_sqr(b) == 4.0 == det(b,b)
  @test conj(a*b) == b*a == reverse(a*b) == ~(a*b)

  @test factor(a*b) == (scalar(a)*scalar(b), e₁, e₂)

  @test basis_1blades(e₁) == basis_kblades(e₁, 1)
  @test length(basis_kblades(e₂₃, 2)) == 6
  @test length(basis_kblades(e₂₃, 3)) == length(basis_kblades(e₂₃, 1))

  g₁, g₂, g₃, g₄, g₁₂, g₁₃, g₁₄, g₂₃, g₂₄, g₃₄, g₁₂₃, g₁₂₄, g₁₃₄, g₂₃₄, g₁₂₃₄ = alle(g4,4)

  @test dual(dual(1g₁)) == ⟂(⟂(1g₁))
  @test subspace(dual(dual(1e₁))) == subspace(⟂(⟂(1g₁)))
  @test subspace(Blades.nonmetric_dual(Blades.nonmetric_dual(3e₁))) == subspace(⟂(⟂(3g₁)))
  @test sign(dual(dual(1pg2.e₁))) == sign(⟂(⟂(1g3.e₁)))
  @test sign(Blades.nonmetric_dual(Blades.nonmetric_dual(3pg2.e₁))) == sign(⟂(⟂(3g3.e₁)))
  @test sign(dual(3pg2.e₁))== sign(⟂(3g3.e₁))
  @test sign(Blades.nonmetric_dual(3pg2.e₁)) == sign(⟂(3g3.e₁))

  # convenience operators on type
  @test e₁∧e₂ == e₁₂ == e₁*e₂
  @test e₂*e₁ == e₂∧e₁ == Nothing
  @test 2(e₁∧e₂) == e₁₂(2)
end

@testset "G2" begin

  e₁, e₂, e₁₂ = alle(g2,2)
  eᵢ  = one.(alle( g2, 2)[1:2])
  #More generally, on any flat space we have
  #⋆(dxi1 ∧dxi2 ∧···∧dxik) = dxik+1 ∧dxik+2 ∧···∧dxin, where (i1,i2,...,in) is any even permutation of (1,2,...,n)

  for p in Iterators.filter(p->parity(p)==0,permutations(1:2))
    for s in 1:1
      @test ⋆(reduce(∧,eᵢ[p[1:s]])) == reduce(∧,eᵢ[p[s+1:end]])
    end
  end
  @test ⋆(1e₁) == 1e₂
  @test ⋆(1e₂) == -1e₁
  ot = -1.0e₁₂
  @test ot∧⋆(ot) == (ot⋅ot)pseudoscalar(ot)
  @test ⋆(2.0e₁₂) == -2.0
  ot = 2.0e₁₂
  @test ot∧⋆(ot) == (ot⋅ot)pseudoscalar(ot)

  @test raise(e₁) == alld(g2,2)[1]
  @test lower(raise(e₁₂)) == e₁₂
end

#==
@testset "DEC course notes G5" begin
  eᵢ  = one.(alle( G5, 5)[1:5])
  #More generally, on any flat space we have
  #⋆(dxi1 ∧dxi2 ∧···∧dxik) = dxik+1 ∧dxik+2 ∧···∧dxin, where (i1,i2,...,in) is any even permutation of (1,2,...,n)

  for p in Iterators.filter(p->parity(p)==0,permutations(1:5))
    for s in 1:4
      @test ⋆(reduce(∧,eᵢ[p[1:s]])) == reduce(∧,eᵢ[p[s+1:end]])
      w = reduce(∧,eᵢ[p[1:s]])
    end
  end
end
==#

@testset "hodge star" begin
  rb = alle(G5, 5)[1:5]
  𝑖 = mapreduce(eᵢ->1.0*eᵢ, ∧, rb[1:4])
  for i in 1:100
    A = rand()rand(rb[1:3])∧rand()rand(rb[1:3])
    B = rand()rand(rb[1:3])∧rand()rand(rb[1:3])
    @test A∧⋆(B, 𝑖) == (A⋅B)*𝑖
  end

  for i in 1:1000
    A = rand()rand(rb[1:5])
    B = rand()rand(rb[1:5])∧rand()rand(rb[1:5])
    @test A∧⋆(B) == (A⋅B)pseudoscalar(A)
    @test B∧⋆(A) == (B⋅A)pseudoscalar(A)
    C = rand()rand(rb[1:5])∧rand()rand(rb[1:5])∧rand()rand(rb[1:5])
    @test A∧⋆(C) == (A⋅C)pseudoscalar(A)
    @test C∧⋆(A) == (C⋅A)pseudoscalar(A)
    @test B∧⋆(C) == (B⋅C)pseudoscalar(A)
  end
end


text:
using PyPlot

function calculate_midpoint(a,b)
  x,y = ((a[1]+b[1])/2), ((a[2]+b[2])/2)
  return (x,y)
end

function sierpinski(maxiter,pts)
  while maxiter > 0
    a = pts[4]
    b = pts[rand(1:3)]
    pts[4] = calculate_midpoint(a,b)
    ax[:scatter](pts[4][1],pts[4][2])
    maxiter = maxiter -1
  end
end

fig,ax = subplots()
p = [(1.0,0.0),(2.5,5.0),(5.0,0.0),(0.0,0.0)]
sierpinski(1500,p)

max_stars_repo_name: krudutta/Sierpinski-Gasket

max_issues_repo_name: krudutta/Sierpinski-Gasket

max_forks_repo_name: krudutta/Sierpinski-Gasket