text:
function get_affine_map(ds::ArchGDAL.IDataset)
    # ArchGDAL fails hard on datasets without
    # an affinemap. GDAL documents that on fail
    # a default affinemap should be returned.
    local gt
    try
        gt = ArchGDAL.getgeotransform(ds)
    catch y
        @warn y.msg
        gt = [0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
    end
    geotransform_to_affine(gt)
end

function geotransform_to_affine(gt::SVector{6,<:AbstractFloat})
    # See https://lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html
    # for an explanation of the geotransform format
    AffineMap(SMatrix{2,2}([gt[2] gt[3]; gt[5] gt[6]]), SVector{2}([gt[1], gt[4]]))
end
geotransform_to_affine(A::Vector{<:AbstractFloat}) = geotransform_to_affine(SVector{6}(A))

function affine_to_geotransform(am::AffineMap)
    l = am.linear
    t = am.translation
    (length(l) == 4 && length(t) == 2) || error("AffineMap has wrong dimensions.")
    [t[1], l[1], l[3], t[2], l[2], l[4]]
end

"""Check wether the AffineMap of a GeoArray contains rotations."""
function is_rotated(ga::GeoArray)
    ga.f.linear[2] != 0.0 || ga.f.linear[3] != 0.0
end

function bbox(ga::GeoArray)
    i, j, _ = size(ga)
    if is_rotated(ga)
        ax, ay = ga.f(SVector(0, 0))
        bx, by = ga.f(SVector(i, 0))
        cx, cy = ga.f(SVector(0, j))
        dx, dy = ga.f(SVector(i, j))
        (min_x=min(ax, bx, cx, dx), min_y=min(ay, by, cy, dy), max_x=max(ax, bx, cx, dx), max_y=max(ay, by, cy, dy))
    else
        ax, ay = ga.f(SVector(0, 0))
        bx, by = ga.f(SVector(i, j))
        (min_x=min(ax, bx), min_y=min(ay, by), max_x=max(ax, bx), max_y=max(ay, by))
    end
end

function unitrange_to_affine(x::StepRangeLen, y::StepRangeLen)
    δx, δy = float(step(x)), float(step(y))
    AffineMap(
        SMatrix{2,2}(δx, 0.0, 0.0, δy),
        SVector(x[1] - δx / 2, y[1] - δy / 2)
    )
end

function bbox_to_affine(size::Tuple{Integer,Integer}, bbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
    AffineMap(
        SMatrix{2,2}(float(bbox.max_x - bbox.min_x) / size[1], 0, 0, float(bbox.max_y - bbox.min_y) / size[2]),
        SVector(float(bbox.min_x), float(bbox.min_y))
    )
end

"""Set geotransform of `GeoArray` by specifying a bounding box.
Note that this only can result in a non-rotated or skewed `GeoArray`."""
function bbox!(ga::GeoArray, bbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
    ga.f = bbox_to_affine(size(ga)[1:2], bbox)
    ga
end

"""Generate bounding boxes for GeoArray cells."""
function bboxes(ga::GeoArray)
    c = coords(ga, Vertex())
    m, n = size(c)
    cellbounds = Matrix{NamedTuple}(undef, (m - 1, n - 1))
    for j in 1:n-1, i in 1:m-1
        v = c[i:i+1, j:j+1]
        minx, maxx = extrema(first.(v))
        miny, maxy = extrema(last.(v))
        cellbounds[i, j] = (min_x=minx, max_x=maxx, min_y=miny, max_y=maxy)
    end
    cellbounds
end

# Extend CoordinateTransformations
CoordinateTransformations.compose(ga::GeoArray, t2::AffineMap) = CoordinateTransformations.compose(ga.f, t2)
CoordinateTransformations.compose(ga::GeoArray, t2::LinearMap) = CoordinateTransformations.compose(ga.f, t2)
CoordinateTransformations.compose(ga::GeoArray, t2::Translation) = CoordinateTransformations.compose(ga.f, t2)

"""Transform an GeoArray by applying a Transformation."""
function compose!(ga::GeoArray, t2::AffineMap)
    ga.f = compose(ga, t2)
end
function compose!(ga::GeoArray, t2::LinearMap)
    ga.f = compose(ga, t2)
end
function compose!(ga::GeoArray, t2::Translation)
    ga.f = compose(ga, t2)
end

"""Function to flip GeoArray upside down to adjust to GDAL ecosystem."""
function flipud!(ga::GeoArray)
    # Flip data upside down
    ga.A = reverse(ga.A, dims=2)

    # Find new corner coordinates
    ux, uy = ga.f(SVector{2}([0, size(ga)[2]]))

    # Define y mirror and compose
    lm = LinearMap(SMatrix{2,2}([1.0 0.0; 0.0 -1.0]))
    am = compose(ga.f, lm)  # AffineMap
    translate = SVector{2}([ux, uy])
    f = AffineMap(am.linear, translate)

    ga.f = f
    ga
end

"""Check bbox overlapping

Return true if two bboxes overlap.
"""
function bbox_overlap(
    bbox_a::NamedTuple{(:min_x, :min_y, :max_x, :max_y)},
    bbox_b::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})

    # Check if bboxes are valid
    if (bbox_a.min_x >= bbox_a.max_x) ||
       (bbox_a.min_y >= bbox_a.max_y)
        error("Invalid bbox (min >= max)", bbox_a)
    end
    if (bbox_b.min_x >= bbox_b.max_x) ||
       (bbox_b.min_y >= bbox_b.max_y)
        error("Invalid bbox (min >= max)", bbox_b)
    end

    # Check if bboxes overlap
    if (bbox_a.max_x < bbox_b.min_x) ||
       (bbox_a.max_y < bbox_b.min_y) ||
       (bbox_a.min_x > bbox_b.max_x) ||
       (bbox_a.min_y > bbox_b.max_y)
        return false
    else
        return true
    end
end

"""
    function crop(ga::GeoArray, cbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y),Tuple{Float64,Float64,Float64,Float64}})

Crop input GeoArray by coordinates (box or another GeoArray). If the coordinates range is larger than the GeoArray, only the overlapping part is given in the result.
"""
function crop(ga::GeoArray, cbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
    # Check if ga and cbox overlap
    if !bbox_overlap(bbox(ga), cbox)
        error("GeoArray and crop box do not overlap")
    end
    # Check extent and get bbox indices
    ga_x, ga_y, = size(ga)
    i_min_x, i_min_y = indices(ga, [cbox.min_x, cbox.max_y])
    i_max_x, i_max_y = indices(ga, [cbox.max_x, cbox.min_y])

    # Determine indices for crop area
    i_min_x = max(i_min_x, 1)
    i_max_x = min(i_max_x, ga_x)
    i_min_y = max(i_min_y, 1)
    i_max_y = min(i_max_y, ga_y)

    # Subset and return GeoArray
    return ga[i_min_x:i_max_x, i_min_y:i_max_y, :]
end

function crop(ga::GeoArray, cga::GeoArray)
    # Check if ga and crop ga CRS are the same
    if ga.crs != cga.crs
        error("GeoArrays have different CRS")
    end

    # Crop
    return crop(ga, bbox(cga))
end

"""
    sample!(ga::GeoArray, ga2::GeoArray)

Sample values from ga2 to ga.
"""
function sample!(ga::GeoArray, ga2::GeoArray)
    if ga.crs != ga2.crs
        error("GeoArrays have different CRS")
    end
    wo, ho, zo = size(ga)
    w, h = size(ga2)[1:2]

    # Compose AffineMap that translates from logical coordinates
    # in `ga` to logical coordinates in `ga2`
    x = inv(ga2.f) ∘ ga.f

    for io in 1:wo, jo in 1:ho
        i, j = round.(Int, x((io, jo)) .+ 0.5)
        if (1 <= i <= w) && (1 <= j <= h)
            # Loop over bands
            for z in 1:zo
                ga[io, jo, z] = ga2[i, j, z]
            end
        end
    end

    ga
end

function _sizeof(ga::GeoArray, bbox)
    x = bbox.max_x - bbox.min_x
    y = bbox.max_y - bbox.min_y
    abs.(round.(Int, (x / ga.f.linear[1], y / ga.f.linear[4], size(ga)[3])))
end

"""
    straighten(ga::GeoArray)

Straighten a rotated GeoArray, i.e. let its AffineMap only scale the coordinates.
"""
function straighten(ga)
    is_rotated(ga) || return ga

    dtype = eltype(ga)
    if !(isa(dtype, Union) && dtype.a == Missing)
        dtype = Union{Missing,dtype}
    end

    A = Array{dtype}(undef, _sizeof(ga, bbox(ga)))
    fill!(A, missing)
    gar = GeoArray(A)

    bbox!(gar, bbox(ga))
    gar.crs = ga.crs
    GeoArrays.sample!(gar, ga)
    gar
end


text:
# Autogenerated wrapper script for HarfBuzz_jll for x86_64-unknown-freebsd
export libharfbuzz, libharfbuzz_gobject, libharfbuzz_icu, libharfbuzz_subset

using Glib_jll
using FreeType2_jll
using Graphite2_jll
using Libffi_jll
using Gettext_jll
using Fontconfig_jll
using Cairo_jll
using ICU_jll
JLLWrappers.@generate_wrapper_header("HarfBuzz")
JLLWrappers.@declare_library_product(libharfbuzz, "libharfbuzz.so.0")
JLLWrappers.@declare_library_product(libharfbuzz_gobject, "libharfbuzz-gobject.so.0")
JLLWrappers.@declare_library_product(libharfbuzz_icu, "libharfbuzz-icu.so.0")
JLLWrappers.@declare_library_product(libharfbuzz_subset, "libharfbuzz-subset.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Glib_jll, FreeType2_jll, Graphite2_jll, Libffi_jll, Gettext_jll, Fontconfig_jll, Cairo_jll, ICU_jll)
    JLLWrappers.@init_library_product(
        libharfbuzz,
        "lib/libharfbuzz.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libharfbuzz_gobject,
        "lib/libharfbuzz-gobject.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libharfbuzz_icu,
        "lib/libharfbuzz-icu.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libharfbuzz_subset,
        "lib/libharfbuzz-subset.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# linear constraints

MOI.supports_constraint(o::Optimizer, ::Type{SAF}, ::Type{<:BOUNDS}) = true

function MOI.add_constraint(o::Optimizer, func::F, set::S) where {F <: SAF, S <: BOUNDS}
    if func.constant != 0.0
        throw(MOI.ScalarFunctionConstantNotZero{Float64, F, S}(func.constant))
    end

    allow_modification(o)

    varrefs = [VarRef(t.variable.value) for t in func.terms]
    coefs = [t.coefficient for t in func.terms]

    lhs, rhs = bounds(set)
    lhs = lhs === nothing ? -SCIPinfinity(o) : lhs
    rhs = rhs === nothing ?  SCIPinfinity(o) : rhs

    cr = add_linear_constraint(o.inner, varrefs, coefs, lhs, rhs)
    ci = CI{F, S}(cr.val)
    register!(o, ci)
    register!(o, cons(o, ci), cr)
    return ci
end

function MOI.delete(o::Optimizer, ci::CI{<:SAF, S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    allow_modification(o)
    delete!(o.constypes[SAF, S], ConsRef(ci.value))
    delete!(o.reference, cons(o, ci))
    delete(o.inner, ConsRef(ci.value))
    return nothing
end

function MOI.set(o::SCIP.Optimizer, ::MOI.ConstraintSet, ci::CI{<:SAF,S}, set::S) where {S <: BOUNDS}
    allow_modification(o)

    lhs, rhs = bounds(set)
    lhs = lhs === nothing ? -SCIPinfinity(o) : lhs
    rhs = rhs === nothing ?  SCIPinfinity(o) : rhs

    @SCIP_CALL SCIPchgLhsLinear(o, cons(o, ci), lhs)
    @SCIP_CALL SCIPchgRhsLinear(o, cons(o, ci), rhs)

    return nothing
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{<:SAF, S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    c = cons(o, ci)
    nvars::Int = SCIPgetNVarsLinear(o, c)
    vars = unsafe_wrap(Array{Ptr{SCIP_VAR}}, SCIPgetVarsLinear(o, c), nvars)
    vals = unsafe_wrap(Array{Float64}, SCIPgetValsLinear(o, c), nvars)

    terms = [AFF_TERM(vals[i], VI(ref(o, vars[i]).val)) for i=1:nvars]
    # can not identify constant anymore (is merged with lhs,rhs)
    return SAF(terms, 0.0)
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{<:SAF, S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    lhs = SCIPgetLhsLinear(o, cons(o, ci))
    rhs = SCIPgetRhsLinear(o, cons(o, ci))
    return from_bounds(S, lhs, rhs)
end

function MOI.modify(o::Optimizer, ci::CI{<:SAF, <:BOUNDS},
                    change::MOI.ScalarCoefficientChange{Float64})
    allow_modification(o)
    @SCIP_CALL SCIPchgCoefLinear(o, cons(o, ci), var(o, change.variable), change.new_coefficient)
    return nothing
end


text:
struct S end
struct G end

"""
Optimization is performed on this model
"""
mutable struct OptimModel{T}
	ntg::Int64
	nt::Int64
	nts::Int64
	nr::Int64
	obs::Conv.P_conv{T,2,2,1} # observed convolutional model
	cal::Conv.P_conv{T,2,2,1} # calculated convolutional model
	dg::Array{T,2} # store gradients
	ds::Vector{T} # store gradients
	ddcal::Array{T,2}
	gfilt::Array{T,2}
	sfilt::Array{T,1}
	dvec::Vector{T}
end

function OptimModel(ntg, nt, nr, nts, T; fftwflag=FFTW.MEASURE, slags=nothing, dlags=nothing, glags=nothing)

	obs=Conv.P_conv(T,ssize=[nts], dsize=[nt,nr], gsize=[ntg,nr],
	slags=slags,
	glags=glags,
	dlags=dlags,
	fftwflag=fftwflag,
	)
	cal=deepcopy(obs)

	# initial values are random
	ds=zero(cal.s)
	dg=zero(cal.g)
	ddcal=zero(cal.d)
	gfilt=zero(cal.g)
	sfilt=zero(cal.s)
	dvec=zeros(nt*nr) # if data needs to be stored as vector

	return OptimModel{T}(ntg, nt, nts, nr, obs, cal, dg, ds, ddcal, gfilt, sfilt, dvec)
end

function replace!(pa::OptimModel, x, fieldmod::Symbol, field::Symbol)
	pao=getfield(pa, fieldmod)
	paoo=getfield(pao, field)
	if(length(paoo)==length(x))
		for i in eachindex(paoo)
			paoo[i]=x[i]
		end
	else
		error("error dimension replacing")
	end
	# perform modelling :<)
	Conv.mod!(pao, Conv.D())
end

"""
Store actual
"""
mutable struct ObsModel{T}
	ntg::Int64
	nt::Int64
	nts::Int64
	nr::Int64
	g::Array{T,2}
	s::Vector{T}
	d::Array{T,2}
end

function ObsModel(ntg, nt, nr, nts, T::DataType;
	       d=nothing,
	       g=nothing,
	       s=nothing,
	       slags=[nts-1, 0],
	       glags=nothing,
	       dlags=nothing
	       )
	(s===nothing) && (s=zeros(nts))
	(g===nothing) && (g=zeros(ntg,nr))
	if(d===nothing)
		(iszero(g) || iszero(s)) && error("need gobs and sobs")
		obstemp=Conv.P_conv(T,ssize=[nts], dsize=[nt,nr], gsize=[ntg,nr], slags=slags, dlags=dlags, glags=glags)
		copyto!(obstemp.g, g)
		copyto!(obstemp.s, s)
		Conv.mod!(obstemp, Conv.D()) # do a convolution to model data
		d=obstemp.d
	end

	return ObsModel{T}(ntg, nt, nts, nr, g, s, d)
end

"""
Inversion variables
"""
mutable struct X{T}
	x::Vector{T}
	last_x::Vector{T}
	lower_x::Vector{T}
	upper_x::Vector{T}
	precon::Vector{T}
	preconI::Vector{T}
	weights::Vector{T}
end

function X(n,T)
	return X{T}(zeros(n), randn(n), fill(-Inf,n), fill(Inf,n), ones(n), ones(n), zeros(n))
end

mutable struct P_bandpass{T<:Real}
	nt::Int
	x::Vector{T}
	y::Vector{T}
	dfreq::Vector{Complex{T}}
	dfftp::FFTW.rFFTWPlan
	difftp::FFTW.Plan
	filter::Vector{T}
end

function P_bandpass(T=Float64;
		    fmin=0.0, fmax=0.5, nt)

	l=nt
	x=zeros(T,l)
	y=zeros(T,l)
	nrfft=div(l,2)+1
	dfftp=plan_rfft(zeros(T, l),[1])
	difftp=plan_irfft(complex.(zeros(T, nrfft)),l,[1])

	filter=ones(T,nrfft)
	ω=FFTW.rfftfreq(l)
	for i in eachindex(filter)
		filter[i]=(((ω[i]-fmin)*(fmax-ω[i])) ≥ 0.0) ? 1.0 : 0.0
	end

	dfreq=complex.(zeros(T,nrfft))
	return P_bandpass(nt,x,y,dfreq,dfftp, difftp, filter)
end


text:
using Spells
using Spells: Context, If, Loop, emit, Func
using SPIRV
using Test

spirv_file(filename) = joinpath(@__DIR__, "resources", filename * ".spv")

@testset "Spells.jl" begin
    ir = IR(SPIRV.Module(spirv_file("unicolor.vert")))

    ctx = Context()
    f = Func()
end


text:
@testset "GridTools" begin
    nt = length(T_all)
    n = nt * nb
    @unpack xgeo = p
    x = xgeo * ones(n)
    @testset "convert to 3D array" for xᵢ in state_to_tracers(x, nb, nt)
        xᵢ3D = rearrange_into_3Darray(xᵢ, grd)
        @test size(xᵢ3D) == size(grd)
    end
    @testset "convert to 3D array" for xᵢ in unpack_tracers(x, grd)
        xᵢ3D = rearrange_into_3Darray(xᵢ, grd)
        @test size(xᵢ3D) == size(grd)
    end
end


text:
function restaurant(p::Player)
    print_dict(p, "restaurant_intro")
    while true
        user_input = readline(p.socket)
        if (user_input == "c")
            counter(p)
            print_dict(p, "restaurant_intro")
        elseif (user_input == "l")
            return
        elseif (user_input == "🧀")
            searchCheese(p)
        end
    end
end

function checkToken(p::Player, code::String)
    token_json = JSON.parse(p.token, inttype=UInt8)
    enc_msg = UInt8.(token_json[1])

    cryptomaterial = generate_cryptomaterial(p)
    cryptomaterial[2] = UInt8.(token_json[2])

    msg = AESCBC(enc_msg, cryptomaterial[1], cryptomaterial[2], false)

    x = abs(p.dimension) % length(code)+1
    code = string(code[x+1:length(code)], code[1:x])

    msg = String(msg)
    msg = rstrip(msg, '\x00')

    return msg == code

end

function checkTokenExistence(p::Player)
    return !(p.token == "")
end

function counter(p::Player)
    print_dict(p, "restaurant_counter")
    print_dict(p, "restaurant_menu")
    while true
        user_input = lowercase(readline(p.socket))
        if (user_input == "l")
            return
        elseif (user_input == "f")
            print_dict(p, "restaurant_flies_and_moths")
            return
        elseif (user_input == "r")
            print_dict(p, "restaurant_robot_choice")
            return
        elseif (user_input == "c")
            if !checkTokenExistence(p)
                print_dict(p, "restaurant_no_token")
                continue
            end
            if !checkToken(p, "CASINO_ROYALE")
                print_dict(p, "restaurant_token_fail")
                continue
            end
            print_dict(p, "restaurant_casino_royale")
            p.diarrhea = true
            break
        elseif (user_input == "d")
            print_dict(p, "restaurant_dips_chips")
        else
            print_dict(p, "restaurant_invalid_input")
        end
    end
end

function searchCheese(p::Player)
    if !(isfile("data/.leaflet"))
        print_dict(p, "restaurant_cheese_fail")
    else
        print_dict(p, "restaurant_cheese")
        f = open_file_try("data/.leaflet", "r")
        write(p.socket, "$(read(f, String))\n")
        close(f)
    end
end


text:
# This file exemplifies the workflow from data input to optimization result generation

using TimeSeriesClustering

# load data
data_path=normpath(joinpath(dirname(@__FILE__),"..","data","TS_GER_18"))
ts_input_data = load_timeseries_data(data_path; T=24, years=[2015])

 # define simple extreme days of interest
 ev1 = SimpleExtremeValueDescr("wind-dena42","max","absolute")
 ev2 = SimpleExtremeValueDescr("solar-dena42","min","integral")
 ev3 = SimpleExtremeValueDescr("el_demand-dena21","max","absolute")
 ev = [ev1, ev2, ev3]
 # simple extreme day selection
 ts_input_data_mod,extr_vals,extr_idcs = simple_extr_val_sel(ts_input_data,ev;rep_mod_method="feasibility")

 # run clustering
ts_clust_res = run_clust(ts_input_data_mod;method="kmeans",representation="centroid",n_init=10,n_clust=5) # default k-means

# representation modification
ts_clust_extr = representation_modification(extr_vals,ts_clust_res.clust_data)


text:
# Adapted from julia/special/trig.jl

if VERSION < v"0.7.0-beta2.199"
    function Base.DomainError(x)
        throw(DomainError())
    end
end

function sintau(x::Real)
    if !isfinite(x)
        isnan(x) && return x
        throw(DomainError(x))
    end

    rx = copysign(float(rem(x,1)),x)
    arx = abs(rx)

    if arx <= oftype(rx,1/8)
        return sin(tau*rx)
    elseif arx <= oftype(rx,3/8)
        arx = oftype(rx,1/4) - arx
        return copysign(cos(tau*arx),rx)
    elseif arx < oftype(rx,5/8)
        rx = (oftype(rx,1/2) - arx)*sign(rx)
        return sin(tau*rx)
    elseif arx <= oftype(rx,7/8)
        arx = oftype(rx,3/4) - arx
        return -copysign(cos(tau*arx),rx)
    else
        rx = rx - copysign(one(rx),rx)
        return sin(tau*rx)
    end
end

function costau(x::Real)
    if !isfinite(x)
        isnan(x) && return x
        throw(DomainError(x))
    end

    rx = abs(float(rem(x,1)))

    if rx <= oftype(rx,1/8)
        return cos(tau*rx)
    elseif rx < oftype(rx,3/8)
        rx = oftype(rx,1/4) - rx
        return sin(tau*rx)
    elseif rx <= oftype(rx,5/8)
        rx = oftype(rx,1/2) - rx
        return -cos(tau*rx)
    elseif rx < oftype(rx,7/8)
        rx = rx - oftype(rx,3/4)
        return sin(tau*rx)
    else
        rx = one(rx) - rx
        return cos(tau*rx)
    end
end

sintau(x::Integer) = zero(x)
costau(x::Integer) = one(x)

# Implementations of complex sintau/costau are adapted from julia/base/special/trig.jl.
function sintau(z::Complex{T}) where T
    F = float(T)
    zr, zi = reim(z)
    if isinteger(zr)
        Complex(zero(F), sinh(tau * zi))
    elseif !isfinite(zr)
        if zi == 0 || isinf(zi)
            Complex(F(NaN), F(zi))
        else
            Complex(F(NaN), F(NaN))
        end
    else
        tauzi = tau * zi
        Complex(sintau(zr) * cosh(tauzi), costau(zr) * sinh(tauzi))
    end
end

function costau(z::Complex{T}) where T
    F = float(T)
    zr, zi = reim(z)
    if isinteger(zr)
        s = copysign(zero(F), zr)
        Complex(cosh(tau * zi), isnan(zi) ? s : -flipsign(s,zi))
    elseif !isfinite(zr)
        if zi == 0
            Complex(F(NaN), isnan(zr) ? zero(F) : -flipsign(F(zi), zr))
        elseif isinf(zi)
            Complex(F(Inf), F(NaN))
        else
            Complex(F(NaN), F(NaN))
        end
    else
        tauzi = tau * zi
        Complex(costau(zr) * cosh(tauzi), -sintau(zr) * sinh(tauzi))
    end
end

const sinτ = sintau
const cosτ = costau


text:
using Cxx
CCompiler = Cxx.new_clang_instance(false, true)

cxx"""
#include "clang/AST/RecordLayout.h"
"""

# Iteration for ArrayRef
import Base: start, next, length, done
const ArrayRef = cxxt"llvm::ArrayRef<$T>" where T
start(AR::ArrayRef) = 0
function next(AR::cxxt"llvm::ArrayRef<$T>", i) where T
    (icxx"""
        // Force a copy, otherwise we'll retain reference semantics in julia
        // which is not what people expect.
        $T element = ($AR)[$i];
        return element;
    """, i+1)
end
length(AR::ArrayRef) = icxx"$AR.size();"
done(AR::ArrayRef, i) = i >= length(AR)

# Working with clang::IdentifierInfo
Base.String(II::pcpp"clang::IdentifierInfo") = unsafe_string(icxx"$II->getName().str();")

# Printing for tokens
function Base.String(Tok::cxxt"clang::Token")
    II = icxx"$Tok.getIdentifierInfo();"
    @assert II != C_NULL
    String(II)
end
tok_is_identifier(Tok) = icxx"$Tok.is(clang::tok::identifier);"
tok_is_comma(Tok) = icxx"$Tok.is(clang::tok::comma);"
tok_is_numeric(Tok) = icxx"$Tok.is(clang::tok::numeric_constant);"
getSpelling(PP, Tok) = unsafe_string(icxx"$PP->getSpelling($Tok);")
function Base.show(io::IO, Tok::cxxt"clang::Token")
    print(io, unsafe_string(icxx"clang::tok::getTokenName($Tok.getKind());"))
    if haskey(io, :preprocessor)
        PP = get(io, :preprocessor, nothing)
        print(io, " '", getSpelling(PP, Tok), "'")
    elseif tok_is_identifier(Tok)
        print(io, " '", String(Tok), "'")
    end
end

# Iteration over IdentifierTable
start(tab::rcpp"clang::IdentifierTable") = icxx"$tab.begin();"
next(tab::rcpp"clang::IdentifierTable", it) = (icxx"$it->second;", icxx"++$it;")
done(tab::rcpp"clang::IdentifierTable", it) = icxx"$it == $tab.end();"
length(tab::rcpp"clang::IdentifierTable") = icxx"$tab.size();"

function AdvanceIfEof(P)
  icxx"""
  if ($P->getPreprocessor().isIncrementalProcessingEnabled() &&
    $P->getCurToken().is(clang::tok::eof))
      $P->ConsumeToken();
  """
end

# Parsing
function EnterTokenStream(PP::pcpp"clang::Preprocessor", tokens::Vector{cxxt"clang::Token"})
  # Vector memory layout is incompatible, convert to clang::Token**
  toks = typeof(tokens[1].data)[x.data for x in tokens]
  icxx"$PP->EnterTokenStream(llvm::ArrayRef<clang::Token>{
    (clang::Token*)$(pointer(toks)),
    (size_t)$(length(toks))
  },false);"
end
function ParseTypeName(P::pcpp"clang::Parser")
  AdvanceIfEof(P)
  res = icxx"$P->ParseTypeName(nullptr, clang::Declarator::TypeNameContext);"
  !icxx"$res.isUsable();" && error("Parsing failed")
  Cxx.QualType(icxx"clang::Sema::GetTypeFromParser($res.get());")
end

PP = icxx"&$(Cxx.active_instances[2].CI)->getPreprocessor();"
P  = Cxx.active_instances[2].Parser

Cxx.cxxinclude(CCompiler, "linux/termios.h")
Cxx.cxxinclude(CCompiler, "linux/usbdevice_fs.h")
map(String, Iterators.filter(II->icxx"$II->hasMacroDefinition();", icxx"$PP->getIdentifierTable();"))

Base.String(F::pcpp"clang::FieldDecl") = unsafe_string(icxx"$F->getName();")
# C structure to julia array of fields
function inspectStruct(CC, S)
    CC = Cxx.instance(CC)
    ASTCtx = icxx"&$(CC.CI)->getASTContext();"
    fields = Any[]
    icxx"""
        $S->dump();
        auto &ARL = $ASTCtx->getASTRecordLayout($S);
        for (auto field : ($S)->fields()) {
          unsigned i = field->getFieldIndex();
          field->dump();
          if (field->getType()->isUnionType())
            continue;
          if (field->getType()->isArrayType())
            continue;
          $:(begin
            QT = Cxx.QualType(icxx"return field->getType();")
            push!(fields, (
              String(icxx"return field;"),
              Cxx.juliatype(QT),
              icxx"return $ASTCtx->toCharUnitsFromBits(ARL.getFieldOffset(i)).getQuantity();"
            ))
          nothing
        end);
        }
    """
    fields
end

function desugar(C, QT)
  C = Cxx.instance(C)
  Cxx.QualType(icxx"$QT.getDesugaredType($(C.CI)->getASTContext());")
end
isPointerType(QT) = icxx"$QT->isPointerType();"
isRecordType(QT) = icxx"$QT->isRecordType();"

println("Hello World")

getMacroInfo(PP, II) = icxx"$PP->getMacroInfo($II);"
tokens(x) = icxx"$x->tokens();"

ioctl_defs = filter(map(II->(II,collect(tokens(getMacroInfo(PP, II)))),
  Iterators.filter(II->icxx"$II->hasMacroDefinition();", icxx"$PP->getIdentifierTable();"))) do x
      II, tokens = x
      isempty(tokens) && return false
      tok_is_identifier(tokens[1]) && String(tokens[1]) in ["_IO","_IOR","_IOW","_IOWR"]
  end;

function argumenttype_for_QT(QT)
  QT = desugar(CCompiler, QT)
  argumenttype = nothing
  # skip pointer/interger types for now
  if !isPointerType(QT) && isRecordType(QT)
      RD = icxx"((clang::RecordType*)$(Cxx.extractTypePtr(QT)))->getDecl();"
      if !icxx"$RD->isCompleteDefinition();"
          warn(string("RecordDecl ",unsafe_string(icxx"$RD->getName();")," used in ioctl but forward defined"))
      else
          argumenttype = inspectStruct(CCompiler, RD)
      end
  end
  argumenttype
end

ioctls = map(ioctl_defs) do x
  II, tokens = x
  # Processes e.g.
  # _IOR('U', 19, struct usbdevfs_hub_portinfo)
  # identifier '_IOR', l_paren, char_constant, comma, numeric_constant, comma, struct, identifier 'usbdevfs_hub_portinfo', r_paren
  category = getSpelling(PP,tokens[3])[2]
  code = parse(Int, getSpelling(PP,tokens[5]))
  argumenttype = nothing
  if String(tokens[1]) != "_IO"
      subt = tokens[findlast(tok_is_comma, tokens)+1:end-1]
      EnterTokenStream(PP, Vector{cxxt"clang::Token"}(subt))
      QT = ParseTypeName(P)
      argumenttype = argumenttype_for_QT(QT)
  end
  (String(II), category, code, 0, argumenttype)
end

getIdentifierInfo(PP, name) = icxx"$PP->getIdentifierInfo($name);"

# Some ioctls are not defined with the _IOR macros. Process them manually
function addLegacyIoctls(names)
  for name in names
    argumenttype = nothing
    if isa(name, Tuple)
        name, structname = name
        Cxx.EnterBuffer(Cxx.instance(CCompiler), structname)
        argumenttype = argumenttype_for_QT(ParseTypeName(P))
    end
    toks = collect(tokens(getMacroInfo(PP, getIdentifierInfo(PP, name))))
    (length(toks) == 1 && tok_is_numeric(toks[])) || continue
    val = parse(Int,getSpelling(PP,toks[]))
    push!(ioctls,(name, Char((val & 0xFF00)>>8), val&0xFF, 0, argumenttype))
  end
end

addLegacyIoctls([
 ("TCGETS", "struct termios"),
 ("TIOCGWINSZ", "struct winsize"),
 "TCSETS", "TCSETSW", "TCSETSF", "TCGETA", "TCSETA", "TCSETAW",
 "TCSETAF", "TCSBRK", "TCXONC", "TCFLSH", "TIOCEXCL", "TIOCNXCL", "TIOCSCTTY",
 "TIOCGPGRP", "TIOCSPGRP", "TIOCOUTQ", "TIOCSTI", "TIOCSWINSZ",
 "TIOCMGET", "TIOCMBIS", "TIOCMBIC", "TIOCMSET", "TIOCGSOFTCAR", "TIOCSSOFTCAR",
 "FIONREAD", "TIOCINQ", "TIOCLINUX", "TIOCCONS", "TIOCGSERIAL", "TIOCSSERIAL",
 "TIOCPKT", "FIONBIO", "TIOCNOTTY", "TIOCSETD", "TIOCGETD", "TCSBRKP",
 "TIOCSBRK", "TIOCCBRK", "TIOCGSID", "TIOCGRS485", "TIOCSRS485", "TCGETX",
 "TCSETX", "TCSETXF", "TCSETXW", "TIOCVHANGUP", "FIONCLEX", "FIOCLEX",
 "FIOASYNC", "TIOCSERCONFIG", "TIOCSERGWILD", "TIOCSERSWILD", "TIOCGLCKTRMIOS",
 "TIOCSLCKTRMIOS", "TIOCSERGSTRUCT", "TIOCSERGETLSR", "TIOCSERGETMULTI",
 "TIOCSERSETMULTI", "TIOCMIWAIT", "TIOCGICOUNT", "FIOQSIZE"])

open(joinpath(dirname(@__FILE__), "../generated/alltheioctls.jl"),"w") do f
    println(f,"const ioctls = Tuple{String,Char,Int64,Int64,Any}[")
    sort!(ioctls, by = x->(x[2],x[3]))
    for ioctl in ioctls
        println(f, "   ", ioctl)
    end
    println(f,"]")
end

#=
cxx"""
#include <iostream>
class IoctlCatcher : public clang::PPCallbacks {
public:
    virtual void MacroExpands(const clang::Token &MacroNameTok,
                 const clang::MacroDefinition &MD, clang::SourceRange Range,
                 const clang::MacroArgs *Args) override
    {
        std::cout << MacroNameTok.getIdentifierInfo()->getName().str() << std::endl;
    }

    virtual void MacroDefined(const clang::Token &MacroNameTok,
                              const clang::MacroDirective *MD) override
    {
        std::cout << "DEFINED " << MacroNameTok.getIdentifierInfo()->getName().str() << std::endl;
    }
};
"""
icxx"$(Cxx.active_instances[2].CI)->getPreprocessor().addPPCallbacks(
  std::unique_ptr<clang::PPCallbacks>{new IoctlCatcher});"
Cxx.cxxinclude(CCompiler, "termios.h")
Cxx.cxxinclude(CCompiler, "linux/usbdevice_fs.h")

function getMacroInfo(PP, name)
  icxx"$PP->getMacroInfo($PP->getIdentifierInfo($name));"
end
function dump_macro(PP, name)
  icxx"$PP->DumpMacro(*$(getMacroInfo(PP, name)));"
end

PP=icxx"&$(Cxx.active_instances[2].CI)->getPreprocessor();"
=#


text:
# NOTE: This file contains many parts that are copied from the file runtests.jl from the Package ParallelStencil.jl.
push!(LOAD_PATH, "../src")

using PseudoTransientStokes

function runtests()
    exename = joinpath(Sys.BINDIR, Base.julia_exename())
    testdir = pwd()
    istest(f) = endswith(f, ".jl") && startswith(basename(f), "test_")
    testfiles = sort(filter(istest, vcat([joinpath.(root, files) for (root, dirs, files) in walkdir(testdir)]...)))

    nfail = 0
    printstyled("Testing PseudoTransientStokes.jl\n"; bold=true, color=:white)

    for f in testfiles
        try
            run(`$exename -O3 --startup-file=no --check-bounds=no $(joinpath(testdir, f))`)
        catch ex
            nfail += 1
        end
    end
    return nfail
end

exit(runtests())


text:
function measure_mean!(tuning_run::Run, x::Configuration)
    next_proc   = @task chooseproc()
    @sync begin
        for i = 1:tuning_run.cost_evaluations
            @async begin
                tuning_run.cost_values[i] = remotecall_fetch(consume(next_proc),
                                                             tuning_run.cost,
                                                             x,
                                                             tuning_run.cost_arguments)
            end
        end
    end
    mean(tuning_run.cost_values)
end

function sequential_measure_mean!(tuning_run::Run, x::Configuration)
   for i = 1:tuning_run.cost_evaluations
        tuning_run.cost_values[i] = tuning_run.cost(x, tuning_run.cost_arguments)
    end
    mean(tuning_run.cost_values)
end


text:
###############################################################################
#
#  NfRelOrd
#
###############################################################################

# I don't like that we have to drag S around
# Still hoping for julia/#18466

mutable struct NfRelOrdSet{T}
  nf::RelativeExtension{T}

  function NfRelOrdSet{T}(K::RelativeExtension{T}) where {T}
    a = new(K)
    return a
  end
end

mutable struct NfRelOrd{T, S} <: Ring
  nf::RelativeExtension{T}
  basis_nf::Vector{RelativeElement{T}}
  basis_mat::Generic.Mat{T}
  basis_mat_inv::Generic.Mat{T}
  basis_pmat::PMat{T, S}
  pseudo_basis::Vector{Tuple{RelativeElement{T}, S}}

  disc_abs::NfOrdIdl # used if T == nf_elem
  disc_rel#::NfRelOrdIdl{T} # used otherwise; is a forward declaration
  parent::NfRelOrdSet{T}

  isequation_order::Bool

  ismaximal::Int                   # 0 Not known
                                   # 1 Known to be maximal
                                   # 2 Known to not be maximal

  trace_mat::Generic.Mat{T}

  inv_coeff_ideals::Vector{S}

  function NfRelOrd{T, S}(K::RelativeExtension{T}) where {T, S}
    z = new{T, S}()
    z.nf = K
    z.parent = NfRelOrdSet{T}(K)
    z.isequation_order = false
    z.ismaximal = 0
    return z
  end

  function NfRelOrd{T, S}(K::RelativeExtension{T}, M::PMat{T, S}) where {T, S}
    z = NfRelOrd{T, S}(K)
    z.nf = K
    z.parent = NfRelOrdSet{T}(K)
    z.basis_pmat = M
    z.basis_mat = M.matrix
    return z
  end

  function NfRelOrd{T, S}(K::RelativeExtension{T}, M::Generic.Mat{T}) where {T, S}
    z = NfRelOrd{T, S}(K)
    z.nf = K
    z.parent = NfRelOrdSet{T}(K)
    z.basis_mat = M
    z.basis_pmat = pseudo_matrix(M)
    return z
  end
end

###############################################################################
#
#  NfRelOrdElem
#
###############################################################################

mutable struct NfRelOrdElem{T} <: RingElem
  parent#::NfRelOrd{T, S} # I don't want to drag the S around
  elem_in_nf::RelativeElement{T}
  elem_in_basis::Vector{T}
  has_coord::Bool

  function NfRelOrdElem{T}(O::NfRelOrd{T}) where {T}
    z = new{T}()
    z.parent = O
    z.elem_in_nf = zero(nf(O))
    z.elem_in_basis = Vector{T}(degree(O))
    z.has_coord = false
    return z
  end

  function NfRelOrdElem{T}(O::NfRelOrd{T}, a::RelativeElement{T}) where {T}
    z = new{T}()
    z.parent = O
    z.elem_in_nf = a
    z.elem_in_basis = Vector{T}(degree(O))
    z.has_coord = false
    return z
  end

  function NfRelOrdElem{T}(O::NfRelOrd{T}, a::RelativeElement{T}, arr::Vector{T}) where {T}
    z = new{T}()
    z.parent = O
    z.elem_in_nf = a
    z.elem_in_basis = arr
    z.has_coord = true
    return z
  end
end

###############################################################################
#
#  NfRelOrdFracIdl
#
###############################################################################

mutable struct NfRelOrdFracIdlSet{T, S}
  order::NfRelOrd{T, S}

  function NfRelOrdFracIdlSet{T, S}(O::NfRelOrd{T, S}) where {T, S}
    a = new(O)
    return a
  end
end

mutable struct NfRelOrdFracIdl{T, S}
  order::NfRelOrd{T, S}
  parent::NfRelOrdFracIdlSet{T, S}
  basis_pmat::PMat{T, S}
  pseudo_basis::Vector{Tuple{RelativeElement{T}, S}}
  basis_mat::Generic.Mat{T}
  basis_mat_inv::Generic.Mat{T}
  den::fmpz

  norm
  has_norm::Bool

  function NfRelOrdFracIdl{T, S}(O::NfRelOrd{T, S}) where {T, S}
    z = new{T, S}()
    z.order = O
    z.parent = NfRelOrdFracIdlSet{T, S}(O)
    z.has_norm = false
    return z
  end

  function NfRelOrdFracIdl{T, S}(O::NfRelOrd{T, S}, M::PMat{T, S}) where {T, S}
    z = NfRelOrdFracIdl{T, S}(O)
    z.basis_pmat = M
    z.basis_mat = M.matrix
    return z
  end

  function NfRelOrdFracIdl{T, S}(O::NfRelOrd{T, S}, a::Array{Tuple{T1, S}}) where {T1 <: RelativeElement{T}, S} where T
    z = NfRelOrdFracIdl{T, S}(O)
    z.pseudo_basis = a
    return z
  end
end

###############################################################################
#
#  NfRelOrdIdl
#
###############################################################################

mutable struct NfRelOrdIdlSet{T, S}
  order::NfRelOrd{T, S}

  function NfRelOrdIdlSet{T, S}(O::NfRelOrd{T, S}) where {T, S}
    a = new(O)
    return a
  end
end

mutable struct NfRelOrdIdl{T, S}
  order::NfRelOrd{T, S}
  parent::NfRelOrdIdlSet{T, S}
  basis_pmat::PMat{T, S}
  pseudo_basis::Vector{Tuple{RelativeElement{T}, S}}
  basis_mat::Generic.Mat{T}
  basis_mat_inv::Generic.Mat{T}

  norm
  has_norm::Bool
  is_prime::Int            # 0: don't know
                           # 1 known to be prime
                           # 2 known to be not prime
  splitting_type::Tuple{Int, Int}

  minimum
  non_index_div_poly::fq_poly # only used if the ideal is a prime ideal not dividing the index
  p_uniformizer::NfRelOrdElem{T}
  anti_uniformizer::RelativeElement{T}

  function NfRelOrdIdl{T, S}(O::NfRelOrd{T, S}) where {T, S}
    z = new{T, S}()
    z.order = O
    z.parent = NfRelOrdIdlSet{T, S}(O)
    z.has_norm = false
    z.is_prime = 0
    z.splitting_type = (0,0)
    return z
  end

  function NfRelOrdIdl{T, S}(O::NfRelOrd{T, S}, M::PMat{T, S}) where {T, S}
    z = NfRelOrdIdl{T, S}(O)
    z.basis_pmat = M
    z.basis_mat = M.matrix
    return z
  end

  function NfRelOrdIdl{T, S}(O::NfRelOrd{T, S}, a::Array{Tuple{T1, S}}) where {T1 <: RelativeElement{T}, S} where T
    z = NfRelOrdIdl{T, S}(O)
    z.pseudo_basis = a
    return z
  end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

struct Pair{A,B}
    first::A
    second::B
end

const => = Pair

start(p::Pair) = 1
done(p::Pair, i) = i>2
next(p::Pair, i) = (getfield(p,i), i+1)
eltype(p::Pair{A,B}) where {A,B} = Union{A,B}

indexed_next(p::Pair, i::Int, state) = (getfield(p,i), i+1)

hash(p::Pair, h::UInt) = hash(p.second, hash(p.first, h))

==(p::Pair, q::Pair) = (p.first==q.first) & (p.second==q.second)
isequal(p::Pair, q::Pair) = isequal(p.first,q.first) & isequal(p.second,q.second)

isless(p::Pair, q::Pair) = ifelse(!isequal(p.first,q.first), isless(p.first,q.first),
                                                             isless(p.second,q.second))
getindex(p::Pair,i::Int) = getfield(p,i)
getindex(p::Pair,i::Real) = getfield(p, convert(Int, i))
reverse{A,B}(p::Pair{A,B}) = Pair{B,A}(p.second, p.first)

endof(p::Pair) = 2
length(p::Pair) = 2
first(p::Pair) = p.first
last(p::Pair) = p.second

convert(::Type{Pair{A,B}}, x::Pair{A,B}) where {A,B} = x
function convert(::Type{Pair{A,B}}, x::Pair) where {A,B}
    Pair{A,B}(convert(A, x[1]), convert(B, x[2]))
end

promote_rule(::Type{Pair{A1,B1}}, ::Type{Pair{A2,B2}}) where {A1,B1,A2,B2} =
    Pair{promote_type(A1, A2), promote_type(B1, B2)}

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
using Base.Test
# spec.jl
import WaveletScattering: checkspec_super, default_ɛ,
    default_max_qualityfactor, default_motherfrequency,
    default_n_filters_per_octave, default_n_octaves
# meta.jl
import WaveletScattering: get_γ, get_scale, ΦMeta, ΨMeta1D
# morlet1d.jl
import WaveletScattering: uncertainty
# spec.jl
import WaveletScattering: AbstractSpec
# spec1d.jl
import WaveletScattering: Spec1D
# waveletclass.jl
import WaveletScattering: MexicanHat, Morlet, RedundantWaveletClass

# default_motherfrequency
immutable NullWaveletClass <: RedundantWaveletClass end
class = NullWaveletClass()
for nfo in [1, 2, 4, 8, 12, 24, 32]
    ξ = default_motherfrequency(class, nfo)
    @test 2ξ ≈ (ξ*2^(-1/nfo) + (1-ξ))
end

# default_n_octaves
numerictypes = [Float16, Float32, Float64]
nfos = [1, 2, 4, 8, 12, 24, 32]
for T in numerictypes, nfo in nfos, max_q in nfos[nfos.<=nfo],
    log2_s in (7+ceil(Int, log2(nfo)):18), max_s in [max_q*exp2.(5:14); Inf]
    spec = Spec1D(
        log2_size = log2_s,
        max_qualityfactor = max_q,
        max_scale = max_s,
        n_filters_per_octave = nfo,
        signaltype = T)
    siglength = 1 << log2_s
    if max_s > siglength
        min_centerfrequency = uncertainty(spec) / siglength * max_q
    else
        min_centerfrequency = uncertainty(spec) / max_s * 1.0
    end
    n_octaves = default_n_octaves(nothing, spec.class, log2_s,
                                Float64(max_q), max_s, spec.motherfrequency,
                                nfo)
    ξs = [ meta.centerfrequency for meta in spec.ψmetas ]
    n_octaves_a = floor(Int, log2(ξs[1] / min_centerfrequency))
    n_octaves_b = log2_s - 1 - ceil(Int, log2(nfo))
    @test n_octaves == min(n_octaves_a, n_octaves_b)
end

# checkspec_super
immutable UncheckedSpec <: AbstractSpec
    ɛ::Float64
    log2_size::Int
    ϕmeta::ΦMeta
    ψmetas::Array{ΨMeta1D,3}
    class::RedundantWaveletClass
    max_qualityfactor::Float64
    max_scale::Float64
    motherfrequency::Float64
    n_filters_per_octave::Int
    n_octaves::Int
end
ɛ = 1e-3
log2_size = 13
max_qualityfactor = 12.0
max_scale = 5e3
motherfrequency = 0.45
n_filters_per_octave = 16
n_octaves = 8
ϕmeta = Spec1D().ϕmeta
ψmetas = Spec1D().ψmetas
class = Morlet()
@test_throws ErrorException checkspec_super(UncheckedSpec(-0.1, log2_size,
    ϕmeta, ψmetas, MexicanHat(),
    2.0, max_scale, motherfrequency, n_filters_per_octave,
    n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(-0.1, log2_size,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, motherfrequency, n_filters_per_octave,
    n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(1.0, log2_size,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, motherfrequency, n_filters_per_octave,
    n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, 1,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, motherfrequency, n_filters_per_octave,
    n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    0.9, max_scale, motherfrequency, n_filters_per_octave, n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, 0.0, n_filters_per_octave, n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, 0.51, n_filters_per_octave, n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, motherfrequency, 0, n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, motherfrequency, n_filters_per_octave, 0))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, motherfrequency, 11, n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    16.1, max_scale, motherfrequency, n_filters_per_octave, n_octaves))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, motherfrequency, n_filters_per_octave, 14))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    max_qualityfactor, max_scale, motherfrequency, n_filters_per_octave, 12))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, log2_size,
    ϕmeta, ψmetas, class,
    1.0, 1e1, motherfrequency, 1, 12))
