text:
using NEOS, Base.Test

using MathProgBase

@testset "MPSWriter" begin
	include(joinpath(dirname(@__FILE__), "MPSWriter.jl"))
end

# Null out this method for testing
NEOS.getobjbound(m::NEOS.MPSModel) = NaN

#
# As part of the NEOS terms of use, some solvers
# 	require a valid email to be supplied.
#
# TESTING_EMAIL = readchomp(`git config --get user.email`)
TESTING_EMAIL = "odow@users.noreply.github.com"
@testset "Old solvers" begin
	s = NEOSServer()
	@test_throws Exception NEOSCPLEXSolver(s)
	@test_throws Exception NEOSMOSEKSolver(s)
	@test_throws Exception NEOSXpressSolver(s)
	@test_throws Exception NEOSSYMPHONYSolver(s)
end
@testset "Test NEOS Server" begin
	s = NEOSServer()
	addemail!(s, TESTING_EMAIL)
	@test s.email == TESTING_EMAIL
	addemail!(s, "")
	s = NEOSServer(email=TESTING_EMAIL)
	@test s.email == TESTING_EMAIL
	v = version(s)
	@test match(r"neos version ([0-9]+) \(.+\)", v) != nothing

	@test contains(neosHelp(s), "class NeosServer")

	_solvers = listAllSolvers(s)
	@test "milp:CPLEX:MPS" in _solvers
	@test "milp:MOSEK:MPS" in _solvers
	@test "lp:MOSEK:MPS" in _solvers
	@test "milp:SYMPHONY:MPS" in _solvers
	@test "milp:FICO-Xpress:MPS" in _solvers

	@test "Mixed Integer Linear Programming" in listCategories(s)

	_solvers = listSolversInCategory(s, "milp")
	@test "CPLEX:MPS" in _solvers
	@test "MOSEK:MPS" in _solvers
	@test "SYMPHONY:MPS" in _solvers
	@test "FICO-Xpress:MPS" in _solvers

	j = NEOS.NEOSJob(3804943, "OfRcoMbp")
	_info = getJobInfo(s, j)
	@test _info[1] == "milp"
	@test _info[2] == "FICO-Xpress"
	@test _info[3] == "MPS"
	@test _info[4] == "Done"

	@test getJobStatus(s, j) == "Done"
	@test killJob(s, j) == "Job #3804943 is finished"
	@test getIntermediateResults(s, j)[1] == getIntermediateResultsNonBlocking(s, j)[1] == "Results for Job #3804943 are no longer available"
end

@testset "Test MPSModel" begin
	m = MathProgBase.LinearQuadraticModel(NEOSSolver(solver=:SYMPHONY))
	@test isa(m.solver, NEOS.NEOSSolver{:SYMPHONY, :MPS})
	@test MathProgBase.getsolution(m) == []
	@test MathProgBase.getobjval(m) == 0.
	@test MathProgBase.getsense(m) == :Min
	MathProgBase.setsense!(m, :Max)
	@test MathProgBase.getsense(m) == :Max
	@test MathProgBase.status(m) == NEOS.NOTSOLVED
	MathProgBase.loadproblem!(m, [1. 2. 3.; 1. 1. 1.], [-1., 0., 0.], [1., 1., Inf], [0., 0., 1.], [1.25, 1.], [1.25, 1.], :Max)
	MathProgBase.setvartype!(m, [:SemiCont, :Cont, :Bin])
	@test_throws Exception NEOS.addCOLS(m, "")
	MathProgBase.setvartype!(m, [:SemiInt, :Cont, :Bin])
	@test_throws Exception NEOS.addCOLS(m, "")
end

function basictest(solver, timelimit)
	@testset "Test basic solver stuff for $(typeof(solver))" begin
		@test isa(solver, NEOS.NEOSSolver)
		fields = fieldnames(solver)
		for sym in [:server, :requires_email, :solves_sos, :provides_duals,
			:template, :params, :gzipmodel, :print_results, :result_file]
			@test sym in fields
		end

		m = MathProgBase.LinearQuadraticModel(solver)

		if solver.provides_duals
			m.nrow, m.ncol = 1, 2
			@test MathProgBase.getreducedcosts(m) == []
			@test MathProgBase.getconstrduals(m) == []
		end

		if !solver.solves_sos
			@test_throws Exception MathProgBase.addsos1!(m, [], [])
			@test_throws Exception MathProgBase.addsos2!(m, [], [])
		end

		addemail!(m, "Test")
		@test m.solver.server.email == "Test"
		addemail!(solver, TESTING_EMAIL)
		@test solver.server.email == TESTING_EMAIL

		addparameter!(solver, "key", 0)
		@test solver.params["key"] == 0
		solver.params = Dict{String, Any}()

		addparameter!(solver, timelimit, 60)
		# _solver = @eval $(s)($(timelimit)=60, email=TESTING_EMAIL)
		@test solver.params[string(timelimit)] == 60
		@test solver.server.email == TESTING_EMAIL

		if solver.requires_email
			addemail!(solver, "")
			@test_throws Exception linprog([-1.0,0.0;],sparse([2.0 -1.0;]),'<',1.5, [-1.0, -Inf], [1.0, 0.0], solver)
			addemail!(solver, TESTING_EMAIL)
		end
	end
end
function feasibleproblem(solver, timelimit)
    @testset "Testing feasible problem $(typeof(solver))" begin
	    sol = linprog([-1.0,0.0;],sparse([2.0 -1.0;]),'<',1.5, [-1.0, -Inf], [1.0, 0.0], solver)
	    @test sol.status == :Optimal
	    @test isapprox(sol.objval, -0.75, atol=1e-5)
	    @test  isapprox(sol.sol, [0.75, 0.;], atol=1e-5)
		if solver.provides_duals
	    	@test isapprox(sol.attrs[:lambda], [-0.5;], atol=1e-5)
	    	@test isapprox(sol.attrs[:redcost], [0., -0.5;], atol=1e-5)
		end

		addparameter!(solver, string(timelimit), 60)
        sol = mixintprog(-[5.0,3.0,2.0,7.0,4.0;],Float64[2.0 8.0 4.0 2.0 5.0;],'<',10.0,:Int,-0.5,1.0,solver)
		@test sol.status == :Optimal
 		@test isapprox(sol.objval, -16.0, atol=1e-6)
		@test isapprox(sol.sol, [1.0, 0.0, 0.0, 1.0, 1.0;], atol=1e-4)
	end
end
function infeasibleproblem(solver, timelimit)
	@testset "Testing infeasible problem $(typeof(solver))" begin
		solver.gzipmodel=false
	    sol = linprog([1.0,0.0;],[-2.0 -1.0;],'>',1.0, solver)
		@test (sol.status == NEOS.INFEASIBLE || sol.status == NEOS.UNBNDORINF)
	end
end
function unboundedproblem(solver, timelimit)
	@testset "Testing unbounded problem $(typeof(solver))" begin
	    solver.result_file = randstring(5)
		sol = linprog([-1.0,-1.0;],[-1.0 2.0;],'<',[0.0;], solver)
	    @test (sol.status == NEOS.UNBOUNDED || sol.status == NEOS.UNBNDORINF)
		@test length(readstring(solver.result_file)) > 0
		rm(solver.result_file)
		solver.result_file = ""
    end
end
function nullproblem(solver, timelimit)
	@testset "Testing null problem $(typeof(solver))" begin
		solver.result_file = randstring(5)
		m = MathProgBase.LinearQuadraticModel(solver)
		MathProgBase.loadproblem!(m, Array{Float64}(0,1), [0.0], [Inf], [1.0], [], [], :Min)
		MathProgBase.optimize!(m)
		rm(solver.result_file)
		solver.result_file = ""
	end
end
function sosproblem(solver, timelimit)
	@testset "Testing SOS problem $(typeof(solver))" begin
		m = MathProgBase.LinearQuadraticModel(solver)
		MathProgBase.loadproblem!(m, [1.0 2.0 3.0; 1.0 1.0 1.0], [-1.0, 0.0, 0.0], [1.0, 1.0, Inf], [0.0, 0.0, 1.0], [1.25, 1.0], [1.25, 1.0], :Max)
		MathProgBase.addsos2!(m, [1, 2, 3], [1.0, 2.0, 3.0])
		@test MathProgBase.optimize!(m) == NEOS.OPTIMAL
		@test isapprox(MathProgBase.getobjval(m), 0., atol=1e-6)
		@test isapprox(MathProgBase.getsolution(m), [0.75, 0.25, 0.0], atol=1e-6)

		m = MathProgBase.LinearQuadraticModel(solver)
		MathProgBase.loadproblem!(m, [1.0 1.0 1.0], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0], [1.0, 3.0, 2.0], [0.0], [1.5], :Max)
		MathProgBase.addsos1!(m, [1, 2, 3], [1.0, 2.0, 3.0])
		@test MathProgBase.optimize!(m) == NEOS.OPTIMAL
		@test isapprox(MathProgBase.getobjval(m), 3.0, atol=1e-6)
		@test isapprox(MathProgBase.getsolution(m), [0.0, 1.0, 0.0], atol=1e-6)
	end
end

SOLVERS = [
	(NEOSSolver(solver=:CPLEX, format=:NL), :timelimit),
	(NEOSSolver(solver=:CPLEX, format=:MPS), :timelimit),
	(NEOSSolver(solver=:MOSEK, format=:MPS), :MSK_DPAR_OPTIMIZER_MAX_TIME),
	(NEOSSolver(solver=:SYMPHONY, format=:MPS), :time_limit),
	(NEOSSolver(solver=:Xpress, format=:MPS), :MAXTIME)
]

nullproblem(NEOSSolver(solver=:CPLEX, format=:NL, print_level=2, email=TESTING_EMAIL), :timelimit)

for (solver, timelimit) in SOLVERS
	basictest(solver, timelimit)
	feasibleproblem(solver, timelimit)
	infeasibleproblem(solver, timelimit)
	unboundedproblem(solver, timelimit)
	nullproblem(solver, timelimit)
	!solver.solves_sos && continue
	sosproblem(solver, timelimit)
end


text:
using Materials, Test
using Materials: Simulator
using PyPlot

num = 100 # Number of steps
E = 200000.0 # Young's modulus [GPa]
poissons_ratio = 0.3 # Poisson's ratio
stress_y = 200.0 # Yield stress
dt = 0.1
n = 0.92
K = 180.0e3

# Simulation 1
mat = Material(ViscoPlastic, (:norton, [K, n]))
mat.properties.youngs_modulus = E
mat.properties.poissons_ratio = poissons_ratio
mat.properties.yield_stress = stress_y
times = [0.0]
strains = [zeros(6)]

L = 1000.0               # Initial length [mm]
dL = range(0, stop=3.1, length=num) # Change of length [mm]

strain = collect(dL / L)
# strain_total = vec([strain; reverse(strain[2:end-1]); -strain; reverse(-strain[2:end-1])])
strain_total = strain
nsteps = length(strain_total)
dt = 0.03

for i=1:nsteps
    str = zeros(6)
    str[1] = strain_total[i]
    str[2] = -strain_total[i] * poissons_ratio
    str[3] = -strain_total[i] * poissons_ratio

    push!(times, times[end] + dt)
    push!(strains, str)
end

println(strains[2][1] / dt)
sim = Simulator(mat)
Materials.initialize!(sim, strains, times)
Materials.run!(sim)

s11s = [s[1] for s in sim.stresses]
e11s = [e[1] for e in strains]

# Simulation 2
mat = Material(ViscoPlastic, (:norton, [K, n]))
mat.properties.youngs_modulus = E
mat.properties.poissons_ratio = poissons_ratio
mat.properties.yield_stress = stress_y
times = [0.0]
strains = [zeros(6)]

dt = 0.08

for i=1:nsteps
    str = zeros(6)
    str[1] = strain_total[i]
    str[2] = -strain_total[i] * poissons_ratio
    str[3] = -strain_total[i] * poissons_ratio

    push!(times, times[end] + dt)
    push!(strains, str)
end
println(strains[3][1] / dt)
sim = Simulator(mat)
Materials.initialize!(sim, strains, times)
Materials.run!(sim)

s11s2 = [s[1] for s in sim.stresses]
e11s2 = [e[1] for e in strains]

# Simulation 3
mat = Material(ViscoPlastic, (:norton, [K, n]))
mat.properties.youngs_modulus = E
mat.properties.poissons_ratio = poissons_ratio
mat.properties.yield_stress = stress_y
times = [0.0]
strains = [zeros(6)]

dt = 0.5

