text:
__precompile__(true)
# Name: ODLCommonTools
# Description: Common Tools for ODL
#   Performs the function of forward declaring abstract types

module ODLCommonTools
using ArrayViews
import ArrayViews.view  # use ArrayViews rather than Base.view
import Base.==

include("topo.jl")

import Base.show
import Base.isless
import Base.copy
import Base.copy!
export AbstractSolutionData, AbstractParamType, Abstract3DArray, Abstract4DArray
export AbstractMesh, AbstractCGMesh, AbstractDGMesh, AbstractOperator,
       AbstractSharedFaceData
export Boundary, Interface, getElementL, getFaceL, NullBoundary, NullInterface
export AbstractFunctional, AbstractIntegralFunctional,
       AbstractBoundaryFunctional, setupFunctional, getParallelData,
       getParallelDataString, getParallelDataEnum, CompositeFunctional,
       PARALLEL_DATA_NONE, PARALLEL_DATA_FACE, PARALLEL_DATA_ELEMENT
export Boundary
export Interface
export BCType, BCType_revm, BCType_revq, SRCType, FluxType, FluxType_revm,
       FluxType_diff, FluxType_revq, FunctionalType
export calcDiffElementArea
export functorThatErrors, functorThatErrors_revm

# topo.jl
export ElementTopology3, ElementTopology2, ElementTopology

# eqn_copy.jl
export copyForMultistage!

# eqn_deepcopy.jl
export eqn_deepcopy, eqn_deepcopy_fields

export ROView, ro_sview, ROVector, ROMatrix, ROArray

# misc.jl
export prepend_path, append_path, split_fname

# io.jl
export write_binary, read_binary!, writeSolutionFiles, readSolutionFiles

# lapack.jl
export getrf!, getrs2!, BlasInt, laswp!, applyIpiv!, solve_suitesparse,
       UMFPACK_A, UMFPACK_At, UMFPACK_Aat

export sview

# null_array.jl
export NullArray

# empty_array.jl
export EmptyArray

"""

  This abstract type is the supertype for all the objects that store the
  solution data. Every physics module should implement its own subtype.

  **Static parameters**

   * Tsol: datatype of solution variables
   * Tres: datatype of the mesh variables

  See the [AbstractSolutionData](@ref) for the description of everything this
  type must implement.

"""
abstract type AbstractSolutionData{Tsol, Tres} end

"""
  This abstract type is the supertype for all mesh objects.  Every interface to
  a mesh software should define its own implementation.

  **Static parameters**

   * Tmsh: datatype of the mesh data (coordinates, mapping to/from parametric
          space, mapping jacobian).

  See the [AbstractMesh](@ref) for the description of everything this
  type must implement.

"""
abstract type AbstractMesh{Tmsh} end

@doc """
### ODLCommonTools.AbstractCGMesh

  The abstrac type is the supertype of all continuous Galerkin meshes
"""->
abstract type AbstractCGMesh{Tmsh} <: AbstractMesh{Tmsh} end

@doc """
### ODLCommonTools.AbstractDGGMesh

  The abstract type is the supertype of all discontinuous Galerkin meshes
"""->
abstract type AbstractDGMesh{Tmsh} <: AbstractMesh{Tmsh} end

"""
  Abstract type for all discretization operators (SBP, FE, etc.)
"""
abstract type AbstractOperator{T} end

"""
  This abstract type is the supertype for all Param objects, which hold values
  needed for the computation in a place that is fast to access.

   The Param type is also useful for dispatching to low level functions which
   the AbstractSolutionData might not be passed (depending on the organization
   of the physics module.

  **Static Parameters**:

   * Tdim: the dimensionality of the equation being solved (2d or 3d usually)
"""
abstract type AbstractParamType{Tdim} end

@doc """
### ODLCommonTools.Abstract3DArray

  A typealias useful for specify a 3 dimensional AbstractArray without
  specifying the element datatype.

"""->
const Abstract3DArray{T} = AbstractArray{T, 3}

"""
  A typealias for any 4D array.  Element is the static parameter
"""
const Abstract4DArray{T} =  AbstractArray{T, 4}

"""
Abstract supertype of all SharedFaceData implemenations, used to storing
the data needed for parallel communication

  **Static Parameters**

   * Tsol: the datatype of the solution variables
"""
abstract type AbstractSharedFaceData{Tsol} end

function show(io::IO, obj::AbstractMesh)
  println(io, "$(typeof(obj)) with $(obj.numEl) p = $(obj.order) elements")
end

function show(io::IO, obj::AbstractSolutionData)
  println(io, "$(typeof(obj)) with $(size(obj.q, 3)) elements")
end

function show(io::IO, obj::AbstractOperator)
  println(io, "degree $(obj.degree) $(typeof(obj))")
end

"""

  Copies only the essential data from one AbstractSolutionData to another.

  Currently copies:
    * q
    * q_vec
    * res
    * res_vec
    * shared_data

  Implementation Notes:
    Avoids double copying when q and q_vec alias
"""
function copyForMultistage!(eqn_dest::AbstractSolutionData, eqn_src::AbstractSolutionData)

  copy!(eqn_dest.q, eqn_src.q)

  # slight optimization- checking if we need to copy q_vec at all,
  #   in case q_vec & q point to the same thing
  if pointer(eqn_src.q) != pointer(eqn_src.q_vec)
    copy!(eqn_dest.q_vec, eqn_src.q_vec)
  end

  # This is needed because we suspect copy! will only work down one level
  num_inner_arrays = length(eqn_src.shared_data)
  for i = 1:length(eqn_src.shared_data)
    copy!(eqn_dest.shared_data[i], eqn_src.shared_data[i])
  end

  copy!(eqn_dest.res, eqn_src.res)
  # same slight optimization as above for q & q_vec
  if pointer(eqn_src.res) != pointer(eqn_src.res_vec)
    copy!(eqn_dest.res_vec, eqn_src.res_vec)
  end

  return nothing

end

@doc """

Not in place copy function for AbstractSolutionData.
    We want this to throw an error because copy! is better.

"""->
function copy(eqn_src::AbstractSolutionData)

  throw(ErrorException("Use copy!(eqn_dest, eqn_src) for AbstractSolution data, not copy"))

  return nothing

end

@doc """
### ODLCommonTools.Boundary

Used to identify boundary faces in a finite-element grid.

**Fields**

* `element` : index of the element to which the boundary face belongs
* `face` : the face index of the boundary (local index to the element)

**Example**

To mark face 2 of element 7 to be a boundary face, use `Boundary(7,2)`

"""->
struct Boundary
  element::UInt32
  face::UInt8
end

"""
  Boundary(0, 0).  Useful as a default value for function arguments
"""
global const NullBoundary = Boundary(0,0)

@doc """
### ODLCommonTools.Interface

Used to identify interfaces between elements in a finite-element grid.

**Fields**

* `elementL` : index of the so-called left element in the pair
* `elementR` : index of the so-called right element in the pair
* `faceL` : the face index of the interface with respect to the left element
* `faceR` : the face index of the interface with respect to the right element
* `orient` : orientation of the 'right' element relative to the 'left'

**Example**

Consider an interface between elements 2 and 5.  Suppose the interface is on
face 1 of element 2 and face 3 of element 5.  Furthermore, suppose element 5 has
orientation 1 relative to element 1 (defintion of orientation TBD).  This can be
indicated as `Interface(2,5,1,3,1)`

"""->
struct Interface
  elementL::UInt32
  elementR::UInt32
  faceL::UInt8
  faceR::UInt8
  orient::UInt8
end

global const NullInterface = Interface(0, 0, 0, 0, 0)

# small interface for Boundary and Interface

"""
  This function returns either the `element` field of a Boundary or the
  `elementL` field of an interface.
"""
function getElementL(bndry::Boundary)
  return bndry.element
end

function getElementL(iface::Interface)
  return iface.elementL
end

"""
  This function returns either the `face` field of a Boundary or the
  `faceL` field of an Interface
"""
function getFaceL(bndry::Boundary)
  return bndry.face
end

function getFaceL(iface::Interface)
  return iface.faceL
end

"""
  Abstract supertype of all boundary condition functors
"""
abstract type BCType end

"""
  Equality test for two BCTypes of the same type.  Tests that all fields
  are equal
"""
function (==)(x::T, y::T) where {T <: BCType}

  fnames = fieldnames(T)
  val = true
  for fname in fnames
    val = val && (getfield(x, fname) == getfield(y, fname))
  end

  return val
end

"""
  Abstract supertype of all boundary condition functors that compute the
  reverse mode with respect to the metrics
"""
abstract type BCType_revm end

"""
  Abstract supertype of all boundary condition functors that compute the
  reverse mode with respect to the solution
"""
abstract type BCType_revq end

"""
  Abstract supertype of all source term functors
"""
abstract type SRCType end

"""
  Abstract supertype of all 2 point numerical flux functions.
"""
abstract type FluxType end # functor DG flux abstract type

"""
  Like [`Fluxtype`](@ref) but for functors that compute the reverse mode
  with respect to the metrics.
"""
abstract type FluxType_revm end # functor type for reverse mode of DG interface fluxes w.r.t mesh metrics

"""
  Like [`FluxType`](@ref), but for flux functions that compute the
  jacobian of the flux with respect to the left and right solutions.
"""
abstract type FluxType_diff end

"""
  Like [`FluxType`](@ref) but for functors that compute the reverse mode
  with respect to the solution q
"""
abstract type FluxType_revq end

abstract type FunctionalType end # functor for functional abstract type

"""
  Show method for Boundary objects
"""
function show(io::IO, object::Boundary)
  print(io, "Boundary element, face = ", object.element, ", ", object.face)
end

"""
  Show method for Interface objects
"""
function show(io::IO, obj::Interface)
  print(io, "Interface elementL, elementR, faceL, faceR, orient = ",
        obj.elementL, ", ",obj.elementR, ", ", obj.faceL, ", ", obj.faceR,
        ", ",obj.orient)
end

"""
  Compare boundaries first by element, then by face
"""
function isless(a::Boundary, b::Boundary)
  if a.element < b.element
    return true
  elseif a.element > b.element
    return false
  elseif a.face > b.face     # elements are same
    return true
  else
    return false
  end

end

"""
  Compare Interfaces, first by elementL, then by elementR
"""
function isless(a::Interface, b::Interface)
  if a.elementL < b.elementL
    return true
  elseif a.elementL > b.elementL
    return false
  end

  if a.elementR < b.elementR
    return true
  else
    return false
  end
end

@doc """
### ODLCommonTools.calcDiffElementArea

  This function calculates the equivalent of the differential area of an element
  at a point in the x and y directions?

  Inputs:
    nrm: a normal vector in the parametric coordinate system of the element
         at the node
    dxidx: the 2x2 matrix dxi/dx at the node

  Inputs/Outputs:
    workvec: [dx, dy]

"""->
function calcDiffElementArea(nrm::AbstractArray{T,1},
                            dxidx::AbstractArray{T2,2},
                            workvec::AbstractArray{T3,1}) where {T, T2, T3}
  fill!(workvec, zero(T3))
  for di1 = 1:size(nrm,1)
    for di2 = 1:size(nrm,1)
      workvec[di2] += nrm[di1]*dxidx[di1,di2]
    end
  end
  return norm(workvec)
end

# it would be better if this used @boundscheck
@doc """
### Utils.safe_views

  This bool value controls whether the function named sview refers to
  view or unsafe_view from the ArrayViews package
"""->
global const safe_views = true
if safe_views
  global const sview = ArrayViews.aview
else
  global const sview = ArrayViews.unsafe_aview
end

"""
  This macro eliminates blocks of code if safe_views is false
"""
macro ifsafeview(ex)

  if safe_views
    return quote
      $(esc(ex))
    end
  else
    return nothing
  end
end  # end macro

@doc """
### ODLCommonTools.functorThatErrors

  This functor will error.

  It is intended to be used to initialize functor fields in eqn to some value
    so that eqn_deepcopy can operate upon them.

  This should never be called; if it is called, then the field of eqn has not been
    properly defined for its desired usage.

  Inputs: none
  Outputs: none

"""->
mutable struct functorThatErrors <: FluxType
end

function (obj::functorThatErrors)(params::AbstractParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F::AbstractVector{Tres}) where {Tsol, Tres, Tmsh}

  error("Default functor has been called. You have not properly initialized something.")
  return nothing
end

@doc """
### ODLCommonTools.functorThatErrors_revm

  This functor will error.

  It is intended to be used to initialize functor fields in eqn to some value
    so that eqn_deepcopy can operate upon them.

  This should never be called; if it is called, then the field of eqn has not been
    properly defined for its desired usage.

  Inputs: none
  Outputs: none

"""->
mutable struct functorThatErrors_revm <: FluxType_revm
end