@test_throws ErrorException checkspec_super(UncheckedSpec(ɛ, 8,
    ϕmeta, ψmetas, class,
    1.0, Inf, motherfrequency, 1, 1))

# default_ɛ
@test default_ɛ(Float16) ≈ 1e-3
@test default_ɛ(Float32) ≈ 1e-7
@test default_ɛ(Float64) ≈ 1e-15
@test default_ɛ(Complex{Float16}) ≈ 1e-3
@test default_ɛ(Complex{Float32}) ≈ 1e-7
@test default_ɛ(Complex{Float64}) ≈ 1e-15

# default_max_qualityfactor
@test default_max_qualityfactor(8.0, nothing) ≈ 8.0
@test default_max_qualityfactor(nothing, 8) ≈ 8.0
@test default_max_qualityfactor(nothing, nothing) ≈ 1.0

# default_n_filters_per_octave
@test default_n_filters_per_octave(12, nothing) == 12
@test default_n_filters_per_octave(nothing, 7.9) == 8
@test default_n_filters_per_octave(nothing, nothing) == 1

# default_n_octaves (fallback)
type WhateverType end
@test default_n_octaves(5, WhateverType) == 5


text:
# This file is part of Kpax3. License is MIT.

function test_settings_exceptions()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, miss=UInt8[63; 0])
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, l=-1)
  @test_throws Kpax3.KInputError  Kpax3.KSettings(ifile, ofile, gamma=[1.0; 0.0])
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, gamma=[1.0; 0.0; -1.0])
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, r=0.0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, r=-1.0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, maxclust=0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, maxunit=0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, verbosestep=-1)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, popsize=-1)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, popsize=0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, popsize=1)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, maxiter=0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, maxgap=-1)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, xrate=-1.0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, xrate=2.0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, mrate=-1.0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, mrate=2.0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, T=0)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, burnin=-1)
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, tstep=-1)
  @test_throws Kpax3.KInputError  Kpax3.KSettings(ifile, ofile, op=[1.0; 0.0])
  @test_throws Kpax3.KDomainError Kpax3.KSettings(ifile, ofile, op=[1.0; -1.0; 0.0])

  nothing
