text:
"""
    list_libraries(conn::Union{LibPQ.Connection, DBInterface.Connection})

Load the list of Postgres schemata
the user has permission to access
"""
function list_libraries(conn::Union{LibPQ.Connection, DBInterface.Connection})
    query = """
        WITH RECURSIVE "names"("name") AS (
            SELECT n.nspname AS "name"
                FROM pg_catalog.pg_namespace n
                WHERE n.nspname !~ '^pg_'
                    AND n.nspname <> 'information_schema')
            SELECT "name"
                FROM "names"
                WHERE pg_catalog.has_schema_privilege(
                    current_user, "name", 'USAGE') = TRUE;
        """

        return run_sql_query(conn, query)
end

"""
    check_schema_perms(conn::Union{LibPQ.Connection, DBInterface.Connection}, library::String)

Verify that the user can access a schema
"""
function check_schema_perms(conn::Union{LibPQ.Connection, DBInterface.Connection}, library::String)
    # Verify that the library exists

    if library in list_libraries(conn).name
        return true
    else
        schemas = run_sql_query(conn, query).schema_name
        if library in schemas
            error("You do not have permission to access the $library library")
        else
            error("The $library library is not found")
        end
    end
end

"""
    list_tables(conn::Union{LibPQ.Connection, DBInterface.Connection}, library::String)

List all of the views/tables/foreign tables within a schema
"""
function list_tables(conn::Union{LibPQ.Connection, DBInterface.Connection}, library::String)
    if check_schema_perms(conn, library)
        query = """SELECT table_name FROM INFORMATION_SCHEMA.views
                    WHERE table_schema IN ('$library');"""
        return run_sql_query(conn, query)
    end
end

"""
    approx_row_count(conn::Union{LibPQ.Connection, DBInterface.Connection}, library::String, table::String)

Get an approximate count of the number of rows in a table
"""
function approx_row_count(conn::Union{LibPQ.Connection, DBInterface.Connection}, library::String, table::String)
    if check_schema_perms(conn, library)
        query = """
            SELECT reltuples::bigint AS estimate
            FROM   pg_class
            WHERE  oid = '$library.$table'::regclass;
        """
        #TODO: We should be able to cast the result directly to int, instead of DF
        return run_sql_query(conn, query)[1, "estimate"]
    end

end

"""
    describe_table(conn::Union{LibPQ.Connection, DBInterface.Connection}, library::String, table::String)

Get a table's description (row count, columns, column types)
"""
function describe_table(conn::Union{LibPQ.Connection, DBInterface.Connection}, library::String, table::String)
    if check_schema_perms(conn, library)
        row_count = approx_row_count(conn, library, table)
        println("There are approximately $row_count rows in $library.$table")

        query = """
            SELECT column_name, data_type, is_nullable
            FROM
                information_schema.columns
            WHERE
                table_name = '$table';
        """
        return run_sql_query(conn, query)
    end
end

"""
    get_table(
                    conn::Union{LibPQ.Connection, DBInterface.Connection},
                    library::String,
                    table::String;
                    obs::Union{Nothing, Int} = nothing,
                    offset::Int = 0,
                    cols = nothing
                )

Create a DataFrame from a table
"""
function get_table(
                    conn::Union{LibPQ.Connection, DBInterface.Connection},
                    library::String,
                    table::String;
                    obs::Union{Nothing, Int} = nothing,
                    offset::Int = 0,
                    cols = nothing
                )
    if check_schema_perms(conn, library)

        limit = obs !== nothing && obs > 0 ? "LIMIT $obs" : ""
        columns = cols === nothing ? "*" : join(cols, ",")

        query = "SELECT $columns FROM $library.$table $limit OFFSET $offset"

        return run_sql_query(conn, query)
    end
end

"""
    raw_sql(
        conn::Union{LibPQ.Connection, DBInterface.Connection},
        query::String
    )

Executes raw sql code, and converts code to a DataFrame
"""
function raw_sql(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    query::String
)
    return run_sql_query(conn, query)
end


text:
@testset "HashDictionary" begin
    @test HashDictionary() isa HashDictionary{Any, Any}

    d = HashDictionary{Int64, Int64}()

    @test eltype(d) === Int64
    @test keytype(d) === Int64
    @test tokentype(d) === Int
    @test Base.IteratorSize(d) === Base.HasLength()
    @test isinsertable(d)
    @test length(d) == 0
    @test isempty(d)
    @test isempty(keys(d))
    @test d == d
    @test d == copy(d)
    @test d == HashDictionary(copy(keys(d)), d)
    @test isequal(d, d)
    @test isequal(copy(d), d)
    @test isequal(HashDictionary(copy(keys(d)), d), d)
    @test !isless(d, d)
    @test !isless(copy(d), d)
    @test !isless(HashDictionary(copy(keys(d)), d), d)
    @test cmp(d, d) == 0
    @test cmp(copy(d), d) == 0
    @test cmp(HashDictionary(copy(keys(d)), d), d) == 0
    @test_throws IndexError d[10]
    @test get(d, 10, 15) == 15
    @test length(unset!(d, 10)) == 0
    io = IOBuffer(); print(io, d); @test String(take!(io)) == "{}"
    io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "0-element HashDictionary{Int64,Int64}"
    @test_throws IndexError d[10] = 11
    @test_throws IndexError delete!(d, 10)

    insert!(d, 10, 11)

    @test d[10] == 11
    @test get(d, 10, 15) == 11
    @test get!(d, 10, 15) == 11
    @test length(d) == 1
    @test keys(d) === keys(keys(d))
    @test unique(d)::HashIndices{Int64} == HashIndices([11])
    @test !isempty(d)
    @test d == d
    @test d == copy(d)
    @test d == HashDictionary(copy(keys(d)), d)
    @test d != empty(d)
    @test empty(d) != d
    @test fill(0, d) != d
    @test d != fill(0, d)
    @test fill(0, copy(keys(d))) != d
    @test d != fill(0, copy(keys(d)))
    @test isequal(d, d)
    @test isequal(copy(d), d)
    @test isequal(HashDictionary(copy(keys(d)), d), d)
    @test !isequal(d, empty(d))
    @test !isequal(empty(d), d)
    @test !isequal(fill(0, d), d)
    @test !isequal(d, fill(0, d))
    @test !isequal(fill(0, copy(keys(d))), d)
    @test !isequal(d, fill(0, copy(keys(d))))
    @test !isless(d, d)
    @test !isless(copy(d), d)
    @test !isless(HashDictionary(copy(keys(d)), d), d)
    @test !isless(d, empty(d))
    @test isless(empty(d), d)
    @test isless(fill(0, d), d)
    @test !isless(d, fill(0, d))
    @test isless(fill(0, copy(keys(d))), d)
    @test !isless(d, fill(0, copy(keys(d))))
    @test cmp(d, d) == 0
    @test cmp(copy(d), d) == 0
    @test cmp(HashDictionary(copy(keys(d)), d), d) == 0
    @test cmp(d, empty(d)) == 1
    @test cmp(empty(d), d) == -1
    @test cmp(fill(0, d), d) == -1
    @test cmp(d, fill(0, d)) == 1
    @test cmp(fill(0, copy(keys(d))), d) == -1
    @test cmp(d, fill(0, copy(keys(d)))) == 1
    @test_throws IndexError insert!(d, 10, 12)
    @test d[10] == 11
    set!(d, 10, 12)
    @test length(d) == 1
    @test d[10] == 12
    d[10] = 13
    @test length(d) == 1
    @test d[10] == 13
    io = IOBuffer(); print(io, d); @test String(take!(io)) == "{10 │ 13}"
    io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "1-element HashDictionary{Int64,Int64}\n 10 │ 13"
    @test !isequal(d, empty(d))
    @test isequal(d, copy(d))
    @test isempty(empty(d))

    delete!(d, 10)
    @test isequal(d, HashDictionary{Int64, Int64}())

    @test get!(d, 10, 14) == 14
    @test d[10] == 14
    delete!(d, 10)

    for i = 2:2:1000
        insert!(d, i, i+1)
        @test d[i] == i + 1 ? true : (@show i; false)
    end
    @test all(in(i, d) == !iseven(i) for i in 2:2:1000)
    @test all(in(i, keys(d)) == iseven(i) for i in 2:2:1000)
    @test isempty(empty!(d))

    # TODO token interface
end

ext: jl


text:
using Printf
using Oceananigans
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.OutputWriters
using Oceananigans.Diagnostics
using Oceananigans.Utils
using Oceananigans.AbstractOperations

arch = GPU()
FT   = Float64

# length of the simulation
days = 0.25
end_time = days*day

# output name
filename_1 = "free_convection_day$(days)"
# output directory
base_dir = "les_output/$(filename_1)/"

# NetCDF output writer can't make new directories or overwrite files by itself
try
    mkdir(base_dir) # make output directory
catch SystemError
    # directory exists, so remove .nc files within it
    cd(base_dir)
    run(`rm -f "*".nc`)
    cd("../../")
end

## grid
Lx = Ly = Lz = 100 # extent
Nx = Ny = Nz = 128 # make sure that gridpoints are a multiple of 16 for the GPU

topology = (Periodic, Periodic, Bounded) # horizontally periodic BCs
grid = RegularCartesianGrid(topology=topology, size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=(-Lz, 0))
Hz = grid.Hz # number of halo points

## rotation
f = 1e-4 # coriolis parameter (positive - northern hemisphere)
coriolis = FPlane(FT, f=f)

# α = 2e-4  # Thermal expansion coefficient [K⁻¹]
# eos = LinearEquationOfState(FT, α=α, β=0)
buoyancy = BuoyancyTracer()

κh = νh = 1e-5   # Horizontal diffusivity and viscosity [m²/s]
κv = νv = 1e-5   # Vertical diffusivity and viscosity [m²/s]

N_s = 1e-2 * 2e-3  # Uniform vertical stratification [s⁻¹]
buoyancy_flux = 1.96e-7 / 4 # bouyancy flux

## closure
# closure = ConstantAnisotropicDiffusivity(FT, νh=νh, νv=νv, κh=κh, κv=κv)
closure = AnisotropicMinimumDissipation(FT)

#  *--*--*--*--*--*--*--*--*
#  |  Boundary Conditions  |
#  *--*--*--*--*--*--*--*--*

top_b_bc = FluxBoundaryCondition(buoyancy_flux)
bottom_b_bc = GradientBoundaryCondition(N_s)
b_bcs = TracerBoundaryConditions(grid, top = top_b_bc, bottom = bottom_b_bc)
top_C_bc = ValueBoundaryCondition(1.0)
C_bcs = TracerBoundaryConditions(grid, top=top_C_bc)
bcs = (b = b_bcs,)

#  *--*--*--*
#  | Model  |
#  *--*--*--*

model = IncompressibleModel(
            architecture = arch,
                float_type = FT,
                    grid = grid,
                coriolis = coriolis,
                buoyancy = buoyancy,
                closure = closure,
                tracers = (:b,),
    boundary_conditions = bcs
)

#  *--*--*--*--*--*--*--*--*
#  |  Initial Conditions   |
#  *--*--*--*--*--*--*--*--*

T_s  = 12.0    # Surface temperature [°C]
ε(σ) = σ * randn()
B₀(x, y, z) = N_s * (z + 100) + ε(1e-8) + 20*2e-3
set!(model, b=B₀) # set ICs

#  *--*--*--*--*
#  | Timestep  |
#  *--*--*--*--*

# Δt= 1.0
Δt_wizard = TimeStepWizard(cfl=0.3, Δt = 1.0, max_change=1.0, max_Δt = 3.0)
cfl = AdvectiveCFL(Δt_wizard)

#  *--*--*--*--*--*
#  | Simulation   |
#  *--*--*--*--*--*

# Take Ni "intermediate" time steps at a time before printing a progress
# statement and updating the time step.
Ni = 50

function print_progress(simulation)
    model = simulation.model
    i, t = model.clock.iteration, model.clock.time

    progress = 100 * (model.clock.time / end_time)

    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)

    @printf("[%05.2f%%] i: %d, t: %.2e days, umax: (%6.3e, %6.3e, %6.3e) m/s, CFL: %6.4e, next Δt: %.1e s\n",
    	    progress, i, t / day, umax, vmax, wmax, cfl(model), Δt_wizard.Δt)
end

simulation = Simulation(model, Δt=Δt_wizard, stop_time=end_time, progress=print_progress, progress_frequency=Ni)

u = model.velocities.u
v = model.velocities.v
w = model.velocities.w
b = model.tracers.b

#  *--*--*--*--*--*--*
#  | Output Writers  |
#  *--*--*--*--*--*--*
#  Set up output writers to write output variables to NetCDF (.nc) files

global_attributes = Dict(
    "f"                 => Dict("longname" => "Coriolis parameter", "units" => "---"),
    "model.closure.ν"   => Dict("longname" => "Coriolis parameter", "units" => "---"),
    "model.closure.ν"   => Dict("longname" => "Viscosity", "units" => "---"),
    "model.closure.κ.b" => Dict("longname" => "diffusivity", "units" => "---"),
    "Q_str"             => Dict("longname" => "Net surface heat flux into the ocean (+=down), >0 increases theta", "units" => "W/m²"),
    "f"                 => Dict("longname" => "Coriolis parameter", "units" => "---"),
    "rho"               => Dict("longname" => "Density", "units" => "---"),
    "ρ₀"                => Dict("longname" => "Density", "units" => "---")
    # file["parameters/viscosity"] = model.closure.ν
    # file["parameters/diffusivity"] = model.closure.κ.b
    # file["parameters/diffusivity_T"] = model.closure.κ.b
    # file["parameters/diffusivity_S"] = model.closure.κ.S
    # file["parameters/surface_cooling"] = Q_str
    # file["parameters/temperature_gradient"] = ∂θ∂z
    # file["boundary_conditions/top/FT"] = Fθ_str
    # file["boundary_conditions/top/Fu"] = Fu_str
    # file["boundary_conditions/bottom/dTdz"] = ∂θ∂z
)

