text:
export AbstractProblem,
       Problem,
       TrackingType,
       AffineTracking,
       ProjectiveTracking,
       homotopy,
       homogenization,
       embed,
       homvars,
       problem_startsolutions

const problem_startsolutions_supported_keywords = [
    [
     :seed,
     :homvar,
     :homvars,
     :variable_groups,
     :variable_ordering,
     :homotopy,
     :system,
     :system_scaling,
     :affine_tracking,
     :projective_tracking,
     :start_system,
     :only_torus,
    ]
    input_supported_keywords
]

const DEFAULT_SYSTEM = @static VERSION < v"1.1" ? FPSystem : SPSystem
const DEFAULT_HOMOTOPY = StraightLineHomotopy
const DEFAULT_SYSTEM_SCALING = :equations

"""
	TrackingType

Abstract type determining how paths are tracked.
"""
abstract type TrackingType end

"""
	AffineTracking <: TrackingType

Indicates that paths should be tracked in affine space.
"""
struct AffineTracking <: TrackingType end

"""
	ProjectiveTracking <: TrackingType

Indicates that paths should be tracked in projective space.
"""
struct ProjectiveTracking <: TrackingType end

abstract type AbstractProblem{T<:TrackingType,VG<:VariableGroups} end
Base.broadcastable(P::AbstractProblem) = Ref(P)

"""
    homotopy(prob::AbstractProblem)

Get the homotopy stored in the problem `prob`.
"""
homotopy(prob::AbstractProblem) = prob.homotopy

"""
    seed(prob::AbstractProblem)

Get the random seed used for the problem `prob`.
"""
seed(prob::AbstractProblem) = prob.seed

"""
    Problem(T::TrackingType, H::AbstractHomotopy, homogenization::AbstractHomogenization, seed::Int)

Construct a `Problem`. If `T <: ProjectiveTracking` then the homotopy `H` needs to be homogeneous.
"""
struct Problem{T<:TrackingType,VG<:VariableGroups} <: AbstractProblem{T,VG}
    tracking_type::T
    homotopy::AbstractHomotopy
    vargroups::VG
    seed::Int
    startsolutions_need_reordering::Bool
    regauging_factors::Union{Nothing,Vector{Float64}}
end
function Problem(
    T::TrackingType,
    H::AbstractHomotopy,
    vargroups::VariableGroups,
    seed::Int;
    startsolutions_need_reordering = false,
    regauging_factors = nothing,
)
    if isa(T, AffineTracking) && startsolutions_need_reordering
        error("Affine tracking doesn't support reodering.")
    end
    Problem(T, H, vargroups, seed, startsolutions_need_reordering, regauging_factors)
end
function Problem(
    T::TrackingType,
    G::AbstractSystem,
    F::AbstractSystem,
    vargroups::VariableGroups,
    seed::Int;
    homotopy = DEFAULT_HOMOTOPY,
    kwargs...,
)
    Problem(T, homotopy(G, F), vargroups, seed; kwargs...)
end
Problem{T}(args...; kwargs...) where {T<:TrackingType} = Problem(T(), args...; kwargs...)

struct PolyhedralProblem{T<:TrackingType,VG<:VariableGroups} <: AbstractProblem{T,VG}
    tracking_type::T
    toric_homotopy::AbstractHomotopy
    generic_homotopy::AbstractHomotopy
    vargroups::VG
    seed::Int
    startsolutions_need_reordering::Bool
    regauging_factors::Union{Nothing,Vector{Float64}}
end
function PolyhedralProblem(
    T::TrackingType,
    toric_homotopy::ToricHomotopy,
    generic_homotopy::AbstractHomotopy,
    vargroups::VariableGroups,
    seed::Int;
    startsolutions_need_reordering = false,
    regauging_factors = nothing,
)
    if isa(T, AffineTracking) && startsolutions_need_reordering
        error("Affine tracking doesn't support reodering.")
    end
    PolyhedralProblem(
        T,
        toric_homotopy,
        generic_homotopy,
        vargroups,
        seed,
        startsolutions_need_reordering,
        regauging_factors,
    )
end
PolyhedralProblem{T}(args...; kwargs...) where {T<:TrackingType} =
    Problem(T(), args...; kwargs...)

struct OverdeterminedProblem{
    T<:TrackingType,
    VG<:VariableGroups,
    AP<:AbstractProblem{T,VG},
} <: AbstractProblem{T,VG}
    problem::AP
    target_system::AbstractSystem
    regauging_factors::Union{Nothing,Vector{Float64}}
end
function OverdeterminedProblem(
    problem::AbstractProblem,
    target_system::AbstractSystem;
    regauging_factors = nothing,
)
    OverdeterminedProblem(problem, target_system, regauging_factors)
end

"""
    homvars(prob::AbstractProblem)

Get the homogenization variables of the problem. Returns `nothing` if there are no.
"""
function homvars(prob::AbstractProblem{ProjectiveTracking})
    if has_dedicated_homvars(prob.vargroups)
        map(last, prob.vargroups.groups)
    else
        nothing
    end
end
homvars(prob::AbstractProblem{AffineTracking}) = nothing
homvars(prob::OverdeterminedProblem) = homvars(prob.problem)

"""
    embed(prob::Problem{ProjectiveTracking}, v)

Embed the vector `v` into projective space if necessary.
"""
function embed!(x, prob::AbstractProblem{ProjectiveTracking}, v)
    if prob.startsolutions_need_reordering
        embed_projective!(x, prob.vargroups, v)
    else
        ProjectiveVectors.embed!(x, v)
    end
end
embed!(x, prob::OverdeterminedProblem{ProjectiveTracking}, v) = embed!(x, prob.problem, v)
embed!(x, prob::AbstractProblem{ProjectiveTracking}, v::PVector) =
    ProjectiveVectors.embed!(x, v)
embed!(x, prob::AbstractProblem{AffineTracking}, v::AbstractVector) = (x .= v; x)

function embed(prob::AbstractProblem{ProjectiveTracking}, v)
    if prob.startsolutions_need_reordering
        embed_projective(prob.vargroups, v)
    else
        ProjectiveVectors.embed(v, projective_dims(prob.vargroups))
    end
end
embed(prob::OverdeterminedProblem{ProjectiveTracking}, v) = embed(prob.problem, v)
embed(prob::AbstractProblem{ProjectiveTracking}, v::PVector) = v
embed(prob::AbstractProblem{AffineTracking}, v::AbstractVector) = v

tracking_vector_type(prob::AbstractProblem{AffineTracking}) = Vector{ComplexF64}
function tracking_vector_type(prob::AbstractProblem{
    ProjectiveTracking,
    <:VariableGroups{N},
}) where {N}
    PVector{ComplexF64,N}
end

"""
    pull_back(prob::Problem, x)

Pull the solution `x` into affine space if necessary. Creates a copy.
"""
function pull_back(prob::AbstractProblem{AffineTracking}, x::AbstractVector; regauge = true)
    if regauge
        regauge!(copy(x), prob)
    else
        copy(x)
    end
end
function pull_back(prob::OverdeterminedProblem{ProjectiveTracking}, x::PVector; kwargs...)
    pull_back(prob.problem, x; kwargs...)
end
function pull_back(prob::AbstractProblem{ProjectiveTracking}, x::PVector; regauge = true)
    if pull_back_is_to_affine(prob)
        λ = prob.regauging_factors
        if λ === nothing || !regauge
            map(kᵢ -> x[kᵢ[1]] / x[kᵢ[2]], prob.vargroups.pull_back_mapping)
        else
            map(
                kᵢ -> (λ[kᵢ[1]] * x[kᵢ[1]]) / (λ[kᵢ[2]] * x[kᵢ[2]]),
                prob.vargroups.pull_back_mapping,
            )
        end
    elseif regauge
        LinearAlgebra.normalize!(regauge!(copy(x), prob))
    else
        LinearAlgebra.normalize(x)
    end
end

function regauge!(x, prob::AbstractProblem)
    if prob.regauging_factors !== nothing
        x .= prob.regauging_factors .* x
    end
    x
end

"""
    pull_back_is_to_affine(prob::ProjectiveProblem)

Returns `true` if [`pull_back`](@ref) would pull a solution `x` into affine space.
"""
pull_back_is_to_affine(prob::AbstractProblem{AffineTracking}) = true
function pull_back_is_to_affine(prob::AbstractProblem{ProjectiveTracking})
    has_dedicated_homvars(prob.vargroups)
end
function pull_back_is_to_affine(prob::OverdeterminedProblem{ProjectiveTracking})
    pull_back_is_to_affine(prob.problem)
end

function construct_system(F::Composition, system_constructor; homvars = nothing, kwargs...)
    CompositionSystem(F, system_constructor; homvars = homvars, kwargs...)
end
function construct_system(F::MPPolys, system_constructor; homvars = nothing, kwargs...)
    system_constructor(F; kwargs...)
end

function apply_system_scaling(F, vars, system_scaling::Union{Nothing,Symbol,Bool})
    if system_scaling == :equations_and_variables || system_scaling == true
        precondition(F, vars)
    elseif system_scaling == :equations
        normalize_coefficients(F), nothing
    elseif system_scaling === nothing || system_scaling == false
        F, nothing
    else
        throw(ArgumentError("Got unsupported argument `system_scaling=$(system_scaling)`." *
                            " Valid values are `nothing`, `:equations` and `:equations_and_variables`."))
    end
end

function default_affine_tracking(F::TargetSystemInput{<:MPPolyInputs}, hominfo)
    !(is_homogeneous(F.system, hominfo))
end
function default_affine_tracking(
    F::TargetSystemInput{<:MPPolyInputs},
    hominfo::HomogenizationInformation,
)
    is_hom = is_homogeneous(F.system, hominfo)
    if !is_hom && !isnothing(hominfo.homvars)
        throw(ArgumentError("Input system is not homogeneous although `homvars=$(hominfo.homvars)` was passed."))
    end

    !(is_hom || (hominfo.vargroups !== nothing && length(hominfo.vargroups) > 1))
end
function default_affine_tracking(input::StartTargetInput{<:MPPolyInputs}, hominfo)
    F, G = input.target, input.start
    !(is_homogeneous(F, hominfo) && is_homogeneous(G, hominfo))
end
function default_affine_tracking(input::HomotopyInput, hominfo)
    !(is_homogeneous(input.H))
end
function default_affine_tracking(F::ParameterSystemInput{<:MPPolyInputs}, hominfo)
    !(is_homogeneous(F.system, hominfo; parameters = F.parameters))
end
function default_affine_tracking(input::ParameterSystemInput{<:AbstractSystem}, hominfo)
    H = ParameterHomotopy(
        input.system,
        p₁ = input.p₁,
        p₀ = input.p₀,
        γ₁ = input.γ₁,
        γ₀ = input.γ₀,
    )
    !(is_homogeneous(H))
end

is_homogeneous(F::AbstractSystem) = compute_numerically_degrees(F) !== nothing
is_homogeneous(H::AbstractHomotopy) = is_homogeneous(FixedHomotopy(H, rand()))

"""
    problem_startsolutions(F; options...)
    problem_startsolutions(G, F, startsolutions; options...)
    problem_startsolutions(H::AbstractHomotopy, startsolutions; options...)
    problem_startsolutions(prob::TotalDegreeProblem; options...)
    problem_startsolutions(prob::StartTargetProblem; options...)

    Construct the problem and if necessary the startsolutions. This steps
    constructs a homotopy and homogenizes the systems if necessary.

    The `options` are
    * `seed::Int`: Random seed used in the construction.
    * `system=SPSystem`: A constructor to assemble a [`AbstractSystem`](@ref). The constructor
    is called with `system(polys, variables)` where `variables` determines the variable ordering.
    * `homotopy=StraightLineHomotopy`: A constructor to construct a [`AbstractHomotopy`](@ref) an `AbstractSystem`. The constructor
    is called with `homotopy(start, target)` where `start` and `target` are systems constructed
    with `system`.
	* `affine_tracking::Bool=true`: Indicate whether path tracking should happen in affine space.
	* `projective_tracking::Bool=false`: Indicate whether path tracking should happen in projective space. The flag `affine_tracking` is dominant.
"""
function problem_startsolutions end

function problem_startsolutions(
    args...;
    seed = randseed(),
    variable_ordering = nothing,
    kwargs...,
)
    Random.seed!(seed)
    supported, rest = splitkwargs(kwargs, input_supported_keywords)
    input, startsolutions = input_startsolutions(
        args...;
        variables = variable_ordering,
        supported...,
    )
    if variable_ordering !== nothing
        problem_startsolutions(
            input,
            startsolutions,
            seed;
            variables = variable_ordering,
            rest...,
        )
    else
        problem_startsolutions(input, startsolutions, seed; rest...)
    end
end

function problem_startsolutions(
    input::AbstractInput,
    startsolutions = nothing;
    seed = randseed(),
    kwargs...,
)
    problem_startsolutions(input, startsolutions, seed; kwargs...)
end
function problem_startsolutions(
    input::AbstractInput,
    startsolutions,
    seed::Int;
    homvar = nothing,
    homvars = nothing,
    variable_groups = nothing,
    affine_tracking = nothing,
    projective_tracking = nothing,
    kwargs...,
)

    homvar_info = HomogenizationInformation(
        ;
        homvar = homvar,
        homvars = homvars,
        variable_groups = variable_groups,
    )

    #projective_tracking is for the frontend
    #internally, we use affine_tracking
    if isnothing(affine_tracking)
        if isnothing(projective_tracking)
            problem_startsolutions(input, startsolutions, homvar_info, seed; kwargs...)
        else
            problem_startsolutions(
                input,
                startsolutions,
                homvar_info,
                seed;
                affine_tracking = !projective_tracking,
                kwargs...,
            )
        end
    else
        problem_startsolutions(
            input,
            startsolutions,
            homvar_info,
            seed;
            affine_tracking = affine_tracking,
            kwargs...,
        )
    end