end

test_settings_exceptions()

function test_settings_constructor()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"
  protein = false
  miss = zeros(UInt8, 0)
  l = 10
  α = 0.0
  θ = 1.0
  γ = [0.3; 0.3; 0.4]
  r = log(0.01) / log(0.9)
  maxclust = 2
  maxunit = 2
  verbose = false
  verbosestep = 1
  popsize = 10
  maxiter = 5
  maxgap = 1
  xrate = 0.8
  mrate = 0.2
  T = 10
  burnin = 100
  tstep = 2
  op = [0.7; 0.3; 0.0]

  settings = Kpax3.KSettings(ifile, ofile, protein=protein, miss=miss, l=l, alpha=α, theta=θ, gamma=γ, r=r, maxclust=maxclust, maxunit=maxunit, verbose=verbose, verbosestep=verbosestep, popsize=popsize, maxiter=maxiter, maxgap=maxgap, xrate=xrate, mrate=mrate, T=T, burnin=burnin, tstep=tstep, op=op)

  @test settings.ifile == ifile
  @test settings.ofile == ofile
  @test settings.l == l
  @test settings.α == α
  @test settings.θ == θ
  @test settings.γ == γ
  @test settings.r == r
  @test settings.maxclust == maxclust
  @test settings.maxunit == maxunit
  @test settings.verbose == verbose
  @test settings.verbosestep == verbosestep
  @test settings.popsize == popsize
  @test settings.maxiter == maxiter
  @test settings.maxgap == maxgap
  @test settings.xrate == xrate
  @test settings.mrate == mrate
  @test settings.T == T
  @test settings.burnin == burnin
  @test settings.tstep == tstep
  @test isa(settings.op, StatsBase.ProbabilityWeights)
  @test values(settings.op) == op

  settings = Kpax3.KSettings(ifile, ofile, protein=true, miss=zeros(UInt8, 0))

  @test settings.protein
  @test settings.miss == UInt8['?', '*', '#', '-', 'b', 'j', 'x', 'z']

  settings = Kpax3.KSettings(ifile, ofile, protein=true, miss=UInt8['?', '*', '#', 'b', 'j', 'x', 'z'])

  @test settings.protein
  @test settings.miss == UInt8['?', '*', '#', 'b', 'j', 'x', 'z']

  settings = Kpax3.KSettings(ifile, ofile, protein=false, miss=zeros(UInt8, 0))

  @test !settings.protein
  @test settings.miss == UInt8['?', '*', '#', '-', 'b', 'd', 'h', 'k', 'm', 'n', 'r', 's', 'v', 'w', 'x', 'y', 'j', 'z']

  settings = Kpax3.KSettings(ifile, ofile, protein=false, miss=UInt8['?', '*', '#', 'b', 'd', 'h', 'k', 'm', 'n', 'r', 's', 'v', 'w', 'x', 'y', 'j', 'z'])

  @test !settings.protein
  @test settings.miss == UInt8['?', '*', '#', 'b', 'd', 'h', 'k', 'm', 'n', 'r', 's', 'v', 'w', 'x', 'y', 'j', 'z']

  settings = Kpax3.KSettings(ifile, ofile, misscsv=Array{String}(undef, 0))

  @test settings.misscsv == [""]

  settings = Kpax3.KSettings(ifile, ofile, misscsv=["", "?", "X", "NA", "-", "."])

  @test settings.misscsv == ["", "?", "X", "NA", "-", "."]

  settings = Kpax3.KSettings(ifile, ofile, misscsv=["?", "X", "NA", "-", "."])

  @test settings.misscsv == ["", "?", "X", "NA", "-", "."]

  nothing
