text:
@time @testset "splinevalue" begin

    # b1 = BSplineBasis(1, 0:5)
    s1 = Spline(b1, 1.0:5.0)
    # b2 = BSplineBasis(2, big(0):1:big(5))
    s2 = Spline(b2, 6.0:-1:1.0)
    # b3 = BSplineBasis(3, 0//1:5//1)
    s3 = Spline(b3, -30:10:30)
    # b4 = BSplineBasis(4, Int128[0,1,2,3,4,5])
    s4 = Spline(b4, 1:8)
    # b5 = BSplineBasis(3, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
    s5 = Spline(b5, ones(Float64, 7))
    # b6 = BSplineBasis(4, Rational{Int64}[-5//1, -7//2, -11//5, 1//1, 2//1, 18//5])
    s6 = Spline(b6, [1//i for i=eachindex(b6)])
    # b7 = BSplineBasis(3, Int8[-6,-4,-2,0,0,1,2])
    s7 = Spline(b7, view(1:10, 2:9))

    @testset "Return types" begin
        function test_splinevalue_returntype(spline, xtype, bspline_type, leftknot=intervalindex(basis(spline), 1))
            x = one(xtype)
            @test @inferred(splinevalue(spline, x, leftknot=leftknot)) isa bspline_type
            @test @inferred(splinevalue(spline, x, Derivative(1), leftknot=leftknot)) isa bspline_type
        end
        function test_splinevalue_returntype_workspace(spline, leftknot=intervalindex(basis(spline), 1))
            for (T,N) = ((Float32,0), (Float64,1), (BigFloat,2), (Rational{Int64},3), (Rational{BigInt},4))
                work = Vector{T}(undef, order(spline))
                @test @inferred(splinevalue(spline, 1, Derivative(N), leftknot=leftknot, workspace=work)) isa T
            end
        end

        # b1 = BSplineBasis(1, 0:5)
        b1_bspline_types = [(Int,Float64), (BigInt,BigFloat),
                            (Float32,Float32), (Float64,Float64), (BigFloat,BigFloat),
                            (Rational{Int},Rational{Int}), (Rational{BigInt},Rational{BigInt})]
        # b2 = BSplineBasis(2, big(0):1:big(5))
        b2_bspline_types = [(Int,BigFloat), (Float32,BigFloat), (Rational{Int},Rational{BigInt})]
        # b3 = BSplineBasis(3, 0//1:5//1)
        b3_bspline_types = [(BigInt,Rational{BigInt}), (Float64,Float64), (Rational{BigInt},Rational{BigInt})]
        # b4 = BSplineBasis(4, Int128[0,1,2,3,4,5])
        b4_bspline_types = [(Int,Float64), (BigFloat,BigFloat), (Rational{Int},Rational{Int128})]
        # b5 = BSplineBasis(3, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
        b5_bspline_types = [(BigInt,BigFloat), (Float32,Float64), (Rational{BigInt},BigFloat)]
        # b6 = BSplineBasis(4, Rational{Int64}[-5//1, -7//2, -11//5, 1//1, 2//1, 18//5])
        b6_bspline_types = [(Int,Rational{Int64}), (Float64,Float64), (Rational{Int},Rational{Int64})]
        # b7 = BSplineBasis(3, Int8[-6,-4,-2,0,0,1,2])
        b7_bspline_types = [(BigInt,BigFloat), (BigFloat,BigFloat), (Rational{Int},Rational{Int})]
        for (basis, types) = [(b1, b1_bspline_types), (b2, b2_bspline_types),
                              (b3, b3_bspline_types), (b4, b4_bspline_types),
                              (b5, b5_bspline_types), (b6, b6_bspline_types),
                              (b7, b7_bspline_types)]
            bspline = basis[1]
            leftknot = intervalindex(basis, 1)
            test_splinevalue_returntype_workspace(bspline, leftknot)
            for (xtype, bspline_type) in types
                test_splinevalue_returntype(bspline, xtype, bspline_type, leftknot)
            end
        end

        # b1 = BSplineBasis(1, 0:5)
        b1_coeffs_and_types = [(1:5,Int,Float64),
                               (1:5,Float32,Float32)]
        # b2 = BSplineBasis(2, big(0):1:big(5))
        b2_coeffs_and_types = [(1//1:6//1,Int8,Rational{BigInt}),
                               (1:6,Int,BigFloat)]
        # b3 = BSplineBasis(3, 0//1:5//1)
        b3_coeffs_and_types = [(1:7,Float32,Float32),
                               (Float32.(1:7),Float64,Float64)]
        # b4 = BSplineBasis(4, Int128[0,1,2,3,4,5])
        b4_coeffs_and_types = [(Float64.(1:8),Rational{Int},Float64),
                               (1:8,Rational{Int},Rational{Int128})]
        # b5 = BSplineBasis(3, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
        b5_coeffs_and_types = [(Float32.(1:7),Float32,Float64),
                               (big(1):big(7),Float64,BigFloat)]
        # b6 = BSplineBasis(4, Rational{Int64}[-5//1, -7//2, -11//5, 1//1, 2//1, 18//5])
        b6_coeffs_and_types = [(1:8,Int,Rational{Int64}),
                               (1:8,Float64,Float64)]
        # b7 = BSplineBasis(3, Int8[-6,-4,-2,0,0,1,2])
        b7_coeffs_and_types = [(Float64.(1:8),Rational{Int},Float64),
                               (BigFloat.(1:8),Int,BigFloat)]
        for (basis, coeffs_and_types) = [(b1, b1_coeffs_and_types), (b2, b2_coeffs_and_types),
                                         (b3, b3_coeffs_and_types), (b4, b4_coeffs_and_types),
                                         (b5, b5_coeffs_and_types), (b6, b6_coeffs_and_types),
                                         (b7, b7_coeffs_and_types)]
            leftknot = intervalindex(basis, 1)
            for (coeffs, xtype, bspline_type) in coeffs_and_types
                spline = Spline(basis, coeffs)
                test_splinevalue_returntype(spline, xtype, bspline_type, leftknot)
            end
        end
        for spline = [s1, s2, s3, s4, s5, s6, s7]
            test_splinevalue_returntype_workspace(spline)
        end
    end

    function test_splinevalue_zero(spline::Spline, x)
        work = zeros(order(spline))
        for N = 0:3
            @test splinevalue(spline, x, Derivative(N)) == 0
            @test splinevalue(spline, x, Derivative(N), workspace=work) == 0
            @test splinevalue(spline, x, Derivative(N), leftknot=nothing) == 0
            @test splinevalue(spline, x, Derivative(N), leftknot=nothing, workspace=work) == 0
            @test_throws ArgumentError splinevalue(spline, x, Derivative(N), leftknot=1)
            @test_throws ArgumentError splinevalue(spline, x, Derivative(N), leftknot=1, workspace=work)
        end
    end
    test_splinevalue_zero(b1[1], Inf16)
    test_splinevalue_zero(b2[2], -Inf32)
    test_splinevalue_zero(b3[3], Inf64)
    test_splinevalue_zero(b4[4], 1//0)
    test_splinevalue_zero(b5[5], -1//0)
    test_splinevalue_zero(b6[6], -10.0)
    test_splinevalue_zero(b7[7], 100)

    function test_splinevalue_nan(spline::Spline, x)
        work = zeros(order(spline))
        for N = 0:3
            @test isnan(splinevalue(spline, x, Derivative(N)))
            @test isnan(splinevalue(spline, x, Derivative(N), workspace=work))
        end
        if isnan(x)
            for N = 0:3
                @test isnan(splinevalue(spline, x, Derivative(N), leftknot=nothing))
                @test isnan(splinevalue(spline, x, Derivative(N), leftknot=nothing, workspace=work))
                @test_throws ArgumentError splinevalue(spline, x, Derivative(N), leftknot=1)
                @test_throws ArgumentError splinevalue(spline, x, Derivative(N), leftknot=1, workspace=work)
            end
        end
    end
    test_splinevalue_nan(b1[1], NaN16)
    test_splinevalue_nan(b2[2], NaN32)
    test_splinevalue_nan(b3[3], NaN64)
    test_splinevalue_nan(b4[4], big(NaN))

    function splinevalue_exact(splines, coeffs, x, nderiv=0)
        if nderiv ≥ size(splines,2)
            0
        elseif coeffs isa BSplines.StandardBasisVector
            (splines[coeffs.index,nderiv+1])(maybebig(x))
        else
            bigx = maybebig(x)
            sum(coeffs[i]*(splines[i,nderiv+1])(bigx) for i=eachindex(coeffs))
        end
    end
    # Enlarge x for calculating exact B-splines:
    # * for Integer/Rational: use 64-bit to avoid overflow in rational arithmetic
    # * for floating-point: use BigFloat to minimize error
    maybebig(x::Integer) = Int64(x)
    maybebig(x::Rational) = Rational{Int64}(x)
    maybebig(x::AbstractFloat) = big(x)

    function test_splinevalue_nonzero(spline::Spline, basis_exact, x, isexact)
        leftknot = intervalindex(basis(spline), x)
        leftknot === nothing && error("this function is only suitable for x which are inside the support of the basis")
        work = Array{isexact ? Rational{Int64} : Float64}(undef, order(spline))
        for N = 0:3
            spl = splinevalue(spline, x, Derivative(N), leftknot=leftknot)
            exact = splinevalue_exact(basis_exact, coeffs(spline), x, N)
            if isexact
                @test splinevalue(spline, x, Derivative(N), leftknot=leftknot) == exact
                @test splinevalue(spline, x, Derivative(N), leftknot=leftknot, workspace=work) == exact
            else
                @test splinevalue(spline, x, Derivative(N), leftknot=leftknot) ≈ₑₗ exact
                @test splinevalue(spline, x, Derivative(N), leftknot=leftknot, workspace=work) ≈ₑₗ exact
            end
        end
    end

    function test_splinevalue_bsplines(basis::BSplineBasis, basis_exact, xtype, isexact)
        for bspline = basis
            a, b = support(bspline)
            for x = xrange(xtype, a, b, 5)
                test_splinevalue_nonzero(bspline, basis_exact, x, isexact)
            end
            a′, b′ = support(basis)
            for x = xrange(xtype, a′, b′, 5)
                test_splinevalue_nonzero(bspline, basis_exact, x, isexact)
            end
        end
    end

    function xrange(xtype, a, b, length)
        if xtype <: Integer
            xs = ceil(xtype, a):floor(xtype, b)
        elseif xtype <: AbstractFloat
            xa = xtype(a)
            xb = xtype(b)
            xa < a && (xa = nextfloat(xa))
            xb > b && (xb = prevfloat(xb))
            xs = range(xa, stop=xb, length=length)
        else
            xs = range(xtype(a), stop=xtype(b), length=length)
        end
        xs
    end

    # b1 = BSplineBasis(1, 0:5)
    test_splinevalue_bsplines(b1, b1_exact, BigInt,           true)
    test_splinevalue_bsplines(b1, b1_exact, Float32,          true)
    # b2 = BSplineBasis(2, big(0):1:big(5))
    test_splinevalue_bsplines(b2, b2_exact, BigFloat,         false)
    test_splinevalue_bsplines(b2, b2_exact, Rational{Int32},  true)
    # b3 = BSplineBasis(3, 0//1:5//1)
    test_splinevalue_bsplines(b3, b3_exact, Float64,          false)
    test_splinevalue_bsplines(b3, b3_exact, Rational{BigInt}, true)
    # b4 = BSplineBasis(4, Int128[0,1,2,3,4,5])
    test_splinevalue_bsplines(b4, b4_exact, BigInt,           false)
    test_splinevalue_bsplines(b4, b4_exact, Rational{Int64},  true)
    # b5 = BSplineBasis(3, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
    test_splinevalue_bsplines(b5, b5_exact, Int,              false)
    test_splinevalue_bsplines(b5, b5_exact, Float32,          false)
    # b6 = BSplineBasis(4, Rational{Int64}[-5//1, -7//2, -11//5, 1//1, 2//1, 18//5])
    test_splinevalue_bsplines(b6, b6_exact, Int,              true)
    test_splinevalue_bsplines(b6, b6_exact, Float64,          false)
    # b7 = BSplineBasis(3, Int8[-6,-4,-2,0,0,1,2])
    test_splinevalue_bsplines(b7, b7_exact, BigFloat,         false)
    test_splinevalue_bsplines(b7, b7_exact, Rational{Int16},  true)

    function test_splinevalue(spline::Spline, basis_exact, xtype, isexact)
        a, b = support(spline)
        for x = xrange(xtype, a, b, 10)
            test_splinevalue_nonzero(spline, basis_exact, x, isexact)
        end
    end

    # b1 = BSplineBasis(1, 0:5)
    # s1 = Spline(b1, 1.0:5.0)
    test_splinevalue(s1, b1_exact, BigInt,           true)
    test_splinevalue(s1, b1_exact, Float64,          true)
    # b2 = BSplineBasis(2, big(0):1:big(5))
    # s2 = Spline(b2, 6.0:-1:1.0)
    test_splinevalue(s2, b2_exact, BigFloat,         false)
    test_splinevalue(s2, b2_exact, Rational{Int32},  false)
    # b3 = BSplineBasis(3, 0//1:5//1)
    # s3 = Spline(b3, -30:10:30)
    test_splinevalue(s3, b3_exact, Float32,          false)
    test_splinevalue(s3, b3_exact, Rational{BigInt}, true)
    # b4 = BSplineBasis(4, Int128[0,1,2,3,4,5])
    # s4 = Spline(b4, 1:8)
    test_splinevalue(s4, b4_exact, BigInt,           false)
    test_splinevalue(s4, b4_exact, Rational{Int64},  true)
    # b5 = BSplineBasis(3, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
    # s5 = Spline(b5, ones(Float64, 7))
    test_splinevalue(s5, b5_exact, Int,              false)
    test_splinevalue(s5, b5_exact, Float32,          false)
    # b6 = BSplineBasis(4, Rational{Int64}[-5//1, -7//2, -11//5, 1//1, 2//1, 18//5])
    # s6 = Spline(b6, [1//i for i=eachindex(b6)])
    test_splinevalue(s6, b6_exact, Int,              true)
    test_splinevalue(s6, b6_exact, Float64,          false)
    # b7 = BSplineBasis(3, Int8[-6,-4,-2,0,0,1,2])
    # s7 = Spline(b7, view(1:10, 2:9))
    test_splinevalue(s7, b7_exact, BigFloat,         false)
    test_splinevalue(s7, b7_exact, Rational{Int16},  true)

    test_splinevalue_zero(s1, Inf16)
    test_splinevalue_zero(s2, -Inf32)
    test_splinevalue_zero(s3, Inf64)
    test_splinevalue_zero(s4, 1//0)
    test_splinevalue_zero(s5, -1//0)
    test_splinevalue_zero(s6, -10.0)
    test_splinevalue_zero(s7, 100)
    test_splinevalue_nan(s1, NaN16)
    test_splinevalue_nan(s2, NaN32)
    test_splinevalue_nan(s3, NaN64)
    test_splinevalue_nan(s4, big(NaN))
    test_splinevalue_nan(Spline(b4, fill(NaN, length(b4))), 1.5)

    k = order(s5);
    @test_throws DimensionMismatch splinevalue(s5, 0, workspace=zeros(k+1))
    @test_throws DimensionMismatch splinevalue(s5, 0, workspace=zeros(k-1))
    @test_throws DimensionMismatch splinevalue(s5, 0, workspace=zeros(1, k))
    @test_throws DimensionMismatch splinevalue(s5, 0, workspace=zeros(k, 1))
    @test_throws DimensionMismatch splinevalue(s5, 0, workspace=zeros(k, k))
    leftknot = intervalindex(basis(s5), 0)
    @test_throws ArgumentError splinevalue(s5, 0, leftknot=leftknot-1)
    @test_throws ArgumentError splinevalue(s5, 0, leftknot=leftknot+1)
    @test_throws ArgumentError splinevalue(s5, 0, leftknot=nothing)
    @test_throws ArgumentError splinevalue(s5, 0, Derivative(0), leftknot=leftknot-1)
    @test_throws ArgumentError splinevalue(s5, 0, Derivative(1), leftknot=leftknot+1)
    @test_throws ArgumentError splinevalue(s5, 0, Derivative(2), leftknot=nothing)
    @test s5(0) == s5(0, leftknot=leftknot) == splinevalue(s5, 0, leftknot=leftknot)
    @test s5(0, Derivative(1)) == s5(0, Derivative(1), leftknot=leftknot) == splinevalue(s5, 0, Derivative(1), leftknot=leftknot)
end


text:
module StartUpDG

using ConstructionBase: ConstructionBase
using LinearAlgebra: diagm, eigvals, Diagonal, I
using StaticArrays: SVector, SMatrix
using Setfield: setproperties, @set    # for "modifying" structs (setproperties)
using Kronecker: kronecker # for Hex element matrix manipulations

using Reexport
@reexport using UnPack            # for getting values in RefElemData and MeshData
@reexport using NodesAndModes     # for basis functions

using NodesAndModes: meshgrid
using SparseArrays: sparse, droptol!
using MAT # read in .mat files
using RecipesBase, Colors

# reference element utility functions
include("RefElemData.jl")
include("RefElemData_polynomial.jl")
include("RefElemData_SBP.jl")
include("ref_elem_utils.jl")
export RefElemData, Polynomial
export SBP, DefaultSBPType, TensorProductLobatto, Hicken, Kubatko  # types for SBP node dispatch
export LobattoFaceNodes, LegendreFaceNodes # type parameters for SBP{Kubatko{...}}
export hybridized_SBP_operators, inverse_trace_constant, face_type

include("MeshData.jl")
export MeshData

include("geometric_functions.jl")
export geometric_factors, estimate_h

# spatial connectivity routines
include("connectivity_functions.jl")
export make_periodic

# for tagging faces on boundaries
include("boundary_utils.jl")
export boundary_face_centroids, tag_boundary_faces

# uniform meshes + face vertex orderings are included
include("mesh/simple_meshes.jl")
export readGmsh2D, uniform_mesh

# Plots.jl recipes for meshes
include("mesh/mesh_visualization.jl")
export VertexMeshPlotter, MeshPlotter

using Requires
function __init__()
    @require Triangulate="f7e6ffb2-c36d-4f8f-a77e-16e897189344" begin
        using Printf
        using .Triangulate: TriangulateIO, triangulate
        include("mesh/triangulate_utils.jl")
        export refine, triangulateIO_to_VXYEToV, get_node_boundary_tags
        export BoundaryTagPlotter
        include("mesh/triangulate_example_meshes.jl")
        export triangulate_domain
        export Scramjet, SquareDomain, RectangularDomain, RectangularDomainWithHole
    end
end

# simple explicit time-stepping included for conveniencea
export ck45 # LSERK 45
include("explicit_timestep_utils.jl")

end

ext: jl


text:
using DataFrames

mutable struct DFBindInfo
    df::DataFrame
    result_types::Vector{Type}
    scan_functions::Vector{Function}

    function DFBindInfo(df::DataFrame, result_types::Vector{Type}, scan_functions::Vector{Function})
        return new(df, result_types, scan_functions)
    end
end

function DFResultType(df, entry)
    column_type = eltype(df[!, entry])
    if typeof(column_type) == Union
        # remove Missing type from the union
        column_type = Core.Compiler.typesubtract(column_type, Missing, 1)
    end
    return column_type
end

function DFJuliaType(column_type)
    if column_type == Date
        column_type = Int32
    elseif column_type == Time
        column_type = Int64
    elseif column_type == DateTime
        column_type = Int64
    end
    return column_type
end

ValueToDuckDB(val::T) where {T <: Date} = convert(Int32, Dates.date2epochdays(val) - ROUNDING_EPOCH_TO_UNIX_EPOCH_DAYS)
ValueToDuckDB(val::T) where {T <: Time} = convert(Int64, Dates.value(val) / 1000)
ValueToDuckDB(val::T) where {T <: DateTime} =
    convert(Int64, (Dates.datetime2epochms(val) - ROUNDING_EPOCH_TO_UNIX_EPOCH_MS) * 1000)
function ValueToDuckDB(val::T) where {T <: AbstractString}
    throw(
        NotImplementedException(
            "Cannot use ValueToDuckDB to convert string values - use DuckDB.AssignStringElement on a vector instead"
        )
    )
end
function ValueToDuckDB(val::T) where {T}
    return val
end

function DFScanColumn(
    df::DataFrame,
    df_offset::Int64,
    col_idx::Int64,
    scan_count::Int64,
    output::DuckDB.DataChunk,
    ::Type{T}
) where {T}
    vector = DuckDB.GetVector(output, col_idx)
    result_array::Vector{T} = DuckDB.GetArray(vector, T)
    validity = DuckDB.GetValidity(vector)
    input_column = df[!, col_idx]
    for i in 1:scan_count
        if input_column[df_offset + i] === missing
            DuckDB.SetInvalid(validity, i)
        else
            result_array[i] = ValueToDuckDB(input_column[df_offset + i])
        end
    end
end

function DFScanStringColumn(
    df::DataFrame,
    df_offset::Int64,
    col_idx::Int64,
    scan_count::Int64,
    output::DuckDB.DataChunk,
    ::Type{T}
) where {T}
    vector = DuckDB.GetVector(output, col_idx)
    validity = DuckDB.GetValidity(vector)
    input_column = df[!, col_idx]
    for i in 1:scan_count
        if input_column[df_offset + i] === missing
            DuckDB.SetInvalid(validity, i)
        else
            DuckDB.AssignStringElement(vector, i, input_column[df_offset + i])
        end
    end
end

function DFScanFunction(df, entry)
    result_type = DFResultType(df, entry)
    if result_type <: AbstractString
        return DFScanStringColumn
    end
    return DFScanColumn
end

function DFBindFunction(info::DuckDB.BindInfo)
    # fetch the df name from the function parameters
    parameter = DuckDB.GetParameter(info, 0)
    name = DuckDB.GetValue(parameter, String)
    # fetch the actual df using the function name
    extra_data = DuckDB.GetExtraData(info)
    df = extra_data[name]

    # register the result columns
    result_types::Vector{Type} = Vector()
    scan_functions::Vector{Function} = Vector()
    for entry in names(df)
        result_type = DFResultType(df, entry)
        scan_function = DFScanFunction(df, entry)
        push!(result_types, DFJuliaType(result_type))
        push!(scan_functions, scan_function)

        DuckDB.AddResultColumn(info, entry, result_type)
    end
    return DFBindInfo(df, result_types, scan_functions)
end

mutable struct DFInitInfo
    pos::Int64

    function DFInitInfo()
        return new(0)
    end
end

function DFInitFunction(info::DuckDB.InitInfo)
    return DFInitInfo()
end

function DFScanFunction(info::DuckDB.FunctionInfo, output::DuckDB.DataChunk)
    bind_info = DuckDB.GetBindInfo(info, DFBindInfo)
    init_info = DuckDB.GetInitInfo(info, DFInitInfo)

    column_count = size(names(bind_info.df), 1)
    row_count = size(bind_info.df, 1)
    scan_count = DuckDB.VECTOR_SIZE
    if init_info.pos + scan_count >= row_count
        scan_count = row_count - init_info.pos
    end

    for col_idx in 1:column_count
        bind_info.scan_functions[col_idx](
            bind_info.df,
            init_info.pos,
            col_idx,
            scan_count,
            output,
            bind_info.result_types[col_idx]
        )
    end
    init_info.pos += scan_count
    DuckDB.SetSize(output, scan_count)
    return
end

function RegisterDataFrame(con::Connection, df::DataFrame, name::AbstractString)
    con.db.data_frames[name] = df
    DBInterface.execute(
        con,
        string("CREATE OR REPLACE VIEW \"", name, "\" AS SELECT * FROM julia_df_scan('", name, "')")
    )
    return
end
RegisterDataFrame(db::DB, df::DataFrame, name::AbstractString) = RegisterDataFrame(db.main_connection, df, name)

function UnregisterDataFrame(con::Connection, name::AbstractString)
    pop!(con.db.data_frames, name)
    DBInterface.execute(con, string("DROP VIEW IF EXISTS \"", name, "\""))
    return
end
UnregisterDataFrame(db::DB, name::AbstractString) = UnregisterDataFrame(db.main_connection, name)

function AddDataFrameScan(db::DB)
    # add the data frame scan function
    DuckDB.CreateTableFunction(
        db.main_connection,
        "julia_df_scan",
        [String],
        DFBindFunction,
        DFInitFunction,
        DFScanFunction,
        db.handle.data_frames
    )
    return
end


text:
using Salsa
using BenchmarkTools

@declare_input in1(rt, ::Int)::Int

@derived function d1(rt, x::Int)
    return in1(rt, x) + 1
end

@derived function d2(rt, x::Int)
    return d1(rt, x) + 1
end

function simple_bench(rt::Runtime; N, iters)
    for j in 1:iters
        Salsa.new_epoch!(rt)
        for i in 1:N
            set_in1!(rt, i, i+j)
        end
        for i in 1:N
            d2(rt, i)
        end
    end
end

function parallel_bench(rt::Runtime; N, iters)
    for j in 1:iters
        Salsa.new_epoch!(rt)
        for i in 1:N
            set_in1!(rt, i, i+j)
        end
        @sync for i in 1:N
            Threads.@spawn begin
                d2(rt, i)
            end
        end
    end
end

# Run all the above benchmarks and report the results to stdout
function run(; N, iters)
    @btime simple_bench(rt; N=$N, iters=$iters) setup=(rt=Runtime())
    @btime parallel_bench(rt; N=$N, iters=$iters) setup=(rt=Runtime())
end


text:
const Tsent = Vector{Tuple{Ti, Ti}}

struct ZhuXian
        words::Vector{String}
        postags::Vector{String}
        sents::Vector{Tsent}
end

function ZhuXian(path)
        wmp, pmp = Dict{String, Ti}(), Dict{String, Ti}()
        function go(file)
                sents = Tsent[]
                for line in eachline(file)
                        sent = Tsent()
                        for wordpos in split(line)
                                word, pos = String.(split(wordpos, '_'))
                                iw = getid!(wmp, word)
                                ip = getid!(pmp, pos)
                                push!(sent, (ip, iw))
                        end
                        push!(sents, sent)
                end
                sents
        end
        sents = Tsent[]
        append!(sents, go(joinpath(path, "train.zhuxian.wordpos")))
        append!(sents, go(joinpath(path, "dev.zhuxian.wordpos")))
        append!(sents, go(joinpath(path, "test.zhuxian.wordpos")))

        ZhuXian(vec(wmp), vec(pmp), sents)
end

function ZhuXian()
        file = KongYiji.dir("zhuxian.jld2")
        zfile = KongYiji.dir("zhuxian.jld2.7z")
        if isfile(file)
                r = load(file)["zx"]
        elseif isfile(zfile)
                r = load(KongYiji.decompress(zfile))["zx"]
        else
                error("zhuxian not downloaded.")
        end
        return r
end

import Base.length
length(zx::ZhuXian) = length(zx.sents)

sents(zx::ZhuXian) = zx.sents
words(zx::ZhuXian) = zx.words
postags(zx::ZhuXian) = zx.postags

function rawsents(zx::ZhuXian)
        r = [join([zx.words[p[2]] for p in sent]) for sent in sents(zx)]
        r
end

function wordsents(zx::ZhuXian)
        r = [[zx.words[p[2]] for p in sent] for sent in sents(zx)]
        r
end

function posidsents(zx::ZhuXian)
        r = [[p[1] for p in sent] for sent in sents(zx)]
        r
end

function foldbatch(zx::ZhuXian, nbatch::Int)
        sents = collect(KongYiji.sents(zx)) |> shuffle!
        ns = length(sents)
        nb = min(ns, nbatch)
        @assert 2 <= nb
        (begin
                te = Tsent[]
                tr = Tsent[]
                npick = div(ns + nb - 1, nb)
                from = min(ns + 1, (ib - 1) * npick + 1)
                to = min(ns, from + npick - 1)
                append!(te, sents[from:to])
                append!(tr, sents[1:from - 1])
                append!(tr, sents[to + 1:ns])
                te = ZhuXian(zx.words, zx.postags, te)
                tr = ZhuXian(zx.words, zx.postags, tr)
                (tr, te)
        end for ib in 1:nb)
end

import Base.show
function show(io::IO, zx::ZhuXian)
        ns = length(zx)
        nw = mapreduce(length, +, sents(zx))
        print(io, "Zhu Xian of $(ns) sents, $(length(zx.words))/$(nw) words")
end

function ==(a::ZhuXian, b::ZhuXian)
        return all(fname -> getfield(a, fname) == getfield(b, fname), fieldnames(ZhuXian))
end

function mask(zx::ZhuXian)
        words = mask.(zx.words) |> sort! |> unique!
        wmp = dict(words)
        ns = length(zx)
        sents = Tsent[]
        for sent in KongYiji.sents(zx)
                sent2 = Tsent()
                for (ip, iw) in sent
                        iw = wmp[mask[zx.words[iw]]]
                        push!(sent2, (ip, iw))
                end
                push!(sents, sent2)
        end
        ZhuXian(words, zx.postags, sents)
end


text:
# Autogenerated wrapper script for Raylib_jll for aarch64-apple-darwin
export libraylib

using GLFW_jll
JLLWrappers.@generate_wrapper_header("Raylib")
JLLWrappers.@declare_library_product(libraylib, "@rpath/libraylib.400.dylib")
function __init__()
    JLLWrappers.@generate_init_header(GLFW_jll)
    JLLWrappers.@init_library_product(
        libraylib,
        "lib/libraylib.4.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct StorageAccountKey <: SwaggerModel
    keyName::Any # spec type: Union{ Nothing, String } # spec name: keyName
    value::Any # spec type: Union{ Nothing, String } # spec name: value
    permissions::Any # spec type: Union{ Nothing, String } # spec name: permissions

    function StorageAccountKey(;keyName=nothing, value=nothing, permissions=nothing)
        o = new()
        validate_property(StorageAccountKey, Symbol("keyName"), keyName)
        setfield!(o, Symbol("keyName"), keyName)
        validate_property(StorageAccountKey, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(StorageAccountKey, Symbol("permissions"), permissions)
        setfield!(o, Symbol("permissions"), permissions)
        o
    end
end # type StorageAccountKey

const _property_map_StorageAccountKey = Dict{Symbol,Symbol}(Symbol("keyName")=>Symbol("keyName"), Symbol("value")=>Symbol("value"), Symbol("permissions")=>Symbol("permissions"))
const _property_types_StorageAccountKey = Dict{Symbol,String}(Symbol("keyName")=>"String", Symbol("value")=>"String", Symbol("permissions")=>"String")
Base.propertynames(::Type{ StorageAccountKey }) = collect(keys(_property_map_StorageAccountKey))
Swagger.property_type(::Type{ StorageAccountKey }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_StorageAccountKey[name]))}
Swagger.field_name(::Type{ StorageAccountKey }, property_name::Symbol) =  _property_map_StorageAccountKey[property_name]

const _allowed_StorageAccountKey_permissions = ["Read", "Full"]

function check_required(o::StorageAccountKey)
    true
end

function validate_property(::Type{ StorageAccountKey }, name::Symbol, val)
    if name === Symbol("permissions")
        Swagger.validate_param(name, "StorageAccountKey", :enum, val, _allowed_StorageAccountKey_permissions)
    end
end


text:
@testset "cuda/msgpass" begin
    T = Float32
    in_channel = 10
    out_channel = 5
    N = 6

    adj = [0. 1. 0. 0. 0. 0.;
        1. 0. 0. 1. 1. 1.;
        0. 0. 0. 0. 0. 1.;
        0. 1. 0. 0. 1. 0.;
        0. 1. 0. 1. 0. 1.;
        0. 1. 1. 0. 1. 0.]

    struct NewCudaLayer{T} <: MessagePassing
        weight::T
    end
    NewCudaLayer(m, n) = NewCudaLayer(randn(T, m,n))
    @functor NewCudaLayer

    function (l::NewCudaLayer)(fg::FeaturedGraph, X::AbstractMatrix)
        _, x, _ = GeometricFlux.propagate(l, fg, edge_feature(fg), X, global_feature(fg), +)
        x
    end

    (l::NewCudaLayer)(fg::FeaturedGraph) = FeaturedGraph(fg, nf = l(fg, node_feature(fg)))

    GeometricFlux.message(n::NewCudaLayer, x_i, x_j, e_ij) = n.weight * x_j
    GeometricFlux.update(::NewCudaLayer, m, x) = m

    X = rand(T, in_channel, N)
    fg = FeaturedGraph(adj, nf=X) |> gpu
    l = NewCudaLayer(out_channel, in_channel) |> gpu

    fg_ = l(fg)
    @test size(node_feature(fg_)) == (out_channel, N)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8699f9e2-20cc-11ec-30f2-89cda443cd19
using PlutoUI

# ╔═╡ f437f2d8-4106-4557-9199-75393e6b7a1d
md"""
# Slider reset method

This slider method was proposed by Gerhard Dorn in a response to my Zulip post of Sep 28 2021.

I had not seen this syntax before.
"""

# ╔═╡ 6def317b-85d5-423f-8d12-98c4a1d3ee14
button_widget = @bind button_variable Button("Reset")

# ╔═╡ 12df3e12-b727-426d-b0ef-ce0e2555d333
begin
	button_variable
	slider_widget = @bind slider_variable Slider(1:10)
end ;

# ╔═╡ 1554a4a8-3854-4118-9043-103c682ab9c0
begin
	md"""
	Use the slider $(slider_widget) or set it: $(button_widget)
	"""
end

# ╔═╡ e832539b-2b86-49fb-a074-80384dbeca92
slider_variable

# ╔═╡ 9219b1e7-0b54-41cc-a489-e93eca18c6be
slider_widget

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.11"
"""

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

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

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

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

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

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

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

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

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

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

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "0c3e067931708fa5641247affc1a1aceb53fff06"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.11"

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

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

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

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

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

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

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

# ╔═╡ Cell order:
# ╠═f437f2d8-4106-4557-9199-75393e6b7a1d
# ╠═8699f9e2-20cc-11ec-30f2-89cda443cd19
# ╠═1554a4a8-3854-4118-9043-103c682ab9c0
# ╠═6def317b-85d5-423f-8d12-98c4a1d3ee14
# ╠═12df3e12-b727-426d-b0ef-ce0e2555d333
# ╠═e832539b-2b86-49fb-a074-80384dbeca92
# ╠═9219b1e7-0b54-41cc-a489-e93eca18c6be
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# Show grid and random walk target trajectory

using SubmodularMaximization
using PyPlot

close("all")

steps = 100

grid = Grid(10, 10)

initial = random_state(grid)
states = Array{typeof(initial)}(undef, steps)
states[1] = initial

plot_state_space(grid)

for ii = 2:steps
  println("Step ", ii)
  states[ii] = target_dynamics(grid, states[ii-1])

  plots = plot_trajectory(states[1:ii])

  sleep(0.1)

  if ii < steps
    foreach(x->x.remove(), plots)
  end
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
#Run analysis with neutral data, input parameters:
# mu = 20.0
# cm = 200
# cellularity = 0.7
println()
println("####################")
println("Checking that inference on neutral simulated data with known input parameters returns ground truth...")
Random.seed!(123)
@time out = fitABCmodels("data/neutral.txt",
  "neutral",
  read_depth = 150,
  minvaf = 0.05,
  resultsdirectory = "output",
  nparticles = 100,
  maxiterations = 10^4,
  Nmax = 10^3,
  maxclones = 1,
  firstpass = false,
  progress = true,
  verbose = true,
  save = true);

# check if we get the correct model
println("")
println("\tTesting posterior model probability returns neutral as most probable model...")
@test argmax(out.ModelProb[:Probability]) == 1

# extract parameters
mu = out.Posterior[1].Parameters[:mu]
clonalmutations = out.Posterior[1].Parameters[:clonalmutations]
cellularity = out.Posterior[1].Parameters[:cellularity]

println("\tChecking true parameters are within the 95% credible interval range...")
# check if parameters are within 80% credible interval
@test quantile(mu, 0.025) < 20.0 < quantile(mu, 0.975)
@test quantile(clonalmutations, 0.025) < 200.0 < quantile(clonalmutations, 0.975)
@test quantile(cellularity, 0.025) < 0.7 < quantile(cellularity, 0.975)

println("\tChecking plotting functions work and are saved...")
saveallplots(out)
@test isfile("output/neutral/finalpopulation/plots/neutral-histogram-0clone.pdf")
@test isfile("output/neutral/finalpopulation/posterior/neutral-histogram-clone0.csv")
rm("output", recursive = true)

println("Check if saving intermediate outputs works correctly")
@time out = fitABCmodels("data/neutral.txt",
  "neutral",
  read_depth = 150,
  minvaf = 0.05,
  resultsdirectory = "output",
  nparticles = 100,
  maxiterations = 10^4,
  Nmax = 10^3,
  maxclones = 1,
  firstpass = false,
  progress = true,
  verbose = true,
  save = true,
  savepopulations = true);
  @test isfile("output/neutral/populations/population_1/posterior/neutral-histogram-clone0.csv")
  rm("output", recursive = true)

println("All tests passed for neutral simulated data.")


text:
module Production
##################Dependencies####################
using Revise

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

if "$(pwd())\\data" ∉ LOAD_PATH
  push!(LOAD_PATH,"$(pwd())\\data")
end

using Finometrics #my personal module
using DataFrames, Distributions, GZip, ForwardDiff, Dates, LinearAlgebra, Statistics,
   Formatting, StatsBase, DataStructures, Distributed,
  Serialization,  CSV, Random#, StatsPlots, UnicodePlots
  #, Measures, Gadfly, Fontconfig, Cairo NOTE: uncomment for graphs

############################Exported functions
export Firm, Problem

#############################Constants
const STATES_PER_VAR = 10

###Firm Specific DefaultParameters
const DEF_τ = [0.2, 0.3]
const DEF_δ = 0.3
const DEF_Θ = (0.8, 0.2)
const DEF_α = 1.05
const DEF_β = 0.98
const DEF_γ = Float64(1/3)
const DEF_ζ = (1.0, 0.0)
const DEF_KSTATES = collect(range(0.0, stop=5.0, length=STATES_PER_VAR))
const DEF_YSTATES = collect(range(0.0, stop=5.0, length=STATES_PER_VAR))

const TOL_CONSTRAINT = 1.0e-10
const N_TOL = 1.0e-10

#Other problem default parameters
const DEF_p = [0.81 0.68; 0.19 0.32]
const DEF_M = 3
const DEF_D_TOL = 10e-6
const DEF_MAX_ITER = 1000

const OUT_NAME = "test"
const OUT_PAPER = "C:\\Users\\Clinton Tepper\\Dropbox\\Apps\\Overleaf\\productionproject\\sub"
const OUT_PATH = "output"

const MAX_WORKERS = Base.Sys.CPU_THREADS ÷ 2 # gets num of physical cores

##########################code files##########################
include("code\\Numerical.jl")
include("code\\Firm.jl")
include("code\\Problem.jl")

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#read_pdf.jl
module Read_PDF
export read_pdf

using JSON

function to_tuple(X) :: Any
    if typeof(X) <: Array
        tmp = [ to_tuple(x) for x in X ]
        x = (tmp[1:end]...)
    else
        return X
    end
end

function array_to_dict(a) :: Dict{Tuple{Any,Any,Any}, Float64}
    local d = Dict{Tuple{Any,Any,Any}, Float64}()
    for p in a
        t = Tuple{Any,Any,Any}( to_tuple( p[1] ) )
        d[t] = Float64( p[2] )
    end
    return d
end

function read_pdf(filename::String) ::  Dict{Tuple{Any,Any,Any}, Float64}
    a = JSON.parsefile(filename)
    return array_to_dict(a)
end

end

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
using Klara

### Define the log-target as a function (generic or anonymous):

plogtarget(z::Vector{Float64}) = -dot(z, z)

### Define the parameter via BasicContMuvParameter (it is a continuous multivariate variable)
### The input arguments for BasicContMuvParameter are:
### 1) the variable key,
### 2) the log-target

p = BasicContMuvParameter(:p, logtarget=plogtarget)

### Define the model using the likelihood_model generator
### The second argument informs the likelihood_model generator that p.index has not been set

model = likelihood_model(p, false)

### Define a Metropolis-Hastings sampler with an identity covariance matrix

sampler = MH(ones(2))

### Set MCMC sampling range

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

### Set initial values for simulation

v0 = Dict(:p=>[5.1, -0.9])

### Specify job to be run

job = BasicMCJob(model, sampler, mcrange, v0)

### Run the simulation

run(job)

### Get simulated values

chain = output(job)

chain.value

### Check that the simulated values are close to the zero-mean target

mean(chain)

reset(job, [3.2, 9.4])

run(job)

chain = output(job)

job = BasicMCJob(model, sampler, mcrange, v0, tuner=VanillaMCTuner(verbose=true))

run(job)

chain = output(job)

outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget])

job = BasicMCJob(model, sampler, mcrange, v0, tuner=VanillaMCTuner(verbose=true), outopts=outopts)

run(job)

chain = output(job)

chain.logtarget

outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget], :diagnostics=>[:accept])

job = BasicMCJob(model, sampler, mcrange, v0, tuner=VanillaMCTuner(verbose=true), outopts=outopts)

filepath = dirname(@__FILE__)
filesuffix = "csv"

outopts = Dict{Symbol, Any}(
  :monitor=>[:value, :logtarget],
  :diagnostics=>[:accept],
  :destination=>:iostream,
  :filepath=>filepath
)

job = BasicMCJob(model, sampler, mcrange, v0, tuner=VanillaMCTuner(verbose=true), outopts=outopts)

run(job)

for fname in (:value, :logtarget, :diagnosticvalues)
  rm(joinpath(filepath, string(fname)*"."*filesuffix))
end

outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget], :diagnostics=>[:accept])

job = BasicMCJob(model, sampler, mcrange, v0, tuner=VanillaMCTuner(verbose=true), outopts=outopts, plain=false)

run(job)

chain = output(job)

reset(job, [-2.8, 3.4])

run(job)

chain = output(job)

plogtarget(z::Vector{Float64}) = -dot(z, z)

pgradlogtarget(z::Vector{Float64}) = -2*z

p = BasicContMuvParameter(:p, logtarget=plogtarget, gradlogtarget=pgradlogtarget)

model = likelihood_model(p, false)

### Set driftstep to 0.9

sampler = MALA(0.9)

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

v0 = Dict(:p=>[5.1, -0.9])

### Save grad-log-target along with the chain (value and log-target)

outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget, :gradlogtarget], :diagnostics=>[:accept])

job = BasicMCJob(model, sampler, mcrange, v0, tuner=VanillaMCTuner(verbose=true), outopts=outopts)

run(job)

chain = output(job)

chain.gradlogtarget

mean(chain)

job = BasicMCJob(
  model,
  sampler,
  mcrange,
  v0,
  tuner=AcceptanceRateMCTuner(0.6, verbose=true),
  outopts=outopts
)

run(job)

chain = output(job)

plogtarget(z::Vector) = -dot(z, z)

p = BasicContMuvParameter(:p, logtarget=plogtarget, autodiff=:forward)

model = likelihood_model(p, false)

sampler = MALA(0.9)

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

v0 = Dict(:p=>[5.1, -0.9])

outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget, :gradlogtarget], :diagnostics=>[:accept])

job = BasicMCJob(model, sampler, mcrange, v0, tuner=VanillaMCTuner(verbose=true), outopts=outopts)

run(job)

chain = output(job)

plogtarget(z::Vector) = -dot(z, z)

p = BasicContMuvParameter(:p, logtarget=plogtarget, autodiff=:reverse, init=[(:z, ones(2))])

model = likelihood_model(p, false)

sampler = MALA(0.9)

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

v0 = Dict(:p=>[5.1, -0.9])

outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget, :gradlogtarget], :diagnostics=>[:accept])

job = BasicMCJob(model, sampler, mcrange, v0, tuner=VanillaMCTuner(verbose=true), outopts=outopts)

run(job)

chain = output(job)

p = BasicContMuvParameter(:p, logtarget=:(-dot(z, z)), autodiff=:reverse, init=[(:z, ones(2))])

model = likelihood_model(p, false)

sampler = MALA(0.9)

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

v0 = Dict(:p=>[5.1, -0.9])

outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget, :gradlogtarget], :diagnostics=>[:accept])

job = BasicMCJob(model, sampler, mcrange, v0, tuner=VanillaMCTuner(verbose=true), outopts=outopts)

run(job)

chain = output(job)


text:
function Fuhrer1(a1,ap,Varey,w,gamm,Varep,Dbig,qy,qpi,qf,bet)
#Fuhrer1 Generate A,B,K,Q,U, and R from parameter values of simplified Fuhrer model.
#
#
#  Purpose:   Generate A,B,K,Q,U, and R from parameter values of simplified
#             Fuhrer model.
#
#  Usage:     (A,B,K,Q,U,R) = Fuhrer1(a1,ap,Varey,w,gamm,Varep,Dbig,qy,qpi,qf,bet)
#
#  The period loss function is
#
#  qy*y^2 +  qpi*pi^2 + qf*i^2, which is rewritten here as
#
#  x'Qx + 2x'Uu + u'Ru.
#
#  x = [x1,x2], with x1(t) = [ep(t),y(t),Dx(t-1)] and x2(t) = [R[t],Dx(t)]
#
#  See lecture notes for further comments.
#
#
#
#
#
#
#  Paul Soderlind, 20 August 1997
#----------------------------------------------------------------------------

                             #model->general setup
  Ax = [ (  -1/(w*(1-w)^2)                        );
         (  -gamm*a1/(w*(1-w)) - gamm/(1-w)^2     );
         (  -1                                    );
         (  -gamm*ap/(w*(1-w))                    );
         (   2   ) ]'

  Ar = [ 0  0  0  ((1+Dbig)/Dbig)  (4/Dbig*(1-w) ) ]
  Ar =  Ar + 4/Dbig*w*Ax

  A = [ zeros(1,5);                    #x(t+1)=A*x(t)+B*u(t)+ e(t+)
        0  a1   0   ap  0;
        0  0    0    0  1;
        Ar;
        Ax ]

  B = [ 0;0;0;(-1/Dbig);0 ]

  K  = [ 0  1  0         0   0      0;   #[y,pi,f] = K*[x,u]
         0  0  4*(1-w)   0   4*w    0;
         0  0  0         0   0      1  ]

  Wbig = K'Diagonal([qy;qpi;qf])*K

  Q    = Wbig[1:5,1:5]
  U    = Wbig[1:5,6]
  R    = Wbig[6,6]

  return A,B,K,Q,U,R

end
#-----------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
# This file was generated, do not modify it. # hide
select!(df, Not([:yr_renovated, :sqft_basement, :zipcode]));


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.query;
#
# import java.util.function.Function;
#
# public interface QueryFuture<T> {
#
#     T get();
#
#     default <U> QueryFuture<U> map(Function<T, U> function) {
#         return new Mapped<>(this, function);
#     }
#
#     class Mapped<T, U> implements QueryFuture<U> {
#
#         private final QueryFuture<T> queryFuture;
#         private final Function<T, U> function;
#
#         public Mapped(QueryFuture<T> queryFuture, Function<T, U> function) {
#             this.queryFuture = queryFuture;
#             this.function = function;
#         }
#
#         @Override
#         public U get() {
#             return function.apply(queryFuture.get());
#         }
#     }
# }


text:
function update!(d::TreeStochastic, m::Model)
    d.distr = d.eval(m)
    d
end

function update!(d::TreeLogical, m::Model)
    d.value = d.eval(m)
    d
end

function names(d::TreeStochastic, nodekey::Symbol)
    AbstractString["Tree_height["*string(nodekey)*"]" , "Tree_length["*string(nodekey)*"]"]
end

function names(d::TreeLogical, nodekey::Symbol)
    AbstractString["Tree_height["*string(nodekey)*"]" , "Tree_length["*string(nodekey)*"]"]
end

function unlist(root::N) where N<:GeneralNode
    [tree_height(root), tree_length(root)]
end

function unlist(d::T) where T <: TreeVariate
    unlist(d.value)
end

function unlist_tree(root::N) where N<:GeneralNode
    [tree_height(root), tree_length(root)]
end

function unlist_tree(s::AbstractTreeStochastic)
    unlist_tree(s.value)
end

function unlist(s::AbstractTreeStochastic, x::N, transform::Bool=false)  where N <: GeneralNode
    s.value
end


text:
# using CUDA
using StaticArrays
using Interpolations
using LinearAlgebra
using DifferentialEquations
using Cubature
using DelimitedFiles
# using Interpolations
# using Plots
using Dierckx
using Distributed
using IntelVectorMath
using LoopVectorization
using BenchmarkTools
# using GR
# import Pkg; Pkg.add("SymPy")
# using SymPy
# const sympy_parsing_mathematica =
#     SymPy.PyCall.pyimport("sympy.parsing.mathematica")
# mathematica2julia(
#     s::AbstractString,
#     substitutions::Pair{<:AbstractString,<:AbstractString}...,
# ) = SymPy.walk_expression(sympy_parsing_mathematica."mathematica"(
#     s,
#     Dict(substitutions...),
# ))
#
# Pkg.add("SymPy")
#
# Pkg.build("PyCall")
#
# ENV["PYTHON"]
#
# ENV["PYTHON"]="python3"
#
# using Conda
# rm(Conda.ROOTENV, recursive=true)
# Conda.add("numpy")
# Conda.runconda(`search numpy=1.15.4`)
# Pkg.add("Conda")
#
# rm(raw"C:\Users\tyy\.juliapro\JuliaPro_v1.4.2-1\conda", recursive=true)
# Pkg.build("PyCall")
# #
# #
# #
# #
# # Pkg.add("Conda")
# #
# # using Conda
# #
# # Conda.add("numpy")

max_stars_repo_name: Yangyang-Tan/FRG-Code-Collection

max_issues_repo_name: Yangyang-Tan/FRG-Code-Collection

max_forks_repo_name: Yangyang-Tan/FRG-Code-Collection


text:
using XLSX, DataFrames
include("functions.jl")
#This file contains the exported functions
################################################################################
##########################  Default wind profile ###############################
df = DataFrame(XLSX.readtable("../cordoba/src/economics/input_data/data.xlsx", "wind_data")...)
wnd=wndF_wndPrf([getproperty(df,Symbol("Norther"))])
wind_module.save_wind4_module(wnd,"Norther")
################################################################################
########################### Base equipment functions  ##########################
################################################################################

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#NOTE EENS is set to zero for flexlan !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#Cost of HVAC cable of capacity mva and length km
function AC_cbl(mva,km)
    cbl=cable();cbl.mva=mva;cbl.wnd="Norther"
    ac_cbl=optimal_ac_cable(cbl,get_220kV_cables(),km,get_Cost_Data())
    if (ac_cbl.num<1); ac_cbl.costs.cpx_p=1e12; end
    return ac_cbl
end
#Cost of HVDC cable of capacity mva and length km

function DC_cbl(mva,km)
    cbl=cable();cbl.mva=mva;cbl.wnd="Norther"
    dc_cbl=optimal_hvdc_cable(cbl,get_300kV_cables(),km,get_Cost_Data())
    return dc_cbl
end
#Cost of AC platform of capacity mva and depth m
function AC_plat(mva, m::Float64=17.1)
    plat=platform();
    plat.mva=mva
    plat.wnd="Norther"
    plat=cost_ac_platform(plat,get_Cost_Data())
    plat=adjust_base_ac_platform(plat,get_Cost_Data())
    #adjusts costs based on depth default and neutral depth is 17.1m
    multiplier=0.0136*abs(m) + 0.7676
    plat.costs.cm=plat.costs.cm*multiplier
    plat.costs.cpx=plat.costs.cpx*multiplier
    plat.costs.ttl=plat.costs.cpx+plat.costs.cm
    return plat
end
#Cost of DC platform of capacity mva and depth m
function DC_plat(mva, m::Float64=17.1)
    plat=platform();
    plat.mva=mva
    plat.wnd="Norther"
    plat=cost_dc_platform(plat,get_Cost_Data())
    plat=adjust_base_dc_platform(plat,get_Cost_Data())
    #adjusts costs based on depth default and neutral depth is 17.1m
    multiplier=0.0136*abs(m) + 0.7676
    plat.costs.cm=plat.costs.cm*multiplier
    plat.costs.cpx=plat.costs.cpx*multiplier
    plat.costs.ttl=plat.costs.cpx+plat.costs.cm
    return plat
end
#offshore transformer of capacity mva
function off_xfm(mva)
    xfo_off=transformer()
    xfo_off.relia=get_offshore_xfo_failure_data(xfo_off)
    xfo_off.wnd="Norther"
    xfo_off.mva=mva
    xfo_off=xfo_oss(xfo_off,get_Cost_Data(),get_Xfo_Data())
    return xfo_off
end
#onshore transformer of capacity mva
function on_xfm(mva)
    xfo_on=transformer()
    xfo_on.relia=get_offshore_xfo_failure_data(xfo_on)
    xfo_on.wnd="Norther"
    xfo_on.mva=mva
    xfo_pcc(xfo_on,get_Cost_Data(),get_Xfo_Data())
    return xfo_on
end
#cost of offshore converter of size mva
function off_conv(mva)
    conv_off=converter()
    conv_off.relia=get_offshore_conv_failure_data(conv_off)
    conv_off.wnd="Norther"
    conv_off.mva=mva
    conv_off.elec.mva=mva
    conv_off.num=1
    conv_off=cost_hvdc_oss(conv_off,get_Cost_Data())
    conv_off=adjust_base_hvdc_offshore_converter(conv_off,get_Cost_Data())
    return conv_off
end
#cost of onshore converter of size mva
function on_conv(mva)
    conv_on=converter()
    conv_on.relia=get_onshore_conv_failure_data(conv_on)
    conv_on.wnd="Norther"
    conv_on.mva=mva
    conv_on.elec.mva=mva
    conv_on.num=1
    conv_on=cost_hvdc_pcc(conv_on,get_Cost_Data())
    conv_on=adjust_base_hvdc_onshore_converter(conv_on,get_Cost_Data())
    return conv_on
end
################################################################################
########################### Finding MV max distance ############################
################################################################################
#**
function find_max_mv_transmission(owpps,database)
    ks=get_Cost_Data()
    for owpp in owpps
        pnt0=0.001
        pnt1=deepcopy(database["cables"]["66.0"][string(owpp.mva)][length(database["cables"]["66.0"][string(owpp.mva)])].length)
        mv_cable_pnt0=deepcopy(mvac_cable(owpp.mva,pnt0,owpp.wnd,database["cables"]["66.0"][string(owpp.mva)],ks))
        mv_cable_pnt1=deepcopy(mvac_cable(owpp.mva,pnt1,owpp.wnd,database["cables"]["66.0"][string(owpp.mva)],ks))
        plat=platform();plat.costs.ttl=0.0
        plat.mva=owpp.mva
        plat.wnd=owpp.wnd
        plat=cost_ac_platform(plat,ks)
        plat=adjust_base_ac_platform(plat,ks)
        ####################################### Addition
        xfo_off=transformer()
        xfo_off.relia=get_offshore_xfo_failure_data(xfo_off)
        xfo_off.wnd=owpp.wnd
        xfo_off.mva=owpp.mva
        xfo_off=xfo_oss(xfo_off,ks,get_Xfo_Data())
        ########################################
        plat.costs.ttl=plat.costs.ttl+xfo_off.costs.ttl
        np=0
        if (database["bits"]["hvac"]==true || database["bits"]["mpc_ac"]==true)
            hv_cable_pnt0_220=deepcopy(hvac_cable(owpp.mva,pnt0,owpp.wnd,database["cables"]["220.0"][string(owpp.mva)],ks))
            hv_cable_pnt0_400=deepcopy(hvac_cable(owpp.mva,pnt0,owpp.wnd,database["cables"]["400.0"][string(owpp.mva)],ks))
            if (hv_cable_pnt0_220.costs.ttl<hv_cable_pnt0_400.costs.ttl)
                hv_cable_pnt0=hv_cable_pnt0_220
                hv_c="220.0"
            else
                hv_cable_pnt0=hv_cable_pnt0_400
                hv_c="400.0"
            end
        elseif (database["bits"]["hvdc"]==true)
            hv_c="300.0"
        end

        if (mv_cable_pnt0.costs.ttl<hv_cable_pnt0.costs.ttl+plat.costs.ttl)
            hv_cable_pnt1=deepcopy(hvac_cable(owpp.mva,pnt1,owpp.wnd,database["cables"][hv_c][string(owpp.mva)],ks))
            #if AC cable range is maxed out this backs it up to a feasible distance
            while (hv_cable_pnt1.costs.ttl==Inf)
                pnt1=pnt1-25
                hv_cable_pnt1=deepcopy(hvac_cable(owpp.mva,pnt1,owpp.wnd,database["cables"][hv_c][string(owpp.mva)],ks))
                mv_cable_pnt1=deepcopy(mvac_cable(owpp.mva,pnt1,owpp.wnd,database["cables"]["66.0"][string(owpp.mva)],ks))
            end
            np=pnt1
            if (mv_cable_pnt1.costs.ttl>hv_cable_pnt1.costs.ttl+plat.costs.ttl)
                while (pnt1-pnt0>1)
                    np=(pnt0+pnt1)/2
                    mv_cable_np=deepcopy(mvac_cable(owpp.mva,np,owpp.wnd,database["cables"]["66.0"][string(owpp.mva)],ks))
                    hv_cable_np=deepcopy(hvac_cable(owpp.mva,np,owpp.wnd,database["cables"][hv_c][string(owpp.mva)],ks))
                    if (mv_cable_np.costs.ttl<hv_cable_np.costs.ttl+plat.costs.ttl)
                        pnt0=deepcopy(np)
                    elseif (mv_cable_np.costs.ttl>hv_cable_np.costs.ttl+plat.costs.ttl)
                        pnt1=deepcopy(np)
                    else
                        break
                    end
                end
            end
        end
        ############### addition #########
        # Account for collector circuit###
        ###### 6MW/km square assumed######
        #### Diameter of area travelled###
        collector=2*sqrt((owpp.mva/6)/pi)
        np_minus_collector=np-collector
        if (np_minus_collector<0)
            #np_minus_collector=0
            np_minus_collector=0.1
        end
        owpp.mv_zone=deepcopy(np_minus_collector)
        ##################################
        #owpp.mv_zone=deepcopy(np)
        owpp.kV=66.0
    end
    return owpps
end

################################################################################
############################## Finding lookup table ############################
################################################################################
#**
function get_equipment_tables(km_mva_set,wnd,kv)
    cable_database=get_cable_table(km_mva_set,wnd)
    xfo_database,conv_database=get_xfoConv_table(km_mva_set,wnd)
    database=Dict{String, Dict{String,Any}}()
    push!(database,"cables"=>cable_database)
    push!(database,"transformers"=>xfo_database)
    push!(database,"converters"=>conv_database)

    #vheck if hvdc is to be included
    bits=Dict([("hvdc", true), ("mpc_ac", false), ("hvac", true)])
    push!(database,"bits"=>bits)
    hvdc_bit=check4hvdc(database,km_mva_set,wnd,kv)

    #check hvac connections
    hvac_bit=true
    if (hvdc_bit==true)
        hvac_bit=check4hvac(database,km_mva_set,wnd,kv)
    end
    mpc_bit=false

    #final bits sets
    bits=Dict([("hvdc", hvdc_bit), ("mpc_ac", mpc_bit), ("hvac", hvac_bit)])
    push!(database,"bits"=>bits)
    database["cables"]=set_cables_per_km_cost(database["cables"],km_mva_set)
    return database
end

################################################################################
######################## Optimal equipment using lookup table ##################
################################################################################
#finds the offshore transformer using the look up table
#when increasing meshed transformer xfo_oss_meshed is used
#**
function xfo_oss(xfo0,ks,xfo_data)
    xfo=transformer()
    xfo.costs.ttl=Inf
    xfo.mva=xfo0.mva
    for xd in xfo_data
        if ((xd+10<=xfo0.mva) && (6*xd>xfo0.mva))
            xfo0.num=1
            xfo0.elec.mva=xd
            while ((xfo0.num*xfo0.elec.mva)<xfo0.mva)
                xfo0.num=xfo0.num+1
            end
            xfo0=cost_xfo_oss(xfo0,ks)
            if (xfo0.costs.ttl<xfo.costs.ttl)
                xfo=deepcopy(xfo0)
            end
        end
    end
    return xfo
end

#finds the onshore transformer using the look up table
#**
function xfo_pcc(xfo0,ks,xfo_data)
    xfo=transformer()
    xfo.costs.ttl=Inf
    xfo.mva=xfo0.mva
    for xd in xfo_data
        if ((xd+10<=xfo0.mva) && ((3*xd>xfo0.mva) || ((xfo0.mva>2250) && (6*xd>xfo0.mva))))
            xfo0.num=1
            xfo0.elec.mva=xd
            while ((xfo0.num*xfo0.elec.mva)<xfo0.mva)
                xfo0.num=xfo0.num+1
            end
            xfo0=cost_xfo_pcc(xfo0,ks)
            if (xfo0.costs.ttl<xfo.costs.ttl)
                xfo=deepcopy(xfo0)
            end
        end
    end
    return xfo
end

#finds the hvdc cable using the look up table
#**
function hvdc_cable(mva,km,wnd,cable_array,ks)
    cbl=deepcopy(cable_array[1])
    cbl.wnd=wnd
    cbl.length=km
    cbl=cost_hvdc_cable(cbl,ks)
    return cbl
end

#finds the hvac cable using the look up table
#**
function hvac_cable(mva,km,wnd,cable_array,ks)
    cbl=cable()
    cbl.costs.ttl=Inf
    cbl.costs.grand_ttl=Inf
    for cbl0 in cable_array
        if (km<=cbl0.length+0.01)
            #+0.01 is just to be certain floating point error does not occur
            cbl=deepcopy(cbl0)
            break
        end
    end
    if (cbl.size!=0.0)
        cbl.wnd=wnd
        cbl.length=km
        cbl.elec.mva=get_newQ_Capacity(cbl.elec.freq,km,cbl.elec.volt,cbl.elec.farrad,cbl.elec.amp)
        cbl=get_cbl_failure_data(cbl)
        cbl=cost_hvac_cable(cbl,ks)
    end
    return cbl
end

#finds the mvac cable using the look up table
#**
function mvac_cable(mva,km,wnd,cable_array,ks)
    cbl=cable()
    cbl.costs.ttl=Inf
    cbl.costs.grand_ttl=Inf
    for cbl0 in cable_array
        if (km<=cbl0.length+0.001)#+0.001 is just to be certain floating point error does not occur
            cbl=deepcopy(cbl0)
            break
        end
    end
    max_in_parallel=12#if changing check economics/database/functions function optimal_ac_cable(cbl0,cbl_data,km,ks)
    if (cbl.num>max_in_parallel && cbl.size!=0.0)#66kv up to 2000MW eception
        cbl.mva=cbl.mva/2
        cbl.num=cbl.num/2
        cbl.costs.perkm_cpx=cbl.costs.perkm_cpx/2
        cbl.elec.mva=get_newQ_Capacity(cbl.elec.freq,km,cbl.elec.volt,cbl.elec.farrad,cbl.elec.amp)
        cbl.wnd=wnd
        cbl.length=km
        cbl=get_cbl_failure_data(cbl)
        cbl=cost_mvac_cable(cbl,ks)
        cbl.num=cbl.num*2
        cbl.mva=cbl.mva*2
        cbl.costs.perkm_cpx=cbl.costs.perkm_cpx*2
        cbl.costs.sg=cbl.costs.sg*2
        cbl.costs.cpx_p=cbl.costs.cpx_p*2
        cbl.costs.cpx_i=cbl.costs.cpx_i*2
        cbl.costs.rlc=cbl.costs.rlc*2
        cbl.costs.cm=cbl.costs.cm*2
        cbl.costs.eens=cbl.costs.eens*2
        cbl.costs.ttl=cbl.costs.ttl*2
        cbl.costs.grand_ttl=cbl.costs.ttl
    elseif (cbl.size!=0.0)
        cbl.elec.mva=get_newQ_Capacity(cbl.elec.freq,km,cbl.elec.volt,cbl.elec.farrad,cbl.elec.amp)
        cbl.wnd=wnd
        cbl.length=km
        cbl=get_cbl_failure_data(cbl)
        cbl=cost_mvac_cable(cbl,ks)
        if (cbl.costs.ttl==Inf)
            cbl.costs.perkm_ttl=10^9
            cbl.costs.grand_ttl=0.0
        else
            cbl.costs.grand_ttl=cbl.costs.ttl
        end
    end
    return cbl
end

#Converter
#calculates the cost of an offshore HVDC converter
#each line needing a converter shares the total before the fixed cost is added at the end
#**
function cost_hvdc_oss(conv,ks)
    #capex converter
    conv.costs.cpx=capex_hvdc(conv,ks)
    #cost of losses
    conv.costs.tlc=cost_tlc(conv,ks)
    #corrective maintenance
    conv.costs.cm=cost_cm(conv.costs.cpx,ks.opx_co)
    #eens calculation
    conv.costs.eens=cost_eens(conv,ks)
    #totals the xfo cost
    conv.costs.ttl=cost_conv_sum(conv)
    return conv
end

#once all incoming lines with converters are combined the fixed cost is added
#**
function adjust_base_hvdc_offshore_converter(conv,ks)
    #adjust to add base cost capex converter
    conv.costs.cpx=conv.costs.cpx+ks.conv_d
    #corrective maintenance
    conv.costs.cm=conv.costs.cm+ks.conv_d*ks.opx_co*npv_years()
    #totals the converter cost
    conv.costs.ttl=cost_conv_sum(conv)
    return conv
end

#calculates the cost of an onshore HVDC converter
#each line needing a converter shares the total before the fixed cost is added at the end
#**
function cost_hvdc_pcc(conv,ks)
    #capex converter
    conv.costs.cpx=capex_hvdc(conv,ks)
    #cost of losses
    conv.costs.tlc=cost_tlc(conv,ks)
    #corrective maintenance
    conv.costs.cm=cost_cm(conv.costs.cpx,ks.opx_cp)
    #eens calculation
    conv.costs.eens=cost_eens(conv,ks)
    #totals the xfo cost
    conv.costs.ttl=cost_conv_sum(conv)
    return conv
end

#once all incoming lines with converters are combined the fixed cost is added
#**
function adjust_base_hvdc_onshore_converter(conv,ks)
    #adjust to add base cost capex converter
    conv.costs.cpx=conv.costs.cpx+ks.conv_d
    #corrective maintenance
    conv.costs.cm=conv.costs.cm+ks.conv_d*ks.opx_cp*npv_years()
    #totals the converter cost
    conv.costs.ttl=cost_conv_sum(conv)
    return conv
end

################################### platform ###################################
#calculates incremental cost of an AC platform NOT including fixed cost
#**
function cost_ac_platform(plat,ks)
    #capex oss
    plat.costs.cpx=capex_plat_ac(plat,ks)
    #corrective maintenance
    plat.costs.cm=cost_cm(plat.costs.cpx,ks.opx_pl)
    #totals the plat cost
    plat.costs.ttl=plat.costs.cpx+plat.costs.cm
    return plat
end

#once all incoming lines established the fixed cost is added to the of incremental sum - AC
#**
function adjust_base_ac_platform(plat,ks)
    #adjust to add base cost capex converter
    plat.costs.cpx=plat.costs.cpx+ks.pac_f
    #corrective maintenance
    plat.costs.cm=plat.costs.cm+ks.pac_f*ks.opx_pl*npv_years()
    #totals the converter cost
    plat.costs.ttl=plat.costs.cpx+plat.costs.cm
    return plat
end

#calculates incremental cost of an DC platform NOT including fixed cost
#**
function cost_dc_platform(plat,ks)
    #capex oss
    plat.costs.cpx=capex_plat_dc(plat,ks)
    #corrective maintenance
    plat.costs.cm=cost_cm(plat.costs.cpx,ks.opx_pl)
    #totals the plat cost
    plat.costs.ttl=plat.costs.cpx+plat.costs.cm
    return plat
end

#once all incoming lines established the fixed cost is added to the of incremental sum - DC
#**
function adjust_base_dc_platform(plat,ks)
    #adjust to add base cost capex converter
    plat.costs.cpx=plat.costs.cpx+ks.pdc_h
    #corrective maintenance
    plat.costs.cm=plat.costs.cm+ks.pdc_h*ks.opx_pl*npv_years()
    #totals the converter cost
    plat.costs.ttl=plat.costs.cpx+plat.costs.cm
    return plat
end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.@schedule
# Date: 12/11/2014

module AbstractGenerativeModelImpl

export AbstractGenerativeModel

abstract AbstractGenerativeModel

end #module

module AbstractGenerativeModelInterfaces

export get,
       isterminal

function get() end

function isterminal() end

end #module


text:
module ExploringOceanData

export @boilerplate, @analysisplate

macro boilerplate()
    boiler_block = :(
        using GLMakie;
        using Pkg;
        using Statistics;
        using NCDatasets;
        using DataStructures;
    )
    return boiler_block
end

macro analysisplate()
    boiler_block = :(
        using Statistics;
        using NCDatasets;
        using DataStructures;
    )
    return boiler_block
end

end # module


text:
"""
    catbyrow(arr)

Convert an array of array to a 2-dimenional matrix, the inner vector
is transposed as concatenated as rows.
"""
catbyrow(arr) = vcat(transpose(arr)...);

"""
    toparameter(parameters)

Convert a parameter vector of length 1 to a single parameter.
"""
function toparameter(parameters::Vector{Symbol})
    length(parameters) == 1 || error("Number of parameters should be 1")
    return parameters[1]
end

effsamplesize(n1, n2) = 1/(1/n1 + 1/n2)

dofpooled(n1, n2) = n1 + n2 - 2

function dofwelch(sd1, sd2, n1, n2)
    se1_sq = sd1^2/n1
    se2_sq = sd2^2/n2
    (se1_sq + se2_sq)^2 / (se1_sq^2/(n1-1) + se2_sq^2/(n2-1))
end

pooledsd(sd1, sd2, n1, n2) = sqrt(((n1-1)*sd1^2 + (n2-1)*sd2^2) / (n1 + n2 - 2))

zstat(x̄, μ0, σ, n) = sqrt(n)*(x̄ - μ0)/σ


text:
module ReinforcementLearning

export RL
const RL = ReinforcementLearning

using Reexport

@reexport using ReinforcementLearningBase
@reexport using ReinforcementLearningCore
@reexport using ReinforcementLearningEnvironments
@reexport using ReinforcementLearningZoo

end


text:
import Random
using Test: @test, @testset, @inferred
using Statistics: mean, std, var
using LogExpFunctions: softmax
using ConvolutionalRBMs: randgumbel, cartesian_sample_from_logits

@testset "randgumbel" begin
    Random.seed!(3)
    data = [@inferred randgumbel() for _ = 1:10^6]
    @test mean(data) ≈ MathConstants.γ rtol=0.01
    @test std(data) ≈ π / √6 rtol=0.01
end

@testset "categorical_sample_from_logits" begin
    Random.seed!(52)
    logits = randn(2, 3, 2)
    @test size(@inferred cartesian_sample_from_logits(logits; dims=2)) == (2,1,2)
    p = mean(cartesian_sample_from_logits(logits; dims=2) .== CartesianIndices(logits) for _ in 1:10^6)
    @test p ≈ softmax(logits; dims=2) rtol=0.01
end


text:
using LinearAlgebra
export prune, prune_parent, rm_redundant_rotation, find_smallest, prune_train

function prune(x::RotationGate; threshold=1e-1)
    if abs(x.theta) < threshold
        return NoGate()
    else
        return x
    end
end

prune(x::AbstractBlock; threshold=1e-1) =
    chsubblocks(x, prune.(subblocks(x); threshold=threshold))

prune_parent(c::NoGate, p::RotationGate) = NoGate()

function prune_parent(c::ChainBlock, p::ChainBlock)
    c[1] isa ControlBlock && return copy(c)
    chsubblocks(p, prune_parent.(subblocks(c), subblocks(p)))
end

prune_parent(c::AbstractBlock, p::AbstractBlock) =
    chsubblocks(p, prune_parent.(subblocks(c), subblocks(p)))

function rm_cnot(x::ChainBlock{N}) where N
    x[1] isa ControlBlock || return x
    deleteat!(x.blocks, rand(1:length(x)))
    return x
end

rm_cnot(x::AbstractBlock) = chsubblocks(x, rm_cnot.(subblocks(x)))

function rm_redundant_rotation(x::ChainBlock{1})
    length(x) === 3 || return x
    # this is an identity
    if tr((mat(x) - IMatrix(2))' * (mat(x) - IMatrix(2))) < 1e-2
        return NoGate()
    end

    if x[2] isa NoGate
        return Rz(x[1].theta + x[3].theta)
    end
    return x
end

rm_redundant_rotation(x::AbstractBlock) =
    chsubblocks(x, rm_redundant_rotation.(subblocks(x)))

"""
    find_smallest(ngates, circuit)

find the first `ngates` smallest absolute value of parameters
of gates in the given circuit.
"""
function find_smallest(ngates, circuit)
    return sort!(parameters(circuit); by=abs)[ngates]
end

"""
    prune_train(opt, circuit, ham; verbose=false, nprune=10, nepochs=10, niteration=100)

# Parameters
- `opt`: optimizer from `Flux.Optimise`.
- `circuit`: the circuit to prune
- `ham`: the Hamiltonian

# Keywords

- `verbose`: print the verbose info in each training process if `true`
- `nprune`: number of gates to prune
- `nepochs`: the number of epochs of pruning procedure
- `niteration`: number of iterations for training
"""
function prune_train(opt, circuit, ham; verbose=false, nprune=10, nepochs=10, niteration=100, relative_error=true, least_prune=10, cnot_prune=true)
    total_history = Float64[]
    n = nqubits(ham)
    original = copy(circuit)
    if relative_error
        E0 = eigmin(Matrix(mat(ham)))/4n
    end

    try
        for k in 1:nepochs
            @info "epoch = $k"
            @info "nparameters = $(nparameters(circuit))"
            history = train!(opt, niteration, ham, circuit; verbose=verbose)
            append!(total_history, history)
            @info "E/n = $(history[end])"

            if relative_error
                re = (history[end] - E0)/E0
                @info "(E - E0)/E0 = $re"

                if k > least_prune && re > 1e-2
                    break
                end
            end

            circuit = prune(circuit; threshold=abs(find_smallest(nprune, circuit)));
            circuit = rm_redundant_rotation(circuit);
            circuit = copy(prune_parent(circuit, original))

            if cnot_prune
                circuit = chsubblocks(circuit, rm_cnot.(subblocks(circuit)))
            end
        end
    catch e
        if e isa InterruptException
            return circuit, total_history
        else
            rethrow(e)
        end
    end
    return circuit, total_history
end


text:
#  Fast, stripped down expression evaluation code for investigating
#  speed of doeval, meval.  In particular, Table takes unevaled args
#  and calls doeval. So we call meval1 instead in Table (for testing)

# This is stripped down evaluation

@inline meval1(x::Complex) = x.im == 0 ? x.re : x
@inline meval1(x) = x
@inline meval1(s::SJSym) = symval(s)

function meval1(mx::Mxpr)
    args = margs(mx)
    len = length(args)
    nargs = newargs(len)
#   evalargs(args,nargs)
    @inbounds for i in 1:len
        res = meval1(args[i])
        nargs[i] = res
    end
#    mx.args = nargs
    nmx = mxpr(mhead(mx),nargs)
    nmx = flatcanon!(nmx)
    nmx = ev_upvalues(nmx)
    nmx = ev_downvalues(nmx)
    merge_args_if_emtpy_syms(nmx)
    apprules(nmx)
end


text:
@inline function ArbFloat(x)
    prec = DEFAULT_PRECISION.x
    res  = ArbFloat{prec}(x)
    return res
end
@inline function ArbReal(x)
    prec = DEFAULT_PRECISION.x
    res  = ArbReal{prec}(x)
    return res
end

@inline function ArbComplex(x)
    prec = DEFAULT_PRECISION.x
    res  = ArbComplex{prec}(x)
    return res
end

@inline function ArbComplex(x, y)
    prec = DEFAULT_PRECISION.x
    res  = ArbComplex{prec}(x, y)
    return res
end

@inline function ArbComplex(x, y::T) where {T<:AbstractFloat}
    prec = DEFAULT_PRECISION.x
    res  = ArbComplex{prec}(x, y)
    return res
end

# IEEEFloat

# rounds up (widens)
function Float64(x::Mag)
    z = ccall(@libarb(mag_get_d), Float64, (Ref{Mag},), x)
    return z
end

function Float64(x::ArbFloat{P}) where {P}
    rounding = match_rounding_mode(RoundNearest)
    z = ccall(@libarb(arf_get_d), Float64, (Ref{ArbFloat}, Cint), x, rounding)
    return z
end
Float32(x::ArbFloat{P}) where {P} = Float32(Float64(x))
Float16(x::ArbFloat{P}) where {P} = Float16(Float64(x))

function Float64(x::ArbReal{P}) where {P}
    x = midpoint(x)
    y = ArbFloat{P}(x)
    Float64(y)
end
Float32(x::ArbReal{P}) where {P} = Float32(Float64(x))
Float16(x::ArbReal{P}) where {P} = Float16(Float64(x))

function Float64(x::ArbComplex{P}) where {P}
    x = midpoint(real(x))
    y = ArbFloat{P}(x)
    Float64(y)
end
Float32(x::ArbComplex{P}) where {P} = Float32(Float64(x))
Float16(x::ArbComplex{P}) where {P} = Float16(Float64(x))

function Float64(x::ArbFloat{P}, roundingmode::RoundingMode) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = ccall(@libarb(arf_get_d), Float64, (Ref{ArbFloat}, Cint), x, rounding)
    return z
end
Float32(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Float32(Float64(x, roundingmode))
Float16(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Float16(Float64(x, roundingmode))

function Float64(x::ArbReal{P}, roundingmode::RoundingMode) where {P}
    x = midpoint(x)
    y = ArbFloat{P}(x)
    Float64(y, roundingmode)
end
Float32(x::ArbReal{P}, roundingmode::RoundingMode) where {P} = Float32(Float64(x, roundingmode))
Float16(x::ArbReal{P}, roundingmode::RoundingMode) where {P} = Float16(Float64(x, roundingmode))

function Float64(x::ArbComplex{P}, roundingmode::RoundingMode) where {P}
    x = midpoint(real(x))
    y = ArbFloat{P}(x)
    Float64(y, roundingmode)
end
Float32(x::ArbComplex{P}, roundingmode::RoundingMode) where {P} = Float32(Float64(x, roundingmode))
Float16(x::ArbComplex{P}, roundingmode::RoundingMode) where {P} = Float16(Float64(x, roundingmode))

# parts

function real(x::ArbComplex{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(acb_get_real), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x)
    return z
end

function imag(x::ArbComplex{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(acb_get_imag), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x)
    return z
end

function arg(x::ArbComplex{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(acb_arg), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end

# retype

ArbFloat(x::ArbReal{P}) where {P} = ArbFloat{P}(x)
ArbFloat(x::ArbComplex{P}) where {P} = ArbFloat{P}(real(x))
ArbReal(x::ArbFloat{P}) where {P} = ArbReal{P}(x)
ArbReal(x::ArbComplex{P}) where {P} = ArbReal{P}(real(x))
# ArbComplex(x::ArbFloat{P}) where {P} = ArbComplex{P}(ArbReal{P}(x))
# ArbComplex(x::ArbReal{P}) where {P} = ArbComplex{P}(x)

ArbFloat{Q}(x::ArbReal{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(x))
ArbFloat{Q}(x::ArbComplex{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(real(x)))
ArbReal{Q}(x::ArbFloat{P}) where {P,Q} = ArbReal{Q}(ArbFloat{Q}(x))
ArbReal{Q}(x::ArbComplex{P}) where {P,Q} = ArbReal{Q}(real(x))
ArbComplex{Q}(x::ArbFloat{P}) where {P,Q} = ArbComplex{Q}(ArbReal{Q}(ArbFloat{Q}(x)))
ArbComplex{Q}(x::ArbReal{P}) where {P,Q} = ArbComplex{Q}(ArbReal{Q}(x))

@inline convert(::Type{ArbFloat{Q}}, x::ArbReal{P}) where {P,Q} = ArbFloat{Q}(x)
@inline convert(::Type{ArbFloat{Q}}, x::ArbComplex{P}) where {P,Q} = ArbFloat{Q}(x)
@inline convert(::Type{ArbReal{Q}}, x::ArbFloat{P}) where {P,Q} = ArbReal{Q}(x)
@inline convert(::Type{ArbReal{Q}}, x::ArbComplex{P}) where {P,Q} = ArbReal{Q}(x)
@inline convert(::Type{ArbComplex{Q}}, x::ArbFloat{P}) where {P,Q} = ArbComplex{Q}(x)
@inline convert(::Type{ArbComplex{Q}}, x::ArbReal{P}) where {P,Q} = ArbComplex{Q}(x)

# change precision

function ArbFloat{P}(x::ArbFloat{Q}, roundingmode::RoundingMode) where {P,Q}
    rounding = match_rounding_mode(roundingmode)
    z = ArbFloat{P}()
    rnd = ccall(@libarb(arf_set_round), Cint,
                (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, P, rounding)
    return z
end
ArbFloat{P}(x::ArbFloat{Q}) where {P,Q} = ArbFloat{P}(x, RoundNearest)

function ArbReal{P}(x::ArbReal{Q}) where {P,Q}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_round), Cvoid,
          (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, P)
    return z
end

function ArbComplex{P}(x::ArbComplex{Q}) where {P,Q}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_round), Cvoid,
          (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end


text:
using Test
using HoJBot: mod_make_regex, mod_check_message

# initialization bad word list since it's lazy
HoJBot.mod_init()

@testset "mod" begin
    function test_matches(regex, strs, expected)
        for i in 1:length(strs)
            @test (match(regex, strs[i]) !== nothing) == expected[i]
        end
    end

    # regular word
    let regex = mod_make_regex("bad")
        strs = ["bad", "notbad", "badness"]
        expected = [true, false, false]
        test_matches(regex, strs, expected)
    end

    # having symbols in the match word
    let regex = mod_make_regex("b.a.d")
        strs = ["b.a.d", "b1a2d"]
        expected = [true, false]
        test_matches(regex, strs, expected)
    end

    @test isempty(mod_check_message("this is good")) == true
    @test isempty(mod_check_message("this is shit")) == false
    @test isempty(mod_check_message("this is ||shit||")) == true
    @test isempty(mod_check_message("this is damn ||shit||")) == false
    @test isempty(mod_check_message("this is ok ||shit")) == true
end


text:
module Spirograph

using Plots

function lastly()
    plot!(aspectratio = 1, legend = false, axis = false, grid = false, ticks = false)
end

@inline inner_angle(a::Int, b::Int, t::Real)::Float64 = b * t / a

@inline function spot(a::Int, b::Int, offset::Real, t::Real)::Complex{Float64}
    θ = inner_angle(a, b, t)
    return im * ((a + b) * exp(θ * im) + offset * exp(t * im))
end

@inline stop_t(a::Int, b::Int) = 2 * π * (lcm(a, b) ÷ abs(b))

function spiro_points(a::Int, b::Int, offset::Real, step::Real)
    if a <= 0 || b == 0
        error(
            "Radius of fixed circle must be positive and radius of the rolling disk must be nonzero",
        )
    end
    T = stop_t(a, b) #+ step
    return [spot(a, b, offset, t) for t = 0:step:T]
end

_DEFAULT_STEP = 0.05

"""
    spirograph(a,b,offset;args...)
Draw a spirograph picture where
+ `a` is the radius of the fixed wheel,
+ `b` is the radius of the moving wheel, and
+ `offset` is the distance from the center of the moving wheel to the pen.

The optional `args` are passed to `plot`.

## Example
`spirograph(20,-9,12,linecolor=:green)`
"""
function spirograph(a::Int, b::Int, offset::Real, step::Real = _DEFAULT_STEP; args...)
    pts = spiro_points(a, b, offset, step)
    plot(real(pts), imag(pts); args...)
    lastly()
end

"""
    spirograph!
Same as `spirograph` but does not clear the drawing screen first. This way spirograph
images can be combined.
"""
function spirograph!(a::Int, b::Int, offset::Real, step::Real = _DEFAULT_STEP; args...)
    pts = spiro_points(a, b, offset, step)
    plot!(real(pts), imag(pts); args...)
    lastly()
end

export spirograph, spirograph!
end # module


text:
#
# Define tests
#

function fill_ghost_cells!(c, κtop, κbottom, model, fieldbcs)
    fill_bottom_ghost_cell!(fieldbcs.bottom, c, κbottom, model)
    fill_top_ghost_cell!(fieldbcs.top, c, κtop, model, c.grid.N)
    return nothing
end

function cell_field_construction(T, N, H)
    grid = UniformGrid(T, N, H)
    c = CellField(grid)
    typeof(c) <: CellField
end

function face_field_construction(T, N, H)
    grid = UniformGrid(T, N, H)
    f = FaceField(grid)
    typeof(f) <: FaceField
end

function field_indexing()
    grid = UniformGrid(3, 4.2)
    c = CellField(grid)
    f = FaceField(grid)
    val = 2.1
    c[2] = val
    f[2] = val
    c[2] == val && f[2] == val
end

function cell_∂z(T)
    grid = UniformGrid(T, 2, 2.0)
    c = CellField([2, 5], grid)
    cz = ∂z(c)
    cz_answer = FaceField([2, 3, -5], grid)
    cz.data[1:3] == cz_answer.data[1:3]
end

function face_∂z(T)
    grid = UniformGrid(T, 2, 2.0)
    f = FaceField([2, 4, 6], grid)
    fz = ∂z(f)
    fz_answer = CellField([2, 2], grid)
    fz.data[1:2] == fz_answer.data[1:2]
end

function cell_plus(T)
    grid = UniformGrid(T, 4, 2.0)
    c1 = CellField([0, 1, 2, 0], grid)
    c2 = CellField([0, 3, 4, 0], grid)
    c3 = CellField([0, 4, 6, 0], grid)
    c1_plus_c2 = c1 + c2
    c1_plus_c2.data == c3.data
end

function cell_times(T)
    grid = UniformGrid(T, 4, 2.0)
    c1 = CellField([0, 1, 2, 0], grid)
    c2 = CellField([0, 3, 4, 0], grid)
    c3 = CellField([0, 3, 8, 0], grid)
    c1_times_c2 = c1 * c2
    c1_times_c2.data == c3.data
end

function face_plus(T)
    grid = UniformGrid(T, 2, 2.0)
    c1 = FaceField([1, 2, 2], grid)
    c2 = FaceField([3, 4, 3], grid)
    c3 = FaceField([4, 6, 5], grid)
    c1_plus_c2 = c1 + c2
    c1_plus_c2.data == c3.data
end

function face_times(T)
    grid = UniformGrid(T, 2, 2.0)
    f1 = FaceField([1, 2, 2], grid)
    f2 = FaceField([3, 4, 3], grid)
    f3 = FaceField([3, 8, 6], grid)
    f1_times_f2 = f1 * f2
    f1_times_f2.data == f3.data
end

function set_scalar_field(loc, T)
    g = UniformGrid(T, 2, 2.0)
    f = Field(loc, g)
    a = 7
    set!(f, a)
    !any(@. !(f.data == a))
end

function set_array_field(loc, T)
    g = UniformGrid(T, 2, 2.0)
    f = Field(loc, g)
    data = rand(1:10, length(f))
    set!(f, data)
    f.data[1:length(f)] == data
end

function set_function_field(loc, T)
    g = UniformGrid(T, 2, 2.0)
    f = Field(loc, g)
    fcn(z) = z^2
    z = nodes(f)
    set!(f, fcn)
    i = 2
    f[i] == fcn(z[i])
end

function integral_is_correct(T)
    grid = UniformGrid(T, 3, 3.0)
    c = CellField([1, 2, 3], grid)
    integral(c) / grid.H == 2
end

function integral_range(T, N, H, z₋, z₊)
    grid = UniformGrid(T, N, H)
    c = CellField(ones(N), grid)
    integral(c, z₋, z₊) == z₊ - z₋
end

struct DummyModel end

function ghost_cell_value(T)
    model = DummyModel()
    grid = UniformGrid(T, 3, 3.0)
    c = CellField([1, 2, 3], grid)
    c_top = T(0.2)
    c_bottom = T(0.11)
    top_bc = ValueBoundaryCondition(c_top)
    bottom_bc = ValueBoundaryCondition(c_bottom)
    fieldbcs = FieldBoundaryConditions(bottom_bc, top_bc)
    fill_ghost_cells!(c, 0, 0, model, fieldbcs)

    onface(c, 1) ≈ c_bottom && onface(c, grid.N+1) ≈ c_top
end

function ghost_cell_gradient(T)
    model = DummyModel()
    grid = UniformGrid(T, 3, 3.0)
    c = CellField([1, 2, 3], grid)
    cz_top = T(0.2)
    cz_bottom = T(0.11)
    top_bc = GradientBoundaryCondition(cz_top)
    bottom_bc = GradientBoundaryCondition(cz_bottom)
    fieldbcs = FieldBoundaryConditions(bottom_bc, top_bc)

    fill_ghost_cells!(c, 0, 0, model, fieldbcs)

    ∂z(c, 1) ≈ cz_bottom && ∂z(c, grid.N+1) ≈ cz_top
end

function ghost_cell_flux(T)
    model = DummyModel()
    grid = UniformGrid(T, 3, 3.0)
    c = CellField([1, 2, 3], grid)
    Fc_top = T(0.2)
    Fc_bottom = T(0.11)
    top_bc = FluxBoundaryCondition(Fc_top)
    bottom_bc = FluxBoundaryCondition(Fc_bottom)
    fieldbcs = FieldBoundaryConditions(bottom_bc, top_bc)

    κ = T(0.47)
    fill_ghost_cells!(c, κ, κ, model, fieldbcs)

    (OceanTurb.flux(0, κ, c, 1) ≈ Fc_bottom
        && OceanTurb.flux(0, κ, c, grid.N+1) ≈ Fc_top)
end

function absolute_error_is_correct(T)
    grid = UniformGrid(T, 3, 3)
    c = CellField([0, 0, 0], grid)
    d = CellField([2, 2, 2], grid)
    absolute_error(c, d) == 2
end

function absolute_error_fine_to_coarse()
    fine_grid = UniformGrid(N=4, H=1)
    coarse_grid = UniformGrid(N=2, H=1)
    fine_field = CellField([2, 4, 6, 8], fine_grid)
    coarse_field = CellField([0, 1], coarse_grid)
    return absolute_error(coarse_field, fine_field) == sqrt(45/2)
end

function absolute_error_coarse_to_fine()
    fine_grid = UniformGrid(N=4, H=1)
    coarse_grid = UniformGrid(N=2, H=1)
    fine_field = CellField([2, 4, 6, 8], fine_grid)
    coarse_field = CellField([0, 1], coarse_grid)
    return absolute_error(fine_field, coarse_field) == sqrt(47/2)
end

function relative_error_is_correct(T)
    grid = UniformGrid(T, 3, 1)
    c = CellField([0, 0, 0], grid)
    d = CellField([2, 2, 2], grid)
    relative_error(c, d) == 1
end

@testset "Fields" begin
    N, H = 3, 4.2
    for T in (Float64,)
        @test cell_field_construction(T, N, H)
        @test face_field_construction(T, N, H)
        @test cell_∂z(T)
        @test face_∂z(T)
        @test cell_plus(T)
        @test cell_times(T)
        for loc in (Face, Cell)
            @test set_scalar_field(loc, T)
            @test set_array_field(loc, T)
            @test set_function_field(loc, T)
        end
        @test integral_is_correct(T)
        @test integral_range(T, 10, 10, -5, -3)
        @test integral_range(T, 10, 10, -5.1, -4.5)
        @test integral_range(T, 17, 3, -2.5, -1.3)
        @test integral_range(T, 3, 30, -2.5, -1.3)
        @test ghost_cell_value(T)
        @test ghost_cell_gradient(T)
        #@test ghost_cell_flux(T)
        @test absolute_error_is_correct(T)
        @test relative_error_is_correct(T)
    end
    @test absolute_error_fine_to_coarse()
    @test absolute_error_coarse_to_fine()
    @test field_indexing()
end


text:
using Compat, Compat.Test
using LightGraphs

function to_adjlist(g::Graph)
    I = Int[]
    J = Int[]
    for e in edges(g)
        push!(I,src(e))
        push!(J,dst(e))
    end
    return gen_adjlist(I,J,length(vertices(g)))
end

import ReverseDiffSparse.Coloring: acyclic_coloring, recovery_preprocess, reverse_topological_sort_by_dfs, gen_adjlist, hessian_color_preprocess, prepare_seed_matrix!, recover_from_matmat!, seed_matrix

# tests for acyclic coloring

g = Graph(10)
color, numcolors = acyclic_coloring(to_adjlist(g))
@test numcolors == 1

add_edge!(g, 2, 4)
color, numcolors = acyclic_coloring(to_adjlist(g))
@test numcolors == 2
@test color[2] != color[4]

add_edge!(g, 2, 3)
color, numcolors = acyclic_coloring(to_adjlist(g))
@test numcolors == 2
@test color[3] == color[4]

add_edge!(g, 3, 4)
color, numcolors = acyclic_coloring(to_adjlist(g))
@test numcolors == 3
recovery_preprocess(to_adjlist(g), color, numcolors, Int[])

g = Graph(3)
add_edge!(g, 1, 3)
add_edge!(g, 2, 3)
color, numcolors = acyclic_coloring(to_adjlist(g))
@test numcolors == 2

g = Graph(4)
add_edge!(g, 1, 2)
add_edge!(g, 2, 3)
add_edge!(g, 3, 4)
add_edge!(g, 4, 1)
color, numcolors = acyclic_coloring(to_adjlist(g))
@test numcolors == 3

# test our topological sort method
#=
g = Graph(6)
add_edge!(g, 1,2)
add_edge!(g, 1,3)
add_edge!(g, 1,6)
add_edge!(g, 2,4)
add_edge!(g, 2,5)
=#
vec = [3,6,2,1,4,5,1,2,2,1]
offset = [1,4,7,8,9,10,11]

v = reverse_topological_sort_by_dfs(vec, offset, 6, zeros(Int,6))
@test v[1] == [3,6,4,5,2,1]
@test v[2] == [0,1,1,2,2,1]

I,J,rinfo = hessian_color_preprocess(Set([(1,2)]),2)
num_colors = rinfo.num_colors
N = length(rinfo.color)
R = seed_matrix(rinfo)
prepare_seed_matrix!(R, rinfo)
@test I == [1,2,2]
@test J == [1,2,1]
@test R == [1.0 0.0; 0.0 1.0]
hess = [3.4 2.1; 2.1 1.3]
matmat = hess*R
V = zeros(3)
recover_from_matmat!(V, matmat, rinfo, zeros(3))
@test V == [3.4,1.3,2.1]


text:
const SymString = Union{AbstractString, Symbol}

macro map(expr...)
    T, ex =
        if length(expr) == 1
            nothing, first(expr)
        else
            first(expr), last(expr)
    end

    args = map(ex.args) do arg
        if Meta.isexpr(arg, :(=))
            key, value = arg.args
            return :( $(QuoteNode(key)) => $value )
        elseif Meta.isexpr(arg, :(:=))
            key, value = arg.args
            return :( $key => $value )
        else
            error("Invalid syntax for @map, expected expression of type `foo = bar` or `foo := bar`")
        end
    end

    result =
        if isnothing(T)
            Expr(:call, :Dict, args...)
        else
            Expr(:call, Expr(:curly, :Dict, T.args...), args...)
        end

    return esc(result)
end

macro widget(expr)
    if @capture(expr, mutable struct name_ <: super_ fields__ end)
        quote
            mutable struct $name <: Widget{Lynx.gtype($super)}
                $(fields...)
            end
            Lynx.gwidget(w::$name) = Lynx.gwidget(getfield(w, :widget))
        end |> esc
    else
        name = gensym("widget")
        quote
            begin
                $(name) = $(expr)
                for (prop, value) in pairs(props)
                    setprop!($(name), prop, value)
                end
                $(name)
            end
        end |> esc
    end
end

macro container(expr)
    name = gensym("widget")
    quote
        begin
            $(name) = $(expr)
            for (prop, value) in pairs(props)
                setprop!($(name), prop, value)
            end
            for child in children
                push!($(name), gwidget(child))
            end
            $(name)
        end
    end |> esc
end

macro secure(expr, msg)
    esc(quote
        try
            $(expr)
        catch e
            @error $(msg) exception = e
        end
    end)
end

macro unpack(expr)
    if @capture(expr, (vars__,) = object_)
        res = Expr[]
        for var in vars
            push!(res, :( $var = $(object).$var ))
        end
        return esc(quote $(res...) end)
    end
end

function Base.convert(::Type{T}, color::GdkRGBA) where {T <: Colorant}
    @unpack r, g, b, a = color
    return convert(T, RGBA(r, g, b, a))
end

function Base.convert(::Type{GdkRGBA}, color::Colorant)
    rgba = convert(RGBA, color)
    @unpack r, g, b, alpha = rgba
    return GdkRGBA(r, g, b, alpha)
end

function middle(r::AbstractRange)
    N = length(r)
    return isodd(N) ? r[N ÷ 2 + 1] : r[N ÷ 2]
end

const Center   = Gtk.GtkAlign.CENTER
const Fill     = Gtk.GtkAlign.FILL
const Start    = Gtk.GtkAlign.START
const End      = Gtk.GtkAlign.END
const Baseline = Gtk.GtkAlign.BASELINE

getprop(w::GtkWidget, prop::SymString, T::DataType) = get_gtk_property(w, prop, T)
setprop!(w::GtkWidget, prop::SymString, value::T) where {T} = set_gtk_property!(w, prop, value)
onevent(callback::Function, event::SymString, w::GtkWidget) = signal_connect(callback, w, event)
disconnect(w::GtkWidget, id) = signal_handler_disconnect(w, id)

function setprop!(widget::GtkWidget, prop::SymString, signal::Observable)
    on(signal) do value
        set_gtk_property!(widget, prop, value)
    end
end

abstract type Iterable end

function Base.iterate(itr::T, state = 1) where T <: Iterable
    if state <= fieldcount(T)
        return getfield(itr, state), state + 1
    else
        return nothing
    end
end

var"@width"(::LineNumberNode, ::Module, widget)  = :( width($( esc(widget) )) )
var"@height"(::LineNumberNode, ::Module, widget) = :( height($( esc(widget) )) )

const _ref_dict = IdDict{Any, Any}()

"""
    gcpreserve(widget::GtkWidget, obj)
Preserve `obj` until `widget` has been destroyed.
"""
function gcpreserve(widget::Union{GtkWidget,GtkCanvas}, obj)
    _ref_dict[obj] = obj
    signal_connect(widget, :destroy) do w
        delete!(_ref_dict, obj)
    end
end

macro gcpreserve(expr)
    if @capture(expr, new_(widget_) | new_(widget_, args__))
        quote
            begin
                this = $(expr)
                gcpreserve($widget, this)
                this
            end
        end |> esc
    end
end

function splitv(f::Function, array::AbstractVector{T}) where {T}
    i = 1
    result = Vector{T}[ T[] ]
    foreach(array) do element
        if f(element)
            i += 1
            push!(result, T[])
        else
            push!(result[i], element)
        end
    end
    return filter(!isempty, result)
end

splitv(x::T, array::AbstractVector{T}) where {T} = splitv(==(x), array)