for i=1:nsteps
    str = zeros(6)
    str[1] = strain_total[i]
    str[2] = -strain_total[i] * poissons_ratio
    str[3] = -strain_total[i] * poissons_ratio

    push!(times, times[end] + dt)
    push!(strains, str)
end
println(strains[3][1] / dt)
sim = Simulator(mat)
Materials.initialize!(sim, strains, times)
Materials.run!(sim)

s11s3 = [s[1] for s in sim.stresses]
e11s3 = [e[1] for e in strains]
legend(bbox_to_anchor=[1.05,1],loc=2,borderaxespad=0)
plot(e11s, s11s, label="de/dt=0.001346")
plot(e11s2, s11s2, label="de/dt=0.00081")
plot(e11s3, s11s3, label="de/dt=0.000008")
grid()
show()

ext: jl


text:
@testset "Globals" begin

@test isempty(lintstr("""
y=1
function f(x)
    x + y
end
"""))

@test isempty(lintstr("""
y = 1
function f(x)
    global y = 2
    x + y
end
"""))

msgs = lintstr("""
yyyyyy = 1
function f(x)
    yyyyyy = 2
    x + yyyyyy
end
""")
@test msgs[1].code == :I341
@test msgs[1].variable == "yyyyyy"
@test occursin("local variable", msgs[1].message)
@test occursin("shadows global variable", msgs[1].message)

s = """
y = 1
function f(x)
    y = 2
    x + y
end
"""
msgs = lintstr(s)
@test_broken isempty(msgs) # short names are grandfathered to be ok

s = """
const y = 1
function f(x)
    y = 2
    x + y
end
"""
msgs = lintstr(s)
@test_broken isempty(msgs) # short names are grandfathered to be ok

s = """
const y
function f(x)
    y = 2
    x + y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E111
@test occursin("expected assignment after \\\"const\\\"", msgs[1].message)

s = """
global 5
"""
msgs = lintstr(s)
@test msgs[1].code == :E134
@test occursin("unknown global pattern", msgs[1].message)

s = """
f() = x
x = 5
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
x
x = 5
"""
msgs = lintstr(s)
@test msgs[1].code == :E321
@test occursin("use of undeclared symbol", msgs[1].message)

# Test gloabls defined in other files
# File in package src
msgs = lintfile("FakePackage/src/subfolder2/file2.jl")
@test_broken isempty(msgs)
# File in package test
msgs = lintfile("FakePackage/test/file2.jl")
@test_broken isempty(msgs)
# File in base julia
msgs = lintfile("FakeJulia/base/file2.jl")
@test_broken isempty(msgs)

msgs = lintfile("filename","something")
@test msgs[1].code == :E321
end


text:
##############################
## by Qin Yu, Nov 2018
## using Julia 1.0.1
##############################

############################## Starting here:
using LinearAlgebra
using Statistics
using Random
using Plots
using MAT  # Windows only
using JLD  # macOS Friendly
push!(LOAD_PATH, ".")
using SuperLearn
gr()

############################## The Boston Housing Dataset:
# Detailed Description from Toronto University:
### http://www.cs.toronto.edu/~delve/data/boston/bostonDetail.html

# MAT File from Prof. Mark Herbster, University College London:
### http://www0.cs.ucl.ac.uk/staff/M.Herbster/boston

# This dataset contains information collected by the U.S Census Service
#   concerning housing in the area of Boston Mass. It was obtained
#   from the StatLib archive (http://lib.stat.cmu.edu/datasets/boston), and
#   has been used extensively throughout the literature to benchmark algorithms.
# However, these comparisons were primarily done outside of Delve
#   and are thus somewhat suspect.
# The dataset is small in size with only 506 cases.
# The data was originally published
#   by Harrison, D. and Rubinfeld, D.L. `Hedonic prices and the demand for clean air',
#   J. Environ. Economics & Management, vol.5, 81-102, 1978.

#------------------------------
# The Boston Housing Dataset:
# There are 14 attributes in each case of the dataset. They are:
# CRIM - per capita crime rate by town
# ZN - proportion of residential land zoned for lots over 25,000 sq.ft.
# INDUS - proportion of non-retail business acres per town.
# CHAS - Charles River dummy variable (1 if tract bounds river; 0 otherwise)
# NOX - nitric oxides concentration (parts per 10 million)
# RM - average number of rooms per dwelling
# AGE - proportion of owner-occupied units built prior to 1940
# DIS - weighted distances to five Boston employment centres
# RAD - index of accessibility to radial highways
# TAX - full-value property-tax rate per $10,000
# PTRATIO - pupil-teacher ratio by town
# B - 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
# LSTAT - % lower status of the population
# MEDV - Median value of owner-occupied homes in $1000's
#------------------------------

############################## Processing From :
#----- (for my Windows 10)
# Import Boston:
matdata = matread("boston.mat")["boston"]
# If MAT package is not available for your version of Julia,
# use JLD to prepare data on another machine:
save("./boston.jld", "boston", matdata)

#----- (for my macOS)
# use JLD to load the transfored data:
matdata = load("boston.jld")["boston"]
# and have a look at all data of y:
plot(matdata[:, 14])

############################## Naive Regression (i.e. Linear Regression with 0 Attribute):
# Generating test/training sets:
nrow, ncol = size(matdata)
nrow_test  = div(nrow, 3)
nrow_train = nrow - nrow_test
ones_test  = ones(nrow_test)
ones_train = ones(nrow_train)

# 20 runs of randomly splitted datasets:
plot(matdata[:, 14], lab="boston[:,14]")
sum_all_20_te = sum_all_20_tse = 0
for i = 1:20
    # Obtain randomly splitted 𝒚_test/training:
    permuted_matdata = matdata[randperm(nrow),:]
    𝒚_train = permuted_matdata[nrow_test+1:nrow, 14]
    𝒚_test  = permuted_matdata[1:nrow_test, 14]

    # Plot each run:
    trl(x_test) = trained_regression_line(x_test, ones_train, 𝒚_train, 1)
    plot!(trl, 1, nrow, lab="$i th run")

    # Obtain MSEs:
    global sum_all_20_te += training_error_k_dim_basis(ones_train, 𝒚_train, 1)
    global sum_all_20_tse += test_error_k_dim_basis(ones_test, 𝒚_test, ones_train, 𝒚_train, 1)
end
current()
savefig("./graph/3.1.pdf")
# Obtain average MSEs for testing/training over 20 runs:
avg_all_20_te  = sum_all_20_te / 20
avg_all_20_tse = sum_all_20_tse / 20

############################## Linear Regression with 1 Attribute:
# Using SuperLearn without basis function (but with integrated bias/intercept term, [ ,1]):
Core.eval(SuperLearn, :(TRANS_BASIS = false))

# the 𝑖th element is for the 𝑖th 𝑥, and the last entry saved for next section: all-attributes-regression
sum_all_20_te = sum_all_20_tse = avg_all_20_te = avg_all_20_tse = zeros(14)
for j = 1:20
    permuted_matdata = matdata[randperm(nrow),:]
    for i = 1:13
        𝒙_test  = permuted_matdata[1:nrow_test, i]
        𝒙_train = permuted_matdata[nrow_test+1:nrow, i]
        𝒚_test  = permuted_matdata[1:nrow_test, 14]
        𝒚_train = permuted_matdata[nrow_test+1:nrow, 14]
        scatter(𝒙_train, 𝒚_train)

        trl(x_test) = trained_regression_line(x_test, 𝒙_train, 𝒚_train, nothing)
        display(plot!(trl, minimum(𝒙_train), maximum(𝒙_train)))

        sum_all_20_te[i]  += training_error_k_dim_basis(𝒙_train, 𝒚_train, nothing)
        sum_all_20_tse[i] += test_error_k_dim_basis(𝒙_test, 𝒚_test, 𝒙_train, 𝒚_train, nothing)
        j == 1 && savefig("./graph/3.2-$i.pdf")
    end
end
avg_all_20_te  = sum_all_20_te / 20
avg_all_20_tse = sum_all_20_tse / 20

############################## Linear Regression with All Attributes
plot(matdata[:, 14])
X_test  = matdata[1:nrow_test, 1:13]
X_train = matdata[nrow_test+1:nrow, 1:13]
𝒚_test  = matdata[1:nrow_test, 14]
𝒚_train = matdata[nrow_test+1:nrow, 14]
trl(X_test) = trained_regression_line_M(X_test, X_train, 𝒚_train, nothing)
scatter!(nrow_test+1:nrow, trl(X_train))

# Have a look at our training result:
sorted_matdata = sort_matrix_by_jth_col(matdata, 14)
plot(sorted_matdata[:, 14])
plot!(1:nrow, trl(sorted_matdata[:, 1:13]))

sum_all_20_te[14] = sum_all_20_tse[14] = 0
sorted_matdata = sort_matrix_by_jth_col(matdata, 14)
plot(sorted_matdata[:, 14], lab="real y value")
for j = 1:20
    permuted_matdata = matdata[randperm(nrow),:]
    X_test  = permuted_matdata[1:nrow_test, 1:13]
    X_train = permuted_matdata[nrow_test+1:nrow, 1:13]
    𝒚_test  = permuted_matdata[1:nrow_test, 14]
    𝒚_train = permuted_matdata[nrow_test+1:nrow, 14]

    trl(X_test) = trained_regression_line_M(X_test, X_train, 𝒚_train, nothing)
    display(plot!(1:nrow, trl(sorted_matdata[:, 1:13]), lab="$j th run"))  # current() outside loop, display() inside!!

    sum_all_20_te[14]  += training_error_k_dim_basis(X_train, 𝒚_train, nothing)
    sum_all_20_tse[14] += test_error_k_dim_basis(X_test, 𝒚_test, X_train, 𝒚_train, nothing)
end
savefig("./graph/3.3.pdf")
avg_all_20_te[14]  = sum_all_20_te[14] / 20
avg_all_20_tse[14] = sum_all_20_tse[14] / 20

############################## Comparing the Result (More Detailed Comparisons in 4.jl)
plot(avg_all_20_te, xlabel="attributes, 14th is all-attribute", xticks=0:14, ylabel="MSE", lab="avg training MSE")
plot!(avg_all_20_tse, lab="avg testing MSE")
savefig("./graph/3.4.pdf")

max_stars_repo_name: qin-yu/julia-regression-boston-housing

max_issues_repo_name: qin-yu/ml-julia-boston-housing

max_forks_repo_name: qin-yu/ml-julia-boston-housing


text:
"""
`SQLite.Source(db, sql, values=[]; rows::Int=0, stricttypes::Bool=true)`

Independently constructs an `SQLite.Source` in `db` with the SQL statement `sql`.
Will bind `values` to any parameters in `sql`.
`rows` is used to indicate how many rows to return in the query result if known beforehand. `rows=0` (the default) will return all possible rows.
`stricttypes=false` will remove strict column typing in the result set, making each column effectively `Vector{Any}`. `nullable::Bool=true` indicates
whether to allow null values when fetching results; if set to `false` and a null value is encountered, a `NullException` will be thrown.

Note that no results are returned; `sql` is executed, and results are ready to be returned (i.e. streamed to an appropriate `Data.Sink` type)
"""
function Source(db::DB, sql::AbstractString, values=[]; rows::Int=-1, stricttypes::Bool=true, nullable::Bool=true)
    stmt = SQLite.Stmt(db, sql)
    bind!(stmt, values)
    status = SQLite.execute!(stmt)
    cols = SQLite.sqlite3_column_count(stmt.handle)
    header = Array(String, cols)
    types = Array(DataType, cols)
    for i = 1:cols
        header[i] = unsafe_string(SQLite.sqlite3_column_name(stmt.handle, i))
        # do better column type inference; query what the column was created for?
        if nullable
            types[i] = Nullable{stricttypes ? SQLite.juliatype(stmt.handle, i) : Any}
        else
            types[i] = stricttypes ? SQLite.juliatype(stmt.handle, i) : Any
        end
    end
    return SQLite.Source(Data.Schema(header, types, rows), stmt, status)
end

"""
`SQLite.Source(sink::SQLite.Sink, sql="select * from \$(sink.tablename)")`