function (obj::functorThatErrors_revm)( params::AbstractParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F::AbstractVector{Tres}) where {Tsol, Tres, Tmsh}

  error("Default functor has been called. You have not properly initialized something.")
  return nothing
end

include("types.jl")
include("abstract_functional.jl")
include("eqn_deepcopy.jl")
include("ro_view.jl")
include("io.jl")
include("misc.jl")
include("getAllTypeParams.jl")
include("lapack.jl")
include("null_array.jl")
include("empty_array.jl")
end     # module


text:
using Printf # Needed to use printf
using Statistics

include("variables.jl")
include("strings.jl")
include("conditionals.jl")
include("looping.jl")
include("arrays.jl")
include("tuples.jl")
include("dictionaries.jl")
include("sets.jl")
include("functions.jl")
include("anonymous_functions.jl")
include("maths.jl")
include("enums.jl")
include("symbols.jl")
include("structs.jl")
include("abstract_types.jl")
include("exceptions.jl")
include("file_io.jl")
include("macros.jl")
include("map_filter_reduce.jl")

function main()
    variable_examples()
    string_examples()
    conditional_examples()
    looping_examples()
    array_examples()
    tuple_examples()
    dictionary_examples()
    set_examples()
    function_examples()
    anonymous_function_examples()
    maths_examples()
    enum_examples()
    symbol_examples()
    struct_examples()
    abstract_type_examples()
    exception_examples()
    file_io_examples()
    macro_examples()
    map_filter_reduce_examples()
end

main()

ext: jl

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
const Layoutable = Union{LAxis, LObject, GridLayout}

defaultlayout(layoutable::Layoutable) = ProtrusionLayout(layoutable)

function align_to_bbox!(layoutable::Layoutable, bbox)
    layoutable.layoutnodes.suggestedbbox[] = bbox
end

computedsizenode(layoutable::Layoutable) = layoutable.layoutnodes.computedsize
protrusionnode(layoutable::Layoutable) = layoutable.layoutnodes.protrusions

function Base.getproperty(layoutable::T, s::Symbol) where T <: Layoutable
    if s in fieldnames(T)
        getfield(layoutable, s)
    else
        layoutable.attributes[s]
    end
end

function Base.setproperty!(layoutable::T, s::Symbol, value) where T <: Layoutable
    if s in fieldnames(T)
        setfield!(layoutable, s, value)
    else
        layoutable.attributes[s][] = value
    end
end

function Base.propertynames(layoutable::T) where T <: Layoutable
    [fieldnames(T)..., keys(layoutable.attributes)...]
end

Base.Broadcast.broadcastable(l::Layoutable) = Ref(l)


text:
using Math
using Test

@testset "Math" begin
    # Write your tests here.
end


text:
mutable struct Agent{T<:AbstractFloat}
    u::Dict
    p::Particle{T}
    cnti::Signed
    cntk::Signed
end

function create_agents(num_agents::Signed, init_width::T,
    init_center::AbstractArray{T, 1}, u) where {T<:AbstractFloat}
    agt = []
    for i = 1:num_agents
        init_pos = init_width*(rand(2)-0.5) .+ init_center
        init_vel = 2π*rand()
        p = Particle(init_pos..., init_vel)
        cnti = div(num_agents*(num_agents-1), 2)
        cntk = num_agents - 1
        a = Agent(deepcopy(u), p, cnti, cntk)
        push!(agt, a)
    end
    return agt
end

function detected_crossing!(agent, color)
    agent.cntk -= 1
    if agent.cntk != 0 && agent.cnti == agent.cntk*(agent.cntk + 1)/2 + 1
        agent.cnti -= 1
        for key in keys(agent.u)
            agent.u[key] = 0
        end
        agent.u[color] = 1
    else
        for key in keys(agent.u)
            agent.u[key] = 1
        end
    end
    return (color, agent.cntk)
end

function received_broadcast!(agent, broadcast)
    agent.cnti -= 1
    if agent.cnti == agent.cntk*(agent.cntk - 1)/2
        for key in keys(agent.u)
            agent.u[key] = 1
        end
    elseif agent.cnti == agent.cntk*(agent.cntk + 1)/2 + 1 && agent.cntk == broadcast[2]
        agent.cnti -= 1
        for key in keys(agent.u)
            agent.u[key] = 0
        end
        agent.u[broadcast[1]] = 1
    # Below to handle cyclic graphs with d ≤ 2
    # (but can be added to chain case too, as it will not make a difference there)
    elseif agent.cnti > agent.cntk*(agent.cntk - 1)/2 && agent.cntk > broadcast[2]
        for key in keys(agent.u)
            agent.u[key] = 1
        end
        agent.u[broadcast[1]] = 0
    end
end

max_stars_repo_name: gbarsih/Representations-in-Robotics

max_issues_repo_name: gbarsih/Representations-in-Robotics

max_forks_repo_name: gbarsih/Representations-in-Robotics


text:
__precompile__()

module SparseRF

using Compat

import Base: length, convert, promote_rule, show, start, next, done

export Leaf, Node, Ensemble, print_tree, depth, build_stump, build_tree,
       prune_tree, apply_tree, apply_tree_proba, nfoldCV_tree, build_forest,
       apply_forest, apply_forest_proba, nfoldCV_forest, build_adaboost_stumps,
       apply_adaboost_stumps, apply_adaboost_stumps_proba, nfoldCV_stumps,
       majority_vote, ConfusionMatrix, confusion_matrix, mean_squared_error,
       R2, _int

# ScikitLearn API
export DecisionTreeClassifier, DecisionTreeRegressor, RandomForestClassifier,
       RandomForestRegressor, AdaBoostStumpClassifier,
       # Should we export these functions? They have a conflict with
       # DataFrames/RDataset over fit!, and users can always
       # `using ScikitLearnBase`.
       predict, predict_proba, fit!, get_classes

#####################################
##### Compatilibity Corrections #####

const Range1 = Range
_int(x) = map(y->round(Integer, y), x)
float(x) = map(Float64, x)

_squeeze(m::Matrix, i::Int) = squeeze(m, i)
_squeeze(v::Vector, i::Int) = v
neg(arr) = map(!, arr)

###########################
########## Types ##########

const NO_BEST=(0,0)

immutable Leaf{S}
    majority::S
    values::Vector{S}
end

immutable Node{S,T}
    featid::Int
    featval::T
    left::Union{Leaf{S},Node{S, T}}
    right::Union{Leaf{S},Node{S, T}}
end

@compat const LeafOrNode{S,T} = Union{Leaf{S},Node{S, T}}

immutable Ensemble{S, T}
    trees::Vector{Node{S, T}}
end

immutable UniqueRanges
    v::AbstractVector
end

start(u::UniqueRanges) = 1
done(u::UniqueRanges, s) = done(u.v, s)
next(u::UniqueRanges, s) = (val = u.v[s];
                            t = searchsortedlast(u.v, val, s, length(u.v), Base.Order.Forward);
                            ((val, s:t), t+1))

##############################
########## Includes ##########

include("measures.jl")
include("classification.jl")
include("regression.jl")
include("scikitlearnAPI.jl")

#############################
########## Methods ##########

length(leaf::Leaf) = 1
length(tree::Node) = length(tree.left) + length(tree.right)
length(ensemble::Ensemble) = length(ensemble.trees)

depth(leaf::Leaf) = 0
depth(tree::Node) = 1 + max(depth(tree.left), depth(tree.right))

function print_tree(leaf::Leaf, depth=-1, indent=0)
    matches = find(leaf.values .== leaf.majority)
    ratio = string(length(matches)) * "/" * string(length(leaf.values))
    println("$(leaf.majority) : $(ratio)")
end

function print_tree(tree::Node, depth=-1, indent=0)
    if depth == indent
        println()
        return
    end
    println("Feature $(tree.featid), Threshold $(tree.featval)")
    print("    " ^ indent * "L-> ")
    print_tree(tree.left, depth, indent + 1)
    print("    " ^ indent * "R-> ")
    print_tree(tree.right, depth, indent + 1)
end

function show(io::IO, leaf::Leaf)
    println(io, "Decision Leaf")
    println(io, "Majority: $(leaf.majority)")
    print(io,   "Samples:  $(length(leaf.values))")
end

function show(io::IO, tree::Node)
    println(io, "Decision Tree")
    println(io, "Leaves: $(length(tree))")
    print(io,   "Depth:  $(depth(tree))")
end

function show(io::IO, ensemble::Ensemble)
    println(io, "Ensemble of Decision Trees")
    println(io, "Trees:      $(length(ensemble))")
    println(io, "Avg Leaves: $(mean([length(tree) for tree in ensemble.trees]))")
    print(io,   "Avg Depth:  $(mean([depth(tree) for tree in ensemble.trees]))")
end

end # module


text:
using BeaconK8sUtilities
using Documenter

DocMeta.setdocmeta!(BeaconK8sUtilities, :DocTestSetup, :(using BeaconK8sUtilities);
                    recursive=true)

makedocs(; modules=[BeaconK8sUtilities],
         repo="https://github.com/beacon-biosignals/BeaconK8sUtilities.jl/blob/{commit}{path}#{line}",
         sitename="BeaconK8sUtilities.jl",
         format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://beacon-biosignals.github.io/BeaconK8sUtilities.jl",
                                assets=String[], ansicolor=true),
         pages=["Home" => "index.md", "Utilities" => "utilities.md",
                "Templates" => "templates.md", "Preferences" => "preferences.md"])

deploydocs(; repo="github.com/beacon-biosignals/BeaconK8sUtilities.jl", push_preview=true,
           devbranch="main")


text:
#THIS IS THE MOST HORRIBLE FUNCTION I HAVE EVER CREATED!!!
#extracts only atoms with projections
function extract_atoms(atoms_block::T, proj_block::T, cell::Mat3{LT}, spinors=false) where {T<:InputData,LT<:Length}
    if atoms_block.name == :atoms_cart
        cell = Mat3(Matrix(1.0Ang*I, 3, 3))
    end
    projections_ = proj_block.data
    atoms = atoms_block.data
    t_start = 1
    out_ats = Atom{Float64, LT}[]
    if projections_ != nothing
        t_ats = Atom{Float64, LT}[]
        for (pos_at, pos) in atoms
            for ps in pos
                for (proj_at, projs) in projections_
                    if proj_at != pos_at
                        continue
                    end
                    for proj in projs
                        size = spinors ? 2*orbsize(proj) : orbsize(proj)
                        push!(t_ats, Atom(pos_at, element(pos_at), cell * ps, ps, projections=[Projection(orbital(proj), t_start, t_start + size - 1)]))
                        t_start += size
                    end
                end
            end
        end
        for (pos_at, pos) in atoms
            for ps in pos
                same_ats = Atom{Float64, LT}[]
                for at in t_ats
                    if position_cart(at) == cell * ps
                        push!(same_ats, at)
                    end
                end
                if isempty(same_ats)
                    continue
                end
                if length(same_ats) > 1
                    for at in same_ats[2:end]
                        push!(projections(same_ats[1]), projections(at)[1])
                    end
                end
                push!(out_ats, same_ats[1])
            end
        end
    else
        for (pos_at, pos) in atoms
            for p in pos
                push!(out_ats, Atom(pos_at, element(pos_at), cell * p, p, projections=:random))
            end
        end
    end
    return out_ats
end

function extract_structure(name, cell_block::T, atoms_block::T, projections_block::T, spinors=false) where T <: Union{InputData, Nothing}
    if atoms_block == nothing || cell_block == nothing
        return nothing
    end
    if cell_block.option == :bohr
        cell = cell_block.data' .* 1a₀
    else
        cell = cell_block.data' .* 1Ang
    end

    atoms = extract_atoms(atoms_block, projections_block, cell, spinors)
    return Structure(name, Mat3(cell), atoms)
end