#every output (each key in profiles) must have an associated attribute
output_attributes = Dict(
     "u"  => Dict("longname" => "Velocity in the x-direction", "units" => "m/s"),
     "v"  => Dict("longname" => "Velocity in the y-direction", "units" => "m/s"),
     "w"  => Dict("longname" => "Velocity in the z-direction", "units" => "m/s"),
     "b" => Dict("longname" => "Buoyancy",                    "units" => "m/s²"),
     "nu" => Dict("longname" => "Eddy viscosity",              "units" => "m²/s"), #####
 "kappaT" => Dict("longname" => "Eddy diffusivity of conservative temperature", "units" => "m²/s"),
     "uu" => Dict("longname" => "Velocity covariance between u and u", "units" => "m²/s²"),
     "vv" => Dict("longname" => "Velocity covariance between v and v", "units" => "m²/s²"),
     "ww" => Dict("longname" => "Velocity covariance between w and w", "units" => "m²/s²"),
     "uv" => Dict("longname" => "Velocity covariance between u and v", "units" => "m²/s²"),
     "uw" => Dict("longname" => "Velocity covariance between u and w", "units" => "m²/s²"),
     "vw" => Dict("longname" => "Velocity covariance between v and w", "units" => "m²/s²"),
     "wT" => Dict("longname" => "Vertical turbulent heat flux", "units" => "K*m/s"),
 "vshear" => Dict("longname" => "Conservative temperature",    "units" => "°C")
)

#  *--*--*--*--*--*--*--*--*--*--*--*
#  |  NetCDF Output Writer: Fields  |
#  *--*--*--*--*--*--*--*--*--*--*--*

fields = Dict(
    "u" => u,
    "v" => v,
    "w" => w,
    "b" => b
)

# field_writer = NetCDFOutputWriter(model, fields;
#                                   filename=base_dir*filename_1*"_fields.nc",
#                                 # global_attributes=global_attributes,
#                                 # output_attributes=output_attributes, # use default output_attributes
#                                 # dimensions = field_dims, # computes field variable dimensions automatically
#                                   interval=2hour)

# simulation.output_writers[:field_writer] = field_writer

#  *--*--*--*--*--*--*--*--*--*--*--*
#  | NetCDF Output Writer: Profiles |
#  *--*--*--*--*--*--*--*--*--*--*--*

# diagnostics
 Up    = HorizontalAverage(u;                        return_type=Array)
 Vp    = HorizontalAverage(v;                        return_type=Array)
 Wp    = HorizontalAverage(w;                        return_type=Array)
 Tp    = HorizontalAverage(b;                        return_type=Array)
 νp    = HorizontalAverage(model.diffusivities.νₑ;   return_type=Array)
κTp    = HorizontalAverage(model.diffusivities.κₑ.b; return_type=Array)
uu     = HorizontalAverage(u*u, model;               return_type=Array)
vv     = HorizontalAverage(v*v, model;               return_type=Array)
ww     = HorizontalAverage(w*w, model;               return_type=Array)
uv     = HorizontalAverage(u*v, model;               return_type=Array)
uw     = HorizontalAverage(u*w, model;               return_type=Array)
vw     = HorizontalAverage(v*w, model;               return_type=Array)
wT     = HorizontalAverage(w*b, model;               return_type=Array)
vshear = HorizontalAverage(∂z(u)^2 + ∂z(v)^2, model; return_type=Array)

# output variables
profiles = Dict(
     "u"  => model ->  Up(model)[1,1,1+Hz:end-Hz],
     "v"  => model ->  Vp(model)[1,1,1+Hz:end-Hz],
     "w"  => model ->  Wp(model)[1,1,1+Hz:end-Hz],
     "b"  => model ->  Tp(model)[1,1,1+Hz:end-Hz],
     "nu" => model ->  νp(model)[1,1,1+Hz:end-Hz],
 "kappaT" => model -> κTp(model)[1,1,1+Hz:end-Hz],
     "uu" => model ->  uu(model)[1,1,1+Hz:end-Hz],
     "vv" => model ->  vv(model)[1,1,1+Hz:end-Hz],
     "ww" => model ->  ww(model)[1,1,1+Hz:end-Hz],
     "uv" => model ->  uv(model)[1,1,1+Hz:end-Hz],
     "uw" => model ->  uw(model)[1,1,1+Hz:end-Hz],
     "vw" => model ->  vw(model)[1,1,1+Hz:end-Hz],
     "wT" => model ->  wT(model)[1,1,1+Hz:end-Hz],
 "vshear" => model -> vshear(model)[1,1,1+Hz:end-Hz]
)

profile_dims = Dict(k => ("zC",) for k in keys(profiles))
profile_dims["w"] = ("zF",)

profile_writer = NetCDFOutputWriter(model, profiles;
                                    filename = base_dir*filename_1*"_profiles.nc",
                                    # global_attributes = global_attributes,
                                    output_attributes = output_attributes,
                                    dimensions = profile_dims,
                                    interval = 1hour)

simulation.output_writers[:profile_writer] = profile_writer

#  *--*--*--*--*--*--*
#  | Run simulation  |
#  *--*--*--*--*--*--*

# time_step!(model, 1)
run!(simulation)

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
######################## CVRPLIB parsers ########################
function parse_cvrplib(file_path::String)
    scan = Scanner.Scan(file_path)
    garbage = Scanner.next(scan, String)
    garbage = Scanner.next(scan, String)
    id = Scanner.next(scan, String)
    garbage = Scanner.nextline(scan)
    garbage = Scanner.nextline(scan)
    n = Scanner.next(scan, Int)
    garbage = Scanner.nextline(scan)
    capacity = Scanner.next(scan, Int)
    garbage = Scanner.nextline(scan)
    points = Scanner.nextmatrix(scan, Int, n, 3, rowmajor = true)[1:end,2:end]
    garbage = Scanner.nextline(scan)
    demands = Scanner.nextmatrix(scan, Int, n, 2, rowmajor = true)[1:end,2:end]
    garbage = Scanner.nextline(scan)
    depot_idx = Scanner.next(scan, Int)
    Scanner.finish_scan(scan)

    xs = [points[i,1] for i in 1:n]
    ys = [points[i,2] for i in 1:n]

    travel_matrix_periods = [Range()]
    period_to_matrix_id = Dict{Range,String}(travel_matrix_periods[1] =>
                                             "unique_mat")
    mat = generate_symmetric_distance_matrix(xs, ys)
    for j in 1:n
        for i in 1:n
            mat[i,j] = round(mat[i,j], RoundingMode{:Nearest}())
        end
    end
    travel_distance_matrices = Dict{String,Array{Float64,2}}("unique_mat" => mat)
    travel_time_matrices = Dict{String,Array{Float64,2}}()
    energy_consumption_matrices = Dict{String,Array{Float64,2}}()

    locations = [Location(
        id = string("loc_", i), index = i,
        long_x = xs[i], lat_y = ys[i]
    ) for i in 1:n]
    locations[depot_idx].id = "depot"
    location_groups = create_singleton_location_groups(locations)

    product_compatibility_classes = ProductCompatibilityClass[]
    product_sharing_classes = ProductSharingClass[]
    product_specification_classes = [ProductSpecificationClass(
        id = "unique_p_spec_c",
        capacity_consumptions =
                Dict{String,Tuple{Float64,Float64}}("unique_measure" =>
                                                    (1.0,1.0))
    )]

    requests = Request[]
    req_idx = 1
    for i in 1:n
        if i != depot_idx
            req = Request(
                request_type = 1,
                id = string("req_", req_idx),
                product_specification_class_id = "unique_p_spec_c",
                product_quantity_range = Range(demands[i]),
                pickup_location_group_id = location_groups[i].id
            )
            push!(requests, req)
            req_idx += 1
        end
    end

    vehicle_categories = [VehicleCategory(
        id = "unique_vehicle_category",
        capacities = VehicleCharacteristics(of_vehicle = Dict{String,Float64}("unique_measure" => capacity))
    )]
    vehicle_cost = CostPeriod(travel_distance_unit_cost = 1.0)
    vehicle_sets = [HomogeneousVehicleSet(
        id = "unique_vehicle_set",
        vehicle_category_id = "unique_vehicle_category",
        cost_periods = [vehicle_cost],
        departure_location_group_id = "depot_loc_group",
        arrival_location_group_id = "depot_loc_group",
        nb_of_vehicles_range = FlexibleRange(soft_range = Range(0, n-1),
                                             hard_range = Range(0, n-1))
    )]

    travel_specifications = [TravelSpecification("unique_period_cat",
            Array{Float64,2}(undef, 0, 0), mat, Array{Float64,2}(undef, 0, 0))]
    travel_periods = [TravelPeriod(Range(),
                                    "unique_period_cat")]

    coordinate_mode = 0
    distance_mode = 0
    data = RvrpInstance(
        id, coordinate_mode, distance_mode, travel_specifications,
        travel_periods, locations, location_groups,
        product_compatibility_classes, product_sharing_classes,
        product_specification_classes, requests, vehicle_categories,
        vehicle_sets
    )
    preprocess_instance(data)
    return data

end


text:
# Error type and decoding functionality

import Base: ==, show, showerror

export CuError

#
# API requirements errors
#

immutable CuVersionError <: Exception
    symbol::Symbol
    minver::VersionNumber
end

function showerror(io::IO, err::CuVersionError)
    @printf(io, "CuVersionError: call to %s requires at least driver v%s",
            err.symbol, err.minver)
end

#
# Runtime API errors
#

immutable CuError <: Exception
    code::Int
    info::Nullable{String}

    CuError(code) = new(code, Nullable{String}())
    CuError(code, info) = new(code, Nullable{String}(info))
end

==(x::CuError,y::CuError) = x.code == y.code

const return_codes = Dict{Int,Tuple{Symbol,String}}(
    0   => (:SUCCESS, "Success"),

    1   => (:ERROR_INVALID_VALUE, "Invalid value"),
    2   => (:ERROR_OUT_OF_MEMORY, "Out of memory"),
    3   => (:ERROR_NOT_INITIALIZED, "Driver not initialized"),
    4   => (:ERROR_DEINITIALIZED, "Driver being shutdown"),
    5   => (:ERROR_PROFILER_DISABLED, "Profiler disabled"),
    6   => (:ERROR_PROFILER_NOT_INITIALIZED, "Profiler not initialized"),
    7   => (:ERROR_PROFILER_ALREADY_STARTED, "Profiler already started"),
    8   => (:ERROR_PROFILER_ALREADY_STOPPED, "Profiler already stopped"),
    100 => (:ERROR_NO_DEVICE, "No CUDA-capable device"),
    101 => (:ERROR_INVALID_DEVICE, "Invalid device ordinal"),
    200 => (:ERROR_INVALID_IMAGE, "Invalid kernel image"),
    201 => (:ERROR_INVALID_CONTEXT, "Invalid context"),
    202 => (:ERROR_CONTEXT_ALREADY_CURRENT, "Context already current"),
    205 => (:ERROR_MAP_FAILED, "Map operation failed"),
    206 => (:ERROR_UNMAP_FAILED, "Unmap operation failed"),
    207 => (:ERROR_ARRAY_IS_MAPPED, "Array mapped"),
    208 => (:ERROR_ALREADY_MAPPED, "Resource already mapped"),
    209 => (:ERROR_NO_BINARY_FOR_GPU, "No kernel image available/suitable for GPU"),
    210 => (:ERROR_ALREADY_ACQUIRED, "Resource already acquired"),
    211 => (:ERROR_NOT_MAPPED, "Resource not mapped"),
    212 => (:ERROR_NOT_MAPPED_AS_ARRAY, "Resource not mapped as array"),
    213 => (:ERROR_NOT_MAPPED_AS_POINTER, "Resource not mapped as pointer"),
    214 => (:ERROR_ECC_UNCORRECTABLE, "Uncorrectable ECC error detected"),
    215 => (:ERROR_UNSUPPORTED_LIMIT, "Unsupported limit"),
    216 => (:ERROR_CONTEXT_ALREADY_IN_USE, "Context already in use"),
    217 => (:ERROR_PEER_ACCESS_UNSUPPORTED, "Peer access not supported"),
    218 => (:ERROR_INVALID_PTX, "Invalid PTX code"),
    300 => (:ERROR_INVALID_SOURCE, "Invalid kernel source"),
    301 => (:ERROR_FILE_NOT_FOUND, "File not found"),
    302 => (:ERROR_SHARED_OBJECT_SYMBOL_NOT_FOUND, "Shared object symbol not found"),
    303 => (:ERROR_SHARED_OBJECT_INIT_FAILED, "Shared object initialization failed"),
    304 => (:ERROR_OPERATING_SYSTEM, "OS call failed"),
    400 => (:ERROR_INVALID_HANDLE, "Invalid handle"),
    500 => (:ERROR_NOT_FOUND, "Named symbol not found"),
    600 => (:ERROR_NOT_READY, "Not ready"),
    700 => (:ERROR_ILLEGAL_ADDRESS, "Illegal memory access"),
    701 => (:ERROR_LAUNCH_OUT_OF_RESOURCES, "Launch out of resources"),
    702 => (:ERROR_LAUNCH_TIMEOUT, "Launch timeout"),
    703 => (:ERROR_LAUNCH_INCOMPATIBLE_TEXTURING, "Incompatible texturing mode"),
    704 => (:ERROR_PEER_ACCESS_ALREADY_ENABLED, "Peer access already enabled"),
    705 => (:ERROR_PEER_ACCESS_NOT_ENABLED, "Peer access not enabled"),
    708 => (:ERROR_PRIMARY_CONTEXT_ACTIVE, "Primary context already active"),
    709 => (:ERROR_CONTEXT_IS_DESTROYED, "Context destroyed"),
    710 => (:ERROR_ASSERT, "Assertion triggered failure"),
    711 => (:ERROR_TOO_MANY_PEERS, "Too many peers"),
    712 => (:ERROR_HOST_MEMORY_ALREADY_REGISTERED, "Host memory already registered"),
    713 => (:ERROR_HOST_MEMORY_NOT_REGISTERED, "Host memory not registered"),
    714 => (:ERROR_HARDWARE_STACK_ERROR, "Hardware stack error"),
    715 => (:ERROR_ILLEGAL_INSTRUCTION, "Illegal instruction"),
    716 => (:ERROR_MISALIGNED_ADDRESS, "Misaligned address"),
    717 => (:ERROR_INVALID_ADDRESS_SPACE, "Invalid address space"),
    718 => (:ERROR_INVALID_PC, "Invalid program counter"),
    719 => (:ERROR_LAUNCH_FAILED, "Kernel launch failed"),
    800 => (:ERROR_NOT_PERMITTED, "Operation not permitted"),
    801 => (:ERROR_NOT_SUPPORTED, "Operation not supported"),
    999 => (:ERROR_UNKNOWN, "Unknown error")
)