end

function problem_startsolutions(
    input::HomotopyInput,
    startsolutions,
    homvar_info,
    seed;
    affine_tracking = default_affine_tracking(input, homvar_info),
    system_scaling = nothing,
    variables = nothing,
    kwargs...,
)
    if !affine_tracking
        is_homogeneous(input.H) || throw(ArgumentError("Input homotopy is not homogeneous by our numerical check."))
        if homvar_info === nothing
            N = size(input.H)[2]
            if start_solution_sample(startsolutions) isa PVector
                vargroups = VariableGroups(N)
            else
                vargroups = VariableGroups(N, N)
            end
            Problem{ProjectiveTracking}(input.H, vargroups, seed), startsolutions
        else
            vargroups = VariableGroups(size(input.H)[2], homvar_info)
            Problem{ProjectiveTracking}(input.H, vargroups, seed), startsolutions
        end
    else
        Problem{AffineTracking}(input.H, VariableGroups(size(input.H)[2]), seed),
        startsolutions
    end
end

##############
# TOTALDEGREE
##############

function problem_startsolutions(
    input::TargetSystemInput{<:MPPolyInputs},
    ::Nothing,
    homvar_info,
    seed;
    start_system = :total_degree,
    affine_tracking = default_affine_tracking(input, homvar_info),
    system_scaling = DEFAULT_SYSTEM_SCALING,
    system = DEFAULT_SYSTEM,
    only_torus = false,
    variables = variables(input.system),
    kwargs...,
)
    if affine_tracking
        vargroups = VariableGroups(variables)
        homvars = nothing
        F = input.system
        tracking_type = AffineTracking()
    else
        F, vargroups, homvars = homogenize_if_necessary(
            input.system,
            homvar_info;
            vars = variables,
        )
        tracking_type = ProjectiveTracking()
    end

    classifcation = classify_system(F, vargroups; affine_tracking = affine_tracking)
    if classifcation == :underdetermined
        throw(ArgumentError("Underdetermined polynomial systems are currently not supported." *
                            " Consider adding linear polynomials to your system in order to reduce your system" *
                            " to a zero dimensional system."))
# The following case is too annoying right now
    elseif classifcation == :overdetermined && ngroups(vargroups) > 1
        error(ArgumentError("Overdetermined polynomial systems with a multi-homogeneous" *
                            " structure are currently not supported."))
    end

    vars = flattened_variable_groups(vargroups)
    target_constructor(f) = construct_system(f, system; variables = vars, homvars = homvars)

    # Scale systems
    f, regauging_factors = apply_system_scaling(F, vars, system_scaling)

    if ngroups(vargroups) > 1 # Multihomogeneous
        affine_tracking && error("Affine tracking is currently not supported for variable groups.")

        D = multidegrees(F, vargroups)
        C = multi_bezout_coefficients(D, projective_dims(vargroups))

        g = MultiHomTotalDegreeSystem(D, C)
        problem = Problem(
            tracking_type,
            g,
            target_constructor(f),
            vargroups,
            seed;
            regauging_factors = regauging_factors,
            kwargs...,
        )
        startsolutions = totaldegree_solutions(g, vargroups)
    # We have ngroups(vargroups) == 1
    elseif start_system == :total_degree
        n = affine_tracking ? length(vars) : length(vars) - 1
        degrees = maxdegrees(F)

        if classifcation == :overdetermined
            perm = sortperm(degrees; rev = true)
            # reorder polynomials by minimal degree
            f = permute(f, perm)
            degrees = degrees[perm]
        end

        g = TotalDegreeSystem(degrees[1:n]; affine = affine_tracking)

        if classifcation == :overdetermined
            A = randn(ComplexF64, n, npolynomials(f) - n)
            f̄ = target_constructor(f)
            problem = OverdeterminedProblem(
                Problem(
                    tracking_type,
                    g,
                    SquaredUpSystem(f̄, A, degrees),
                    vargroups,
                    seed;
                    kwargs...,
                ),
                f̄;
                regauging_factors = regauging_factors,
            )
        # square case
        else
            problem = Problem(
                tracking_type,
                g,
                target_constructor(f),
                vargroups,
                seed;
                regauging_factors = regauging_factors,
                kwargs...,
            )
        end
        startsolutions = totaldegree_solutions(g, vargroups)
    elseif start_system == :polyhedral
        n = affine_tracking ? length(vars) : length(vars) - 1
        degrees = maxdegrees(F)
        if f isa Composition
            f = expand(f)
        end

        if classifcation == :overdetermined
            # square the system up
            # We do not automatically sort the polynomials by degree since it is not
            # clear that this is actually always advantagous
            f̂ = [LinearAlgebra.I randn(ComplexF64, n, length(f) - n)] * f
        else
            f̂ = f
        end

        support, coeffs = support_coefficients(f̂, vars)
        if !affine_tracking
            affine_support = map(A -> A[1:end-1, :], support)
        else
            affine_support = support
        end

        if !only_torus
            for i = 1:length(support)
                !iszero(@view affine_support[i][:, end]) || continue
                affine_support[i] = [affine_support[i] zeros(Int32, n)]
                push!(coeffs[i], zero(eltype(coeffs[i])))
                if !affine_tracking
                    support[i] = [support[i] [zeros(Int32, n); degrees[i]]]
                end
            end
        end
        cell_iter = PolyhedralStartSolutionsIterator(affine_support)
        toric_homotopy = ToricHomotopy(affine_support, cell_iter.start_coefficients)
        generic_homotopy = CoefficientHomotopy(
            support,
            cell_iter.start_coefficients,
            coeffs,
        )

        if classifcation == :overdetermined
            # always use fp system for the target system since
            problem = OverdeterminedProblem(
                PolyhedralProblem(
                    tracking_type,
                    toric_homotopy,
                    generic_homotopy,
                    vargroups,
                    seed;
                    kwargs...,
                ),
                FPSystem(f),
                regauging_factors = regauging_factors,
            )
        else
            problem = PolyhedralProblem(
                tracking_type,
                toric_homotopy,
                generic_homotopy,
                vargroups,
                seed;
                regauging_factors = regauging_factors,
                kwargs...,
            )
        end
        startsolutions = cell_iter
    else
        throw(ArgumentError("Unsupported argument `start_system=$start_system`. " *
                            "Possible values are `:total_degree` and `:polyhedral`"))
    end

    problem, startsolutions
end

function problem_startsolutions(
    input::TargetSystemInput{<:AbstractSystem},
    ::Nothing,
    hominfo::Union{Nothing,HomogenizationInformation},
    seed;
    system = DEFAULT_SYSTEM,
    system_scaling = nothing,
    variables = nothing,
    kwargs...,
)
    n, N = size(input.system)

    degrees, is_homogeneous = degrees_ishomogeneous(input.system)
    G = TotalDegreeSystem(degrees; affine = !is_homogeneous)
    variable_groups = VariableGroups(N, hominfo)
    classification = classify_system(
        input.system,
        variable_groups;
        affine_tracking = !is_homogeneous,
    )

    # overdetermined
    if classification == :overdetermined
        m = N - !is_homogeneous
        A = randn(ComplexF64, m, n - m)
        f = SquaredUpSystem(input.system, A, degrees)
        problem = OverdeterminedProblem(
            Problem{is_homogeneous ? ProjectiveTracking : AffineTracking}(
                G,
                f,
                variable_groups,
                seed;
                kwargs...,
            ),
            input.system,
        )
        starts = totaldegree_solutions(max.(degrees[1:m], maximum(degrees[m+1:end])))
    # underdetermined
    elseif classification == :underdetermined
        throw(ArgumentError("Underdetermined polynomial systems are currently not supported." *
                            " Consider adding linear polynomials to your system in order to reduce your system" *
                            " to a zero dimensional system."))
    else # square
        problem = Problem{is_homogeneous ? ProjectiveTracking : AffineTracking}(
            G,
            input.system,
            variable_groups,
            seed;
            kwargs...,
        )
        starts = totaldegree_solutions(degrees)
    end
    # Check overdetermined case
    problem, starts
end

function degrees_ishomogeneous(F)
    n, N = size(F)

    degs = compute_numerically_degrees(F)
    is_homogeneous = !isnothing(degs)
    if isnothing(degs)
        degs = degrees(F)
    end
    isnothing(degs) && throw(ArgumentError("Cannot compute degrees of the input system. Consider overloading `system_degree(F)."))

    degs, is_homogeneous
end

###############
# START TARGET
###############

function problem_startsolutions(
    input::StartTargetInput{<:MPPolyInputs,<:MPPolyInputs},
    startsolutions,
    homvar,
    seed;
    affine_tracking = default_affine_tracking(input, homvar),
    system_scaling = DEFAULT_SYSTEM_SCALING,
    system = DEFAULT_SYSTEM,
    variables = variables(input.target),
    kwargs...,
)
    vars = variables
    F, G = input.target, input.start
    F_ishom, G_ishom = is_homogeneous.((F, G))
    if !affine_tracking && F_ishom && G_ishom
        f, regauging_factors = apply_system_scaling(F, vars, system_scaling)
        vargroups = VariableGroups(vars, homvar)
        F̄ = construct_system(f, system; variables = vars, homvars = homvar)
        Ḡ = construct_system(G, system; variables = vars, homvars = homvar)
        prob = Problem{ProjectiveTracking}(
            Ḡ,
            F̄,
            vargroups,
            seed;
            regauging_factors = regauging_factors,
            startsolutions_need_reordering = true,
            kwargs...,
        )
        prob, startsolutions
    elseif F_ishom || G_ishom
        throw(ArgumentError("One of the input polynomials is homogeneous and the other not. Make either both homogeneous or non-homogeneous!"))
    elseif !affine_tracking
        homvar === nothing || throw(ArgumentError("Input system is not homogeneous although `homvar` was passed."))

        h = uniquevar(F)
        push!(vars, h)
        sort!(vars, rev = true)
        g, f = homogenize(G, h), homogenize(F, h)
        f, regauging_factors = apply_system_scaling(f, vars, system_scaling)
        F̄ = construct_system(f, system; variables = vars, homvars = homvar)
        Ḡ = construct_system(g, system; variables = vars, homvars = homvar)
        vargroups = VariableGroups(vars, h)
        prob = Problem{ProjectiveTracking}(
            Ḡ,
            F̄,
            vargroups,
            seed;
            regauging_factors = regauging_factors,
            startsolutions_need_reordering = true,
            kwargs...,
        )
        prob, startsolutions
    else # affine_tracking
        f, regauging_factors = apply_system_scaling(F, vars, system_scaling)
        F̂ = construct_system(f, system; variables = vars)
        Ĝ = construct_system(G, system; variables = vars)
        prob = Problem{AffineTracking}(
            Ĝ,
            F̂,
            VariableGroups(vars),
            seed;
            regauging_factors = regauging_factors,
            kwargs...,
        )
        prob, startsolutions
    end
end

#####################
# Parameter homotopy
#####################

function problem_startsolutions(
    input::ParameterSystemInput{<:MPPolyInputs},
    startsolutions,
    hominfo,
    seed;
    affine_tracking = default_affine_tracking(input, hominfo),
    system = SPSystem,
    system_scaling = nothing,
    variables = variables(input.system; parameters = input.parameters),
    kwargs...,
)
    Prob = affine_tracking ? Problem{AffineTracking} : Problem{ProjectiveTracking}
    if affine_tracking
        variable_groups = VariableGroups(variables)
        homvars = nothing
        F = input.system
    else
        F, variable_groups, homvars = homogenize_if_necessary(
            input.system,
            hominfo;
            vars = variables,
            parameters = input.parameters,
        )
    end
    vars = flattened_variable_groups(variable_groups)
    F̂ = construct_system(
        F,
        system;
        homvars = homvars,
        variables = vars,
        parameters = input.parameters,
    )
    H = ParameterHomotopy(F̂, p₁ = input.p₁, p₀ = input.p₀, γ₁ = input.γ₁, γ₀ = input.γ₀)
    Prob(H, variable_groups, seed; startsolutions_need_reordering = !affine_tracking),
    startsolutions
end

function problem_startsolutions(
    input::ParameterSystemInput{<:AbstractSystem},
    startsolutions,
    hominfo,
    seed;
    affine_tracking = default_affine_tracking(input, hominfo),
    system = nothing,
    system_scaling = nothing,
    kwargs...,
)
    n, N = size(input.system)
    H = ParameterHomotopy(
        input.system,
        p₁ = input.p₁,
        p₀ = input.p₀,
        γ₁ = input.γ₁,
        γ₀ = input.γ₀,
    )
    variable_groups = VariableGroups(N, hominfo)
    tracking_type = affine_tracking ? AffineTracking() : ProjectiveTracking()
    Problem(
        tracking_type,
        H,
        variable_groups,
        seed;
        startsolutions_need_reordering = false,
    ),
    startsolutions
end


text: ## experiments with variational Bayes

ext: jl