constructs an SQLite.Source from an SQLite.Sink; selects all rows/columns from the underlying Sink table by default
"""
Source(sink::SQLite.Sink,sql::AbstractString="select * from $(sink.tablename)") = Source(sink.db, sql::AbstractString)

function juliatype(handle,col)
    x = SQLite.sqlite3_column_type(handle, col)
    if x == SQLITE_BLOB
        val = sqlitevalue(Any,handle, col)
        return typeof(val)
    else
        return juliatype(x)
    end
end
juliatype(x) = x == SQLITE_INTEGER ? Int : x == SQLITE_FLOAT ? Float64 : x == SQLITE_TEXT ? String : Any

sqlitevalue{T<:Union{Signed,Unsigned}}(::Type{T}, handle, col) = convert(T, sqlite3_column_int64(handle, col))
const FLOAT_TYPES = Union{Float16, Float32, Float64} # exclude BigFloat
sqlitevalue{T<:FLOAT_TYPES}(::Type{T}, handle, col) = convert(T, sqlite3_column_double(handle, col))
#TODO: test returning a WeakRefString instead of calling `bytestring`
sqlitevalue{T<:AbstractString}(::Type{T}, handle, col) = convert(T,unsafe_string(sqlite3_column_text(handle, col)))
function sqlitevalue{T}(::Type{T}, handle, col)
    blob = convert(Ptr{UInt8}, sqlite3_column_blob(handle, col))
    b = sqlite3_column_bytes(handle, col)
    buf = zeros(UInt8, b) # global const?
    unsafe_copy!(pointer(buf), blob, b)
    r = sqldeserialize(buf)::T
    return r
end

# DataStreams interface
function Data.isdone(s::Source, row, col)
    (s.status == SQLITE_DONE || s.status == SQLITE_ROW) || sqliteerror(s.stmt.db)
    return s.status == SQLITE_DONE
end
# resets an SQLite.Source, ready to read data from at the start of the resultset
Data.reset!(io::SQLite.Source) = (sqlite3_reset(io.stmt.handle); execute!(io.stmt))
Data.streamtype{T<:SQLite.Source}(::Type{T}, ::Type{Data.Field}) = true

# `T` might be Int, Float64, String, WeakRefString, any Julia type, Any, NullType
# `t` (the actual type of the value we're returning), might be SQLITE_INTEGER, SQLITE_FLOAT, SQLITE_TEXT, SQLITE_BLOB, SQLITE_NULL
# `SQLite.getfield` returns the next `Nullable{T}` value from the `SQLite.Source`
function Data.getfield{T}(source::SQLite.Source, ::Type{Nullable{T}}, row, col)
    handle = source.stmt.handle
    t = SQLite.sqlite3_column_type(handle, col)
    if t == SQLite.SQLITE_NULL
        val = Nullable{T}()
    else
        TT = SQLite.juliatype(t) # native SQLite Int, Float, and Text types
        val = Nullable{T}(sqlitevalue(ifelse(TT === Any && !isbits(T), T, TT), handle, col))
    end
    col == source.schema.cols && (source.status = sqlite3_step(handle))
    return val
end
function Data.getfield{T}(source::SQLite.Source, ::Type{T}, row, col)
    handle = source.stmt.handle
    t = SQLite.sqlite3_column_type(handle, col)
    if t == SQLite.SQLITE_NULL
        throw(NullException)
    else
        TT = SQLite.juliatype(t) # native SQLite Int, Float, and Text types
        val::T = sqlitevalue(ifelse(TT === Any && !isbits(T), T, TT), handle, col)
    end
    col == source.schema.cols && (source.status = sqlite3_step(handle))
    return val
end

"""
`SQLite.query(db, sql::String, sink=DataFrame, values=[]; rows::Int=0, stricttypes::Bool=true)`

convenience method for executing an SQL statement and streaming the results back in a `Data.Sink` (DataFrame by default)

Will bind `values` to any parameters in `sql`.
`rows` is used to indicate how many rows to return in the query result if known beforehand. `rows=0` (the default) will return all possible rows.
`stricttypes=false` will remove strict column typing in the result set, making each column effectively `Vector{Any}`
"""
function query(db::DB, sql::AbstractString, sink=DataFrame, args...; append::Bool=false, values=[], rows::Int=-1, stricttypes::Bool=true, nullable::Bool=true)
    source = Source(db, sql, values; rows=rows, stricttypes=stricttypes, nullable=nullable)
    sink = Data.stream!(source, sink, append, args...)
    Data.close!(sink)
    return sink
end

function query{T}(db::DB, sql::AbstractString, sink::T; append::Bool=false, values=[], rows::Int=-1, stricttypes::Bool=true, nullable::Bool=true)
    source = Source(db, sql, values; rows=rows, stricttypes=stricttypes, nullable=nullable)
    sink = Data.stream!(source, sink, append)
    Data.close!(sink)
    return sink
end
query(source::SQLite.Source, sink=DataFrame, args...; append::Bool=false) = (sink = Data.stream!(source, sink, append, args...); Data.close!(sink); return sink)
query{T}(source::SQLite.Source, sink::T; append::Bool=false) = (sink = Data.stream!(source, sink, append); Data.close!(sink); return sink)

"""
`SQLite.tables(db, sink=DataFrame)`

returns a list of tables in `db`
"""
tables(db::DB, sink=DataFrame) = query(db, "SELECT name FROM sqlite_master WHERE type='table';", sink)

"""
`SQLite.indices(db, sink=DataFrame)`

returns a list of indices in `db`
"""
indices(db::DB, sink=DataFrame) = query(db, "SELECT name FROM sqlite_master WHERE type='index';", sink)

"""
`SQLite.columns(db, table, sink=DataFrame)`

returns a list of columns in `table`
"""
columns(db::DB,table::AbstractString, sink=DataFrame) = query(db, "PRAGMA table_info($(esc_id(table)))", sink)


text:
function svd_finite(A::Array{T,2}) where T<:ForwardDiff.Dual
    m = size(A,1)
    n = size(A,2)

    Av = map(ForwardDiff.value,A)
    num_A = length(Av)

    Uv,sv,Vv = svd(Av)

    dU = zeros(size(Uv,1),size(Uv,2),num_A)
    ds = zeros(length(sv),num_A)
    dV = zeros(size(Vv,1),size(Vv,2),num_A)

    ϵ = sqrt(eps(1.))
    for i = 1:num_A
        δ = zeros(size(Av))
        δ[i] = ϵ
        U,s,V = svd(Av + δ)
        dU[:,:,i] = (U - Uv) ./ ϵ
        ds[:,i] = (s - sv) ./ ϵ
        dV[:,:,i] = (V - Vv) ./ ϵ
    end

    Up = zeros(typeof(A[1].partials),size(Uv))
    sp = zeros(typeof(A[1].partials),size(sv))
    Vp = zeros(typeof(A[1].partials),size(Vv))
    for i = 1:size(Up,1)
        for j = 1:size(Up,2)
            Up[i,j] = sum([dU[i,j,k] * A[k].partials for k = 1:length(Av)])
        end
    end
    for i = 1:length(sp)
        sp[i] = sum([ds[i,k] * A[k].partials for k = 1:length(Av)])
    end
    for i = 1:size(Vp,1)
        for j = 1:size(Vp,2)
            Vp[i,j] = sum([dV[i,j,k] * A[k].partials for k = 1:length(Av)])
        end
    end

    U = map(T,Uv,Up)
    s = map(T,sv,sp)
    V = map(T,Vv,Vp)

    return U, s, V
end


text:
using Test, Random
function main()
    # OLS
    x = [ones(10,1) (1:10)]
    Random.seed!(1)
    y = rand(10,1)
    b = ols(y,x,silent=true)[1]
    @testset "ols" begin
        @test b[1,1] ≈ 0.07655377367377458
    end
    # samin
    xopt, fopt = samin()
    @testset "samin" begin
    @test xopt[1] ≈ 0.0 atol=1e-5
    @test fopt ≈ 2.0 atol=1e-5
    end
    # npreg
    @testset "npreg" begin
        yhat = npreg();
        @test yhat[1,1] ≈ -0.8108308248103517
    end
end
main()

max_stars_repo_name: manuelmontesinos/Econometrics

max_issues_repo_name: manuelmontesinos/Econometrics

max_forks_repo_name: manuelmontesinos/Econometrics


text:
function define_packed(ty::DataType)
    packed_offsets = cumsum([sizeof(x) for x in ty.types])
    sz = pop!(packed_offsets)
    pushfirst!(packed_offsets, 0)

    @eval begin
        function Base.read(io::IO, ::Type{$ty})
            bytes = read!(io, Array{UInt8}(undef, $sz...))
            hdr = $(Expr(:new, ty, [:(unsafe_load(convert(Ptr{$(ty.types[i])}, pointer(bytes)+$(packed_offsets[i])))) for i = 1:length(packed_offsets)]...,))
            if hdr.sizeof_hdr == ntoh(Int32(348))
                return byteswap(hdr), true
            end
            hdr, false
        end
        function Base.write(io::IO, x::$ty)
            bytes = UInt8[]
            for name in fieldnames($ty)
                append!(bytes, reinterpret(UInt8, [getfield(x,name)]))
            end
            write(io, bytes)
            $sz
        end
    end
    nothing
end

mutable struct NIfTI1Header
    sizeof_hdr::Int32

    data_type::NTuple{10,UInt8}
    db_name::NTuple{18,UInt8}
    extents::Int32
    session_error::Int16
    regular::Int8

    dim_info::Int8
    dim::NTuple{8,Int16}
    intent_p1::Float32
    intent_p2::Float32
    intent_p3::Float32
    intent_code::Int16
    datatype::Int16
    bitpix::Int16
    slice_start::Int16
    pixdim::NTuple{8,Float32}
    vox_offset::Float32
    scl_slope::Float32
    scl_inter::Float32
    slice_end::Int16
    slice_code::Int8
    xyzt_units::Int8
    cal_max::Float32
    cal_min::Float32
    slice_duration::Float32
    toffset::Float32

    glmax::Int32
    glmin::Int32

    descrip::NTuple{80,UInt8}
    aux_file::NTuple{24,UInt8}

    qform_code::Int16
    sform_code::Int16
    quatern_b::Float32
    quatern_c::Float32
    quatern_d::Float32
    qoffset_x::Float32
    qoffset_y::Float32
    qoffset_z::Float32

    srow_x::NTuple{4,Float32}
    srow_y::NTuple{4,Float32}
    srow_z::NTuple{4,Float32}

    intent_name::NTuple{16,UInt8}

    magic::NTuple{4,UInt8}
end
define_packed(NIfTI1Header)

# byteswapping

function byteswap(hdr::NIfTI1Header)
    for fn in fieldnames(typeof(hdr))
        val = getfield(hdr, fn)
        if isa(val, Number) && sizeof(val) > 1
            setfield!(hdr, fn, ntoh(val))
        elseif isa(val, NTuple) && sizeof(eltype(val)) > 1
            setfield!(hdr, fn, map(ntoh, val))
        end
    end
    hdr
end

"""
    freqdim(x)::Int

Returns the frequency dimension.
"""
freqdim(x::NIfTI1Header) = Int(getfield(x, :dim_info) & 0x03)
freqdim(x) = freqdim(header(x))

"""
    phasedim(x)::Int

Returns the phase dimension.
"""
phasedim(x::NIfTI1Header) = Int((getfield(x, :dim_info) >> 2) & 0x03)
phasedim(x) = phasedim(header(x))

"""
    slicedim(x)::Int

Returns the slice dimension.
"""
slicedim(x::NIfTI1Header) = Int((getfield(x, :dim_info) >> 4) & 0x03)
slicedim(x) = slicedim(header(x))

"""
    slice_start(x)::Int

Which slice corresponds to the first slice acquired during MRI acquisition (i.e. not padded slices).
"""
slice_start(x::NIfTI1Header) = Int(getfield(x, :slice_start)) + 1
slice_start(x) = slice_start(header(x))

"""
    slice_end(x)::Int

Which slice corresponds to the last slice acquired during MRI acquisition (i.e. not padded slices).
"""
slice_end(x::NIfTI1Header) = Int(getfield(x, :slice_end)) + 1
slice_end(x) = slice_end(header(x))

"""
    slice_duration(x)