end

test_settings_constructor()


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LibOSXUnwind_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LibOSXUnwind")
JLLWrappers.@generate_main_file("LibOSXUnwind", UUID("a83860b7-747b-57cf-bf1f-3e79990d037f"))
end  # module LibOSXUnwind_jll


text:
module Laplacians2DSymmetric

using ..Lattices
using ..Points
using SparseArrays

import ..Symmetric, ..Unsymmetric
import LinearAlgebra: I
import ..Laplacian

################################################################################
# CARTESIAN
"""
	Laplacian(::Lattice{2}, α::Vector, α_half::Vector) -> lap
"""
function Laplacian{Symmetric}(
			lattice::Lattice{2,Cartesian},
			αx::Vector{ComplexF64},
			αy::Vector{ComplexF64},
			α_halfx::Vector{ComplexF64},
			α_halfy::Vector{ComplexF64})

	dx = lattice.dx
	dy = lattice.dy

	Nx = length(αx)
	Ny = length(αy)
	N = Nx*Ny

	α_halfx⁻¹ = 1 ./α_halfx
	α_halfy⁻¹ = 1 ./α_halfy

	# note that x_half[1] = x[1] - dx/2, so that D₋½=D(x_half[1])
	# and that D½=D(x_half[2])
	rows0 = 1:Nx	; cols0 = 1:Nx	; vals0 = - α_halfx⁻¹[rows0] - α_halfx⁻¹[rows0.+1]
	rowsm = 2:Nx	; colsm = 1:Nx-1; valsm = α_halfx⁻¹[1 .+ colsm]
	rowsp = 1:Nx-1	; colsp = 2:Nx	; valsp = α_halfx⁻¹[colsp]
	∂ₓα⁻¹∂ₓ = sparse(vcat(rowsm,rows0,rowsp),vcat(colsm,cols0,colsp),vcat(valsm,vals0,valsp)/lattice.dx^2,Nx,Nx)

	rows0 = 1:Ny	; cols0 = 1:Ny	; vals0 = - α_halfy⁻¹[rows0] - α_halfy⁻¹[rows0.+1]
	rowsm = 2:Ny	; colsm = 1:Ny-1; valsm = α_halfy⁻¹[1 .+ colsm]
	rowsp = 1:Ny-1	; colsp = 2:Ny	; valsp = α_halfy⁻¹[colsp]
	∂ᵤα⁻¹∂ᵤ = sparse(vcat(rowsm,rows0,rowsp),vcat(colsm,cols0,colsp),vcat(valsm,vals0,valsp)/lattice.dy^2,Ny,Ny)

	Ax = spdiagm(0=>αx)
	Ay = spdiagm(0=>αy)

	∇² = kron(Ay, ∂ₓα⁻¹∂ₓ) + kron(∂ᵤα⁻¹∂ᵤ, Ax)

	return Laplacian{2,Symmetric}(∇², kron(Ay,Ax))