text:
import Statistics: mean
using Printf
using DelimitedFiles
import Images
import ImageFiltering
using Base.Threads
using MultivariateStats
using Statistics

"""
    swap_xy(lms)

Swaps x and y coordinates of a landmark array.
"""
function swap_xy(lms)
  new = deepcopy(lms)
  for i in 1:size(lms, 2)
    for cor in 1:3:size(lms,1)
      new[cor, i] = lms[cor+1, i]
      new[cor+1, i] = lms[cor, i]
    end
  end
  return new
end

"""
    empty_3rd_dimension(lms)

returns a new 2D array with 3 coordinates per landmark
with the third (z-) coordinate being 0 for all landmarks.
"""
function empty_3rd_dimension(lms)
    new_lms = zeros(Float32, Int64(size(lms,1)*3/2), size(lms, 2))
    counter = 0
    for i in 1:3:60
        new_lms[i:i+1,:] = lms[i-counter:i-(counter-1), :]
        counter+=1
    end
    return new_lms
end

"""
    change_values!(m, from, to, rel)

Quickly change elements that are `rel`(`from`) to `to`.
"""
function change_values!(m, from, to, rel)
  for index in 1:length(m)
    if rel(m[index], from)
      m[index] = to
    end
  end
end

"""
    reshape_4d_to_5d(tensor)

Takes a 4d tensor with 3d volumes and returns them in a 5d tensor
with one channel per volume.
"""
function reshape_4d_to_5d(tensor)
  img_data = reshape(tensor, (size(tensor)[1],size(tensor)[2],size(tensor)[3], 1, size(tensor)[4]))
  return img_data
end

"""
    resize_2d_images(imgs, lms, out_size)

Resizes 4D tensor of 2D images to size(outsize[1], outsize[2]) and
adjusts landmarks accordingly. out_size is a tuple with sizes in
dim1 and dim2.
"""
function resize_2d_images(imgs, lms, out_size)
    n_inds = size(imgs, 4)
    out = zeros(Float32, out_size[1], out_size[2], size(imgs, 3), n_inds)
    out_lms = deepcopy(lms)
    ratio_x = out_size[1]/size(imgs, 1)
    ratio_y = out_size[2]/size(imgs, 2)
    for i in 1:2:size(lms, 1)
        out_lms[i,:] .= lms[i,:] .* ratio_x
        out_lms[i+1, :] .= lms[i+1, :] .* ratio_y
    end
    for ind in 1:n_inds
        for img in 1:size(imgs, 3)
            out[:,:,img, ind] .= Images.imresize(imgs[:,:,img,ind], out_size)
        end
    end
    return out, out_lms
end

# need to build custom function, the scikit one doesn't work due to dimension problems
"""
    train_test_split_3d(X, y, train_ratio)

Randomly split features and labels with ratio 'train_ratio' and
return training and testing set. Takes a 5d tensor as X.
"""
function train_test_split_3d(X, y, train_ratio)
  indices_train = []
  indices_test = []
  for i in 1:size(X)[5]
    if rand() < train_ratio
      push!(indices_train, i)
    else
      push!(indices_test, i)
    end
  end
  x_train = X[:,:,:,:,indices_train]
  y_train = y[:,indices_train]
  x_test = X[:,:,:,:,indices_test]
  y_test = y[:,indices_test]

  return x_train, x_test, y_train, y_test
end

"""
    regular_train_test_split_3d(X, y)

split features and labels so that every fifth entry in the original
dataset is in the testin set (80/20 split) and
return training and testing set. Takes a 5d tensor as X.
"""
function regular_train_test_split_3d(X, y)
  indices_train = []
  indices_test = []
  for i in 1:size(y)[2]
    if i % 5 != 0
      push!(indices_train, i)
    else
      push!(indices_test, i)
    end
  end
  if length(size(X)) == 5
    x_train = X[:,:,:,:,indices_train]
    y_train = y[:,indices_train]
    x_test = X[:,:,:,:,indices_test]
    y_test = y[:,indices_test]
  else
    x_train = X[:,:,:,indices_train]
    y_train = y[:,indices_train]
    x_test = X[:,:,:,indices_test]
    y_test = y[:,indices_test]
  end
  return x_train, x_test, y_train, y_test
end

"""
    regular_train_test_split_3d(X, y)

split features and labels so that every fifth entry in the original
dataset is in the testin set (80/20 split) and
return training and testing set. Takes a 4d tensor as X.
"""
function regular_train_test_split_2d(X, y)
  indices_train = []
  indices_test = []
  for i in 1:size(X)[4]
    if i % 5 != 0
      push!(indices_train, i)
    else
      push!(indices_test, i)
    end
  end
  x_train = X[:,:,:,indices_train]
  y_train = y[:,indices_train]
  x_test = X[:,:,:,indices_test]
  y_test = y[:,indices_test]

  return x_train, x_test, y_train, y_test
end

"""
    landmark_to_surface(volumes, landmarks, radius)

Moves landmarks that are not already on the volume to the closest
point to them on the surface of the volume.
"""
function landmark_to_surface(volumes, landmarks, radius)
  inds = size(landmarks, 2)
  lm_out = deepcopy(landmarks)
  for ind in 1:inds
    image = volumes[:,:,:,ind]
    maxi = maximum(volumes[:,:,:,ind])
    landm = landmarks[:, ind]
    n_points = floor(Int, length(landm)/3)
    coordinates_lm = zeros(n_points, 3)
    for i in 1:n_points
        for j in 1:3
            coordinates_lm[i,j] = landm[3*i-(3-j)]
        end
    end
    lm_coord = floor.(Int, round.(coordinates_lm .*10))
    for point in 1:size(lm_coord, 1)
      if image[Base.max(1, lm_coord[point, 1]), Base.max(1, lm_coord[point, 2]), Base.max(1, lm_coord[point, 3])] != maxi
        continue
      end
      points_in_area = findall(x->x!=maxi, image[Base.max(1, lm_coord[point, 1]-radius):Base.min(size(image, 1), lm_coord[point, 1]+radius),
      Base.max(1, lm_coord[point, 2]-radius):Base.min(size(image, 2), lm_coord[point, 2]+radius),
      Base.max(1, lm_coord[point, 3]-radius):Base.min(size(image, 3), lm_coord[point, 3]+radius)])
      if length(points_in_area) == 0
        continue
      end
      nearest = [0, 10000000]
      for p in 1:length(points_in_area)
        dist = sum([(radius-points_in_area[p][1])^2, (radius-points_in_area[p][2])^2,
        (radius-points_in_area[p][3])^2])
        if dist < nearest[2]
          nearest[2] = dist
          nearest[1] = p
        end
      end
      lm_out[point*3-2, ind] = (points_in_area[nearest[1]][1] + lm_coord[point, 1]-radius)/10
      lm_out[point*3-1, ind] = (points_in_area[nearest[1]][2] + lm_coord[point, 2]-radius)/10
      lm_out[point*3, ind] = (points_in_area[nearest[1]][3] + lm_coord[point, 3]-radius)/10
    end
  end
  return lm_out
end

"""
    image_gradients(x)

Adds the sum of the image gradients in x and y direction for each channel
of each image in a 4-dimensional tensor as a new channel.
"""
function image_gradients(x)
  copyx = deepcopy(x)
  for ind in 1:size(x, 4)
    for chan in 1:size(x, 3)
      gradix, gradiy = ImageFiltering.imgradients(x[:,:,chan,ind], ImageFiltering.Kernel.ando3)
      copyx[:,:,chan,ind] = abs.(gradix) .+ abs.(gradiy)
    end
  end
  out = cat(x, copyx, dims=3)
  return out
end