Time to acquire one slice
"""
slice_duration(x::NIfTI1Header) = getfield(x, :slice_duration)
slice_duration(x) = slice_duration(header(x))


text:
println("Test misc")

@testset "misc" begin
  @testset "wsample_logprob" begin
    p = [.3, .5, .2]
    logp = log.(p)

    Random.seed!(1234)
    x = wsample_logprob(logp, 100_000)
    @test isapprox(mean(x .== 1), p[1], atol=1e-2)
    @test isapprox(mean(x .== 2), p[2], atol=1e-2)
    @test isapprox(mean(x .== 3), p[3], atol=1e-2)
  end

  @testset "quantiles" begin
    Random.seed!(0)
    X = randn(1_000_000, 5)
    @test all(isapprox.(quantiles(X, .025, dims=1), -1.96, atol=1e-2))
    @test all(isapprox.(quantiles(X, .975, dims=1), 1.96, atol=1e-2))
  end
end


text:
module MatrixOptInterface

using LinearAlgebra, SparseArrays
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

# export MatrixOptimizer

@enum ConstraintSense EQUAL_TO GREATER_THAN LESS_THAN INTERVAL
_sense(::Type{<:MOI.EqualTo}) = EQUAL_TO
_sense(::Type{<:MOI.LessThan}) = LESS_THAN
_sense(::Type{<:MOI.GreaterThan}) = GREATER_THAN
_sense(::Type{<:MOI.Interval}) = INTERVAL

_no_upper(bound) = bound != typemax(bound)
_no_lower(bound) = bound != typemin(bound)

function _bound_sense(lb::T, ub::T) where {T}
    if _no_upper(ub)
        if _no_lower(lb)
            if lb == ub
                return EQUAL_TO
            else
                return INTERVAL
            end
        else
            return LESS_THAN
        end
    else
        if _no_lower(lb)
            return GREATER_THAN
        else
            return
        end
    end
end

function _bound_set(lb::T, ub::T) where T
    if _no_upper(ub)
        if _no_lower(lb)
            if ub == lb
                return MOI.EqualTo(lb)
            else
                return MOI.Interval(lb, ub)
            end
        else
            return MOI.LessThan(ub)
        end
    else
        if _no_lower(lb)
            return MOI.GreaterThan(lb)
        else
            return
        end
    end
end

@enum VariableType CONTINUOUS INTEGER BINARY

include("sparse_matrix.jl")
include("conic_form.jl")
include("matrix_input.jl")
include("change_form.jl")

end


text:
const g = 9.80665

# ODE Functions
@taylorize function lander_thrust_on!(dx, x, p, t)
    # Unpack variables
    rx, ry, vx, vy, m_prop, Isp, m_dry, T, θ = x

    # Intermediate parameters
    m = m_dry + m_prop
    Fx, Fy = T .* (sin(θ), cos(θ))

    # Needed because numbers can't be converted into Taylor?
    _0 = zero(x[1])

    # State derivatives
    dx[1] = vx              # rx
    dx[2] = vy              # ry
    dx[3] = Fx/m            # vx
    dx[4] = -g + Fy/m       # vy
    dx[5] = -T/(Isp * g)    # m_prop
    dx[6] = _0              # Isp
    dx[7] = _0              # m_dry
    dx[8] = _0              # T
    dx[9] = _0              # θ
    return nothing
end

@taylorize function lander_thrust_off!(dx, x, p, t)
    # Unpack variables
    rx, ry, vx, vy, m_prop, Isp, m_dry, T = x

    # Needed because numbers can't be converted into Taylor?
    _1 = one(x[1])
    _0 = zero(x[1])

    # State derivatives
    dx[1] = vx      # rx
    dx[2] = vy      # ry
    dx[3] = _0      # vx
    dx[4] = -_1*g   # vy
    dx[5] = _0      # m_prop
    dx[6] = _0      # Isp
    dx[7] = _0      # m_dry
    dx[8] = _0      # T
    dx[9] = _0      # θ
    return nothing
end

@taylorize function lander_terminated!(dx, x, p, t)
    _0 = zero(x[1])
    dx[1] = _0  # rx
    dx[2] = _0  # ry
    dx[3] = _0  # vx
    dx[4] = _0  # vy
    dx[5] = _0  # m_prop
    dx[6] = _0  # Isp
    dx[7] = _0  # m_dry
    dx[8] = _0  # T
    dx[9] = _0  # θ
    return nothing
end


text:
function visualize(problem::FluidFlowProblem, quadrature::Quadrature, f_in, time, stats)
    return
    q = quadrature

    # Pre-allocate Macroscopic Variables
    ρ = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
    u = Array{Float64}(undef, size(f_in, 1), size(f_in, 2), dimension(q))
    p = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
    T = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
    σ_xx = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
    σ_xy = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))

    Nx = size(f_in, 1)
    Ny = size(f_in, 2)

    density_field = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
    velocity_field = Array{Float64}(undef, size(f_in, 1), size(f_in, 2), dimension(q))
    pressure_field = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
    σ_xx_field = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
    σ_xy_field = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))

    x_range, y_range = range(problem)

    f = Array{Float64}(undef, size(f_in, 3))
    u_ = zeros(dimension(q))
    @inbounds for x_idx in 1:Nx, y_idx in 1:Ny
        @inbounds for f_idx in 1:size(f_in, 3)
            f[f_idx] = f_in[x_idx, y_idx, f_idx]
        end

        x = x_range[x_idx]
        y = y_range[y_idx]

        density_field[x_idx, y_idx] =
            LatticeBoltzmann.density(quadrature, problem, x, y, time)
        pressure_field[x_idx, y_idx] =
            LatticeBoltzmann.pressure(quadrature, problem, x, y, time)
        velocity_field[x_idx, y_idx, :] = LatticeBoltzmann.velocity(problem, x, y, time)
        σ_exact = deviatoric_tensor(q, problem, x, y, time)
        σ_xx_field[x_idx, y_idx] = σ_exact[1, 1]
        σ_xy_field[x_idx, y_idx] = σ_exact[1, 2]

        ρ_ = density(q, f)
        velocity!(q, f, ρ_, u_)
        T_ = temperature(q, f, ρ_, u_)
        p_ = pressure(q, f, ρ_, u_)

        τ = q.speed_of_sound_squared * LatticeBoltzmann.lattice_viscosity(problem)
        σ_ = deviatoric_tensor(q, τ, f, ρ_, u_)

        ρ_ = dimensionless_density(problem, ρ_)
        u_ = dimensionless_velocity(problem, u_)
        T_ = dimensionless_temperature(q, problem, T_)
        p_ = dimensionless_pressure(q, problem, p_)
        ρ[x_idx, y_idx] = ρ_
        u[x_idx, y_idx, :] = u_
        p[x_idx, y_idx] = p_
        T[x_idx, y_idx] = T_

        σ_ = dimensionless_stress(problem, σ_)
        σ_xx[x_idx, y_idx] = σ_[1, 1]
        σ_xy[x_idx, y_idx] = σ_[1, 2]
    end

    s = (1000, 500)

    domain = (2:(problem.NY - 1)) ./ (problem.NY - 2)

    if (typeof(problem) != DecayingShearFlow)
        x_pos = round(Int, problem.NX / 2)
        domain = y_range[1:Ny]

        velocity_profile_x = plot(domain, u[x_pos, :, 1], label = "solution", title = "u_x")
        plot!(velocity_profile_x, domain, velocity_field[x_pos, :, 1], label = "exact")

        velocity_profile_y =
            plot(u[x_pos, 1:(problem.NY), 2], domain, label = "solution", title = "u_y")
        plot!(
            velocity_profile_y,
            velocity_field[x_pos, 1:(problem.NY), 2],
            domain,
            label = "exact",
        )

        pressure_profile =
            plot(domain, p[x_pos, 1:(problem.NY)], label = "solution", title = "p")
        plot!(
            pressure_profile,
            domain,
            pressure_field[x_pos, 1:(problem.NY)],
            label = "exact",
        )

        temperature_profile =
            plot(domain, T[x_pos, 1:(problem.NY), 1], label = "solution", title = "T")

        sigma_xx_profile = plot(
            domain,
            σ_xx[x_pos, 1:(problem.NY)],
            label = "solution",
            title = "sigma_xx",
        )
        plot!(sigma_xx_profile, domain, σ_xx_field[x_pos, 1:(problem.NY)], label = "exact")

        sigma_xy_profile = plot(
            domain,
            σ_xy[x_pos, 1:(problem.NY)],
            label = "solution",
            title = "sigma_xy",
        )
        plot!(sigma_xy_profile, domain, σ_xy_field[x_pos, 1:(problem.NY)], label = "exact")
    else
        y_pos = round(Int, problem.NY / 2)
        domain = x_range[1:Nx]

        velocity_profile_x = plot(domain, u[:, y_pos, 1], label = "solution", title = "u_x")
        plot!(velocity_profile_x, domain, velocity_field[:, y_pos, 1], label = "exact")

        velocity_profile_y = plot(u[:, y_pos, 2], domain, label = "solution", title = "u_y")
        plot!(velocity_profile_y, velocity_field[:, y_pos, 2], domain, label = "exact")

        pressure_profile = plot(domain, p[:, y_pos], label = "solution", title = "p")
        plot!(pressure_profile, domain, pressure_field[:, y_pos], label = "exact")

        temperature_profile = plot(domain, T[:, y_pos], label = "solution", title = "T")

        sigma_xx_profile =
            plot(domain, σ_xx[:, y_pos], label = "solution", title = "sigma_xx")
        plot!(sigma_xx_profile, domain, σ_xx_field[:, y_pos], label = "exact")

        sigma_xy_profile =
            plot(domain, σ_xy[:, y_pos], label = "solution", title = "sigma_xy")
        plot!(sigma_xy_profile, domain, σ_xy_field[:, y_pos], label = "exact")
    end

    plot(
        contour(ρ, title = "Density", fill = true, cbar = true),
        contour(p, title = "Pressure", fill = true),
        contour(T, title = "Temperature", fill = true),
        # contour(pressure_field', title="pressure analytical", fill=true),
        plot!(streamline(u), title = "Computed"),
        plot!(streamline(velocity_field), title = "exact"),
        plot(getfield.(stats, :error_u), legend = false, title = "U_e"),
        plot(getfield.(stats, :error_p), legend = false, title = "P_e"),
        plot(getfield.(stats, :error_σ_xx), legend = false, title = "sigma_xx_e"),
        plot(getfield.(stats, :error_σ_yx), legend = false, title = "sigma_yx_e"),
        plot(getfield.(stats, :error_σ_xy), legend = false, title = "sigma_xy_e"),
        plot(getfield.(stats, :error_σ_yy), legend = false, title = "sigma_yy_e"),
        # plot(getfield.(stats, :error_u), legend=false, title="U_e"),
        # plot(getfield.(stats, :error_p), legend=false, title="P_e"),
        # kinetic_energy_profile,
        pressure_profile,
        temperature_profile,
        velocity_profile_x,
        velocity_profile_y,
        sigma_xx_profile,
        sigma_xy_profile,
        size = (1000, 600),
        legend = nothing,
    )
    gui()
end

function streamline(
    j;
    amount_of_arrows = 5,
    step = round(Int, size(j, 1) / amount_of_arrows),
)
    s = (1000, 500)
    velocity_field =
        contour((j[:, :, 1] .^ 2 .+ j[:, :, 2] .^ 2)', fill = true, title = "Momentum")
    N = size(j, 1)
    X = [
        i
        for i in range(2, size(j, 1), step = step), j in range(1, size(j, 2), step = step)
    ]
    Y = [
        j
        for i in range(2, size(j, 1), step = step), j in range(1, size(j, 2), step = step)
    ]

    quiver!(
        velocity_field,
        X,
        Y,
        quiver = (x, y) -> (
            0.1 * amount_of_arrows * j[Int(x), Int(y), 1] /
                sqrt(j[Int(x), Int(y), 1]^2 + j[Int(x), Int(y), 2]^2),
            0.1 * amount_of_arrows * j[Int(x), Int(y), 2] /
                sqrt(j[Int(x), Int(y), 1]^2 + j[Int(x), Int(y), 2]^2),
        ),
        color = "white",
    )

    return velocity_field
end


text:
module OrdinalRegr

using Convex, SCS

struct OrdinalSVM
    w
    b
    L
end

function fitordinalsvm(y, X, L; C=1e-2)
    n, p = size(X)
    l = length(L)
    w = Variable(p)
    b = Variable(l - 1)
    # get indexes of the labels
    inds = convert(Vector{Int}, indexin(y, L))
    # compute objective
    obj = C * sum(w .* w)
    f = X * w
    for (j, lj) in enumerate(L[1:end-1])
        for (i, yi) in enumerate(y)
            obj += yi ≤ lj ? max(0.0, 1.0 + (f[i] - b[j])) : max(0.0, 1.0 - (f[i] - b[j]))
        end
    end
    problem = minimize(obj, (b[i] <= b[i+1] for i in 1:l-2)...)
    solve!(problem, () -> SCS.Optimizer(verbose=false))
    return OrdinalSVM(w.value[:], b.value[:], L)
end

fitordinalsvm(y, X; C=1e-2) = fitordinalsvm(y, X, sort!(unique(y)); C=1e-2)

"""
    score(model::OrdinalSVM, X::AbstractMatrix)

Compute the scoring using an `OrdinalSVM` model.
"""
function score(model::OrdinalSVM, X::AbstractMatrix)
    X * model.w
end

"""
predict(model::OrdinalSVM, X::AbstractMatrix)