end

################################################################################
# POLAR

function Laplacian{Symmetric}(lattice::Lattice{2,Polar}, α::Vector{ComplexF64}, α_half::Vector{ComplexF64})#, nnm, nnp, indices, interior, surface)
	throw("not implemented for 2D Polar yet")
	N = length(α)

	dr = lattice.dr
	r = (1/2 .+ (0:N-1))*dr

	α_half⁻¹ = 1 ./α_half
	r_half = (0:N)*dr

	∂ₓα⁻¹∂ₓ_bulk = _bulk_laplacian(lattice,N,interior,surface,nnm,nnp,α_half⁻¹)
	∂ₓα⁻¹∂ₓ_surface = _surface_laplacian(lattice,N,interior,surface,nnm,nnp,α_half⁻¹,indices)

	∂ₓ_bulk = _bulk_derivative(lattice,N,interior,surface,nnm,nnp)
	∂ₓ_surface = _surface_derivative(lattice,N,interior,surface,nnm,nnp,indices)
	∂ₓ = ∂ₓ_bulk + ∂ₓ_surface
	i∂ₓ = 1im*∂ₓ

	return Laplacian{2}(∂ₓα⁻¹∂ₓ)
end

function _bulk_laplacian(lattice::Lattice{2},N,interior,surface,nnm,nnp,α_half⁻¹)
	rows = findall(interior .& .!surface)
	colsm = nnm[1][interior .& .!surface]
	colsp = nnp[1][interior .& .!surface]
	valsm = α_half⁻¹[rows] # note that x_half[1] = x[1] - dx/2, so that D₋½=D(x_half[1])
	valsp = α_half⁻¹[colsp] #  and D½=D(x_half[2])
	cols0 = rows
	vals0 = - valsm - valsp
	return sparse(vcat(rows,rows,rows),vcat(colsm,cols0,colsp),vcat(valsm,vals0,valsp)/lattice.dx^2,N,N)
end
function _surface_laplacian(lattice::Lattice{2},N,interior,surface,nnm,nnp,α_half⁻¹,indices)
	rows_sf = findall(surface .& interior)
	cols_sfm = map(i->findfirst(isequal(indices[i]-CartesianIndex(1,)),indices),rows_sf)
	rows_sfm = [rows_sf[findfirst(.!isnothing.(cols_sfm))]]
	cols_sfm = [cols_sfm[findfirst(.!isnothing.(cols_sfm))]]
	vals_sfm = α_half⁻¹[rows_sfm] # note that x_half[1] = x[1] - dx/2, so that D₋½=D(x_half[1])
	cols_sfp = map(i->findfirst(isequal(indices[i]+CartesianIndex(1,)),indices),rows_sf)
	rows_sfp = [rows_sf[findfirst(.!isnothing.(cols_sfp))]]
	cols_sfp = [cols_sfp[findfirst(.!isnothing.(cols_sfp))]]
	vals_sfp = α_half⁻¹[cols_sfp] # and D½=D(x_half[2])
	cols_sf0 = rows_sf
	vals_sf0 = - α_half⁻¹[[0,1] .+ rows_sf] - α_half⁻¹[[0,1] .+ rows_sf]
	return sparse(vcat(rows_sfm,rows_sf,rows_sfp),vcat(cols_sfm,cols_sf0,cols_sfp),vcat(vals_sfm,vals_sf0,vals_sfp)/lattice.dx^2,N,N)
end

function _bulk_derivative(lattice::Lattice{2},N,interior,surface,nnm,nnp)
	rows = findall(interior .& .!surface)
	colsm = nnm[1][interior .& .!surface]
	valsm = -ones(size(colsm))
	colsp = nnp[1][interior .& .!surface]
	valsp = ones(size(colsm))
	return sparse(vcat(rows,rows),vcat(colsm,colsp),vcat(valsm,valsp)/lattice.dx,N,N)
end
function _surface_derivative(lattice::Lattice{2},N,interior,surface,nnm,nnp,indices)
	rows_sf = findall(surface .& interior)
	cols_sfm = map(i->findfirst(isequal(indices[i]-CartesianIndex(1,)),indices),rows_sf)
	rows_sfm = [rows_sf[findfirst(.!isnothing.(cols_sfm))]]
	cols_sfm = [cols_sfm[findfirst(.!isnothing.(cols_sfm))]]
	vals_sfm = -ones(size(rows_sfm))
	cols_sfp = map(i->findfirst(isequal(indices[i]+CartesianIndex(1,)),indices),rows_sf)
	rows_sfp = [rows_sf[findfirst(.!isnothing.(cols_sfp))]]
	cols_sfp = [cols_sfp[findfirst(.!isnothing.(cols_sfp))]]
	vals_sfp = ones(size(rows_sfm))
	return sparse(vcat(rows_sfm,rows_sfp),vcat(cols_sfm,cols_sfp),vcat(vals_sfm,vals_sfp)/lattice.dx,N,N)
end

################################################################################
# SPHERICAL
Laplacian{Symmetric}(lattice::Lattice{2,Spherical}, α::Vector{ComplexF64}, α_half::Vector{ComplexF64}, nnm, nnp, indices, interior, surface) = throw("not implemented for 2D Spherical yet")

end

using .Laplacians2DSymmetric


text:
# gpu versn = 20000
x = randn(10, n)
using GPUArrays
using CUBLAS

function gauss{T}(X::Matrix{T}, alpha::T)
    n = size(X, 2)
    xy = LinAlg.BLAS.syrk('U', 'T', one(T), X)
    x2 = [ sum(X[:, i] .^ 2) for i in 1:n ]

    LinAlg.BLAS.syr2k!('U', 'N', one(T), x2, ones(T, n), T(-2), xy)

    @inbounds for i in 1:n
        for j in 1:i
            xy[j, i] = exp(-alpha * xy[j, i])
        end
    end
    LinAlg.copytri!(xy, 'U')
end
gauss(x, 1.0)ions

type KRR_GPU{T <: AbstractFloat} <: AbstractKRR{T}
    λ :: T
    X :: Matrix{T}
    α :: Vector{T}
    ϕ :: MLKernels.MercerKernel{T}

    function KRR(λ, X, α, ϕ)
        @assert λ > zero(λ)
        @assert size(X, 2) == length(α)
        new(λ, X, α, ϕ)
    end
end

function KRR{T <: AbstractFloat}(
    λ :: T,
    X :: Matrix{T},
    α :: Vector{T},
    ϕ :: MLKernels.MercerKernel{T}
)
    KRR{T}(λ, X, α, ϕ)
end

function StatsBase.fit{T <: AbstractFloat}(
      :: Type{KRR},
    X :: Matrix{T},
    y :: Vector{T},
    λ :: T,
    ϕ :: MLKernels.Kernel{T}
)
    d, n = size(X)
    K = MLKernels.kernelmatrix!(MLKernels.ColumnMajor(),
                                Matrix{T}(n, n),
                                ϕ, X, true)
    for i = 1:n
        # the n is important to make things comparable between fast and normal
        # KRR
        @inbounds K[i, i] += n * λ
    end

    α = cholfact!(K) \ y

    KRR(λ, X, α, ϕ)
end

function StatsBase.predict!{T <: AbstractFloat}(
    KRR :: KRR{T},
    X   :: Matrix{T},
    y   :: Vector{T},
    K   :: Matrix{T}
)
    n, n_new = (size(KRR.X, 2), size(X, 2))
    @assert (n_new, n) == size(K)
    @assert length(y) == n_new

    MLKernels.kernelmatrix!(
        MLKernels.ColumnMajor(),
        K, KRR.ϕ, X, KRR.X
    )
    LinAlg.BLAS.gemv!('N', one(T), K, KRR.α, zero(T), y)
    return y
end

function StatsBase.predict!{T <: AbstractFloat}(KRR::KRR{T}, X::Matrix{T}, y::Vector{T})
    StatsBase.predict!(
        KRR, X, y, MLKernels.kernelmatrix!(
            MLKernels.ColumnMajor(),
            Matrix{T}(size(X, 2), size(KRR.X, 2)),
            KRR.ϕ, X, KRR.X
    ))
end

function StatsBase.predict{T <: AbstractFloat}(KRR::KRR{T}, X::Matrix{T})
    StatsBase.predict!(KRR, X, Vector{T}(size(X, 2)))
end

# predict the KRR with the data in X and add the result to y, used to speedup
# predict(fast_krr::FastKRR)
function predict_and_add!{T <: AbstractFloat}(
    KRR :: KRR{T},
    X   :: Matrix{T},
    y   :: Vector{T},
    K   :: Matrix{T}
)
    n, n_new = (size(KRR.X, 2), size(X, 2))
    @assert (n_new, n) == size(K)
    @assert length(y) == n_new

    MLKernels.kernelmatrix!(
        MLKernels.ColumnMajor(),
        K, KRR.ϕ, X, KRR.X
    )
    LinAlg.BLAS.gemv!('N', one(T), K, KRR.α, one(T), y)
    return y
end


text:
"""
    x0_volume_liquid(model,T,z)

Returns an initial guess to the liquid volume, dependent on temperature and composition. by default is 1.25 times the [lower bound volume](@ref lb_volume).
"""
function x0_volume_liquid(model,T,z)
    v_lb = lb_volume(model,z)
    return v_lb*1.25
end

"""
    x0_volume_gas(model,p,T,z)

Returns an initial guess to the gas volume, depending of pressure, temperature and composition. by default uses a [virial aproximation](@ref volume_virial)
"""
function x0_volume_gas(model,p,T,z)
    return volume_virial(model,p,T,z)
end

"""
    x0_volume(model,p,T,z; phase = :unknown)

Returns an initial guess of the volume at a pressure, temperature, composition and suggested phase.

If the suggested phase is `:unkwown` or `:liquid`, calls [`x0_volume_liquid`](@ref).

If the suggested phase is `:gas`, calls [`x0_volume_gas`](@ref).

"""
function x0_volume(model,p,T,z; phase = :unknown)
    phase = Symbol(phase)
    if phase === :unknown || is_liquid(phase)
        return x0_volume_liquid(model,T,z)
    elseif is_vapour(phase)
        return x0_volume_gas(model,p,T,z)
    elseif is_supercritical(phase)
     else
        error("unreachable state on x0_volume")
    end
end

"""
    lb_volume(model::EoSModel,z=SA[1.0])

Returns the lower bound volume.

It has different meanings depending on the Equation of State, but symbolizes the minimum allowable volume at a certain composition:

- SAFT EoS: the packing volume
- Cubic EoS, covolume (b) parameter

On empiric equations of state, the value is chosen to match the volume of the conditions at maximum pressure and minimum temperature
, but the equation itself normally can be evaluated at lower volumes.

On SAFT and Cubic EoS, volumes lower than `lb_volume` will likely error.

The lower bound volume is used for guesses of liquid volumes at a certain pressure, saturated liquid volumes and critical volumes.
"""
function lb_volume end

"""
    x0_sat_pure(model::EoSModel,T,z=SA[1.0])

Returns a 2-tuple corresponding to `(log10(Vₗ),log10(Vᵥ))`, where Vₗ and Vᵥ are the liquid and vapor initial guesses.
Used in [`sat_pure`](@ref).
"""
function x0_sat_pure(model,T,z=SA[1.0])
    #=theory as follows
    #given T = Teos:
    #calculate B(Teos)
    PB = Peos = -2B

    veos = volume_compress(model,PB,T)
    with a (P,V,T,B) pair, change to
    (P,V,a,b)
    =#
    B = second_virial_coefficient(model,T,SA[1.0])
    lb_v = lb_volume(model,SA[1.0])*one(T)
    _0 = zero(B)
    #virial volume below lower bound volume.
    #that means that we are way over the critical point
    if -2B < lb_v
        _nan = _0/_0
        return (_nan,_nan)
    end
    p = -0.25*R̄*T/B
    vl = x0_volume(model,p,T,z,phase=:l)
    vl = _volume_compress(model,p,T,SA[1.0],vl)
    #=the basis is that p = RT/v-b - a/v2
    we have a (p,v,T) pair
    and B = 2nd virial coefficient = b-a/RT
    with that, we solve for a and b
    as a and b are vdW, Pc and Tc can be calculated
    with Tc and Pc, we use [1] to calculate vl0 and vv0
    with Tc, we can also know in what regime we are.
    in near critical pressures, we use directly vv0 = -2B
    and vl0 = 4*lb_v

    [1]
    DOI: 10.1007/s10910-007-9272-4
    Journal of Mathematical Chemistry, Vol. 43, No. 4, May 2008 (© 2007)
    The van der Waals equation: analytical and approximate solutions
    =#
    γ = p*vl*vl/(R̄*T)
    _c = vl*vl + B*vl - γ*vl
    _b = γ - B - vl
    Δ = _b*_b - 4*_c
    if isnan(vl) | (Δ < 0)
        #fails on two ocassions:
        #near critical point, or too low.
        #old strategy
        x0l = 4*lb_v
        x0v = -2*B + 2*lb_v
        return (log10(x0l),log10(x0v))
    end
    Δsqrt = sqrt(Δ)
    b1 = 0.5*(-_b + Δsqrt)
    b2 = 0.5*(-_b - Δsqrt)
    if b1 < 0
        b = b2
    elseif b1 > vl
        b = b2
    else
        b = b1
    end
    a = -R̄*T*(B-b)
    Vc = 3*b
    Ωa =  27/64
    Ωb =  1/8
    ar = a/Ωa
    br = b/Ωb
    Tc = ar/br/R̄
    Pc = ar/(br*br)
    Tr = T/Tc
    #Tr(vdW approx) > Tr(model)
    if Tr >= 1
        _nan = _0/_0
        return (_nan,_nan)
    end
    # if b1/b2 < -0.95, then T is near Tc.
    #if b<lb_v then we are in trouble
    #critical regime or something horribly wrong happened
    if (b1/b2 < -0.95) | (b<lb_v) | (Tr>0.99)
        x0l = 4*lb_v
        x0v = -2*B #gas volume as high as possible
        return (log10(x0l),log10(x0v))
    end
    Vl0,Vv0 = vdw_x0_xat_pure(T,Tc,Pc,Vc)
    x0l = min(Vl0,vl)
    x0v = min(1e4*one(Vv0),Vv0) #cutoff volume
    return (log10(x0l),log10(x0v))