"""
    align_principal(volumes, landmarks, output_size)

Aligns all volumes in tensor "volumes" along their principal
axis (new z-axis) and transforms respective landmark data alongside.
The output size can be specified with output_size
(a cube with size output_size^3 will be returned), but must be
large enough to contain the largest aligned volume. a median filter
will be applied to handle caused by transformation. Also returns
a reconstruction array used to translate predicted landmarks back
to their original volume.
"""
function align_principal(volumes, landmarks, output_size::Int)
  inds = size(landmarks, 2)
  out = []
  lm_out = deepcopy(landmarks)
  reconstruction = zeros(3, 5, inds)
  for ind in 1:inds
    image = volumes[:,:,:,ind]
    maxi = maximum(volumes[:,:,:,ind])
    points = findall(x->x!=maxi, image)
    coordinates = zeros(length(points), 3)
    for i in 1:length(points)
        for j in 1:3
            coordinates[i,j] = points[i][j]
        end
    end
    pcs = fit(PCA, coordinates', maxoutdim=5)
    tra = transform(pcs, coordinates')
    positive = tra .- (minimum(tra)-1)
    maxi = round(maximum(positive))
    pos = floor.(Int, round.(positive))
    vol2 = zeros(output_size, output_size, output_size, 1)
    for pnt in 1:size(pos,2)
        vol2[pos[1, pnt], pos[2, pnt], pos[3, pnt], 1] = 1.0
    end
    vol3 = deepcopy(vol2)
    # for i in 1:output_size
    #   vol3[:,:,i,1] = ImageFiltering.mapwindow(median!, vol2[:,:,i,1], (3,3))
    # end
    change_values!(vol3, 0.0, 2.0, ==)
    if ind == 1
      out = vol3
    else
      out = cat(out, vol3, dims = 4)
    end
    landm = landmarks[:, ind]
    n_points = floor(Int, length(landm)/3)
    coordinates_lm = zeros(n_points, 3)
    for i in 1:n_points
        for j in 1:3
            coordinates_lm[i,j] = landm[3*i-(3-j)]
        end
    end
    lm_coord = coordinates_lm .*10
    new_lm = transform(pcs, lm_coord')
    pos_lm = new_lm .- (minimum(tra)-1)
    for i in 1:n_points
        for j in 1:3
            lm_out[3*i-(3-j), ind] = pos_lm[j, i]
        end
    end
    reconstruction[:,1:3, ind] = projection(pcs)
    reconstruction[1,4, ind] = minimum(tra)
    reconstruction[2,4, ind] = 1
    reconstruction[:,5,ind] = mean(pcs)
  end
  return out, relu.(lm_out ./ 10), reconstruction
end

"""
    resize_relevant(vols, lms, out_size)

Takes the output from align_principal and resizes the relevant
part (part containing the actual item) to a cube of specified size.
Also returns the ratio between resized and original volume size.
Adjusts landmarks.
"""
function resize_relevant(vols, lms, out_size)
  inds = size(vols, 4)
  out = zeros(out_size, out_size, out_size, inds)
  scales = []
  lms_out = deepcopy(lms)
  for ind in 1:inds
    maxi = maximum(vols[:,:,:,ind])
    coords = findall(x->x!=maxi, vols[:,:,:,ind])
    max_coord = maximum(coords)
    ar = zeros(1,3)
    for i in 1:3
      ar[i] = max_coord[i]
    end
    maximu = Int(maximum(ar))
    scale = out_size/maximu
    push!(scales, scale)
    out[:,:,:,ind] = Images.imresize(vols[1:maximu,1:maximu,1:maximu,ind], (out_size, out_size, out_size))
    lms_out[:, ind] = scale .* lms[:, ind]
  end
  return out, lms_out, scales
end

"""
    translate_lms_back(lms, reconstruction_array)

Takes an array of predicted landmarks and a reconstruction array
as returned by align_principal() and translates the landmarks
back to their respective original volume.
"""
function translate_lms_back(lms, reconstruction_array)
  lm_out = deepcopy(lms)
  for ind in 1:size(lms, 2)
    landm = lms[:,ind]
    n_points = floor(Int, length(landm)/3)
    coordinates_lm = zeros(n_points, 3)
    for i in 1:n_points
        for j in 1:3
            coordinates_lm[i,j] = landm[3*i-(3-j)]
        end
    end
    coord = coordinates_lm .* 10
    coord_scaled = coord ./ reconstruction_array[2,4,ind]
    coord_center = coord_scaled .+ (reconstruction_array[1,4,ind] - 1)
    new_coord = (reconstruction_array[1:3,1:3,ind] * coord_center')

    for i in 1:n_points
        for j in 1:3
            lm_out[3*i-(3-j), ind] = new_coord[j, i] + reconstruction_array[j, 5, ind]
        end
    end
  end
  return lm_out ./ 10
end

"""
    choose_dims(y, dims)

choose which dimensions of x,y and z should be included in the labels and return
a new array containing only the specified dims.
"""
function choose_dims(y, dims)
  newy = []
  for i in 1:3:size(y)[1]
    for o in 1:3
      if dims[o] == 1
        if size(newy)[1] == 0
          newy = y[i+(o-1), :]
        else
          newy = hcat(newy, y[i+(o-1), :])
        end
      else
        continue
      end
    end
  end
  return newy'
end

"""
    give_z_value(x)

sets every value outside of the actual volume to 0 and changes the voxel-values
of every voxel within the volume to its corresponding z-coordinate. Returns
the new tensor with all values normalized between -1 and 1.
"""
function give_z_value(x)
  copyx = deepcopy(x)
  Threads.@threads for ind in 1:size(x)[4]
    siz = size(x)[3]
    for i in 1:siz
      img = x[:,:,i, ind]
      ma = maximum(img)
      img[findall(x->x==ma, img)] .= -1
      img[findall(x->x!=-1, img)] .= (2*((i-1)/(siz-1))-1)*-1
      copyx[:,:,i,ind] = img
    end
  end
  return copyx
end

"""
    make_elevation_map(x)

Returns a 2d image for every volume in tensor x and sets the pixel values
to the maximum value within each column of voxels in dimension z as returned
by the function give_z_value().
"""
function depth_map(x)
  max_val = maximum(x)
  x_out = zeros(Float32, size(x)[1], size(x)[2], 1, size(x)[4])
  Threads.@threads for inds in 1:size(x)[4]
    for dim1 in 1:size(x)[1]
      for dim2 in 1:size(x)[2]
        x_out[dim1, dim2, 1, inds] = maximum(x[dim1,dim2,:,inds])/max_val
      end
    end
  end
  return x_out
end

"""
    depth_map_all_sides(x)

returns a 4d image tensor with 6 channels per image, each representing
a depthmap from one side of the volume.
"""
function depth_map_all_sides(x)
  mock_y = zeros(3, size(x, 4))
  first = give_z_value(x)
  second = -1 .* first
  change_values!(second, 1.0, -1.0, ==)
  out = depth_map(first)
  out = cat(out, depth_map(second), dims=3)
  first = nothing
  second = nothing
  front_back, y_ = flip_volume_front(x, mock_y)
  third = give_z_value(front_back)
  forth = -1 .* third
  change_values!(forth, 1.0, -1.0, ==)
  out = cat(out, depth_map(third), dims=3)
  out = cat(out, depth_map(forth), dims=3)
  third = nothing
  forth = nothing
  side_side, y_ = flip_volume_side(front_back, mock_y)
  front_back = nothing
  fifth = give_z_value(side_side)
  sixth = -1 .* fifth
  change_values!(sixth, 1.0, -1.0, ==)
  out = cat(out, depth_map(fifth), dims=3)
  out = cat(out, depth_map(sixth), dims=3)
  return out
end

max_stars_repo_name: imi-diener/AutoLM

max_issues_repo_name: imi-diener/AutoLM

max_forks_repo_name: imi-diener/AutoLM


text:
"""
    @rec(rng, eqn)

Computes recursive calculations for the given range `rng` and equation `eqn`.

### Examples
```julia-repl
julia> s = TSeries(1U, zeros(1))
1-element Unit TSeries from 1U:
      1U : 0.0

julia> s[1U] = 0
0

julia> s[2U] = 1
1

julia> s
2-element Unit TSeries from 1U:
      1U : 0.0
      2U : 1.0

julia> @rec(3U:10U, s[t] = s[t-1] + s[t-2])

julia> s
10-element Unit TSeries from 1U:
      1U : 0.0
      2U : 1.0
      3U : 1.0
      4U : 2.0
      5U : 3.0
      6U : 5.0
      7U : 8.0
      8U : 13.0
      9U : 21.0
     10U : 34.0
```
"""
macro rec(rng, eqn)
    eqn isa Expr && eqn.head == :(=) || error("Expression must be an assignment.")
    lhs, rhs = eqn.args
    lhs isa Expr && lhs.head == :ref || error("Left hand side of assignment must be an indexing expression.")
    vn, inds = lhs.args[1], lhs.args[2:end]
    tn = let # try to find the indexing variable
        tn = Symbol[]
        for i in 1:length(inds)
            isa(inds[i], Symbol) || continue
            inds[i] == :(:) && continue
            if inds[i] == :t # prefer :t if it's in there
                tn = [:t]
                break
            else
                push!(tn, inds[i])
            end
        end
        isempty(tn) && error("Time index not found.")
        # :t not there and multiple other symbols.
        length(tn) > 1 && error("Ambiguous time index. Use `t`.")
        tn[1]
    end
    return quote
        for $(tn) in $(rng)
            $(eqn)
        end
    end |> esc
end

export @rec


text:
#=
Model: PDE (84 variables, 1 input)
=#
using Reachability, MAT

pde(o::Pair{Symbol, <:Any}...) = pde(Options(Dict{Symbol, Any}(o)))

function pde(input_options::Options)
    # =====================
    # Problem specification
    # =====================
    file = matopen(@relpath "pde.mat")
    A = float(read(file, "A")) # this matrix has Int entries

    # initial set
    n = size(A, 1)
    center0 = zeros(n)
    radius0 = zeros(n)
    center0[65:80] .= 0.00125
    center0[81:84] .= -0.00175
    radius0[65:84] .= 0.00025
    X0 = Hyperrectangle(center0, radius0)

    # input set
    B = read(file, "B")
    U = BallInf([0.75], .25)

    # instantiate continuous LTI system
    S = InitialValueProblem(
        ConstrainedLinearControlContinuousSystem(A, B, nothing, U), X0)

    # property: y < 12 for linear combination y
    property = SafeStatesProperty(HalfSpace(
        read(matopen(@relpath "out.mat"), "M")[1, :], 12.))

    # =======================
    # Problem default options
    # =======================
    partition = [(2*i-1:2*i) for i in 1:42] # 2D blocks

    if input_options[:mode] == "reach"
        problem_options = Options(:vars => [1],
                                  :partition => partition,
                                  :plot_vars => [0, 1])
        # :projection_matrix => sparse(read(matopen(@relpath "out.mat"), "M"))
    elseif input_options[:mode] == "check"
        problem_options = Options(:vars => 1:84,
                                  :partition => partition,
                                  :property => property)
    end

    return (S, merge(problem_options, input_options))
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


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

# Compute y += a*x using @simd for vectors x and y
function simd_axpy( a, x, y )
    # LLVM's auto-vectorizer typically vectorizes this loop even without @simd
    @simd for i=1:length(x)
        @inbounds y[i] += a*x[i]
    end
end

# Run axpy(a,x,y) m times
function flog_axpy( m, a, x, y )
    for j=1:m
        simd_axpy(a,x,y)
    end
end

# Run axpy for Float32 and Float64
for t in [Float32,Float64]
    n = 1000
    x = rand(t,n)
    y = rand(t,n)
    a = convert(t,0.5)
    bits = 8*sizeof(t)
    @timeit(flog_axpy(100,a,x,y), "simd_axpy_$bits", "SIMD BLAS axpy for type $t", "SIMD")
end


text:
# Autogenerated wrapper script for AprilTags_jll for i686-w64-mingw32
export libapriltag

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using OrdinaryDiffEq
using DynamicalSystemsBase
using ChaosTools

ds = Systems.lorenz()

tim = time()

lyapunovs(ds, 2000)

tim = time() - tim

if typeof(ds) <: DiscreteDynamicalSystem
    println("Discrete system. ")
else
    println("Default solver: ", DynamicalSystemsBase.CDS_KWARGS.alg)
end
println("Time to first run of `lyapunov` ( ≡ compile time):")
println(tim, " seconds.\n")

# Results:
# Default solver: SimpleDiffEq.SimpleATsit5()
# Time to first run of `lyapunov` ( ≡ compile time):
# 7.906000137329102 seconds.
#
# Default solver: Tsit5()
# Time to first run of `lyapunov` ( ≡ compile time):
# 13.582000017166138 seconds.
#
# Default solver: Vern9(true)
# Time to first run of `lyapunov` ( ≡ compile time):
# 24.50499987602234 seconds.


text:
module EllipticEquationMod

using PDESolver
using ArrayViews
import ArrayViews.view
using SolverCommon
using ODLCommonTools
using SummationByParts
using PdePumiInterface
#using NonlinearSolvers
using Jacobian
using LinearSolvers
using Input
# using ForwardDiff
using Utils
using MPI
using Input  # input file processing
using PETSc2

import ODLCommonTools.sview

export AbstractEllipticData, EllipticData, EllipticData_, run_elliptic

abstract type AbstractEllipticData{Tsol, Tres} <: AbstractSolutionData{Tsol, Tres} end
abstract type EllipticData{Tsol, Tres, Tdim} <: AbstractEllipticData{Tsol, Tres} end

# now that EllipticData is declared, include other files that use it
include(joinpath(Pkg.dir("PDESolver"), "src/solver/debug.jl"))  # debug macro
include("elliptic.jl")
# include("time_advance.jl")
include("bc.jl")
include("ic.jl")
include("flux.jl")
include("check_options.jl")
include("source.jl")
#include("output.jl")
include("exactSolution.jl")
include("diffusion.jl")
include("functional.jl")
include("types.jl")
include("startup_func.jl")  # function for invoking the solver

global const PhysicsName = "Elliptic"
register_physics(PhysicsName, EllipticEquationMod, createObjects, checkOptions)
end # end of module


text:
using Documenter, MyExample

makedocs(;
    modules=[MyExample],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/Prabu Annamalai/MyExample.jl/blob/{commit}{path}#L{line}",
    sitename="MyExample.jl",
    authors="Prabu Annamalai",
    assets=String[],
)

deploydocs(;
    repo="github.com/Prabu Annamalai/MyExample.jl",
)


text:
# Some helper functions
function get_cellset(cell_type, cells)
    return Set(findall(c-> typeof(c) == cell_type, cells))
end

function create_field(;name, field_dim, order, spatial_dim, cellshape)
    interpolation = Lagrange{spatial_dim, cellshape, order}()
    return Field(name, interpolation, field_dim)
end

function get_2d_grid()
    # GIVEN: two cells, a quad and a triangle sharing one face
    cells = [
        Quadrilateral((1, 2, 3, 4)),
        Triangle((3, 2, 5))
        ]
    coords = zeros(Vec{2,Float64}, 5)
    nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 5)]
    return Grid(cells,nodes)
end

function test_1d_bar_beam()
    # Something like a truss and a Timoshenko beam.
    # two line-cells with 2 dofs/node -> "bars"
    # one line-cell with a mixed field: one vector field and one scalar -> "beam"
    cells = [
        Line((1, 2)),
        Line((2, 3)),
        Line((1, 3)),
        ]
    nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 3)]
    grid = Grid(cells, nodes)

    field1 = create_field(name=:u, spatial_dim=1, field_dim=2, order=1, cellshape=RefCube)
    field2 = create_field(name=:u, spatial_dim=1, field_dim=2, order=1, cellshape=RefCube)
    field3 = create_field(name=:θ, spatial_dim=1, field_dim=1, order=1, cellshape=RefCube)

    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field2, field3], Set(3)));
    push!(dh, FieldHandler([field1], Set((1, 2))));
    close!(dh)
    @test ndofs(dh) == 8
    @test celldofs(dh, 3) == collect(1:6)
    @test celldofs(dh, 1) == [1, 2, 7, 8]
    @test celldofs(dh, 2) == [7, 8, 3, 4]

    #        7,8
    #       /   \
    #     /      \
    #   /_________\
    # 1,2,5       3,4,6
end

function test_2d_scalar()

    grid = get_2d_grid()
    # WHEN: adding a scalar field for each cell and generating dofs
    field1 = create_field(name=:u, spatial_dim=2, field_dim=1, order=1, cellshape=RefCube)
    field2 = create_field(name=:u, spatial_dim=2, field_dim=1, order=1, cellshape=RefTetrahedron)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1], Set(1)));
    push!(dh, FieldHandler([field2], Set(2)));
    close!(dh)

    # THEN: we expect 5 dofs and dof 2 and 3 being shared
    @test ndofs(dh) == 5
    @test dh.cell_dofs.values == [1, 2, 3, 4, 3, 2, 5]
    @test celldofs(dh, 1) == [1, 2, 3, 4]
    @test celldofs(dh, 2) == [3, 2, 5]
end

function test_2d_error()

    grid = get_2d_grid()
    # the refshape of the field must be the same as the refshape of the elements it is added to
    field1 = create_field(name=:u, spatial_dim=2, field_dim=1, order=1, cellshape=RefTetrahedron)
    field2 = create_field(name=:u, spatial_dim=2, field_dim=1, order=1, cellshape=RefCube)
    dh = MixedDofHandler(grid);
    @test_throws ErrorException push!(dh, FieldHandler([field1], Set(1)));
    @test_throws ErrorException push!(dh, FieldHandler([field2], Set(2)));
    # all cells within a FieldHandler should be of the same celltype
    @test_throws ErrorException push!(dh, FieldHandler([field1], Set([1,2])));

end

function test_2d_vector()

    grid = get_2d_grid()
    ## vector field
    field1 = create_field(name = :u, spatial_dim=2, field_dim = 2, order = 1, cellshape = RefCube)
    field2 = create_field(name = :u, spatial_dim=2, field_dim = 2, order = 1, cellshape = RefTetrahedron)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1], Set(1)));
    push!(dh, FieldHandler([field2], Set(2)));
    close!(dh)

    # THEN: we expect 10 dofs and dof 3-6 being shared
    @test ndofs(dh) == 10
    @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8]
    @test celldofs(dh, 2) == [5, 6, 3, 4, 9, 10]

end

function test_2d_mixed_1_el()
    grid = get_2d_grid()
    ## mixed field of same order
    field1 = create_field(name = :u, spatial_dim=2, field_dim = 2, order = 1, cellshape = RefCube)
    field2 = create_field(name = :p, spatial_dim=2, field_dim = 1, order = 1, cellshape = RefCube)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1, field2], Set(1)));
    close!(dh)

    # THEN: we expect 12 dofs
    @test ndofs(dh) == 12
    @test ndofs_per_cell(dh, 1) == 12
    @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
end