name(err::CuError)        = return_codes[err.code][1]
description(err::CuError) = return_codes[err.code][2]

function showerror(io::IO, err::CuError)
    if isnull(err.info)
        @printf(io, "%s (CUDA error #%d, %s)",
                    description(err), err.code, name(err))
    else
        @printf(io, "%s (CUDA error #%d, %s)\n%s",
                    description(err), err.code, name(err), get(err.info))
    end
end

show(io::IO, err::CuError) =
    @printf(io, "%s(%d)", name(err), err.code)

for code in return_codes
    @eval $(code[2][1]) = CuError($(code[1]))
end


text:
@testset "Interval matrix methods" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.8..4.2 0.0..0.9])
    @test opnorm(m) == opnorm(m, Inf) ≈ 5.2
    @test opnorm(m, 1) ≈ 5.3

    l = inf(m)
    r = sup(m)
    c = mid(m)
    d = diam(m)
    rad = radius(m)
    m2 = copy(m)
    @test m2 isa IntervalMatrix && m.mat == m2.mat
    @test l == inf.(m) && r == sup.(m) && c == mid.(m)
    @test d ≈ r - l
    @test rad ≈ d/2

    sm = scale(m, 2.0)
    @test sm ==  2.0 .* m
    @test sm ≠ m
    scale!(m, 2.0) # in-place
    @test sm == m

    m3 = IntervalMatrix([-2.0..2.0 -2.0..0.0; 0.0..2.0 -1.0..1.0])
    m4 = IntervalMatrix([-1.0..1.0 -1.0..1.0; -1.0..1.0 -2.0..2.0])
    @test m3 ∩ m4 == IntervalMatrix([-1.0..1.0 -1.0..0.0; 0.0..1.0 -1.0..1.0])
    @test m3 ∪ m4 == IntervalMatrix([-2.0..2.0 -2.0..1.0; -1.0..2.0 -2.0..2.0])
    @test diam_norm(m3) ≈ 6.0 # default diameter p-norm is Inf
    @test diam_norm(m3, 1) ≈ 6.0
end

@testset "Interval matrix membership" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    a1 = [0. -4.;  4. 0.]
    a2 = [0. -3.; 4. 0.]
    @test a1 ∈ m
    @test a2 ∉ m
end

@testset "Interval matrix inclusion" begin
    m1 = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    m2 = IntervalMatrix([-1.2..1.0 -4.1.. -3.9; 3.9..4.2 -1.2..0.9])
    @test m1 ⊆ m2 && !(m2 ⊆ m1)

end

@testset "Interval matrix rand" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    for i in 1:3
        @test rand(m) ∈ m
    end
end


text:
using JuMP, Ipopt, LinearAlgebra

## Example from Nocedal 16.2 (p 452-453)
x̄ = [2;-1;1]
λ = [3;-2]

# set up model
m = JuMP.Model(solver=IpoptSolver(print_level=0))

# set up primal variables
n = 3
@variable(m, x[1:n])

# set up constraints
G = [6. 2. 1.; 2. 5. 2.; 1. 2. 4.]
c = [-8.; -3.; -3]
A = [1. 0. 1.; 0. 1. 1.]
b = [3.; 0.]

@objective(m, Min, 0.5*x'*G*x + x'*c)

@constraint(m, con, A*x .== b)

print(m)

status = JuMP.solve(m)

# Solution
println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("λ = ", getdual(con))

##

# set up model
m = JuMP.Model(solver=IpoptSolver())

# set up primal variables
n = 3
@variable(m, x[1:n])

# set up constraints
G = [6. 2. 1.; 2. 5. 2.; 1. 2. 4.]
c = [-8.; -3.; -3]
A = [1. 0. 1.; 0. 1. 1.]
b = [3.; 0.]