Predict the label using an `OrdinalSVM` model.
"""
function predict(model::OrdinalSVM, X::AbstractMatrix)
    l = [model.L[1] for i in 1:size(X, 1)]
    for (i, s) in enumerate(score(model, X))
        for (Lj, bj) in zip(model.L[2:end], model.b)
            s - bj ≥ 0 && (l[i] = Lj)
        end
    end
    return l

end

end


text:
abstract type Dataset end

abstract type Mode end
struct Train <: Mode end
struct Dev <: Mode end
struct Test <: Mode end

function dataset( ::Type{M}, d::D, args...; kwargs...) where M <: Mode where D <: Dataset
    ds = datafile(M, d, args...; kwargs...)
    rds = reader.(ds)
    rds
end

datafile(::Type{Train}, d::D, args...; kwargs...) where D <: Dataset = trainfile(d, args...; kwargs...)
datafile(::Type{Dev}, d::D, args...; kwargs...) where D <: Dataset = devfile(d, args...; kwargs...)
datafile(::Type{Test}, d::D, args...; kwargs...) where D <: Dataset = testfile(d, args...; kwargs...)

trainfile(d, args...; kwargs...) = (println("Trainset not found"); nothing)
devfile(d, args...; kwargs...) = (println("Devset not found"); nothing)
testfile(d, args...; kwargs...) = (println("Testset not found"); nothing)

function reader(file::AbstractString)
    ch = Channel{String}(0)
    task = @async foreach(x->put!(ch, x), eachline(open(file)))
    bind(ch, task)
    ch
end

function reader(iter)
    ch = Channel{String}(0)
    task = @async foreach(x->put!(ch, x), iter)
    bind(ch, task)
    ch
end

function batched(xs)
    s = length(xs)
    sx = s != 0 ? length(xs[1]) : 0
    res = [Vector{typeof(xs[1][i])}() for i = 1:sx]
    for x ∈ xs
        for (i, xi) ∈ enumerate(x)
            push!(res[i], xi)
        end
    end
    res
end

function get_batch(c::Channel, n=1)
    res = Vector(undef, n)
    for (i, x) ∈ enumerate(c)
        res[i] = x
        i >= n && break
    end
    isassigned(res, n) ? batched(res) : nothing
end

function get_batch(cs::Container{C}, n=1) where C <: Channel
    res = Vector(undef, n)
    for (i, xs) ∈ enumerate(zip(cs...))
        res[i] = xs
        i >= n && break
    end
    isassigned(res, n) ? batched(res) : nothing
end

get_vocab(::D, args...; kwargs...) where D <: Dataset = (println("No prebuild vocab"); nothing)

function token_freq(files...; vocab::Dict{String, Int} = Dict{String,Int}(), min_freq::Int = 3)
    for f ∈ files
        open(f) do fd
            for line ∈ eachline(fd)
                for token ∈ tokenize(line)
                    token = intern(token)
                    vocab[token] = get(vocab, token, 0) + 1
                end
            end
        end
    end

    for key ∈ keys(vocab)
        if vocab[key] < min_freq
            delete!(vocab, key)
        end
    end
    vocab
end


text:
abstract type AbstractCallback end
abstract type SafeCallback <: AbstractCallback end
const Callback = SafeCallback
abstract type UnsafeCallback <: AbstractCallback end

abstract type AbstractMetric <: Callback end
abstract type AbstractLogger <: Callback end

stateaccess(::Callback) = (;)
runafter(::AbstractCallback) = ()

abstract type ConflictResolution end
struct NotDefined <: ConflictResolution end
struct Unresolvable <: ConflictResolution end
struct RunFirst <: ConflictResolution cb end
struct NoConflict <: ConflictResolution end

function _resolveconflict(cb1, cb2)
    r = resolveconflict(cb1, cb2)
    if r === NotDefined()
        return resolveconflict(cb2, cb1)
    else
        return r
    end
end
resolveconflict(::AbstractCallback, ::AbstractCallback) = NotDefined()

abstract type FitException <: Exception end
struct CancelBatchException <: FitException
    msg::String
end
struct CancelEpochException <: FitException
    msg::String
end
struct CancelFittingException <: FitException
    msg::String
end

# Callback hook

"""
    on(event::FitEvent, phase::AbstractFittingPhase, callback::AbstractCallback, learner)

Handle `event` with `callback`. Can dispatch on an `AbstractFittingPhase` and
receives `learner` as an additional argument.

If not overwritten with a more specific method, does nothing.

To see events which an `AbstractCallback` handles, use

    `methods(Training.on, (Any, Any, MyCallbackType, Any)`
"""
on(::FitEvent, ::AbstractFittingPhase, ::AbstractCallback, learner) = return

_on(e, p, cb, learner) = on(e, p, cb, learner)
_on(e, p, cb::SafeCallback, learner) = on(e, p, cb, protect(learner, stateaccess(cb)))
#_on(e, p, cb::SafeCallback, learner) = on(e, p, cb, learner)


text:
# ____________________________________________________________________
# ironuv
#
function ironuv_read()
    ia   = readdlm(qsfit_data() * "/VW2001/Fe_UVtemplt_A.asc")
    ib   = readdlm(qsfit_data() * "/VW2001/Fe_UVtemplt_B.asc")
    i34  = readdlm(qsfit_data() * "/VW2001/Fe3UV34modelB2.asc")
    i47  = readdlm(qsfit_data() * "/VW2001/Fe3_UV47.asc")
    i191 = readdlm(qsfit_data() * "/VW2001/Fe2_UV191.asc")

    lmin, lmax = extrema([ia[:,1]; ib[:,1]; i34[:,1]; i47[:,1]; i191[:,1]])
    lmax = 3.3e3 # avoid overshooting to negative values

    λ = lmin:(ia[2,1]-ia[1,1]):lmax
    ia   = Spline1D(  ia[:,1],   ia[:,2], k=1, bc="zero")(λ)
    ib   = Spline1D(  ib[:,1],   ib[:,2], k=1, bc="zero")(λ)
    i34  = Spline1D( i34[:,1],  i34[:,2], k=1, bc="zero")(λ)
    i47  = Spline1D( i47[:,1],  i47[:,2], k=1, bc="zero")(λ)
    i191 = Spline1D(i191[:,1], i191[:,2], k=1, bc="zero")(λ)

    L = ib .+ i47 .+ i191
    L ./= sum(L .* (λ[2] - λ[1]))
    return (λ, L)
end

mutable struct ironuv <: AbstractComponent
    norm::Parameter
    fwhm::Float64
    L::Vector{Float64}  # pre-computed

    function ironuv(fwhm::Number)
        out = new(Parameter(1), float(fwhm), Vector{Float64}())
        out.norm.low = 0
        return out
    end
end

function prepare!(comp::ironuv, domain::Domain{1})
    @assert comp.fwhm > 900
    λ, L = QSFit.ironuv_read()
    L = QSFit.conv_gauss(λ, L, sqrt(comp.fwhm^2. - 900^2) / 2.355)
    L ./= int_tabulated(λ, L)[1]
    comp.L = Spline1D(λ, L, k=1, bc="zero")(domain[:])
    return fill(NaN, length(domain))
end

function evaluate!(buffer, comp::ironuv, domain::Domain{1},
                   norm)
    buffer .= norm .* comp.L
end

#=
    x = Domain(1e3:1:3.5e3)
    comp = QSFit.ironuv(3000.)
    ceval = GFit.CompEval(comp, x)
    GFit.evaluate_cached(ceval)
    @gp x[:] ceval.buffer ./ maximum(ceval.buffer) "w l"
=#


text:
# Functions for text printing (setting font ect. is done in setup.jl

function text(gr::Graph, x::Real, y::Real, z::Real, txt::String; fnt::String=":AC", size::Real=-1.)
    cmgl.puts(gr.ptr, x, y, z, txt, fnt, size)
end

function text(gr::Graph, x::Real, y::Real, txt::String; fnt::String=":AC", size::Real=-1.)
    cmgl.puts(gr.ptr, x, y, 0., txt, fnt, size)
end

function text(gr::Graph, p::AbstractArray, txt::String; fnt::String=":AC", size::Real=-1.)
    if length(p) == 3
        cmgl.puts(gr.ptr, p..., txt, fnt, size)
    elseif length(p) == 2
        cmgl.puts(gr.ptr, p..., 0., txt, fnt, size)
    else
        error("Lengths of the input vectors must match (2 or 3 dimensional)")
    end
end

#=function text(gr::Graph, x::Real, y::Real, z::Real, txt::UTF8String; fnt::String=":AC", size::Real=-1.)=#
    #=cmgl.putsw(gr.ptr, x, y, z, txt, fnt, size)=#
#=end=#

#=function text(gr::Graph, x::Real, y::Real, txt::UTF8String; fnt::String=":AC", size::Real=-1.)=#
    #=cmgl.putsw(gr.ptr, x, y, 0., txt, fnt, size)=#
#=end=#

#=function text(gr::Graph, p::AbstractArray, txt::UTF8String; fnt::String=":AC", size::Real=-1.)=#
    #=if length(p) == 3=#
        #=cmgl.putsw(gr.ptr, p..., txt, fnt, size)=#
    #=elseif length(p) == 2=#
        #=cmgl.putsw(gr.ptr, p..., 0., txt, fnt, size)=#
    #=else=#
        #=error("Lengths of the input vectors must match (2 or 3 dimensional)")=#
    #=end=#
#=end=#

function text(gr::Graph, x::Real, y::Real, z::Real, dx::Real, dy::Real, dz::Real, txt::String; fnt::String=":AC", size::Real=-1.)
    cmgl.puts_dir(gr.ptr, x, y, z, dx, dy, dz, txt, fnt, size)
end

function text(gr::Graph, x::Real, y::Real, dx::Real, dy::Real, txt::String; fnt::String=":AC", size::Real=-1.)
    cmgl.puts_dir(gr.ptr, x, y, 0., dx, dy, 0., txt, fnt, size)
end

function text(gr::Graph, p::AbstractArray, dp::AbstractArray, txt::String; fnt::String=":AC", size::Real=-1.)
    if length(p) == length(dp) == 3
        cmgl.puts_dir(gr.ptr, p..., dp..., txt, fnt, size)
    elseif length(2) == length(dp) == 2
        cmgl.puts_dir(gr.ptr, p..., 0., dp..., 0., txt, fnt, size)
    else
        error("Lengths of the input vectors must match (2 or 3 dimensional)")
    end
end

#=function text(gr::Graph, x::Real, y::Real, z::Real, dx::Real, dy::Real, dz::Real, txt::UTF8String; fnt::String=":AC", size::Real=-1.)=#
    #=cmgl.putsw_dir(gr.ptr, x, y, z, dx, dy, dz, txt, fnt, size)=#
#=end=#

#=function text(gr::Graph, x::Real, y::Real, dx::Real, dy::Real, txt::UTF8String; fnt::String=":AC", size::Real=-1.)=#
    #=cmgl.putsw_dir(gr.ptr, x, y, 0., dx, dy, 0., txt, fnt, size)=#
#=end=#

#=function text(gr::Graph, p::AbstractArray, dp::AbstractArray, txt::UTF8String; fnt::String=":AC", size::Real=-1.)=#
    #=if length(p) == length(dp) == 3=#
        #=cmgl.putsw_dir(gr.ptr, p..., dp..., txt, fnt, size)=#
    #=elseif length(2) == length(dp) == 2=#
        #=cmgl.putsw_dir(gr.ptr, p..., 0., dp..., 0., txt, fnt, size)=#
    #=else=#
        #=error("Lengths of the input vectors must match (2 or 3 dimensional)")=#
    #=end=#
#=end=#

#TODO: Strings along curves


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   J2 osculating orbit propagator algorithm.
#
#       This algorithm propagates the orbit considering the secular and
#       short-period perturbations introduced by the J2 gravitational term. The
#       algorithm is based on Kwok version as indicated in [1, p. 708-710].
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export j2osc_init, j2osc!

################################################################################
#                                  Functions
################################################################################

"""
    j2osc_init(epoch::Tepoch, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number, dn_o2::Number, ddn_o6::Number; j2_gc::J2_GravCte{T} = j2_gc_egm08) where {Tepoch, T}

Initialize the data structure of J2 osculating orbit propagator algorithm.

!!! note
    The type used in the propagation will be the same as used to define the
    gravitational constants in the structure `j2_gc`.

# Args