function test_2d_mixed_2_el()

    grid = get_2d_grid()
    ## mixed field of same order
    field1_quad = create_field(name = :u, spatial_dim=2, field_dim = 2, order = 1, cellshape = RefCube)
    field2_quad = create_field(name = :p, spatial_dim=2, field_dim = 1, order = 1, cellshape = RefCube)
    field1_tri = create_field(name = :u, spatial_dim=2, field_dim = 2, order = 1, cellshape = RefTetrahedron)
    field2_tri = create_field(name = :p, spatial_dim=2, field_dim = 1, order = 1, cellshape = RefTetrahedron)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1_quad, field2_quad], Set(1)));
    push!(dh, FieldHandler([field1_tri, field2_tri], Set(2)));
    close!(dh)

    # THEN: we expect 15 dofs
    @test ndofs(dh) == 15
    @test ndofs_per_cell(dh, 1) == 12
    @test ndofs_per_cell(dh, 2) == 9
    @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
    @test celldofs(dh, 2) == [5, 6, 3, 4, 13, 14, 11, 10, 15]
end

function test_face_dofs_2_tri()
    cells = [
        Triangle((1, 2, 3)),
        Triangle((2, 4, 3))
        ]
    nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 4)]
    grid = Grid(cells, nodes);
    field1 = create_field(name = :u, spatial_dim = 2, field_dim = 2, order = 2, cellshape = RefTetrahedron)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1], Set((1, 2))));
    #push!(dh, FieldHandler([field2], Set(2)));
    close!(dh)

    # THEN:
    @test ndofs(dh) == 18
    @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
    @test celldofs(dh, 2) == [3, 4, 13, 14, 5, 6, 15, 16, 17, 18, 9, 10]
end

function test_3d_tetrahedrons()
    cells = [
        Tetrahedron((1, 2, 3, 7)),
        Tetrahedron((1, 5, 2, 7)),
        Tetrahedron((2, 4, 3, 7)),
        Tetrahedron((2, 8, 4, 7)),
        Tetrahedron((2, 5, 6, 7)),
        Tetrahedron((2, 6, 8, 7)),
        ]
    nodes = [Node(coord) for coord in zeros(Vec{3,Float64}, 8)]
    grid = Grid(cells, nodes)
    field = create_field(name = :u, spatial_dim=3,  field_dim = 3, order = 2, cellshape = RefTetrahedron)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field], Set((1, 2, 3, 4, 5, 6))));
    close!(dh)

    # reference using the regular DofHandler
    tet_grid = generate_grid(Tetrahedron, (1, 1,1))
    tet_dh = DofHandler(tet_grid)
    push!(tet_dh, :u, 3, Lagrange{3,RefTetrahedron,2}())
    close!(tet_dh)

    for i in 1:6
        @test celldofs(dh, i) == celldofs(tet_dh, i)
    end
end

function test_face_dofs_quad_tri()
    # quadratic quad and quadratic triangle
    grid = get_2d_grid()
    field1 = create_field(name = :u, spatial_dim = 2, field_dim = 2, order = 2, cellshape = RefCube)
    field2 = create_field(name = :u, spatial_dim = 2, field_dim = 2, order = 2, cellshape = RefTetrahedron)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1], Set(1)));
    push!(dh, FieldHandler([field2], Set(2)));
    close!(dh)

    # THEN:
    @test ndofs(dh) == 24
    @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]
    @test celldofs(dh, 2) == [5, 6, 3, 4, 19, 20, 11, 12, 21, 22, 23, 24]
end

function test_serendipity_quad_tri()
    # bi-quadratic quad (Serendipity) and quadratic triangle
    grid = get_2d_grid()
    interpolation = Serendipity{2, RefCube, 2}()
    field1 = Field(:u, interpolation, 2)
    field2 = create_field(name = :u, spatial_dim = 2, field_dim = 2, order = 2, cellshape = RefTetrahedron)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1], Set(1)));
    push!(dh, FieldHandler([field2], Set(2)));
    close!(dh)

    # THEN:
    @test ndofs(dh) == 22
    @test celldofs(dh, 1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
    @test celldofs(dh, 2) == [5, 6, 3, 4, 17, 18, 11, 12, 19, 20, 21, 22]
end

function test_2d_mixed_field_triangles()
    # Mixed field (Taylor-Hood)
    # quadratic vector field + linear scalar field on
    # celltypes: 2 Triangles
    cells = [
        Triangle((1, 2, 3)),
        Triangle((2, 4, 3))
        ]
    nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 4)]
    grid = Grid(cells, nodes)
    field1 = create_field(name=:u, spatial_dim=2, field_dim=2, order=2, cellshape=RefTetrahedron)
    field2 = create_field(name=:p, spatial_dim=2, field_dim=1, order=1, cellshape=RefTetrahedron)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1, field2], Set((1, 2))));
    close!(dh)
    @test ndofs(dh) == 22
    @test celldofs(dh, 1) == collect(1:15)
    @test celldofs(dh, 2) == [3, 4, 16, 17, 5, 6, 18, 19, 20, 21, 9, 10, 14, 22, 15]
end

function test_2d_mixed_field_mixed_celltypes()
    # Mixed field (Taylor-Hood)
    # quadratic vector field + linear scalar field on
    # celltypes: 1 Quadrilateral and 1 Triangle

    grid = get_2d_grid()
    field1 = create_field(name=:u, spatial_dim=2, field_dim=2, order=2, cellshape=RefCube)
    field2 = create_field(name=:p, spatial_dim=2, field_dim=1, order=1, cellshape=RefCube)
    field3 = create_field(name=:u, spatial_dim=2, field_dim=2, order=2, cellshape=RefTetrahedron)
    field4 = create_field(name=:p, spatial_dim=2, field_dim=1, order=1, cellshape=RefTetrahedron)
    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1, field2], Set(1)));
    push!(dh, FieldHandler([field3, field4], Set(2)));
    close!(dh)
    @test ndofs(dh) == 29
    @test celldofs(dh, 1) == collect(1:22)
    @test celldofs(dh, 2) == [5, 6, 3, 4, 23, 24, 11, 12, 25, 26, 27, 28, 21, 20, 29]
end

function test_3d_mixed_field_mixed_celltypes()
    # Mixed celltypes in 3d
    # one hex-cell with a vector field (3 dofs/node)
    # one quad-cell (connected along bottom edge of the hex) with two mixed vector fields (3+3 dofs/node)
    # One field type is common for both cells
    cells = [
        Hexahedron((1, 2, 3, 4, 5, 6, 7, 8)),
        Quadrilateral((3, 2, 9, 10)),
        ]
    nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 10)]
    grid = Grid(cells, nodes)

    # E.g. 3d continuum el -> 3dofs/node
    field1 = create_field(name=:u, spatial_dim=3, field_dim=3, order=1, cellshape=RefCube)
    # E.g. displacement field + rotation field -> 6 dofs/node
    field2 = create_field(name=:u, spatial_dim=3, field_dim=3, order=1, cellshape=RefCube)
    field3 = create_field(name=:θ, spatial_dim=3, field_dim=3, order=1, cellshape=RefCube)

    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([field1], Set(1)));
    push!(dh, FieldHandler([field2, field3], Set(2)));
    close!(dh)
    @test ndofs(dh) == 42
    @test celldofs(dh, 1) == collect(1:24)
    @test celldofs(dh, 2) == [7, 8, 9, 4, 5, 6, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42]
end

function test_2_element_heat_eq()
    # Regression test solving the heat equation.
    # The grid consists of two quad elements treated as two different cell types to test the Grid and all other necessary functions

    temp_grid = generate_grid(Quadrilateral, (2, 1))
    grid = Grid(temp_grid.cells, temp_grid.nodes)  # regular grid -> mixed grid
    grid.facesets = temp_grid.facesets;

    # Create two identical fields
    f1 = create_field(name=:u, spatial_dim=2, field_dim=1, order=1, cellshape=RefCube)
    f2 = create_field(name=:u, spatial_dim=2, field_dim=1, order=1, cellshape=RefCube)

    dh = MixedDofHandler(grid);
    push!(dh, FieldHandler([f1], Set(1)));  # first field applies to cell 1
    push!(dh, FieldHandler([f2], Set(2)));  # second field applies to cell 2
    close!(dh)

    # Create two Dirichlet boundary conditions - one for each field.
    ch = ConstraintHandler(dh);
    ∂Ω1 = getfaceset(grid, "left")
    ∂Ω2 = getfaceset(grid, "right")
    dbc1 = Dirichlet(:u, ∂Ω1, (x, t) -> 0)
    dbc2 = Dirichlet(:u, ∂Ω2, (x, t) -> 0)
    add!(ch, dh.fieldhandlers[1], dbc1);
    add!(ch, dh.fieldhandlers[2], dbc2);
    close!(ch)

    function doassemble(cellset, cellvalues, assembler, dh)

        n = ndofs_per_cell(dh, first(cellset))
        n_basefuncs = getnbasefunctions(cellvalues)
        Ke = zeros(n, n)
        fe = zeros(n)
        eldofs = zeros(Int, n)

        for cellnum in cellset
            celldofs!(eldofs, dh, cellnum)
            xe = getcoordinates(dh.grid, cellnum)
            reinit!(cellvalues, xe)
            fill!(Ke, 0)
            fill!(fe, 0)

            for q_point in 1:getnquadpoints(cellvalues)
                dΩ = getdetJdV(cellvalues, q_point)

                for i in 1:n_basefuncs
                    v  = shape_value(cellvalues, q_point, i)
                    ∇v = shape_gradient(cellvalues, q_point, i)
                    fe[i] += v * dΩ
                    for j in 1:n_basefuncs
                        ∇u = shape_gradient(cellvalues, q_point, j)
                        Ke[i, j] += (∇v ⋅ ∇u) * dΩ
                    end
                end
            end
            assemble!(assembler, eldofs, fe, Ke)
        end
    end

    K = create_sparsity_pattern(dh)
    f = zeros(ndofs(dh));
    assembler = start_assemble(K, f);
    # Use the same assemble function since it is the same weak form for both cell-types
    for fh in dh.fieldhandlers
        qr = QuadratureRule{2, RefCube}(2)
        interp = fh.fields[1].interpolation
        cellvalues = CellScalarValues(qr, interp)
        doassemble(fh.cellset, cellvalues, assembler, dh)
    end

    update!(ch, 0.0);
    apply!(K, f, ch)
    u = K \ f;

    # tested against heat_equation.jl (in the examples folder) using 2x1 cells and no
    # dbc on top and bottom boundary
    @test u == [0.0, 0.5, 0.5, 0.0, 0.0, 0.0]

    gridfilename = "mixed_grid"
    addcellset!(grid, "cell-1", [1,])
    addcellset!(grid, "cell-2", [2,])
    vtk_grid(gridfilename, grid) do vtk
        vtk_cellset(vtk, grid, "cell-1")
        vtk_cellset(vtk, grid, "cell-2")
        vtk_point_data(vtk, dh, u)
        # vtk_point_data(vtk, ch)  #FIXME
    end
    sha = bytes2hex(open(SHA.sha1, gridfilename*".vtu"))
    @test sha == "e96732c000b0b385db7444f002461468b60b3b2c"

end

function test_element_order()
    # Check that one can have non-contigous ordering of cells in a Grid
    # Something like this:
    #        ______
    #      /|     |\
    # (1) / | (2) | \ (3)
    #    /__|_____|__\
    cells = [
        Triangle((1, 2, 3)),
        Quadrilateral((2, 4, 5, 3)),
        Triangle((4, 6, 5))
        ]
    nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 6)]
    grid = Grid(cells, nodes)
    field1_tri = create_field(name=:u, spatial_dim=3, field_dim=2, order=1, cellshape=RefTetrahedron)
    field1_quad = create_field(name=:u, spatial_dim=3, field_dim=2, order=1, cellshape=RefCube)

    dh = MixedDofHandler(grid);
    # Note the jump in cell numbers
    push!(dh, FieldHandler([field1_tri], Set((1, 3))));
    push!(dh, FieldHandler([field1_quad], Set(2)));
    # Dofs are first created for cell 1 and 3, thereafter cell 2
    close!(dh)

    @test ndofs(dh) == 12
    @test celldofs(dh, 1) == collect(1:6)
    @test celldofs(dh, 2) == [3, 4, 7, 8, 11, 12, 5, 6]
    @test celldofs(dh, 3) == [7,8, 9, 10, 11, 12]

end

function test_field_on_subdomain()
    grid = get_2d_grid() # cell 1: quad, cell2: triangle
    dh = MixedDofHandler(grid)

    # assume two fields: a scalar field :s and a vector field :v
    # :v lives on both cells, :s lives only on the triangle
    ip_tri = Lagrange{2,RefTetrahedron,1}()
    ip_quad = Lagrange{2,RefCube,1}()
    v_tri = Field(:v, ip_tri, 2)
    v_quad = Field(:v, ip_quad, 2)
    s = Field(:s, ip_tri, 1)

    push!(dh, FieldHandler([v_quad], Set((1,))))
    push!(dh, FieldHandler([v_tri, s], Set((2,))))

    close!(dh)

    # retrieve field dimensions
    @test Ferrite.getfielddim(dh, :v) == 2
    @test Ferrite.getfielddim(dh, :s) ==1

    # find field in FieldHandler
    @test Ferrite.find_field(dh.fieldhandlers[1], :v) == 1
    @test Ferrite.find_field(dh.fieldhandlers[2], :v) == 1
    @test Ferrite.find_field(dh.fieldhandlers[2], :s) == 2
    @test_throws ErrorException Ferrite.find_field(dh.fieldhandlers[1], :s)
end