end

function vdw_x0_xat_pure(T,T_c,P_c,V_c)
    Tr = T/T_c
    Trm1 = 1.0-Tr
    Trmid = sqrt(Trm1)
    if Tr >= 0.7
        c_l = 1.0+2.0*Trmid + 0.4*Trm1 - 0.52*Trmid*Trm1 +0.115*Trm1*Trm1 #Eq. 29
    else
        c_l = 1.5*(1+sqrt(1-(32/27)*Tr)) #Eq. 32
    end

    if Tr >= 0.46
        #Eq. 30, valid in 0.46 < Tr < 1
        c_v = 1.0-2.0*Trmid + 0.4*Trm1 + 0.52*Trmid*Trm1 +0.207*Trm1*Trm1
    elseif Tr <= 0.33
        #Eq. 33, valid in 0 < Tr < 0.33
        c_v = (3*c_l/(ℯ*(3-c_l)))*exp(-(1/(1-c_l/3)))

    else
        #Eq. 31 valid in 0.25 < Tr < 1
        mean_c = 1.0 + 0.4*Trm1 + 0.161*Trm1*Trm1
        c_v = 2*mean_c - c_l
    end
    #volumes predicted by vdW
    Vl0 = (1/c_l)*V_c
    Vv0 = (1/c_v)*V_c
    return (Vl0,Vv0)
end

function scale_sat_pure(model,z=SA[1.0])
    p    = 1/p_scale(model,z)
    μ    = 1/R̄/T_scale(model,z)
    return p,μ
end

"""
    x0_crit_pure(model::SAFTModel)

Returns a 2-tuple corresponding to
    `(k,log10(Vc0))`, where `k` is `Tc0/T_scale(model,z)`
"""
function x0_crit_pure end

function x0_crit_pure(model::EoSModel)
    lb_v = lb_volume(model)
    (1.5, log10(lb_v/0.3))
end

"""
    T_scale(model::EoS,z=SA[1.0])

Represents a temperature scaling factor.

On any EoS based on Critical parameters (Cubic or Empiric EoS), the temperature scaling factor is chosen to be the critical temperature.

On SAFT or other molecular EoS, the temperature scaling factor is chosen to be a function of the potential depth ϵ.

Used as scaling factors in [`saturation_pressure`](@ref) and as input for solving [`crit_pure`](@ref)
"""
function T_scale end

"""
    p_scale(model::SAFTModel,z=SA[1.0])

Represents a pressure scaling factor

On any EoS based on Critical parameters (Cubic or
Empiric EoS), the pressure scaling factor is
chosen to be a function of the critical pressure.

On SAFT or other molecular EoS, the temperature
scaling factor is chosen to a function of ∑(zᵢ*ϵᵢ*(σᵢᵢ)³)

Used as scaling factors in [`saturation_pressure`](@ref) and as input for solving [`crit_pure`](@ref)

"""
function p_scale end

#=
the following methods are fallbacks,
that require just the definition of T_scale,p_scale and lb_volume
respectively. if possible, each eos should define those
=#
function T_scales(model,z)
    n = length(z)
    x = zeros(n)
    res = zeros(n)
    for i = 1:n
        x[i] = 1.0
        res[i] = T_scale(model,x)
        x[i] = 0.0
    end
    return res
end


text:
d = @destruct [a, b] = Dict(:a => 1, :b => 2)
@test d == Dict(:a => 1, :b => 2)
@test (a, b) == (1, 2)

@destruct [a] = Dict("a" => "foo")
@test a == "foo"

@destruct [foo = :a || 5, b = :b || 6, c || 7] = Dict(:a => 1)
@test (foo, b, c) == (1, 6, 7)

@destruct _.data = "foo"
@test _ == "foo"
@test data == "foo".data

@destruct _.(re, im) = Complex(1,2)
@test (re, im) == (1, 2)

@destruct [s.data = :a] = Dict(:a => "foo")
@test s == "foo"
@test data === s.data

@destruct x[a, [c, d] = b] = Dict(:a => 1, :b => Dict(:c => 2, :d => 3))
@test x == Dict(:a => 1, :b => Dict(:c => 2, :d => 3))
@test (a, c, d) == (1, 2, 3)
@test b == Dict(:c => 2, :d => 3)


text:
using DiffEqFlux, Test

@testset "DiffEqFlux" begin

include("layers.jl")
include("utils.jl")
include("neural_de.jl")

end


text:
struct LatLonSeparableDriver <: DriverModel{LatLonAccel}
    mlat::LateralDriverModel
    mlon::LaneFollowingDriver
end

function observe!(model::LaneFollowingDriver, scene::Scene, roadway::Roadway, egoid::Int)

    vehicle_index = findfirst(egoid, scene)

    fore = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())

    v_ego = scene[vehicle_index].state.v
    v_oth = NaN
    headway = NaN

    if fore.ind != nothing
        v_oth = scene[fore.ind].state.v
        headway = fore.Δs
    end

    track_longitudinal!(model, v_ego, v_oth, headway)

    return model
end

get_name(model::LatLonSeparableDriver) = @sprintf("%s + %s", get_name(model.mlat), get_name(model.mlon))
function set_desired_speed!(model::LatLonSeparableDriver, v_des::Float64)
    set_desired_speed!(model.mlon, v_des)
    model
end
function observe!(model::LatLonSeparableDriver, scene::Scene, roadway::Roadway, egoid::Int)
    observe!(model.mlat, scene, roadway, egoid)
    observe!(model.mlon, scene, roadway, egoid)
    model
end
function Base.rand(model::LatLonSeparableDriver)
    alat = rand(model.mlat)
    alon = rand(model.mlon).a
    LatLonAccel(alat, alon)
end
Distributions.pdf(model::LatLonSeparableDriver, a::LatLonAccel) = pdf(model.mlat, a.a_lat) * pdf(model.mlon, a.a_lon)
Distributions.logpdf(model::LatLonSeparableDriver, a::LatLonAccel) = logpdf(model.mlat, a.a_lat) + logpdf(model.mlon, a.a_lon)

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
module CbcMathProgSolverInterface

using Cbc.CbcCInterface

using SparseArrays

import MathProgBase
const MPB = MathProgBase

export CbcMathProgModel, CbcSolver

mutable struct CbcMathProgModel <: MPB.AbstractLinearQuadraticModel
    inner::CbcModel
    check_warmstart::Bool
    binaries::Vector{Int} # indices of binary variables
end

mutable struct CbcSolver <: MPB.AbstractMathProgSolver
    options
end
CbcSolver(;kwargs...) = CbcSolver(kwargs)

function CbcMathProgModel(;check_warmstart::Bool=true,options...)
    c = CbcModel()
    setParameter(c, "logLevel", "0")
    for (optname, optval) in options
        setParameter(c, string(optname), string(optval))
    end
    return CbcMathProgModel(c, check_warmstart, Int[])
end

MPB.LinearQuadraticModel(s::CbcSolver) = CbcMathProgModel(;s.options...)

MPB.ConicModel(s::CbcSolver) = MPB.LPQPtoConicBridge(MPB.LinearQuadraticModel(s))
MPB.supportedcones(s::CbcSolver) = [:Free,:Zero,:NonNeg,:NonPos]

function MPB.setparameters!(s::CbcSolver; mpboptions...)
    opts = collect(Any, s.options)
    silent = false
    for (optname, optval) in mpboptions
        if optname == :TimeLimit
            @static if VERSION >= v"0.7-"
                push!(opts, :seconds => optval)
            else
                push!(opts, (:seconds, optval))
            end
        elseif optname == :Silent
            if optval == true
                @static if VERSION >= v"0.7-"
                    push!(opts, :logLevel => 0)
                else
                    push!(opts, (:logLevel, 0))
                end
            end
        else
            error("Unrecognized parameter $optname")
        end
    end
    s.options = opts
end

function MPB.setparameters!(s::CbcMathProgModel; mpboptions...)
    for (optname, optval) in mpboptions
        if optname == :TimeLimit
            setParameter(s.inner, "seconds", string(optval))
        elseif optname == :Silent
            if optval == true
                setParameter(s.inner, "logLevel","0")
            end
        else
            error("Unrecognized parameter $optname")
        end
    end
end

function MPB.loadproblem!(m::CbcMathProgModel, A, collb, colub, obj, rowlb, rowub, sense)
    loadProblem(m.inner, A, collb, colub, obj, rowlb, rowub)
    MPB.setsense!(m, sense)
end

function MPB.writeproblem(m::CbcMathProgModel, filename::AbstractString)
    if endswith(filename,".mps")
        writeMps(m.inner, filename)
    else
        error("Only MPS output supported")
    end
end

function MPB.setsense!(m::CbcMathProgModel,sense)
    @assert sense == :Min || sense == :Max
    if sense == :Min
       setObjSense(m.inner, 1)
    else
       setObjSense(m.inner, -1)
    end
end

function MPB.getsense(m::CbcMathProgModel)
    s = getObjSense(m.inner)
    if s == 1
        return :Min
    elseif s == -1
        return :Max
    else
        error("Internal error: Unknown sense $s")
    end
end

MPB.numvar(m::CbcMathProgModel) = getNumCols(m.inner)
MPB.numconstr(m::CbcMathProgModel) = getNumRows(m.inner)

function MPB.getvartype(m::CbcMathProgModel)
    ncol = MPB.numvar(m)
    vartype = fill(:Cont,ncol)
    for i in 1:ncol
        if isInteger(m.inner, i-1)
            vartype[i] = :Int
        end
    end
    for k in m.binaries
        vartype[k] = :Bin
    end
    return vartype
end

function MPB.setvartype!(m::CbcMathProgModel,vartype::Vector{Symbol})
    ncol = MPB.numvar(m)
    @assert length(vartype) == ncol
    m.binaries = Int[]
    for i in 1:ncol
        if vartype[i] == :Int
            setInteger(m.inner, i-1)
        elseif vartype[i] == :Cont
            setContinuous(m.inner, i-1)
        elseif vartype[i] == :Bin
            setInteger(m.inner, i-1)
            push!(m.binaries, i)
        else
            error("Unrecognized variable type $(vartype[i])")
        end
    end
end

function MPB.optimize!(m::CbcMathProgModel)
    lb = getColLower(m.inner)
    ub = getColUpper(m.inner)
    # tighten bounds on binaries
    for i in m.binaries
        setColLower(m.inner, i-1, max(lb[i],0.0))
        setColUpper(m.inner, i-1, min(ub[i],1.0))
    end
    solve(m.inner)
end

function MPB.status(m::CbcMathProgModel)
    if isProvenOptimal(m.inner)
        return :Optimal
    elseif isProvenInfeasible(m.inner)
        return :Infeasible
    elseif isContinuousUnbounded(m.inner)
        return :Unbounded # is this correct?
    elseif isNodeLimitReached(m.inner) || isSecondsLimitReached(m.inner) || isSolutionLimitReached(m.inner)
        return :UserLimit
    elseif isAbandoned(m.inner)
        return :Error
    else
        error("Internal error: Unrecognized solution status")
    end
end

function MPB.getconstrmatrix(m::CbcMathProgModel)
    starts = getVectorStarts(m.inner)
    rowval = getIndices(m.inner)
    nzval = getElements(m.inner)
    return SparseMatrixCSC(MPB.numconstr(m), MPB.numvar(m), starts .+ 1,
                           rowval .+ 1, nzval)
end

MPB.getobjval(m::CbcMathProgModel) = getObjValue(m.inner)

MPB.getobjbound(m::CbcMathProgModel) = getBestPossibleObjValue(m.inner)

function MPB.getobjgap(m::CbcMathProgModel)
    b = getBestPossibleObjValue(m.inner)
    f = getObjValue(m.inner)
    return abs(b-f)/abs(f)
end

MPB.getnodecount(m::CbcMathProgModel) = getNodeCount(m.inner)

MPB.getsolution(m::CbcMathProgModel) = getColSolution(m.inner)

MPB.getrawsolver(m::CbcMathProgModel) = m.inner

function MPB.setwarmstart!(m::CbcMathProgModel, v)
    if any(isnan, v)
        @warn("Ignoring partial starting solution. Cbc requires a feasible " *
              "value to be specified for all variables.", maxlog = 1)
        return
    end

    # ignore if not feasible
    if m.check_warmstart
        @assert length(v) == MPB.numvar(m)
        l = getColLower(m.inner)
        u = getColUpper(m.inner)
        for i in 1:length(v)
            if !(l[i] - 1e-6 <= v[i] <= u[i] + 1e-6)
                return
            end
            if isInteger(m.inner, i-1) && !isinteger(l[i])
                return
            end
        end
        lb = getRowLower(m.inner)
        ub = getRowUpper(m.inner)
        A = MPB.getconstrmatrix(m)
        rowval = A*v
        for i in 1:MPB.numconstr(m)
            if !(lb[i] - 1e-6 <= rowval[i] <= ub[i] + 1e-6)
                return
            end
        end
    end
    setInitialSolution(m.inner, v)
end

function MPB.addsos1!(m::CbcMathProgModel, idx, weight::Vector{Float64})
    idxc = convert(Vector{Cint}, idx)
    addSOS(m.inner, 1, Cint[1,length(idx)+1], idxc, weight, 1)
end

function MPB.addsos2!(m::CbcMathProgModel, idx, weight::Vector{Float64})
    idxc = convert(Vector{Cint}, idx)
    addSOS(m.inner, 1, Cint[1,length(idx)+1], idxc, weight, 2)
end

end


text:
module CompositeTest

using Test
using ForneyLab
using ForneyLab: @composite, outboundType, isApplicable, messagePassingSchedule
using ForneyLab: SPClamp, SPGaussianMeanVarianceOutNPP, Product, condense, flatten, step!, setTargets!

# Define new node type called StateTransition, with exposed variables called (y, x_prev, x):
@composite StateTransition (y, x_prev, x) begin
    @RV x ~ GaussianMeanVariance(x_prev, constant(1.0))
    @RV y ~ GaussianMeanVariance(x, constant(1.0))
end

@testset "@composite" begin
    @test StateTransition <: FactorNode
end