- `epoch::Number`: Epoch of the initial mean orbital elements [Julian Day].
- `a_0::Number`: Initial mean semi-major axis [m].
- `e_0::Number`: Initial mean eccentricity.
- `i_0::Number`: Initial mean inclination [rad].
- `Ω_0::Number`: Initial mean right ascension of the ascending node [rad].
- `ω_0::Number`: Initial mean argument of perigee [rad].
- `f_0::Number`: Initial mean true anomaly [rad].
- `dn_o2::Number`: First time derivative of the mean motion divided by two
    [rad/s^2].
- `ddn_o6::Number`: Second time derivative of the mean motion divided by six
    [rad/s^3].

# Keywords

- `j2_gc::J2_GravCte`: J2 orbit propagator gravitational constants (see
    [`J2_GravCte`](@ref)). (**Default** = `j2_gc_egm08`)

# Returns

The structure [`J2osc_Structure`](@ref) with the initialized parameters.

# Remarks

The inputs are the mean orbital elements.
"""
function j2osc_init(
    epoch::Tepoch,
    a_0::Number,
    e_0::Number,
    i_0::Number,
    Ω_0::Number,
    ω_0::Number,
    f_0::Number,
    dn_o2::Number,
    ddn_o6::Number;
    j2_gc::J2_GravCte{T} = j2_gc_egm08
) where {Tepoch, T}
    j2d = j2_init(epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, dn_o2, ddn_o6; j2_gc = j2_gc)

    # Initialize the structure.
    j2oscd = J2osc_Structure{Tepoch, T}(j2d, 0, 0, 0, 0, 0, 0, 0, 0)

    # Call the propagation one first time to update the osculating elements.
    j2osc!(j2oscd, 0)

    return j2oscd
end

"""
    j2osc!(j2oscd::J2osc_Structure{Tepoch, T}, t::Number) where {Tepoch, T}

Propagate the orbit defined in `j2oscd` (see [`J2osc_Structure`](@ref)) until
the time `t` [s].

!!! note
    The internal values in `j2oscd` will be modified.

# Returns

- The position vector represented in the inertial frame at time `t` [m].
- The velocity vector represented in the inertial frame at time `t` [m/s]

# Remarks