function test_mixed_grid_show()
    grid = get_2d_grid()
    str = sprint(show, MIME("text/plain"), grid)
    @test occursin("2 Quadrilateral/Triangle cells", str)
end

# regression tests for https://github.com/KristofferC/JuAFEM.jl/issues/315
function test_subparametric_quad()
    #linear geometry
    grid = generate_grid(Quadrilateral, (1,1))
    ip      = Lagrange{2,RefCube,2}()

    field = Field(:u, ip, 2)
    fh = FieldHandler([field], Set(1:getncells(grid)))

    dh = MixedDofHandler(grid)
    push!(dh, fh)
    close!(dh)

    ch = ConstraintHandler(dh)
    dbc1 = Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 0.0, 2)
    add!(ch, dbc1)
    close!(ch)
    update!(ch, 1.0)
    @test getnbasefunctions(Ferrite.getfieldinterpolation(dh,1)) == 9 # algebraic nbasefunctions
    @test celldofs(dh, 1) == [i for i in 1:18]
end

function test_subparametric_triangle()
    #linear geometry
    grid = generate_grid(Triangle, (1,1))

    ip = Lagrange{2,RefTetrahedron,2}()

    field = Field(:u, ip, 2)
    fh = FieldHandler([field], Set(1:getncells(grid)))

    dh = MixedDofHandler(grid)
    push!(dh, fh)
    close!(dh)

    ch = ConstraintHandler(dh)
    dbc1 = Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 0.0, 2)
    add!(ch, dbc1)
    close!(ch)
    update!(ch, 1.0)
    @test getnbasefunctions(Ferrite.getfieldinterpolation(dh,1)) == 6 # algebraic nbasefunctions
    @test celldofs(dh, 1) == [i for i in 1:12]
end

function test_separate_fields_on_separate_domains()
    # 5_______6
    # |\      |
    # |   \   |
    # 3______\4
    # |       |
    # |       |
    # 1_______2
    # Given: a vector field :q defined on the quad and a scalar field :t defined on the triangles
    nodes = [Node((0.0, 0.0)),
            Node((1.0, 0.0)),
            Node((0.0, 1.0)),
            Node((1.0, 1.0)),
            Node((0.0, 2.0)),
            Node((1.0, 2.0))]
    cells = Ferrite.AbstractCell[Quadrilateral((1,2,4,3)),
            Triangle((3,4,5)),
            Triangle((4,6,5))]
    mesh = Grid(cells, nodes)
    addcellset!(mesh, "quads", Set{Int}((1,)))
    addcellset!(mesh, "tris", Set{Int}((2, 3)))

    dh = MixedDofHandler(mesh)
    ip_tri = Lagrange{2,RefTetrahedron,1}()
    ip_quad = Lagrange{2,RefCube,1}()
    field = Field(:q, ip_quad, 2) # vector field :q only on quad
    fh_quad = FieldHandler([field], getcellset(mesh, "quads"))
    push!(dh, fh_quad)

    field = Field(:t, ip_tri, 1) # scalar field :t only on tris
    fh_tri = FieldHandler([field], getcellset(mesh, "tris"))
    push!(dh, fh_tri)
    close!(dh)

    # Expect: 8 dofs for the quad and 4 new dofs for the triangles
    @test ndofs(dh) == 12
    @test celldofs(dh, 1) == [i for i in 1:8]
    @test celldofs(dh, 2) == [9, 10, 11]
    @test celldofs(dh, 3) == [10, 12, 11]

end

@testset "MixedDofHandler" begin

    test_1d_bar_beam();
    test_2d_scalar();
    test_2d_error();
    test_2d_vector();
    test_2d_mixed_1_el();
    test_2d_mixed_2_el();
    test_face_dofs_2_tri();
    test_face_dofs_quad_tri();
    test_serendipity_quad_tri();
    test_3d_tetrahedrons();
    test_2d_mixed_field_triangles();
    test_2d_mixed_field_mixed_celltypes();
    test_3d_mixed_field_mixed_celltypes();
    test_2_element_heat_eq();
    test_element_order();
    test_field_on_subdomain();
    test_mixed_grid_show();
    test_subparametric_quad();
    test_subparametric_triangle();
    test_separate_fields_on_separate_domains();
end


text:
# Test starting and different position
@test check_string(UInt32[0x110000, 0x1f596], 2) == (1,0x10,1,0,0,0,0)

# Test valid sequences
for (seq, res) in (
    ([0x0],                (1,0,0,0,0,0,0)),  # Nul byte, beginning of ASCII range
    ([0x7f],               (1,0,0,0,0,0,0)),  # End of ASCII range
    ([0xc2,0x80],          (1,2,0,0,0,1,0)),  # \u80, beginning of Latin1 range
    ([0xc3,0xbf],          (1,2,0,0,0,1,0)),  # \uff, end of Latin1 range
    ([0xc4,0x80],          (1,4,0,0,1,0,0)),  # \u100, beginning of non-Latin1 2-byte range
    ([0xdf,0xbf],          (1,4,0,0,1,0,0)),  # \u7ff, end of non-Latin1 2-byte range
    ([0xe0,0xa0,0x80],     (1,8,0,1,0,0,0)),  # \u800, beginning of 3-byte range
    ([0xed,0x9f,0xbf],     (1,8,0,1,0,0,0)),  # \ud7ff, end of first part of 3-byte range
    ([0xee,0x80,0x80],     (1,8,0,1,0,0,0)),  # \ue000, beginning of second part of 3-byte range
    ([0xef,0xbf,0xbf],     (1,8,0,1,0,0,0)),  # \uffff, end of 3-byte range
    ([0xf0,0x90,0x80,0x80],(1,16,1,0,0,0,0)), # \U10000, beginning of 4-byte range
    ([0xf4,0x8f,0xbf,0xbf],(1,16,1,0,0,0,0)), # \U10ffff, end of 4-byte range
    ([0x0000],             (1,0,0,0,0,0,0)),  # Nul byte, beginning of ASCII range
    ([0x007f],             (1,0,0,0,0,0,0)),  # End of ASCII range
    ([0x0080],             (1,2,0,0,0,1,0)),  # Beginning of Latin1 range
    ([0x00ff],             (1,2,0,0,0,1,0)),  # End of Latin1 range
    ([0x0100],             (1,4,0,0,1,0,0)),  # Beginning of non-Latin1 2-byte range
    ([0x07ff],             (1,4,0,0,1,0,0)),  # End of non-Latin1 2-byte range
    ([0x0800],             (1,8,0,1,0,0,0)),  # Beginning of 3-byte range
    ([0xd7ff],             (1,8,0,1,0,0,0)),  # End of first part of 3-byte range
    ([0xe000],             (1,8,0,1,0,0,0)),  # Beginning of second part of 3-byte range
    ([0xffff],             (1,8,0,1,0,0,0)),  # End of 3-byte range
    ([0xd800,0xdc00],      (1,16,1,0,0,0,0)), # \U10000, beginning of 4-byte range
    ([0xdbff,0xdfff],      (1,16,1,0,0,0,0)), # \U10ffff, end of 4-byte range
    ([0x00000],            (1,0,0,0,0,0,0)),  # Nul byte, beginning of ASCII range
    ([0x0007f],            (1,0,0,0,0,0,0)),  # End of ASCII range
    ([0x00080],            (1,2,0,0,0,1,0)),  # Beginning of Latin1 range
    ([0x000ff],            (1,2,0,0,0,1,0)),  # End of Latin1 range
    ([0x00100],            (1,4,0,0,1,0,0)),  # Beginning of non-Latin1 2-byte range
    ([0x007ff],            (1,4,0,0,1,0,0)),  # End of non-Latin1 2-byte range
    ([0x00800],            (1,8,0,1,0,0,0)),  # Beginning of 3-byte range
    ([0x0d7ff],            (1,8,0,1,0,0,0)),  # End of first part of 3-byte range
    ([0x0e000],            (1,8,0,1,0,0,0)),  # Beginning of second part of 3-byte range
    ([0x0ffff],            (1,8,0,1,0,0,0)),  # End of 3-byte range
    ([0x10000],            (1,16,1,0,0,0,0)), # \U10000, beginning of 4-byte range
    ([0x10ffff],           (1,16,1,0,0,0,0))) # \U10ffff, end of 4-byte range
    @test check_string(seq) == res
end


text:
function interpolate_parameters_gdp(p, v, d, t)
    # interpolation of parameters, see notes in run_timestep for why these parameters
    if is_first(t)
        for annual_year = 2015:(gettime(t))
            yr = annual_year - 2015 + 1
            for r in d.region

                v.pop_population_ann[yr, r] = p.pop_population[t, r]
                v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r]

            end
        end
    else
        for annual_year = (gettime(t - 1) + 1):(gettime(t))
            yr = annual_year - 2015 + 1
            frac = annual_year - gettime(t - 1)
            fraction_timestep = frac / ((gettime(t)) - (gettime(t - 1)))

            for r in d.region
                if use_linear
                    v.pop_population_ann[yr, r] = p.pop_population[t, r] * fraction_timestep + p.pop_population[t - 1, r] * (1 - fraction_timestep)
                    v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r] * fraction_timestep + p.grw_gdpgrowthrate[t - 1, r] * (1 - fraction_timestep)
                elseif use_logburke
                    ## fully linear (partially because everything except for pop_population causes imaginary numbers (due to negative numbers)).
                    # linear
                    v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r] * fraction_timestep + p.grw_gdpgrowthrate[t - 1, r] * (1 - fraction_timestep)
                    v.pop_population_ann[yr, r] = p.pop_population[t, r] * fraction_timestep + p.pop_population[t - 1, r] * (1 - fraction_timestep)
                elseif use_logpopulation
                    # linear
                    v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r] * fraction_timestep + p.grw_gdpgrowthrate[t - 1, r] * (1 - fraction_timestep)

                    # log
                    v.pop_population_ann[yr, r] = p.pop_population[t, r]^fraction_timestep * p.pop_population[t - 1, r]^(1 - fraction_timestep)
                elseif use_logwherepossible
                    # linear
                    v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r] * fraction_timestep + p.grw_gdpgrowthrate[t - 1, r] * (1 - fraction_timestep)

                    # log
                    v.pop_population_ann[yr, r] = p.pop_population[t, r]^fraction_timestep * p.pop_population[t - 1, r]^(1 - fraction_timestep)
                else
                    error("NO INTERPOLATION METHOD SELECTED! Specify linear or logarithmic interpolation.")
                end

            end
        end
    end
end