@testset "Composite node construction" begin
    g = FactorGraph()
    nd = StateTransition(Variable(), Variable(), Variable())

    # Node fields should be of correct types
    @test isa(nd.id, Symbol)
    @test isa(nd.interfaces, Vector{Interface})
    @test isa(nd.i, Dict)
    @test isa(nd.inner_graph, FactorGraph)
    @test isa(nd.terminals, Vector{Terminal})

    # Node constructor should automatically assign an id
    @test nd.id == :statetransition_1

    # Node constructor should assign interfaces to itself
    for iface in nd.interfaces
        @test ===(iface.node, nd)
    end

    # Node constructor should add node to graph
    @test g.nodes[:statetransition_1] == nd
end

#-------------
# Update rules
#-------------

# Define custom rule for sum-product message towards x
@sumProductRule(:node_type     => StateTransition,
                :outbound_type => Message{Gaussian},
                :inbound_types => (Message{PointMass}, Message{Gaussian}, Nothing),
                :name          => SPStateTransitionX)

@testset "Custom SPStateTransitionX" begin
    @test SPStateTransitionX <: SumProductRule{StateTransition}
    @test outboundType(SPStateTransitionX) == Message{Gaussian}
    @test isApplicable(SPStateTransitionX, [Message{PointMass}, Message{Gaussian}, Nothing])
    @test !isApplicable(SPStateTransitionX, [Message{Gaussian}, Message{PointMass}, Nothing])
end

@testset "Composite node scheduling and algorithm compilation" begin
    fg = FactorGraph()

    x_prev = Variable(id=:x_prev)
    nd = GaussianMeanVariance(x_prev, constant(0.0), constant(1.0))
    x = Variable(id=:x)
    y = Variable(id=:y)
    cnd = StateTransition(placeholder(y, :y), x_prev, x)
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    algo = InferenceAlgorithm(pfz)

    # Build SP schedule
    schedule = messagePassingSchedule(pf)
    @test length(schedule) == 5
    @test ScheduleEntry(nd.i[:m].partner, SPClamp{Univariate}) in schedule
    @test ScheduleEntry(nd.i[:v].partner, SPClamp{Univariate}) in schedule
    @test ScheduleEntry(nd.i[:out], SPGaussianMeanVarianceOutNPP) in schedule
    @test ScheduleEntry(cnd.i[:y].partner, SPClamp{Univariate}) in schedule
    @test ScheduleEntry(cnd.i[:x], SPStateTransitionX) in schedule
    pf.schedule = condense(flatten(schedule))

    # Build marginal schedule
    marginal_table = marginalTable(x)
    @test length(marginal_table) == 1
    @test marginal_table[1].target == x
    @test marginal_table[1].interfaces[1] == cnd.i[:x]
    @test marginal_table[1].marginal_update_rule == Nothing
    pf.marginal_table = marginal_table

    # Build SP algorithm for Julia execution
    ForneyLab.assembleInferenceAlgorithm!(algo)
    code = ForneyLab.algorithmSourceCode(algo)

    @test occursin("Array{Message}(undef, 2)", code)
    @test occursin("messages[1] = ruleSPGaussianMeanVarianceOutNPP(nothing, Message(Univariate, PointMass, m=0.0), Message(Univariate, PointMass, m=1.0))", code)
    @test occursin("messages[2] = ruleSPStateTransitionX(Message(Univariate, PointMass, m=data[:y]), messages[1], nothing)", code)
    @test occursin("marginals[:x] = messages[2].dist", code)
end

@testset "Composite node algorithm execution" begin
    # Implement custom rule for Julia execution
    ruleSPStateTransitionX(::Message{PointMass, Univariate}, ::Message{F, Univariate}, ::Nothing) where F<:Gaussian = Message(Univariate, GaussianMeanVariance, m=2.0, v=3.0) # Send some dummy message

    # Resulting algorithm ---
    function step!(marginals::Dict, data::Dict)

    messages = Array{Message}(undef, 2)

    messages[1] = ruleSPGaussianMeanVarianceOutNPP(nothing, Message(Univariate, PointMass, m=0.0), Message(Univariate, PointMass, m=1.0))
    messages[2] = ruleSPStateTransitionX(Message(Univariate, PointMass, m=data[:y]), messages[1], nothing)

    marginals[:x] = messages[2].dist

    end
    # ---

    marginals = Dict()
    data = Dict(:y => 1.0)
    step!(marginals, data)

    @test marginals[:x] == ProbabilityDistribution(Univariate, GaussianMeanVariance, m=2.0, v=3.0)
end

# Composite node definition without custom rule
@composite StateTransition2 (y, x_prev, x) begin
    @RV x ~ GaussianMeanVariance(x_prev, 1.0)
    @RV y ~ GaussianMeanVariance(x, 1.0)
end

@testset "Composite node algorithm compilation without custom rule" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(1.0, 0.0)
    @RV y ~ GaussianMeanVariance(1.0, 0.0)
    @RV z ~ StateTransition2(x, y)
    placeholder(z, :z)

    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)
    code = algorithmSourceCode(algo)

    @test occursin("messages::Vector{Message}=Array{Message}(undef, 5)", code)
    @test !occursin("StateTransition", code)
end

end #module


text:
using BasisFunctions: AbstractSolverOperator, SamplingOperator
import BasisFunctions: apply, operator

export GenericAZSolver
struct GenericAZSolver <: AbstractSolverOperator
    A           ::  SynthesisOperator
    Zt          ::  AnalysisOperator
    plunge_op   ::  AbstractOperator     # The plunge operator P = (I-A*Zt)
    S           ::  SamplingOperator     # The sampling operator S
    psolver     ::  DictionaryOperator   # The low rank plunge solver for S*P*A
end

function GenericAZSolver(ap::ApproximationProblem, A;
            REG = default_regularization,
            rankestimate = 40,
            threshold = default_threshold(A),
            options...)
    pstyle = GenericOperatorStyle()
    Zt = AZ_Zt(pstyle, ap; options...)
    P = plungeoperator(pstyle, ap; options...)
    M = plungematrix(pstyle, ap; options...) # TODO some inner products are calculated twice.
    S = samplingoperator(pstyle, ap; options...)

    psolver = REG(M; threshold = threshold, rankestimate = rankestimate, options...)

    GenericAZSolver(A, Zt, P, S, psolver)
end

export default_threshold
default_threshold(A::SynthesisOperator) = regularization_threshold(coefficienttype(src(A)))

operator(op::GenericAZSolver) = op.A

function apply(op::GenericAZSolver, fun; options...)
    Pb = apply(op.S, apply(op.plunge_op, fun; options...); options...)
    x1 = op.psolver*Pb
    fun1 = x->fun(x...)-(op.A*x1)(x...)
    x2 = apply(op.Zt, fun1; options...)
    x1+x2
end


text:
module EconLH

# ToDo
# - Factor out production functions into a separate module

using DocStringExtensions
using Distributions, Random, StatsFuns

include("vectorLH.jl");
include("extreme_value_decision.jl");
include("production/productionFunctionsLH.jl");

"""
    $(SIGNATURES)

Present value of a function of time.
First value is not discounted.
"""
function present_value(f, R :: T1, T :: T2) where
    {T1 <: AbstractFloat, T2 <: Integer}

    return sum(f.(1 : T) ./ (R .^ (1:T))) .* R;
end

end # module


text:
# using Revise, Plots
using OrdinaryDiffEq, ForwardDiff, Test
	using BifurcationKit, LinearAlgebra, Parameters, Setfield
	const BK = BifurcationKit
	const FD = ForwardDiff

norminf = x -> norm(x, Inf)

function Fsl!(f, u, p, t)
	@unpack r, μ, ν, c3, c5 = p
	u1 = u[1]
	u2 = u[2]

	ua = u1^2 + u2^2

	f[1] = r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1
	f[2] = r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2
	return f
end

Fsl(x, p) = Fsl!(similar(x), x, p, 0.)

Fode(f, x, p, t) = Fsl!(f, x, p, t)
dFsl(x, dx, p) = FD.derivative(t -> Fsl(x .+ t .* dx, p), 0.)

par_sl = (r = 0.5, μ = 0., ν = 1.0, c3 = 1.0, c5 = 0.0,)
u0 = [.001, .001]

function FslMono!(f, x, p, t)
	u = x[1:2]
	du = x[3:4]
	Fsl!(f[1:2], u, p, t)
	f[3:4] .= dFsl(u, du, p)
end

####################################################################################################
par_hopf = (@set par_sl.r = 0.1)
####################################################################################################
# continuation
optconteq = ContinuationPar(ds = -0.01, detectBifurcation = 3, pMin = -0.5, nInversion = 4)
br, = continuation(Fsl, u0, par_hopf, (@lens _.r), optconteq)
####################################################################################################
prob = ODEProblem(Fode, u0, (0., 100.), par_hopf)
probMono = ODEProblem(FslMono!, vcat(u0, u0), (0., 100.), par_hopf)
####################################################################################################
sol = solve(probMono, KenCarp4(), abstol =1e-9, reltol=1e-6)
sol = solve(prob, KenCarp4(), abstol =1e-9, reltol=1e-6)
# plot(sol[1,:], sol[2,:])
####################################################################################################
section(x, T) = x[1] #* x[end]
# standard simple shooting
M = 1
dM = 1
_pb = ShootingProblem(Fsl, par_hopf, prob, KenCarp4(), 1, section; abstol =1e-10, reltol=1e-9)

initpo = [0.13, 0., 6.]
res = _pb(initpo, par_hopf)

# test of the differential of the shooting method

_dx = rand(3)
resAD = FD.derivative(z -> _pb(initpo .+ z .* _dx, par_hopf), 0.)
resFD = (_pb(initpo .+ 1e-8 .* _dx, par_hopf) - _pb(initpo, par_hopf)) .* 1e8
resAN = _pb(initpo, par_hopf, _dx; δ = 1e-8)
@test norm(resAN - resFD, Inf) < 4e-6
@test norm(resAN - resAD, Inf) < 4e-6
####################################################################################################
# test shooting interface M = 1
_pb = ShootingProblem(Fsl, par_hopf, prob, KenCarp4(), [initpo[1:end-1]]; abstol =1e-10, reltol=1e-9)
res = _pb(initpo, par_hopf)
res = _pb(initpo, par_hopf, initpo)

# test the jacobian of the functional in the case M=1
_Jad = FD.jacobian( x -> _pb(x, par_hopf), initpo)
_Jana = _pb(Val(:JacobianMatrix), initpo, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-7

_pb2 = ShootingProblem(Fsl, par_hopf, prob, Rodas4(), probMono, Rodas4(autodiff=false), [initpo[1:end-1]]; abstol =1e-10, reltol=1e-9)
res = _pb2(initpo, par_hopf)
res = _pb2(initpo, par_hopf, initpo)
BK.isSimple(_pb2)

# we test this using Newton - Continuation
ls = GMRESIterativeSolvers(reltol = 1e-5, N = length(initpo))
optn = NewtonPar(verbose = false, tol = 1e-9,  maxIter = 20, linsolver = ls)
deflationOp = BK.DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]), 1.0, [zeros(3)])
outpo, _, flag,_ = newton(_pb,
	initpo, par_hopf,
	optn,
	normN = norminf)
	@test flag

BK.getPeriod(_pb, outpo, par_hopf)
BK.getAmplitude(_pb, outpo, par_hopf)
BK.getMaximum(_pb, outpo, par_hopf)
BK.getTrajectory(_pb, outpo, par_hopf)

opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.01, ds= -0.01, pMax = 4.0, maxSteps = 30, detectBifurcation = 2, nev = 2, newtonOptions = @set optn.tol = 1e-7)
	br_pok2, upo , _= continuation(
		_pb,
		outpo, par_hopf, (@lens _.r),
		opts_po_cont;
		tangentAlgo = BorderedPred(),
		verbosity = 0,
		plot = false,
		# plotSolution = (x, p; kwargs...) -> plot!(x[1:end-1]; kwargs...),
		printSolution = (u, p) -> norm(u[1:2]),
		normC = norminf)
# plot(br_pok2)
####################################################################################################
# test automatic branch switching
JFsl(x, p) = FD.jacobian(t -> Fsl(t, p), x)
d1Fsl(x, p, dx) = FD.derivative(t -> Fsl(x .+ t .* dx, p), 0.)
d2Fsl(x, p, dx1, dx2) = FD.derivative(t -> d1Fsl(x .+ t .* dx2, p, dx1), 0.)
d3Fsl(x, p, dx1, dx2, dx3) = FD.derivative(t -> d2Fsl(x .+ t .* dx3, p, dx1, dx2), 0.)
jet = (Fsl, JFsl, d2Fsl, d3Fsl)

br_pok2, = continuation(jet..., br, 1, opts_po_cont, ShootingProblem(1, par_hopf, prob, KenCarp4());printSolution = (u, p) -> norm(u[1:2]), normC = norminf)