function wan_read_input(::Type{T}, f::IO) where T
    flags = Dict{Symbol,Any}()
    data  = Vector{InputData}()
    atoms_block = nothing
    cell_block  = nothing
    proj_block  = nothing
    line = readline(f)
    while !eof(f)
        @label start_label

        if occursin("!", line) || line == "" || occursin("end", lowercase(line)) || occursin("#", line)
            line = readline(f)
            continue
        end

        if occursin("begin", lowercase(line))
            block_name = Symbol(split(lowercase(line))[end])

            if block_name == :projections
                proj_dict = Dict{Symbol,Array{Symbol,1}}()
                line      = readline(f)
                while !occursin("end", lowercase(line))
                    if occursin("!", line) || line == ""
                        line = readline(f)
                        continue
                    end
                    if occursin("random", line)
                        proj_block = InputData(:projections, :random, nothing)
                        line = readline(f)
                        break
                    else
                        split_line      = strip_split(line, ':')
                        atom            = Symbol(split_line[1])
                        projections     = [Symbol(proj) for proj in strip_split(split_line[2], ';')]
                        proj_dict[atom] = projections
                        line = readline(f)
                    end
                end
                proj_block = InputData(:projections, :none, proj_dict)
                @goto start_label

            elseif block_name == :kpoint_path
                line = readline(f)
                k_path_array = Array{Tuple{Symbol,Array{T,1}},1}()
                while !occursin("end", lowercase(line))
                if occursin("!", line) || line == ""
                    line = readline(f)
                    continue
                end
                split_line = split(line)
                push!(k_path_array, (Symbol(split_line[1]), parse_string_array(T, split_line[2:4])))
                push!(k_path_array, (Symbol(split_line[5]), parse_string_array(T, split_line[6:8])))
                line = readline(f)
            end
            push!(data, InputData(:kpoint_path, :none, k_path_array))
            @goto start_label

            elseif block_name == :unit_cell_cart
                line = readline(f)
                if length(split(line)) == 1
                    option = Symbol(strip(lowercase(line)))
                    line = readline(f)
                else
                    option = :ang
                end
                cell_param = Matrix{T}(undef, 3, 3)
                for i = 1:3
                    cell_param[i, :] = parse_line(T, line)
                    line = readline(f)
                end
                cell_block = InputData(:unit_cell_cart, option, Mat3(cell_param))
                # line = readline(f)
                @goto start_label

            elseif block_name == :atoms_frac || block_name == :atoms_cart
                line   = readline(f)
                atoms  = Dict{Symbol, Array{Point3{T}, 1}}()
                option = :ang
                while !occursin("end", lowercase(line))
                    if occursin("!", line) || line == ""
                        line = readline(f)
                        continue
                    end
                    if length(split(line)) == 1
                        option = Symbol(strip(line))
                        line = readline(f)
                        continue
                    end
                    split_line = strip_split(line)
                    atom       = Symbol(split_line[1])
                    position   = Point3(parse_string_array(T, split_line[2:4]))
                    if !haskey(atoms,atom)
                        atoms[atom] = [position]
                    else
                        push!(atoms[atom], position)
                    end
                    line = readline(f)
                end
                atoms_block = InputData(block_name, option, atoms)
                @goto start_label

            elseif block_name == :kpoints
                line     = readline(f)
                k_points = Array{Array{T,1},1}()
                while !occursin("end", lowercase(line))
                    if line == ""
                        line = readline(f)
                        continue
                    end
                    push!(k_points, parse_line(T, line))
                    line = readline(f)
                end
                push!(data, InputData(:kpoints, :none, k_points))
                @goto start_label
            end

        else
            if occursin("mp_grid", line)
                flags[:mp_grid] = parse_string_array(Int, split(line)[end-2:end])
            else
                flag, val = wan_parse_flag_line(line)
                flags[flag] = val
            end
        end
        line = readline(f)
    end
    return flags, data, atoms_block, cell_block, proj_block
end

wan_read_input(f::IO) = wan_read_input(Float64, f)

"""
    wan_read_input(filename::String, T=Float64; runcommand= Exec(""), run=true, exec=Exec("wannier90.x"), structure_name="NoName")

Reads a `DFInput{Wannier90}` and the included `Structure` from a WANNIER90 input file.
"""
function wan_read_input(filename::String, T=Float64; execs=[Exec("wannier90.x")], run=true, structure_name="NoName")
    flags = Dict{Symbol,Any}()
    data  = Vector{InputData}()
    atoms_block = nothing
    cell_block  = nothing
    proj_block  = nothing
    open(filename,"r") do f
    	flags, data, atoms_block, cell_block, proj_block = wan_read_input(T, f)
	end
    structure = extract_structure(structure_name, cell_block, atoms_block, proj_block, get(flags, :spinors, false))
    dir, file = splitdir(filename)
    flags[:preprocess] = hasflag(getfirst(x->x.exec == "wannier90.x", execs), :pp)  ? true : false
    return DFInput{Wannier90}(splitext(file)[1], dir, flags, data, execs, run), structure
end

function wan_parse_flag_line(line::String)
    split_line = strip_split(line, '=')
    flag       = Symbol(split_line[1])
    flagtyp   = flagtype(Wannier90, flag)
    value      = strip(lowercase(split_line[2]), '.')
    if flagtyp != String
	    value = replace(value, "d" => "e")
    end
    spl_val    = split(value)
    if eltype(flagtyp) == Bool
        for i = 1:length(spl_val)
            if spl_val[i] == "t"
                spl_val[i] = "true"
            elseif spl_val[i] == "f"
                spl_val[i] = "false"
            end
        end
    end
    parsed_val = flagtyp == String ? spl_val : parse.(eltype(flagtyp), spl_val)
    if length(parsed_val) == 1
        val = parsed_val[1]
    else
        val = parsed_val
    end
    return flag, val
end

function projections_string(at::AbstractAtom)
	prjs = projections(at)
	str = "$(name(at)): $(prjs[1].orb.name)"
    if length(prjs) > 1
        for proj in prjs[2:end]
            str *= ";$(proj.orb.name)"
        end
    end
    return str
end

function wan_write_projections(f::IO, atoms::Vector{<:AbstractAtom})
    write(f, "begin projections\n")
    uniats = unique(atoms)
    projs = projections.(uniats)
    # projs = projections.(unique(atoms(structure)))
    if all(isempty.(projs))
        write(f, "random\n")
    else
        for (at, prjs) in zip(uniats, projs)
            if isempty(prjs)
                continue
            end
            write(f, projections_string(at))
            write(f, "\n")
        end
    end
    write(f, "end projections\n")
end

