text:
import Base: convert, real, imag, float, eps
import Base: isfinite, isnan, isinf, isless
import Base: trunc, ceil, floor, round

function evalf(b::Basic, bits::Integer=53, real::Bool=false)
    c = Basic()
    bits > 53 && real && (have_mpfr || throw(ArgumentError("libsymengine has to be compiled with MPFR for this feature")))
    bits > 53 && !real && (have_mpc || throw(ArgumentError("libsymengine has to be compiled with MPC for this feature")))
    status = ccall((:basic_evalf, libsymengine), Cint, (Ref{Basic}, Ref{Basic}, Culong, Cint), c, b, Culong(bits), Int(real))
    if status == 0
        return c
    else
        throw(ArgumentError("symbolic value cannot be evaluated to a numeric value"))
    end
end

## Conversions from SymEngine -> Julia at the ccall level
function convert(::Type{BigInt}, b::BasicType{Val{:Integer}})
    a = BigInt()
    c = Basic(b)
    ccall((:integer_get_mpz, libsymengine), Nothing, (Ref{BigInt}, Ref{Basic}), a, c)
    return a
end

function convert(::Type{BigFloat}, b::BasicType{Val{:RealMPFR}})
    c = Basic(b)
    a = BigFloat()
    ccall((:real_mpfr_get, libsymengine), Nothing, (Ref{BigFloat}, Ref{Basic}), a, c)
    return a
end

function convert(::Type{Cdouble}, b::BasicType{Val{:RealDouble}})
    c = Basic(b)
    return ccall((:real_double_get_d, libsymengine), Cdouble, (Ref{Basic},), c)
end

function real(b::BasicType{Val{:ComplexDouble}})
    c = Basic(b)
    a = Basic()
    ccall((:complex_double_real_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c)
    return a
end

function imag(b::BasicType{Val{:ComplexDouble}})
    c = Basic(b)
    a = Basic()
    ccall((:complex_double_imaginary_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c)
    return a
end

function real(b::BasicType{Val{:Complex}})
    c = Basic(b)
    a = Basic()
    ccall((:complex_real_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c)
    return a
end

function imag(b::BasicType{Val{:Complex}})
    c = Basic(b)
    a = Basic()
    ccall((:complex_imaginary_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c)
    return a
end

function real(b::BasicType{Val{:ComplexMPC}})
    c = Basic(b)
    a = Basic()
    ccall((:complex_mpc_real_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c)
    return a
end

function imag(b::BasicType{Val{:ComplexMPC}})
    c = Basic(b)
    a = Basic()
    ccall((:complex_mpc_imaginary_part, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}), a, c)
    return a
end

##################################################
# N
"""

Convert a SymEngine numeric value into a Julian number

"""
N(a::Integer) = a
N(a::Rational) = a
N(a::Complex) = a
N(b::Basic) = N(BasicType(b))

function N(b::BasicType{Val{:Integer}})
    a = convert(BigInt, b)
    if (a.size > 1 || a.size < -1)
        return a
    elseif (a.size == 0)
        return 0
    else
        u = unsafe_load(a.d, 1)
        w = signed(u)
        if w > 0
            return w * a.size
        elseif a.size == 1
            return u
        else
            return a
        end
    end
end

N(b::BasicType{Val{:Rational}}) = Rational(N(numerator(b)), N(denominator(b))) # TODO: conditionally wrap rational_get_mpq from cwrapper.h
N(b::BasicType{Val{:RealDouble}}) = convert(Cdouble, b)
N(b::BasicType{Val{:RealMPFR}}) = convert(BigFloat, b)
N(b::BasicType{Val{:NaN}}) = NaN
N(b::BasicType{Val{:Infty}}) = (string(b) == "-inf") ? -Inf : Inf

## Mapping of SymEngine Constants into julia values
constant_map = Dict("pi" => π, "eulergamma" => γ, "exp(1)" => e, "catalan" => catalan,
                    "goldenratio" => φ)

N(b::BasicType{Val{:Constant}}) = constant_map[toString(b)]

N(b::BasicComplexNumber) = complex(N(real(b)), N(imag(b)))
function N(b::BasicType)
    b = convert(Basic, b)
    fs = free_symbols(b)
    if length(fs) > 0
        throw(ArgumentError("Object can have no free symbols"))
    end
    out = evalf(b)
    imag(out) == Basic(0.0) ? real(out) : out
end

##  Conversions SymEngine -> Julia
function as_numer_denom(x::Basic)
    a, b = Basic(), Basic()
    ccall((:basic_as_numer_denom, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}, Ref{Basic}), a, b, x)
    return a, b
end

as_numer_denom(x::BasicType) = as_numer_denom(Basic(x))
denominator(x::SymbolicType) = as_numer_denom(x)[2]
numerator(x::SymbolicType)   = as_numer_denom(x)[1]

## Complex
real(x::Basic) = Basic(real(SymEngine.BasicType(x)))
real(x::SymEngine.BasicType) = x

imag(x::Basic) = Basic(imag(SymEngine.BasicType(x)))
imag(x::BasicType{Val{:Integer}}) = Basic(0)
imag(x::BasicType{Val{:RealDouble}}) = Basic(0)
imag(x::BasicType{Val{:RealMPFR}}) = Basic(0)
imag(x::BasicType{Val{:Rational}}) = Basic(0)
imag(x::SymEngine.BasicType) = throw(InexactError())

## define convert(T, x) methods leveraging N()
convert(::Type{Float64}, x::Basic)           = convert(Float64, N(evalf(x, 53, true)))
convert(::Type{BigFloat}, x::Basic)          = convert(BigFloat, N(evalf(x, precision(BigFloat), true)))
convert(::Type{Complex{Float64}}, x::Basic)  = convert(Complex{Float64}, N(evalf(x, 53, false)))
convert(::Type{Complex{BigFloat}}, x::Basic) = convert(Complex{BigFloat}, N(evalf(x, precision(BigFloat), false)))
convert(::Type{Number}, x::Basic)            = x
convert(::Type{T}, x::Basic) where {T <: Real}      = convert(T, N(x))
convert(::Type{Complex{T}}, x::Basic) where {T <: Real}    = convert(Complex{T}, N(x))

# Constructors no longer fall back to `convert` methods
Base.Int64(x::Basic) = convert(Int64, x)
Base.Float64(x::Basic) = convert(Float64, x)
Base.BigInt(x::Basic) = convert(BigInt, x)
Base.Real(x::Basic) = convert(Real, x)

## For generic programming in Julia
float(x::Basic) = float(N(x))

# trunc, flooor, ceil, round, rem, mod, cld, fld,
isfinite(x::Basic) = x-x == 0
isnan(x::Basic) = isnan(N(x))
isinf(x::Basic) = !isnan(x) & !isfinite(x)
isless(x::Basic, y::Basic) = isless(N(x), N(y))

## These should have support in symengine-wrapper, but currently don't
trunc(x::Basic, args...) = Basic(trunc(N(x), args...))
trunc(::Type{T},x::Basic, args...) where {T <: Integer} = convert(T, trunc(x,args...))

ceil(x::Basic) = Basic(ceil(N(x)))
ceil(::Type{T},x::Basic) where {T <: Integer} = convert(T, ceil(x))

floor(x::Basic) = Basic(floor(N(x)))
floor(::Type{T},x::Basic) where {T <: Integer} = convert(T, floor(x))

round(x::Basic) = Basic(round(N(x)))
round(::Type{T},x::Basic) where {T <: Integer} = convert(T, round(x))

prec(x::BasicType{Val{:RealMPFR}}) = ccall((:real_mpfr_get_prec, libsymengine), Clong, (Ref{Basic},), x)

# eps
eps(x::Basic) = eps(BasicType(x))
eps(x::BasicType{T}) where {T} = eps(typeof(x))
eps(::Type{T}) where {T <: BasicType} = 0
eps(::Type{T}) where {T <: Basic} = 0
eps(::Type{BasicType{Val{:RealDouble}}}) = 2^-52
eps(::Type{BasicType{Val{:ComplexDouble}}}) = 2^-52
eps(x::BasicType{Val{:RealMPFR}}) = evalf(Basic(2), prec(x), true) ^ (-prec(x)+1)
eps(x::BasicType{Val{:ComplexMPFR}}) = eps(real(x))


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

using Markdown
using InteractiveUtils

# ╔═╡ 3909a581-7d5d-49da-a7b6-7d574a9a8e67
begin
	using CSV
	using DataFrames
	using MLDataUtils
end

# ╔═╡ 973690c7-001f-4110-b1be-5066123f067d
begin
	using ScikitLearn
	using DecisionTree
	using Statistics
end

# ╔═╡ 7bc884cc-f825-11eb-05ac-9d97ec1c19b5
md"# Housing Prices Competition for Kaggle Learn Users"

# ╔═╡ 8fc7b9a9-5edb-445f-acfc-efbf80289567
md"## Download datasets"

# ╔═╡ f0bb2843-6351-46e6-a9af-6001a2c01bec
md"Please make sure `kaggle` command is installed. [kaggle-api github](https://github.com/Kaggle/kaggle-api)"

# ╔═╡ dfa20ffd-442d-4fe9-915f-76f176f9ee3a
begin
	run(`kaggle competitions download -c home-data-for-ml-course -f train.csv`)
	run(`kaggle competitions download -c home-data-for-ml-course -f test.csv`)
end

# ╔═╡ 88ed2deb-b2b8-41fc-adc7-0b79fb345bf6
md"## Prepare Datasets"

# ╔═╡ 1ad6c7fa-a9f6-4886-8424-64bc7c44117c
md"### Import Modules"

# ╔═╡ edb88725-c559-4dc9-b8b1-e45e728800a3
md"### Load Training Dataset"

# ╔═╡ 34ba5692-1a67-4fe8-9727-af2df3a0e776
begin
	home_data = CSV.File("train.csv") |> DataFrame
	train_home_data, val_home_data = splitobs(home_data, at = 0.9)
end

# ╔═╡ 029acc00-ca97-4c77-b954-7f818988ac65
md"### Select Features"

# ╔═╡ ea1795ae-946e-4579-b684-9608d7734c3c
features = [
	"LotArea",
	"YearBuilt",
	"1stFlrSF",
	"2ndFlrSF",
	"FullBath",
	"BedroomAbvGr",
	"TotRmsAbvGrd"
]

# ╔═╡ b910f0d5-afa2-464a-bd9b-33ec9aee5c20
md"### Training and Validation Datasets"

# ╔═╡ eb9d7236-785b-40ed-ae42-97f9a65814ee
begin
	X = home_data[!, features]
	y = home_data.SalePrice
	train_X = train_home_data[!, features]
	train_y = train_home_data.SalePrice
	val_X = val_home_data[!, features]
	val_y = val_home_data.SalePrice
end

# ╔═╡ 0280809c-fe3c-4f4f-b88e-f6ebae7d4008
md"### Testing Datasets"

# ╔═╡ 11bfdd61-a0ae-4018-9dfc-65135c72b286
begin
	test_home_data = CSV.File("test.csv") |> DataFrame
	test_X = test_home_data[!, features]
end

# ╔═╡ c3e2c05f-fe5d-41e3-bd77-ae67b055c6f1
md"## Train Models"

# ╔═╡ f69630ee-d18d-4cce-a7ed-2b8dc5decd3f
md"### Import Modules"

# ╔═╡ 290266cc-4aa5-4fc0-8b5c-8a7f5ad9b6e8
md"### Mean Absolute Error"

# ╔═╡ 4d8eea32-5062-425b-bd4c-0d78f1a9a7f5
function get_mae(val, pred)
	return mean(abs(sum(val .- pred)))
end

# ╔═╡ 0df54ead-3090-44c2-9612-514f88500264
md"### Decision Tree Regressor"

# ╔═╡ 09251e92-1ced-497c-8aa2-9c18c90ac688
begin
	dtr_model = DecisionTreeRegressor(rng=1)
	fit!(dtr_model, Matrix(train_X), train_y)
	dtr_pred_y = ScikitLearn.predict(dtr_model, Matrix(val_X))
	dtr_mae = get_mae(val_y, dtr_pred_y)
end

# ╔═╡ c7ec23d9-354c-44c1-a24e-e3b06f77bd26
begin
	full_dtr_model = DecisionTreeRegressor(rng=1)
	fit!(full_dtr_model, Matrix(X), y)
	dtr_result = ScikitLearn.predict(full_dtr_model, Matrix(test_X))
end

# ╔═╡ a11012ea-469e-4140-8e53-fc06d24e6fcb
begin
	dtr_df = DataFrame(
		Id = test_home_data.Id,
		SalePrice = dtr_result
	)
	CSV.write("dtr_result.csv", dtr_df)
end

# ╔═╡ 27932df1-5f3a-4ec7-bb99-3a3e8f94b582
md"### Random Forest Regressor"

# ╔═╡ 2d1f6590-cb93-4773-9e01-54cbbfece0c9
begin
	rfr_model = RandomForestRegressor(rng=1)
	fit!(rfr_model, Matrix(train_X), train_y)
	rfr_pred_y = ScikitLearn.predict(rfr_model, Matrix(val_X))
	rfr_mae = get_mae(val_y, rfr_pred_y)
end

# ╔═╡ 0106bd03-e5a8-40ab-8cf9-a995a011ad69
begin
	full_rfr_model = RandomForestRegressor(rng=1)
	fit!(full_rfr_model, Matrix(X), y)
	rfr_result = ScikitLearn.predict(full_rfr_model, Matrix(test_X))
end

# ╔═╡ ceda4142-51a5-4bfa-a26e-b619ec6c3754
begin
	rfr_df = DataFrame(
		Id = test_home_data.Id,
		SalePrice = rfr_result
	)
	CSV.write("rfr_result.csv", rfr_df)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DecisionTree = "7806a523-6efd-50cb-b5f6-3fa6f1930dbb"
MLDataUtils = "cc2ba9b6-d476-5e6d-8eaf-a92d5412d41d"
ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
CSV = "~0.8.5"
DataFrames = "~1.2.2"
DecisionTree = "~0.10.10"
MLDataUtils = "~0.5.4"
ScikitLearn = "~0.6.4"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "344f143fa0ec67e47917848795ab19c6a455f32c"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.32.0"

[[Conda]]
deps = ["JSON", "VersionParsing"]
git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.5.2"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.9"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DecisionTree]]
deps = ["DelimitedFiles", "Distributed", "LinearAlgebra", "Random", "ScikitLearnBase", "Statistics", "Test"]
git-tree-sha1 = "8b58db7954a6206399d9f66ef1a328da8c0f1d19"
uuid = "7806a523-6efd-50cb-b5f6-3fa6f1930dbb"
version = "0.10.10"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
deps = ["Test"]
git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.0.0"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LearnBase]]
deps = ["LinearAlgebra", "StatsBase"]
git-tree-sha1 = "47e6f4623c1db88570c7a7fa66c6528b92ba4725"
uuid = "7f8f8fb0-2700-5f03-b4bd-41f8cfc144b6"
version = "0.3.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MLDataPattern]]
deps = ["LearnBase", "MLLabelUtils", "Random", "SparseArrays", "StatsBase"]
git-tree-sha1 = "e99514e96e8b8129bb333c69e063a56ab6402b5b"
uuid = "9920b226-0b2a-5f5f-9153-9aa70a013f8b"
version = "0.5.4"

[[MLDataUtils]]
deps = ["DataFrames", "DelimitedFiles", "LearnBase", "MLDataPattern", "MLLabelUtils", "Statistics", "StatsBase"]
git-tree-sha1 = "ee54803aea12b9c8ee972e78ece11ac6023715e6"
uuid = "cc2ba9b6-d476-5e6d-8eaf-a92d5412d41d"
version = "0.5.4"

[[MLLabelUtils]]
deps = ["LearnBase", "MappedArrays", "StatsBase"]
git-tree-sha1 = "3211c1fdd1efaefa692c8cf60e021fb007b76a08"
uuid = "66a33bbf-0c2b-5fc8-a008-9da813334f0a"
version = "0.5.6"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[MappedArrays]]
git-tree-sha1 = "e8b359ef06ec72e8c030463fe02efe5527ee5142"
uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
version = "0.4.1"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.0"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.2"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "bfd7d8c7fd87f04543810d9cbd3995972236ba1b"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.2.1"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "169bb8ea6b1b143c5cf57df6d34d022a7b60c6db"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.92.3"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[ScikitLearn]]
deps = ["Compat", "Conda", "DataFrames", "Distributed", "IterTools", "LinearAlgebra", "MacroTools", "Parameters", "Printf", "PyCall", "Random", "ScikitLearnBase", "SparseArrays", "StatsBase", "VersionParsing"]
git-tree-sha1 = "ccb822ff4222fcf6ff43bbdbd7b80332690f168e"
uuid = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
version = "0.6.4"