The inertial frame in which the output is represented depends on which frame it
was used to generate the orbit parameters. Notice, however, that the
perturbation theory requires an inertial frame with true equator.
"""
function j2osc!(j2oscd::J2osc_Structure{Tepoch, T}, t::Number) where {Tepoch, T}
    # First, we need to propagate the mean elements since they are necessary to
    # compute the short-periodic perturbations.
    j2d = j2oscd.j2d
    j2!(j2d, t)

    @unpack_J2_GravCte j2d.j2_gc

    # Obtain the mean elements at this time instant.
    al_k = j2d.al_k
    e_k  = j2d.e_k
    e_k² = e_k * e_k
    i_k  = j2d.i_k
    Ω_k  = j2d.Ω_k
    ω_k  = j2d.ω_k
    M_k  = j2d.M_k
    f_k  = j2d.f_k
    p_k  = al_k * R0 * (1 - e_k²)
    p_k² = p_k * p_k
    u_k  = ω_k + f_k

    # Auxiliary variables to reduce the computational burden.
    KJ2 = J2 * R0 * R0
    sin_i_k, cos_i_k = sincos(i_k)
    sin_f_k, cos_f_k = sincos(f_k)
    sin_2u_k, cos_2u_k = sincos(2u_k)
    sin_2ω_f_k, cos_2ω_f_k = sincos(2ω_k + f_k)
    sin_2ω_3f_k, cos_2ω_3f_k = sincos(2ω_k + 3f_k)
    sin_i_k² = sin_i_k * sin_i_k
    cos_i_k² = cos_i_k * cos_i_k
    e_cos_f_k = e_k * cos_f_k
    e_sin_f_k = e_k * sin_f_k

    aux1 = 3cos(2u_k) + 3e_k * cos_2ω_f_k + e_k * cos_2ω_3f_k
    aux2 = sqrt(1 - e_k²)
    aux3 = 3cos_i_k² - 1

    # Compute the short-periodic perturbations considering only the J2
    # gravitational term.
    δisp_k = +KJ2 * sin_i_k * cos_i_k / (4p_k²) * aux1

    δpsp_k = +KJ2 * sin_i_k² / (2p_k) * aux1

    δΩsp_k = -KJ2 * cos_i_k / (4p_k²) * (
        6 * (f_k - M_k + e_sin_f_k) - 2sin_2u_k - 3e_k * sin_2ω_f_k - e_k * sin_2ω_3f_k
    )

    δrsp_k = -KJ2 / (4p_k) * (
        aux3 * (2aux2 / (1 + e_cos_f_k) + e_cos_f_k / (1 + aux2)) - sin_i_k² * cos_2u_k
    )

    δṙsp_k = +KJ2 * sqrt(μm) / (4p_k^T(2.5)) * (
        aux3 * e_sin_f_k * (aux2 + (1 + e_cos_f_k)^2 / (1 + aux2)) -
        sin_i_k² * (1 - e_cos_f_k)^2 * sin_2u_k
    )

    δusp_k = +KJ2 / (8p_k²) * (
        (6 - 30cos_i_k²) * (f_k - M_k) +
        4e_k * sin_f_k * ((1 - 6cos_i_k²) - aux3 / (1 + aux2)) -
        aux3 / (1 + aux2) * e_k² * sin(2f_k) +
        (5cos_i_k² - 2) * (2e_k) * sin_2ω_f_k +
        (7cos_i_k² - 1) * sin_2u_k +
        2cos_i_k² * e_k * sin_2ω_3f_k
    )

    r_k = p_k / (1 + e_k * cos_f_k)
    ṙ_k = sqrt(μm / p_k) * e_k * cos_f_k

    r_osc_k = r_k + δrsp_k
    ṙ_osc_k = ṙ_k + δṙsp_k
    p_osc_k = p_k + δpsp_k

    A_k = p_osc_k / r_osc_k - 1
    B_k = sqrt(p_osc_k / μm) * ṙ_osc_k

    e_osc_k² = A_k^2 + B_k^2
    e_osc_k  = sqrt(e_osc_k²)
    a_osc_k  = p_osc_k / (1 - e_osc_k²)
    i_osc_k  = i_k + δisp_k
    Ω_osc_k  = Ω_k + δΩsp_k
    u_osc_k  = u_k + δusp_k
    f_osc_k  = atan(B_k, A_k)
    ω_osc_k  = u_osc_k - f_osc_k
    M_osc_k  = f_to_M(e_osc_k, f_osc_k)

    # Compute the position and velocity considering the osculating elements.
    r_i_k, v_i_k = kepler_to_rv(a_osc_k, e_osc_k, i_osc_k, Ω_osc_k, ω_osc_k, f_osc_k)

    # Update the J2 orbit propagator structure.
    j2oscd.Δt  = T(t)
    j2oscd.a_k = a_osc_k
    j2oscd.e_k = e_osc_k
    j2oscd.i_k = i_osc_k
    j2oscd.Ω_k = Ω_osc_k
    j2oscd.ω_k = ω_osc_k
    j2oscd.M_k = M_osc_k
    j2oscd.f_k = f_osc_k

    return r_i_k, v_i_k
end


text:
using Test
using LinearAlgebra
using UnPack
using MOONS.Grids
using MOONS.Fields

#####
##### CellCenter <-> CellFace
#####

function interp_simple!(cent::CellCenter, face::CellFace{1}, grid::Grid)
    for k in over_points(grid.c[3], All(), Center1D())
        for j in over_points(grid.c[2], All(), Center1D())
            for i in over_points(grid.c[1], All(), Center1D())
                D = 1//2
                U = 1//2
                cent[i, j, k] = D*face[i, j, k] + U*face[i+1, j, k]
            end
        end
    end
end

function interp_simple!(face::CellFace{1}, cent::CellCenter, grid::Grid)
    face .= 0 # boundaries are not populated.
    for k in over_points(grid.c[3], All(), Center1D())
        for j in over_points(grid.c[2], All(), Center1D())
            for i in over_points(grid.c[1], Boundary(), Vertex1D())
                D = grid.c[1].interp.diag[i-1]
                U = grid.c[1].interp.upper[i-1]
                face[i, j, k] = D*cent[i-1, j, k] + U*cent[i, j, k]
            end
        end
    end
end

@testset "Interpolations - Simple implementation Center <-> Face" begin
    FT = Float64
    n, a, b, β = 10, FT(0.0), FT(1.0), FT(1.01)
    grid = Cube(a, b, n; warpfun=Roberts_both, args=(β,))

    # Face -> Center
    face = CellFace(grid, 1)
    cent_correct = CellCenter(grid)
    cent = CellCenter(grid)
    @inbounds for local_grid in GridField(grid,face)
        @unpack ijk,x,y,z = local_grid
        face[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
    end
    interp_simple!(cent_correct, face, grid)
    interp!(cent, face, grid)
    @test all(cent .≈ cent_correct)

    # Center -> Face
    cent = CellCenter(grid)
    face_correct = CellFace(grid, 1)
    face = CellFace(grid, 1)
    @inbounds for local_grid in GridField(grid,cent)
        @unpack ijk,x,y,z = local_grid
        cent[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
    end
    interp_simple!(face_correct, cent, grid)
    interp!(face, cent, grid)
    @test all(face .≈ face_correct)
end

#####
##### CellCorner <-> CellEdge
#####

function interp_simple!(corn::CellCorner, edge::CellEdge{1}, grid::Grid)
    corn .= 0 # boundaries are not populated.
    for k in over_points(grid.c[3], All(), Vertex1D())
        for j in over_points(grid.c[2], All(), Vertex1D())
            for i in over_points(grid.c[1], Boundary(), Vertex1D())
                D = grid.c[1].interp.diag[i-1]
                U = grid.c[1].interp.upper[i-1]
                corn[i, j, k] = D*edge[i-1, j, k] + U*edge[i, j, k]
            end
        end
    end
end

function interp_simple!(edge::CellEdge{1}, corn::CellCorner, grid::Grid)
    for k in over_points(grid.c[3], All(), Vertex1D())
        for j in over_points(grid.c[2], All(), Vertex1D())
            for i in over_points(grid.c[1], All(), Center1D())
                D = 1//2
                U = 1//2
                edge[i, j, k] = D*corn[i, j, k] + U*corn[i+1, j, k]
            end
        end
    end
end

@testset "Interpolations - Simple implementation Corner <-> Edge" begin
    FT = Float64
    n, a, b, β = 10, FT(0.0), FT(1.0), FT(1.01)
    grid = Cube(a, b, n; warpfun=Roberts_both, args=(β,))

    # Edge -> Corner
    edge = CellEdge(grid, 1)
    corn_correct = CellCorner(grid)
    corn = CellCorner(grid)
    @inbounds for local_grid in GridField(grid,edge)
        @unpack ijk,x,y,z = local_grid
        edge[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
    end
    interp_simple!(corn_correct, edge, grid)
    interp!(corn, edge, grid)
    @test all(corn .≈ corn_correct)

    # corner -> edge
    corn = CellCorner(grid)
    edge_correct = CellEdge(grid, 1)
    edge = CellEdge(grid, 1)
    @inbounds for local_grid in GridField(grid,corn)
        @unpack ijk,x,y,z = local_grid
        corn[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
    end
    interp_simple!(edge_correct, corn, grid)
    interp!(edge, corn, grid)
    @test all(edge .≈ edge_correct)
end

#####
##### CellFace <-> CellEdge
#####

function interp_simple!(face::CellFace{1}, edge::CellEdge{2}, grid::Grid)
    for k in over_points(grid.c[3], All(), Center1D())
        for j in over_points(grid.c[2], All(), Center1D())
            for i in over_points(grid.c[1], All(), Vertex1D())
                D = 1//2
                U = 1//2
                face[i, j, k] = D*edge[i, j, k] + U*edge[i, j, k+1]
            end
        end
    end
end

function interp_simple!(edge::CellEdge{2}, face::CellFace{1}, grid::Grid)
    edge .= 0 # boundaries are not populated.
    for k in over_points(grid.c[3], Boundary(), Vertex1D())
        for j in over_points(grid.c[2], All(), Center1D())
            for i in over_points(grid.c[1], All(), Vertex1D())
                D = grid.c[3].interp.diag[k-1]
                U = grid.c[3].interp.upper[k-1]
                edge[i, j, k] = D*face[i, j, k-1] + U*face[i, j, k]
            end
        end
    end
end

@testset "Interpolations - Simple implementation Face <-> Edge" begin
    FT = Float64
    n, a, b, β = 10, FT(0.0), FT(1.0), FT(1.01)
    grid = Cube(a, b, n; warpfun=Roberts_both, args=(β,))

    # Edge -> Face
    edge = CellEdge(grid, 2)
    face_correct = CellFace(grid, 1)
    face = CellFace(grid, 1)
    @inbounds for local_grid in GridField(grid,edge)
        @unpack ijk,x,y,z = local_grid
        edge[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
    end
    interp_simple!(face_correct, edge, grid)
    interp!(face, edge, grid)
    @test all(face .≈ face_correct)

    # Face -> edge
    face = CellFace(grid, 1)
    edge_correct = CellEdge(grid, 2)
    edge = CellEdge(grid, 2)
    @inbounds for local_grid in GridField(grid,face)
        @unpack ijk,x,y,z = local_grid
        face[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
    end
    interp_simple!(edge_correct, face, grid)
    interp!(edge, face, grid)
    @test all(edge .≈ edge_correct)
end

#####
##### nested / non-base interpolations
#####

@testset "Interpolations - nested / non-base interps" begin
    FT = Float64
    n, a, b, β = 10, FT(0.0), FT(1.0), FT(1.01)
    grid = Cube(a, b, n; warpfun=Roberts_both, args=(β,))

    edge = (CellEdge(grid, 1), CellEdge(grid, 2), CellEdge(grid, 3))
    face = (CellFace(grid, 1), CellFace(grid, 2), CellFace(grid, 3))
    cent = CellCenter(grid)
    corn = CellCorner(grid)
    @inbounds for local_grid in GridField(grid, cent)
        @unpack ijk,x,y,z = local_grid
        cent[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
    end

    # Just making sure these are callable:
    interp!(face[1], edge[1], grid)
    interp!(edge[1], face[1], grid)

    interp!(cent, edge[1], grid)
    interp!(edge[1], cent, grid)

    interp!(corn, face[1], grid)
    interp!(face[1], corn, grid)

    interp!(cent, corn, grid)
    interp!(corn, cent, grid)
end


text:
using QuantEcon

B = 10
M = 5
alpha = 0.5
beta = 0.9
u(c) = c^alpha
n = B + M + 1
m = M + 1

s_indices = Int64[]
a_indices = Int64[]
Q = Array{Float64}(0, n)
R = Float64[]

b = 1.0 / (B + 1)

for s in 0:(M + B)
    for a in 0:min(M, s)
        s_indices = [s_indices; s + 1]
        a_indices = [a_indices; a + 1]
        q = zeros(Float64, 1, n)
        q[(a + 1):((a + B) + 1)] = b
        Q = [Q; q]
        R = [R; u(s-a)]
    end
end

ddp = DiscreteDP(R, Q, beta, s_indices, a_indices);
results = solve(ddp, PFI)


text:
include("../common.jl")
build_julia(ARGS, v"1.4.2")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using Distributions, ForwardDiff, Plots, LinearAlgebra

include(string(pwd(),"/hmc/rkhs.jl"))

replics = 3000

x = randn(500)

y = 1.5*x + randn(500)

prior = Normal(0.5,1)

likelihood(y,x,par) = sum(logpdf.(Normal(0,1),y - x .*par))

loglike(y,x,par) = likelihood(y,x,par)+logpdf.(prior,par)

unif = Uniform(0,1)

moment_dist = Normal(0,2)

N = 200

dados = zeros(N,2)

for i in 1:N
    xx = rand(prior)
    dados[i,:] = [xx;loglike(y,x,xx)]
end

sig = 1

rkhs_approx = rkhs_gaussian(dados[:,1],dados[:,2],sig,0)

f_hat(x) = rkhs_approx(x)

dif_aprox(x) = ForwardDiff.derivative(f_hat,x)
dif_exact(xx) = ForwardDiff.derivative(z->loglike(y,x,z),xx)

gridd = -2:0.05:2

diff_exact_vals = dif_exact.(gridd)
diff_approx_vals = dif_aprox.(gridd)

plot(gridd,diff_approx_vals)
plot!(gridd,diff_exact_vals)

par = zeros(replics)
par[1] = rand(prior)

moments = zeros(replics,500)

for i in 2:replics
    moment = rand(moment_dist)
    old_moment = moment
    par_aux = par[i-1]

    L = ceil(rand(Uniform(100,500)))
    leapfrog_step = 1/L

    for j in 1:L
        new_moment = moment + 1/2*leapfrog_step*ForwardDiff.derivative(pp->f_hat(pp),par_aux)
         par_aux = par_aux .+ leapfrog_step*1*new_moment
         moment = new_moment + 1/2*leapfrog_step*ForwardDiff.derivative(pp->f_hat(pp),par_aux)
         moments[i,Int64(j)] = new_moment
    end

    num = loglike(y,x,par_aux) + logpdf(moment_dist,moment)#+ logpdf(kernel_novo,par[i-1])

    dem = loglike(y,x,par[i-1]) + logpdf(moment_dist,old_moment)#+ logpdf(kernel,novo_par)
    alpha = exp(min(0,num - dem))
    p = rand(unif)
    if alpha < p
        par[i] = par[i-1]
    else
        par[i] = par_aux
    end
    println(i)
    if i % 400 == 0
        novos_dados = sample(par[1:i],N)
        novas_ll = [loglike(y,x,novos_dados[j]) for j in 1:N]
        rkhs_approx = rkhs_gaussian(novos_dados,novas_ll,sig,0)
        global f_hat(x) = rkhs_approx(x)
        dif_aprox(x) = ForwardDiff.derivative(f_hat,x)
        dif_exact(xx) = ForwardDiff.derivative(z->loglike(y,x,z),xx)

        gridd = -2:0.05:2

        diff_exact_vals = dif_exact.(gridd)
        diff_approx_vals = dif_aprox.(gridd)
        #println("Dif Error", sqrt(sum((diff_exact_vals - diff_approx_vals).^2)))
    end
end

histogram(par[100:3000])

par_rwmh = zeros(replics)
par_rwmh[1] = par[1]

for i in 2:replics
    kernel = Normal(par_rwmh[i-1],0.5)
    novo_par = rand(kernel)
    kernel_novo = Normal(novo_par,0.5)

    num = likelihood(y,x,novo_par) + logpdf(prior,novo_par) + logpdf(kernel_novo,par_rwmh[i-1])

    dem = likelihood(y,x,par_rwmh[i-1]) + logpdf(prior,par_rwmh[i-1]) + logpdf(kernel,novo_par)
    alpha = exp(min(0,num - dem))
    p = rand(unif)
    if alpha < p
        par_rwmh[i] = par_rwmh[i-1]
    else
        par_rwmh[i] = novo_par
    end
    println(i)
end

using StatsPlots

StatsPlots.density(par[100:1000])
StatsPlots.density!(par_rwmh[100:1000])
vline!([1.5])

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using Test;

function run(fname)
    # read fil
    f = open(fname)
    lines = readlines(f)
    close(f)

    # parse input
    earliest_depart = parse(Int, lines[1])
    nums = split(lines[2], ",")
    filter!(x -> x != "x", nums)
    bus_ids = map(x -> parse(Int, x), nums)

    println("Earliest departure = ", earliest_depart)
    println("Bus IDs = ", bus_ids)

    result_bus = -1
    result_wait = +Inf
    for bid in bus_ids
        possible_departure = (div(earliest_depart, bid) + 1) * bid
        possible_wait_time = possible_departure - earliest_depart
        println("Possible = ", possible_wait_time, " (Bus ID = ", bid, " Departure = ", possible_departure, ")")
        if possible_wait_time < result_wait
            result_bus = bid
            result_wait = possible_wait_time
        end
    end

    result = result_bus * result_wait
    println("Result = ", result, " (Bus ID = ", result_bus, ", Wait Time = ", result_wait, ")")

    return result
end

# 59*5 = 295
@test run("13ex.txt") == 295
run("13.txt")

# Chinese Remainder Theorem
# https://en.wikipedia.org/wiki/Chinese_remainder_theorem#Computation
# I first tried to use mods.jl CRT() function but hit error (https://github.com/scheinerman/Mods.jl)
# Now following approach here https://www.freecodecamp.org/news/how-to-implement-the-chinese-remainder-theorem-in-java-db88a3f1ffe0/
function CRT(nums, rems)
    # Step 1: Find the product of all the numbers
    product = reduce(*, nums)
    # Step 2: Find the partial product of each number.
    partial_products = map((x) -> Int(product / x), nums)
    # Step 3: Find the modular multiplicative inverse of number[i] modulo partialProduct[i].
    # TODO: should this order be flipped ^
    inverses = []
    for i in 1:length(nums)
        push!(inverses, invmod(partial_products[i], nums[i]))
    end

    # Step 4: Final total
    total = 0
    for i in 1:length(nums)
        total += partial_products[i] * inverses[i] * rems[i];
    end

    return mod(total, product)
end

function run2(l)
    # Parse the input
    nums = []
    rems = []
    for (idx, num) in enumerate(l)
        if num == 0
            # skip 0's, these are the X's in the puzzle input
            continue
        end

        rem = 0
        if idx != 1
            rem = num - (idx - 1)
        end
        push!(nums, num)
        push!(rems, rem)
    end

    result = CRT(nums, rems)
    println("[part 2] Result = ", result)
    return result
end

# NOTE: replaced x's with 0's
@test run2([17,0,13,19]) == 3417.
# 3417 mod 17 = (17-0) # where 0 is offset in list
# 3417 mod 13 = (13-2) # where -2 is offset in list
# 3417 mod 19 = (19-3) # where -3 is offset in list
# i.e.
# (3417 + 0) mod 17 = 0    # where 0 is offset in list
# (3417 + 2) mod 13 = 0    # where -2 is offset in list
# (3417 + 3) mod 19 = 0    # where -3 is offset in list
# i.e.
# trying to solve:
#     (X mod 17) =  (X+2 mod 13) = (X+3 mod 19)

@test run2([67,7,59,61]) == 754018.
@test run2([67,0,7,59,61]) == 779210.
@test run2([67,7,0,59,61]) == 1261476.
@test run2([1789,37,47,1889]) == 1202161486.

run2([17,0,0,0,0,0,0,0,0,0,0,37,0,0,0,0,0,439,0,29,0,0,0,0,0,0,0,0,0,0,13,0,0,0,0,0,0,0,0,0,23,0,0,0,0,0,0,0,787,0,0,0,0,0,0,0,0,0,41,0,0,0,0,0,0,0,0,19])


text:
using Fast12thRoot
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@test fast12throot(0.0) == 0.0
#@test fast12throot(-0.0) == 0.0
@test fast12throot(1.0) == 1.0
@test isnan(fast12throot(NaN)) == true
@test isinf(fast12throot(Inf)) == true

# test smallest and largest normal and subnormal doubles
const vsmallestsubnorm = 2.0^(-1074)
const vlargestsubnorm = 2.0^(-1022)*(1.0-2.0^-52)
const vsmallestnorm = 2.0^(-1023)*(2.0 - 2.0^-52)
const vlargestnorm = 2.0^(1022)*(2.0-2.0^-52)

@test abs(big(vsmallestsubnorm)^(one(BigFloat)/big(12))/fast12throot(vsmallestsubnorm)-1.0)<0.5eps(Float64)
@test abs(big(vlargestsubnorm)^(one(BigFloat)/big(12))/fast12throot(vlargestsubnorm)-1.0)<0.5eps(Float64)
@test abs(big(vsmallestnorm)^(one(BigFloat)/big(12))/fast12throot(vsmallestnorm)-1.0)<0.5eps(Float64)
@test abs(big(vlargestnorm)^(one(BigFloat)/big(12))/fast12throot(vlargestnorm)-1.0)<0.5eps(Float64)

for ii=1:100
    @test fast12throot(Float64(ii)^(12.0)) == Float64(ii)
end


text:
@testset "@test_logs" begin
    function foo(n)
        @info "Doing foo with n=$n"
        for i=1:n
            @debug "Iteration $i"
        end
    end

    @test_logs (Info,"Doing foo with n=2") foo(2)

    # Log pattern matching
    # Regex
    @test_logs (Info,r"^Doing foo with n=[0-9]+$") foo(10)
    @test_logs (Info,r"^Doing foo with n=[0-9]+$") foo(1)
    # Level symbols
    @test_logs (:debug,) min_level=Debug @debug "foo"
    @test_logs (:info,)  @info  "foo"
    @test_logs (:warn,)  @warn  "foo"
    @test_logs (:error,) @error "foo"

    # Pass through so the value of the expression can also be tested
    @test (@test_logs (Info,"blah") (@info "blah"; 42)) == 42

    # Debug level log collection
    @test_logs (Info,"Doing foo with n=2") (Debug,"Iteration 1") (Debug,"Iteration 2") min_level=Debug foo(2)

    @test_logs (Debug,"Iteration 5") min_level=Debug match_mode=:any foo(10)

    # Test failures
    fails = @testset NoThrowTestSet "check that @test_logs detects bad input" begin
        @test_logs (Warn,) foo(1)
        @test_logs (Warn,) match_mode=:any @info "foo"
        @test_logs (Debug,) @debug "foo"
    end
    @test length(fails) == 3
    @test fails[1] isa LogTestFailure
    @test fails[2] isa LogTestFailure
    @test fails[3] isa LogTestFailure
end

if MicroLogging.core_in_base
# Can't meaningfully use @test_deprecated, except with support in Base
function newfunc()
    42
end
@deprecate oldfunc newfunc

@testset "@test_deprecated" begin
    @test_deprecated oldfunc()

    # Expression passthrough
    @test (@test_deprecated oldfunc()) == 42

    fails = @testset NoThrowTestSet "check that @test_deprecated detects bad input" begin
        @test_deprecated newfunc()
        @test_deprecated r"Not found in message" oldfunc()
    end
    @test length(fails) == 2
    @test fails[1] isa LogTestFailure
    @test fails[2] isa LogTestFailure
end
end


text:
using Plots, Colors, ColorTypes # for plotting

export bar, dimplot

################################################################
function invert(c::RGB, distinguish::Bool = false)
     cout = RGB{Float64}(1-c.r,1-c.g,1-c.b)
     if distinguish && cout == c #gray
         cout = RGB{Float64}(0,0,0)
     end
     return cout
end
################################################################
"""
    bar(config::StructArray{Funb}; uniqueType::Bool = false, dimAnnotation::Bool = true, uniqueDim::Bool = false)