function calc_gdp(p, v, d, t, annual_year)
    # setting the year for entry in lists.
    yr = annual_year - 2015 + 1 # + 1 because of 1-based indexing in Julia

    for r in d.region
        # eq.28 in Hope 2002
        if is_first(t)
            v.gdp_ann[yr, r] = p.gdp_0[r] * (1 + (v.grw_gdpgrowthrate_ann[yr,r] / 100))^(p.y_year_ann[yr] - p.y_year_0)

            v.grwnet_realizedgdpgrowth_ann[yr,r] = v.grw_gdpgrowthrate_ann[yr,r]
            v.gdp_leveleffect_ann[yr,r] = v.gdp_ann[yr,r]

            v.cons_consumption_ann[yr, r] = v.gdp_ann[yr, r] * (1 - p.save_savingsrate / 100)
            v.cons_percap_consumption_ann[yr, r] = v.cons_consumption_ann[yr, r] / v.pop_population_ann[yr, r]

        elseif isequal(gettime(t), 2030)
            v.grwnet_realizedgdpgrowth_ann[yr,r] = v.grw_gdpgrowthrate_ann[yr,r] - p.ge_growtheffects * p.isat_ImpactinclSaturationandAdaptation_ann[6,r]
            v.gdp_ann[yr, r] = v.gdp_ann[yr - 1, r] * (1 + (v.grwnet_realizedgdpgrowth_ann[yr,r] / 100))^(p.y_year_ann[yr] - p.y_year_ann[yr - 1])
            v.gdp_leveleffect_ann[yr,r] = v.gdp_leveleffect_ann[yr - 1, r] *  (1 + (v.grw_gdpgrowthrate_ann[yr,r] / 100))^(p.y_year_ann[yr] - p.y_year_ann[yr - 1])

            v.cons_consumption_ann[yr, r] = v.gdp_ann[yr, r] * (1 - p.save_savingsrate / 100)
            v.cons_percap_consumption_ann[yr, r] = v.cons_consumption_ann[yr, r] / v.pop_population_ann[yr, r]

            # let boundary take effect if pc consumption is in the neighbourhood of the boundary
            if p.use_convergence == 1.
                if v.cons_percap_consumption_ann[yr,r] >= v.cbabsn_pcconsumptionbound_neighbourhood
                    v.cbreg_regionsatbound_ann[yr,r] = 0.
                    v.cons_consumption_noconvergence_ann[yr,r] = v.cons_consumption_ann[yr,r]
                    v.cons_percap_consumption_noconvergence_ann[yr,r] = v.cons_consumption_noconvergence_ann[yr,r] / v.pop_population_ann[yr,r]
                else
                    # calculate the consumption level if there was no convergence system
                    v.cons_consumption_noconvergence_ann[yr,r] = v.cons_consumption_noconvergence_ann[yr - 1, r] * (1 + (v.grwnet_realizedgdpgrowth_ann[yr,r] / 100))^(p.y_year_ann[yr] - p.y_year_ann[yr - 1])
                    v.cons_percap_consumption_noconvergence_ann[yr,r] = v.cons_consumption_noconvergence_ann[yr,r] / v.pop_population_ann[yr,r]

                    # send the pc cconsumption on a logistic path convergenging against the bound
                    v.cons_percap_consumption_ann[yr,r] = v.cbabsn_pcconsumptionbound_neighbourhood - 0.5 * v.cbaux1_pcconsumptionbound_auxiliary1 +
                                        v.cbaux1_pcconsumptionbound_auxiliary1 * exp(v.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                                (v.cons_percap_consumption_noconvergence_ann[yr,r] - v.cbabsn_pcconsumptionbound_neighbourhood)) /
                                                                            (1 + exp(v.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                                (v.cons_percap_consumption_noconvergence_ann[yr,r] - v.cbabsn_pcconsumptionbound_neighbourhood)))

                # recalculate all variables accordingly
                    v.cons_consumption_ann[yr, r] = v.cons_percap_consumption_ann[yr,r] * v.pop_population_ann[yr,r]
                    v.gdp_ann[yr, r] = v.cons_consumption_ann[yr, r] / (1 - p.save_savingsrate / 100)
                    v.grwnet_realizedgdpgrowth_ann[yr,r] = 100 * ((v.gdp_ann[yr, r] / v.gdp_ann[yr - 1, r])^(1 / (p.y_year_ann[yr] - p.y_year_ann[yr - 1])) - 1)

                    v.cbreg_regionsatbound_ann[yr,r] = 1.

                end
            else
                if v.cons_percap_consumption_ann[yr,r] < p.cbabs_pcconsumptionbound
                    v.cons_percap_consumption_ann[yr,r] = p.cbabs_pcconsumptionbound

                    # recalculate all variables accordingly
                    v.cons_consumption_ann[yr, r] = v.cons_percap_consumption_ann[yr,r] * v.pop_population_ann[yr,r]
                    v.gdp_ann[yr, r] = v.cons_consumption_ann[yr, r] / (1 - p.save_savingsrate / 100)
                    v.grwnet_realizedgdpgrowth_ann[yr,r] = 100 * ((v.gdp_ann[yr, r] / v.gdp_ann[yr - 1, r])^(1 / (p.y_year_ann[yr] - p.y_year_ann[yr - 1])) - 1)

                    v.cbreg_regionsatbound_ann[yr,r] = 1.
                end
            end

        else
            lag = gettime(t) - gettime(t - 1)
            v.grwnet_realizedgdpgrowth_ann[yr,r] = v.grw_gdpgrowthrate_ann[yr,r] - p.ge_growtheffects * p.isat_ImpactinclSaturationandAdaptation_ann[yr - lag,r]
            v.gdp_ann[yr, r] = v.gdp_ann[yr - 1, r] * (1 + (v.grwnet_realizedgdpgrowth_ann[yr,r] / 100))^(p.y_year_ann[yr] - p.y_year_ann[yr - 1])
            v.gdp_leveleffect_ann[yr,r] = v.gdp_leveleffect_ann[yr - 1, r] *  (1 + (v.grw_gdpgrowthrate_ann[yr,r] / 100))^(p.y_year_ann[yr] - p.y_year_ann[yr - 1])

            v.cons_consumption_ann[yr, r] = v.gdp_ann[yr, r] * (1 - p.save_savingsrate / 100)
            v.cons_percap_consumption_ann[yr, r] = v.cons_consumption_ann[yr, r] / v.pop_population_ann[yr, r]

            # let boundary take effect if pc consumption is in the neighbourhood of the boundary
            if p.use_convergence == 1.
                if v.cons_percap_consumption_ann[yr,r] >= v.cbabsn_pcconsumptionbound_neighbourhood
                    v.cbreg_regionsatbound_ann[yr,r] = 0.
                    v.cons_consumption_noconvergence_ann[yr,r] = v.cons_consumption_ann[yr,r]
                    v.cons_percap_consumption_noconvergence_ann[yr,r] = v.cons_consumption_noconvergence_ann[yr,r] / v.pop_population_ann[yr,r]
                else
                    # calculate the consumption level if there was no convergence system
                    v.cons_consumption_noconvergence_ann[yr,r] = v.cons_consumption_noconvergence_ann[yr - 1, r] * (1 + (v.grwnet_realizedgdpgrowth_ann[yr,r] / 100))^(p.y_year_ann[yr] - p.y_year_ann[yr - 1])
                    v.cons_percap_consumption_noconvergence_ann[yr,r] = v.cons_consumption_noconvergence_ann[yr,r] / v.pop_population_ann[yr,r]

                    # send the pc cconsumption on a logistic path convergenging against the bound
                    v.cons_percap_consumption_ann[yr,r] = v.cbabsn_pcconsumptionbound_neighbourhood - 0.5 * v.cbaux1_pcconsumptionbound_auxiliary1 +
                                        v.cbaux1_pcconsumptionbound_auxiliary1 * exp(v.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                                (v.cons_percap_consumption_noconvergence_ann[yr,r] - v.cbabsn_pcconsumptionbound_neighbourhood)) /
                                                                            (1 + exp(v.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                                (v.cons_percap_consumption_noconvergence_ann[yr,r] - v.cbabsn_pcconsumptionbound_neighbourhood)))

                # recalculate all variables accordingly
                    v.cons_consumption_ann[yr, r] = v.cons_percap_consumption_ann[yr,r] * v.pop_population_ann[yr,r]
                    v.gdp_ann[yr, r] = v.cons_consumption_ann[yr, r] / (1 - p.save_savingsrate / 100)
                    v.grwnet_realizedgdpgrowth_ann[yr,r] = 100 * ((v.gdp_ann[yr, r] / v.gdp_ann[yr - 1, r])^(1 / (p.y_year_ann[yr] - p.y_year_ann[yr - 1])) - 1)

                    v.cbreg_regionsatbound_ann[yr,r] = 1.

                end
            else
                if v.cons_percap_consumption_ann[yr,r] < p.cbabs_pcconsumptionbound
                    v.cons_percap_consumption_ann[yr,r] = p.cbabs_pcconsumptionbound

                    # recalculate all variables accordingly
                    v.cons_consumption_ann[yr, r] = v.cons_percap_consumption_ann[yr,r] * v.pop_population_ann[yr,r]
                    v.gdp_ann[yr, r] = v.cons_consumption_ann[yr, r] / (1 - p.save_savingsrate / 100)
                    v.grwnet_realizedgdpgrowth_ann[yr,r] = 100 * ((v.gdp_ann[yr, r] / v.gdp_ann[yr - 1, r])^(1 / (p.y_year_ann[yr] - p.y_year_ann[yr - 1])) - 1)

                    v.cbreg_regionsatbound_ann[yr,r] = 1.
                end
            end
        end

        v.lgdp_gdploss_ann[yr,r] = v.gdp_leveleffect_ann[yr,r] - v.gdp_ann[yr,r]
    end

end

@defcomp GDP begin
    # GDP: Gross domestic product $M
    # GRW: GDP growth rate %/year
    region            = Index()
    year              = Index()

    # Variables
    gdp               = Variable(index=[time, region], unit="\$M")
    gdp_ann               = Variable(index=[year, region], unit="\$M")
    cons_consumption  = Variable(index=[time, region], unit="\$million")
    cons_consumption_ann  = Variable(index=[year, region], unit="\$million")
    cons_percap_consumption = Variable(index=[time, region], unit="\$/person")
    cons_percap_consumption_ann = Variable(index=[year, region], unit="\$/person")
    cons_percap_consumption_0 = Variable(index=[region], unit="\$/person")
    yagg_periodspan = Variable(index=[time], unit="year")

    # Parameters
    y_year_0          = Parameter(unit="year")
    y_year            = Parameter(index=[time], unit="year")
    y_year_ann        = Parameter(index=[year], unit="year")
    grw_gdpgrowthrate = Parameter(index=[time, region], unit="%/year") # From p.32 of Hope 2009
    grw_gdpgrowthrate_ann = Variable(index=[year, region], unit="%/year")                                   # interpolation
    gdp_0             = Parameter(index=[region], unit="\$M") # GDP in y_year_0
    save_savingsrate  = Parameter(unit="%", default=15.00) # pp33 PAGE09 documentation, "savings rate".
    pop0_initpopulation = Parameter(index=[region], unit="million person")
    pop_population    = Parameter(index=[time,region], unit="million person")
    pop_population_ann = Variable(index=[year, region], unit="million person")                              # interpolation

    # Saturation, used in impacts
    isat0_initialimpactfxnsaturation = Parameter(unit="unitless", default=20.0) # pp34 PAGE09 documentation
    isatg_impactfxnsaturation = Variable(unit="unitless")

    ###############################################
    # Growth Effects - additional variables and parameters
    ###############################################
    # parameters and variables for growth effects
    gdp_leveleffect   = Variable(index=[time, region], unit="\$M")
    gdp_leveleffect_ann   = Variable(index=[year, region], unit="\$M")
    isat_ImpactinclSaturationandAdaptation = Parameter(index=[time,region])
    isat_ImpactinclSaturationandAdaptation_ann = Parameter(index=[year,region])                              # from MarketDamagesBurke (if default, otherwise MarketDamages)
    lgdp_gdploss =  Variable(index=[time, region], unit="\$M")
    lgdp_gdploss_ann =  Variable(index=[year, region], unit="\$M")
    ge_growtheffects = Parameter(unit="none", default=0.)
    grwnet_realizedgdpgrowth = Variable(index=[time, region], unit="%/year")
    grwnet_realizedgdpgrowth_ann = Variable(index=[year, region], unit="%/year")
    # bound variables
    use_convergence = Parameter(unit="none", default=1.)
    cbabs_pcconsumptionbound = Parameter(unit="\$/person", default=740.65)
    cbabsn_pcconsumptionbound_neighbourhood = Variable(unit="\$/person")
    cbaux1_pcconsumptionbound_auxiliary1 = Variable(unit="none")
    cbaux2_pcconsumptionbound_auxiliary2 = Variable(unit="none")
    cbreg_regionsatbound = Variable(index=[time, region], unit="regions")
    cbreg_regionsatbound_ann = Variable(index=[year, region], unit="regions")
    cons_consumption_noconvergence = Variable(index=[time, region], unit="\$million")
    cons_consumption_noconvergence_ann = Variable(index=[year, region], unit="\$million")
    cons_percap_consumption_noconvergence = Variable(index=[time, region], unit="\$/person")
    cons_percap_consumption_noconvergence_ann = Variable(index=[year, region], unit="\$/person")
    # switch to overwrite the growth effects parameter with empirical distribution
    ge_empirical_distribution = Parameter(index=[draw], unit="none")
    ge_use_empiricaldistribution = Parameter(unit="none", default=0.)
    ge_seed_empiricaldistribution = Parameter(unit="none", default=1.)
    ###############################################

    function init(p, v, d)

        v.isatg_impactfxnsaturation = p.isat0_initialimpactfxnsaturation * (1 - p.save_savingsrate / 100)
        for rr in d.region
            v.cons_percap_consumption_0[rr] = (p.gdp_0[rr] / p.pop0_initpopulation[rr]) * (1 - p.save_savingsrate / 100)
        end
    end

    function run_timestep(p, v, d, t)

        interpolate_parameters_gdp(p, v, d, t)

        # Analysis period ranges - required for abatemnt costs and equity weighting, from Hope (2006)
        if is_first(t)
            ylo_periodstart = p.y_year_0
        else
            ylo_periodstart = (p.y_year[t] + p.y_year[t - 1]) / 2
        end

        if t.t == length(p.y_year)
            yhi_periodend = p.y_year[t]
        else
            yhi_periodend = (p.y_year[t] + p.y_year[t + 1]) / 2
        end

        v.yagg_periodspan[t] = yhi_periodend - ylo_periodstart

        if is_first(t)
            # calculate the lower consumption bound and the consumption level which triggers the convergence system
            v.cbabsn_pcconsumptionbound_neighbourhood = p.cbabs_pcconsumptionbound * 1.5

            # define auxiliary parameters for the convergence system to make formulas easier to read
            v.cbaux1_pcconsumptionbound_auxiliary1 = 2 * (p.cbabs_pcconsumptionbound - v.cbabsn_pcconsumptionbound_neighbourhood)
            v.cbaux2_pcconsumptionbound_auxiliary2 = 4 / v.cbaux1_pcconsumptionbound_auxiliary1

            # if the switch is set to one, overwrite the growth effects parameter with an empirical distribution, using the seed parameter
            if p.ge_use_empiricaldistribution == 1.
                Random.seed!(trunc(Int, p.ge_seed_empiricaldistribution))
                p.ge_growtheffects = p.ge_empirical_distribution[Random.rand(1:10^6)]
            elseif p.ge_use_empiricaldistribution == 2. # set negative rho values to zero if switch set to two
                Random.seed!(trunc(Int, p.ge_seed_empiricaldistribution))
                p.ge_growtheffects = max(p.ge_empirical_distribution[Random.rand(1:10^6)], 0)
            end
        end

        for r in d.region
            # eq.28 in Hope 2002
            if is_first(t)
                v.gdp[t, r] = p.gdp_0[r] * (1 + (p.grw_gdpgrowthrate[t,r] / 100))^(p.y_year[t] - p.y_year_0)

                v.grwnet_realizedgdpgrowth[t,r] = p.grw_gdpgrowthrate[t,r]
                v.gdp_leveleffect[t,r] = v.gdp[t,r]

                v.cons_consumption[t, r] = v.gdp[t, r] * (1 - p.save_savingsrate / 100)
                v.cons_percap_consumption[t, r] = v.cons_consumption[t, r] / p.pop_population[t, r]
            else
                v.grwnet_realizedgdpgrowth[t,r] = p.grw_gdpgrowthrate[t,r] - p.ge_growtheffects * p.isat_ImpactinclSaturationandAdaptation[t - 1,r]
                v.gdp[t, r] = v.gdp[t - 1, r] * (1 + (v.grwnet_realizedgdpgrowth[t,r] / 100))^(p.y_year[t] - p.y_year[t - 1])
                v.gdp_leveleffect[t,r] = v.gdp_leveleffect[t - 1, r] *  (1 + (p.grw_gdpgrowthrate[t,r] / 100))^(p.y_year[t] - p.y_year[t - 1])

                v.cons_consumption[t, r] = v.gdp[t, r] * (1 - p.save_savingsrate / 100)
                v.cons_percap_consumption[t, r] = v.cons_consumption[t, r] / p.pop_population[t, r]

                # let boundary take effect if pc consumption is in the neighbourhood of the boundary
                if p.use_convergence == 1.
                    if v.cons_percap_consumption[t,r] >= v.cbabsn_pcconsumptionbound_neighbourhood
                        v.cbreg_regionsatbound[t,r] = 0.
                        v.cons_consumption_noconvergence[t,r] = v.cons_consumption[t,r]
                        v.cons_percap_consumption_noconvergence[t,r] = v.cons_consumption_noconvergence[t,r] / p.pop_population[t,r]
                    else
                        # calculate the consumption level if there was no convergence system
                        v.cons_consumption_noconvergence[t,r] = v.cons_consumption_noconvergence[t - 1, r] * (1 + (v.grwnet_realizedgdpgrowth[t,r] / 100))^(p.y_year[t] - p.y_year[t - 1])
                        v.cons_percap_consumption_noconvergence[t,r] = v.cons_consumption_noconvergence[t,r] / p.pop_population[t,r]

                        # send the pc cconsumption on a logistic path convergenging against the bound
                        v.cons_percap_consumption[t,r] = v.cbabsn_pcconsumptionbound_neighbourhood - 0.5 * v.cbaux1_pcconsumptionbound_auxiliary1 +
                                            v.cbaux1_pcconsumptionbound_auxiliary1 * exp(v.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                                (v.cons_percap_consumption_noconvergence[t,r] - v.cbabsn_pcconsumptionbound_neighbourhood)) /
                                                                                (1 + exp(v.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                                (v.cons_percap_consumption_noconvergence[t,r] - v.cbabsn_pcconsumptionbound_neighbourhood)))

                    # recalculate all variables accordingly
                        v.cons_consumption[t, r] = v.cons_percap_consumption[t,r] * p.pop_population[t,r]
                        v.gdp[t, r] = v.cons_consumption[t, r] / (1 - p.save_savingsrate / 100)
                        v.grwnet_realizedgdpgrowth[t,r] = 100 * ((v.gdp[t, r] / v.gdp[t - 1, r])^(1 / (p.y_year[t] - p.y_year[t - 1])) - 1)

                        v.cbreg_regionsatbound[t,r] = 1.

                    end
                else
                    if v.cons_percap_consumption[t,r] < p.cbabs_pcconsumptionbound
                        v.cons_percap_consumption[t,r] = p.cbabs_pcconsumptionbound

                        # recalculate all variables accordingly
                        v.cons_consumption[t, r] = v.cons_percap_consumption[t,r] * p.pop_population[t,r]
                        v.gdp[t, r] = v.cons_consumption[t, r] / (1 - p.save_savingsrate / 100)
                        v.grwnet_realizedgdpgrowth[t,r] = 100 * ((v.gdp[t, r] / v.gdp[t - 1, r])^(1 / (p.y_year[t] - p.y_year[t - 1])) - 1)

                        v.cbreg_regionsatbound[t,r] = 1.
                    end
                end
            end

            v.lgdp_gdploss[t,r] = v.gdp_leveleffect[t,r] - v.gdp[t,r]

        end

        # calculate  for this specific year
        if is_first(t)
            for annual_year = 2015:(gettime(t))
                calc_gdp(p, v, d, t, annual_year)
            end
        else
            for annual_year = (gettime(t - 1) + 1):(gettime(t))
                calc_gdp(p, v, d, t, annual_year)
            end
        end

    end