@objective(m, Min, x'*x)

@constraint(m, con, G*x - c .== 0)
# @constraint(m, con2, A*x .== b)

print(m)

status = JuMP.solve(m)

# Solution
println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("λ = ", getdual(con))


text:
# find the longest edge of a cell
function longest_edge(grid, cellid)

    max_length = 0.0
    faceid = 0

    cell = getcells(grid, cellid)
    for f in 1:nfaces(cell)
        # assumes linear elements
        node_idx1, node_idx2 = Ferrite.faces(cell)[f]
        l = norm(grid.nodes[node_idx1].x - grid.nodes[node_idx2].x)
        if l > max_length
            max_length = l
            faceid = f
        end
    end

    return FaceIndex(cellid, faceid)
end

# create a new node on the mid-point the the face given by faceindex
function new_node(grid, faceindex)
    cellid, faceid = faceindex
    cell = getcells(grid, cellid)
    facenodes = Ferrite.faces(cell)[faceid]
    X = mean([grid.nodes[i].x for i in facenodes]) # hard-coded linear elements
    return Node(X)
end

# add nodes to new_nodes Vector and to split_edges dict
function add_nodes!(new_nodes, split_edges, grid, face::FaceIndex, nodeid)
    faces = Ferrite.full_neighborhood(grid, face, true)
    if !haskey(split_edges, first(faces)) # this face hasn't been assigned a new node yet
        for f in faces[2:end]
            @assert !haskey(split_edges, f) # the neighboring face should not be stored either
        end
        push!(new_nodes, new_node(grid, first(faces)))
        nodeid += 1
        merge!(split_edges, Dict(faces .=> nodeid))
    end
    return nodeid
end

"""
    refine(grid::Grid, cells_to_split)

Refine the given `grid` by applying the [Rivara algorithm](https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620200412) where `cells_to_split` is the list of cells to define (should contain cell numbers).
The `grid` must include a `Topology` for the refinement to work. Construct the first grid as follows, all grids that are returned by `refine` include the topology information:

```julia
grid = generate_grid(Triangle, (nel_x, nel_y), lower_left, upper_right; ; build_topology=true)
```

This implementation is restricted to linear 2D Triangles.

In order to use this function, you need to use the CA4 branch of Ferrite.jl. You can add it as follows:
```
pkg> add Ferrite#CA4
```

!!! warning
    `nodesets`, `facesets` and `cellsets` are currently lost after refinement and must be reconstructed. They can be added to a grid in the following manner:
    ```julia
    addfaceset!(mesh, "left", x -> x[1] ≈ 0.0)
    addfaceset!(mesh, "bottom", x -> x[2] ≈ 0.0)
    addfaceset!(mesh, "Γₗ", x -> (x[1] <= h && x[2] ≈ h))
    ```
"""
function refine(grid::Grid{2,Triangle}, cells_to_split)

    # compute new nodes
    nnodes = getnnodes(grid)
    nodeid = nnodes
    new_nodes = Node{2,Float64}[]
    split_edges = Dict{FaceIndex, Int}()

    old_cell_list = Int[]
    cell_list = cells_to_split
    while !(old_cell_list == cell_list)
        for cellid in cell_list
            face = longest_edge(grid, cellid)
            nodeid = add_nodes!(new_nodes, split_edges, grid, face, nodeid)
        end
        old_cell_list = copy(cell_list)
        cell_list = unique!([f[1] for f in keys(split_edges)])
    end

    new_cells = Triangle[]
    # construct new cells
    for cellid in cell_list
        longest_face = longest_edge(grid, cellid)
        # n_face = Ferrite.full_neighborhood(grid, longest_face)[1] # assumes only one neighboring face
        node_longest_face = split_edges[longest_face]

        cell = getcells(grid, cellid)
        for f in 1:nfaces(cell)
            if f == longest_face[2]
                continue
            end
            face = Ferrite.faces(cell)[f]
            node_f = get(split_edges, FaceIndex(cellid, f), nothing)
            if node_f === nothing # no node on the face
                c = Triangle((face[1], face[2], node_longest_face))
                push!(new_cells, c)
            else
                c1 = Triangle((face[1], node_f, node_longest_face))
                c2 = Triangle((node_f, face[2], node_longest_face))
                push!(new_cells, c1)
                push!(new_cells, c2)
            end
        end
    end

    # construct new grid
    cells = vcat(deleteat!(grid.cells, sort!(cell_list)), new_cells)
    nodes = vcat(grid.nodes, new_nodes)

    new_grid = Grid(cells, nodes; topology=Ferrite.ExclusiveTopology(cells)) # currently looses facesets etc.
    return new_grid
end

"""
    linear_to_quadratic(dh_lin::DofHandler{2,Triangle}, dh_quad::DofHandler{2,Triangle}, a_lin)

Interpolate the linear solution `a_lin` (associated with `dh_lin`) to the degrees of freedom of the quadratic `DofHandler` `dh_quad`.
Return the interpolated values as a vector ordered according to the dofs in `dh_quad`. Restricted to vector valued fields.
"""
function linear_to_quadratic(dh_lin::DofHandler{2,Triangle}, dh_quad::DofHandler{2,Triangle}, a_lin)
    dim = 2

    qr = QuadratureRule{dim,RefTetrahedron,Float64}([NaN for i=1:6], Vec.([(1.0, 0.0), (0.0, 1.0), (0.0, 0.0), (0.5, 0.5), (0.0, 0.5), (0.5, 0.0)]))
    ip = Lagrange{dim, RefTetrahedron, 1}()
    cv_nodes = CellVectorValues(qr, ip)

    a_quad = Vector{Float64}(undef, ndofs(dh_quad))

    for cell in CellIterator(dh_quad)
        dofs = celldofs(cell)
        dofs_lin = celldofs(dh_lin, cellid(cell))
        for (node, dof) in enumerate(1:dim:length(dofs))
            dof1 = dofs[dof]
            a_quad[dof1:(dof1+dim-1)] = function_value(cv_nodes, node, a_lin[dofs_lin])
        end
    end

    return a_quad
end


text:
module TestHash
using Compat
using Compat.Test
using CategoricalArrays

@testset "hash() for CategoricalPool{Int} and CategoricalPool{Float64} and its values" begin
    pool1 = CategoricalPool([1, 2, 3])
    pool2 = CategoricalPool([2.0, 1.0, 3.0])

    @test (hash(pool1) == hash(pool1)) === true
    @test (hash(pool1) == hash(pool2)) === false
    @test (hash(pool2) == hash(pool2)) === true

    opool1 = CategoricalPool([1, 2, 3], true)
    opool2 = CategoricalPool([2.0, 1.0, 3.0], true)

    @test (hash(opool1) == hash(opool1)) === true
    @test (hash(opool1) == hash(opool2)) === false
    @test (hash(opool2) == hash(opool2)) === true

    nv1a = CategoricalArrays.catvalue(1, pool1)
    nv2a = CategoricalArrays.catvalue(1, pool2)
    nv1b = CategoricalArrays.catvalue(2, pool1)
    nv2b = CategoricalArrays.catvalue(2, pool2)

    @test (hash(nv1a) == hash(nv1a)) === true
    @test (hash(nv1a) == hash(nv2a)) === false
    @test (hash(nv1a) == hash(nv1b)) === false
    @test (hash(nv1a) == hash(nv2b)) === true

    @test (hash(nv1b) == hash(nv1a)) === false
    @test (hash(nv1b) == hash(nv2a)) === true
    @test (hash(nv1b) == hash(nv1b)) === true
    @test (hash(nv1b) == hash(nv2b)) === false

    @test (hash(nv2a) == hash(nv1a)) === false
    @test (hash(nv2a) == hash(nv2a)) === true
    @test (hash(nv2a) == hash(nv1b)) === true
    @test (hash(nv2a) == hash(nv2b)) === false

    @test (hash(nv2b) == hash(nv1a)) === true
    @test (hash(nv2b) == hash(nv2a)) === false
    @test (hash(nv2b) == hash(nv1b)) === false
    @test (hash(nv2b) == hash(nv2b)) === true

    ov1a = CategoricalArrays.catvalue(1, opool1)
    ov2a = CategoricalArrays.catvalue(1, opool2)
    ov1b = CategoricalArrays.catvalue(2, opool1)
    ov2b = CategoricalArrays.catvalue(2, opool2)

    @test (hash(ov1a) == hash(ov1a)) === true
    @test (hash(ov1a) == hash(ov2a)) === false
    @test (hash(ov1a) == hash(ov1b)) === false
    @test (hash(ov1a) == hash(ov2b)) === true

    @test (hash(ov1b) == hash(ov1a)) === false
    @test (hash(ov1b) == hash(ov2a)) === true
    @test (hash(ov1b) == hash(ov1b)) === true
    @test (hash(ov1b) == hash(ov2b)) === false

    @test (hash(ov2a) == hash(ov1a)) === false
    @test (hash(ov2a) == hash(ov2a)) === true
    @test (hash(ov2a) == hash(ov1b)) === true
    @test (hash(ov2a) == hash(ov2b)) === false

    @test (hash(ov2b) == hash(ov1a)) === true
    @test (hash(ov2b) == hash(ov2a)) === false
    @test (hash(ov2b) == hash(ov1b)) === false
    @test (hash(ov2b) == hash(ov2b)) === true

    @testset "ordered and non-ordered values hash equal" begin
    @test (hash(ov1a) == hash(nv1a)) === true
    @test (hash(ov1a) == hash(nv2a)) === false
    @test (hash(ov1a) == hash(nv1b)) === false
    @test (hash(ov1a) == hash(nv2b)) === true

    @test (hash(ov1b) == hash(nv1a)) === false
    @test (hash(ov1b) == hash(nv2a)) === true
    @test (hash(ov1b) == hash(nv1b)) === true
    @test (hash(ov1b) == hash(nv2b)) === false

    @test (hash(ov2a) == hash(nv1a)) === false
    @test (hash(ov2a) == hash(nv2a)) === true
    @test (hash(ov2a) == hash(nv1b)) === true
    @test (hash(ov2a) == hash(nv2b)) === false

    @test (hash(ov2b) == hash(nv1a)) === true
    @test (hash(ov2b) == hash(nv2a)) === false
    @test (hash(ov2b) == hash(nv1b)) === false
    @test (hash(ov2b) == hash(nv2b)) === true
    end
end

end


text:
function multitestH0(projMatrix::Array{Float64,2},dataVar::Array{Float64,1},thrs::Nothing)
    projMatrix = transpose(projMatrix)./dataVar
    (nc,nr) = size(projMatrix)
    kMax = floor(Integer,log2(nr))-1
    inds = round.(Int,exp2.(0:kMax))
    pval = zeros(nc)
    projMatrix = abs2.(projMatrix)
    cs = Array(cumsum(projMatrix,dims=2))[:,(inds[2:end].-1)]
    pval = maximum(cdf.(Chisq.(inds[1:(end-1)]),transpose(hcat(cs[:,inds[2]-1],diff(cs,2)))),dims=1)

   # for (inds,inds2) in zip(inds[1:end-1],inds[2:end] .- 1)

    #    norm2proj = vec(sum(projMatrix[:,inds:inds2],2))

    #    @inbounds   pval = max.(pval,cdf.(Chisq(inds),norm2proj))
    #end
    return((1-pval)[:])
end

function multitestH0(projMatrix::Array{Float64},dataVar,thrs::Array{Float64,1})
    projMatrix = transpose(projMatrix)./dataVar
    (nc,nr) = size(projMatrix)
    kMax = floor(Integer,log2(nr))-1
    inds = round.(Int,exp2.(0:kMax))
    #@show inds
   # test = trues(nc)
    projMatrix = abs2.(projMatrix)
    cs = Array(cumsum(projMatrix,dims=2))[:,inds[2:end].-1]
    norm2proj = permutedims(hcat(cs[:,inds[2]-1],diff(cs,dims=2)))
    #@show size(cs)
    #norm2proj = diff(cs,2)
    #@show size(norm2proj)
    #@show length(thrs)
    test = all(norm2proj.<=thrs,dims = 1)

    #for (inds,inds2,i) in zip(inds[1:end-1],inds[2:end] .- 1,collect(1:length(inds) .- 1))

    #    norm2proj = vec(sum(projMatrix[:,inds:inds2],2))

    #    @inbounds test = test .& (norm2proj.<=thrs[i])
    #end
    return(test[:])
end


text:
push!(LOAD_PATH, "../models")
push!(LOAD_PATH, "../")
using SparseInverseProblems
using SuperResModels
using TestCases
using Utils

using PyCall
@pyimport numpy as np
@pyimport matplotlib.pyplot as plt

# Containers to handle parameters
type StaticParameters
    # allowed frequencies
    f_c::Int64
    # spatial max val.
    x_max::Float64
    # ???
    num_x::Int64
    # Weights on the allowed frequencies.
    filter_x::Array{Float64}
end
type DynamicParameters
    # 2K+1 is the number of time samples.
    K::Int64
    # sampling rate
    tau::Float64
    # maxium allowed speed
    v_max::Float64
    # ????
    num_v::Int64
end

### Static parameters and container for use.
f_c = 20
x_max = 1.0
num_x = 10*f_c
filter_x =  ones(2*f_c+1);
static_parameters = StaticParameters(f_c, x_max, num_x, filter_x)

model_static = SuperResModels.Fourier1d(static_parameters.f_c, static_parameters.x_max, static_parameters.filter_x, static_parameters.num_x)
### Dynamic parameters and container for use.
K = 2
tau = 1.0/(K*1.0)
# the maximum velocity on the domain is given by x_max/2 + K tau v_max = 1
v_max = 0.5
num_v = 10
dynamic_parameters = DynamicParameters(K, tau, v_max, num_v)
model_dynamic = SuperResModels.DynamicFourier1d(model_static, dynamic_parameters.v_max, dynamic_parameters.tau, dynamic_parameters.K, dynamic_parameters.num_v*dynamic_parameters.K)
### Number of points
@everywhere begin
    model_static = $model_static
    model_dynamic = $model_dynamic
    x_max = $x_max
    v_max = $v_max
    test_case = () -> TestCases.cloud_1d_full(x_max, v_max, 0.1,10,K, tau, 4)
    noises_data = [0.01, 0.05, 0.1]
    #noises_data = linspace(0, 0.1, 5)
    #noises_data = noises_data[2:end]
    noises_position = [0.0001, 0.001]
    #noises_position = linspace(0, 2.5e-3, 5)
    #noises_position = noises_position[2:end]
end

# Replicating the process of Utils.generate_and_reconsturct_all step by step.
#results = [Utils.generate_and_reconstruct_all(model_static, model_dynamic, test_case, noises_data, noises_position)]
#println(" the results array is:", results )
#println(" with size : ", size(results))

(thetas,weights) = test_case()
# We plot the generated particles and their respective weights.
# This is plotted in space x velocity plane, with the size of the particles
# indicating the weights of the respective particles.
sizeAmp = 50
plt.scatter(thetas[1,:], thetas[2,:], s = weights*sizeAmp, alpha = 0.8)
plt.xlim((0, x_max))
plt.ylim((-v_max, v_max))
#plt.show()

results = zeros(1 + length(noises_data) + length(noises_position), 7)
# Obtain the dynamical norm of the configuration of particles
norm = minimum([abs(thetas[1, i] - thetas[1, j]) + model_dynamic.times[end]* abs(thetas[2, i] - thetas[2, j]) + (i==j) for i in 1:size(thetas, 2), j in 1:size(thetas, 2)])
norm2 = Utils.separation_val(thetas,K,tau,x_max)
# Obtain the minimum separation in space
dx = minimum([abs(thetas[1, i] - thetas[1, j]) + (i==j) for i in 1:size(thetas, 2), j in 1:size(thetas, 2)])
# Obtain the minimum separation in velocity
dv = minimum([abs(thetas[2, i] - thetas[2, j]) + (i==j) for i in 1:size(thetas, 2), j in 1:size(thetas, 2)])

println(" Results: ",results)
println(" Norm: ", norm)
println(" Norm2: ", norm2)
println(" dx infinite norm: ", dx)
println(" dv infinite norm: ", dv)

# Dismembering the function Utils.generate_and_reconstruct_dynamic(model_dynamic, thetas, weights, 0.0,0.0)
(thetas_est, weights_est) = Utils.run_simulation(model_dynamic, thetas, weights)
# See the estimated particles and weights, compared to the real ones
plt.scatter(thetas[1,:], thetas[2,:], s = weights*sizeAmp, alpha = 0.8)
plt.scatter(thetas_est[1,:], thetas_est[2,:], s=weights_est*sizeAmp/3, alpha = 0.8)
plt.xlim((0, x_max))
plt.ylim((-v_max, v_max))
plt.show()

# Now we obtain the difference with the actual solution

# We set a threshold weight in which we discard the generated particles.
(dist_x, dist_v, dist_w) = Utils.generate_and_reconstruct_dynamic(model_dynamic, thetas, weights, 0.0, 0.0)
println("Recons Dx = ", dist_x)
println("Recons Dv = ", dist_v)
println("Recons Dw = ", dist_w)

(dx, dv, norm, results) = Utils.generate_and_reconstruct_all(model_static, model_dynamic, test_case, noises_data, noises_position)
println("The noiseless results are: ", results[1,:])
println("Measurement noise results :", results[2,:])
println("position noise results :", results[3:end,:])
# We test the whole results function


text:
module annsp
    values = Dict()
    #inits
    __init__ = function(data, labels)
                    ver,len_ = annsp.ver__dim(data)
                    if ver
                        values["x"] = data
                        values["y"] = labels
                        values["eta"] = 0.01
                        values["weights"] = randn(len_+1)
                        annsp.train()
                    else
                        print("no")
                    end
                end

    #verifies data structure and dimensions
    ver__dim = function(x)
                    random_pick = x[rand(1:length(x))]
                    len_k = 0
                    for a in x
                        len_k = len_k + length(a)
                    end
                    if int(len_k/length(x)) == length(random_pick)
                        return true, length(random_pick)
                    else
                        return false
                    end
               end

     #binary signum
     signum = function(val)
                if val >= 0
                    return 1
                else
                    return -1
                end
              end

     #adder function
     adder = function(w,v)
                return dot(w,v)
             end

     #weight vector adaptation method for updating weights
     vector_adapt = function(current_weights, x, expected_decision, decision)
                        return current_weights + x*(values["eta"]*(expected_decision-decision))
                    end

     #train method
     train = function()
                   n_instances = length(values["x"])
                   for i in linspace(1,n_instances, n_instances)
                        insert!(values["x"][int(i)], 1, values["y"][int(i)])
                   end

                   for i in linspace(1,20,20)
                        for j in linspace(1, n_instances, n_instances)
                            expected_decision = values["y"][int(j)]
                            decision = annsp.signum(annsp.adder(values["weights"], values["x"][int(j)]))

                            if decision  != expected_decision
                                values["weights"] = annsp.vector_adapt(values["weights"], values["x"][int(j)], expected_decision, decision)
                            end

                        end
                   end
             end

      #to get the final calculated weights
      get_weights = function()
                        return values["weights"]
                    end

      #simpel classifier
      classify = function(v)
                    j = insert!(v, 1, 1)
                    decision = annsp.signum(annsp.adder(values["weights"], j))
                    return decision
                 end

    return {__init__ = __init__, classify=classify, get_weights = get_weights}
end


text:
"""
  The ring ZZ<x> := {c f/g | c in ZZ, f, g in ZZ[x], primitive, pos leading coeff}
  is a PID, even euclidean

  The key interest is
  IntCls(Z[x], F) = IntCls(Z<x>, F) cap IntCls(Q[x], F)
  Since Z<x> and Q[x] are PID (even euclidean) the last 2 can be
  computed using the Round2

  The name is bad, but stems from Florian Hess' PhD thesis,
  Chapter 2.10
  (Actually he covers that for a Dedekind ring R we have that
   R<x> is also Dedekind. The Round2, fundamentally would work
   for Dedekind rings, using PMats)
"""
module HessQRModule
using Hecke
import AbstractAlgebra, Nemo
import Base: +, -, *, gcd, lcm, divrem, div, rem, mod, ^, ==
export HessQR
import AbstractAlgebra: expressify

struct HessQR <: AbstractAlgebra.Ring
  R::FmpzPolyRing
  Qt::Generic.RationalFunctionField{fmpq}

  function HessQR(R::FmpzPolyRing, Qt::Generic.RationalFunctionField)
    new(R, Qt)
  end
end

function Base.show(io::IO, R::HessQR)
  print(io, AbstractAlgebra.obj_to_string(R, context = io))
end

function expressify(R::HessQR; context = nothing)
  return Expr(:sequence, Expr(:text, "extended valuation ring over "),
                         expressify(R.R, context = context),
                         Expr(:text, " in "),
                         expressify(R.Qt, context = context))
end

mutable struct HessQRElem <: RingElem
  parent::HessQR
  c::fmpz
  f::fmpz_poly
  g::fmpz_poly

  function HessQRElem(P::HessQR, c::fmpz, f::fmpz_poly, g::fmpz_poly)
    if iszero(c) || iszero(f)
      r = new(P, fmpz(0), zero(P.R), one(P.R))
      @assert parent(r.f) == P.R
      @assert parent(r.g) == P.R
      return r
    end
    if parent(f) != P.R
      f = map_coefficients(FlintZZ, f, parent = P.R)
    end
    if parent(g) != P.R
      g = map_coefficients(FlintZZ, g, parent = P.R)
    end
    gc = gcd(f, g)
    f = divexact(f, gc)
    g = divexact(g, gc)
    cf = content(f)*sign(leading_coefficient(f))
    cg = content(g)*sign(leading_coefficient(g))
    @assert (c*cf) % cg == 0
    r = new(P, divexact(c*cf, cg), divexact(f, cf), divexact(g, cg))
    @assert parent(r.f) == P.R
    @assert parent(r.g) == P.R
    return r
  end
  function HessQRElem(P::HessQR, q::Generic.Rat{fmpq})
    f = numerator(q)
    g = denominator(q)
    return HessQRElem(P, f, g)
  end
  function HessQRElem(P::HessQR, f::fmpq_poly)
    return HessQRElem(P, f, one(parent(f)))
  end
  function HessQRElem(P::HessQR, f::fmpq_poly, g::fmpq_poly)
    df = reduce(lcm, map(denominator, coefficients(f)), init = fmpz(1))
    dg = reduce(lcm, map(denominator, coefficients(g)), init = fmpz(1))
    ff = map_coefficients(FlintZZ, df*f, parent = P.R)
    gg = map_coefficients(FlintZZ, dg*g, parent = P.R)
    #ff/df//gg/dg = dg/df * ff/gg
    return HessQRElem(P, divexact(dg, df), ff, gg)
  end
  function HessQRElem(P::HessQR, c::fmpz)
    r = new(P, c, one(P.R), one(P.R))
    @assert parent(r.f) == P.R
    @assert parent(r.g) == P.R
    return r
  end
  function HessQRElem(P::HessQR, c::fmpz_poly)
    d = content(c)*sign(leading_coefficient(c))
    r = new(P, d, divexact(c, d), one(P.R))
    @assert parent(r.f) == P.R
    @assert parent(r.g) == P.R
    return r
  end
end

function Base.show(io::IO, a::HessQRElem)
  print(io, AbstractAlgebra.obj_to_string(a, context = io))
end

function expressify(a::HessQRElem; context = nothing)
  return  Expr(:call, :*, expressify(a.c, context = context),
             Expr(:call, ://, expressify(a.f, context = context),
                              expressify(a.g, context = context)))
end

check_parent(a::HessQRElem, b::HessQRElem) = parent(a) == parent(b) || error("Incompatible rings")

function Hecke.integral_split(a::Generic.Rat{fmpq}, S::HessQR)
  if iszero(a)
    return zero(S), one(S)
  end
  n = numerator(a)
  d = denominator(a)
  dn = reduce(lcm, map(denominator, coefficients(n)), init = fmpz(1))
  dd = reduce(lcm, map(denominator, coefficients(d)), init = fmpz(1))
  zn = S.R(n*dn)
  zd = S.R(d*dd)
  cn = content(zn)
  cd = content(zd)
  zn = divexact(zn, cn)
  zd = divexact(zd, cd)
  cn *= dd
  cd *= dn
  g = gcd(cn, cd)
  cn = divexact(cn, g)
  cd = divexact(cd, g)
  return HessQRElem(S, cn, zn, zd), S(cd)
end

function Hecke.numerator(a::Generic.Rat, S::HessQR)
  return integral_split(a, S)[1]
end

function Hecke.denominator(a::Generic.Rat, S::HessQR)
  return integral_split(a, S)[2]
end

Nemo.elem_type(::HessQR) = HessQRElem
Nemo.elem_type(::Type{HessQR}) = HessQRElem
Nemo.parent_type(::HessQRElem) = HessQR
Nemo.parent_type(::Type{HessQRElem}) = HessQR
Nemo.isdomain_type(::Type{HessQRElem}) = true

Base.parent(a::HessQRElem) = a.parent

(R::HessQR)(a::Generic.Rat{fmpq}) = HessQRElem(R, a)
(R::HessQR)(a::fmpz) = HessQRElem(R, a)
(R::HessQR)(a::Integer) = HessQRElem(R, fmpz(a))
(R::HessQR)(a::fmpz_poly) = HessQRElem(R, a)
(R::HessQR)(a::fmpq_poly) = HessQRElem(R, a)
(R::HessQR)(a::HessQRElem) = a
(R::HessQR)() = R(0)

(F::Generic.RationalFunctionField)(a::HessQRElem) = a.c*F(a.f)//F(a.g)

Nemo.iszero(a::HessQRElem) = iszero(a.c)
Nemo.isone(a::HessQRElem) = isone(a.c) && isone(a.f) && isone(a.g)

Nemo.zero(R::HessQR) = R(0)
Nemo.one(R::HessQR) = R(1)
Nemo.canonical_unit(a::HessQRElem) = HessQRElem(parent(a), fmpz(sign(a.c)), a.f, a.g)

Base.deepcopy_internal(a::HessQRElem, dict::IdDict) = HessQRElem(parent(a), Base.deepcopy_internal(a.c, dict), Base.deepcopy_internal(a.f, dict), Base.deepcopy_internal(a.g, dict))

Base.hash(a::HessQRElem, u::UInt=UInt(12376599)) = hash(a.g, hash(a.f, hash(a.c, u)))

+(a::HessQRElem, b::HessQRElem) = check_parent(a, b) && HessQRElem(parent(a), fmpz(1), a.c*a.f*b.g+b.c*b.f*a.g, a.g*b.g)
-(a::HessQRElem, b::HessQRElem) = check_parent(a, b) && HessQRElem(parent(a), fmpz(1), a.c*a.f*b.g-b.c*b.f*a.g, a.g*b.g)
-(a::HessQRElem) = HessQRElem(parent(a), -a.c, a.f, a.g)
*(a::HessQRElem, b::HessQRElem) = check_parent(a, b) && HessQRElem(parent(a), a.c*b.c, a.f*b.f, a.g*b.g)

==(a::HessQRElem, b::HessQRElem) = check_parent(a, b) && a.c*a.f == b.c *b.f && a.g == b.g

Base.:^(a::HessQRElem, n::Int) = HessQRElem(parent(a), a.c^n, a.f^n, a.g^n)

function Hecke.mul!(a::HessQRElem, b::HessQRElem, c::HessQRElem)
  d = b*c
  @assert parent(a.f) == parent(d.f)
  @assert parent(a.g) == parent(d.g)
  a.c = d.c
  a.f = d.f
  a.g = d.g
  return a
end

function Hecke.add!(a::HessQRElem, b::HessQRElem, c::HessQRElem)
  d = b+c
  @assert parent(a.f) == parent(d.f)
  @assert parent(a.g) == parent(d.g)
  a.c = d.c
  a.f = d.f
  a.g = d.g
  return a
end

function Hecke.addeq!(a::HessQRElem, b::HessQRElem)
  d = a+b
  @assert parent(a.f) == parent(d.f)
  @assert parent(a.g) == parent(d.g)
  a.c = d.c
  a.f = d.f
  a.g = d.g
  return a
end

function divrem(a::HessQRElem, b::HessQRElem)
  check_parent(a, b)
  if iszero(b)
    error("div by 0")
  end
  if iszero(a)
    return a, a
  end
  #= af/g mod b:
    it is enough to figure this out for b in Z>0, the rest is units
    that will disappaear into the quotient q

    - c = cont(g mod d), then since cont(g) is 1, c is a unit mod d
    - cd = 1 mod d, then
    - cont(g(cx^?+1) mod d) = 1 if ? = deg(g)+1 (or larger):
      g(cx?+1) = cg x^? + g and cont(cg mod d = 1)...
    - af/g -d/(cx^?+1) now has a denominator with cont( mod d) = 1
    - af/g - (af - du)/(g - dv) =
      (af(g-dv) -g(af-du))/(g(g-dv)) = d*..., so once cont(g %d) =1,
      we can replace f and g mod d (and figure out the quotient afterwards)

    - for d = p a prime, the rep is unqiue, thus F_p(t)

    - other approach:
      af/g = (af mod b)/(g mod b) mod b:

      (af + bA)  af    (af + bA)g - af       bAg
      -------- - -- =  --------------- =  ---------- = b * ....
        g + bB    g     (g + bB) g        (g + bB) g

      the tricky thing is that reducing af and g mod b they may no longer be
      coprime... and gcd mod b is non-trivial to compute. gcd_sircana might work
  =#
  r = rem(a,b)
  return divexact(a-r, b), r

  return HessQRElem(parent(a), q, a.f*b.g, a.g*b.f), HessQRElem(parent(a), r, a.f, a.g)
end

function div(a::HessQRElem, b::HessQRElem)
  check_parent(a, b)
  if iszero(a)
    return a
  end
  return divexact(a-rem(a, b), b)
  q = div(a.c, b.c)
  return HessQRElem(parent(a), q, a.f*b.g, a.g*b.f)
end

function rem(a::HessQRElem, b::HessQRElem)
#  @show :rem, a, b
  #NOT unique...., e.g. gcd(f,g mod b) might be <> 1....
  check_parent(a, b)
  if iszero(a)
    return a
  end
  #see above for reasoning
  d = abs(b.c)
  if isone(d)
    return zero(parent(a))
  end
  R = parent(a).R
  F, mF = quo(ZZ, d)
  aa = map_coefficients(mF, a.c*a.f)
  if iszero(aa)
    z = mF(one(domain(mF)))
  else
    z, aa = Hecke.primsplit!(aa)
  end
  bb = map_coefficients(mF, a.g, parent = parent(aa))
  _, f, g = Hecke.gcd_sircana(aa, bb)
  f *= z
  gg = lift(R, g)
  cg = content(gg)
  if !isone(cg)
    c = inv(canonical_unit(mF(cg)))
    gg = lift(R, c*g)
    @assert isone(content(gg))
    f *= c
  end
  ff = lift(R, f)
  cf = content(ff)
  if !iszero(cf)
    ff = divexact(ff, cf)
  end
  r = HessQRElem(parent(a), cf, ff, gg)
  @assert r.c < b.c
#  divexact(a-r, b)
  return r
end

function Nemo.divexact(a::HessQRElem, b::HessQRElem; check::Bool = false)
  check_parent(a, b)
#  @show :divexact, a, b
  q = HessQRElem(parent(a), divexact(a.c, b.c; check = true), a.f*b.g, a.g*b.f)
  @assert q*b == a
  return q
end

function gcd(a::HessQRElem, b::HessQRElem)
  return HessQRElem(parent(a), gcd(a.c, b.c))
end

function Nemo.gcdx(a::HessQRElem, b::HessQRElem)
  check_parent(a, b)
  R = parent(a)
  g, u, v = Nemo.gcdx(a.c, b.c)
  q,w, e =  R(g), HessQRElem(R, u, a.g, a.f), HessQRElem(R, v, b.g, b.f)
  @assert q == w*a+e*b
  return q, w, e
end

function lcm(a::HessQRElem, b::HessQRElem)
  check_parent(a, b)
  return HessQRElem(parent(a), lcm(a.c, b.c))
end

Hecke.isunit(a::HessQRElem) = isunit(a.c)

Nemo.dense_poly_type(::Type{gfp_fmpz_elem}) = gfp_fmpz_poly

function Nemo.ResidueField(a::HessQR, b::HessQRElem)
  @assert parent(b) == a
  @assert isprime(b.c)
  F = GF(b.c)
  Ft, t = RationalFunctionField(F, String(var(a.R)), cached = false)
  R = parent(numerator(t))
  return Ft, MapFromFunc(x->F(x.c)*Ft(map_coefficients(F, x.f, parent = R))//Ft(map_coefficients(F, x.g, parent = R)),
                         y->HessQRElem(a, fmpz(1), map_coefficients(lift, numerator(y)), map_coefficients(lift, denominator(y))), a, Ft)
end

function Nemo.ResidueRing(a::HessQR, b::HessQRElem)
  F = ResidueRing(FlintZZ, b.c)
  Fx, x = PolynomialRing(F, cached = false)
  Q = FractionField(Fx, cached = false)
  return Q, MapFromFunc(
     x->map_coefficients(F, x.c*x.f, parent = Fx)//map_coefficients(F, x.g, parent = Fx),
     y->HessQRElem(a, fmpz(1), lift(parent(b.f), numerator(y, false)), lift(parent(b.f), denominator(y, false))),
     a, Q)
end

function +(a::FracElem{T}, b::FracElem{T}) where T <: PolyElem{<:ResElem{S}} where S <: Hecke.IntegerUnion
  na = numerator(a, false)
  da = denominator(a, false)

  nb = numerator(b, false)
  db = denominator(b, false)

  nc = na*db + da*nb
  dc = da*db
  _, da, db = Hecke.gcd_sircana(nc, dc)
  return parent(a)(da, db)
end

function -(a::FracElem{T}, b::FracElem{T}) where T <: PolyElem{<:ResElem{S}} where S <: Hecke.IntegerUnion
  na = numerator(a, false)
  da = denominator(a, false)

  nb = numerator(b, false)
  db = denominator(b, false)

  nc = na*db - da*nb
  dc = da*db
  _, da, db = Hecke.gcd_sircana(nc, dc)
  return parent(a)(da, db)
end

function *(a::FracElem{T}, b::FracElem{T}) where T <: PolyElem{<:ResElem{S}} where S <: Hecke.IntegerUnion
  na = numerator(a, false)
  da = denominator(a, false)

  nb = numerator(b, false)
  db = denominator(b, false)

  nc = na*nb
  dc = da*db
  _, da, db = Hecke.gcd_sircana(nc, dc)
  return parent(a)(da, db)
end

function Hecke.factor(a::HessQRElem)
  f = factor(a.c)
  R = parent(a)
  return Fac(R(a.f), Dict((R(p),k) for (p,k) = f.fac))
end

function Hecke.factor(a::Generic.Rat, R::HessQR)
  d1 = reduce(lcm, map(denominator, coefficients(numerator(a))), init = fmpz(1))
  f1 = factor(R(d1*numerator(a)))
  d2 = reduce(lcm, map(denominator, coefficients(denominator(a))), init = fmpz(1))
  f2 = factor(R(d1*denominator(a)))

  for (p,k) = f2.fac
    if haskey(f1.fac, p)
      f1.fac[p] -= k
    else
      f1.fac[p] = k
    end
  end
  f1.unit = divexact(f1.unit, f2.unit)
  return f1
end

function Hecke.isconstant(a::HessQRElem)
  return iszero(a) || (isconstant(a.f) && isconstant(a.g))
end

end

using .HessQRModule


text:
include("../src/mesh.jl")
include("../src/advection.jl")
include("../src/lagrange.jl")
# include("../src/lagrange2d.jl")
# include("../src/advection2d.jl")
include("../src/bspline.jl")
include("../src/bsplinelu.jl")
include("../src/bsplinefft.jl")
include("../src/interpolation.jl")

using LinearAlgebra
using Plots

function exact(tf::T, mesh1::UniformMesh{T}, mesh2::UniformMesh{T}) where {T}

    f = zeros(T, (mesh1.length, mesh2.length))
    for (i, x) in enumerate(mesh1.points), (j, y) in enumerate(mesh2.points)
        xn = cos(tf) * x + sin(tf) * y
        yn = -sin(tf) * x + cos(tf) * y
        f[i, j] = exp(-6 * ((xn)^2 + (yn + big"1.8")^2))
    end

    return f

end

function fctadv(interp, mesh1, mesh2, v1, v2, refdeb, refend, dt)
    adv_x1 = Advection1d(mesh1, interp)
    adv_x2 = Advection1d(mesh2, interp)

    f = copy(refdeb)
    ft = copy(refdeb)
    advection!(adv_x1, f, v1, tan(dt / 2))
    transpose!(ft, f)
    advection!(adv_x2, ft, v2, sin(dt))
    transpose!(f, ft)
    advection!(adv_x1, f, v1, tan(dt / 2))

    return norm(f - refend, Inf)
end

function fctmain(sz, dt::T, ordmax) where {T}

    mesh1 = UniformMesh(T(-5), T(5), sz)
    mesh2 = UniformMesh(T(-5), T(5), sz)

    refdeb = exact(zero(T), mesh1, mesh2)
    refend = exact(dt, mesh1, mesh2)

    v1 = -collect(mesh2.points)
    v2 = +collect(mesh1.points)

    indmax = div(ordmax - 3, 2) + 1

    y = ones(indmax, 3)
    x = zeros(Int64, indmax)

    labels = Array{String,2}(undef, 1, 3)
    labels[1, 1] = "Lagrange"
    labels[1, 2] = "B-Spline LU"
    labels[1, 3] = "B-Spline FFT"

    ind = 1

    for order = 3:2:ordmax
        println("order=$order")
        lag = Lagrange(BigFloat, order)
        splu = B_SplineLU(order, mesh1.length, zero(T))
        spfft = B_SplineFFT(order, mesh1.length, zero(T))

        y[ind, 1] = fctadv(lag, mesh1, mesh2, v1, v2, refdeb, refend, dt)
        y[ind, 2] = fctadv(splu, mesh1, mesh2, v1, v2, refdeb, refend, dt)
        y[ind, 3] = fctadv(spfft, mesh1, mesh2, v1, v2, refdeb, refend, dt)

        println("ind=$ind order=$order y[ind,:]=$(y[ind,:])")

        x[ind] = order

        ind += 1
    end

    p = Plots.plot(
        x,
        log10.(y),
        xlabel = "order",
        ylabel = "error",
        legend = :bottomleft,
        label = labels,
        marker = 2,
    )
    prec = precision(BigFloat)

    Plots.savefig(p, "out/result2_$(sz)_$(prec)_$(ordmax).pdf")
end

function fctinter(interp::InterpolationType, sz, dec)

    fct(x) = cos(2big(pi) * x)
    mesh = big.(1:sz) / sz
    deb = fct.(mesh)
    ref = fct.(mesh .+ dec / sz)
    res = zeros(BigFloat, sz)
    interpolate!(res, deb, dec, interp)
    return Float64(norm(ref - res, Inf))
end

function fctmain2(sz, ordmax)

    indmax = ordmax - 2

    y = ones(indmax, 2)
    x = zeros(Int64, indmax)

    labels = Array{String,2}(undef, 1, 2)
    labels[1, 1] = "Lagrange"
    labels[1, 2] = "B-Spline LU"
    #    labels[1,3] = "B-Spline FFT"

    ind = 1

    for order = 3:ordmax
        println("order=$order")
        lag = Lagrange(BigFloat, order)
        splu = B_SplineLU(order, sz, big"0.")
        #        spfft = B_SplineFFT(order, sz, big"0.")

        y[ind, 1] = fctinter(lag, sz)
        y[ind, 2] = fctinter(splu, sz)
        #        y[ind,3] = fctinter(spfft, sz)

        println("ind=$ind order=$order y[ind,:]=$(y[ind,:])")

        x[ind] = order

        ind += 1
    end

    p = Plots.plot(
        x,
        log10.(y),
        xlabel = "order",
        ylabel = "error",
        legend = :bottomleft,
        label = labels,
        marker = 2,
    )
    prec = precision(BigFloat)

    Plots.savefig(p, "out/resinterp_$(sz)_$(prec)_$(ordmax).pdf")
end

function formule(order, sz, dec)
    return Float64(
        (big(pi) / sz)^(order + 1) * 4 * dec * (1 - dec) / sqrt(big(pi) * (order + 1) / 2),
    )
end

function fctmain_gnuplot(sz, ordmax)
    prec = precision(BigFloat)
    println("# begin gnuplot size=$sz julia precision=$prec")
    println("# order\tLagrange\tformule\t(formule-lag)\tB_SplineLU")
    dec = big"0.351726155665655665187291927162514231451"
    for order = 3:ordmax
        lag = Lagrange(BigFloat, order)
        splu = B_SplineLU(order, sz, big"0.")
        reslag = fctinter(lag, sz, dec)
        ressplu = fctinter(splu, sz, dec)
        resformule = formule(order, sz, dec)
        resdiff = resformule - reslag
        println("$order\t$reslag\t$resformule\t$resdiff\t$ressplu")
    end
end
function fctmainsize_gnuplot(szlist, orderlist)

    prec = precision(BigFloat)
    println("# begin gnuplot order=$orderlist julia precision=$prec")
    print("# size")
    for order in orderlist
        print("\tLagrange($order)\tB_Spline($order)")
    end
    println("")
    dec = big"0.351726155665655665187291927162514231451"
    for sz in szlist
        print("$sz")
        for order in orderlist
            lag = Lagrange(BigFloat, order)
            splu = B_SplineLU(order, sz, big"0.")
            reslag = fctinter(lag, sz, dec)
            ressplu = fctinter(splu, sz, dec)
            print("\t$reslag\t$ressplu")
        end
        println("")
    end
end

setprecision(1024) do
    fctmain_gnuplot(1001, 51)
    #    fctmainsize_gnuplot(101:100:10001,[3,11,31])
end


text:
module BoundaryTriangulationsTests

using Test
using Gridap.Helpers
using Gridap.Fields
using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry
using FillArrays

domain = (0,4,0,4)
partition = (2,2)
model = simplexify(CartesianDiscreteModel(domain,partition))

btrian = BoundaryTriangulation(model,tags=[7,8])
test_triangulation(btrian)
@test get_background_model(btrian) === model
glue = get_glue(btrian,Val(1))
@test glue.tface_to_mface === btrian.glue.face_to_bgface
glue = get_glue(btrian,Val(2))
@test glue.tface_to_mface === btrian.glue.face_to_cell

face_s_q = glue.tface_to_mface_map

s1 = Point(0.0)
s2 = Point(0.5)
s = [s1,s2]
face_to_s = Fill(s,length(face_s_q))

face_to_q = lazy_map(evaluate,face_s_q,face_to_s)
@test isa(face_to_q,Geometry.FaceCompressedVector)

cell_grid = get_grid(get_background_model(btrian))
cell_shapefuns = get_cell_shapefuns(cell_grid)
cell_grad_shapefuns = lazy_map(Broadcasting(∇),cell_shapefuns)

face_shapefuns = lazy_map(Reindex(cell_shapefuns),glue.tface_to_mface)
face_grad_shapefuns = lazy_map(Reindex(cell_grad_shapefuns),glue.tface_to_mface)

face_shapefuns_q = lazy_map(evaluate,face_shapefuns,face_to_q)
test_array(face_shapefuns_q,collect(face_shapefuns_q))
@test isa(face_shapefuns_q,Geometry.FaceCompressedVector)

face_grad_shapefuns_q = lazy_map(evaluate,face_grad_shapefuns,face_to_q)
test_array(face_grad_shapefuns_q,collect(face_grad_shapefuns_q))
@test isa(face_grad_shapefuns_q,Geometry.FaceCompressedVector)

face_to_nvec = get_facet_normal(btrian)
face_to_nvec_s = lazy_map(evaluate,face_to_nvec,face_to_s)
test_array(face_to_nvec_s,collect(face_to_nvec_s))

##print_op_tree(face_shapefuns)
#print_op_tree(face_to_q)
#print_op_tree(face_shapefuns_q)
#print_op_tree(face_grad_shapefuns_q)
#print_op_tree(face_to_nvec_s)

Ω = Triangulation(model)
Γ = btrian
∂Γ = BoundaryTriangulation(Γ)
test_triangulation(∂Γ)
@test ∂Γ.rtrian === Γ
@test isa(∂Γ.dtrian,BoundaryTriangulation)
glue = get_glue(∂Γ,Val(2))
glue = get_glue(∂Γ,Val(1))
@test is_change_possible(Ω,∂Γ)
@test is_change_possible(Γ,∂Γ)

domain = (0,4,0,4)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

btrian = BoundaryTriangulation(model)
test_triangulation(btrian)
@test get_background_model(btrian) === model
glue = get_glue(btrian,Val(1))
@test glue.tface_to_mface === btrian.glue.face_to_bgface
glue = get_glue(btrian,Val(2))
@test glue.tface_to_mface === btrian.glue.face_to_cell

@test isa(get_cell_map(btrian)[1],AffineMap)

face_s_q = glue.tface_to_mface_map

s1 = Point(0.0)
s2 = Point(0.5)
s = [s1,s2]
face_to_s = Fill(s,length(face_s_q))

face_to_q = lazy_map(evaluate,face_s_q,face_to_s)
@test isa(face_to_q,Geometry.FaceCompressedVector)

cell_grid = get_grid(get_background_model(btrian))
cell_shapefuns = get_cell_shapefuns(cell_grid)
cell_grad_shapefuns = lazy_map(Broadcasting(∇),cell_shapefuns)

face_shapefuns = lazy_map(Reindex(cell_shapefuns),glue.tface_to_mface)
face_grad_shapefuns = lazy_map(Reindex(cell_grad_shapefuns),glue.tface_to_mface)

face_shapefuns_q = lazy_map(evaluate,face_shapefuns,face_to_q)
test_array(face_shapefuns_q,collect(face_shapefuns_q))
@test isa(face_shapefuns_q,Geometry.FaceCompressedVector)

face_grad_shapefuns_q = lazy_map(evaluate,face_grad_shapefuns,face_to_q)
test_array(face_grad_shapefuns_q,collect(face_grad_shapefuns_q))
@test isa(face_grad_shapefuns_q,Geometry.FaceCompressedVector)

face_to_nvec = get_facet_normal(btrian)
face_to_nvec_s = lazy_map(evaluate,face_to_nvec,face_to_s)
test_array(face_to_nvec_s,collect(face_to_nvec_s))
@test isa(face_to_nvec_s,Geometry.FaceCompressedVector)

#print_op_tree(face_shapefuns_q)
#print_op_tree(face_grad_shapefuns_q)
#print_op_tree(face_to_nvec_s)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))

s2q = glue.tface_to_mface_map

q = lazy_map(evaluate,s2q,s)
r = Vector{Point{2,Float64}}[
  [(0.25,0.0),(0.75,0.0)],[(0.0,0.25),(0.0,0.75)],
  [(0.25,0.0),(0.75,0.0)],[(1.0,0.25),(1.0,0.75)],
  [(0.25,1.0),(0.75,1.0)],[(0.0,0.25),(0.0,0.75)],
  [(0.25,1.0),(0.75,1.0)],[(1.0,0.25),(1.0,0.75)]]
test_array(q,r)

q2x = get_cell_map(cell_grid)
s2x = lazy_map(∘,lazy_map(Reindex(q2x),glue.tface_to_mface),s2q)

x = lazy_map(evaluate,s2x,s)
r = Vector{Point{2,Float64}}[
  [(0.5,0.0),(1.5,0.0)],[(0.0,0.5),(0.0,1.5)],
  [(2.5,0.0),(3.5,0.0)],[(4.0,0.5),(4.0,1.5)],
  [(0.5,4.0),(1.5,4.0)],[(0.0,2.5),(0.0,3.5)],
  [(2.5,4.0),(3.5,4.0)],[(4.0,2.5),(4.0,3.5)]]
test_array(x,r)

nvec = get_facet_normal(btrian)

nvec_s = lazy_map(evaluate,nvec,s)
r = Vector{Point{2,Float64}}[
  [(0.0,-1.0),(0.0,-1.0)],[(-1.0,0.0),(-1.0,0.0)],
  [(0.0,-1.0),(0.0,-1.0)],[(1.0,-0.0),(1.0,-0.0)],
  [(0.0,1.0),(0.0,1.0)],[(-1.0,0.0),(-1.0,0.0)],
  [(0.0,1.0),(0.0,1.0)],[(1.0,-0.0),(1.0,-0.0)]]
test_array(nvec_s,r)

cellids = collect(1:num_cells(model))

face_to_cellid = lazy_map(Reindex(cellids),glue.tface_to_mface)
@test face_to_cellid == glue.tface_to_mface

q2x = get_cell_map(cell_grid)
s2q = glue.tface_to_mface_map
s2x = lazy_map(∘,lazy_map(Reindex(q2x),glue.tface_to_mface),s2q)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))
x = lazy_map(evaluate,s2x,s)
r = Vector{Point{2,Float64}}[
  [(0.5,0.0),(1.5,0.0)],[(0.0,0.5),(0.0,1.5)],
  [(2.5,0.0),(3.5,0.0)],[(4.0,0.5),(4.0,1.5)],
  [(0.5,4.0),(1.5,4.0)],[(0.0,2.5),(0.0,3.5)],
  [(2.5,4.0),(3.5,4.0)],[(4.0,2.5),(4.0,3.5)]]