# test matrix-free computation of floquet coefficients
eil = EigKrylovKit(dim = 2, x₀=rand(2))
opts_po_contMF = @set opts_po_cont.newtonOptions.eigsolver = eil
opts_po_contMF = @set opts_po_cont.detectBifurcation=0
br_pok2, = continuation(jet...,br,1, opts_po_contMF, ShootingProblem(1, par_hopf, prob, Rodas4());printSolution = (u, p) -> norm(u[1:2]), normC = norminf, plot=false)
####################################################################################################
# test shooting interface M > 1
_pb = ShootingProblem(Fsl, par_hopf, prob, KenCarp4(), [initpo[1:end-1],initpo[1:end-1],initpo[1:end-1]]; abstol =1e-10, reltol=1e-9)
initpo = [0.13, 0, 0, 0.13, 0, 0.13 , 6.3]
res = _pb(initpo, par_hopf)
res = _pb(initpo, par_hopf, initpo)
# test the jacobian of the functional in the case M=1
_Jad = FD.jacobian( x -> _pb(x, par_hopf), initpo)
_Jana = _pb(Val(:JacobianMatrix), initpo, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-7
####################################################################################################
# Single Poincaré Shooting with hyperplane parametrization
normals = [[-1., 0.]]
centers = [zeros(2)]

probPsh = PoincareShootingProblem(2, par_hopf,prob, Rodas4(), probMono, Rodas4(); abstol=1e-10, reltol=1e-9)
probPsh = PoincareShootingProblem(2, par_hopf,prob, Rodas4(); rtol = abstol=1e-10, reltol=1e-9)

probPsh = PoincareShootingProblem(Fsl, par_hopf,
		prob, Rodas4(),
		probMono, Rodas4(),
		normals, centers; rtol = abstol =1e-10, reltol=1e-9)

initpo_bar = BK.R(probPsh, [0, 0.4], 1)

BK.E(probPsh, [1.0,], 1)
initpo_bar = [0.4]

probPsh(initpo_bar, par_hopf)

# test of the analytical formula for jacobian of the functional
_Jad = BifurcationKit.finiteDifferences( x-> probPsh(x, par_hopf), initpo_bar)
_Jana = probPsh(Val(:JacobianMatrix), initpo_bar, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-4

ls = GMRESIterativeSolvers(reltol = 1e-7, N = length(initpo_bar), maxiter = 500, verbose = false)
	eil = EigKrylovKit(dim = 1, x₀=rand(1))
	optn = NewtonPar(verbose = false, tol = 1e-8,  maxIter = 140, linsolver = ls, eigsolver = eil)
	deflationOp = BK.DeflationOperator(2.0, dot, 1.0, [zero(initpo_bar)])
	outpo, _, flag, = newton(probPsh,
			initpo_bar, par_hopf,
			optn; normN = norminf)
	@test flag

BK.getPeriod(probPsh, outpo, par_hopf)
BK.getAmplitude(probPsh, outpo, par_hopf)
BK.getMaximum(probPsh, outpo, par_hopf)
BK.getTrajectory(probPsh, outpo, par_hopf)

probPsh = PoincareShootingProblem(Fsl, par_hopf,
		prob, Rodas4(),
		# probMono, Rodas4(autodiff=false),
		normals, centers; rtol = 1e-8)

opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.015, ds= -0.01, pMax = 4.0, maxSteps = 50, newtonOptions = setproperties(optn;tol = 1e-9, eigsolver = eil), detectBifurcation = 1)
	br_pok2, = continuation(
		probPsh, outpo, par_hopf, (@lens _.r),
		opts_po_cont; verbosity = 0,
		tangentAlgo = BorderedPred(),
		plot = false,
		# plotSolution = (x, p;kwargs...) -> plot!(x; kwargs...),
		printSolution = (u, p) -> norm(u), normC = norminf)
# plot(br_pok2)
####################################################################################################
# normals = [[-1., 0.], [1, -1]]
# centers = [zeros(2), zeros(2)]
# initpo_bar = [1.04, -1.04/√2]

normals = [[-1., 0.], [1, 0]]
centers = [zeros(2), zeros(2)]
initpo_bar = [0.2, -0.2]

probPsh = PoincareShootingProblem(Fsl, par_hopf, prob, KenCarp4(), normals, centers; abstol=1e-10, reltol=1e-9)
# version with analytical jacobian
probPsh2 = PoincareShootingProblem(Fsl, par_hopf, prob, KenCarp4(), normals, centers; abstol=1e-10, reltol=1e-9, δ = 0)

# test of the analytical formula for jacobian of the functional
_Jad = BifurcationKit.finiteDifferences( x-> probPsh(x, par_hopf), initpo_bar)
# _Jphifd = BifurcationKit.finiteDifferences(x->probPsh.flow(x, par_hopf, Inf64), [0,0.4]; δ=1e-8)
_Jana = probPsh(Val(:JacobianMatrix), initpo_bar, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-5

ls = GMRESIterativeSolvers(reltol = 1e-5, N = length(initpo_bar), maxiter = 500, verbose = false)
	eil = EigKrylovKit(dim = 1, x₀=rand(1))
	optn = NewtonPar(verbose = false, tol = 1e-9,  maxIter = 140, linsolver = ls, eigsolver = eil)
	deflationOp = DeflationOperator(2.0, dot, 1.0, [zero(initpo_bar)])
	outpo, = newton(probPsh2, initpo_bar, par_hopf, optn; normN = norminf)
	outpo, = newton(probPsh, initpo_bar, par_hopf, optn; normN = norminf)

getPeriod(probPsh, outpo, par_hopf)

opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.025, ds= -0.01, pMax = 4.0, maxSteps = 50, newtonOptions = (@set optn.tol = 1e-9), detectBifurcation = 1, nev = 2)
	br_pok2, = continuation(probPsh, outpo, par_hopf, (@lens _.r),
		opts_po_cont; verbosity = 0,
		tangentAlgo = BorderedPred(),
		plot = false,
		# plotSolution = (x, p;kwargs...) -> plot!(x; kwargs...),
		printSolution = (u, p) -> norm(u), normC = norminf)

####################################################################################################
normals = [[-1., 0.], [1, 0], [0, 1]]
centers = [zeros(2), zeros(2), zeros(2)]
initpo = [[0., 0.4], [0, -.3], [0.3, 0]]

probPsh = PoincareShootingProblem(Fsl, par_hopf, prob, KenCarp4(), normals, centers; abstol=1e-10, reltol=1e-9)

initpo_bar = reduce(vcat, [BK.R(probPsh, initpo[ii], ii) for ii in eachindex(centers)])
# same with projection function
initpo_bar = reduce(vcat, BK.projection(probPsh, initpo))

# test of the other projection function
BK.projection(probPsh, reduce(hcat, initpo)')

probPsh(initpo_bar, par_hopf; verbose = true)

# test of the analytical formula for jacobian of the functional
_Jad = BifurcationKit.finiteDifferences( x-> probPsh(x, par_hopf), initpo_bar)
_Jana = probPsh(Val(:JacobianMatrix), initpo_bar, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-5

ls = GMRESIterativeSolvers(reltol = 1e-7, N = length(initpo_bar), maxiter = 10, verbose = false)
	optn = NewtonPar(verbose = false, tol = 1e-9,  maxIter = 50, linsolver = ls, eigsolver = eil)
	outpo, = newton(probPsh, initpo_bar, par_hopf, optn; normN = norminf)

for ii=1:length(normals)
	@show BK.E(probPsh	, [outpo[ii]], ii)
end

getPeriod(probPsh, outpo, par_hopf)

opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.025, ds= -0.005, pMax = 4.0, maxSteps = 10, newtonOptions = setproperties(optn; tol = 1e-8), detectBifurcation = 1)
	br_hpsh, upo , _= continuation(probPsh, outpo, par_hopf, (@lens _.r),
		opts_po_cont; normC = norminf)
# plot(br_hpsh)
####################################################################################################
# test automatic branch switching with most possible options
# calls with analytical jacobians
br_psh, = continuation(jet..., br,1, (@set opts_po_cont.ds = 0.005), PoincareShootingProblem(1, par_hopf, prob, KenCarp4(), probMono, KenCarp4(); abstol=1e-10, reltol=1e-9); normC = norminf)

opts_po_cont = @set opts_po_cont.detectBifurcation = 0
for M in [1,2], linearPO in [:autodiffMF, :MatrixFree, :autodiffDense, :FiniteDifferencesDense]
	@show M, linearPO
	br_psh, = continuation(jet..., br, 1, (@set opts_po_cont.ds = 0.005), PoincareShootingProblem(M, par_hopf, prob, Rodas4P(); abstol=1e-10, reltol=1e-9, parallel = true); normC = norminf, updateSectionEveryStep = 2, linearPO = linearPO == :autodiffMF ? :FiniteDifferencesDense : linearPO, verbosity = 0)

	br_ssh, = continuation(jet..., br, 1, (@set opts_po_cont.ds = 0.005), ShootingProblem(M, par_hopf, prob, Rodas4P(); abstol=1e-10, reltol=1e-9, parallel = true); normC = norminf, updateSectionEveryStep = 2, linearPO = linearPO, verbosity = 0)
end

br_psh, = continuation(jet..., br,1, (@set opts_po_cont.ds = 0.005), PoincareShootingProblem(2, par_hopf, prob, KenCarp4(), probMono, KenCarp4(); abstol=1e-10, reltol=1e-9); normC = norminf)

# plot(br_hpsh, vars = (:step, :p))
# plot(br_hpsh)


text:
# Extends PowerModels/src/core/types.jl

##### Linear Approximations #####

"""
Linearized AC branch flow model for radial networks.

Variables:
- squared voltage magnitude;
- branch active power;
- branch reactive power.

Properties:
- same voltage angle for all buses;
- lossless.

Differences with respect to `BFAPowerModel`:
- shunt admittances of the branches are neglected;
- the complex power in the thermal limit constraints of the branches is limited by an octagon
  instead of a circle, so as to keep the model linear.
"""
mutable struct BFARadPowerModel <: _PM.AbstractBFAModel _PM.@pm_fields end


text:
using Base:Symbol
"""
    read_weather(file[,args...];
        date_format = DateFormat("yyyy-mm-ddTHH:MM:SS.s"),
        hour_format = DateFormat("HH:MM:SS")
    )

Read a meteo file. The meteo file is a CSV, and optionnaly with metadata in a header formatted
as a commented YAML. The column names should match exactly the fields of [`Atmosphere`](@ref), or
the user should provide their transformation as arguments (`args`) to help mapping the two. The
transformations are given as for `DataFrames`.

# Arguments

- `file::String`: path to a meteo file
- `var_names = Dict()`: A Dict to map the file variable names to the Atmosphere variable names
- `date_format = DateFormat("yyyy-mm-ddTHH:MM:SS.s")`: the format for the `DateTime` columns
- `hour_format = DateFormat("HH:MM:SS")`: the format for the `Time` columns (*e.g.* `hour_start`)

# Examples

```julia
using Dates

file = joinpath(dirname(dirname(pathof(PlantBiophysics))),"test","inputs","meteo.csv")

meteo = read_weather(
    file,
    :temperature => :T,
    :relativeHumidity => (x -> x ./100) => :Rh,
    :wind => :Wind,
    :atmosphereCO2_ppm => :Cₐ,
    date_format = DateFormat("yyyy/mm/dd")
)
```
"""
function read_weather(
    file,args...;
    date_format = DateFormat("yyyy-mm-ddTHH:MM:SS.s"),
    hour_format = DateFormat("HH:MM:SS")
    )

    arguments = (args...,)
    data, metadata = read_weather(file, DataFrame)

    # Clean-up the variable names:
    length(arguments) > 0 && select!(data, arguments...)

    # If there's a "use" field in the YAML, parse it and rename it:
    if haskey(metadata, "use")
        splitted_use = split(metadata["use"], r"[,\s]")
        metadata["use"] = Symbol.(splitted_use[findall(x -> length(x) > 0, splitted_use)])

        orig_names = [i.first for i in arguments]
        new_names = [isa(i.second, Pair) ? i.second.second : i.second for i in arguments]
        length(arguments) > 0 && replace!(metadata["use"], Pair.(orig_names, new_names)...)
    end
    # NB: the "use" field is not used in PlantBiophysics, but it is still correctly parsed.

    compute_date!(data, date_format, hour_format)
    compute_duration!(data, hour_format)

    cols = fieldnames(Atmosphere)
    select!(data, names(data, x -> Symbol(x) in cols))

    Weather(data, metadata)
end

function read_weather(file, ::Type{DataFrame})
    yaml_data = open(file, "r") do io
        yaml_data = ""
        is_yaml = true
        while is_yaml
            line = readline(io, keep = true)
            if line[1:2] == "#'"
                yaml_data *= lstrip(line[3:end])
            else
                is_yaml = false
            end
        end
        return yaml_data
    end

    metadata = YAML.load(yaml_data)
    push!(metadata, "file" => file)

    met_data = CSV.read(file, DataFrame; comment = "#")

    (data = met_data, metadata = metadata)
end

"""
    compute_date!(df, date_format, hour_format)

Compute the `date` column depending on several cases:

- If it is already in the DataFrame and is a `DateTime`, does nothing.
- If it is a `String`, tries and parse it using a user-input `DateFormat`
- If it is a `Date`, return it as is, or try to make it a `DateTime` if there's a column named
`hour_start`
"""
function compute_date!(df, date_format, hour_format)
        if hasproperty(df, :date) && typeof(df.date[1]) != DateTime
        # There's a "date" column but it is not a DateTime
        # Trying to parse it with the user-defined format:
        try
            df.date = Dates.DateTime.(df.date, date_format)
        catch
            error(
                "The values in the `date` column cannot be parsed.",
                " Please check the format of the dates or provide the format as argument."
            )
        end

        if typeof(df.date[1]) == Dates.Date && hasproperty(df, :hour_start)
            # The `date` column is of Date type, we have to add the Time if there's a column named
            # `hour_start`:
            if typeof(df.hour_start[1]) != Dates.Time
                # There's a "hour_start" column but it is not of Time type
                # If it is a String, it did not parse at reading with CSV, so trying to use
                # the user-defined format:
                try
                    df.hour_start = Dates.Time.(df.hour_start, hour_format)
                catch
                    error(
                        "The values in the `hour_start` column cannot be parsed.",
                        " Please check the format of the hours or provide the format as argument."
                    )
                end
            end
            # Adding the Time to the Date to make a DateTime:
            df.date = df.date .+ df.hour_start
        end
    end
end

"""
    compute_duration!(df, hour_format)

Compute the `duration` column depending on several cases:

- If it is already in the DataFrame, does nothing.
- If it is not, but there's a column named `hour_end` and another one either called `hour_start`
or `date`, compute the duration from the period between `hour_start` (or `date`) and `hour_end`.
"""
function compute_duration!(df, hour_format)
        # `duration` is not in the df but there is an `hour_end` column:
    if hasproperty(df, :hour_end) && !hasproperty(df, :duration)
        if typeof(df.hour_end[1]) != Dates.Time
            # There's a `hour_end` column but it is not of Time type
            # If it is a String, it did not parse at reading with CSV, so trying to use
            # the user-defined format:
            if typeof(df.hour_end[1]) != String
                try
                    df.hour_end = Dates.Time.(df.hour_end, hour_format)
                catch
                    error(
                            "The values in the `hour_end` column cannot be parsed.",
                            " Please check the format of the hours or provide the format as argument."
                        )
                end
            end

            # If it is of Time type, transform it into a DateTime:
            if typeof(df.hour_end[1]) == dates.DateTime
                df.hour_end = Dates.Time(df.hour_end)
            end
        end

        if hasproperty(df, :hour_start)
            df.duration = Minute.(df.hour_end .- df.hour_start)
        elseif hasproperty(df, :date)
            df.duration = Minute.(df.hour_end .- Time.(df.date))
        end
    end
end