end


text:
module eQuilibrator

using PyCall
using Unitful
using Measurements
using Crayons
using DocStringExtensions

include("utils.jl")
include("equilibrator_type.jl")
include("univariate_reaction.jl")
include("multivariate_reaction.jl")
include("multicompartment_reaction.jl")
include("formation.jl")

# export everything that isn't prefixed with _ (inspired by JuMP.jl, thanks!)
for sym in names(@__MODULE__, all = true)
    if sym in [Symbol(@__MODULE__), :eval, :include] || startswith(string(sym), ['_', '#'])
        continue
    end
    @eval export $sym
end

end # module


text:
# Packages required
using Hwloc
using .Threads
using Parameters
using LinearAlgebra
using Calculus
using DifferentialEquations
using ForwardDiff
using Distributions
using StatsBase
using Random
using SymPy
using PyPlot
using DataFrames
using CSV


text:
# Categorical distribution that works with arbitrary objects (Rather than consecutive ints)
struct ObjectCategorical
    cat::Categorical
    objs
    ObjectCategorical(objs::T) where {T <: AbstractArray} = new(Categorical(length(objs)), objs)
    ObjectCategorical(objs::T) where {T <: Tuple} = new(Categorical(length(objs)), [objs...])
end

Base.rand(d::ObjectCategorical, sz::Int...) = d.objs[rand(d.cat, sz...)]
function Distributions.logpdf(d::ObjectCategorical, x)
    if x isa AbstractArray
        x = mapslices(findfirst, d.objs .== reshape(x, 1,:), dims=1)
        return logpdf.(d.cat, x)
    else
        x = findfirst(d.objs .== x)
        return logpdf(d.cat, x)
    end

end

# Constant Layer
struct ConstantLayer{T}
    vec::T
end

Flux.@functor ConstantLayer
(m::ConstantLayer)(x::AbstractArray) = m.vec

## Useful functions
whiten(v) = whiten(v, mean(v), std(v))
whiten(v, μ, σ) = (v .- μ) ./ σ

to2D(W) = reshape(W, :, size(W, ndims(W))) # convert a multidimensional weight matrix to 2D

# Weighted mean aggregator
weighted_mean(weights) = (y) -> mean(y .* weights)

function LinearAlgebra.norm(grads::Flux.Zygote.Grads; p::Real = 2)
    v = []
    for θ in grads.params
        !isnothing(grads[θ]) && push!(v, norm(grads[θ] |> cpu, p))
    end
    norm(v, p)
end

## Early stopping
# Early stopping function that terminates training on validation error increase
function stop_on_validation_increase(π, 𝒫, 𝒟_val, loss; window=5)
    k = "validation_error"
    (infos) -> begin
        ve = loss(π, 𝒫, 𝒟_val) # Compute the validation error
        infos[end][k] = ve # store it
        N = length(infos)
        if length(infos) >= 2*window
            curr_window = mean([infos[i][k] for i=N-window+1:N])
            old_window = mean([infos[i][k] for i=N-2*window+1:N-window])
            return curr_window >= old_window # check if the error has gone up
        end
        false
    end
end

## Losses
function td_loss(;loss=Flux.mse, name=:Qavg, s_key=:s, a_key=:a, weight=nothing)
    (π, 𝒫, 𝒟, y; info=Dict()) -> begin
        Q = value(π, 𝒟[s_key], 𝒟[a_key])

        # Store useful information
        ignore() do
            info[name] = mean(Q)
        end

        loss(Q, y, agg = isnothing(weight) ? mean : weighted_mean(𝒟[weight]))
    end
end

function double_Q_loss(;name1=:Q1avg, name2=:Q2avg, kwargs...)
    l1 = td_loss(;name=name1, kwargs...)
    l2 = td_loss(;name=name2, kwargs...)

    (π, 𝒫, 𝒟, y; info=Dict()) -> begin
        .5f0*(l1(π.C.N1, 𝒫, 𝒟, y, info=info) + l2(π.C.N2, 𝒫, 𝒟, y, info=info))
    end
end

function multi_td_loss(;names, indices=1:length(names), kwargs...)
    ls = [td_loss(;name=name, kwargs...) for name in names]

    (π, 𝒫, 𝒟, ys; info=Dict()) -> begin
        mean([loss(net, 𝒫, 𝒟, y, info=info) for (loss, net, y) in zip(ls, π.networks[indices], ys)])
    end
end

td_error(π, 𝒟, y) = abs.(value(π, 𝒟[:s], 𝒟[:a])  .- y)

## Scheduling
struct LinearDecaySchedule{R<:Real} <: Function
    start::R
    stop::R
    steps::Int
end

function (schedule::LinearDecaySchedule)(i)
    rate = (schedule.start - schedule.stop) / schedule.steps
    val = schedule.start - i*rate
    val = max(schedule.stop, val)
end

function MultitaskDecaySchedule(steps::Int, task_ids; start = 1.0, stop = 0.1)
    schedule = LinearDecaySchedule(start, stop, steps)
    function val(i)
        taskindex = ceil(Int, i / steps)
        taskindex < 1 && return start
        taskindex > length(task_ids) && return stop
        taskid = task_ids[taskindex]
        used_steps = steps*sum(task_ids[1:taskindex-1] .== taskid)
        schedule(used_steps + mod1(i, steps))
    end
end

function logcompσ(x::Real)
    # Computes log(1 - sigmoid(x)) in a numerically stable way
    return logσ(x) - x
end


text:
using Test
using LogicCircuits

@testset "CNF file parser tests" begin

   circuit = zoo_cnf("8.cnf")

   @test circuit isa PlainLogicCircuit
   @test circuit isa Plain⋀Node
   @test circuit.children[1] isa Plain⋁Node
   @test circuit.children[end] isa Plain⋁Node
   @test num_variables(circuit) == 227
   @test num_nodes(circuit) == 1168

   circuit = zoo_cnf("iscas89/s386.scan.cnf")

   @test circuit isa PlainLogicCircuit
   @test circuit isa Plain⋀Node
   @test circuit.children[1] isa Plain⋁Node
   @test circuit.children[end] isa Plain⋁Node
   @test num_variables(circuit) == 172
   @test num_children(circuit) == 506

end

@testset "DNF file parser tests" begin

   circuit = zoo_dnf("8.dnf")

   @test circuit isa PlainLogicCircuit
   @test circuit isa Plain⋁Node
   @test circuit.children[1] isa Plain⋀Node
   @test circuit.children[end] isa Plain⋀Node
   @test num_variables(circuit) == 227
   @test num_nodes(circuit) == 1168

end


text:
macro inject(x::Symbol)
	t = eval(x)
	@assert isa(t, DataType)
	returnValue = Expr[]
	for i =1:length(t.types)
		push!(returnValue, Expr(:(::), t.names[i], t.types[i].name.name))
	end
	esc((Expr(:block, returnValue...)))
end
type LOL
	x::Int
	y::Int
end
type Test
	c::String
end
type LOL3
	@inject Test
	foo::Array{Float32, 1}
end
type LOL2
	@inject LOL
	@inject LOL3
end
println(names(LOL2))

type FunctionT{T}
	func::Function
	args::T
end

function simon(a::Int, b::Float32, c::Array{Int,1})
	println(a)
	println(b)
	println(c)
end
function callF(f::FunctionT)
	apply(f.func, f.args)
end
function simon(a::Int, b::Float32, c::Int)
	println(a*b*c)
end
function foo(f::FunctionT{(Int, Float32, Int)})
	callF(f)
end

function foo(f::FunctionT{(Int, Float32, Array{Int, 1})})
	callF(f)
end
a = FunctionT(simon, (1, 1.0f0, [1,2,3,4]))
b = FunctionT(simon, (1, 1.0f0, 1))
foo(b)

max_stars_repo_name: dadarakt/jamon

max_issues_repo_name: dadarakt/jamon

max_forks_repo_name: dadarakt/jamon


text:
using SVG

using Images

bw = load(ARGS[1])

height, width = size(bw)
hscale = 5
wscale = 5

svg = Svg(2wscale * width, 2hscale * height)

plusminus = repeat([1,-1], wscale)

function row(yo, pxs)
	tps = Vector{Float64}(undef, 2wscale * width)
	xs = 1
	xe = 2wscale
	for x in 1:width
		tps[xs:xe] = 1.2*hscale*(1-pxs[x]) * rand(2wscale) .* plusminus
		xs += 2wscale
		xe += 2wscale
	end

	line = SVG.Polyline(length(tps), style=SVG.Style(strokewidth=0.5))
	for x in 1:length(tps)
		line.points[x] = (x, yo+tps[x])
	end
	line
end

for h in 1:height
	push!(svg.objects, row(2hscale*h, map(Float64, map(Gray, bw[h, :]))))
end

write(stdout, svg)

max_stars_repo_name: lawless-m/jsvg


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastSurface
using GeometricExamples

const Δt     = 2.
const ntime  = 12500
const nx     = 200
const ny     = 200
const nsave  = 25
const nplot  = 5
const run_id = "poincare_2nd_fast_vprk_dt2"

include("guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_2nd.jl")


text:
using Example

Example.hello("CORE 2021")
Example.domath(10)