test_array(x,r)

nvec = get_facet_normal(btrian)
nvec_x = lazy_map(evaluate,nvec,s)

domain = (0,4,0,4)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

btrian = BoundaryTriangulation(model,tags="tag_8")
test_triangulation(btrian)
@test all(1 .== btrian.glue.face_to_lcell)

btrian = BoundaryTriangulation(model,get_face_labeling(model),tags="tag_8")
test_triangulation(btrian)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))
nvec = get_facet_normal(btrian)
nvec_x = lazy_map(evaluate,nvec,s)
s2x = get_cell_map(btrian)
x = lazy_map(evaluate,s2x,s)

btrian = BoundaryTriangulation(model,tags=["tag_6","tag_7"])
test_triangulation(btrian)

ids = [1,3,4]
vtrian = view(btrian,ids)
test_triangulation(vtrian)
@test num_cells(vtrian) == 3
bglue = get_glue(btrian,Val(1))
vglue = get_glue(vtrian,Val(1))
@test vglue.tface_to_mface == bglue.tface_to_mface[ids]
bglue = get_glue(btrian,Val(2))
vglue = get_glue(vtrian,Val(2))
@test vglue.tface_to_mface == bglue.tface_to_mface[ids]

end # module


text:
"""
    Layers(layers::Vector{Layer})

Algebraic object encoding a list of [`AlgebraOfGraphics.Layer`](@ref) objects.
`Layers` objects can be added or multiplied, yielding a novel `Layers` object.
"""
struct Layers
    layers::Vector{Layer}