"""
    save(input::DFInput{Wannier90}, structure, filename::String=inpath(input))

Writes the `DFInput{Wannier90}` and `structure` to a file, that can be interpreted by WANNIER90.
The atoms in the structure must have projections defined.
"""
function save(input::DFInput{Wannier90}, structure, filename::String=inpath(input))
    open(filename, "w") do f
        for (flag, value) in input.flags
            write_flag_line(f, flag, value)
        end
        write(f, "\n")

        if structure != nothing
            write(f,"begin unit_cell_cart\n")
            write_cell(f, ustrip.(uconvert.(Ang, cell(structure)')))
            write(f,"end unit_cell_cart\n")
            write(f, "\n")
        end
        wan_write_projections(f, atoms(structure))

        write(f, "\n")
        write(f, "begin atoms_frac\n")
        for at in atoms(structure)
            pos = round.(position_cryst(at), digits=5)
            write(f, "$(name(at))  $(pos[1]) $(pos[2]) $(pos[3])\n")
        end
        write(f, "end atoms_frac\n")
        write(f, "\n")

        for block in input.data
            write(f, "begin $(block.name)\n")
            if block.name == :kpoint_path
                for i = 1:2:length(block.data)
                    letter1, k_points1 = block.data[i]
                    letter2, k_points2 = block.data[i+1]
                    write(f, "$letter1 $(k_points1[1]) $(k_points1[2]) $(k_points1[3]) $letter2 $(k_points2[1]) $(k_points2[2]) $(k_points2[3])\n")
                end

            elseif block.name == :kpoints
                for k in block.data
                    write(f, "$(k[1]) $(k[2]) $(k[3])\n")
                end
            end
            write(f, "end $(block.name)\n\n")
        end
    end
end

"""
    wan_read_output(filename)

Reads an outputfile for wannier.
Parsed info:
    :disentanglement,
    :wannierise,
    :final_state,
"""
function wan_read_output(filename)
    DisTuple = NamedTuple{(:Iter, :Ω_i_1, :Ω_i, :δ, :Time), Tuple{Int, Float64, Float64, Float64, Float64}}
    WanTuple = NamedTuple{(:Iter, :δ_spread, :∇RMS, :Spread, :Time), Tuple{Int, Float64, Float64, Float64, Float64}}
    WfcTuple = NamedTuple{(:index, :center, :Spread), Tuple{Int, Point3{Float64}, Float64}}
    data = Dict{Symbol, Any}()
    open(filename, "r") do f
        sreadline() = readline(f) |> strip
        line = sreadline()
        while !eof(f)
            if line == "Extraction of optimally-connected subspace"
                data[:disentanglement] = DisTuple[]
                for i = 1:4 readline(f) end
                line = sreadline()
                while !isempty(line)
                    push!(data[:disentanglement], parse(DisTuple, split(line)[1:end-2]))
                    line = sreadline()
                end
            elseif line == "Initial State"
                data[:wannierise] = WanTuple[]
                line = sreadline()
                while !occursin("Final State", line)
                    if occursin("CONV", line)
                        push!(data[:wannierise], parse(WanTuple, split(line)[1:end-2]))
                    end
                    line = sreadline()
                end
                data[:final_state] = WfcTuple[]
                line = sreadline()
                while !occursin("Sum", line)
                    line = replace(replace(replace(line, "(" => ""), ")" => ""), "," => "")
                    push!(data[:final_state], parse(WfcTuple, split(line)[5:end]))
                    line = sreadline()
                end
            end
            line = sreadline()
        end
    end
    return data
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Tiago Os��io
#* $Id$
#*-------------------------------------------------------------------
type StepSignal
	StepSignal()=begin
		new(
			positive(Dict{Symbol,Any}(
				:Unit=>"s"
			)),
			control_signal(),
			control_signal(),
			control_signal(Dict{Symbol,Any}(
				:PosX=>1,
				:PosY=>0.5
			)),
			[
				:(OutSignal = StartValue),
				:(OutSignal = FinalValue),
			],
			[
				"","",
			],
			[:StepTime,:StartValue,:FinalValue,],
			[:OutSignal,]
		)
	end
	StepTime::positive
	StartValue::control_signal
	FinalValue::control_signal
	OutSignal::control_signal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export StepSignal
function setEquationFlow(in::StepSignal)
	if(time < StepTime)
		addEquation(1)
	else
		addEquation(2)
	end
end
function atributes(in::StepSignal,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/StepSignal"
	drive!(fields,_)
	return fields
end
StepSignal(_::Dict{Symbol,Any})=begin
	newModel=StepSignal()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
"""
    get_EUR_to_USD(region::String)
convert Euro to US dollars
introduced because the clusters generated by the python script are in EUR for GER
"""
function get_EUR_to_USD(region::String)
   if region =="GER"
     ret = 1.109729
   else
     ret =1
   end
   return ret
end

"""
    sort_centers(centers::Array,weights::Array)
- centers: hours x days e.g.[24x9]
- weights: days [e.g. 9], unsorted
sorts the centers by weights from largest to smallest
"""
function sort_centers(centers::Array,
                      weights::Array
                      )
  i_w = sortperm(-weights)   # large to small (-)
  weights_sorted = weights[i_w]
  centers_sorted = centers[:,i_w]
  return centers_sorted, weights_sorted
end # function

"""
    z_normalize(data::ClustData;scope="full")
scope: "full", "sequence", "hourly"
"""
function z_normalize(data::ClustData;
                    scope="full"
                    )
 data_norm = Dict{String,Array}()
 mean= Dict{String,Array}()
 sdv= Dict{String,Array}()
 for (k,v) in data.data
   data_norm[k],mean[k],sdv[k] = z_normalize(v,scope=scope)
 end
 return ClustData(data.region,data.years,data.K,data.T,data_norm,data.weights,data.deltas,data.k_ids;mean=mean,sdv=sdv)
end

"""
     z_normalize(data::Array;scope="full")
z-normalize data with mean and sdv by hour
data: input format: (1st dimension: 24 hours, 2nd dimension: # of days)
scope: "full": one mean and sdv for the full data set; "hourly": univariate scaling: each hour is scaled seperately; "sequence": sequence based scaling
"""
function z_normalize(data::Array;
                    scope="full"
                    )
  if scope == "sequence"
    seq_mean = zeros(size(data)[2])
    seq_sdv = zeros(size(data)[2])
    data_norm = zeros(size(data))
    for i=1:size(data)[2]
      seq_mean[i] = mean(data[:,i])
      seq_sdv[i] = std(data[:,i])
      isnan(seq_sdv[i]) &&  (seq_sdv[i] =1)
        data_norm[:,i] = data[:,i] .- seq_mean[i]
      # handle edge case sdv=0
      if seq_sdv[i]!=0
        data_norm[:,i] = data_norm[:,i]./seq_sdv[i]
      end
    end
    return data_norm,seq_mean,seq_sdv
  elseif scope == "hourly"
    hourly_mean = zeros(size(data)[1])
    hourly_sdv = zeros(size(data)[1])
    data_norm = zeros(size(data))
    for i=1:size(data)[1]
      hourly_mean[i] = mean(data[i,:])
      hourly_sdv[i] = std(data[i,:])
      data_norm[i,:] = data[i,:] .- hourly_mean[i]
      # handle edge case sdv=0
      if hourly_sdv[i] !=0
        data_norm[i,:] = data_norm[i,:]./hourly_sdv[i]
      end
    end
    return data_norm, hourly_mean, hourly_sdv
  elseif scope == "full"
    hourly_mean = mean(data)*ones(size(data)[1])
    hourly_sdv = std(data)*ones(size(data)[1])
    # handle edge case sdv=0
    if hourly_sdv[1] != 0
      data_norm = (data.-hourly_mean[1])/hourly_sdv[1]
    else
      data_norm = (data.-hourly_mean[1])
    end
    return data_norm, hourly_mean, hourly_sdv #TODO change the output here to an immutable struct with three fields - use struct - "composite type"
  else
    @error("scope _ ",scope," _ not defined.")
  end
end # function z_normalize

"""
    undo_z_normalize(data_norm_merged::Array,mn::Dict{String,Array},sdv::Dict{String,Array};idx=[])
provide idx should usually be done as default within function call in order to enable sequence-based normalization, even though optional.
"""
function undo_z_normalize(data_norm_merged::Array,mn::Dict{String,Array},sdv::Dict{String,Array};idx=[])
  T = div(size(data_norm_merged)[1],length(keys(mn))) # number of time steps in one period. div() is integer division like in c++, yields integer (instead of float as in normal division)
  0 != rem(size(data_norm_merged)[1],length(keys(mn))) && @error("dimension mismatch") # rem() checks the remainder. If not zero, throw error.
  data_merged = zeros(size(data_norm_merged))
  i=0
  for (attr,mn_a) in mn
    i+=1
    data_merged[(1+T*(i-1)):(T*i),:]=undo_z_normalize(data_norm_merged[(1+T*(i-1)):(T*i),:],mn_a,sdv[attr];idx=idx)
  end
  return data_merged
end

"""
    undo_z_normalize(data_norm, mn, sdv; idx=[])
undo z-normalization data with mean and sdv by hour
normalized data: input format: (1st dimension: 24 hours, 2nd dimension: # of days)
hourly_mean ; 24 hour vector with hourly means
hourly_sdv; 24 hour vector with hourly standard deviations
"""
function undo_z_normalize(data_norm::Array, mn::Array, sdv::Array; idx=[])
  if size(data_norm,1) == size(mn,1) # hourly and full- even if idx is provided, doesn't matter if it is hourly
    data = data_norm .* sdv + mn * ones(size(data_norm)[2])'
    return data
  elseif !isempty(idx) && size(data_norm,2) == maximum(idx) # sequence based
    # we obtain mean and sdv for each day, but need mean and sdv for each centroid - take average mean and sdv for each cluster
    summed_mean = zeros(size(data_norm,2))
    summed_sdv = zeros(size(data_norm,2))
    for k=1:size(data_norm,2)
      mn_temp = mn[idx.==k]
      sdv_temp = sdv[idx.==k]
      summed_mean[k] = sum(mn_temp)/length(mn_temp)
      summed_sdv[k] = sum(sdv_temp)/length(sdv_temp)
    end
    data = data_norm * Diagonal(summed_sdv) +  ones(size(data_norm,1)) * summed_mean'
    return data
  elseif isempty(idx)
    @error("no idx provided in undo_z_normalize")
  end
end

"""
     sakoe_chiba_band(r::Int,l::Int)
calculates the minimum and maximum allowed indices for a lxl windowed matrix
for the sakoe chiba band (see Sakoe Chiba, 1978).
Input: radius r, such that |i(k)-j(k)| <= r
length l: dimension 2 of the matrix
"""
function sakoe_chiba_band(r::Int,l::Int)
  i2min = Int[]
  i2max = Int[]
  for i=1:l
    push!(i2min,max(1,i-r))
    push!(i2max,min(l,i+r))
  end
  return i2min, i2max
end

"""
     calc_SSE(data::Array,centers::Array,assignments::Array)
calculates Sum of Squared Errors between cluster representations and the data
"""
function calc_SSE(data::Array,centers::Array,assignments::Array)
  k=size(centers,2) # number of clusters
  n_periods =size(data,2)
  SSE_sum = zeros(k)
  for i=1:n_periods
    SSE_sum[assignments[i]] += sqeuclidean(data[:,i],centers[:,assignments[i]])
  end
  return sum(SSE_sum)
end # calc_SSE

"""
    calc_SSE(data::Array,centers::Array,assignments::Array)
calculates Sum of Squared Errors between cluster representations and the data
"""
function calc_SSE(data::Array,assignments::Array)
  centers=calc_centroids(data, assignments)
  k=size(centers,2) # number of clusters
  n_periods =size(data,2)
  SSE_sum = zeros(k)
  for i=1:n_periods
    SSE_sum[assignments[i]] += sqeuclidean(data[:,i],centers[:,assignments[i]])
  end
  return sum(SSE_sum)
end # calc_SSE

"""
    calc_centroids(data::Array,assignments::Array)
Given the data and cluster assignments, this function finds
the centroid of the respective clusters.
"""
function calc_centroids(data::Array,assignments::Array)
  K=maximum(assignments) #number of clusters
  n_per_period=size(data,1)
  n_periods =size(data,2)
  centroids=zeros(n_per_period,K)
  for k=1:K
    centroids[:,k]=sum(data[:,findall(assignments.==k)];dims=2)/length(findall(assignments.==k))
  end
  return centroids
end

"""
     calc_medoids(data::Array,assignments::Array)
Given the data and cluster assignments, this function finds
the medoids that are closest to the cluster center.
"""
function calc_medoids(data::Array,assignments::Array)
  K=maximum(assignments) #number of clusters
  n_per_period=size(data,1)
  n_periods =size(data,2)
  SSE=Float64[]
  for i=1:K
    push!(SSE,Inf)
  end
  centroids=calc_centroids(data,assignments)
  medoids=zeros(n_per_period,K)
  # iterate through all data points
  for i=1:n_periods
    d = sqeuclidean(data[:,i],centroids[:,assignments[i]])
    if d < SSE[assignments[i]] # if this data point is closer to centroid than the previously visited ones, then make this the medoid
      medoids[:,assignments[i]] = data[:,i]
      SSE[assignments[i]]=d
    end
  end
  return medoids
end

"""
     resize_medoids(data::Array,centers::Array,weights::Array,assignments::Array)
Takes in centers (typically medoids) and normalizes them such that for all clusters the average of the cluster is the same as the average of the respective original data that belongs to that cluster.
In order to use this method of the resize function, add assignments to the function call (e.g. clustids[5,1]).
"""
function resize_medoids(data::Array,centers::Array,weights::Array,assignments::Array)
    new_centers = zeros(centers)
    for k=1:size(centers)[2] # number of clusters
       is_in_k = assignments.==k
       n = sum(is_in_k)
       new_centers[:,k]=resize_medoids(reshape(data[:,is_in_k],:,n),reshape(centers[:,k] , : ,1),[1.0])# reshape is used for the side case with only one vector, so that resulting vector is 24x1 instead of 24-element
    end
    return new_centers
end

"""
    resize_medoids(data::Array,centers::Array,weights::Array)
This is the DEFAULT resize medoids function
Takes in centers (typically medoids) and normalizes them such that the yearly average of the clustered data is the same as the yearly average of the original data.
"""
function resize_medoids(data::Array,centers::Array,weights::Array)
    mu_data = sum(data)
    mu_clust = 0
    w_tot=sum(weights)
    for k=1:size(centers)[2]
      mu_clust += weights[k]/w_tot*sum(centers[:,k]) # weights[k]>=1
    end
    mu_clust *= size(data)[2]
    mu_data_mu_clust = mu_data/mu_clust
    new_centers = centers* mu_data_mu_clust
    return new_centers
end

"""
     resize_medoids(data::Array,centers::Array,weights::Array)
This is the DEFAULT resize medoids function
Takes in centers (typically medoids) and normalizes them such that the yearly average of the clustered data is the same as the yearly average of the original data.
"""
function resize_medoids(data::ClustData,centers::Array,weights::Array)
    (data.T * length(keys(data.data)) != size(centers,1) ) && @error("dimension missmatch between full input data and centers")
    centers_res = zeros(size(centers))
    # go through the attributes within data
    i=0
    for (k,v) in data.data
      i+=1
      # calculate resized centers for each attribute
      centers_res[(1+data.T*(i-1)):(data.T*i),:] = resize_medoids(v,centers[(1+data.T*(i-1)):(data.T*i),:],weights)
    end
    return centers_res
end

"""
    calc_weights(clustids::Array{Int}, n_clust::Int)
Calculates weights for clusters, based on clustids that are assigned to a certain cluster. The weights are absolute:    weights[i]>=1
"""
function calc_weights(clustids::Array{Int}, n_clust::Int)
    weights = zeros(n_clust)
    for j=1:length(clustids)
        weights[clustids[j]] +=1
    end
    return weights
end

"""
    set_clust_config(;kwargs...)
Add kwargs to a new Dictionary with the variables as entries
"""
function set_clust_config(;kwargs...)
  #Create new Dictionary
  config=Dict{String,Any}()
  # Loop through the kwargs and write them into Dictionary
  for kwarg in kwargs
    config[String(kwarg[1])]=kwarg[2]
  end
  # Return Directory with the information of kwargs
  return config
end


text:
#
# Correctness Tests
#

fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
quiet = length(ARGS) > 0 && ARGS[1] == "-q"
anyerrors = false

using Base.Test
using DataFrames

my_tests = ["utils.jl",
            "data.jl",
            "index.jl",
            "dataframe.jl",
            "io.jl",
            # "formula.jl",
            "constructors.jl",
            # "indexing.jl",
            "RDA.jl",
            "sort.jl",
            "grouping.jl",
            "join.jl",
            "vcat.jl",
            "iteration.jl",
            "duplicates.jl",
            "show.jl"]

println("Running tests:")

for my_test in my_tests
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch e
        anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        if fatalerrors
            rethrow(e)
        elseif !quiet
            showerror(STDOUT, e, backtrace())
            println()
        end
    end
end

if anyerrors
    throw("Tests failed")
end

stdin = Pkg.dir("DataFrames", "test", "stdin.sh")
run(`bash $stdin`)


text:
import Pkg
import Pkg: RegistrySpec, PackageSpec

## ------------------------------------------------------------------
function _info(msg; kwargs...)
    println()
    msg = rpad(string(msg, " "), 60, "-")
    @info(msg, kwargs...)
    println()
    println()
end

## ------------------------------------------------------------------
# desable autoprecompilation
ENV["JULIA_PKG_PRECOMPILE_AUTO"]=0

## ------------------------------------------------------------------
_info("Add registries")
Pkg.Registry.add(RegistrySpec(url = "https://github.com/josePereiro/CSC_Registry.jl"))

# ----------------------------------------------------------------------------
_info("Install fundamentals")
Pkg.add(name="ProjAssistant", version="0.2.1")
Pkg.add(name="ArgParse")

# ----------------------------------------------------------------------------
# ARGS
using ArgParse

set = ArgParseSettings()
@add_arg_table! set begin
    "--git-cmd"
        help = "the git command to use"
        default = "git"
end

if isinteractive()
    # Dev values
    GIT_CMD = "git"
else
    parsed_args = parse_args(set)
    GIT_CMD = parsed_args["git-cmd"]
end

## ------------------------------------------------------------------
# activate
using ProjAssistant
@quickactivate

## ------------------------------------------------------------------
_info("instantiate Chemostat_EColi")
Pkg.instantiate()

## ------------------------------------------------------------------
import Chemostat_EColi
const ChE = Chemostat_EColi

## ------------------------------------------------------------------
# globals
PROJ_ROOT = dirname(projdir(ChE))
cd(PROJ_ROOT)

_info("Globals"; PROJ_ROOT, GIT_CMD)

## ------------------------------------------------------------------
# git utils
function _checkout_tag(tag)
    cmd_strs = [GIT_CMD, "checkout", tag]
    run(Cmd(cmd_strs), wait=true)
end
function _clone_tag(outname, url, tag)
    cmd_strs = [ GIT_CMD, "clone", "--depth", "1", "--branch", tag, "--single-branch", url, outname]
    run(Cmd(cmd_strs), wait=true)
end

## ------------------------------------------------------------------
# Clone and dev dependencies
# TODO: (Chemostat_InSilico) pkg> add BenchmarkTools
let
    _info("Adding .jl devs")
    for (pkgname, url, tag) in [
            (
                "Chemostat_InSilico",
                "https://github.com/josePereiro/Chemostat_InSilico.jl",
                "v0.2.3"
            ),
            (
                "Chemostat_Kayser2005",
                "https://github.com/josePereiro/Chemostat_Kayser2005.jl",
                "v0.2.2"
            ),
            (
                "Chemostat_Nanchen2006",
                "https://github.com/josePereiro/Chemostat_Nanchen2006.jl",
                "v0.2.2"
            ),
            (
                "Chemostat_Folsom2014",
                "https://github.com/josePereiro/Chemostat_Folsom2014.jl",
                "v0.2.2"
            ),
        ]

        pkg_path = joinpath(PROJ_ROOT, pkgname)

        _info("clone and checkout"; pkgname, tag)
        if !isdir(pkg_path)
            cd(PROJ_ROOT)
            _clone_tag(pkgname, url, tag)
        end

        cd(pkg_path)
        _checkout_tag(tag)

        # develop
        _info("develop"; pkgname, tag)
        Pkg.activate(projdir(ChE))
        Pkg.develop(;path=pkg_path)

        _info("instanciate"; pkgname, tag)
        Pkg.activate(pkg_path)
        Pkg.instantiate()

        _info("precompile"; pkgname, tag)
        Pkg.precompile()

        cd(PROJ_ROOT)
    end
end

## ------------------------------------------------------------------
# precompile
_info("precompile Chemostat_EColi")
Pkg.activate(projdir(ChE))
Pkg.precompile()

## ------------------------------------------------------------------
# Test import
_info("Testing import")
try
    import Chemostat_Kayser2005
    import Chemostat_Nanchen2006
    import Chemostat_Folsom2014
    import Chemostat_InSilico
catch err
    @warn("Error loading source package, see scripts/0.1_install.jl")
    rethrow(err)
end

# ## ------------------------------------------------------------------
# # TODO: make repo public
# # clone paper repo
# _info("cloning tex project")
# let
#     pkgname = "MaxEnt_EColi_paper"
#     pkg_path = joipath(PROJ_ROOT, pkgname)
#     tag = ""
#     url="https://github.com/josePereiro/MaxEnt_EColi_paper"
#
#     if !isdir(pkg_path)
#         cd(PROJ_ROOT)
#         _clone_tag(pkgname, url, tag)
#     end
#     cd(pkg_path)
#     _checkout_tag(tag)
# end

_info("All good")


text:
# Specialize broadcasting of all of the non-modifying interfaces
import Base.broadcasted

@inline broadcasted(norm::T, l, m, x) where {T<:AbstractLegendreNorm} =
    broadcasted(legendre, norm, l, m, x)

# N.B. The broadcasting functions are near replicas of legendre() functions in
#      calculation.jl. Done so that bounds checking can occur before any array allocations.

function broadcasted(::typeof(legendre),
        norm::AbstractLegendreNorm, l::Integer, m::Integer, x)
    _chkdomain(l, m)
    boundscheck_hook(norm, l, m)
    z = Broadcast.materialize(x)
    Λ = _similar(z)
    unsafe_legendre!(norm, Λ, l, m, z)
    return (ndims(Λ) == 0) ? Λ[] : Λ
end

function broadcasted(::typeof(legendre), norm::AbstractLegendreNorm, l::DimOrInd, m::DimOrInd, x)
    if l isa AbstractUnitRange
        first(l) == 0 || throw(ArgumentError("Range of degrees l must start at 0"))
    end
    if m isa AbstractUnitRange
        if !(l isa AbstractUnitRange)
            throw(ArgumentError("Range of orders m requires range of degrees l"))
        end
        first(m) == 0 || throw(ArgumentError("Range of orders m must start at 0"))
    end
    lmax = l isa AbstractUnitRange ? last(l) : l
    mmax = m isa AbstractUnitRange ? last(m) : m
    _chkdomain(lmax, mmax)
    boundscheck_hook(norm, lmax, mmax)

    z = Broadcast.materialize(x)
    Λ = zeros(eltype(z), axes(z)..., size(l)..., size(m)...)
    unsafe_legendre!(norm, Λ, lmax, mmax, z)
    return Λ
end


text:
struct Features
    """
    A type for holding all the fixed quantities during a simulation. This involves
    simulation constants for detection probabiliites, constraints for consistency
    with data and some data based measures like the local cases population size,
    duration of sims and the length of the observation perod.
    """
    # constraints for the simulation
    max_forecast_cases::Int
    N::Int
    T_observed::Int
    T_end::Int
    # number of days for each consistency check
    τ::Int
    min_cases::Vector{Int}
    max_cases::Vector{Int}
    idxs_limits::Vector{UnitRange{Int64}}
    state::String
    # truncations for the various bits and pieces
    omicron_start_day::Int
    omicron_only_day::Int
end

mutable struct Realisation
    """
    A type for holding the realisations of the simulations. This is a cleaner way of
    holding the information for the three different matrices used. Z is for infections
    D is for observed cases and U is for undetected cases.
    """
    Z::Array{Int}
    Z_historical::Array{Int}
    D::Array{Int}
    U::Array{Int}

    function Realisation(
        sim_duration,
    )
        """
        Initialising a state array object with arrays of zeros.
        We pad the infection array Z with 35 days to account for
        infections occuring prior to the simulation period. We do this
        separately (and not in a struct) as the arrays are large.
        """

        Z = zeros(Int, sim_duration + 30, 3)
        Z_historical = zeros(Int, sim_duration + 30)
        D = zeros(Int, sim_duration, 3)
        U = zeros(Int, sim_duration, 3)

        return new(Z, Z_historical, D, U)
    end
end

mutable struct Results
    """
    A type for holding the realisations of the simulations. This is a cleaner way of
    holding the information for the three different matrices used. Z is for infections
    D is for observed cases and U is for undetected cases.
    """
    Z::Array{Int}
    Z_historical::Array{Int}
    D::Array{Int}
    U::Array{Int}

    function Results(
        sim_duration,
        nsims,
    )
        """
        Initialising a state array object with arrays of zeros.
        We pad the infection array Z with 35 days to account for
        infections occuring prior to the simulation period. We do this
        separately (and not in a struct) as the arrays are large.
        """

        Z = zeros(Int, sim_duration + 30, 3, nsims)
        Z_historical = zeros(Int, sim_duration + 30, nsims)
        D = zeros(Int, sim_duration, 3, nsims)
        U = zeros(Int, sim_duration, 3, nsims)

        return new(Z, Z_historical, D, U)
    end
end

struct Constants{S, T}
    """
    A type for the dynamical / simulation constants
    """
    p_detect::S
    p_symp::S
    p_detect_given_symp::S
    p_detect_given_asymp::S
    p_symp_given_detect::S
    α_s::S
    α_a::S
    p_detect_import::S
    prior_alpha::T
    prior_beta::T
    ϕ::T
    consistency_multiplier::T

    function Constants(start_date, end_date)
        """
        Constructor for the older version of the model. This assumes a simpler form of the CAR and
        doesn't have jurisdiction dependent assumptions in relation to the CAR.
        """

        # date vectors to indicate which period we're in. Noting that this is +5 days (mean inc
        # period) compared to the ranges used in the fitting. This is to deal with detection being
        # related to detection of cases and not infection dates.
        CAR_normal_before = Dates.Date(start_date):Dates.Day(1):Dates.Date("2021-12-09")
        CAR_low = Dates.Date("2021-12-10"):Dates.Day(1):Dates.Date(end_date)
        # all dates
        CAR_dates = [
            CAR_normal_before;
            CAR_low
        ]

        # initialise the
        CAR = 0.75 * ones(Float64, length(CAR_dates))
        # idx = CAR_dates .>= Dates.Date("2021-12-10")
        idx = CAR_dates .>= CAR_low[begin]
        CAR[idx] .= 0.5

        p_symp = zeros(Float64, length(CAR_dates))
        p_detect_given_asymp = zeros(Float64, length(CAR_dates))
        p_detect_given_symp = zeros(Float64, length(CAR_dates))
        p_symp_given_detect = zeros(Float64, length(CAR_dates))
        p_detect_import = zeros(Float64, length(CAR_dates))

        for (i, d) in enumerate(CAR_dates)
            p_detect = CAR[i]

            if d <= Dates.Date("2021-11-01")
                p_symp[i] = 0.7
            else
                p_symp[i] = 0.4
            end

            if d <= Dates.Date("2021-11-01")
                # p_symp = 0.7 and p_d = 0.75
                p_detect_given_symp[i] = 0.9375
                p_detect_given_asymp[i] = 0.3125
                p_detect_import[i] = 0.98
            elseif d <= Dates.Date("2021-12-10")
                # p_symp = 0.4 and p_d = 0.75
                p_detect_given_symp[i] = 0.9
                p_detect_given_asymp[i] = 0.65
                p_detect_import[i] = p_detect
            else
                # p_symp = 0.4 and p_d = 0.5 noting that the line defined by
                # p(d|s) = -p(a)/p(s) p(d|a) + p(d)/p(s)
                # translates vertically when the CA changes and so we just find the point on
                # the reduced CA line which is closest to the previous values
                # (p(d|a), p(d|s)) = (0.65, 0.9)
                p_detect_given_asymp[i] = 0.3615
                p_detect_given_symp[i] = 0.7076
                p_detect_import[i] = p_detect
            end

            p_symp_given_detect[i] = p_detect_given_symp[i] * p_symp[i] / p_detect

        end

        p_detect = CAR

        γ = 0.5     # relative infectiousness of asymptomatic
        α_s = 1 ./ (p_symp .+ γ * (1 .- p_symp))    # contribution of TP attributable to symptomatic
        α_a = γ * α_s   # contribution of TP attributable to asymptomatic

        # for checking consistency when we miss an outbreak
        consistency_multiplier = 5.0
        # prior parametes for the import model
        prior_alpha = 0.5
        prior_beta = 0.1
        # ema smoothing factor for the import model
        ϕ = 0.5

        # use this to ensure the Constants type is correct
        S = typeof(p_symp)
        T = typeof(prior_alpha)

        return new{S, T}(
            p_detect,
            p_symp,
            p_detect_given_symp,
            p_detect_given_asymp,
            p_symp_given_detect,
            α_s,
            α_a,
            p_detect_import,
            prior_alpha,
            prior_beta,
            ϕ,
            consistency_multiplier,
        )
    end

    function Constants(start_date, end_date, state)
        """
        Newer version constructor. This is distinguished from the other by the extra parameter
        `state` which enables us to produce jurisdiction specific scenarios.
        """
        # date vectors to indicate which period we're in
        CAR_normal_before = Dates.Date(start_date):Dates.Day(1):Dates.Date("2021-12-07")
        CAR_low = Dates.Date("2021-12-08"):Dates.Day(1):Dates.Date("2022-01-17")
        CAR_normal_after = Dates.Date("2022-01-18"):Dates.Day(1):Dates.Date(end_date)
        # all dates
        CAR_dates = [
            CAR_normal_before;
            CAR_low;
            CAR_normal_after
        ]

        # initialise the
        CAR = 0.75 * ones(Float64, length(CAR_dates))

        # if we are one of the jurisdictions that got overwhelmed, apply time varying CAR
        if state in ("NSW", "VIC", "ACT", "QLD")
            idx = (CAR_dates .>= CAR_low[begin]) .& (CAR_dates .<= CAR_low[end])
            CAR[idx] .= 0.333
            idx = CAR_dates .> CAR_low[end]
            CAR[idx] .= 0.75
        end

        p_symp = zeros(Float64, length(CAR_dates))
        p_detect_given_asymp = zeros(Float64, length(CAR_dates))
        p_detect_given_symp = zeros(Float64, length(CAR_dates))
        p_symp_given_detect = zeros(Float64, length(CAR_dates))
        p_detect_import = zeros(Float64, length(CAR_dates))

        for (i, d) in enumerate(CAR_dates)
            p_detect = CAR[i]

            if d <= Dates.Date("2021-11-01")
                p_symp[i] = 0.7
            else
                p_symp[i] = 0.4
            end

            if d <= Dates.Date("2021-11-01")
                # p_symp = 0.7 and p_d = 0.75
                p_detect_given_symp[i] = 0.9375
                p_detect_given_asymp[i] = 0.3125
                p_detect_import[i] = 0.98
            elseif d <= Dates.Date("2021-12-07") || d >= Dates.Date("2022-01-18")
                # p_symp = 0.4 and p_d = 0.75
                p_detect_given_symp[i] = 0.9
                p_detect_given_asymp[i] = 0.65
                p_detect_import[i] = p_detect
            else
                # p_symp = 0.4 and p_d = 0.333
                p_detect_given_asymp[i] = 0.16885
                p_detect_given_symp[i] = 0.57923
                p_detect_import[i] = p_detect
            end

            p_symp_given_detect[i] = p_detect_given_symp[i] * p_symp[i] / p_detect

        end

        p_detect = CAR

        γ = 0.5     # relative infectiousness of asymptomatic
        α_s = 1 ./ (p_symp .+ γ * (1 .- p_symp))
        α_a = γ * α_s

        # other parameters
        consistency_multiplier = 5.0
        # prior parametes for the import model
        prior_alpha = 0.5
        prior_beta = 0.1
        # ema smoothing factor
        ϕ = 0.5

        S = typeof(p_symp)
        T = typeof(prior_alpha)

        return new{S, T}(
            p_detect,
            p_symp,
            p_detect_given_symp,
            p_detect_given_asymp,
            p_symp_given_detect,
            α_s,
            α_a,
            p_detect_import,
            γ,
            prior_alpha,
            prior_beta,
            ϕ,
            consistency_multiplier,
        )
    end
end

struct IndividualTypeMap
    S::Int
    A::Int
    I::Int

    function IndividualTypeMap()
        """
        Constructor for the mapping between type symbols S, A, I and integers.
        """
        S = 1
        A = 2
        I = 3

        return new(S, A, I)

    end
end

struct Forecast
    """
    This type wraps all the information relating to the forecast in a single accessible
    object. This drastically simplifies the layout of the code at no performance hit.
    """
    sim_features::Features
    sim_realisation::Realisation
    sim_results::Results
    sim_constants::Constants
    individual_type_map::IndividualTypeMap
end

struct JurisdictionAssumptions
    """
    This type holds all the assumptions for each jurisdiction. This
    involves the assumed starting dates, population sizes and initial conditions. It also
    includes key dates we don't have stored elsewhere.
    """
    simulation_start_dates::Dict{String, String}
    pop_sizes::Dict{String, Int}
    initial_conditions::Dict{
        String,
        NamedTuple{(:S, :A, :I), Tuple{Int64, Int64, Int64}},
    }
    omicron_start_date::Date
    omicron_dominant_date::Date

    function JurisdictionAssumptions()
        simulation_start_dates = Dict{String, String}(
            "NSW" => "2021-06-01",
            "QLD" => "2021-11-20",
            "SA" => "2021-11-01",
            "TAS" => "2021-11-01",
            "WA" => "2021-12-15",
            "ACT" => "2021-12-01",
            "NT" => "2021-12-01",
            "VIC" => "2021-08-01",
        )

        # date we want to apply increase in cases due to Omicron
        # put a small delay on this as we use the infection dates to index things but in
        # the fitting we used the confirmation dates
        omicron_start_date = Dates.Date("2021-11-15")
        omicron_dominant_date = Dates.Date("2021-12-15") - Dates.Day(5)

        pop_sizes = Dict{String, Int}(
            "NSW" => 8189266,
            "QLD" => 5221170,
            "SA" => 1773243,
            "TAS" => 541479,
            "VIC" => 6649159,
            "WA" => 2681633,
            "ACT" => 432266,
            "NT" => 246338,
        )

        initial_conditions = Dict{
            String, NamedTuple{(:S, :A, :I), Tuple{Int64, Int64, Int64}}
        }(
            "NSW" => (S = 0, A = 0, I = 0),
            "QLD" => (S = 1, A = 1, I = 0),
            "SA" => (S = 0, A = 0, I = 0),
            "TAS" => (S = 0, A = 0, I = 0),
            "VIC" => (S = 5, A = 5, I = 0),
            "WA" => (S = 0, A = 0, I = 0),
            "ACT" => (S = 12, A = 10, I = 0),
            "NT" => (S = 0, A = 0, I = 0),
        )

        return new(
            simulation_start_dates,
            pop_sizes,
            initial_conditions,
            omicron_dominant_date,
        )

    end
end


text:
# cudnnDataType_t
const CUDNN_DATA_FLOAT = Cint(0)
const CUDNN_DATA_DOUBLE = Cint(1)
const CUDNN_DATA_HALF = Cint(2)
const CUDNN_DATA_INT8 = Cint(3)
const CUDNN_DATA_INT32 = Cint(4)
const CUDNN_DATA_INT8x4 = Cint(5)

const Cptr = Ptr{Cvoid}
datatype(::Type{Float32}) = CUDNN_DATA_FLOAT
datatype(::Type{Float64}) = CUDNN_DATA_DOUBLE
datatype(::Type{Float16}) = CUDNN_DATA_HALF
datatype(::Type{Int8}) = CUDNN_DATA_INT8
datatype(::Type{Int32}) = CUDNN_DATA_INT32

# cudnnNanPropagation_t
const CUDNN_NOT_PROPAGATE_NAN = Cint(0)
const CUDNN_PROPAGATE_NAN = Cint(1)

# cudnnMathType_t
const CUDNN_DEFAULT_MATH = Cint(0)
const CUDNN_TENSOR_OP_MATH = Cint(1)

# cudnnDeterminism_t
const CUDNN_NON_DETERMINISTIC = Cint(0)
const CUDNN_DETERMINISTIC = Cint(1)

# cudnnOpTensorOp_t
const CUDNN_OP_TENSOR_ADD = Cint(0)
const CUDNN_OP_TENSOR_MUL = Cint(1)
const CUDNN_OP_TENSOR_MIN = Cint(2)
const CUDNN_OP_TENSOR_MAX = Cint(3)
const CUDNN_OP_TENSOR_SQRT = Cint(4)
const CUDNN_OP_TENSOR_NOT = Cint(5)

# cudnnIndicesType_t
const CUDNN_32BIT_INDICES = Cint(0)
const CUDNN_64BIT_INDICES = Cint(1)
const CUDNN_16BIT_INDICES = Cint(2)
const CUDNN_8BIT_INDICES = Cint(3)

# cudnnPoolingMode_t
const CUDNN_POOLING_MAX = Cint(0)
const CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING = Cint(1)
const CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING = Cint(2)
const CUDNN_POOLING_MAX_DETERMINISTIC = Cint(3)

# cudnnLRNMode_t
const CUDNN_LRN_CROSS_CHANNEL_DIM1 = Cint(0)

# cudnnDivNormMode_t
const CUDNN_DIVNORM_PRECOMPUTED_MEANS = Cint(0)

# cudnnBatchNormMode_t
const CUDNN_BATCHNORM_PER_ACTIVATION = Cint(0)
const CUDNN_BATCHNORM_SPATIAL = Cint(1)
const CUDNN_BATCHNORM_SPATIAL_PERSISTENT = Cint(2)

# cudnnSamplerType_t
const CUDNN_SAMPLER_BILINEAR = Cint(0)

# cudnnCTCLossAlgo_t
const CUDNN_CTC_LOSS_ALGO_DETERMINISTIC = Cint(0)
const CUDNN_CTC_LOSS_ALGO_NON_DETERMINISTIC = Cint(1)


text:
################################################################################
#
#  Valuation
#
################################################################################

function valuation(a::UInt, b::UInt)
  return ccall((:n_remove, libflint), Int, (Ref{UInt}, UInt), a, b)
end

# CF:
# The idea is that valuations are mostly small, eg. in the class group
# algorithm. So this version computes the completion and the embedding into it
# at small precision and can thus compute (small) valuation at the effective
# cost of an mod(nmod_poly, nmod_poly) operation.
# Isn't it nice?
function val_func_no_index_small(p::NfOrdIdl)
  P = p.gen_one
  @assert P <= typemax(UInt)
  K = nf(order(p))
  Rx = PolynomialRing(GF(UInt(P), cached=false), cached=false)[1]
  Zx = PolynomialRing(FlintZZ, cached = false)[1]
  gR = Rx(p.gen_two.elem_in_nf)
  f = Rx(K.pol)
  gR = gcd!(gR, gR, f)
  g = lift(Zx, gR)
  k = flog(fmpz(typemax(UInt)), P)
  g = hensel_lift(Zx(K.pol), g, P, k)
  Sx = PolynomialRing(ResidueRing(FlintZZ, UInt(P)^k, cached=false), cached=false)[1]
  g = Sx(g)
  h = Sx()
  uP = UInt(P)
  local vfunc
  let h = h, g = g, P = P, uP = uP
    function vfunc(x::nf_elem, no::fmpq = fmpq(0))
      d = denominator(x)
      nf_elem_to_nmod_poly!(h, x, false) # ignores the denominator
      h = rem!(h, h, g)
      c = Nemo.coeff_raw(h, 0)
      v = c==0 ? typemax(Int) : valuation(c, uP)
      for i=1:degree(h)
        c = Nemo.coeff_raw(h, i)
        v = min(v, c==0 ? typemax(Int) : valuation(c, uP))
      end
      return v-valuation(d, P)::Int
    end
  end
  return vfunc
end

function val_func_no_index(p::NfOrdIdl)
  P = p.gen_one
  K = nf(order(p))
  Rx, g = PolynomialRing(GF(P, cached=false), cached=false)
  Zx = PolynomialRing(FlintZZ, cached = false)[1]
  nf_elem_to_gfp_fmpz_poly!(g, p.gen_two.elem_in_nf, false)
  f = Rx(K.pol)
  g = gcd(g, f)
  g = lift(Zx, g)
  g = hensel_lift(Zx(K.pol), g, P, 10)
  Sx = PolynomialRing(ResidueRing(FlintZZ, P^5, cached=false), cached=false)[1]
  g = Sx(g)
  h = Sx()
  c = fmpz()
  local vfunc
  let h = h, g = g, P = P
    function vfunc(x::nf_elem, no::fmpq = fmpq(0))
      d = denominator(x)
      nf_elem_to_fmpz_mod_poly!(h, x, false) # ignores the denominator
      h = rem!(h, h, g)
      _coeff_as_fmpz!(c, h, 0)
      v = iszero(c) ? 100 : valuation(c, P)
      for i=1:degree(h)
        _coeff_as_fmpz!(c, h, i)
        v = min(v, iszero(c) ? 100 : valuation(c, P))
      end
      return v-valuation(d, P)::Int
    end
  end
  return vfunc
end

function _coeff_as_fmpz!(c::fmpz, f::fmpz_mod_poly, i::Int)
  ccall((:fmpz_mod_poly_get_coeff_fmpz, libflint), Nothing,
        (Ref{fmpz}, Ref{fmpz_mod_poly}, Int, Ref{fmpz_mod_ctx_struct}), c, f, i, f.parent.base_ring.ninv)
  return nothing
end

function val_func_index(p::NfOrdIdl)
  # We are in the index divisor case. In larger examples, a lot of
  # time is spent computing denominators of order elements.
  # By using the representation matrix to multiply, we can stay in the order
  # and still be fast (faster even than in field).
  O = order(p)
  pi = inv(p)
  P = p.gen_one
  pi_2 = mod(pi.num.gen_two.elem_in_nf, P^2)
  M = representation_matrix(O(pi_2, false))

  local val
  let P = P, O = O, M = M, p = p
    function val(x::nf_elem, no::fmpq = fmpq(0))
      v = 0
      d, x_mat = integral_split(x, O)
      Nemo.mul!(x_mat, x_mat, M)
      c = content(x_mat)
      vc = valuation(c, P)
      while vc > 0  # should divide and test in place
        #Carlo : How? Not clear how to improve this code.
	      divexact!(x_mat, x_mat, c)
        mul!(x_mat, x_mat, M)
        v += 1 + (vc-1)*p.splitting_type[1]
        c = content(x_mat)
        vc = valuation(c, P)
      end
      return v-Int(valuation(d, P))*p.splitting_type[1]
    end
  end
  return val
end

# CF:
# Classical algorithm of Cohen, but take a valuation element with smaller (?)
# coefficients. Core idea is that the valuation element is, originally, den*gen_two(p)^-1
# where gen_two(p) is "small". Actually, we don't care about gen_two, we
# need gen_two^-1 to be small, hence this version.
function val_fun_generic_small(p::NfOrdIdl)
  P = p.gen_one
  K = nf(order(p))
  O = order(p)
  e = anti_uniformizer(p)
  local val
  let e = e, P = P, p = p, O = O
    function val(x::nf_elem, no::fmpq = fmpq(0))
      nn = fmpz(0)
      v = 0
      p_mod = fmpz(0)
      d = denominator(x, O)
      x *= d
      if !iszero(no)
        nn = numerator(no*(d^degree(O)))
        vp = min(valuation(nn, norm(p)), 10)
        p_mod = P^(div(vp, ramification_index(p))+1)
      else
        p_mod = P^(div(10, ramification_index(p))+1)
      end
      x = mod(x, p_mod)
      mul!(x, x, e)
      while x in O && v < 10
        v += 1
        if !iszero(no)
          nn = divexact(nn, norm(p))
          if !divisible(nn, norm(p))
            break
          end
        end
        x = mod(x, p_mod)
        mul!(x, x, e)
      end
      if v == 10
        return 100
      end
      return v-valuation(d, P)*p.splitting_type[1]
    end
  end
  return val
end

function val_func_generic(p::NfOrdIdl)
  @hassert :NfOrd 1 isprime(p)
  P = p.gen_one
  K = nf(order(p))
  O = order(p)
  e = anti_uniformizer(p)
  local val
  let e = e, P = P, p = p, O = O
    function val(x::nf_elem, no::fmpq = fmpq(0))
      nn = fmpz(0)
      v = 0
      p_mod = fmpz(0)
      d = denominator(x, O)
      if !iszero(no)
        nn = numerator(no*(d^degree(O)))
        p_mod = P^(div(valuation(nn, norm(p)), ramification_index(p))+1)
	      x = mod(x, p_mod)
      end
      x *= d
      x = x*e
      while x in O
        v += 1
        if !iszero(no)
          nn = divexact(nn, norm(p))
          if !divisible(nn, norm(p))
            break
          end
          x = mod(x, p_mod)
        end
        mul!(x, x, e)
      end
      return v-valuation(d, P)*p.splitting_type[1]
    end
  end
  return val
end

function valuation_with_anti_uni(a::nf_elem, anti_uni::nf_elem, I::NfOrdIdl)
  O = order(I)
  b = a*anti_uni
  if !(b in O)
    return 0
  end
  v = 1
  mul!(b, b, anti_uni)
  while b in O
    v += 1
    mul!(b, b, anti_uni)
  end
  return v
end

function _isindex_divisor(O::NfOrd, P::NfOrdIdl)
  @assert isprime_known(P) && isprime(P)
  if !isone(denominator(P.gen_two.elem_in_nf))
    return true
  end
  R = GF(Int(minimum(P)), cached = false)
  Rt, t = PolynomialRing(R, "x", cached = false)
  f = Rt(nf(P).pol)
  g = Rt(P.gen_two.elem_in_nf)
  d = gcd(f, g)
  if !divides(f, d^2)[1] && isirreducible(d)
    return false
  else
    return true
  end
end

#Function that chooses the valuation depending on the properties of the ideal
function assure_valuation_function(p::NfOrdIdl)
  if isdefined(p, :valuation)
    return nothing
  end
  O = order(p)
  K = nf(O)
  # for generic ideals
  if p.splitting_type[2] == 0
    assure_2_normal(p)
    anti_uni = anti_uniformizer(p)
    local val2
    let O = O, p = p, anti_uni = anti_uni, K = K
      function val2(s::nf_elem, no::fmpq = fmpq(0))
        d = denominator(s, O)
        x = d*s
        if gcd(d, minimum(p, copy = false)) == 1
          return valuation_with_anti_uni(x, anti_uni, p)::Int
        else
          return valuation_with_anti_uni(x, anti_uni, p)::Int - valuation_with_anti_uni(K(d), anti_uni, p)::Int
        end
      end
    end
    p.valuation = val2
    return nothing
  end
  P = minimum(p)
  if degree(O) < 40 && p.splitting_type[1]*p.splitting_type[2] == degree(O)
    local val3
    let P = P, p = p
      function val3(s::nf_elem, no::fmpq = fmpq(0))
        return divexact(valuation(iszero(no) ? norm(s) : no, P)[1], p.splitting_type[2])::Int
      end
    end
    p.valuation = val3
  elseif mod(index(O), P) != 0 && ramification_index(p) == 1
    if fits(UInt, P^2)
      f1 = val_func_no_index_small(p)
      f2 = val_func_generic(p)
      local val1
      let f1 = f1, f2 = f2
        function val1(x::nf_elem, no::fmpq = fmpq(0))
          v = f1(x, no)
          if v > 100  # can happen ONLY if the precision in the .._small function
                      # was too small.
            return f2(x, no)::Int
          else
            return v::Int
          end
        end
      end
      p.valuation = val1
    else
      f8 = val_func_no_index(p)
      f9 = val_func_generic(p)
      local val1
      let f8 = f8, f9 = f9
        function val_large_non_index(x::nf_elem, no::fmpq = fmpq(0))
          v = f8(x, no)
          if v > 10  # can happen ONLY if the precision in the .._small function
                      # was too small.
            return f9(x, no)::Int
          else
            return v::Int
          end
        end
      end
      p.valuation = val_func_generic(p)
    end
  elseif ramification_index(p) == 1 && fits(UInt, P^2) && !_isindex_divisor(O, p)
    f3 = val_func_no_index_small(p)
    if degree(O) > 80
      f5 = val_func_generic(p)
      local val5
      let f3 = f3, f5 = f5
        function val5(x::nf_elem, no::fmpq = fmpq(0))
          v = f3(x, no)
          if v > 100  # can happen ONLY if the precision in the .._small function
                      # was too small.
            return f5(x, no)::Int
          else
            return v::Int
          end
        end
      end
      p.valuation = val5
    else
      f4 = val_func_index(p)
      local val4
      let f3 = f3, f4 = f4
        function val4(x::nf_elem, no::fmpq = fmpq(0))
          v = f3(x, no)
          if v > 100  # can happen ONLY if the precision in the .._small function
                      # was too small.
            return f4(x, no)::Int
          else
            return v::Int
          end
        end
      end
      p.valuation = val4
    end
  elseif degree(O) > 80
    f6 = val_func_generic(p)
    f7 = val_fun_generic_small(p)
    local val_gen
    let f7 = f7, f6 = f6
      function val_gen(x::nf_elem, no::fmpq = fmpq(0))
        vv = f7(x, no)
        if vv == 100
          return f6(x, no)
        else
          return vv
        end
      end
    end
    p.valuation = val_gen
  else
    p.valuation = val_func_index(p)
  end
  return nothing
end

function valuation(a::NfAbsNSElem, p::NfAbsOrdIdl, n::fmpq = fmpq(0))
  return valuation_naive(a, p)
end

function valuation(a::nf_elem, p::NfOrdIdl, no::fmpq = fmpq(0))
  if parent(a) !== nf(order(p))
    throw(error("Incompatible parents"))
  end
  if !isdefining_polynomial_nice(parent(a)) || order(p).ismaximal != 1
    return valuation_naive(a, p)::Int
  end
  @hassert :NfOrd 0 !iszero(a)
  assure_valuation_function(p)
  if p.is_prime != 1
    return Int(p.valuation(a, no))::Int
  end
  #First, check the content of a as a polynomial.
  #We remove the numerator of the content, as the
  #valuation for integers is much easier.
  O = order(p)
  K = nf(O)
  Zx = PolynomialRing(FlintZZ, "x")[1]
  pol_a = Zx(denominator(a)*a)
  c = content(pol_a)
  valnum = Int(valuation(c, p))
  b = divexact(a, c)

  nno = no
  if !iszero(nno)
    nno = divexact(nno, c^degree(K))
  end
  res = Int(p.valuation(b, nno))::Int
  res += valnum
  return res
end

@doc Markdown.doc"""
    valuation(a::nf_elem, p::NfOrdIdl) -> fmpz
    valuation(a::NfOrdElem, p::NfOrdIdl) -> fmpz
    valuation(a::fmpz, p::NfOrdIdl) -> fmpz

Computes the $\mathfrak p$-adic valuation of $a$, that is, the largest $i$
such that $a$ is contained in $\mathfrak p^i$.
"""
valuation(a::NfOrdElem, p::NfOrdIdl) = valuation(a.elem_in_nf, p)

@doc Markdown.doc"""
    valuation(a::nf_elem, p::NfOrdIdl) -> fmpz
    valuation(a::NfOrdElem, p::NfOrdIdl) -> fmpz
    valuation(a::fmpz, p::NfOrdIdl) -> fmpz

Computes the $\mathfrak p$-adic valuation of $a$, that is, the largest $i$
such that $a$ is contained in $\mathfrak p^i$.
"""
function valuation(a::fmpz, p::NfAbsOrdIdl)
  if p.splitting_type[1] == 0
    return valuation_naive(order(p)(a), p)
  end
  P = minimum(p, copy = false)
  return valuation(a, P)* p.splitting_type[1]
end
@doc Markdown.doc"""
    valuation(a::Integer, p::NfOrdIdl) -> fmpz
Computes the $\mathfrak p$-adic valuation of $a$, that is, the largest $i$
such that $a$ is contained in $\mathfrak p^i$.
"""
valuation(a::Integer, p::NfAbsOrdIdl) = valuation(fmpz(a), p)

#TODO: some more intelligence here...
function valuation_naive(A::NfAbsOrdIdl, B::NfAbsOrdIdl)
  @assert !isone(B)
  Bi = inv(B)
  i = 0
  C = simplify(A* Bi)
  while denominator(C) == 1
    C = simplify(Bi*C)
    i += 1
  end
  return i
end

#TODO: some more intelligence here...
#      in non-maximal orders, interesting ideals cannot be inverted
#      maybe this needs to be checked...
function valuation_naive(x::NfAbsOrdElem, B::NfAbsOrdIdl)
  @assert !isone(B)
  i = 0
  C = B
  while x in C
    i += 1
    C *= B
  end
  return i
end

function valuation_naive(x::T, B::NfAbsOrdIdl) where T <: Union{nf_elem, NfAbsNSElem}
  @assert !isone(B)
  i = 0
  C = B
  O = order(B)
  d = denominator(x, O)
  return valuation_naive(O(x*d), B) - valuation_naive(O(d), B)
end

@doc Markdown.doc"""
    valuation(A::NfOrdIdl, p::NfOrdIdl) -> fmpz

Computes the $\mathfrak p$-adic valuation of $A$, that is, the largest $i$
such that $A$ is contained in $\mathfrak p^i$.
"""
function valuation(A::NfAbsOrdIdl, p::NfAbsOrdIdl)
  if has_minimum(A) && has_minimum(p) && !divisible(minimum(A, copy = false), minimum(p, copy = false))
    return 0
  end
  if has_princ_gen_special(A)
    gen = princ_gen_special(A)
    return valuation(gen, p)
  end
  if A.is_principal == 1 && isdefined(A, :princ_gen)
    return valuation(A.princ_gen.elem_in_nf, p, fmpq(norm(A)))
  end
  _assure_weakly_normal_presentation(A)
  if !isdefined(p, :splitting_type) || p.splitting_type[1] == 0 #ie. if p is non-prime...
    return valuation_naive(A, p)
  end
  if iszero(A.gen_two)
    return valuation(A.gen_one, p)
  end
  v1 = valuation(A.gen_one, p)
  return min(v1, valuation(A.gen_two.elem_in_nf, p, fmpq(norm(A))))
end


text:
# use this to time regular mcmc(), but without plots.
using Distributions, Random
function main()
    y = rand(Exponential(3.0),30)
    # set prior, likelihood and proposal
    Prior = θ -> pdf.(Ref(LogNormal(1.0,1.0)), θ)
    lnL = θ -> sum(logpdf.(Ref(Exponential(θ)), y))
    tuning = 0.5
    Proposal = θ -> rand(LogNormal(log(θ),tuning))
    # get the chain, plot posterior, and descriptive stats
    chain = mcmc(1.0, 100000, 10000, Prior, lnL, Proposal, false) # start value, chain length, and burnin
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
# TLC = Three Letters Code

const TLC = Dict(
    "ACX" => :linear_acceleration_x,
    "ACY" => :linear_acceleration_y,
    "ACZ" => :linear_acceleration_z,
    "ANX" => :angular_acceleration_x,
    "ANY" => :angular_acceleration_y,
    "ABZ" => :angular_acceleration_z,  # ABZ? not ANZ?
    "ATS" => :altimeter_pressure,
    "BEI" => :gnss_beidou,
    "BFI" => :blind_flying_instrument,
    "CCL" => :competition_class,
    "CGD" => :change_geodetic_datum,
    "CID" => :competition_id,
    "CLB" => :club,
    "CM2" => :second_crew_member,
    "COT" => :controller_temperature,
    "CUR" => :electrical_current,
    "DAE" => :displacement_east,
    "DAN" => :displacement_north,
    "DB1" => :date_of_birth_pilot_in_charge,
    "DB2" => :date_of_birth_second_crew_member,
    "DTE" => :date,
    "DTM" => :geodetic_datum,
    "EGT" => :exhaust_gas_temperature,
    "FLE" => :fuel_level,
    "FFL" => :fuel_flow,
    "ENL" => :environmental_noise_level,
    "FIN" => :finish,
    "FLP" => :flap_position,
    "FRS" => :flight_recorder_security,
    "FTY" => :flight_recorder_type,
    "FXA" => :fix_accuracy,
    "GAL" => :gnss_galileo,
    "GID" => :glider_id,
    "GLO" => :gnss_glonass,
    "GPS" => :gnss_gps,
    "GSP" => :ground_speed,
    "GTY" => :glider_type,
    "HDM" => :heading_magnetic,
    "HDT" => :heading_true,
    "IAS" => :airspeed,
    "JPT" => :jet_pipe_temperature,
    "LEB" => :battery_left,
    "LAD" => :last_places_of_decimal_minutes_of_latitude,
    "LCU" => :data_seeyou_l,
    "LOD" => :last_places_of_decimal_minutes_of_longitude,
    "LOV" => :low_voltage,
    "MAC" => :mac_ready_setting,
    "MCU" => :data_seeyou_m,
    "MP2" => :additional_means_of_propulsion_2,
    "MP3" => :additional_means_of_propulsion_3,
    "MP4" => :additional_means_of_propulsion_4,
    "MP5" => :additional_means_of_propulsion_5,
    "MP6" => :additional_means_of_propulsion_6,
    "MP7" => :additional_means_of_propulsion_7,
    "MP8" => :additional_means_of_propulsion_8,
    "MP9" => :additional_means_of_propulsion_9,
    "MOT" => :motor_temperature,
    "MOP" => :means_of_propulsion,
    "OA1" => :position_other_aircraft_1,
    "OA2" => :position_other_aircraft_2,
    "OA3" => :position_other_aircraft_3,
    "OA4" => :position_other_aircraft_4,
    "OA5" => :position_other_aircraft_5,
    "OA6" => :position_other_aircraft_6,
    "OA7" => :position_other_aircraft_7,
    "OA8" => :position_other_aircraft_8,
    "OA9" => :position_other_aircraft_9,
    "OAT" => :outside_air_temperature,
    "ONT" => :on_task,
    "OOI" => :oo_id,
    "PEV" => :pilot_event,
    "PFC" => :post_flight_claim,
    "PLT" => :pilot_in_charge,
    "PRS" => :pressure_altitude_sensor,
    "RAI" => :receiver_autonomous_integrity_monitoring,
    "REX" => :record_addition,
    "RFW" => :firmware_version,  # of flight recorder
    "RHW" => :hardware_version,  # of flight recorder
    "RPM" => :revolution_per_minute,
    "SEC" => :security,
    "SIT" => :site,
    "SIU" => :satellites_in_use,
    "STA" => :start_event,
    "TAS" => :airspeed_true,
    "TDS" => :utc_time_decimal_seconds,
    "TEN" => :total_energy_altitude,
    "TPC" => :turn_point_confirmation,
    "TRT" => :true_task,
    "TZN" => :timezone_offset,
    "UND" => :undercarriage,  # landing gear
    "VAR" => :uncompensated_variometer,
    "VOL" => :electric_voltage,
    "VAT" => :compensated_variometer,
    "VXA" => :vertical_fix_accuracy,
    "WDI" => :wind_direction,
    "WSP" => :wind_speed,
)

#=
"XN*"
A manufacturer code where N is the manufacturer's single-character IGC name (para A3.5.6) and * is any character.
The manufacturer must specify its meaning and us in the documentation for the recorder and its use must be
approved by GFAC before IGC-approval. The X prefix is intended to allow a trial with a provisional new code
before deciding whether it is worthwhile adding to the full list.
=#

const OBSOLETE_TLC = Dict(
    "CCN" => :camera_connect,
    "CDC" => :camera_disconnect,
    "DOB" => :date_of_birth_pilot,  # now use DB1
    "PHO" => :photo_taken,  # shutter-press
    "SCM" => :second_crew_member,  # now use CM2
)

const OBSOLETE_ENGINE_CODE_TLC = Dict(
    "EDN" => :engine_down,
    "EOF" => :engine_off,
    "EON" => :engine_on,
    "EUP" => :engine_up,
)


text:
# Run tests on block Newton

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@bench @testset "blocknewton" begin
    nep=nep_gallery("dep0",4);
    V = Matrix(1.0I, size(nep,1), 3)
    S0=zeros(3,3);
    S,V=blocknewton(nep,S=S0,X=V,logger=displaylevel,
                    armijo_factor=0.5,
                    maxit=20)

    λv=eigvals(S);
    for λ=λv
        @test minimum(svdvals(compute_Mder(nep,λ)))<sqrt(eps())
    end

end


text:
###############
# Status type #
###############

type Status <: GitLabType
    id::Nullable{Int}
    total_count::Nullable{Int}
    state::Nullable{GitLabString}
    description::Nullable{GitLabString}
    context::Nullable{GitLabString}
    sha::Nullable{GitLabString}
    url::Nullable{HttpCommon.URI}
    target_url::Nullable{HttpCommon.URI}
    created_at::Nullable{Dates.DateTime}
    updated_at::Nullable{Dates.DateTime}
    creator::Nullable{Owner}
    repository::Nullable{Repo}
    statuses::Nullable{Vector{Status}}

    ## For commit status
    status::Nullable{GitLabString}
    name::Nullable{GitLabString}
    author::Nullable{Owner}
    ref::Nullable{GitLabString}
    started_at::Nullable{Dates.DateTime}
    finished_at::Nullable{Dates.DateTime}
    allow_failure::Nullable{Bool}
end

Status(data::Dict) = json2gitlab(Status, data)
Status(id::Real) = Status(Dict("id" => id))

namefield(status::Status) = status.id

###############
# API Methods #
###############

function create_status(repo, sha; options...)
    result = gh_post_json("/api/v3/projects/$(get(repo.id))/statuses/$(name(sha))"; options...)
    return Status(result)
end

function statuses(repo, ref; options...)
    results, page_data = gh_get_paged_json("/api/v3/projects/$(get(repo.id))/repository/commits/$(name(ref))/statuses"; options...)
    return map(Status, results), page_data
end

#= TODO: no equivalent API
function status(repo, ref; options...)
    result = gh_get_json("/api/v3/projects/$(get(repo.id))/commits/$(name(ref))/status"; options...)
    return Status(result)
end
=#


text:
using Printf

function calculate_area(a, b)
    π * a * b
end

function read_input(input_dir, v)
    open(@sprintf "%s/%s" input_dir v) do file
        parse.(Float64, read(file, String))
    end
end

function write_output(output_dir, output_file, v)
    output_path = @sprintf "%s/%s" output_dir output_file
    open(output_path, "w") do file
        write(file, string(v))
    end
end

function main(input_dir, output_dir)
    a = read_input(input_dir, 'a')
    b = read_input(input_dir, 'b')

    area = calculate_area(a, b)

    write_output(output_dir, "area", area)
    write_output(output_dir, "metadata", "[from julia rawcontainer]")
end

# the keyword ARGS is a special value that contains the command-line arguments
# julia arrays are 1-indexed
input_dir = ARGS[1]
output_dir = ARGS[2]

main(input_dir, output_dir)

max_stars_repo_name: kennyworkman/flytesnacks

max_issues_repo_name: kennyworkman/flytesnacks

max_forks_repo_name: kennyworkman/flytesnacks


text:
"""
    test_attribute_NumberThreads(model::MOI.AbstractOptimizer, config::Config)

Test that the [`MOI.NumberOfThreads`](@ref) attribute is implemented for
`model`.
"""
function test_attribute_NumberThreads(model::MOI.AbstractOptimizer, ::Config)
    @requires MOI.supports(model, MOI.NumberOfThreads())
    # Get the current value to restore it at the end of the test
    value = MOI.get(model, MOI.NumberOfThreads())
    MOI.set(model, MOI.NumberOfThreads(), 1)
    @test MOI.get(model, MOI.NumberOfThreads()) == 1
    MOI.set(model, MOI.NumberOfThreads(), 3)
    @test MOI.get(model, MOI.NumberOfThreads()) == 3
    MOI.set(model, MOI.NumberOfThreads(), value)
    @test value == MOI.get(model, MOI.NumberOfThreads())
    _test_attribute_value_type(model, MOI.NumberOfThreads())
    return
end
test_attribute_NumberThreads(::MOI.ModelLike, ::Config) = nothing

function setup_test(
    ::typeof(test_attribute_NumberThreads),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOI.set(model, MOI.NumberOfThreads(), nothing)
    return
end

"""
    test_attribute_RawStatusString(model::MOI.AbstractOptimizer, config::Config)

Test that the [`MOI.RawStatusString`](@ref) attribute is implemented for
`model`.
"""
function test_attribute_RawStatusString(
    model::MOI.AbstractOptimizer,
    config::Config,
)
    @requires _supports(config, MOI.optimize!)
    @requires _supports(config, MOI.RawStatusString)
    MOI.add_variable(model)
    MOI.optimize!(model)
    _test_attribute_value_type(model, MOI.RawStatusString())
    return
end
test_attribute_RawStatusString(::MOI.ModelLike, ::Config) = nothing

function setup_test(
    ::typeof(test_attribute_RawStatusString),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(
                mock,
                MOI.RawStatusString(),
                "Mock solution set by `mock_optimize!`.",
            )
        end,
    )
    return
end

"""
    test_attribute_Silent(model::MOI.AbstractOptimizer, config::Config)

Test that the [`MOI.Silent`](@ref) attribute is implemented for `model`.
"""
function test_attribute_Silent(model::MOI.AbstractOptimizer, ::Config)
    @requires MOI.supports(model, MOI.Silent())
    # Get the current value to restore it at the end of the test
    value = MOI.get(model, MOI.Silent())
    MOI.set(model, MOI.Silent(), !value)
    @test !value == MOI.get(model, MOI.Silent())
    # Check that `set` does not just take `!` of the current value
    MOI.set(model, MOI.Silent(), !value)
    @test !value == MOI.get(model, MOI.Silent())
    MOI.set(model, MOI.Silent(), value)
    @test value == MOI.get(model, MOI.Silent())
    _test_attribute_value_type(model, MOI.Silent())
    return
end
test_attribute_Silent(::MOI.ModelLike, ::Config) = nothing

function setup_test(
    ::typeof(test_attribute_Silent),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOI.set(model, MOI.Silent(), true)
    return
end

"""
    test_attribute_SolverName(model::MOI.AbstractOptimizer, config::Config)

Test that the [`MOI.SolverName`](@ref) attribute is implemented for `model`.
"""
function test_attribute_SolverName(model::MOI.AbstractOptimizer, config::Config)
    if _supports(config, MOI.SolverName)
        _test_attribute_value_type(model, MOI.SolverName())
    end
    return
end
test_attribute_SolverName(::MOI.ModelLike, ::Config) = nothing

"""
    test_attribute_SolveTimeSec(model::MOI.AbstractOptimizer, config::Config)

Test that the [`MOI.SolveTimeSec`](@ref) attribute is implemented for `model`.
"""
function test_attribute_SolveTimeSec(
    model::MOI.AbstractOptimizer,
    config::Config,
)
    @requires _supports(config, MOI.optimize!)
    @requires _supports(config, MOI.SolveTimeSec)
    MOI.add_variable(model)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.SolveTimeSec()) >= 0.0
    _test_attribute_value_type(model, MOI.SolveTimeSec())
    return
end
test_attribute_SolveTimeSec(::MOI.ModelLike, ::Config) = nothing

function setup_test(
    ::typeof(test_attribute_SolveTimeSec),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOI.set(mock, MOI.SolveTimeSec(), 0.0),
    )
    return
end

"""
    test_attribute_TimeLimitSec(model::MOI.AbstractOptimizer, config::Config)

Test that the [`MOI.TimeLimitSec`](@ref) attribute is implemented for `model`.
"""
function test_attribute_TimeLimitSec(model::MOI.AbstractOptimizer, ::Config)
    @requires MOI.supports(model, MOI.TimeLimitSec())
    # Get the current value to restore it at the end of the test
    value = MOI.get(model, MOI.TimeLimitSec())
    MOI.set(model, MOI.TimeLimitSec(), 0.0)
    @test MOI.get(model, MOI.TimeLimitSec()) == 0.0
    MOI.set(model, MOI.TimeLimitSec(), 1.0)
    @test MOI.get(model, MOI.TimeLimitSec()) == 1.0
    MOI.set(model, MOI.TimeLimitSec(), value)
    @test value == MOI.get(model, MOI.TimeLimitSec()) # Equality should hold
    _test_attribute_value_type(model, MOI.TimeLimitSec())
    return
end
test_attribute_TimeLimitSec(::MOI.ModelLike, ::Config) = nothing

function setup_test(
    ::typeof(test_attribute_TimeLimitSec),
    model::MOIU.MockOptimizer,
    ::Config,
)
    MOI.set(model, MOI.TimeLimitSec(), nothing)
    return
end

"""
    test_attribute_after_empty(model::MOI.AbstractOptimizer, config::Config)

Test that optimizer attributes such as `Silent` are not cleared by `MOI.empty!`.
"""
function test_attribute_after_empty(model::MOI.AbstractOptimizer, ::Config)
    @requires MOI.supports(model, MOI.Silent())
    current = MOI.get(model, MOI.Silent())
    for value in (true, false)
        MOI.set(model, MOI.Silent(), value)
        @test MOI.get(model, MOI.Silent()) == value
        MOI.empty!(model)
        @test MOI.get(model, MOI.Silent()) == value
    end
    # Make sure to reset the value before leaving this function!
    MOI.set(model, MOI.Silent(), current)
    return
end

test_attribute_after_empty(::MOI.ModelLike, ::Config) = nothing