[[ScikitLearnBase]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "7877e55c1523a4b336b433da39c8e8c08d2f221f"
uuid = "6e75b9c4-186b-50bd-896f-2d2496a4843e"
version = "0.5.0"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "a3a337914a035b2d59c9cbe7f1a38aaba1265b02"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.6"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "fed1ec1e65749c4d96fc20dd13bea72b55457e62"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.9"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "d0c690d37c73aeb5ca063056283fde5585a41710"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.2.0"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─7bc884cc-f825-11eb-05ac-9d97ec1c19b5
# ╟─8fc7b9a9-5edb-445f-acfc-efbf80289567
# ╟─f0bb2843-6351-46e6-a9af-6001a2c01bec
# ╠═dfa20ffd-442d-4fe9-915f-76f176f9ee3a
# ╟─88ed2deb-b2b8-41fc-adc7-0b79fb345bf6
# ╟─1ad6c7fa-a9f6-4886-8424-64bc7c44117c
# ╠═3909a581-7d5d-49da-a7b6-7d574a9a8e67
# ╟─edb88725-c559-4dc9-b8b1-e45e728800a3
# ╠═34ba5692-1a67-4fe8-9727-af2df3a0e776
# ╟─029acc00-ca97-4c77-b954-7f818988ac65
# ╠═ea1795ae-946e-4579-b684-9608d7734c3c
# ╟─b910f0d5-afa2-464a-bd9b-33ec9aee5c20
# ╠═eb9d7236-785b-40ed-ae42-97f9a65814ee
# ╟─0280809c-fe3c-4f4f-b88e-f6ebae7d4008
# ╠═11bfdd61-a0ae-4018-9dfc-65135c72b286
# ╟─c3e2c05f-fe5d-41e3-bd77-ae67b055c6f1
# ╟─f69630ee-d18d-4cce-a7ed-2b8dc5decd3f
# ╠═973690c7-001f-4110-b1be-5066123f067d
# ╟─290266cc-4aa5-4fc0-8b5c-8a7f5ad9b6e8
# ╠═4d8eea32-5062-425b-bd4c-0d78f1a9a7f5
# ╟─0df54ead-3090-44c2-9612-514f88500264
# ╠═09251e92-1ced-497c-8aa2-9c18c90ac688
# ╠═c7ec23d9-354c-44c1-a24e-e3b06f77bd26
# ╠═a11012ea-469e-4140-8e53-fc06d24e6fcb
# ╟─27932df1-5f3a-4ec7-bb99-3a3e8f94b582
# ╠═2d1f6590-cb93-4773-9e01-54cbbfece0c9
# ╠═0106bd03-e5a8-40ab-8cf9-a995a011ad69
# ╠═ceda4142-51a5-4bfa-a26e-b619ec6c3754
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl

max_stars_repo_name: even-wei/kaggle-with-julia

max_issues_repo_name: even-wei/kaggle-with-julia

max_forks_repo_name: even-wei/kaggle-with-julia


text:
using DFTK
using LinearAlgebra
using PyCall
using Test
using Unitful
include("testcases.jl")

@testset "bzmesh_uniform agrees with spglib" begin
    function test_against_spglib(kgrid_size; kshift=[0, 0, 0])
        kgrid_size = Vec3(kgrid_size)
        is_shift = ifelse.(kshift .== 0, 0, 1)
        diagonal = Matrix{Int64}(I, 3, 3)
        n_kpts, _, grid =
            DFTK.spglib_get_stabilized_reciprocal_mesh(kgrid_size, [diagonal]; is_shift)

        kcoords_spglib = [(kshift .+ grid[ik]) .// kgrid_size for ik in 1:n_kpts]
        kcoords_spglib = DFTK.normalize_kpoint_coordinate.(kcoords_spglib)
        sort!(kcoords_spglib)

        kcoords, _ = bzmesh_uniform(kgrid_size, kshift=kshift)
        sort!(kcoords)

        @test kcoords == kcoords_spglib
    end

    test_against_spglib([ 2,  3,  2])
    test_against_spglib([ 3,  3,  3])
    test_against_spglib([ 3,  3,  3], kshift=[1//2, 0, 0])
    test_against_spglib([ 2,  3,  4])
    test_against_spglib([ 9, 11, 13])
end

@testset "bzmesh_ir_wedge is correct reduction" begin
    function test_reduction(testcase, kgrid_size, kirredsize;
                            supercell=[1, 1, 1], kshift=[0, 0, 0])
        lattice   = testcase.lattice
        atoms     = testcase.atoms
        positions = testcase.positions

        if supercell != (1, 1, 1)  # Make a supercell
            ase       = ase_atoms(lattice, atoms, positions) * supercell
            lattice   = load_lattice(ase)
            atoms     = load_atoms(ase)
            positions = load_positions(ase)
        end

        red_kcoords, _ = bzmesh_uniform(kgrid_size; kshift)
        symmetries = DFTK.symmetry_operations(lattice, atoms, positions)
        irred_kcoords, ksymops = bzmesh_ir_wedge(kgrid_size, symmetries; kshift)

        @test length(irred_kcoords) == kirredsize

        # Try to reproduce all kcoords from irred_kcoords
        all_kcoords = Vector{Vec3{Rational{Int}}}()
        sym_preserving_grid = DFTK.symmetries_preserving_kgrid(symmetries, red_kcoords)
        for (ik, k) in enumerate(irred_kcoords)
            append!(all_kcoords, [symop.S * k for symop in sym_preserving_grid])
        end

        # Normalize the obtained k-points and test for equality
        red_kcoords = unique(sort([mod.(k .* kgrid_size, kgrid_size) for k in red_kcoords]))
        all_kcoords = unique(sort([mod.(k .* kgrid_size, kgrid_size) for k in all_kcoords]))
        @test all_kcoords == red_kcoords
    end

    test_reduction(silicon, [ 2,  3,  2],   6)
    test_reduction(silicon, [ 3,  3,  3],   4)
    test_reduction(silicon, [ 2,  3,  4],  14)
    test_reduction(silicon, [ 9, 11, 13], 644)

    test_reduction(silicon, [ 3,  3,  3], 6, kshift=[1//2, 1//2, 1//2])
    test_reduction(silicon, [ 3,  3,  3], 6, kshift=[1//2, 0, 1//2])
    test_reduction(silicon, [ 3,  3,  3], 6, kshift=[0, 1//2, 0])

    test_reduction(silicon, [ 1,  4,  4],    7, supercell=[2, 1, 1])
    test_reduction(silicon, [ 1,  16,  16], 73, supercell=[4, 1, 1])

    test_reduction(magnesium, [ 2,  3,  2],   8)
    test_reduction(magnesium, [ 3,  3,  3],   6)
    test_reduction(magnesium, [ 2,  3,  4],  12)
    test_reduction(magnesium, [ 9, 11, 13], 350)

    test_reduction(platinum_hcp, [5, 5, 5], 63)
end

@testset "standardize_atoms" begin
    # Test unperturbed structure
    std = standardize_atoms(silicon.lattice, silicon.atoms, silicon.positions, primitive=true)
    @test length(std.atoms) == 2
    @test std.atoms[1].symbol == :Si
    @test std.atoms[2].symbol == :Si
    @test length(std.positions) == 2
    @test std.positions[1] - std.positions[2] ≈ 0.25ones(3)
    @test std.lattice ≈ silicon.lattice

    # Perturb structure
    plattice   = silicon.lattice .+ 1e-8rand(3, 3)
    patoms     = silicon.atoms
    ppositions = [p + 1e-8rand(3) for p in silicon.positions]
    std = standardize_atoms(plattice, patoms, ppositions, primitive=true)

    # And check we get the usual silicon primitive cell back:
    a = std.lattice[1, 2]
    @test std.lattice == [0  a  a; a 0 a; a a 0]
    @test length(std.atoms) == 2
    @test std.atoms[1].symbol == :Si
    @test std.atoms[2].symbol == :Si
    @test length(std.positions) == 2
    @test std.positions[1] - std.positions[2] ≈ 0.25ones(3)
end

@testset "kgrid_from_minimal_spacing" begin
    # Test that units are stripped from both the lattice and the spacing
    lattice = [[-1.0 1 1]; [1 -1  1]; [1 1 -1]]
    @test kgrid_from_minimal_spacing(lattice * u"angstrom", 0.5 / u"angstrom") == [9; 9; 9]
end

@testset "kgrid_from_minimal_n_kpoints" begin
    lattice = [[-1.0 1 1]; [1 -1  1]; [1 1 -1]]
    @test kgrid_from_minimal_n_kpoints(lattice * u"Å", 1000) == [10, 10, 10]

    @test kgrid_from_minimal_n_kpoints(magnesium.lattice, 1) == [1, 1, 1]
    for n_kpt in [10, 20, 100, 400, 900, 1200]
        @test prod(kgrid_from_minimal_n_kpoints(magnesium.lattice, n_kpt)) ≥ n_kpt
    end

    lattice = diagm([4., 10, 0])
    @test kgrid_from_minimal_n_kpoints(lattice, 1000) == [50, 20, 1]
    @test kgrid_from_minimal_n_kpoints(diagm([10, 0, 0]), 913) == [913, 1, 1]
end


text:
"""
    occlusion_checker(veh_a, veh_b, obs::ConvexPolygon)

returns true if vehicle b is occluded by `obs` from vehicle a, false otherwise.
Occlusion checking relies on ray tracing and a the parallel axis theorem to check for collision between the ray and the obstacle.

"""
function occlusion_checker(veh_a, veh_b, obs::ConvexPolygon)
    line = @SMatrix [posg(veh_a).x posg(veh_a).y;
                     posg(veh_b).x posg(veh_b).y]
    pts = obs.pts
    obstacle = @SMatrix [pts[1].x pts[1].y;
                         pts[2].x pts[2].y;
                         pts[3].x pts[3].y;
                         pts[4].x pts[4].y]

    return AutomotiveSimulator.overlap(line, obstacle)
end

function occlusion_checker(veh_a, veh_b, veh_b_def::VehicleDef)
    ray = posg(veh_b) - posg(veh_a)
    line = @SMatrix [ray.x ray.y]
    obstacle = polygon(veh_b, veh_b_def)
    return AutomotiveSimulator.overlap(line, obstacle)
end

function occlusion_checker(veh_a, veh_b, obstacles::Vector{ConvexPolygon})
    for obs in obstacles
        if occlusion_checker(veh_a, veh_b, obs)
            return true
        end
    end
    return false
end

## For rendering
@with_kw struct OcclusionOverlay
    obstacles::Vector{ConvexPolygon} = ConvexPolygon[]
    egoid::Int64 = 1
    visible_color::Colorant = COLOR_CAR_EGO
    occluded_color::Colorant = RGB(1.0, 0., 0.)
end

function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel, overlay::OcclusionOverlay, scene::Scene, roadway::R) where R
    ego = scene[findfirst(overlay.egoid, scene)]
    #Display ray
    for veh in scene
        if veh.id == overlay.egoid
            continue
        end
        line_pts = [posg(veh.state).x posg(veh.state).y; posg(ego.state).x posg(ego.state).y]'
        if occlusion_checker(ego.state, veh.state, overlay.obstacles)
            add_instruction!( rendermodel, render_line, (line_pts, overlay.occluded_color, 0.1))
        else
            add_instruction!( rendermodel, render_line, (line_pts, overlay.visible_color, 0.1))
        end
    end
end


text:
module FeedbackTrees

using Flux
import Flux: children, mapchildren

using ..Splitters
using ..Mergers
export FeedbackTree

struct FeedbackTree{T<:Tuple}
    layers::T
    FeedbackTree(xs...) = new{typeof(xs)}(xs)
end

# empty feedback tree just returns the input
function (c::FeedbackTree{Tuple{}})(h, x)
    return h, x
end # function (c::FeedbackTree{Tuple{}})

"""
    (c::FeedbackTree)(h, x)

Apply a `FeedbackTree` to input `x` with hidden state `h`. `h` should take the
form of a dictionary mapping `Splitter` names to states.
"""
function (c::FeedbackTree)(h, x)
    newh = Dict{String, Any}()
    for layer ∈ c.layers
        if layer isa Splitter
            newh[layer.name] = x
            x = h[layer.name]
        elseif layer isa Merger
            x = layer(x, h[layer.splitname])
        else
            x = layer(x)
        end
    end
    return newh, x
end # function (c::FeedbackTree)

children(c::FeedbackTree) = c.layers
mapchildren(f, c::FeedbackTree) = FeedbackTree(f.(c.layers)...)

end # module FeedbackTrees


text:
## ------------------------------------------------------------------
# Meta

get_sec(q::RBNote) = get_parent(q)
set_sec!(q::RBNote, sec::RBSection) = set_parent!(q, sec)

## ------------------------------------------------------------------
# Data

get_text(q::RBNote) = getproperty!(q, :text, "")
set_text!(q::RBNote, text::String) = setproperty!(q, :text, strip(text))

## ------------------------------------------------------------------
# show
function Base.show(io::IO, p::RBNote)
    print(io, "RBNote(\"", get_label(p), "\")")
    print(io, "\n - path: '", localrelpath(p), "'")
    text = get_text(p)
    !isempty(text) && print(io, "\n - text: \"", _preview(io, text), "\"")
    tags = join(get_tags(p), ", ")
    !isempty(tags) && print(io, "\n - tags: \"", _preview(io, tags), "\"")
    print(io, "\n", _preview(io, "-"^70))
end


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.2
#   kernelspec:
#     display_name: Julia 1.5.0
#     language: julia
#     name: julia-1.5
# ---

using StaticArrays

"""
AbstractGeometry in R{Dim} with Number type T <: Number
"""
abstract type AbstractGeometry{Dim, T <: Number} end
abstract type GeometryPrimitive{Dim, T} <: AbstractGeometry{Dim, T} end
Base.ndims(x::AbstractGeometry{Dim}) where Dim = Dim

"""
Geometry made of N connected points. Connected as one flat geometry, it makes a Ngon/Polygon.
Connected as volume it will be a Simplex/Tri/Cube.
Note That `Polytope{N} where N == 3` denotes a Triangle both as a Simplex or Ngon.
"""
abstract type Polytope{Dim, T} <: AbstractGeometry{Dim, T} end
abstract type AbstractPolygon{Dim, T} <: Polytope{Dim, T} end

"""
Point and Face as a subtype of StaticVector
"""
abstract type AbstractPoint{Dim ,T} <: StaticVector{Dim, T} end
abstract type AbstractFace{N, T} <: StaticVector{N, T} end
abstract type AbstractSimplex{Dim, N, T} <: StaticVector{Dim, T} end
abstract type AbstractSimplexFace{N, T} <: AbstractFace{N, T} end
abstract type AbstractNgonFace{N, T} <: AbstractFace{N, T} end

"""
Face index, connecting points to form a simplex
"""

@fixed_vector Simplex AbstractSimplexFace
const TetrahedronFace{T} = SimplexFace{4, T}
Face(::Type{<: SimplexFace{N}}, ::Type{T}) where {N, T} = SimplexFace{N, T}

@fixed_vector NgonFace AbstractNgonFace
const LineFace{T} = NgonFace{2, T}
const TriangleFace{T} = NgonFace{3, T}
const QuadFace{T} = NgonFace{4,T}


text:
# ucasemap.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files ucasemap.h
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
   Unicode case mapping functions using a UCaseMap service object.

   The service object takes care of memory allocations, data loading, and setup for the attributes, as usual.

   Currently, the functionality provided here does not overlap with uchar.h and ustring.h, except for ucasemap_toTitle().
"""
module ucasemap end

macro libcasemap(s) ; _libicu(s, iculib,     "ucasemap_") ; end

const casemap  = Ptr{Cvoid}[C_NULL]
const collator = [UCollator()]

for f in (:ToLower, :ToUpper, :FoldCase, :ToTitle)
    lf = Symbol(lowercase(string(f)))
    ff = Symbol(string("utf8", f))
    uf = Symbol(string('_',lf))
    @eval begin
        ($uf)(dest, destsiz, src, srclen, err) =
            ccall(@libcasemap($ff), Int32,
                  (Ptr{Cvoid}, Ptr{UInt8}, Int32, Ptr{UInt8}, Int32, Ptr{UErrorCode}),
                  casemap[], dest, destsiz, src, srclen, err)
        function ($lf)(str::ByteStr)
            dstlen = srclen = sizeof(str)
            dest = _allocate(dstlen)
            err = Ref{UErrorCode}(0)
            siz = ($uf)(dest, dstlen, str, srclen, err)
            # Retry with large enough buffer if got buffer overflow
            if err[] == U_BUFFER_OVERFLOW_ERROR
                err[] = 0
                dstlen = siz
                dest = _allocate(dstlen)
                siz = ($uf)(dest, dstlen, str, srclen, err)
            end
            FAILURE(err[]) && error("failed to map case: $(err[])")
            siz != dstlen ? cvt_utf8(dest[1:siz]) : cvt_utf8(dest)
        end
    end
end

get_break_iterator() = ccall(@libcasemap(getBreakIterator), Ptr{Cvoid}, (Ptr{Cvoid},), casemap[])

function set_locale!(loc::ASCIIStr)
    if casemap[] != C_NULL
        ccall(@libcasemap(close), Cvoid, (Ptr{Cvoid},), casemap[])
        casemap[] = C_NULL
    end
    collator[] = UCollator(loc)
    err = Ref{UErrorCode}(0)
    casemap[] = (loc == ""
                 ? ccall(@libcasemap(open), Ptr{Cvoid}, (Ptr{UInt8}, Int32, Ptr{UErrorCode}),
                         C_NULL, 0, err)
                 : ccall(@libcasemap(open), Ptr{Cvoid}, (Cstring, Int32, Ptr{UErrorCode}),
                         loc, 0, err))
    FAILURE(err[]) && error("could not set casemap $(err[])")
    locale[] = loc
end

set_locale!(loc::AbstractString) = set_locale!(ASCIIStr(loc))


text:
#!/usr/bin/env julia

# The main purpose of this file is to read the config file and copy all of the relevant static
# files into this directory, write out the proper wavelength files, and make sure we can continue
# with the run properly.

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--test"
    help = "Print a test message to see you linked things correctly"
    action = :store_true
    "--new-config"
    help = "Copy a stock configuration file to this directory."
    action = :store_true
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
    # Options for specifying only model velocities.
    "--vstart", "-s"
    help = "Starting velocity in km/s."
    arg_type = Float64
    "--vend", "-e"
    help = "Ending velocity in km/s."
    arg_type = Float64
    "--nvel", "-n"
    help = "Number of velocity channels."
    arg_type = Int
    "--vsys"
    help = "Shift the systemic velocity"
    arg_type = Float64
end

parsed_args = parse_args(ARGS, s)

assets_dir = Pkg.dir("DiskJockey") * "/assets/"

if parsed_args["test"]
    println("Your DiskJockey scripts are successfully linked.")
    println("Exiting")
    quit()
end

if parsed_args["new-config"]
    cp(assets_dir * "config.yaml", pwd() * "/config.yaml")
    println("Copied default config.yaml file to current working directory.")
    println("Exiting")
    quit()
end

import YAML
config = YAML.load(open(parsed_args["config"]))

using HDF5

using DiskJockey.constants
using DiskJockey.hmodel

# To initialize, we basically want to do everything that would be necessary to then go in and
# run a RADMC-3D model.

# Create an output directory that will be useful for the results.
out_base = config["out_base"]
if !ispath(out_base)
    println("Creating ", out_base)
    mkdir(out_base)
end

# Are we using dust or gas?

# gas
if config["gas"]
    ff = assets_dir * "radmc3d.inp.gas"
    run(`cp $ff radmc3d.inp`)

    # Need this as a dummy file, apparently.
    ff = assets_dir * "wavelength_micron.inp"
    run(`cp $ff .`)

    # Copy appropriate molecule files.

    # Which species?
    species = config["species"]
    transition = config["transition"]
    if species == "12CO"
        ff = assets_dir * "molecule_co.inp"
        run(`cp $ff .`)
        ff = assets_dir * "lines_co.inp"
        run(`cp $ff lines.inp`)
    elseif species ==  "13CO"
        ff = assets_dir * "molecule_13co.inp"
        run(`cp $ff .`)
        ff = assets_dir * "lines_13co.inp"
        run(`cp $ff lines.inp`)
    elseif species == "C18O"
        ff = assets_dir * "molecule_c18o.inp"
        run(`cp $ff .`)
        ff = assets_dir * "lines_c18o.inp"
        run(`cp $ff lines.inp`)
    end

    println("Copied over RADMC-3D gas input files.")

    # Now run the correct setup for the grid

    # If we want just want to generate a model without fitting to data, we can specify
    # the wavelengths arbitrarily and then use `plot_model.jl`

    # Load the starting parameters
    pp = config["parameters"]
    params = ["M_star", "r_c", "T_10m", "q_m", "T_10a", "q_a", "T_freeze", "X_freeze", "gamma", "h", "delta", "logM_gas", "ksi", "dpc", "incl", "PA", "vel", "mu_RA", "mu_DEC"]

    nparam = length(params)
    starting_param = Array(Float64, nparam)

    for i=1:nparam
        starting_param[i] = pp[params[i]][1]
    end

    # Convert logM_gas to M_gas
    starting_param[12] = 10^starting_param[12]

    pars = Parameters(starting_param...)

    vel = pars.vel # [km/s]
    # RADMC conventions for inclination and PA
    incl = pars.incl # [deg]
    PA = pars.PA # [deg] Position angle runs counter clockwise
    npix = config["npix"] # number of pixels

    vstart = parsed_args["vstart"]
    vend = parsed_args["vend"]
    nvel = parsed_args["nvel"]
    if vstart != nothing
        # If we have specified the velocities, ignore the data files and Doppler shift
        # and create an evenly spaced array of velocities
        vels = linspace(vstart, vend, nvel) # [km/s]

        lam0 = lam0s[species*transition]

        # convert velocities to wavelengths
        shift_lams = lam0 * (vels/c_kms + 1)
    else
        # read the wavelengths for all data channels
        fid = h5open(config["data_file"], "r")
        lams = read(fid["lams"]) # [μm]
        close(fid)

        # Doppler shift the dataset wavelength according to the velocity in the parameter file
        beta = vel/c_kms # relativistic Doppler formula
        shift_lams =  lams .* sqrt((1. - beta) / (1. + beta)) # [microns]
    end

    grd = config["grid"]
    grid = Grid(grd["nr"], grd["ntheta"], grd["r_in"], grd["r_out"], true)

    # Write everything to the current working directory
    write_grid("", grid)
    tic()
    write_model(pars, "", grid, species)
    println("Model writing time")
    toc()
    write_lambda(shift_lams, "")

    println("Wrote gas model input files for RADMC-3D.")

# dust
else
    ff = assets_dir * "radmc3d.inp.dust"
    run(`cp $ff radmc3d.inp`)
    ff = assets_dir * "stars.inp"
    run(`cp $ff .`)
    ff = assets_dir * "wavelength_micron.inp"
    run(`cp $ff .`)
    ff = assets_dir * "dustkappa_silicate.inp"
    run(`cp $ff .`)
    ff = assets_dir * "dustopac.inp"
    run(`cp $ff .`)

    println("Copied over RADMC-3D dust input files.")
end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
function exampleCategorical()
    nsteps = 300#10000#300
    nseeds = 3

    specs = SpecsMultinomial(stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01],#[1., 0.25, 0.14, 0.1, 0.01],
                           changeprobabilities = [0.1, 0.05, 0.01, 0.005, 0.001, 0.0001],
                           nseeds = nseeds, seedenvpolicy = collect(1:nseeds),
                           seedlearner = collect(1:nseeds) * 99,
                           ns = 5, nsteps = nsteps)
    iseed = 1
    is = 1
    ich = 2
    learnertype = TParticleFilter
    learnerparam = (20, specs.stochasticities[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = TMPN
    # learnerparam = (20, specs.stochasticities[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = TVarSmile
    # learnerparam = (0.01, specs.stochasticities[is])
    # ---------------
    # learnertype = TBayesFilter
    # learnerparam = (specs.stochasticities[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = TSOR
    # learnerparam = (20, specs.stochasticities[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = TSmileOriginal
    # learnerparam = (0.01, specs.stochasticities[is])
    # ---------------
    callback = run!(specs.envtype, learnertype, learnerparam,
                    ns = specs.ns,
                    stochasticity = specs.stochasticities[is],
                    changeprobability = specs.changeprobabilities[ich],
                    nsteps = nsteps,
                    iseedenvpolicy = specs.seedenvpolicy[iseed],
                    iseedlearner = specs.seedlearner[iseed]
                    )
    terror = computesquareddiff(callback.true_history,
                                callback.estimated_history)
    fig = figure(); ax = gca()
    ax.plot(terror)
    iswitch = findall(callback.switchflag)
    if !isempty(iswitch)
        for i in 1:length(iswitch)
            ax.axvline(iswitch[i], color="red", linestyle="--", linewidth=0.5)
        end
    end
    ax.grid(true)
    ax.set_ylabel("average error")
    ax.set_xlabel("steps")

    switches = callback.switchflag
    errswitches = geterrorsafterswitches(terror, switches,
                    Int(1/specs.changeprobabilities[ich]))
    errswitches = hcat(errswitches...)
    merrswitches, stderrswitches, stderrorrerrswitches = calc_stats(errswitches)
    fig = figure(); ax = gca()
    ax.plot(merrswitches, color = "r")
    ax.grid(true)
    ax.set_ylim(0., 0.1);
    ax.set_ylabel("average error")
    ax.set_xlabel("steps after change")
    callback
end

function exampleGaussian()
    nsteps = 300
    nseeds = 3
    specs = SpecsGaussian(; mu_0 = 0., sigma_0 = 1.,
                    sigmas = [5., 2., 1., 0.5, 0.1],
                    changeprobabilities = [0.1, 0.05, 0.01, 0.005, 0.001, 0.0001],
                    nseeds = nseeds, seedenvpolicy = collect(1:nseeds),
                    seedlearner =  collect(1:nseeds) * 99,
                    nsteps = nsteps)
    iseed = 1
    is = 2
    ich = 1
    # ---------------
    learnertype = GParticleFilter
    learnerparam = (20, specs.sigmas[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = GVarSmile
    # learnerparam = (0.1, specs.sigmas[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = GMPN
    # learnerparam = (1, specs.sigmas[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = GNas10
    # learnerparam = (specs.sigmas[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = GNas12
    # learnerparam = (specs.sigmas[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = GBayesFilter
    # learnerparam = (specs.sigmas[is], specs.changeprobabilities[ich])
    # ---------------
    # learnertype = GSmileOriginal
    # learnerparam = (0.1, specs.sigmas[is])
    # ---------------
    # learnertype = GSOR
    # learnerparam = (20, specs.sigmas[is], specs.changeprobabilities[ich])
    # ---------------
    @show learnertype

    callback = run!(specs.envtype, learnertype, learnerparam,
                    sigma = specs.sigmas[is],
                    changeprobability = specs.changeprobabilities[ich],
                    nsteps = nsteps,
                    iseedenvpolicy = specs.seedenvpolicy[iseed],
                    iseedlearner = specs.seedlearner[iseed]
                    )

    upTo = nsteps#150
    fig = figure(); ax = gca()
    ax.plot(callback.true_history[2:upTo], label="true")
    ax.plot(callback.estimated_history[2:upTo], color = "r", label = "estimation")
    ax.scatter(0:(upTo-2), callback.state_history[2:upTo], color = "g", s = 4)
    #
    ax.plot(log10.(callback.Sgm_history[2:upTo]), color = "magenta", label = "Bayes Factor Surprise")
    # ax.plot(callback.Py_0_history[2:upTo], color = "black", )
    # ax.plot(log10.(callback.var_history[2:upTo]), color = "black", linestyle = "--")
    # ax.plot(log10.(callback.Ssh_history[2:upTo]), color = "black", linestyle = "--")
    # ax.plot(callback.γ_history[2:upTo], color = "black", linestyle = "--")

    ax.grid(true, which="major")
    ax.minorticks_on()
    ax.grid(true, which="minor")
    ax,legend()

    merror = computesquareddiff(callback.true_history[2:end],
                                callback.estimated_history[2:end])

    switches = callback.switchflag
    errswitches = geterrorsafterswitches(merror, switches,
                    # Int(1/specs.changeprobabilities[ich])) #ERROR: InexactError: Int64(24.244620170823282)
                    Int(round(1/specs.changeprobabilities[ich])))
    errswitches = hcat(errswitches...)
    merrswitches, stderrswitches, stderrorrerrswitches = calc_stats(errswitches)
    # fig = figure(); ax = gca()
    # ax.plot(merrswitches, color = "r")
    # ax.grid(true)
    # ax.set_ylim(0., 0.014);
    # plot_shadedline_core(merrswitches, stderrorrerrswitches, "error after switches", "b", ax)

    callback, merror
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 38144072-3fab-478c-aa25-947c7086d7d9
begin
	using Pkg
	using Revise
	Pkg.activate("..")
	using LRHousing
end

# ╔═╡ c060061b-1b8d-45f6-98c0-c4d9e28c16f7
begin
	struct DataLine
		idx_numeric
		log_numeric
		categorical
		ignore
		function DataLine(df, ignore)
			ignore = @chain begin
				names(df)
				filter(name -> (count(ismissing, df[!, name]) / length(df[!, name])) ≥ 0.5, _)
				vcat(_, ignore)
				Symbol.(_)
			end
			allnames = @chain begin
				names(df)
				Symbol.(_)
				filter(!in(ignore), _)
				filter(!isequal(ignore), _)
			end
			numeric = @chain begin
				allnames
				filter(name -> Bool != eltype(df[!, name]), _)
				filter(name -> eltype(df[!, name]) <: Number, _)
			end
			log_numeric = filter(name -> skewness(df[!, name]) ≥ 0.75, numeric)
			categorical = filter(!in(numeric), allnames)
			new(numeric, log_numeric, categorical, ignore)
		end
	end

	function (l::DataLine)(df)
		df = select(df, Not(l.ignore))
		function znorm(X)
			X = replace(X, missing => mean(filter(!ismissing, X)))
			standardize(ZScoreTransform, convert.(Float32, X))
		end
		function categorical_(X)
			X = replace(X, missing => mode(filter(!ismissing, X)))
			categorical(X)
		end
		df = @chain begin
			l.idx_numeric
			map(name -> (name => znorm => name), _)
			transform(df, _...)
		end
		df = @chain begin
			l.categorical
			map(name -> (name => categorical_ => name), _)
			transform(df, _...)
		end
		return df
	end
end

# ╔═╡ 34a4eb8d-b724-4f28-9d5c-159a3c2b1a6d
function read_csv(file)
	df = @chain begin
		CSV.read(file, DataFrame, missingstrings=["NA"], normalizenames=true)
# 		map(names(_)) do name
# 			col = _[!, name]
# 			col = if eltype(col) <: AbstractString
# 				categorical(col)
# 			else
# 				col
# 			end
# 			name => col
# 		end
# 		DataFrame(_)
# 	end
	end
end

# ╔═╡ 91eaf5c4-909f-4d62-877c-2f807a872fc0
train_data = read_csv("../data/train.csv")

# ╔═╡ ea0bab7a-6260-4995-a775-857f66651da3
preprocess = DataLine(train_data, ["Id"])

# ╔═╡ 51afcb1f-f749-496c-acc5-4cecb4e4d1df
preprocess(train_data)

# ╔═╡ 70d6defb-8557-43e1-8a9a-639640d2dd67
preprocess(train_data)

# ╔═╡ 5b242981-a0ef-4638-a3b3-d35617686744
test_data = read_csv("../data/test.csv")

# ╔═╡ 15c9cc08-3834-4922-9789-52023a9a2489
function single_functional_term(f, x)
	symbol = Expr(:call, Symbol(f), x)
	fanon = eval(:(($x,) -> log($x)))
	StatsModels.capture_call(f, fanon, (x,), symbol, [Term(x)])
end

# ╔═╡ 5b3a0677-9c64-4669-b90e-8ad7286b0531
e = :(f(x))

# ╔═╡ 5718f7ce-328a-4705-85ad-052494be8e93
predictor_names = setdiff(names(train_data), ["Id", "SalePrice"])

# ╔═╡ 25981a26-7762-4daa-8b69-936a027a88dc
numeric_columns = @chain begin
	predictor_names
	filter(x -> eltype(train_data[!, x]) <: Number, _)
end

# ╔═╡ 8fe249bd-5102-4808-baf7-e769e4b1056a
FormulaTerm(Term(:y), (Term(:x),))

# ╔═╡ 0a9b7ae3-bd8f-4b12-95a9-c132715d7f1b
f = let res = Term(:SalePrice)
	pres = map(numeric_columns) do name
		name = Symbol(name)
		if skewness(train_data[!, name]) > 0.75
			single_functional_term(log, name)
		else
			Term(name)
		end
	end
	FormulaTerm(res, (pres...,))
end

# ╔═╡ 65ed182a-fea7-472d-86c9-e51958ebf6d2
lm(f, train_data)

# ╔═╡ febc052d-82f2-455e-b35b-3fc230899657
ex = Expr(:call, :~, :y, Expr(:call, :+, 1, :x1, :x2))

# ╔═╡ d96899e6-3484-4b0e-b902-b1c709d9c6f1
dump(ex)

# ╔═╡ 4ceaa564-ee7a-47c6-ae2f-efe0f94f968b
ex2 = :(y ~ 1 + x + x2)

# ╔═╡ 424a6f2d-4aa0-46fb-b2a6-527cbedef366
function build_formular(y, xs)
	e = Expr(:call, :~, Symbol(y), Expr(:call, :+, 1, Symbol.(xs)...))
	eval(:(@formula $e))
end

# ╔═╡ 37d737c9-6be9-40f0-8434-53eae3d5936b
f2 = build_formular(:SalePrice, numeric_columns)

# ╔═╡ f3e33cfe-4a73-4033-942e-2dc0b85fb6d7
lm(f2, train_data)

# ╔═╡ ab1753c5-8672-45e2-ac68-f5bf1d6374fc
eval(:(@formula $ex))

# ╔═╡ d3e5677e-a491-4c58-9f7a-2299339fc299
InteractionTerm((Term(:x), Term(:y))

# ╔═╡ e8f21844-31b1-45a0-a4cd-a36c32a9292e
let i = :a
	StatsModels.capture_call(log, (a,)->log(a), (:a,), :(log(a)), [Term(:a)])
	# [StatsModels.Term($i)]
	# )
	# ~ Term(:AX))
	# end
end

# ╔═╡ Cell order:
# ╠═38144072-3fab-478c-aa25-947c7086d7d9
# ╠═c060061b-1b8d-45f6-98c0-c4d9e28c16f7
# ╠═ea0bab7a-6260-4995-a775-857f66651da3
# ╠═51afcb1f-f749-496c-acc5-4cecb4e4d1df
# ╠═70d6defb-8557-43e1-8a9a-639640d2dd67
# ╠═34a4eb8d-b724-4f28-9d5c-159a3c2b1a6d
# ╠═91eaf5c4-909f-4d62-877c-2f807a872fc0
# ╠═5b242981-a0ef-4638-a3b3-d35617686744
# ╠═15c9cc08-3834-4922-9789-52023a9a2489
# ╠═5b3a0677-9c64-4669-b90e-8ad7286b0531
# ╠═5718f7ce-328a-4705-85ad-052494be8e93
# ╠═25981a26-7762-4daa-8b69-936a027a88dc
# ╠═8fe249bd-5102-4808-baf7-e769e4b1056a
# ╠═0a9b7ae3-bd8f-4b12-95a9-c132715d7f1b
# ╠═65ed182a-fea7-472d-86c9-e51958ebf6d2
# ╠═febc052d-82f2-455e-b35b-3fc230899657
# ╠═d96899e6-3484-4b0e-b902-b1c709d9c6f1
# ╠═4ceaa564-ee7a-47c6-ae2f-efe0f94f968b
# ╠═424a6f2d-4aa0-46fb-b2a6-527cbedef366
# ╠═37d737c9-6be9-40f0-8434-53eae3d5936b
# ╠═f3e33cfe-4a73-4033-942e-2dc0b85fb6d7
# ╠═ab1753c5-8672-45e2-ac68-f5bf1d6374fc
# ╠═d3e5677e-a491-4c58-9f7a-2299339fc299
# ╠═e8f21844-31b1-45a0-a4cd-a36c32a9292e


text:
@doc raw"""
    gradient_descent(M, F, ∇F, x)

perform a gradient_descent ``x_{k+1} = \mathrm{retr}_{x_k} s_k\nabla f(x_k)`` with
different choices of ``s_k`` available (see `stepsize` option below).

# Input
* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F\colon\mathcal M\to\mathbb R`` to minimize
* `∇F` – the gradient ``\nabla F\colon\mathcal M\to T\mathcal M`` of F
* `x` – an initial value ``x ∈ \mathcal M``

# Optional
* `stepsize` – ([`ConstantStepsize`](@ref)`(1.)`) specify a [`Stepsize`](@ref)
  functor.
* `retraction_method` – (`ExponentialRetraction()`) a `retraction(M,x,ξ)` to use.
* `stopping_criterion` – ([`StopWhenAny`](@ref)`(`[`StopAfterIteration`](@ref)`(200), `[`StopWhenGradientNormLess`](@ref)`(10.0^-8))`)
  a functor inheriting from [`StoppingCriterion`](@ref) indicating when to stop.
* `return_options` – (`false`) – if activated, the extended result, i.e. the
    complete [`Options`](@ref) are returned. This can be used to access recorded values.
    If set to false (default) just the optimal value `x_opt` if returned
...
and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output
* `x_opt` – the resulting (approximately critical) point of gradientDescent
OR
* `options` - the options returned by the solver (see `return_options`)
"""
function gradient_descent(M::Manifold, F::TF, ∇F::TDF, x; kwargs...) where {TF,TDF}
    x_res = allocate(x)
    copyto!(x_res, x)
    return gradient_descent!(M, F, ∇F, x_res; kwargs...)
end
@doc raw"""
    gradient_descent!(M, F, ∇F, x)

perform a gradient_descent ``x_{k+1} = \mathrm{retr}_{x_k} s_k\nabla f(x_k)`` inplace of `x`
with different choices of ``s_k`` available.

# Input
* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F\colon\mathcal M\to\mathbb R`` to minimize
* `∇F` – the gradient ``\nabla F\colon\mathcal M\to T\mathcal M`` of F
* `x` – an initial value ``x ∈ \mathcal M``

For more options, especially [`Stepsize`](@ref)s for ``s_k``, see [`gradient_descent`](@ref)
"""
function gradient_descent!(
    M::Manifold,
    F::TF,
    ∇F::TDF,
    x;
    stepsize::Stepsize=ConstantStepsize(1.0),
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    stopping_criterion::StoppingCriterion=StopWhenAny(
        StopAfterIteration(200), StopWhenGradientNormLess(10.0^-8)
    ),
    direction=IdentityUpdateRule(),
    return_options=false,
    kwargs..., #collect rest
) where {TF,TDF}
    p = GradientProblem(M, F, ∇F)
    o = GradientDescentOptions(
        x;
        stopping_criterion=stopping_criterion,
        stepsize=stepsize,
        direction=direction,
        retraction_method=retraction_method,
    )
    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
#
# Solver functions
#
function initialize_solver!(p::P, o::O) where {P<:GradientProblem,O<:GradientDescentOptions}
    o.∇ = get_gradient(p, o.x)
    return o
end
function step_solver!(p::P, o::O, iter) where {P<:GradientProblem,O<:GradientDescentOptions}
    s, o.∇ = o.direction(p, o, iter)
    o.x = retract(p.M, o.x, -s .* o.∇, o.retraction_method)
    return o
end
get_solver_result(o::GradientDescentOptions) = o.x


text:
module OpenStreetMapPlotter
using LightXML
using HTTP
using Winston

include("structs.jl")
include("styles.jl")
include("MapCSSParser.jl")
function open_file(filepath::String)
	xdoc = parse_file(filepath)
	xroot = root(xdoc)  # an instance of XMLElement
	return xroot
end

function open_bbox(bbox::Tuple)
	minlon = bbox[1]
	maxlon = bbox[3]
	minlat = bbox[2]
	maxlat = bbox[4]
	url = "http://overpass-api.de/api/map?bbox=$(minlon),$(minlat),$(maxlon),$(maxlat)"
	r = HTTP.request("GET", url)
	return root(parse_string(String(r.body)))
end

function parse_nodes(xroot::XMLElement)
    node_arr = Node[]
    for node in xroot["node"]
        id = attribute(node, "id")
        lat = parse(Float64, attribute(node, "lat"))
        lon = parse(Float64, attribute(node, "lon"))
        push!(node_arr, Node(lon, lat, id))
    end
    return node_arr
end

function find_node(id::String, node_arr::Array{Node})
    return node_arr[findfirst(x -> x.id == id, node_arr)]
end

function find_way(id::String, way_arr::Array{Way})
	try
    	return way_arr[findfirst(x -> x.id == id, way_arr)]
    catch
    	return false
    end
end
function parse_ways(xroot::XMLElement)
    way_arr = Way[]
    node_arr = parse_nodes(xroot)
    for way in xroot["way"]
        cur_way = Way()
        cur_way.id = attribute(way, "id")
        cur_way.visible = attribute(way, "visible") == "true"
        cur_way.version = parse(Int, attribute(way, "version"))
        cur_way.changeset = attribute(way, "changeset")
        cur_way.timestamp = attribute(way, "timestamp")
        cur_way.user = attribute(way, "user")
        cur_way.uid = attribute(way, "uid")
        cur_way.nodes = []
        for node in way["nd"]
            push!(cur_way.nodes, find_node(attribute(node, "ref"), node_arr))
        end
        cur_way.tags = Dict()
        for tag in way["tag"]
        	cur_way.tags[attribute(tag, "k")] = attribute(tag, "v")
        end
        push!(way_arr, cur_way)
    end
    bounds = xroot["bounds"][1]
    bbox = (parse(Float64, attribute(bounds, "minlon")), parse(Float64, attribute(bounds, "minlat")), parse(Float64, attribute(bounds, "maxlon")), parse(Float64, attribute(bounds, "maxlat")))
    return way_arr, bbox
end

function parse_relations(xroot::XMLElement, way_arr::Array{Way}, node_arr::Array{Node})
	rel_arr = []
    for rel in xroot["relation"]
        cur_rel = Relation()
        cur_rel.id = attribute(rel, "id")
        cur_rel.visible = lowercase(attribute(rel, "visible")) == "true"
        cur_rel.version = parse(Int, attribute(rel, "version"))
        cur_rel.changeset = attribute(rel, "changeset")
        cur_rel.timestamp = attribute(rel, "timestamp")
        cur_rel.user = attribute(rel, "user")
        cur_rel.uid = attribute(rel, "uid")
        cur_rel.members = []
        for member in rel["member"]
        	ref = attribute(member, "ref")
        	if attribute(member, "type") == "node"
            	push!(cur_rel.members, find_node(ref, node_arr))
            elseif attribute(member, "type") == "way"
            	try
            		temp_way = find_way(ref, way_arr)
            		if temp_way != false
            			push!(cur_rel.members, temp_way)
            		end
            	catch
            		continue
            	end
            end
        end
        cur_rel.tags = Dict()
        for tag in rel["tag"]
            cur_rel.tags[attribute(tag, "k")] = attribute(tag, "v")
        end
        push!(rel_arr, cur_rel)
    end
    return rel_arr
end

function plot_ways(way_arr::Array{Way}, bbox::Tuple; width::Int64=900, css_file_name::String="")
	minlon = bbox[1]
	maxlon = bbox[3]
	minlat = bbox[2]
	maxlat = bbox[4]
	c_adj = cosd((minlat + maxlat) / 2)
    range_y = maxlat - minlat
    range_x = maxlon - minlon
	aspect_ratio = range_x * c_adj / range_y
	fignum = Winston.figure(name="OpenStreetMap Plot", width=width, height=round(Int, width/aspect_ratio))
	p = FramedPlot(xrange = (minlon, maxlon), yrange = (minlat, maxlat))
	p.attr[:gutter] = 0
	p.attr[:title_style][:fontsize] = 1
	p.x1.attr[:ticklabels_style][:fontsize] = 1
	p.y1.attr[:ticklabels_style][:fontsize] = 1
	p.x1.attr[:draw_subticks] = true
	p.y1.attr[:draw_subticks] = true
	p.x1.attr[:ticks] = 10
	p.y1.attr[:ticks] = 10
	Winston._winston_config.defaults["fontsize_min"] = ".1"
	layers = Vector{Way}[[],[],[],[],[],[],[],[],[],[],[]]
	draw_later = []
	cascade = []
	if css_file_name != ""
		cascade = parse_css(css_file_name)
	end
	for way in way_arr
		if haskey(way.tags, "layer")
			l = tryparse(Int, way.tags["layer"])
			if l == nothing
				println("invalid layer tag")
			else
				push!(layers[l+6], way)
			end
		end
		if haskey(way.tags, "highway")
				if way.tags["highway"] in ["motorway", "trunk", "primary", "secondary", "tertiary"]
					push!(layers[8], way)
				else
					push!(layers[7], way)
				end
		else
			push!(layers[6], way)
		end
	end
	labels = []
	for layer in layers
		for way in layer
			style=get_way_style(way.tags, cascade)
			if style == nothing
				continue
			end
			if way.nodes[1] == way.nodes[end]
				style["is_polygon"] = true
			else
				style["is_polygon"] = false
			end
			if haskey(style, "text")
				for i in keys(way.tags)
					if occursin(i, style["text"])
						c = center_of_points(way.nodes)
						c = ((c[1]-minlon)/range_x, (c[2]-minlat)/range_y)
						cur_label = PlotLabel(c[1], c[2], replace(way.tags[i], "&"=>"&"), fontsize = style["font-size"])
						push!(labels, cur_label)
						break
					end
				end
			end
			f = nothing
			if style["is_polygon"] == true
				split = split_polygon(way.nodes)
				topside = way.nodes[1:split]
				bottomside = way.nodes[split:end]
				f = FillBetween([i.x for i in topside], [i.y for i in topside], [i.x for i in bottomside], [i.y for i in bottomside], fillcolor = style["color"], linewidth=style["width"])
			else
				f = Curve([i.x for i in way.nodes], [i.y for i in way.nodes], color=style["color"], linewidth=style["width"])
			end
	    	if f != nothing
	    		Winston.add(p, f)
	    	end
		end
	end
	println(labels)
	for label in labels
		Winston.add(p, label)
	end
    display(p)
end
function split_polygon(nodes::Array{Node})
	return findmax([i.x for i in nodes])[2]
end
function sort_counterclockwise(nodes::Array{Node})
	center = center_of_points(nodes)
	return sort(nodes, lt=(a, b)->!is_less(a, b, center))
end
#taken from here: https://rosettacode.org/wiki/Shoelace_formula_for_polygonal_area#Julia
shoelacearea(x, y) = abs(sum(i * j for (i, j) in zip(x, append!(y[2:end], y[1]))) - sum(i * j for (i, j) in zip(append!(x[2:end], x[1]), y))) / 2
function get_area(nodes::Array{Node})
	x = [i.x for i in nodes]
	y = [i.y for i in nodes]
	return shoelacearea(x, y)
end
function is_less(a, b, center)
	if a.x-center[1] >= 0 && b.x-center[1] < 0
		return true
	elseif a.x -center[1] == 0 && b.x-center[1] == 0
		return a.y >b.y
	end
	det = (a.x - center[1]) * (b.y - center[2]) - (b.x - center[2]) * (a.y - center[2])
	if det < 0
		return true
	elseif det > 0
		return false
	end
	d1 = (a.x - center[1]) * (a.x - center[1]) + (a.y - center[2]) * (a.y - center[2])
	d2 = (b.x - center[1]) * (b.x - center[1]) + (b.y - center[2]) * (b.y - center[2])
	return d1 > d2
end
function mean(array::AbstractArray)
	return sum(array)/length(array)
end
function center_of_points(nodes::Array{Node})
	return (mean([i.x for i in nodes]), mean([i.y for i in nodes]))
end
function save_json(way_arr::Array{Way}, filepath::String)
	f=open(filepath, "w")
	write(f, "{\n\"type\": \"FeatureCollection\",\n\"features\": [\n")
	flags = [true, true, true]
	for way in way_arr
		if flags[1] == true
			flags[1]=false
		else
			write(f, ",\n")
		end
		write(f, "{ \"type\": \"Feature\", \"properties\": {")
		flags[1] = false
		flags[2] = true
		for key in keys(way.tags)
			key = key
			val = way.tags[key]
			if flags[2] == true
				flags[2] = false
			else
				write(f, ",")
			end
			if in('"', val)
				write(f, "\"$key\":\"$(replace(val, r"\"" => "\\\""))\"")
			else
				write(f, "\"$key\":\"$val\"")
			end
		end
		startarraystr = ""
		endarraystr = ""
		sorted_arr = []
		if way.nodes[1] == way.nodes[end]
			write(f, "}, \"geometry\": { \"type\": \"Polygon\", \"coordinates\":")
			startarraystr = "[ ["
			endarraystr = "]] }}\n"
			sorted_arr = sort_counterclockwise(way.nodes)
			sorted_arr = vcat(sorted_arr, sorted_arr[1])
		elseif length(way.nodes) == 1
			write(f, "}, \"geometry\": { \"type\": \"Point\", \"coordinates\":")
			startarraystr = ""
			endarraystr = "}}\n"
			sorted_arr = way.nodes
		else
			write(f, "}, \"geometry\": { \"type\": \"LineString\", \"coordinates\":")
			startarraystr = "["
			endarraystr = "] }}\n"
			sorted_arr = sort_counterclockwise(way.nodes)
		end
		write(f, startarraystr)
		flags[3]=true

		for coord in sorted_arr
			lon = coord.x
			lat = coord.y
			if flags[3] == true
				flags[3] = false
			else
				write(f, ", ")
			end
			flags[3] = false
			write(f, "[ $lon, $lat]")
		end
		write(f, endarraystr)
	end

	write(f, "\n]\n}")
	close(f)
end
export open_file, open_bbox, parse_nodes, parse_ways, parse_relations, plot_ways, save_json, Node, Tag, Way, highway_styles, building_styles, waterway_styles

end


text:
import DelayEmbeddings: genembed, Dataset
import Statistics: mean

export SymbolicWeightedPermutation

"""
    SymbolicWeightedPermutation(; τ = 1, m = 3, lt = Entropies.isless_rand) <: PermutationProbabilityEstimator

See docstring for [`SymbolicPermutation`](@ref).
"""
struct SymbolicWeightedPermutation{F} <: PermutationProbabilityEstimator
    τ::Int
    m::Int
    lt::F
end
function SymbolicWeightedPermutation(; τ::Int = 1, m::Int = 3, lt::F = isless_rand) where {F <: Function}
    m >= 2 || error("Need m ≥ 2, otherwise no dynamical information is encoded in the symbols.")
    SymbolicWeightedPermutation{F}(τ, m, lt)
end

function weights_from_variance(x, m::Int)
    sum((x .- mean(x)) .^ 2)/m
end

function probabilities(x::AbstractDataset{m, T}, est::SymbolicWeightedPermutation) where {m, T}
    m >= 2 || error("Need m ≥ 2, otherwise no dynamical information is encoded in the symbols.")
    πs = symbolize(x, SymbolicPermutation(m = m, lt = est.lt))  # motif length controlled by dimension of input data
    wts = weights_from_variance.(x.data, m)

    Probabilities(symprobs(πs, wts, normalize = true))
end

function probabilities(x::AbstractVector{T}, est::SymbolicWeightedPermutation) where {T<:Real}
    τs = tuple([est.τ*i for i = 0:est.m-1]...)
    emb = genembed(x, τs)
    πs = symbolize(emb, SymbolicPermutation(m = est.m, lt = est.lt)) # motif length controlled by estimator m
    wts = weights_from_variance.(emb.data, est.m)

    Probabilities(symprobs(πs, wts, normalize = true))
end


text:
export BrickletAmbientLightIlluminanceCallbackThreshold
struct BrickletAmbientLightIlluminanceCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletAmbientLightAnalogValueCallbackThreshold
struct BrickletAmbientLightAnalogValueCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletAmbientLightIdentity
struct BrickletAmbientLightIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAmbientLight
"""
Measures ambient light up to 900lux
"""
mutable struct BrickletAmbientLight <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAmbientLight(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_ambient_light")
        deviceInternal = package.BrickletAmbientLight(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_illuminance
"""
    $(SIGNATURES)

Returns the illuminance of the ambient light sensor.

If you want to get the illuminance periodically, it is recommended to use the
:cb:`Illuminance` callback and set the period with
:func:`Set Illuminance Callback Period`.
"""
function get_illuminance(device::BrickletAmbientLight)
    return device.deviceInternal.get_illuminance()
end

export get_analog_value
"""
    $(SIGNATURES)

Returns the value as read by a 12-bit analog-to-digital converter.

.. note::
 The value returned by :func:`Get Illuminance` is averaged over several samples
 to yield less noise, while :func:`Get Analog Value` gives back raw
 unfiltered analog values. The only reason to use :func:`Get Analog Value` is,
 if you need the full resolution of the analog-to-digital converter.

 Also, the analog-to-digital converter covers three different ranges that are
 set dynamically depending on the light intensity. It is impossible to
 distinguish between these ranges with the analog value.

If you want the analog value periodically, it is recommended to use the
:cb:`Analog Value` callback and set the period with
:func:`Set Analog Value Callback Period`.
"""
function get_analog_value(device::BrickletAmbientLight)
    return device.deviceInternal.get_analog_value()
end

export set_illuminance_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Illuminance` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Illuminance` callback is only triggered if the illuminance has changed
since the last triggering.
"""
function set_illuminance_callback_period(device::BrickletAmbientLight, period)
    device.deviceInternal.set_illuminance_callback_period(period)
end

export get_illuminance_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Illuminance Callback Period`.
"""
function get_illuminance_callback_period(device::BrickletAmbientLight)
    return device.deviceInternal.get_illuminance_callback_period()
end

export set_analog_value_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Analog Value` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Analog Value` callback is only triggered if the analog value has
changed since the last triggering.
"""
function set_analog_value_callback_period(device::BrickletAmbientLight, period)
    device.deviceInternal.set_analog_value_callback_period(period)
end

export get_analog_value_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Analog Value Callback Period`.
"""
function get_analog_value_callback_period(device::BrickletAmbientLight)
    return device.deviceInternal.get_analog_value_callback_period()
end

export set_illuminance_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Illuminance Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the illuminance is *outside* the min and max values"
 "'i'",    "Callback is triggered when the illuminance is *inside* the min and max values"
 "'<'",    "Callback is triggered when the illuminance is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the illuminance is greater than the min value (max is ignored)"
"""
function set_illuminance_callback_threshold(device::BrickletAmbientLight, option, min, max)
    device.deviceInternal.set_illuminance_callback_threshold(option, min, max)
end

export get_illuminance_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Illuminance Callback Threshold`.
"""
function get_illuminance_callback_threshold(device::BrickletAmbientLight)
    return device.deviceInternal.get_illuminance_callback_threshold()
end

export set_analog_value_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Analog Value Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the analog value is *outside* the min and max values"
 "'i'",    "Callback is triggered when the analog value is *inside* the min and max values"
 "'<'",    "Callback is triggered when the analog value is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the analog value is greater than the min value (max is ignored)"
"""
function set_analog_value_callback_threshold(device::BrickletAmbientLight, option, min, max)
    device.deviceInternal.set_analog_value_callback_threshold(option, min, max)
end

export get_analog_value_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Analog Value Callback Threshold`.
"""
function get_analog_value_callback_threshold(device::BrickletAmbientLight)
    return device.deviceInternal.get_analog_value_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Illuminance Reached`,
* :cb:`Analog Value Reached`

are triggered, if the thresholds

* :func:`Set Illuminance Callback Threshold`,
* :func:`Set Analog Value Callback Threshold`

keep being reached.
"""
function set_debounce_period(device::BrickletAmbientLight, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletAmbientLight)
    return device.deviceInternal.get_debounce_period()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAmbientLight)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletAmbientLight, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
# Randomized interpolative decomposition

using LowRankApprox

function rid(X,
             k::Int,
             p::Int=10,
             q::Int=0)

    Q, B = rqb(X, k, p, q)
    Z, j = id(B, k)
    C = A[:, j]

    return A[:, j], Z, j
end


text:
module Simulate_DoublePendulumWithDampers

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
# massProperties = Modia3D.MassProperties()
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

@forceElement Damper(; d=1.0) begin
   w   = ModiaMath.RealScalar("w",   causality=ModiaMath.Input,  numericType=ModiaMath.WR)
   tau = ModiaMath.RealScalar("tau",  causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeTorque(damper::Damper, sim::ModiaMath.SimulationState)
   damper.tau.value = -damper.d*damper.w.value
end

@assembly Bar(;Lx = 0.1, Ly=Lx/5, Lz=Ly, m=1.0) begin
   frame0 = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
   frame1 = Modia3D.Object3D(frame0; r=[-Lx/2, 0.0, 0.0])
   frame2 = Modia3D.Object3D(frame0; r=[ Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2)
   cyl1   = Modia3D.Object3D(frame1, cyl)
   cyl2   = Modia3D.Object3D(frame2, cyl)
end

@assembly DoublePendulumWithDampers(;Lx = 1.0, m=1.0) begin
   world = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))
   bar1  = Bar(Lx=Lx, m=m)
   bar2  = Bar(Lx=Lx, m=m)
   rev1  = Modia3D.Revolute(world, bar1.frame1)
   rev2  = Modia3D.Revolute(bar1.frame2, bar2.frame1)
   d1     = Damper(d=0.2)
   damper1 = Modia3D.AdaptorForceElementToFlange(w=d1.w, tau=d1.tau)
   Modia3D.connect(damper1, rev1)
   d2     = Damper(d=0.3)
   damper2 = Modia3D.AdaptorForceElementToFlange(w=d2.w, tau=d2.tau)
   Modia3D.connect(damper2, rev2)
end

model = Modia3D.SimulationModel( DoublePendulumWithDampers())
result = ModiaMath.simulate!(model, stopTime=5.0, tolerance=1e-6,interval=0.001,log=true)

ModiaMath.plot(result, [("rev1.phi", "rev2.phi"),
                        ("rev1.w"  , "rev2.w"),
                        ("rev1.a"  , "rev2.a")])

println("... success of Simulate_DoublePendulumWithDampers.jl!")
end


text:
using BinaryBuilder

include("../common.jl")

name = "Glibc"
version = v"2.19"

sources = glibc_sources(version)
script = glibc_script()
platforms = glibc_platforms(version)
products = glibc_products()
dependencies = glibc_dependencies()

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", lock_microarchitecture=false)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module datagen

using Random, Distributions
using OrdinaryDiffEq, Flux

export generate_data

clip(x, vmin, vmax) = max(min(x, vmax), vmin)

################################################################################
##                                                                            ##
##          Drawing images of (filled) circles on 32×32 canvases              ##
##                                                                            ##
################################################################################

# Generate quantized circles via Monte Carlo
circs = Matrix{Matrix}(undef, 10, 10)
@info "Pre-calculating circle discretization via Monte Carlo... "
for xx in -0.5:0.1:0.4
    for yy in -0.5:0.1:0.4
        _circmat = zeros(Int, 6, 6)
        _circmat_n = zeros(Int, 6, 6)
        for i in -3:0.01:3-0.01, j in -3:0.01:3-0.01
            _circmat_n[floor(Int, i)+4, floor(Int, j)+4] += 1
            (((i-xx)^2 + (j-yy)^2) <= 1.6^2) &&
                (_circmat[floor(Int, i)+4, floor(Int, j)+4] += 1)
        end
        _circmat = _circmat ./ _circmat_n;
        circs[round(Int, xx/0.1) + 6, round(Int, yy/0.1) + 6] = _circmat
    end
end
println("Done.\n")

"""
    gen_pic_circ!(out::AbstractMatrix, x, y; circ_ims=circs)

Generate a 1-channel image of a circle (via pre-computed `circ_ims`) on a
32×32 canvas supplied as the matrix `out` at the location (`x`, `y`), where
both coordinates are presumed to be in ``0 \\le (x,y) \\le 32``. (If the coods
are outside this range, the circle will be largely clipped, or absent from the
result.) Note that by "circle", I really mean "ball", although this may in
principle be changed in the pre-calculation function.

This is a mutating function which writes directly onto the supplied `out`
matrix. For a non-mutating version, use `gen_pic_circ(x,y)`.
"""
function gen_pic_circ!(out::AbstractMatrix, x, y; circ_ims=circs)
    ox, oy = x, y
    x, y = round(Int, x), round(Int, y)
    rx, ry = ox-x, oy-y
    rx, ry = ceil(Int, rx/0.1 + 5 +eps()), ceil(Int, ry/0.1 + 5 +eps())
    im = circ_ims[rx, ry]
    for xx in -3:2, yy in -3:2
        _x, _y = x+xx, y+yy
        _x = 32 - _x
        !(1 <= _x <= 32) && continue
        !(1 <= _y <= 32) && continue
        out[_x, _y] = max(out[_x, _y], im[xx+4, yy+4])
    end
    return out
end

gen_pic_circ(x, y) = let out = zeros(Float32, 32,32); o=gen_pic_circ!(out, x, y); end

"""
    gen_pic_circ2(x, y, x2, y2)

Draws TWO circles onto the same 32×32 canvas via the function `gen_pic_circ!`
(see docstring of that function). Circle centers are (x,y) and (x2, y2),
presumed to be in range ``0 \\le (x,y) \\le 32``.

Output: 32×32 matrix.
"""
function gen_pic_circ2(x,y,x2,y2)
    out = zeros(Float32, 32,32)
    gen_pic_circ!(out, x, y)
    gen_pic_circ!(out, x2, y2)
end

"""
    gen_pic_circ2(x, y, x2, y2)

Draws one circle onto each of a 2-channel 32×32 canvas via the function
`gen_pic_circ!` (see docstring of that function). Circle centers are (x,y) and
(x2, y2), presumed to be in range ``0 \\le (x,y) \\le 32``.

Output: 32×32×2 tensor.
"""
function gen_pic_circ2_2chan(x,y,x2,y2)
    out = zeros(Float32, 32,32,2)
    @views gen_pic_circ!(out[:,:,1], x, y)
    @views gen_pic_circ!(out[:,:,2], x2, y2)
    return out
end

imcood_tf(x,y) = let xy=7*[x+2, y+2].+2; [xy[1],xy[2]]; end

"""
    gen_pic_circ2_cood_tf(xy1, xy2)

A version of `gen_pic_circ2` (see docstring), but where centers xy1 and xy2 are
specified in the range [-2,2]×[-2,2] and transformed to the [0,32]×[0,32] range.
"""
function gen_pic_circ2_cood_tf(xy1, xy2)
    xy1, xy2 = imcood_tf(xy1...), imcood_tf(xy2...)
    gen_pic_circ2(xy1[2], xy1[1], xy2[2], xy2[1])
end

"""
    gen_pic_circ2_2chan_cood_tf(xy1, xy2)

A version of `gen_pic_circ2_2chan` (see docstring), but where centers xy1 and
xy2 are specified in the range [-2,2]×[-2,2] and transformed to the
[0,32]×[0,32] range.
"""
function gen_pic_circ2_2chan_cood_tf(xy1, xy2)
    xy1, xy2 = imcood_tf(xy1...), imcood_tf(xy2...)
    gen_pic_circ2_2chan(xy1[2], xy1[1], xy2[2], xy2[1])
end

################################################################################
##                                                                            ##
##            Solving and manipulating the Double Pendulum ODE                ##
##                                                                            ##
################################################################################

"""
    double_pendulum(du,u,p,t)

Definition of the double pendulum ODE. Started life from github user `zaman13`.
https://github.com/zaman13/Double-Pendulum-Motion-Animation/blob/master/Code/Double_pendulum-uniformly_spaced_time.ipynb.
See wikipedia article for equations of motion via a Lagrangian approach.
"""
function double_pendulum(du,u,p,t)
    # du = derivatives
    # u = variables
    # p = parameters
    # t = time variable

    m1, m2, L1, L2, g, α₁, α₂ = p;    # α = damping

    c = cos(u[1]-u[3]);  # intermediate variables
    s = sin(u[1]-u[3]);  # intermediate variables

    du[1] = u[2];   # d(theta 1)
    du[2] = ( m2*g*sin(u[3])*c - m2*s*(L1*c*u[2]^2 + L2*u[4]^2) - (m1+m2)*g*sin(u[1]) -α₁*u[2]) /(L1 *(m1+m2*s^2));
    du[3] = u[4];   # d(theta 2)
    du[4] = ((m1+m2)*(L1*u[2]^2*s - g*sin(u[3]) + g*sin(u[1])*c) + m2*L2*u[4]^2*s*c -α₂*u[4]) /(L2 * (m1+m2*s^2));
end

"""
    solve_dblpend(T=20; masses=[1,1], lengths=[1,1], g=9.8, dt=0.05,
    u0 = [2π/3; 0; π/3; 0.0], α₁=0, α₂=0)

Solving the double pendulum ODE defined in `double_pendulum`. The only
positional argument `T` defines the total timespan in terms of seconds, but *not
the length of the output*. The motion is sampled at the interval `dt`, and hence
the total output length is `T/dt`. The ODE parameters are given by the `masses`
and `lengths`, with elements corresponding to the first (inner) and second
(outer) pendulum respectively. The damping factor (which is applied to the
speed) is specified via `α₁` and `α₂` for each pendulum respecitvely. Finally,
`u0` is the initial condition (θ₁, θ̇₁, θ₂, θ̇₂), and `g` is the gravitational
constant.

Output: returns a tuple, the first of which is a matrix of size (T/dt)×4
where each row gives the state (θ₁, θ̇₁, θ₂, θ̇₂) for each point in time; the
second gives the solution in Euclidean coordinates (y₁, x₁, y₂, x₂) for
modelling purposes (and without velocity information).
"""
function solve_dblpend(T=20; masses=[1,1], lengths=[1,1], g=9.8, dt=0.05, u0 = [2π/3; 0; π/3; 0.0],
        α₁=0, α₂=0)
    m1, m2 = masses   # mass of pendulum 1,2 (in kg)
    L1, L2 = lengths  # length of pendulum 1,2 (in meters)
    # u0 == (pend1 θ, θ̇, pend2 θ, θ̇)

    # Solving the system
    pars = [m1, m2, L1, L2, g, α₁, α₂];
    tlims = (0.0, T);
    prob = ODEProblem(double_pendulum, u0, tlims, pars);
    ode_solver = Vern7() # Verner's "Most Efficient" 7/6 Runge-Kutta method
    sol = solve(prob, ode_solver, reltol=1e-6);

    z = reduce(hcat, sol(0:dt:(T-dt)).u)

    xy1 = vcat(L1*sin.(z[1,:]'), -L1*cos.(z[1,:]'))                         # First Pendulum
    xy2 = vcat(xy1[1,:]' + L2*sin.(z[3,:]'), xy1[2,:]' - L2*cos.(z[3,:]'))         # Second Pendulum

    return z, vcat(xy1, xy2)
end

# Various kinetic and potential energy calculations for a pendulum. The
# potential energy is a simple sum of the PE of the bobs, but the KE of the
# combined system requires the extra term, which is reflected in the total
# energy `energy_dblpend`.
ke(θ, θ̇, r=1) = ke(θ̇, r)
ke(θ̇, r=1) = 0.5*(r*θ̇)^2
extra_dbl_ke(θ₁, θ̇₁, θ₂, θ̇₂) = θ̇₁*θ̇₂*cos(θ₁-θ₂)  # extra term from the coupling
pe(θ, θ̇, r=1) = pe(θ)
pe(θ) = 9.8*(1-cos(θ))
energy_dblpend(θ₁, θ̇₁, θ₂, θ̇₂, r₁=1,r₂=1) = sum(2*f(θ₁,θ̇₁,r₁)+f(θ₂,θ̇₂,r₂) for f in [pe,ke]) +θ̇₁*θ̇₂*cos(θ₁-θ₂)

################################################################################
##                                                                            ##
##                            Generating the dataset                          ##
##                                                                            ##
################################################################################

rsgn(; rng=Random.GLOBAL_RNG) = rand(rng, [-1,+1])   # random sign

"""
    sample_init(target)

Sample the initial conditions (θ₁, θ̇₁, θ₂, θ̇₂) of the double pendulum system
with an energy specified by `target`. My initial experiments used the value
`3.5*9.8` which is the default if called without arguments.

It is not entirely uniformly at random over (θ₁, θ̇₁, θ₂, θ̇₂), since some
combinations led to fairly boring behaviour (~simple harmonic motion). I didn't
have time to analyze exactly what initial conditions led to this, but I put a
couple of intuitive (to me) changes in: requiring that the initial height of the
first pendulum (θ₁) is bounded above 0, as is θ₂. Notice also that the ± in the
quadratic formula for θ̇₂ is fixed as the *opposite* sign of θ̇₁. This appears to
avoid some of the synchronization between the pendula that I noticed from
sampling uniformly at random.
"""
function sample_init(target=3.5*9.8; rng=Random.GLOBAL_RNG)
    H = target/9.8
    ul = acos(max(1 - H/2, -1))
    θ₁ = rsgn(rng=rng)*rand(rng, Uniform(ul/2, ul))
    θ₂ = let a=acos(clip(3-H-2*cos(θ₁), -1, 1)); rsgn(rng=rng)*rand(rng, Uniform(a/2, a)); end
    κ = 9.8*((H-3) + 2*cos(θ₁) + cos(θ₂))  # remaining energy
    θ̇₁ = rand(rng, Uniform(0, sqrt(κ)))
    b, c = 2*θ̇₁*cos(θ₁-θ₂), 2*(θ̇₁^2 - κ)
    θ̇₂ = 0.5*(-b -sign(θ̇₁)*sqrt(b^2-4*c))
    return [θ₁, θ̇₁, θ₂, θ̇₂]
end

"""
    generate_data(;N=10, Nvalid=1, Ntest=1, tT=140, _seed=1240,
        max_energy=40, damping_vals=nothing, dt=0.05)

Generate sequences of length `tT` from the cartesian product of the damping_vals
(defaults to
    [0.0,  0.01,  0.02,  0.05,  0.1,  0.20,  0.30,  0.50,  0.75,  1.00]
for each bob). Each element in this cart. prod generates `N` sequences (def 10)
resulting in (10×10)×10 sequences of length `tT`. Also a validation and test
set obtain `Nvalid` and `Ntest` seqeunces from this cartesian product too. The
initial conditions are sampled from `sample_init` with a random energy with
a maximum specified as `max_energy`, with a bias towards higher energy. The
default behaviour is to sample the sequence every `dt=0.05` seconds, and hence
a `tT=140` sequence covers only 7 seconds of ODE rollout.

Output:
(1) `data_xy`: a Dict with keys `[:train]`, `[:valid]`, `[:test]`, each of which
contains a list of all (10×10)×N (default) sequences of length `tT`.

(2) `data_θ`: as above, but now with each row of each sequence containing the
values (θ₁, θ̇₁, θ₂, θ̇₂).

(3) `data_meta`: A Dict of the metadata of each sequence. E.g. the `[:train]`
key holds a list of the [length, %total length, initial energy, damping values,
and initial conditions] of each sequence.
"""
function generate_data(;N=10, Nvalid=1, Ntest=1, tT=140, _seed=1240,
        max_energy=40, damping_vals=nothing, dt=0.05)

    rng = MersenneTwister()
    Random.seed!(rng, _seed)

    damping_vals = something(damping_vals,
        [0.0,  0.01,  0.02,  0.05,  0.1,  0.20,  0.30,  0.50,  0.75,  1.00])

    data_xy = Dict(:train=>[], :valid=>[], :test=>[])
    data_θ = Dict(:train=>[], :valid=>[], :test=>[])
    data_meta = Dict(:train=>[], :valid=>[], :test=>[])

    for α₁ in damping_vals, α₂ in damping_vals
        lens = [repeat([tT], N), repeat([tT], Nvalid), repeat([tT], Ntest)]
        for (ls, tvt) in zip(lens, [:train, :valid, :test])
            for l in ls
                energy = max_energy*rand(rng, Beta(3,1))  # heuristically chosen: bias towards larger energy
                u0 = sample_init(energy; rng=rng);   # original: [2π/3 + 0.00001; 0; π/3 + 0.00001; 0.0]
                lc = ceil(Int, l*dt)
                states, traj = solve_dblpend(lc; u0=u0, α₁=α₁, α₂=α₂, dt=dt)
                data_pend = convert(Array{Float32}, Matrix(traj'))
                thetas = convert(Array{Float32}, Matrix(states'))
                push!(data_xy[tvt], data_pend)
                push!(data_θ[tvt], thetas)
                push!(data_meta[tvt], [lc÷dt, energy, α₁, α₂, u0...])
            end
        end
    end

    data_meta = Dict(k=>Matrix(reduce(hcat, v)) for (k,v) in data_meta)
    data_meta = Dict(k=>(convert(Vector{Int}, v[1,:]), v[1,:]/sum(v[1,:]), v[2,:], v[3:end,:]) for (k,v) in data_meta)

    return data_xy, data_θ, data_meta
end

rand_timeseries(data_xy, data_meta, tT, nbatch; variable_len=false,
                time_slices=true, entire_seq=false, entire_step=1) =
    rand_timeseries(Random.GLOBAL_RNG, data_xy, data_meta, tT, nbatch;
                    variable_len=variable_len, time_slices=time_slices,
                    entire_seq=entire_seq, entire_step=entire_step)

"""
    rand_timeseries([rng::Random.AbstractRNG,] data_xy, data_meta, tT, nbatch;
                         variable_len=false, time_slices=true, entire_seq=false, entire_step=1)

Generate `nbatch` time series from the double pendulum dataset represented by
the (train/test/valid) split given by `data_xy` and `data_meta`. These are
generated as in the `generate_data` function (see docstring). Each time series
will be of length `tT`.

Options: if dataset has variable length time series (I think it no longer can!)
then use `variable_len=true` to weight each time series by its length during
random sampling. For inputs to recurrent models it is often helpful to have the
data in timeslice form, i.e. a vector where each element corresponds to the
input at a time `t` for a given batch. See `time_slices=true/false`. Some models
encode the entire sequence length to infer the latent variable for which the
`entire_seq=true` option may be helpful. This may also be returned for every kth
value, for which specify `entire_step=k`.
"""
function rand_timeseries(rng::Random.AbstractRNG, data_xy, data_meta, tT, nbatch;
                         variable_len=false, time_slices=true, entire_seq=false, entire_step=1)
    if variable_len
        lens = data_meta[1] .- (tT - 1)
        seqixs = rand(rng, Categorical(lens/sum(lens)), nbatch)
        lens = lens[seqixs]
    else
        N = length(data_meta[1])
        seqixs = rand(rng, 1:N, nbatch)
        lens = fill(data_meta[1][1], nbatch)
    end
    ixs = [rand(1:l-tT+1) for l in lens]
    y = Flux.batch([data_xy[s][i .+ (0:tT-1), :] for (s,i) in zip(seqixs,ixs)])
    (time_slices) && (y = [y[tt,:,:] for tt in 1:tT])

    if entire_seq
        variable_len && @warn "Cannot batch variable_len sequences in general. This may error."
        fullseq = Flux.batch([data_xy[s] for s in seqixs])
        (time_slices) && (fullseq = [fullseq[tt,:,:] for tt in 1:entire_step:tT])
        return y, fullseq, seqixs, ixs
    end
    return y, seqixs, ixs
end

end

max_stars_repo_name: ornithos/mtds-dblpend

max_issues_repo_name: ornithos/mtds-dblpend

max_forks_repo_name: ornithos/mtds-dblpend


text:
"""
    File containing all the conventions used to name objects.
"""
bus_name(id_bus::Int64) = "BUS_$id_bus"
get_busid(busname::String) = busname[5:end]
basecase_scenario_name() = "BaseCase"

function variable_name(varname::String, bus::String, elemid::String, scenario::String)
    if elemid == ""
        return "$(scenario)_$(get_busid(bus))_$(varname)"
    else
        return "$(scenario)_$(get_busid(bus))_$(elemid)_$(varname)"
    end
end

function variable_name(varname::String, link::Link, elemid::String, scenario::String)
    if elemid == ""
        return "$(scenario)_$(get_busid(link.orig))-$(get_busid(link.dest))_$(varname)"
    else
        return "$(scenario)_$(get_busid(link.orig))-$(get_busid(link.dest))_$(elemid)_$(varname)"
    end
end

function cstrname_nodal_balance(cstrnames::Set{String}, scenario::String, bus::String)
    cstr = ""
    for cstrname in sort(collect(cstrnames))
        cstr = "$(cstr)_$cstrname"
    end

    if cstr == ""
      cstr = "_LOAD"
    end

    return "$(scenario)_$(get_busid(bus))_BALANCE$(cstr)"
end

function get_cstrname(scenario::String, bus::String, elemid::String, cstrname::String)
    if elemid == ""
        return "$(scenario)_$(get_busid(bus))_$(cstrname)"
    else
        return "$(scenario)_$(get_busid(bus))_$(elemid)_$(cstrname)"
    end
end

function get_cstrname(scenario::String, link::Link, elemid::String, cstrname::String)
    if elemid == ""
        return "$(scenario)_$(get_busid(link.orig))-$(get_busid(link.dest))_$(cstrname)"
    else
        return "$(scenario)_$(get_busid(link.orig))-$(get_busid(link.dest))_$(elemid)_$(cstrname)"
    end
end

### GOC conventions
scenarioname(id_contingency::Int) = "Scen$id_contingency"

# nodal element names
volt_name() = "Volt"
generator_name(gen_id) = "Gen$gen_id"
load_name(id_load) = "Load$id_load"
shunt_name(id_shunt) = "Shunt$id_shunt"

get_delta_varname(scenario::String) = "$(scenario)_Delta"
get_binInf_varname(scenario1::String, scenario2::String, bus::String) = "BinVolt_$(get_busid(bus))_$(scenario1)_inf_$(scenario2)"

# constraint names
get_VoltM_cstrname() = "VOLTM"
get_ShuntDef_cstrname() = "ShuntDef"
get_LoadDef_cstrname() = "LoadDef"
get_GenBounds_cstrname() = "GenBounds"

get_VoltBinDef_upper() = "BinDef_upper"
get_VoltBinDef_lower() = "BinDef_lower"
get_VoltBinDef_complement() = "BinDef_compl"

get_CC_active_cstrname() = "CC_Pgen"
get_CC_reactiveupper_cstrname() = "CC_Qgen_upper"
get_CC_reactivelower_cstrname() = "CC_Qgen_lower"

get_Smax_orig_cstrname() = "Smax_orig"
get_Smax_dest_cstrname() = "Smax_dest"

get_GOC_Volt_ϵ() = 1e-3
get_GOC_BigM() = 1e1

##node elems

##link elems
nullimpedance_notransformer_name(branch_id) = "NullImp_notransfo_$branch_id"
linepi_withtransformer_name(branch_id) = "Lineπ_transfo_$branch_id"
nullimpedance_withtransformer_name(branch_id) = "NullImp_transfo_$branch_id"
linepi_notransformer_name(branch_id) = "Lineπ_notransfo_$branch_id"

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
module EMLsplitter
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("splitter/splitter_n.jl")
	include("splitter/splitter2.jl")
	include("splitter/splitter3.jl")
	include("splitter/splitter_column.jl")
end


text:
# Autogenerated wrapper script for TinyXML_jll for i686-w64-mingw32-cxx03
export libtinyxml

JLLWrappers.@generate_wrapper_header("TinyXML")
JLLWrappers.@declare_library_product(libtinyxml, "libtinyxml.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtinyxml,
        "bin\\libtinyxml.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
Max-affine neural network [1].

# Note
If you specify `n` or `m`, it can also be regarded as bivariate function.

# Variables
x ∈ ℝ^n
u ∈ ℝ^m
z = [xᵀ, uᵀ]ᵀ ∈ ℝ^(n+m)  # l = n+m
α_is: a vector of subgradients, i.e., α_is[i] ∈ ℝ^(n+m).
β_is: a vector of bias terms, i.e., β_is[i] ∈ ℝ.

# References
[1] G. C. Calafiore, S. Gaubert, and C. Possieri, “Log-Sum-Exp Neural Networks and Posynomial Models for Convex and Log-Log-Convex Data,” IEEE Transactions on Neural Networks and Learning Systems, vol. 31, no. 3, pp. 827–838, Mar. 2020, doi: 10.1109/TNNLS.2019.2910417.
"""
struct MA <: ConvexApproximator
    l::Int  # n+m
    n::Int  # the first variable for bivariate function
    m::Int  # the second variable for bivariate function
    i_max::Int
    _α_is::Matrix
    _β_is::Matrix
    function MA(α_is::Vector, β_is::Vector; n=nothing, m=nothing)
        l, i_max, _α_is, _β_is = _construct_convex_approximator(α_is, β_is)
        if n != nothing && m != nothing
            @assert l == n + m
        elseif n != nothing
            m = l - n
            @assert m > 0
        elseif m != nothing
            n = l - m
            @assert n > 0
        end
        new(l, n, m, i_max, _α_is, _β_is)
    end
end
Flux.@functor MA (_α_is, _β_is,)

function (nn::MA)(z::Array)
    is_vector = length(size(z)) == 1
    z_affine = affine_map(nn, z)
    _res = maximum(z_affine; dims=1)
    res = is_vector ? reshape(_res, 1) : _res
end

function (nn::MA)(z::Convex.AbstractExpr)
    z_affine = affine_map(nn, z)
    _res = [maximum(z_affine)]
end

"""
Considering bivariate function approximator
"""
function (nn::MA)(x::AbstractArray, u::AbstractArray)
    nn(vcat(x, u))
end

function (nn::MA)(x::AbstractArray, u::Convex.AbstractExpr)
    nn(vcat(Convex.Constant(x), u))
end


text:
using Documenter
push!(LOAD_PATH,"../src/")
using Check3DPackings, Chen1995, Kurpel2018

makedocs(
  modules = [Check3DPackings, Chen1995, Kurpel2018],
  format = :html,
  sitename = "HBDModules.jl",
  doctest = true,
  html_prettyurls = false
)