end

Base.convert(::Type{Layers}, s::Layer) = Layers([s])
Base.convert(::Type{Layers}, s::Layers) = s

Base.getindex(v::Layers, i::Int) = v.layers[i]
Base.length(v::Layers) = length(v.layers)
Base.eltype(::Type{Layers}) = Layer
Base.iterate(v::Layers, args...) = iterate(v.layers, args...)

const OneOrMoreLayers = Union{Layers, Layer}

function Base.:+(s1::OneOrMoreLayers, s2::OneOrMoreLayers)
    l1::Layers, l2::Layers = s1, s2
    return Layers(vcat(l1.layers, l2.layers))
end

function Base.:*(s1::OneOrMoreLayers, s2::OneOrMoreLayers)
    l1::Layers, l2::Layers = s1, s2
    return Layers([el1 * el2 for el1 in l1 for el2 in l2])
end

function compute_processedlayers_grid(processedlayers, categoricalscales)
    indices = CartesianIndices(compute_grid_positions(categoricalscales))
    pls_grid = map(_ -> ProcessedLayer[], indices)
    for processedlayer in processedlayers
        append_processedlayers!(pls_grid, processedlayer, categoricalscales)
    end
    return pls_grid
end

function compute_attributes(attributes, primary, named)
    attrs = NamedArguments()
    merge!(attrs, attributes)
    merge!(attrs, primary)
    merge!(attrs, named)

    # implement alpha transparency
    alpha = get(attrs, :alpha, automatic)
    color = get(attrs, :color, automatic)
    (color !== automatic) && (alpha !== automatic) && (color = (color, alpha))

    # opt out of the default cycling mechanism
    cycle = nothing

    merge!(attrs, Dictionary(valid_options(; color, cycle)))

    # remove unnecessary information
    return filterkeys(!in((:col, :row, :layout, :alpha)), attrs)