Plots bars for each dimension set.

It is assumed that number of dimension sets are the same.

- To have a same color for the same types, set true as the 2nd argument.
- To turn off dimension annotations pass false as the 3rd argument.
- In case of unique dimensions pass true as 4th argument to print dimension in title instead.

# Examples
```julia
bar(configs)
bar(configs, uniqueType = true, dimAnnotation = false, uniqueDim = true)
```
"""
function bar(config::StructArray{Funb}; uniqueType::Bool = false, dimAnnotation::Bool = true, uniqueDim::Bool = false)

    mkpath("bar")

    local titleText

    numFun = length(config.fun)

    bar_width = numFun/3.0*0.2
    bar_text_font = Int64(bar_width*40)
    dim_text_font = Int64(bar_width*30)
    xticks_font = Int64(bar_width*5)

    if uniqueType
        uniqueTypes = uniqueflatten(config.types[:])
        numUniqueTypes =length(uniqueTypes)
        colors = distinguishable_colors(numUniqueTypes, [RGB(1,1,1)])
    end

    _, numDimsSets = numArgsDims(config.dims[1])

    for iDimSet = 1:numDimsSets
        plt = plot()  # different figure for different dims

        for iFun = 1:numFun

            fname = [string(config.fun[iFun])]

            for iType = 1:length(config.types[iFun])

                x = [(iFun-1)+(bar_width*(iType-1))]

                if uniqueType
                    iTypeUnique = findfirst(x -> x==config.types[iFun][iType], uniqueTypes)
                    fillcolor = colors[iTypeUnique]
                    type_text_color = invert(fillcolor, true)'

                    # labels
                    if iFun == 1
                        labels = string(uniqueTypes[iType])
                    else
                        labels = ""
                    end
                    legend = :topleft
                    barText = ""

                else
                    fillcolor = :auto
                    type_text_color = :black

                    barText = Plots.text(string(config.types[iFun][iType]), pointsize = bar_text_font, :center, rotation = 90, color = type_text_color )

                    labels = ""
                    legend = false
                end

                if uniqueDim
                    dimText = ""
                    titleText = "Dimension $(string(config[1].dims[iFun][1,iDimSet]))"
                else
                    if dimAnnotation
                        dimText = Plots.text(stringMatrix(config.dims[iFun][:,iDimSet]), pointsize = dim_text_font, :center)
                    else
                        dimText = ""
                    end
                    titleText = "Dimension set $iDimSet"
                end

                y = [config.median[iFun][iType][iDimSet]]
                # adding bar
                bar!(plt,
                    x,
                    y,
                    labels = labels,
                    legend = legend,
                    fillcolor = fillcolor,
                    bar_width = bar_width,
                    annotations = ([x, x], [y./2, y .+ bar_width/8], [barText, dimText]),
                    dpi = 600
                )
            end
        end

        xticks!(0:numFun-1, string.(config.fun), rotation = 70, fontsize = xticks_font)

        title!(titleText)
        ylabel!("Time [micro seconds]")
        if uniqueType
            filename = "bar/bench-dims-set$iDimSet-unique.png"
        else
            filename = "bar/bench-dims-set$iDimSet.png"
        end
        savefig(filename)
    end
end

bar(config::Funb; uniqueType::Bool = false, dimAnnotation::Bool = true, uniqueDim::Bool = false) = bar(FunbArray(config), uniqueType=uniqueType, dimAnnotation=dimAnnotation, uniqueDim = uniqueDim)

"""
    bar(config::Pair{StructArray{Funb}}; uniqueType::Bool = false, dimAnnotation::Bool = false, uniqueDim::Bool =false, configName::Pair{String,String} = "1" => "2")

Gets a pair of StructArrays and calculates relative speed of coresponding elements and plots them.

Uses the first element of the pair for the configurations and only uses runtimes from the 2nd element of the pair. The relative speed is calculated like this:

1st element runtimes / 2nd element runtimes

Give the configsets names as a pair for the ylabel.
```julia
bar(configs => configs, uniqueType = true, dimAnnotation = false, uniqueDim =true, "group 1" => "group 2")
```
"""
function bar(config::Pair{StructArray{Funb,T1,T2,T3}, StructArray{Funb,T1,T2,T3}}; uniqueType::Bool = false, dimAnnotation::Bool = false, uniqueDim::Bool =false, configName::Pair{String,String} = "1" => "2") where {T1,T2,T3}

    mkpath("bar")
    local titleText

    bar_width = 0.2
    bar_text_font = Int64(bar_width*40)
    dim_text_font = Int64(bar_width*30)
    xticks_font = Int64(bar_width*5)

    if uniqueType
        uniqueTypes = uniqueflatten(config[1].types[:])
        numUniqueTypes =length(uniqueTypes)
        colors = distinguishable_colors(numUniqueTypes, [RGB(1,1,1)])
    end

    numFun = length(config[1].fun)
    _, numDimsSets = numArgsDims(config[1].dims[1])

    for iDimSet = 1:numDimsSets
        plt = plot()  # different figure for different dims

        for iFun = 1:numFun

            fname = [string(config[1].fun[iFun])]

            for iType = 1:length(config[1].types[iFun])

                x = [(iFun-1)+(bar_width*(iType-1))]

                if uniqueType
                    iTypeUnique = findfirst(x -> x==config[1].types[iFun][iType], uniqueTypes)
                    fillcolor = colors[iTypeUnique]
                    type_text_color = invert(fillcolor, true)

                    # labels
                    if iFun == 1
                        labels = string(uniqueTypes[iType])
                    else
                        labels = ""
                    end
                    legend = :topleft
                    barText = ""
                else
                    fillcolor = :auto
                    type_text_color = :black

                    barText = Plots.text(string(config[1].types[iFun][iType]), pointsize = bar_text_font, :center, rotation = 90, color = type_text_color )

                    labels = ""
                    legend = false
                end

                if uniqueDim
                    dimText = ""
                    titleText = "Dimension $(string(config[1].dims[iFun][1,iDimSet]))"
                else
                    if dimAnnotation
                        dimText = Plots.text(stringMatrix(config[1].dims[iFun][:,iDimSet]), pointsize = dim_text_font, :center)
                    else
                        dimText = ""
                    end
                    titleText = "Dimension set $iDimSet"
                end

                y = [config[1].median[iFun][iType][iDimSet] / config[2].median[iFun][iType][iDimSet]]
                # adding bar
                bar!(plt,
                    x,
                    y,
                    labels = labels,
                    legend = legend,
                    fillcolor = fillcolor,
                    bar_width = bar_width,
                    annotations = ([x, x], [y./2, y .+ bar_width/8], [barText, dimText]),
                    dpi = 600
                )
            end
        end

        xticks!(0:numFun-1, string.(config[1].fun), rotation = 70, fontsize = xticks_font)

        title!(titleText)
        ylabel!("Relative Speed ($(configName[2])/$(configName[1]))")
        if uniqueType
            hline!([1], line=(4, :dash, 0.6, [:green]), labels = 1)
        end
        savefig("bar/bench-dims-set$iDimSet-relative.png")
    end
end

################################################################
"""
    dimplot(config::StructArray{Funb})

Shows runtime of a function for different dimension sets. The axes are logarithmic.

# Examples
```julia
configs2 = Funb( sin, [(-1,1)],[Float32, Float64], [10 20 30 40] );

benchmark!(configs2)

dimplot(configs2)
```
"""
function dimplot(config::StructArray{Funb})

mkpath("dimplot")

bar_width = 0.2
bar_text_font = Int64(bar_width*40)
dim_text_font = Int64(bar_width*30)
xticks_font = Int64(bar_width*5)

numFun = length(config.fun)
_, numDimsSets = numArgsDims(config.dims[1])

    for iFun = 1:numFun

        plt = plot()  # different figure for different dims

        fname = [string(config.fun[iFun])]

        for iType = 1:length(config.types[iFun])

            x = 1:numDimsSets
            xticks = string.(dropdims(flatten([config.dims[iFun][1,iDimSet] for iDimSet = 1:numDimsSets]), dims=1))
            y = [config.median[iFun][iType][iDimSet] for iDimSet = 1:numDimsSets]
            # adding bar
            plot!(plt,
                x,
                y,
                labels = string(config.types[iFun][iType]),
                legend = :bottomright,
                xscale = :log10,
                yscale = :log10,
                dpi = 600
            )
            xticks!(x, xticks, fontsize = xticks_font)
        end

        title!("$(fname[1])")
        xlabel!("Dimension")
        ylabel!("Time [micro seconds]")

        savefig("dimplot/bench-$(fname[1]).png")
    end
end

dimplot(config::Funb) = dimplot(FunbArray(config))

"""
    dimplot(config::Vector{StructArray{Funb}}, labels::Vector{String})

By passing a vector different benchmark sets can be grouped together to be shown in a single plot. The axes are logarithmic.

# Examples
```julia
dimplot([config1,config2, config3], ["group 1", "group 2", "group3"])
```
"""
function dimplot(config::Vector{StructArray{Funb,T1,T2,T3}}, labels::Vector{String}) where {T1,T2,T3}

mkpath("dimplot")
numLabels = length(config)

bar_width = 0.2
bar_text_font = Int64(bar_width*40)
dim_text_font = Int64(bar_width*30)
xticks_font = Int64(bar_width*5)

numFun = length(config[1].fun)
_, numDimsSets = numArgsDims(config[1].dims[1])

    for iFun = 1:numFun

        fname = [string(config[1].fun[iFun])]

        for iType = 1:length(config[1].types[iFun])

            plt = plot()  # different figure for different dims

            x = 1:numDimsSets
            xticks = string.(dropdims(flatten([config[1].dims[iFun][1,iDimSet] for iDimSet = 1:numDimsSets]), dims=1))

            for iLabel = 1:numLabels

                y = [config[iLabel].median[iFun][iType][iDimSet] for iDimSet = 1:numDimsSets]

                # adding bar
                plot!(plt,
                    x,
                    y,
                    # labels = string(config[1].types[iFun][iType]),
                    labels = labels[iLabel],
                    legend = :right,
                    xscale = :log10,
                    yscale = :log10,
                    dpi = 600
                )
                xticks!(x, xticks, fontsize = xticks_font)

            end

            title!("$(fname[1]) - Type $(config[1].types[iFun][iType])")
            xlabel!("Dimension")
            ylabel!("Time [micro seconds]")

            savefig("dimplot/bench-$(fname[1])-Type-$(config[1].types[iFun][iType]).png")
        end
    end
end


text:
function visit_internal_node(obj::BoundaryBuilder, arg0::BSPTree)
    return jcall(obj, "visitInternalNode", void, (BSPTree,), arg0)
end

function visit_leaf_node(obj::BoundaryBuilder, arg0::BSPTree)
    return jcall(obj, "visitLeafNode", void, (BSPTree,), arg0)
end

function visit_order(obj::BoundaryBuilder, arg0::BSPTree)
    return jcall(obj, "visitOrder", BSPTreeVisitor_Order, (BSPTree,), arg0)
end