end

function compute_entries_continuousscales(pls_grid)
    # Here processed layers in `pls_grid` are "sliced",
    # the categorical scales have been applied, but not
    # the continuous scales

    entries_grid = map(_ -> Entry[], pls_grid)
    continuousscales_grid = map(_ -> MixedArguments(), pls_grid)

    for idx in eachindex(pls_grid), pl in pls_grid[idx]
        # Apply continuous transformations
        positional = map(contextfree_rescale, pl.positional)
        plottype = Makie.plottype(pl.plottype, positional...)

        # Compute continuous scales with correct plottype, to figure out role of color
        continuousscales = AlgebraOfGraphics.continuousscales(ProcessedLayer(pl; plottype))
        mergewith!(mergescales, continuousscales_grid[idx], continuousscales)

        # Compute `Entry` with rescaled columns
        named = compute_attributes(pl.attributes, pl.primary, map(contextfree_rescale, pl.named))
        entry = Entry(plottype, positional, named)
        push!(entries_grid[idx], entry)
    end

    # Compute merged continuous scales, as it may be needed to use global extrema
    merged_continuousscales = reduce(mergewith!(mergescales), continuousscales_grid, init=MixedArguments())

    colorscale = get(merged_continuousscales, :color, nothing)
    if !isnothing(colorscale)
        # TODO: might need to change to support temporal color scale
        colorrange = colorscale.extrema
        for entries in entries_grid, entry in entries
            # Safe to do, as each entry has a separate `named` dictionary
            set!(entry.named, :colorrange, colorrange)
        end
    end

    return entries_grid, continuousscales_grid, merged_continuousscales
end

function compute_axes_grid(fig, s::OneOrMoreLayers;
                           axis=NamedTuple(), palettes=NamedTuple())

    axes_grid = compute_axes_grid(s; axis, palettes)
    sz = size(axes_grid)
    if sz !== (1, 1) && fig isa Axis
        error("You can only pass an `Axis` to `draw!`, if the calculated layout only contains one element. Elements: $(sz)")
    end

    return map(ae -> AxisEntries(ae, fig), axes_grid)
end

function compute_axes_grid(s::OneOrMoreLayers;
                           axis=NamedTuple(), palettes=NamedTuple())
    layers::Layers = s
    processedlayers = map(ProcessedLayer, layers)

    theme_palettes = NamedTuple(Makie.current_default_theme()[:palette])
    palettes = merge((layout=wrap,), map(to_value, theme_palettes), palettes)

    categoricalscales = MixedArguments()
    for processedlayer in processedlayers
        mergewith!(
            mergescales,
            categoricalscales,
            AlgebraOfGraphics.categoricalscales(processedlayer, palettes)
        )
    end
    # fit categorical scales (compute plot values using all data values)
    map!(fitscale, categoricalscales, categoricalscales)

    pls_grid = compute_processedlayers_grid(processedlayers, categoricalscales)
    entries_grid, continuousscales_grid, merged_continuousscales =
        compute_entries_continuousscales(pls_grid)

    indices = CartesianIndices(pls_grid)
    axes_grid = map(indices) do c
        return AxisSpecEntries(
            AxisSpec(c, axis),
            entries_grid[c],
            categoricalscales,
            continuousscales_grid[c]
        )
    end

    # Axis labels and ticks
    for ae in axes_grid
        ndims = isaxis2d(ae) ? 2 : 3
        for (i, var) in zip(1:ndims, (:x, :y, :z))
            scale = get(ae.categoricalscales, i) do
                return get(ae.continuousscales, i, nothing)
            end
            isnothing(scale) && continue
            label = getlabel(scale)
            # Use global scales for ticks for now
            # TODO: requires a nicer mechanism that takes into account axis linking
            (scale isa ContinuousScale) && (scale = merged_continuousscales[i])
            for (k, v) in pairs((label=to_string(label), ticks=ticks(scale)))
                keyword = Symbol(var, k)
                # Only set attribute if it was not present beforehand
                get!(ae.axis.attributes, keyword, v)
            end
        end
    end

    return axes_grid
end


text:
module HomotopyContinuation
    include("HomConBase/HomConBase.jl")

    using .HomConBase

    # Exports defined in HomConBase
    export AbstractHomotopy
    export HomConAlgorithm
    export init, solve
    export evaluate, substitute, startsystem, targetsystem,
        differentiate, dt, ∂t, nvariables, degrees,
        nequations, ishomogenous, homogenize
    export affine, projective, totaldegree, randomsystem

    # Exports from FixedPolySystem in HomConBase
    export Poly, PolySystem,
        nvariables, variables, polynomials,
        evaluate, evaluate!, substitute, removepoly, differentiate,
        ishomogenous, homogenize, homogenized, dehomogenize,
        weyldot, weylnorm

    include("Homotopy/Homotopy.jl")

    using .Homotopy
    # Homotopy Exports
    export substitute, removepoly, valuate, startsystem, targetsystem, differentiate,
        dt, homogenize, homogenized, degrees, weylnorm, nvars, nequations
    export StraightLineHomotopy, GammaTrickHomotopy, FixedHomotopy

    include("PredictorCorrector/PredictorCorrector.jl")

    using .PredictorCorrector
    #PredictorCorrector exports
    export PredictorCorrector
    export correct!, predict
    export Spherical, Affine, trackpath, solve, cauchyendgame,
        Result, CauchyEndgameResult, PathResult, ConvergentCluster

    #
    include("TestSystems/TestSystems.jl")
    export TestSystems
end # module


text:
using SymPy: Sym, diff
using ForwardDiff: jacobian
include("kalmanflt.jl")

struct StateVariable

    pred_func::Any

    prev_state_sym::Sym

    pred_func_syms::Array{Sym}

end

struct ObservationVariable

    obs_func::Any

    obs_func_syms::Array{Sym}

end

function eval_jacobian_sym_mat(
        mat::Matrix{Sym},
        par_values_dict::Dict
)::Matrix{Float64}
   jacobian = zeros(Float64, size(mat))
    for row in 1:size(jacobian)[1]
        for col in 1:size(jacobian)[2]
            jacobian[row, col] = float(
                mat[row, col].subs(par_values_dict).evalf()
            )
        end
    end
    return jacobian
end

function eval(
        mat::Matrix{Sym},
        par_values_dict::Dict
)
    return eval_jacobian_sym_mat(mat, par_values_dict)
end

function observation_matrix_sym(
        observation_variables::Vector{ObservationVariable},
        state_variables::Vector{StateVariable}
)::Matrix{Sym}

    unevaluated_jacobian = Matrix{Sym}(undef, (length(observation_variables), length(state_variables)))
    for row in 1:length(observation_variables)
        obs_variable = observation_variables[row].obs_func
        obs_fcn_variables = observation_variables[row].obs_func_syms
        for col in 1:length(state_variables)
            state_n_pred_variable = state_variables[col].prev_state_sym
            deriv = diff(obs_variable(obs_fcn_variables...), state_n_pred_variable, evaluate=false)
            unevaluated_jacobian[row, col] = deriv
        end
    end
    return unevaluated_jacobian
end

function transition_matrix_sym(
        state_variables::Vector{StateVariable}
)::Matrix{Sym}

    unevaluated_jacobian = Matrix{Sym}(
        undef, (length(state_variables), length(state_variables))
    )
    for row in 1:length(state_variables)
        state_n_variable = state_variables[row].pred_func
        state_fcn_variables = state_variables[row].pred_func_syms
        for col in 1:length(state_variables)
            state_nm1_variable = state_variables[col].prev_state_sym
            deriv = diff(
                state_n_variable(state_fcn_variables...),
                state_nm1_variable,
                evaluate=false
            )
            unevaluated_jacobian[row, col] = deriv
        end
    end
    return unevaluated_jacobian
end

@inline function predict_state_cov(
        transition_mat::Matrix,
        prev_state_cov::Matrix,
        noise_mat::Matrix
)::Matrix
    transition_mat * prev_state_cov * transpose(transition_mat) + noise_mat
end

function predict_state(
        state_vars_vec::Vector{StateVariable},
        prev_state::KalmanState,
        model_vars::Dict,
        noise_mat::Matrix
)::KalmanState

    eval_dict = model_vars
    pred_state_vec = zeros(eltype(prev_state.state), size(state_vars_vec))
    for index in 1:length(state_vars_vec)
        key::Sym = state_vars_vec[index].prev_state_sym
        value = prev_state.state[index]
        eval_dict[key] = value
    end

    trans_mat::Matrix{eltype(pred_state_vec)} = eval(
        transition_matrix_sym(state_vars_vec),
        eval_dict
    )
    for index in 1:length(state_vars_vec)
        state_var::StateVariable = state_vars_vec[index]
        pred_func_sym::Sym = state_var.pred_func(state_var.pred_func_syms...)
        pred_state_vec[index] = pred_func_sym.subs(eval_dict).evalf()
    end
    pred_state_cov = predict_state_cov(trans_mat, prev_state.cov, noise_mat)
    return KalmanState{eltype(pred_state_cov)}(pred_state_vec, pred_state_cov)
end

function predict_state(
        state_function_vec::Any,
        prev_state::KalmanState,
        noise_mat::Matrix
)::KalmanState

    predicted_state_vec::Vector{eltype(prev_state.state)} = state_function_vec(prev_state.state)
    trans_mat::Matrix{eltype(prev_state.state)} = jacobian(
        state_function_vec,
        prev_state.state
    )
    pred_state_cov = predict_state_cov(trans_mat, prev_state.cov, noise_mat)
    return KalmanState{eltype(pred_state_cov)}(predicted_state_vec, pred_state_cov)
end

function correct_state(
        pred_state::KalmanState,
        meas_residual::Vector,
        measurement::KalmanObservation
)::KalmanState

    residual_cov = (
        measurement.obs * pred_state.cov * transpose(measurement.obs)
        + measurement.obsCov
    )
    gain = pred_state.cov * transpose(measurement.obs) * residual_cov^(-1)

    corr_state_vec = pred_state.state + gain * meas_residual
    corr_cov = (
        (I(length(corr_state_vec)) - gain * measurement.obs) * pred_state.cov
    )
    return KalmanState{eltype(corr_cov)}(corr_state_vec, corr_cov)
end

function correct_state(
        state_vars_vec::Vector{StateVariable},
        obs_vars_vec::Vector{ObservationVariable},
        model_vars::Dict,
        pred_state::KalmanState,
        measurement::KalmanObservation
)::KalmanState

    eval_dict = model_vars
    for index in 1:length(state_vars_vec)
        key::Sym = state_vars_vec[index].prev_state_sym
        value = pred_state.state[index]
        eval_dict[key] = value
    end

    obs_matrix::Matrix = measurement.obs
    if all(els -> els == 0, obs_matrix)
        obs_matrix = eval(
            observation_matrix_sym(
                obs_vars_vec,
                state_vars_vec
            ), eval_dict
        )
        measurement.obs = obs_matrix
    end
    meas_residual = zeros(eltype(measurement.meas), size(measurement.meas))
    for index in 1:length(obs_vars_vec)
        obs_var = obs_vars_vec[index]
        obs_pred = obs_var.obs_func(obs_var.obs_func_syms...).subs(eval_dict).evalf()
        meas_residual[index] = measurement.meas[index] - obs_pred
    end

    return correct_state(
        pred_state,
        meas_residual,
        measurement
    )
end

function correct_state(
        meas_function_vec::Any,
        pred_state::KalmanState,
        measurement::KalmanObservation
)::KalmanState

    meas_residual::Vector{eltype(pred_state.state)} = (
        measurement.meas - meas_function_vec(pred_state.state)
    )

    obs_matrix::Matrix = measurement.obs
    if all(els -> els == 0, obs_matrix)
        obs_matrix = jacobian(meas_function_vec, pred_state.state)
        measurement.obs = obs_matrix
    end

    return correct_state(
        pred_state,
        meas_residual,
        measurement
    )
end;

max_stars_repo_name: hoganman/kalmanfilterjl

max_issues_repo_name: hoganman/kalmanfilterjl

max_forks_repo_name: hoganman/kalmanfilterjl


text:
##
##  d e r i v a t i v e . j l  Numerical Derivatives
##

#-- Numerical gradient of multivariate function ------------------------------
#
function fd_gradient{T<:Real}(f::Function, x0::Array{T,1}; h::Real = 0.0)
    local heps::Real = h
    if heps == 0.0; heps = eps()^(1.0/3.0); end
    local n::Int = length(x0)
    if length(f(x0)) > 1
        error("Function 'f' must be a scalar function.")
    end

    local hh = zeros(n), gr = zeros(n)
    for i = 1:n
        hh[i] = heps
        gr[i] = ((f(x0 + hh) - f(x0 - hh)) / (2.0*heps))
        hh[i] = 0.0
    end
    return gr
end

function fd_gradient(f::Function, x0::Real; h::Real = 0.0)
    local heps::Real = h
    if heps == 0.0; heps = eps()^(1.0/3.0); end
    if length(f(x0)) > 1
        error("Function 'f' must be a scalar function.")
    end
    return (f(x0 + heps) - f(x0 - heps)) / (2.0*heps)
end

function fd_jacobian{T<:Real}(f::Function, x0::Array{T,1}; h::Real = 0.0)
    local heps::Real = h
    if heps == 0.0; heps = eps()^(1.0/3.0); end
	local n = length(x0), m = length(f(x0))
	local hh = zeros(n)
	jacob = zeros(m, n)
	for i = 1:n
		hh[i] = heps
		jacob[:, i] = (f(x0 + hh) - f(x0 - hh)) / (2.0*heps)
		hh[i] = 0.0
	end
	return jacob
end

function fd_hessian{T<:Real}(f::Function, x0::Array{T,1}; h::Real = 0.0)
    local heps::Real = h
    if heps == 0.0; heps = eps()^(1.0/4.0); end
    local n::Int = length(x0)
    if length(f(x0)) != 1
        error("Function 'f' must a scalar function of n variables.")
    end

    if n == 1
        H = (f(x0.-heps) - 2.0*f(x0) + f(x0.+heps)) / heps^2

    else
        H = zeros(n, n)
        hh = diagm(heps*ones(n))
        for i = 1:(n-1)
            hi = hh[:, i]
            H[i, i] = (f(x0-hi) - 2.0*f(x0) + f(x0+hi)) / heps^2
            for j = (i+1):n
                hj = hh[:, j]
                H[i, j] = (f(x0+hi+hj) - f(x0+hi-hj) - f(x0-hi+hj) + f(x0-hi-hj)) / (4.0*heps^2)
                H[j, i] = H[i, j]
            end
        end
        hi = hh[:, n]
        H[n, n] = (f(x0-hi) - 2.0*f(x0) + f(x0+hi)) / heps^2
    end

    return H
end

function fd_laplacian{T<:Real}(f::Function, x0::Array{T,1}; h::Real = 0.0)
    local heps::Real = h
    if heps == 0.0; heps = eps()^(1.0/4.0); end
    local n::Int = length(x0)
    local hh::Real = zeros(n)

    L = 0.0
    for i = 1:n
        hh[i] = heps
        L += (f(x0+hh) + f(x0-hh) - 2.0*f(x0)) / heps^2
        hh[i] = 0.0
    end

    return L
end

#-- Central difference combined with Richardson approximation ----------------
#
function numderiv(f::Function, x0::Real; n::Int = 16, h::Real = 0.1)

    local epsilon = eps(), err = Inf
    local hstep::Real = h

    j::Int = 1
    D = zeros(n, n)
    D[1, 1] = (f(x0+hstep) - f(x0-hstep))/(2*hstep)

    while j < n
        hstep = hstep / 2.0
        D[j+1, 1] = (f(x0+hstep) - f(x0-hstep)) / (2*hstep)
        for k = 1:j
            D[j+1,k+1] = D[j+1,k] + (D[j+1,k] - D[j,k]) / (4^k - 1)
        end

        err_new = 2 * abs(D[j+1,j+1] - D[j,j]) /
                  (abs(D[j+1,j+1]) + abs(D[j,j]) + epsilon)

        if err_new >= err; break; end
        err = err_new
        j += 1
    end

    return D[j,j]
end

#-- Complex-step derivative approach -----------------------------------------
#
function complex_step(f::Function, x0::Real; h::Real = 1e-20)
    # Check whether function accepts and returns complex values
    if !isa(f(x0), Real)
        error("Function 'f' must take on a real(!) value at 'x0'.")
    end

    fx0hi = try
        f(x0 + h*1im)
    catch err
        error("Function 'f' does not appear to accept complex arguments.")
    end

    # f(x0) must be real, fx0hi complex
    if !isa(fx0hi, Complex)
        error("'f(x0 + h*1im)' must be complex for 'complex-step' to work.")
    end

    return imag(fx0hi) / h
end


text:
using ReverseDiff

####################################################################

function grad_benchmark_driver!(out, f, x)
    println("benchmarking ∇$(f)...")

    cfg = ReverseDiff.GradientConfig(x)
    tp = ReverseDiff.GradientTape(f, x)

    # warmup
    ReverseDiff.gradient!(out, f, x, cfg)
    ReverseDiff.gradient!(out, tp, x)
    ReverseDiff.forward_pass!(tp)
    ReverseDiff.reverse_pass!(tp)

    # actual
    gc()
    print("  gradient! (no precord): ")
    @time ReverseDiff.gradient!(out, f, x, cfg)
    gc()
    print("  gradient!    (precord): ")
    @time ReverseDiff.gradient!(out, tp, x)
    gc()
    print("  forward pass: ")
    @time ReverseDiff.forward_pass!(tp)
    gc()
    print("  reverse pass: ")
    @time ReverseDiff.reverse_pass!(tp)
    gc()

    if length(tp.tape) <= 10000
        @eval begin
            out, x = $out, $x
            ctp = ReverseDiff.compile($tp)

            # warmup
            ReverseDiff.gradient!(out, ctp, x)
            ReverseDiff.forward_pass!(ctp)
            ReverseDiff.reverse_pass!(ctp)

            # actual
            print("  gradient! (compiled): ")
            @time ReverseDiff.gradient!(out, ctp, x)
            gc()
            print("  forward pass (compiled): ")
            @time ReverseDiff.forward_pass!(ctp)
            gc()
            print("  reverse pass (compiled): ")
            @time ReverseDiff.reverse_pass!(ctp)
            gc()
        end
    else
        println("skipping compiled GradientTape benchmark because the tape is too long ($(length(tp.tape)) elements)")
    end
end

####################################################################

rosenbrock(x) = sum(map(ReverseDiff.@forward((i, j) -> (1 - j)^2 + 100*(i - j^2)^2), x[2:end], x[1:end-1]))

# function rosenbrock(x)
#     i = x[2:end]
#     j = x[1:end-1]
#     return sum((1 .- j).^2 + 100*(i - j.^2).^2)
# end
#
# function rosenbrock(x::AbstractVector)
#     a = one(eltype(x))
#     b = 100 * a
#     result = zero(eltype(x))
#     for i in 1:length(x)-1
#         result += (a - x[i])^2 + b*(x[i+1] - x[i]^2)^2
#     end
#     return result
# end

x = rand(100000)
out = zeros(x)
grad_benchmark_driver!(out, rosenbrock, x)

####################################################################

function ackley(x::AbstractVector)
    a, b, c = 20.0, -0.2, 2.0*π
    len_recip = inv(length(x))
    sum_sqrs = zero(eltype(x))
    sum_cos = sum_sqrs
    for i in x
        sum_cos += cos(c*i)
        sum_sqrs += i^2
    end
    return (-a * exp(b * sqrt(len_recip*sum_sqrs)) -
            exp(len_recip*sum_cos) + a + e)
end

x = rand(100000)
out = zeros(x)
grad_benchmark_driver!(out, ackley, x)

####################################################################

function generate_matrix_test(n)
    return x -> begin
        @assert length(x) == 2n^2 + n
        a = reshape(x[1:n^2], n, n)
        b = reshape(x[n^2 + 1:2n^2], n, n)
        return trace(log.((a * b) + a - b))
    end
end

n = 100
matrix_test = generate_matrix_test(n)
x = collect(1.0:(2n^2 + n))
out = zeros(x)
grad_benchmark_driver!(out, matrix_test, x)

####################################################################

relu(x) = log.(1.0 .+ exp.(x))

ReverseDiff.@forward sigmoid(n) = 1. / (1. + exp(-n))

function neural_net(w1, w2, w3, x1)
    x2 = relu(w1 * x1)
    x3 = relu(w2 * x2)
    return sigmoid(dot(w3, x3))
end

xs = (randn(10,10), randn(10,10), randn(10), rand(10))
outs = map(similar, xs)
grad_benchmark_driver!(outs, neural_net, xs)

####################################################################

println("done")


text:
using JuliaDB

using Dates

function search_region_with_julia(vcf_dagger_file_path, chromosome, start_position::Int64, end_position::Int64)

    start_time = now()

    if isfile(vcf_dagger_file_path) == false

        println("Dagger file doesnt exist, fist make_vcf_indexedtable.")

    elseif isfile(vcf_dagger_file_path) == true

        vcf_table = load(vcf_dagger_file_path)

        println("Dagger file loaded.")

    end

    vcf_table_filtered = filter_vcf_table(vcf_table, chromosome, start_position, end_position)

    genotypes = get_genotype(vcf_table_filtered)

    vcf_table_filtered_gt = transform(vcf_table_filtered, :GT => genotypes)

    cols = columns(vcf_table_filtered_gt)

    vcf_table_filtered_gt_cut = table(cols[1], cols[2], cols[4], cols[5], cols[11], names = [:CHROM, :POS, :REF, :ALT, :GT])

    return vcf_table_filtered_gt_cut

    end_time = now()

    println("\nTook $(canonicalize(Dates.CompoundPeriod(end_time - start_time))).\n")

end

export search_region_with_julia


text:
using AudioSchedules: AudioSchedules
using Documenter: doctest

doctest(AudioSchedules)


text:
using ArbNumerics
using Test

setprecision(BigFloat, 512)

include("arithmetic.jl")