text:
"""
    bea_table(dataset::String, TableName::String, frequency::String,
        startyear::Int, endyear::Int; user_id = USER_ID) -> BeaTable

Return a [`BeaTable`](@ref) with data and metadata for `TableName`.  Pass
integer value `0` for both `startyear` and `endyear` to retrieve all available years for
the table. The data values are returned in a `DataFrame` accessed through the
`data_values` field of the `BeaTable` struct.

The `TableName` argument refers to the API `TableName` parameter value for the
requested table.  These can be found using the [`bea_parametervalues`](@ref) method.

Example: to get Table 1.1.6 in the NIPA dataset, quarterly, from 2015 to 2018.
```julia
# User ID stored in ~/.beadatarc
julia> nipa116 = bea_table("NIPA", "T10106", "Q", 2015, 2018)
BEA Table
Dataset:   NIPA
Table No.: 1.1.6
Title:     Real Gross Domestic Product, Chained Dollars
Metric:    Chained Dollars
Units:     Billions of chained (2012) dollars
Frequency: Q
Dates:     2015 - 2018
Revised:   August 29, 2019
```
"""
function bea_table(dataset::String, TableName::String, frequency::String,
    startyear::Int, endyear::Int; user_id = USER_ID)

    supported = ["NIPA", "NIUnderlyingDetail", "FixedAssets"]
    if dataset ∉ supported
        error("Dataset $dataset not supported.")
    end

    api_method = "GetData"
    years = (startyear == 0 && endyear == 0) ? "X" : join(collect(startyear:1:endyear), ",")

    querydict = Dict("UserID" => user_id,
                     "Method" => api_method,
                     "DatasetName" => dataset,
                     "TableName" => TableName,
                     "Frequency" => frequency,
                     "Year" => years,
                     "ResultFormat" => "JSON")

     r_json = bea_query(API_URL, querydict)
     data_dicts = r_json["Data"]
     table_notes_dict = r_json["Notes"]

    # Extract metadata for the table
    tablenum, tabledesc, title_units, notes, api_tablename, last_revised =
        parse_table_metadata(table_notes_dict)
    api_units  = data_dicts[1]["CL_UNIT"]
    metric = data_dicts[1]["METRIC_NAME"]
    units = api_units == "Level" ? title_units : api_units
    unit_mult = data_dicts[1]["UNIT_MULT"]

    # Extract the data
    all_data = [parse_data_dict(d) for d in data_dicts]
    linenums = [tup[1] for tup in all_data]
    linenums = [length(num) == 1 ? "0"*num : num for num in linenums]
    varnames = [tup[2] for tup in all_data]
    sercodes = [tup[3] for tup in all_data]
    dates    = [tup[4] for tup in all_data]

    # Resolve difference in reported units and actual units:
    # billions vs. millions
    datavals = [tup[5] for tup in all_data]
    if unit_mult == "6" && occursin("Billions", units)
        datavals ./= 1000
    end

    dflong = DataFrame(TimePeriod = dates, LineNumber = linenums, value = datavals)
    # Convert to "wide" format: unstack(dfname, :id, :variable, :value)
    # - :id identifies the rows (here it's the date)
    # - :variable identifies what will become columns (here, the line number)
    # - :value identifies the data values
    dfwide = unstack(dflong, :TimePeriod, :LineNumber, :value)
    newnames = [Symbol(string("Line", name)) for name in names(dfwide)[2:end]]
    rename!(dfwide, [:TimePeriod; newnames])
    actual_start = Dates.year(dfwide[!, :TimePeriod][1])
    actual_end   = Dates.year(dfwide[!, :TimePeriod][end])

    linekeys = unique(DataFrame(LineNumber = linenums, SeriesCode = sercodes,
        LineDescription = varnames))

    fullfreq = frequency == "A" ? "Annual" :
        (frequency == "Q" ? "Quarterly" : "Monthly")

    return BeaTable(dataset, tablenum, tabledesc, metric, units, linekeys,
        notes, fullfreq, actual_start, actual_end, api_tablename, last_revised, dfwide)
end

function parse_table_metadata(d)
    tableinfo = d[1]["NoteText"]
    m1 = match(r"\s(?<tablenum>.*)\.\s", tableinfo)
    tablenum = m1[:tablenum] != nothing ? String(m1[:tablenum]) : ""
    m2 = match(r"\.\s(?<tabledesc>.*)\s\[", tableinfo)
    m2 = m2 == nothing ? match(r"\.\s(?<tabledesc>.*)\s\-", tableinfo) : m2
    tabledesc = m2 != nothing ? String(m2[:tabledesc]) : ""
    m3 = match(r"\[(?<units>.*)\]", tableinfo)
    title_units = m3 != nothing ? String(m3[:units]) : ""
    m4 = match(r"(?=LastRevised: (?<revised>.*))", tableinfo)
    last_revised = m4 != nothing ? m4[:revised] : ""

    api_tablename = d[1]["NoteRef"]

    noteref = [note["NoteRef"] for note in d]
    notetext = [note["NoteText"] for note in d]
    notes = DataFrame(noteref = noteref, notetext = notetext)

    return tablenum, tabledesc, title_units, notes, api_tablename, last_revised
end

function parse_data_dict(d::Dict)
    #-Line number and variable name
    linenum = d["LineNumber"]
    varname = d["LineDescription"]
    seriescode = d["SeriesCode"]
    # Add footnote indicator to varname if present
    if haskey(d, "NoteRef") && d["NoteRef"] != d["TableName"]
        fn = []
        for i in split(d["NoteRef"], ",")
             for j in split(i, ".")
                 if j != d["TableName"]
                     push!(fn, j)
                 end
             end
        end
        varname = string(varname, " (", join(fn, ","), ")")
    end

    #--Get the date of the observation--#
    # Convert from String to Date() type
    timeperiod = d["TimePeriod"]
    obs_year = parse(Int, timeperiod[1:4])
    #- Monthly data
    if occursin("M", timeperiod)
        obs_month = parse(Int, timeperiod[6:end])
        date  = Date(obs_year, obs_month, 1)
    #- Quarterly data
    elseif occursin("Q", timeperiod)
        obs_quarter = parse(Int, timeperiod[end])
        date = Date(obs_year, (obs_quarter*3), 1)
    #- Annual data
    else
        date = Date(obs_year, 12, 1)
    end

    #-Get data value-#
    # Remove commas and convert to float
    dataval = parse(Float64, replace(d["DataValue"], "," => ""))

    return (linenum, varname, seriescode, date, dataval)
end

"""
    struct BeaTable

A BEA table with data and metadata returned from a [`bea_table`](@ref) call.

Fields
------
* `dataset`: dataset the table was retrieved from
* `table_number`: the non-API table number
* `table_description`: description of the data contained in the table
* `metric`: measurement metric for the data, e.g. index, dollars, etc.
* `units`: billions, millions, etc.
* `line_descriptions`: `DataFrame` containing line number descriptions
* `table_notes`: `DataFrame` containing any notes for the table
* `frequency`: (M)onthly, (Q)uarterly, or (A)nnual
* `data_startyear`: first year of data returned (may differ from what was requested)
* `data_endyear`: last year of data returned (may differ from what was requested)
* `api_tablename`: `TableName` parameter value for the table
* `last_revised`: date the table was last revised
* `data_values`: `DataFrame` containing the table data values
"""
struct BeaTable
    dataset::String
    table_number::String
    table_description::String
    metric::String
    units::String
    line_descriptions::DataFrame
    table_notes::DataFrame
    frequency::String
    data_startyear::Int
    data_endyear::Int
    api_tablename::String
    last_revised::String
    data_values::DataFrame
end

function Base.show(io::IO, b::BeaTable)
    println(io, "BEA Table")
    println(io, "Dataset:     $(b.dataset)")
    println(io, "Table No.:   $(b.table_number)")
    println(io, "Description: $(b.table_description)")
    println(io, "Metric:      $(b.metric)")
    println(io, "Units:       $(b.units)")
    println(io, "Frequency:   $(b.frequency)")
    println(io, "Dates:       $(b.data_startyear) - $(b.data_endyear)")
    println(io, "Revised:     $(b.last_revised)")
end

struct BeaDataError
    dataset::String
    api_tablename::String
    frequency::String
    startyear::Int
    endyear::Int
    emsg::String
end

function Base.show(io::IO, b::BeaDataError)
    println(io, "BEA Data Error")
    println(io, "$(b.emsg)")
    println(io, "Dataset:   $(b.dataset)")
    println(io, "TableName: $(b.api_tablename)")
    println(io, "Frequency: $(b.frequency)")
    println(io, "Dates:     $(b.startyear) - $(b.endyear)")
end


text:
"""
  read_priors(pathToCsv::String)

Function to load priors stored in a csv file. Returns a dictionary. The csv
file should have the following columns "name", "value", "upper_bound", "lower_bound"
"""
function read_priors(pathToCsv::String)

    # initialize and empty Dictionary:
    #---------------------------------
    dictionary = OrderedDict{String,Array{Float64,1}}()

    # Check the existence of the file
    #----------------------------------
    if isfile(pathToCsv) == false
        error(string(pathToCsv," not found."))
    end

    # If file exists
    #----------------------------------
    #dataFrame = CSV.read(pathToCsv)
    dataFrame = CSV.File(open(read, pathToCsv)) |> DataFrame

    # Check whether the csv file has the appropriate columns
    #-------------------------------------------------------
    list_to_check = ["name", "value", "upper_bound", "lower_bound"]
    col_names = names(dataFrame)

    for col_name in list_to_check
        if in(col_name, col_names) == false
            error("column $(col_name) does not exist.")
        end
    end

    # store the number of rows in the dataframe:
    #------------------------------------------
    number_rows = size(dataFrame, 1)

    # Append the dictionary
    #----------------------
    for i=1:number_rows
        dictionary[dataFrame[i,:name]] = [dataFrame[i,:value],  dataFrame[i,:lower_bound], dataFrame[i,:upper_bound]]
    end

    return dictionary

end

"""
  read_empirical_moments(pathToCsv::String)

Function to load empirical moments stored in a csv file. Returns a dictionary.
The csv file should have the following columns "name", "value", "weight"
"""
function read_empirical_moments(pathToCsv::String)

    # initialize and empty Dictionary:
    #---------------------------------
    dictionary = OrderedDict{String,Array{Float64,1}}()

    # Check the existence of the file
    #----------------------------------
    if isfile(pathToCsv) == false
        error(string(pathToCsv," not found"))
    end

    # If file exists:
    #----------------
    #dataFrame =  CSV.read(pathToCsv)
    dataFrame = CSV.File(open(read, pathToCsv)) |> DataFrame

    #check whether the csv file has the columns "name" and "value"
    #-------------------------------------------------------------
    list_to_check = ["name", "value", "weight"]
    col_names = names(dataFrame)

    for col_name in list_to_check
        if in(col_name, col_names) == false
            error("column $(col_name) does not exist.")
        end
    end

    # store the number of rows in the dataframe:
    #------------------------------------------
    number_rows = size(dataFrame,1)

    # Append the dictionary
    #----------------------
    for i=1:number_rows
        dictionary[dataFrame[i,:name]] = [dataFrame[i,:value],  dataFrame[i,:weight]]
    end

    return dictionary

end

ext: jl


text:
struct ColouredScheme{T}
    opacity::T
end
function (scheme::ColouredScheme)(v)
    0.0 <= v <= 1.0 || println("Warning: color value must be scaled between 0 and 1.")
    if 0.0 <= v < 0.25
        b = 1.0
        g = v/0.25
        r = 0.0
    elseif 0.25 <= v < 0.5
        b = 1.0 - (v - 0.25)/0.25
        g = 1.0
        r = 0.0
    elseif 0.5 <= v < 0.75
        b = 0.0
        g = 1.0
        r = (v - 0.5)/0.25
    else
        b = 0.0
        g = 1.0 - (v - 0.75)/0.25
        r = 1.0
    end
    return r, g, b, scheme.opacity
end

struct BlackScheme end
function (::BlackScheme)(v)
    0.0 <= v <= 1.0 || println("Warning: color value must be scaled between 0 and 1.")
    return 0.7, 0.7, 0.7, v
end

function get_jl_mapped_colors(dataset::T, color_variable_name, component, opacity, color_scheme = BlackScheme(), scale = true) where {T<:AbstractVTKSimpleData}
    if color_variable_name in keys(dataset.point_data)
        point_color = true
        _size = num_of_points(dataset)
        _color_variable = dataset.point_data[color_variable_name]
        _var_dim = var_dim(dataset, color_variable_name, "Point")
    elseif color_variable_name in keys(dataset.cell_data)
        point_color = false
        _size = num_of_cells(dataset)
        _color_variable = dataset.cell_data[color_variable_name]
        _var_dim = var_dim(dataset, color_variable_name, "Cell")
    else
        _size = num_of_points(dataset)
        _color_variable = nothing
    end

    if T <: AbstractVTKStructuredData
        if _var_dim == 1
            color_variable = reshape(_color_variable, (_size,))
        else
            color_variable = reshape(_color_variable, (_var_dim, _size))
        end
    else
        color_variable = _color_variable
    end

    if color_variable !== nothing
        if component == -1 || component == 1 && _var_dim == 1
            if _var_dim == 1
                _values = color_variable
                cmin = minimum(_values)
                cmax = maximum(_values)
            else
                @views _values = [norm(color_variable[:,i]) for i in 1:_size]
                cmin = minimum(_values)
                cmax = maximum(_values)
            end
        else
            if 1 <= component <= _var_dim
                @views _values = color_variable[component,:]
                cmin = minimum(_values)
                cmax = maximum(_values)
            else
                throw("Cannot use component $component of a $color. $color only has $_var_dim components.")
            end
        end
    else
        cmin = cmax = 1.0
        _values = ones(_size)
    end
    jl_mapped_colors = Vector{RGBA{Float32}}(undef, _size)
    if !(isapprox(cmax, cmin, atol = eps(Float32)))
        scaled_values = (_values .- cmin) ./ (cmax - cmin)
    else
        scaled_values = (similar(_values) .= clamp(cmin, 0, 1))
    end
    function fill_jl_mapped_colors(i)
        jl_mapped_colors[i] = RGBA{Float32}((color_scheme(scaled_values[i]))...)
    end
    map(fill_jl_mapped_colors, 1:_size)
    return jl_mapped_colors, cmin, cmax
end

function GLMesh(rectilinear::AbstractVTKStructuredData; kwargs...)
    GLMesh(VTKUnstructuredData(rectilinear); kwargs...)
end
function GLMesh(dataset::AbstractVTKUnstructuredData; color::String="", component::Int=-1, opacity::Float64=1.0, color_scheme = BlackScheme())
    filter_cells!(dataset, [POINT_CELLS; LINE_CELLS])
    tri_cell_data = color in keys(dataset.cell_data)
    tri_dataset = triangulate(dataset, tri_cell_data, GLTriangle)
    tri_dataset = tri_dataset |> remove_unused_vertices |> duplicate_vertices
    if haskey(tri_dataset.cell_data, color)
        celldata_to_pointdata!(tri_dataset)
    end
    color_vec, cmin, cmax = get_jl_mapped_colors(   tri_dataset,
                                                    color,
                                                    component,
                                                    opacity,
                                                    color_scheme
                                                )
    T = eltype(dataset.point_coords)
    vertices = map(1:num_of_points(tri_dataset)) do i
        @views coord = tri_dataset.point_coords[:,i]
        if length(coord) == 2
            Point{3, T}(coord..., zero(T))
        else
            Point{3, T}(coord...)
        end
    end
    faces = tri_dataset.cell_connectivity
    return GLNormalVertexcolorMesh(vertices=vertices, faces=faces, color=color_vec)
end


text:
module FastRationals

export FastRational

using Base: BitInteger, BitSigned, BitUnsigned

using Base.Checked: add_with_overflow, sub_with_overflow, mul_with_overflow,
    checked_neg, checked_abs, checked_add, checked_sub, checked_mul,
    checked_div, checked_rem, checked_fld, checked_mod, checked_cld

import Base: hash, show, repr, string, tryparse,
    zero, one, iszero, isone,
    numerator, denominator, eltype, convert, promote_rule, decompose,
    isinteger, typemax, typemin, sign, signbit, copysign, flipsign, abs, float,
    ==, !=, <, <=, >=, >,
    +, -, *, /, ^, //,
    inv, div, fld, cld, rem, mod, trunc, floor, ceil, round

const FastSigned   = Union{Int8, Int16, Int32, Int64}
const FastUnsigned = Union{UInt8, UInt16, UInt32, UInt64}
const FastInteger  = Union{FastSigned, FastUnsigned}

"""
    RationalState

a state applicable to Rational values
"""
abstract type RationalState end

"""
    IsReduced <: RationalState

This state holds for rational values that are known to have been reduced to lowest terms.
"""
struct IsReduced  <: RationalState end

"""
    MayReduce <: RationalState

This state holds for rational values that may or may not be expressed in lowest terms.
"""
struct MayReduce  <: RationalState end

"""
    FastRational <: Real

### fields

`num`: the rational numerator      (zero and signed integer values)
`den`: the rational denominator    (strictly positive integer values)

"""
struct FastRational{T<:BitInteger, H<:RationalState} <: Real
    num::T
    den::T
end

#=
    FastRational constructors *never*  process 2Tupled pairs
        `num`, `den` already must be reduced to lowest terms
         ( (num::T, den::T) ) ↦ (num, den)

    FastRational constructors *always* process untupled pairs
        `num`, `den` may not yet be reduced to lowest terms
         (num, den) ↦ canonical(num, den)  ↦ (num, den)
=#

FastRational(x::NTuple{2,T}) where {T<:BitInteger} = FastRational{T,IsReduced}(x[1], x[2])
FastRational(num::T, den::T) where {T<:BitInteger} = FastRational(canonical(num, den))

numerator(x::FastRational{T,H})   where {T,H} = x.num
denominator(x::FastRational{T,H}) where {T,H} = x.den
numerdenom(x::FastRational{T,H})  where {T,H} = x.num, x.den

basistype(::Type{FastRational{T,IsReduced}}) where T = T
basistype(::Type{FastRational{T,MayReduce}}) where T = T
basistype(::Type{Rational{T}}) where T = T

FastRational(x::FastRational{T,IsReduced}) where {T} = x
FastRational(x::FastRational{T,MayReduce}) where {T} = FastRational(canonical(x.num, x.den))

FastRational(x::Rational{T}) where {T} = FastRational{T,IsReduced}(x.num, x.den)
FastRational(x::T) where {T<:Integer} = FastRational(x, one(T))
FastRational(x::T) where {T<:Rational} = FastRational(x)
FastRational(x::T) where {T<:Real} = FastRational(rationalize(x))
FastRational{T}(x::T) where {T} = FastRational{T,IsReduced}(x)

FastRational{T,IsReduced}(x::Rational{T}) where {T} = FastRational{T,IsReduced}(x.num, x.den)
FastRational{T,MayReduce}(x::Rational{T}) where {T} = FastRational{T,IsReduced}(x.num, x.den)
FastRational{T1,IsReduced}(x::Rational{T2}) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(T1(x.num), T1(x.den))
FastRational{T1,MayReduce}(x::Rational{T2}) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(T1(x.num), T1(x.den))

FastRational{T,IsReduced}(x::T) where {T<:BitInteger} = FastRational{T,IsReduced}(x, one(T))
FastRational{T,MayReduce}(x::T) where {T<:BitInteger} = FastRational{T,IsReduced}(x, one(T))
FastRational{T1,IsReduced}(x::T2) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(T1(x), one(T1))
FastRational{T1,MayReduce}(x::T2) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(T1(x), one(T1))

FastRational{T1,IsReduced}(x::FastRational{T2,IsReduced},y::FastRational{T2,IsReduced}) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(T1(x.num), T1(x.den))
FastRational{T1,MayReduce}(x::FastRational{T2,IsReduced},y::FastRational{T2,IsReduced}) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(T1(x.num), T1(x.den))
FastRational{T1,IsReduced}(x::FastRational{T2,MayReduce},y::FastRational{T2,MayReduce}) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(canonical(T1(x.num), T1(x.den)))
FastRational{T1,MayReduce}(x::FastRational{T2,MayReduce},y::FastRational{T2,MayReduce}) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(canonical(T1(x.num), T1(x.den)))
FastRational{T1,IsReduced}(x::FastRational{T2,IsReduced},y::FastRational{T2,MayReduce}) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(T1(x.num), T1(x.den))
FastRational{T1,MayReduce}(x::FastRational{T2,MayReduce},y::FastRational{T2,IsReduced}) where {T1<:BitInteger,T2} = FastRational{T1,IsReduced}(T1(x.num), T1(x.den))

FastRational{T,H}(x::FastRational{T,IsReduced},y::FastRational{T,MayReduce}) where {T<:BitInteger,H} = FastRational{T,IsReduced}(canonical(T(x.num), T(x.den)))
FastRational{T,H}(x::FastRational{T,MayReduce},y::FastRational{T,IsReduced}) where {T<:BitInteger,H} = FastRational{T,IsReduced}(canonical(T(x.num), T(x.den)))

#=
FastRational{Int64,H}(x::FastRational{Int64,IsReduced}, y::FastRational{Int64,MayReduce}) where {H<:Main.FastRationals.RationalState} =
    x * inv(FastRational(canonical(y)))
FastRational{Int64,H}(x::FastRational{Int64,MayReduce}, y::FastRational{Int64,IsReduced}) where {H<:Main.FastRationals.RationalState} =
    FastRational(canonical(x) * inv(y)
=#
FastRational{T,H}(x::FastRational{T,IsReduced}, y::FastRational{T,MayReduce}) where {T,H<:RationalState} =
    x * inv(FastRational(canonical(y)))
FastRational{T,H}(x::FastRational{T,MayReduce}, y::FastRational{T,IsReduced}) where {T,H<:RationalState} =
    FastRational(canonical(x) * inv(y))

Rational(x::FastRational{T,IsReduced}) where {T} = Rational{T}(x.num, x.den)
Rational(x::FastRational{T,MayReduce}) where {T} = Rational{T}(canonical(x.num, x.den))
Rational{T}(x::FastRational{T,IsReduced}) where {T} = Rational{T}(x.num, x.den)
Rational{T}(x::FastRational{T,MayReduce}) where {T} = Rational{T}(canonical(x.num, x.den))

//(x::FastRational{T,H}, y::Integer) where {T,H} = x / FastRational(y)
//(x::Integer, y::FastRational{T,H}) where {T,H} = FastRational(x) / y
//(x::FastRational, y::FastRational) = x / y
//(x::FastRational{T,H}, y::Rational) where {T,H} = x / FastRational(y)
//(x::Rational, y::FastRational{T,H}) where {T,H} = FastRational(x) / y

float(x::FastRational{T,H}) where {T,H} = float(Rational(x))
Base.Float64(x::FastRational{T,H}) where {T,H} = Float64(Rational(x))
Base.Float32(x::FastRational{T,H}) where {T,H} = Float32(Rational(x))
Base.BigFloat(x::FastRational{T,H}) where {T,H} = BigFloat(Rational(x))
Base.BigInt(x::FastRational{T,H}) where {T,H} = BigInt(Rational(x))

include("canonical.jl")

include("promote_convert.jl")
# ------------- promotion rules and conversion logic above has been reviewed

# convert(::Type{Rational{T}}, x::FastRational{T,H}) where {T,H<:RationalState} = Rational(canonical(x.num, x.den))

convert(::Type{FastRational{T,H}}, x::T) where {T,H} = FastRational{T,IsReduced}(x, one(T))
convert(::Type{FastRational{T1,H}}, x::T2) where {T1,T2<:Integer,H} = FastRational{T1,IsReduced}(T1(x), one(T1))
convert(::Type{FastRational{T1,H}}, x::Rational{T2}) where {T1,T2,H} =
    FastRational{T1,IsReduced}(T1(x.num), T1(x.den))
# disambiguate
convert(::Type{FastRational{T,H}}, x::FastRational{T,H}) where {T,H} = x
convert(::Type{FastRational{T,IsReduced}}, x::FastRational{T,MayReduce}) where {T} = FastRational(x)
convert(::Type{FastRational{T,MayReduce}}, x::FastRational{T,IsReduced}) where {T} = x

convert(::Type{FastRational{T,H}}, x::AbstractFloat) where {T,H} = FastRational(convert(Rational{T}, x))
convert(::Type{F}, x::FastRational{T,H}) where {T,H,F<:AbstractFloat} = F(convert(Rational{T}, x))

zero(::Type{FastRational{T,H}}) where {T,H} = FastRational{T,IsReduced}(zero(T), one(T))
one(::Type{FastRational{T,H}}) where {T,H} = FastRational{T,IsReduced}(one(T), one(T))

signbit(x::FastRational{T,H}) where {T<:Signed, H} = signbit(x.num) !== signbit(x.den)
sign(x::FastRational{T,H}) where {T<:Signed, H} = FastRational{T,IsReduced}(signbit(x) ? -one(T) : one(T))
abs(x::FastRational{T,IsReduced}) where {T<:Signed} = FastRational{T,IsReduced}(abs(x.num), x.den)
abs(x::FastRational{T,MayReduce}) where {T<:Signed} = FastRational{T,IsReduced}(abs(x.num), abs(x.den))

signbit(x::FastRational{T,H}) where {T<:Unsigned, H} = false
sign(x::FastRational{T,H}) where {T<:Unsigned, H} = FastRational{T,IsReduced}(one(T), one(T))
abs(x::FastRational{T,H}) where {T<:Unsigned, H} = x

# optimize for FastRational{Int32} using Int64 cross-multiplication for comparisons
@inline mulwider(x::T, y::T) where {T<:Integer} = widemul(x,y)

(==)(x::FastRational{T,IsReduced}, y::FastRational{T,IsReduced}) where {T} =
    x.num === y.num && x.den === y.den
(!=)(x::FastRational{T,IsReduced}, y::FastRational{T,IsReduced}) where {T} =
    x.num !== y.num || x.den !== y.den

for F in (:(==), :(!=), :(<), :(<=), :(>=), :(>))
  @eval $F(x::FastRational{T,H1}, y::FastRational{T,H2}) where {T,H1, H2} =
    $F(mulwider(x.num, y.den), mulwider(x.den, y.num))
end

copysign(x::FastRational, y::Real) = y >= 0 ? abs(x) : negate(abs(x))
copysign(x::FastRational{T,H}, y::FastRational) where {T,H} = FastRational{T,H}(copysign(x.num, y.num), x.den)
copysign(x::FastRational{T,H}, y::Rational) where {T,H} = FastRational{T,H}(copysign(x.num,y.num), x.den)

flipsign(x::FastRational{T,H}, y::Real) where {T,H} = FastRational{T,H}(flipsign(x.num,y), x.den)
flipsign(x::FastRational{T,H}, y::FastRational) where {T,H} = FastRational{T,H}(flipsign(x.num,y.num), x.den)
flipsign(x::FastRational{T,H}, y::Rational) where {T,H} = FastRational{T,H}(flipsign(x.num,y.num), x.den)

abs(x::FastRational{T,H}) where {T,H} = FastRational{T,H}(abs(x.num), x.den)

typemin(::Type{FastRational{T,H}}) where {T<:Signed,H} = FastRational{T,IsReduced}(typemin(T),one(T))
typemin(::Type{FastRational{T,H}}) where {T<:Unsigned,H} = FastRational{T,IsReduced}(zero(T),one(T))
typemax(::Type{FastRational{T,H}}) where {T<:Integer,H} = FastRational{T,IsReduced}(typemax(T),one(T))

isinteger(x::FastRational{T,IsReduced}) where {T} = x.den == 1
isinteger(x::FastRational{T,MayReduce}) where {T} = canonical(x.num,x.den)[2] == 1

+(x::FastRational{T,H}) where {T,H} = FastRational{T,H}(+x.num, x.den)
-(x::FastRational{T,H}) where {T,H} = FastRational{T,H}(-x.num, x.den)

inv(x::FastRational{T,IsReduced}) where {T<:BitInteger} =
   FastRational{T,IsReduced}(x.den, x.num)
inv(x::FastRational{T,MayReduce}) where {T<:BitInteger} =
   FastRational{T,MayReduce}(x.den, x.num)

function -(x::FastRational{T,IsReduced}) where {T<:BitSigned}
    x.num == typemin(T) && throw(OverflowError("rational numerator is typemin(T)"))
    FastRational{T,IsReduced}(-x.num, x.den)
end
-(x::FastRational{T,MayReduce}) where {T<:BitSigned} =
    -FastRational{T,IsReduced}(canonical(x.num, x.den)...,)

function -(x::FastRational{T,H}) where {T<:Unsigned,H}
    x.num != zero(T) && throw(OverflowError("cannot negate unsigned number"))
    x
end

# core parts of add, sub

@inline function addovf(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(numerator(x), denominator(y)) # here, numer is a temp
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), numerator(y)) # here, denom is a temp
    ovf |= ovfl
    numer, ovfl = add_with_overflow(numer, denom) # numerator of sum
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), denominator(y)) # denominator of sum
    ovf |= ovfl
    return numer, denom, ovf
end

@inline function subovf(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(numerator(x), denominator(y)) # here, numer is a temp
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), numerator(y)) # here, denom is a temp
    ovf |= ovfl
    numer, ovfl = sub_with_overflow(numer, denom) # numerator of difference
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), denominator(y)) # denominator of difference
    ovf |= ovfl
    return numer, denom, ovf
end

@inline function mulovf(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(numerator(x), numerator(y))
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), denominator(y))
    ovf |= ovfl
    return numer, denom, ovf
end

@inline function divovf(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(numerator(x), denominator(y))
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), numerator(y))
    ovf |= ovfl
    return numer, denom, ovf
end

for (F,G) in ((:(+), :addovf), (:(-), :subovf), (:(*), :mulovf), (:(/), :divovf))
  @eval begin
    function $F(x::FastRational{T,IsReduced}, y::FastRational{T,IsReduced}) where {T}
      num, den, ovf = $G(x, y)
      ovf && throw(OverflowError("$x + $y overflowed"))
      return FastRational{T,MayReduce}(num, den)
    end
    function $F(x::FastRational{T,IsReduced}, y::FastRational{T,MayReduce}) where {T}
      num, den, ovf = $G(x, y)
      if ovf
         num, den = canonical(y.num, y.den)
         return $F(x, FastRational{T,IsReduced}(num, den))
      end
      return FastRational{T,MayReduce}(num, den)
    end
    function $F(x::FastRational{T,MayReduce}, y::FastRational{T,IsReduced}) where {T}
      num, den, ovf = $G(x, y)
      if ovf
         num, den = canonical(x.num, x.den)
         return $F(FastRational{T,IsReduced}(num, den), y)
      end
      return FastRational{T,MayReduce}(num, den)
    end
    function $F(x::FastRational{T,MayReduce}, y::FastRational{T,MayReduce}) where {T}
      num, den, ovf = $G(x, y)
      if ovf
         xnum, xden = canonical(x.num, x.den)
         ynum, yden = canonical(y.num, y.den)
         return $F(FastRational{T,IsReduced}(xnum, xden), FastRational{T,IsReduced}(ynum, yden))
      end
      return FastRational{T,MayReduce}(num, den)
    end
  end
end

const Solidus = "//"
const Solunos = "╱"

function string(x::FastRational{T,IsReduced}) where {T}
    return string(x.num,Solidus,x.den)
end

function string(x::FastRational{T,MayReduce}) where {T}
    num, den = canonical(x.num, x.den)
    return string(num,Solidus,den)
end

function show(io::IO, x::FastRational{T,H}) where {T,H}
    print(io, string(x))
end

function repr(x::FastRational{T,H}) where {T,H}
    return repr(Rational(x))
end

function tryparse(::Type{FastRational{T,H}}, s::String) where {T,H}
    if !occursin(Solunos, s)
        if !occursin(Solidus, s)
            num = tryparse(T, s)
            num = isnothing(num) : one(T) : num
            FastRational{T,IsReduced}(canonical(num, one(T)))
        else
            num, den = String.(split(s,Solidus))
            num = isempty(num) ? zero(T) : parse(T,num)
            den = isempty(den) ? one(T)  : parse(T,den)
            FastRational{T,IsReduced}(canonical(num, den))
        end
    else
        num, den = String.(split(s,Solunos))
        num = isempty(num) ? zero(T) : parse(T,num)
        den = isempty(den) ? one(T)  : parse(T,den)
        FastRational{T,IsReduced}(canonical(num, den))
    end
end

rem(x::FastRational{T,H1}, y::FastRational{T,H2}) where {T,H1,H2} = FastRational{T,IsReduced}(rem(Rational(x), Rational(y)))
mod(x::FastRational{T,H1}, y::FastRational{T,H2}) where {T,H1,H2} = FastRational{T,IsReduced}(mod(Rational(x), Rational(y)))

fld(x::FastRational{T,H1}, y::FastRational{T,H2}) where {T,H1,H2} = fld(Rational(x), Rational(y))
cld(x::FastRational{T,H1}, y::FastRational{T,H2}) where {T,H1,H2} = cld(Rational(x), Rational(y))
trunc(::Type{T}, x::FastRational) where {T} = trunc(Rational(x))
floor(::Type{T}, x::FastRational) where {T} = floor(Rational(x))
ceil(::Type{T}, x::FastRational) where {T} = ceil(Rational(x))
round(::Type{T}, x::FastRational, r::RoundingMode=RoundNearest) where {T} = round(T,Rational(x), r)
round(x::FastRational, r::RoundingMode) = round(Rational(x), r)

decompose(x::FastRational{T,IsReduced}) where {T} = x.num, zero(T), x.den
decompose(x::FastRational{T,MayReduce}) where {T} = decompose(FastRational(x))

end # FastRationals


text:
# * Variation of orbital overlaps

"""
    diff(ab::OrbitalOverlap, o::O)

Vary the orbital overlap ⟨a|b⟩ with respect to |o⟩.
"""
Base.diff(ab::OrbitalOverlap, o::O) where O =
    o == ab.b ? NBodyEquation(Conjugate(ab.a),
                              IdentityOperator{1}()) : 0

"""
    diff(ab::OrbitalOverlap, o::Conjugate{O})

Vary the orbital overlap ⟨a|b⟩ with respect to ⟨o|.
"""
Base.diff(ab::OrbitalOverlap, o::Conjugate{O}) where O =
    conj(o) == ab.a ? NBodyEquation(ab.b,
                                    IdentityOperator{1}()) : 0

# * Variation of orbital matrix elements

"""
    diff(ome::OrbitalMatrixElement, o::O)

Vary the orbital overlap ⟨abc...|Ω|xyz...⟩ with respect to |o⟩.
"""
function Base.diff(ome::OrbitalMatrixElement{N}, o::O) where {N,O}
    i = findfirst(isequal(o), ome.b)
    (isnothing(i) || iszero(ome)) && return 0
    NBodyEquation(Conjugate(ome.a[i]), contract(ome, complement(N, i)...))
end

"""
    diff(ome::OrbitalMatrixElement, o::Conjugate{O})

Vary the orbital overlap ⟨abc...|Ω|xyz...⟩ with respect to ⟨o|.
"""
function Base.diff(ome::OrbitalMatrixElement{N}, o::Conjugate{O}) where {N,O}
    i = findfirst(isequal(conj(o)), ome.a)
    (isnothing(i) || iszero(ome)) && return 0
    NBodyEquation(ome.b[i], contract(ome, complement(N, i)...))
end

# * Variation of NBodyMatrixElements

"""
    diff(me::NBodyMatrixElement, o::O)

Vary the [`NBodyMatrixElement`](@ref) `me` with respect to the orbital `o`.
"""
function Base.diff(me::NBodyMatrixElement, o::O) where O
    equations = Vector{NBodyEquation}()
    for t in me.terms
        for (i,f) in enumerate(t.factors)
            v = diff(f, o)
            iszero(v) && continue
            v = v * NBodyTerm(vcat(t.factors[1:(i-1)]...,
                                   t.factors[(i+1):end]...),
                                   t.coeff)
            push!(equations, v)
        end
    end
    LinearCombinationEquation(equations)
end

# * Variation of energy expression matrices

"""
    diff(E::Matrix{NBodyMatrixElement}, o::O)

Vary the matrix of [`NBodyMatrixElement`](@ref)s with respect to the
orbital `o`.

# Examples

```jldoctest
julia> E = Matrix(OneBodyHamiltonian()+CoulombInteraction(),
                  SlaterDeterminant.([[:a, :b], [:c, :d]]))
2×2 Array{EnergyExpressions.NBodyMatrixElement,2}:
 (a|a) + (b|b) - G(a,b) + F(a,b)  - [a b|d c] + [a b|c d]
 - [c d|b a] + [c d|a b]          (c|c) + (d|d) - G(c,d) + F(c,d)

julia> diff(E, :a)
2×2 SparseArrays.SparseMatrixCSC{LinearCombinationEquation,Int64} with 2 stored entries:
  [1, 1]  =  ⟨a|ĥ + -⟨b|[a|b] + ⟨a|[b|b]
  [2, 1]  =  -⟨d|[c|b] + ⟨c|[d|b]

julia> diff(E, Conjugate(:b))
2×2 SparseArrays.SparseMatrixCSC{LinearCombinationEquation,Int64} with 2 stored entries:
  [1, 1]  =  ĥ|b⟩ + -[a|b]|a⟩ + [a|a]|b⟩
  [1, 2]  =  -[a|d]|c⟩ + [a|c]|d⟩
```
"""
function Base.diff(E::EM, o::O) where {EM<:EnergyExpression,O}
    m,n = size(E)
    eqm = spzeros(LinearCombinationEquation, m, n)
    for i = 1:m
        for j = 1:n
            eq = diff(E[i,j], o)
            iszero(eq) && continue
            eqm[i,j] = eq
        end
    end
    eqm
end


text:
using Dates

mutable struct DSA3217 <: AbstractScanivalve
    "Device name"
    devname::String
    "IP adress"
    ipaddr::IPv4
    "TCP/IP port"
    port::Int
    "Configurations paramenters"
    params::Dict{Symbol,Any}
    "Data acquisition buffer"
    buffer::CircMatBuffer{UInt8}
    "Data aquisition task handling"
    task::DAQTask
    "Index of channels"
    chans::Vector{Int}
    "Names of channels"
    channames::Vector{String}
    "Map from channel name to channel index"
    chanidx::Dict{String,Int}
    "Scanivalve configuration"
    conf::DAQConfig
    "Use threading"
    usethread::Bool
end

"Returns the IP address of the device"
ipaddr(dev::DSA3217) = dev.ipaddr

"Returns the port number used for TCP/IP communication"
portnum(dev::DSA3217) = dev.port

"Is DSA3217 acquiring data?"
AbstractDAQs.isreading(dev::DSA3217) = isreading(dev.task)

"How many samples have been read?"
AbstractDAQs.samplesread(dev::DSA3217) = samplesread(dev.task)

"Convert number to string justifying to the right by padding with zeros"
numstring(x::Integer, n=2) = string(10^n+x)[2:end]

function Base.show(io::IO, dev::DSA3217)
    println(io, "Scanivalve DSA3217")
    println(io, "    Dev Name: $(devname(dev))")
    println(io, "    IP: $(string(dev.ipaddr))")
end

"""
`openscani(ipaddr::IPv4, port=23,  timeout=5)`
`openscani(dev::DSA3217,  timeout=5)`
`openscani(fun::Function, ipaddr::IPv4, port=23, timeout=5)`
`openscani(fun::Function, dev::DSA3217, timeout=5)`

Open a TCP/IP connection to scanivalve.

## Arguments

 * `ipaddr` IP Address in `IPv4` format
 * `port` TCP/IP port
 * `timeout` Timeout to wait while trying to connect to scanivalve
 * `fun` Function that is executed after opening the function.

The best way to use this function is to use `do` construct:

```julia
open(dev) do sock
    # Do stuff with the socket
end
```
"""
function openscani(ipaddr::IPv4, port=23,  timeout=5)

    sock = TCPSocket()
    t = Timer(_ -> close(sock), timeout)
    try
        connect(sock, ipaddr, port)
    catch e
        if isa(e, InterruptException)
            throw(InterruptException())
        else
            error("Could not connect to $ipaddr ! Turn on the device or set the right IP address!")
        end
    finally
        close(t)
    end

    return sock
end

openscani(dev::DSA3217,  timeout=5) = openscani(ipaddr(dev), portnum(dev), timeout)

function openscani(fun::Function, ip, port=23, timeout=5)
    io = openscani(ip, port, timeout)
    try
        fun(io)
    catch e
        throw(e)
    finally
        close(io)
    end
end

function openscani(fun::Function, dev::DSA3217, timeout=5)
    io = openscani(ipaddr(dev), portnum(dev), timeout)
    try
        fun(io)
    catch e
        throw(e)
    finally
        close(io)
    end
end

"""
```
DSA3217(devname="Scanivalve", ipaddr="191.30.80.131";
                 timeout=5, buflen=100_000, tag="", sn="", usethread=true)
```

Create an `AbstractDAQ` device that is used to communicate with DSA3217.

## Arguments

 * `devname`: String with name assigned to the specific `DSA3217` device.
 * `ipaddr`: IP address of pressure scanner.
 * `timeout`: timeout to wait while trying to connect to scanivalve.
 * `buflen`: size of buffer used to store data acquired.
 * `tag`: string with device tag.
 * `sn`: string with serial number of the device.
 * `usethread`: `Bool` specifying whether to use thread when acquiring data asynchronously.

If `usethread` is `false`, `@async` is used. Remember that to use threads, `Julia` should be started with `-t N` (N is the number of threads):
```
julia -t 4    # Starts julia with 4 threads.
```

## Input channels

After establishing connection with the DSA3217, all 16 channels are available and named
`P01` through `P16`. This can be changed using the method [`daqaddinput`](@ref).

## Data acquisition configuration

To configure the data acquisition, there are two possibilities:

 - [`daqconfigdev`](@ref) Uses notation and terminology from DSA3217's manual to configure
 - [`daqconfig`](@ref) Uses parameters such as sampling frequency and total acquisition time

## Acquiring data

To do a synchronous data acquisition, method [`daqacquire`](@ref) should be used. In this case, the function blocks until data acquisition is finalized.

For asynchronous data acquisition, the method [`daqstart`](@ref) starts data aquisition. To read the data use method [`daqread`](@ref). This method will wait until data acquisition is over.

Method [`isreading`](@ref) checks whether data acquisition is going on. [`samplesread`](@ref) returns.

[`daqstop`](@ref) interrupts an asynchronous data acquisition task.

## Examples
```jldoctest
julia> using Scanivalve

julia> scani = DSA3217("dpref", "191.30.80.131")
Scanivalve DSA3217
    Dev Name: dpref
    IP: 191.30.80.131

julia> daqaddinput(scani, 1:4, names=["ptot", "pest", "prtot", "prest"])

julia> daqconfigdev(scani, PERIOD=500, AVG=10, FPS=10)

julia> p = daqacquire(scani);

julia> daqstart(scani)
Task (runnable) @0x00007fa61d769e40

julia> p = daqread(scani);

```
"""
function DSA3217(devname="Scanivalve", ipaddr="191.30.80.131";
                 timeout=5, buflen=100_000, tag="", sn="", usethread=true)

    ip = IPv4(ipaddr)
    port = 23

    openscani(ip, port, timeout) do s
        buffer = zeros(Int8, 0, 0)
        println(s, "SET EU 1")
        println(s, "SET AVG 16")
        println(s, "SET PERIOD 500")
        println(s, "SET FPS 1")
        println(s, "SET BIN 1")
        println(s, "SET XSCANTRIG 0")
        println(s, "SET UNITSCAN PA")
        println(s, "SET TIME 1")
    end

    params = Dict{Symbol,Any}(:FPS=>1, :AVG=>16, :PERIOD=>500, :TIME=>1,:XSCANTRIG=>0, :EU=>1, :UNITSCAN=>"PA")

    ipars = Dict{String,Int}("FPS"=>1, "AVG"=>16, "PERIOD"=>500, "TIME"=>1,
                             "EU"=>1, "XSCANTRIG"=>0)
    spars = Dict{String,String}("UNITSCAN"=>"PA")
    fpars = Dict{String,Float64}()

    conf = DAQConfig(ipars, fpars, spars,
                   devname=devname, ip=ipaddr, model="DSA3217", sn=sn,tag=tag)
    task = DAQTask()
    buf = CircMatBuffer{UInt8}(112, buflen)
    chn = "P" .* numstring.(1:16)
    chanidx = Dict{String,Int}()
    for i in 1:16
        chanidx[chn[i]] = i
    end

    return DSA3217(devname, ip, port, params, buf, task, collect(1:16),
                   chn, chanidx, conf, usethread)

end

"""
`daqpacketsize(::Type{DSA3217}, eu=1, time=1)`
"""
function daqpacketsize(::Type{DSA3217}, eu=1, time=1)
    if time==0
        return (eu==0) ? 72 : 104
    else
        return (eu==0) ? 80 : 112
    end
end

"""
`daqparam(dev, param)`

Returns de value configuration parameter `param` (s `Symbol`)
"""
daqparam(dev, param) = dev.params[param]

"""
`framesize(dev::DSA3217)`

Return the size of a data acquisition frame. A frame contains a data from all channels.
"""
framesize(dev::DSA3217) = daqpacketsize(DSA3217, daqparam(dev,:EU), daqparam(dev,:TIME))

"""
`deltat(dev)`

Calculate the sampling interval calculated from data acquisition paramters:

 * `PERIOD` Time in μs between data acquisition in *each channel*
 * `AVG` Number of frames that are read before averaging.

"""
deltat(dev) = 16*daqparam(dev,:PERIOD)*1e-6 * daqparam(dev,:AVG)

"""
`stopscan(dev::DSA3217)`

Change the `stop` field of `task::DAQTask` field of `DSA3217` object to `true`.
This will flag that data acquisition should stop.

"""
function stopscan(dev::DSA3217)
    tsk = dev.task

    if isreading(tsk)
        tsk.stop = true
    end
end

"""
`scan!(dev::DSA3217)`

Execute data acquisition from `DSA3217` devices. Present configuration should be used.

"""
function scan!(dev::DSA3217)

    # Check if data acquisition is already going on.
    tsk = dev.task
    isreading(tsk) && error("DSA is already reading!")

    # Clear stuff and adjust buffer.
    cleartask!(tsk)

    buf = dev.buffer
    fps = daqparam(dev, :FPS)
    # For long data acquisition, buffer might have to be increased.
    if fps > capacity(buf)
        resize!(buf, fps)
    end
    empty!(buf)

    fsize = framesize(dev)
    nfrs = length(buf)

    if fps == 0 # Read continously
        fps1 = typemax(Int32)
    else
        fps1 = fps
    end

    # Calculate time interval between frames
    δt = deltat(dev) # Time per frame
    stopped = false

    # Register time and date of data acquisition
    tsk.time = now()

    exthrown = false # No exception thrown!

    # Open socket, send SCAN command and acquire data.
    openscani(dev) do sock
        println(sock, "SCAN") # Flags scanivalve to start data acquisition
        tsk.isreading = true
        # Initial data acquisition instant
        t0 = time_ns()
        # Get memory to store first frame
        b = nextbuffer(buf)
        readbytes!(sock, b, fsize)
        # Get the time before second frame. Remember that
        # there might be some latency. But for scanivalve this is not very important
        # since the calculated sampling time (`deltat`) is fairly accurate.
        t1 = time_ns()
        settiming!(tsk, t0, t1, 1)
        tsk.nread += 1
        for i in 2:fps1
            # Check if data acquisition should stop ([`daqstop`]($ref) or
            # [`stopscan`](@ref) commands).
            try
                if tsk.stop
                    stopped = true
                    println(sock, "STOP")
                    sleep(0.5)
                    break
                end
                # Read next frame
                b = nextbuffer(buf)
                readbytes!(sock, b, fsize)
                tn = time_ns()
                tsk.nread += 1
                settiming!(tsk, t1, tn, i-1)
            catch e
                if isa(e, InterruptException)
                    # Ctrl-C captured!
                    # We want to stop the data acquisition safely and then rethwrow it!
                    tsk.stop = true
                    exthrown = true
                else
                    throw(e)
                end
            end

        end

        tsk.isreading = false

        # If we stopped, try to read another frame.
        # there might be something in the buffer. But we will ignore it
        if stopped
            delay = max(2δt, 1.0) # Let's give it some time.
            ev = Base.Event()
            Timer(_ -> begin
                      timeout=true
                      notify(ev)
                  end, delay)
            @async begin
                read(sock, fsize)
                read!(sock, buffer(tsk, idx))
                notify(ev)
            end
            wait(ev)
            if exthrown
                # Data acquisition stopped because of a Ctrl-C
                # We handled stopping DAQ smoothly and now we will rethrow
                throw(InterruptException())
            end

        end
    end
    return
end

"""
Valid units accepted by DSA3217.
"""
const validunits = ["ATM", "BAR", "CMHG", "CMH2O", "DECIBAR", "FTH2O", "GCM2",
                    "INHG", "INH2O", "KPA", "KGCM2", "KGM2", "KIPIN2", "KNM2",
                    "MH2O", "MMHG", "MPA", "NCM2", "MBAR", "OZFT2", "OZIN2",
                    "PA", "PSF", "NM2", "PSI", "TORR"]

"""
`daqconfig(dev::DSA3217; kw...)`

Generic configuration of data acquisition. Here we specify how many samples will be read
and how fast it should be read.

## Arguments

 * `dt`: period of time in seconds between samples
 * `rate`: sampling rate in Hz between samples
 * `nsamples`: number of samples that should be read
 * `time`: time in seconds that data should be acquired.

If `dt` is specified,
```
rate = 1/dt
```

Either the number of samples should be configured or the acquisition time:
```
time = nsamples * dt
```

Samples can be averaged. The `avg` parameter specifies how many samples should be averaged.
Attention: the sampling time (`dt`) is includes the averaging.

Trigger can also be specified.

 * `trigger == 0` Internal trigger
 * `trigger == 1` Use external trigger.

**Important** If a parameter is not used, previous configuration will be used.
Remember to set the value of `avg`. The best option is to specify all the options.

## Example
```jldoctest
julia> daqconfig(scani, rate=100, nsamples=100, avg=1)

julia> p = daqacquire(scani);

julia> samplingrate(p)
100.0

julia> daqconfig(scani, dt=0.1, time=1, avg=1)

julia> p = daqacquire(scani);

julia> samplingrate(p)
10.0

```
"""
function AbstractDAQs.daqconfig(dev::DSA3217; kw...)
    cmds = String[]

    if haskey(kw, :avg)
        avg = round(Int, kw[:avg])
        daqconfigdev(dev, AVG=avg)
    else
        avg = dev.params[:AVG]
    end

    if haskey(kw, :rate) && haskey(kw, :dt)
        error("Parameters `rate` and `dt` can not be specified simultaneously!")
    elseif haskey(kw, :rate) || haskey(kw, :dt)
        if haskey(kw, :rate)
            rate = kw[:rate]
            period = round(Int, 1.0 / (rate*16e-6*avg))
        else
            dt = kw[:dt]
            period = round(Int, dt / (16e-6*avg))
        end
        daqconfigdev(dev, PERIOD=period)
    else
        period = dev.params[:PERIOD]
    end

    if haskey(kw, :nsamples) && haskey(kw, :time)
        error("Parameters `nsamples` and `time` can not be specified simultaneously!")
    elseif haskey(kw, :nsamples) || haskey(kw, :time)
        if haskey(kw, :nsamples)
            nsamples = kw[:nsamples]
        else
            tt = kw[:time]
            dt = period * 16e-6 * avg
            nsamples = round(Int, tt / dt)
        end
        daqconfigdev(dev, FPS=nsamples)
    else
        nsamples = dev.params[:FPS]
    end

    if haskey(kw, :trigger)
        trigger = kw[:trigger]
        daqconfigdev(dev, XSCANTRIG=trigger)
    end

end

"""
`daqstart(dev::DSA3217)`

Start an asynchrounous data acquisition.

If `dev.usethread == true`, a thread will be started with `@spawn`. Othwerwise, `@async` is used. This should be specified when creating the DSA3217 object ([`DSA3217`](@ref)).

## Example

```jldoctest
julia> daqconfig(scani, rate=100, time=10, avg=1) # 10s, 100 Hz

julia> @time begin
             daqstart(scani) # Start data
             sleep(5) # Wait 5s
             daqstop(scani) # Interrupt daq
             p = daqread(scani) # Read data
             size(p.data)
end
  6.531171 seconds (40.08 k allocations: 2.049 MiB, 1.86% compilation time)
(4, 487)

```
"""
function AbstractDAQs.daqstart(dev::DSA3217)
    if isreading(dev)
        error("Scanivalve already reading!")
    end
    if dev.usethread
        tsk = Threads.@spawn scan!(dev)
    else
        tsk = @async scan!(dev)
    end
    dev.task.task = tsk
    return tsk
end

"""
`daqaddinput(dev::DSA3217, chans=1:16; names="P")`

Specify input channels.

The DSA3217 always reads all 16 channels. Choosing a smaller
number of channels will not make data acquisition faster. But it may
be convenient to forget some of them.

## Arguments

 * `chans`: an abstract vector with the index of channels that should be used.
 * `names`: names of channels (see next)

## Channel names

The `names` keyword argument is used to specify the names of the channels used.

If a single string is specified, the channel names will be this string followed by
the index (padded by 1 zero at most). See example.

Individual channels can be named but in this case, every single channel should be named.

## See also

[`DSA3217`](@ref), [`daqconfig`](@ref), [`daqconfigdev`](@ref), [`daqacquire`](@ref),
[`daqstart`](@ref), [`daqread`](@ref)

## Example
```jldoctest
julia> daqaddinput(scani, 8:12, names="press")

julia> daqchannels(scani)
5-element Vector{String}:
 "press08"
 "press09"
 "press10"
 "press11"
 "press12"

julia> daqaddinput(scani, 1:4, names=["ptot", "pest", "prtot", "prest"])

julia> daqchannels(scani)
4-element Vector{String}:
 "ptot"
 "pest"
 "prtot"
 "prest"
```
"""
function AbstractDAQs.daqaddinput(dev::DSA3217, chans=1:16; names="P")

    cmin, cmax = extrema(chans)
    if cmin < 1 || cmax > 16
        throw(ArgumentError("Only channels 1-16 are available to DSA3217"))
    end

    if isa(names, AbstractString) || isa(names, Symbol)
        chn = string(names) .* numstring.(chans)
    elseif length(names) == length(chans)
        chn = string.(names)
    else
        throw(ArgumentError("Argument `names` should have length 1 or the length of `chans`"))
    end

    dev.chans = collect(chans)
    dev.channames = chn
    n = length(chans)
    chanidx = Dict{String,Int}()
    for i in 1:n
        chanidx[chn[i]] = i
    end
    dev.chanidx = chanidx

    return
end

"""
`daqstop(dev::DSA3217)`

Interrupt asynchronous data acquisition.

## Example

```jldoctest
julia> daqconfig(scani, rate=100, time=10, avg=1) # 10s, 100 Hz

julia> @time begin
             daqstart(scani) # Start data
             sleep(5) # Wait 5s
             daqstop(scani) # Interrupt daq
             p = daqread(scani) # Read data
             size(p.data)
end
  6.531171 seconds (40.08 k allocations: 2.049 MiB, 1.86% compilation time)
(4, 487)

```

"""
function AbstractDAQs.daqstop(dev::DSA3217)

    tsk = dev.task.task
    if !istaskdone(tsk) && istaskstarted(tsk)
        stopscan(dev)
        wait(tsk)
    end
    dev.task.stop = false
    dev.task.isreading = false

end

"""
`readpressure(dev::DSA3217)`

Reads pressure, sampling rate and acquisition start time from buffer.

"""
function readpressure(dev::DSA3217)
    isreading(dev) && error("Scanivalve still acquiring data!")

    tsk = dev.task
    buf = dev.buffer
    nsamples = length(buf)
    δt = getdaqtime(dev, nsamples)

    if daqparam(dev, :EU) > 0
        press = read_eu_press(buf, dev.chans)
    else
        error("Reading data without engineering units (EU=1) not implemented yet!")
    end

    return press, 1.0/δt, tsk.time
end

export meastime, samplingrate, measdata, measinfo

"""
`daqread(dev::DSA3217)`

Wait until data acquisition ends and read data.

If continous data acqisition is used, this function will stop data acquisition and return data.

Data acquisition starts with [`daqstart`](@ref) command.

To just to get the data already acquired, there will be function [`daqpeek`](@ref) that has
not been implemented yet.

## Return value

This function returns a [`AbstractDAQs.MeasData`](@ref) object that contains the data as
well as sampling rate and other meta data.

## Usage
See example for [`daqstart`](@ref) command.

"""
function AbstractDAQs.daqread(dev::DSA3217)

    # Check if the reading is continous
    if daqparam(dev, :FPS) == 0
        # Stop reading!
        daqstop(dev)
        sleep(0.1)
    end
    # Wait for task to end
    wait(dev.task.task)

    # Get the data:
    P, fs, t = readpressure(dev)
    unit = dev.params[:UNITSCAN]
    return MeasData{Matrix{Float32},String}(devname(dev), devtype(dev),
                                             t, fs, P, unit, dev.chanidx)

end

"""
`daqacquire(dev::DSA3217)`

Execute a **s**ynchronous data acquisition.

## Return value

This function returns a [`AbstractDAQs.MeasData`](@ref) object that contains the data as
well as sampling rate and other meta data.

## Usage
See example for [`DSA3217`](@ref) command.

"""
function AbstractDAQs.daqacquire(dev::DSA3217)
    scan!(dev)
    P, fs, t = readpressure(dev)
    unit = dev.params[:UNITSCAN]
    return MeasData{Matrix{Float32},String}(devname(dev), devtype(dev),
                                             t, fs, P, unit, dev.chanidx)
end

"""
`read_eu_press(buf, chans)`

Reads pressure in engineering units (EU) from the data acquisition buffer.

"""
function read_eu_press(buf, chans)

    nt = length(buf)
    nch = length(chans)
    P = zeros(Float32, nch, nt)

    for i in 1:nt
        p1 = reinterpret(Float32, view(buf[i], 9:72))
        for k in 1:nch
            P[k,i] = p1[chans[k]]
        end
    end
    return P
end

"""
`getdaqtime(dev, nfr)`

Returns the mean sampling time. If no time is specified (paramter `TIME==0`),
just calculate from acquisition parameters `PERIOD` and `AVG`.

Otherwises, get the value from the frames.
Reads pressure in engineering units (EU) from the data acquisition buffer.

"""
function getdaqtime(dev, nfr)
    buf = dev.buffer
    b = buf[1,1] # Identify the packet

    # Sampling time from scan configuration and fallback value
    δt₀ = daqparam(dev, :AVG) * daqparam(dev, :PERIOD)*1e-6 * 16
    btype = daqparam(dev, :TIME)
    if b == 0x04 || b == 0x05 || nfr==1 # No time
        # Calculate from scan configuration
        return δt₀
    else
        if b == 0x07  # EU with time
            bt1 = reinterpret(Int32, buf[105:108,1])[1]
            bt2 = reinterpret(Int32, buf[105:108,nfr])[1]
        elseif b == 0x06
            bt1 = reinterpret(Int32, buf[73:76,1])[1]
            bt2 = reinterpret(Int32, buf[73:76,nfr])[1]
        else # This shouldn't happen! Fallback
            return δt₀
        end
        ft = (btype==1) ? 1e6 : 1e3
        return  (bt2 - bt1) / (ft*(nfr-1))
    end

end

function scan2press(dev::DSA3217, buf, info)

    nframes = info.FPS
    press = zeros(Float32, nframes, 16)

    if info.EU > 0
        for i in 1:nframes
            press[i,:] = reinterpret(Float32, view(buf, 9:72,i))
        end
    else
        for i in 1:nframes
            press[i,:] = Float32.(reinterpret(Int16, view(buf, 9:40, i)))
        end
    end
    return press

end

import Base.close

"""
`close(scani::DSA3217)`

Return scanivalve to a 'sane' configuration.
"""
function close(scani::DSA3217)
    openscani(scani) do s
        println(s, "SET EU 1")
        println(s, "SET AVG 100")
        println(s, "SET PERIOD 500")
        println(s, "SET FPS 1")
        println(s, "SET BIN 0")
        println(s, "SET XSCANTRIG 0")
        println(s, "SET UNITSCAN PA")
        println(s, "SET TIME 0")
    end

end

"Number of data acquisition channels"
numchans(scani::DSA3217) = 16
"Number of data acquisition channels"
AbstractDAQs.numchannels(scani::DSA3217) = length(scani.chans)
"Name of data acquisition channels"
AbstractDAQs.daqchannels(scani::DSA3217) = scani.channames

#socket(scani) = scani.socket

function scanpacksize(dev::DSA3217, TIME=0, EU=1)
    if TIME==0
        if EU==1
            return 104
        else
            return 72
        end
    else
        if EU==1
            return 112
        else
            return 80
        end
    end

end

"Paramaters valid for `daqconfigdev`"
const validparameters = [:FPS, :PERIOD, :AVG, :TIME, :EU, :UNITSCAN, :XSCANTRIG]

"""
`daqconfigdev(dev::DSA3217; kw...)`

Configure data acquisition with parameters described in the manual.

The following parameters are available:

 * `FPS`: Integer specifying the number of frames (samples) to acquire
 * `PERIOD`: time spent in each pressure sensor in μs
 * `AVG`: Number of frames to average before outputing a mean frame
 * `TIME`: 0 - No time in frame, 1 or 2 (ms or μs)
 * `EU`: 1 if engineering units should be used (not implemented for non EU)
 * `UNITSCAN`: String with pressure units that should be used. Check `validunits` vector.
 * `XSCANTRIG`: 0 for internal trigger or 1 for external trigger.

See DSA3217 manual for further information.

## Examples

```jldoctest
julia> daqconfigdev(scani, PERIOD=500, AVG=10, FPS=10)

julia> p = daqacquire(scani);

julia> samplingrate(p)
12.5

julia> size(p.data)
(16, 10)
```

"""
function AbstractDAQs.daqconfigdev(dev::DSA3217; kw...)

    k = keys(kw)
    cmds = String[]
    if :FPS ∈ k
        fps = kw[:FPS]
        if 0 ≤ fps < 1_000_000
            push!(cmds, "SET FPS $fps")
            dev.params[:FPS] = fps
        else
            throw(DomainError(fps, "FPS outside range (0 - 1000000)"))
        end
    end

    if :PERIOD ∈ k
        period = kw[:PERIOD]
        if 125 ≤ period ≤ 65535
            push!(cmds, "SET PERIOD $period")
            dev.params[:PERIOD] = period
        else
            throw(DomainError(period, "PERIOD outside range (126 - 65535)"))
        end
    end

    if :AVG ∈ k
        avg = kw[:AVG]
        if 1 ≤ avg ≤ 240
            push!(cmds, "SET AVG $avg")
            dev.params[:AVG] = avg
        else
            throw(DomainError(avg, "AVG outside range (1 - 240)"))
        end
    end

    if :TIME ∈ k
        tt = kw[:TIME]
        if 0 ≤ tt ≤ 2
            push!(cmds, "SET TIME $tt")
            dev.params[:TIME] = tt
        else
            throw(DomainError(tt, "TIME outside range (0-2)"))
        end
    end

    if :EU ∈ k
        eu = kw[:EU]
        if 0 ≤ eu ≤ 1
            push!(cmds, "SET EU $eu")
            dev.params[:EU] = eu
        else
            throw(DomainError(eu, "EU outside range (0 or 1)"))
        end
    end

    if :UNITSCAN ∈ k # User should check manually if the correct unit was used
        unitscan = kw[:UNITSCAN]
        if unitscan ∈ validunits
            push!(cmds, "SET UNITSCAN $unitscan")
            dev.params[:UNITSCAN] = unitscan
        else
            throw(DomainError(unitscan, "Invalid unit!"))
        end
    end

    if :XSCANTRIG ∈ k
        xscantrig = kw[:XSCANTRIG]
        if 0 ≤ eu ≤ 1
            push!(cmds, "SET SCANTRIG $xscantrig")
            dev.params[:XSCANTRIG] = xscantrig
        else
            throw(DomainError(xscantrig, "XSCANTRIG should be either 0 or 1"))
        end
    end

    # Send commands to Scanivalve
    openscani(dev) do sock
        for c in cmds
            println(sock, c)
        end

        # Update conf field:
        updateconf!(dev::DSA3217)
    end
end

"""
`updateconf!(dev::DSA3217)`

Update configuration from `params` field of `DSA3217` object.
"""
function updateconf!(dev::DSA3217)
    p = dev.params
    ipars = dev.conf.ipars
    spars = dev.conf.spars

    ipars["FPS"] = p[:FPS]
    ipars["AVG"] = p[:AVG]
    ipars["PERIOD"] = p[:PERIOD]
    ipars["TIME"] = p[:TIME]
    ipars["XSCANTRIG"] = p[:XSCANTRIG]
    ipars["EU"] = p[:EU]

    spars["UNITSCAN"] = p[:UNITSCAN]

    return

end

"""
`daqzero(dev::DSA3217; time=15)`

Execute a hardware zero. The time parameter specifies how long the function should sleep
before returning.

"""
function AbstractDAQs.daqzero(dev::DSA3217; time=15)
    openscani(dev) do io
        println(io, "CALZ")
        sleep(time)
    end
end

"""
`readmanylines(dev, cmd, delay=0.5)`

Several commands return data as a sequence of text lines. Read as many lines as possible
waiting for `delay` seconds before ending any attempt to read.

"""
function readmanylines(dev, cmd, delay=0.5)

    lst = String[]
    openscani(dev) do socket
        println(socket, cmd)
        timeout = false
        while !timeout
            ev = Base.Event()
            Timer(_ -> begin
                      timeout=true
                      notify(ev)
                  end, delay)
            @async begin
                line = readline(socket)
                push!(lst, line)
                notify(ev)
            end
            wait(ev)
        end
    end

    return lst

end

function readstatus(dev::DSA3217)
    openscani(dev) do s
        println(s, "STATUS")
        string(Char.(read(s, 180)[81:100])...)
    end
end

"""
`readnlines(s, n=1)`

Read `n` lines from open socket `s`.
"""
function readnlines(s, n=1)
    lst = String[]
    for i in 1:n
        push!(lst, readline(s))
    end
    return lst
end

"""
`parse_dsa_set(lst)`

Usually configuration parameters from scanivalve are read as a collection of lines with
the following format (each line):
```
SET PARAMETER VALUE
```

This function will create a `Dict{String,String}` from the list of lines with the format
`d["PARAMETER"] = "VALUE"`
"""
function parse_dsa_set(lst)
    p = Dict{String,String}()
    for l in lst
        s = split(l, " ")
        p[s[2]] = s[3]
    end
    return p
end

"""
`listany(dev::DSA3217, cmd, nparams)`

List any configuration parameter by sending command LIST cmd to scanivalve.
"""
function listany(dev::DSA3217, cmd, nparams)
    openscani(dev) do sock
        println(sock, "LIST $cmd\r")
        readnlines(sock, nparams)
    end
end

"""
`listanydict(dev::DSA3217, cmd, nparams)`

List any configuration parameter by sending command LIST cmd to scanivalve.
Build a dictionary with the output.
"""
listanydict(scani::AbstractScanivalve, cmd, nparams) = parse_dsa_set(listany(scani, cmd, nparams))

function listanyval(scani::AbstractScanivalve,
                    cmd, nparams, ::Type{T}=Int) where {T<:Real}
    lst = listany(scani, cmd, nparams)
    nlst = length(lst)
    vals = zeros(T, nlst)

    for i in 1:nlst
        s = split(lst[i], " ")
        vals[i] = parse(T, s[3])
    end

    return vals
end

"Execute command LIST SCAN"
listscan(scani::DSA3217) = listanydict(scani, "S", 14)
"Execute command LIST I"
listident(scani::DSA3217) = listanydict(scani, "I", 4)
"Execute command LIST Z"
listzero(scani::DSA3217) = listanyval(scani, "Z", numchans(scani), Int)
"Execute command LIST O"
listoffset(scani::DSA3217) = listanyval(scani, "O", numchans(scani), Float64)
"Execute command LIST D"
listdelta(scani::DSA3217) = listanyval(scani, "O", numchans(scani), Float64)
"Execute command LIST G"
listgain(scani::DSA3217) = listanyval(scani, "G", numchans(scani), Float64)

"""
`setparam(scani, param, val)`

Set the value of a configuration paramenter on the scanivalve.
"""
setparam(scani, param, val) =
    openscani(scani) do sock
        println(sock, "SET $param $val\r")
    end


text:
const LINE_BREAK = 120

mutable struct GAMSTranslateStream
   io::IOStream
   n_line::Int
end

GAMSTranslateStream(io::IOStream) = GAMSTranslateStream(io, 0)

function write(
   io::GAMSTranslateStream,
   str::String
)
   n = length(str)
   if io.n_line + n > LINE_BREAK
      idx = n+1
      while true
         idx = findprev(isequal(' '), str, idx-1)
         if idx == nothing || io.n_line + idx <= LINE_BREAK
            break
         end
      end
      if idx == nothing
         Base.write(io.io, "\n  ")
         io.n_line = Base.write(io.io, str)
      else
         Base.write(io.io, str[1:idx-1])
         Base.write(io.io, "\n ")
         Base.write(io.io, str[idx+1:end])
         io.n_line = n-idx
      end
   else
      Base.write(io.io, str)
      io.n_line += n
   end
end

function writeln(
   io::GAMSTranslateStream,
   str::String
)
   write(io, str)
   Base.write(io.io, "\n")
   io.n_line = 0
end

function dbl2str(
   value::Number
)
   @sprintf("%.16g", value)
end

function translate_header(
   io::GAMSTranslateStream
)
   writeln(io, "*\n* GAMS Model generated by GAMS.jl\n*\n")
   writeln(io, "\$offlisting")
end

function translate_defsets(
   io::GAMSTranslateStream,
   model::Optimizer
)
   n = length(model.sos1_constraints) + length(model.sos2_constraints)
   if n == 0
      return
   elseif n == 1
      writeln(io, "Set")
   else
      writeln(io, "Sets")
   end
   write(io, "  ")

   first = true
   for (i, con) = enumerate(model.sos1_constraints)
      translate_defsets(io, model, i, con.func, con.set, first=first)
      first = false
   end
   for (i, con) = enumerate(model.sos2_constraints)
      translate_defsets(io, model, i, con.func, con.set, first=first)
      first = false
   end
   writeln(io, "\n");
end

function translate_defsets(
   io::GAMSTranslateStream,
   model::Optimizer,
   idx::Int,
   func::MOI.VectorOfVariables,
   set::Union{MOI.SOS1{Float64}, MOI.SOS2{Float64}};
   first::Bool=true
)
   if ! first
      write(io, ", ")
   end

   if typeof(set) == MOI.SOS1{Float64}
      write(io, "s1s$idx / ")
   elseif typeof(set) == MOI.SOS2{Float64}
      write(io, "s2s$idx / ")
   end

   first_elem = true
   for vi in func.variables
      if ! first_elem
         write(io, ", ")
      end
      write(io, "$(vi.value)")
      first_elem = false
   end

   write(io, " /")
end

function translate_defvars(
   io::GAMSTranslateStream,
   model::Optimizer
)
   translate_defvars(io, model, nothing)
   if model.n_binary > 0
      translate_defvars(io, model, VARTYPE_BINARY)
   end
   if model.n_integer > 0
      translate_defvars(io, model, VARTYPE_INTEGER)
   end
   if model.n_semicont > 0
      translate_defvars(io, model, VARTYPE_SEMICONT)
   end
   if model.n_semiint > 0
      translate_defvars(io, model, VARTYPE_SEMIINT)
   end
   if length(model.sos1_constraints) > 0
      translate_defvars(io, model, VARTYPE_SOS1)
   end
   if length(model.sos2_constraints) > 0
      translate_defvars(io, model, VARTYPE_SOS2)
   end
end

function translate_defvars(
   io::GAMSTranslateStream,
   model::Optimizer,
   filter::Union{GAMSVarType, Nothing}
)
   if filter == VARTYPE_SOS1
      n = length(model.sos1_constraints)
   elseif filter == VARTYPE_SOS2
      n = length(model.sos2_constraints)
   else
      n = length(model.variable_info)
   end

   if filter == VARTYPE_BINARY
      write(io, "Binary ")
   elseif filter == VARTYPE_INTEGER
      write(io, "Integer ")
   elseif filter == VARTYPE_SEMICONT
      write(io, "SemiCont ")
   elseif filter == VARTYPE_SEMIINT
      write(io, "SemiInt ")
   elseif filter == VARTYPE_SOS1
      write(io, "SOS1 ")
   elseif filter == VARTYPE_SOS2
      write(io, "SOS2 ")
   end
   if n == 1
      writeln(io, "Variable")
   else
      writeln(io, "Variables")
   end
   write(io, "  ")

   first = true
   if model.objvar && (isnothing(filter) || filter == VARTYPE_FREE)
      if n == 0
         write(io, "objvar;")
         return
      end
      write(io, "objvar")
      first = false
   end

   for i in 1:n
      if ! isnothing(filter) && model.variable_info[i].type != filter
         continue
      end
      if ! first
         write(io, ", ")
      end
      translate_variable(io, model, i)
      first = false
   end

   # add sos1 variables
   if isnothing(filter) || filter == VARTYPE_SOS1
      for i in 1:length(model.sos1_constraints)
         if ! first
            write(io, ", ")
         end
         write(io, "s1x$i(s1s$(i))")
         first = false
      end
   end

   # add sos2 variables
   if isnothing(filter) || filter == VARTYPE_SOS2
      for i in 1:length(model.sos2_constraints)
         if ! first
            write(io, ", ")
         end
         write(io, "s2x$i(s2s$(i))")
         first = false
      end
   end
   writeln(io, ";\n")
end

function translate_defequs(
   io::GAMSTranslateStream,
   model::Optimizer
)

   m = model.m + length(model.sos1_constraints) + length(model.sos2_constraints) + length(model.complementarity_constraints)

   if m == 0 && ! model.objvar
      return
   end

   if m == 1
      writeln(io, "Equation")
   else m > 1
      writeln(io, "Equations")
   end
   write(io, "  ")

   if model.objvar && m == 0
      writeln(io, "obj;\n")
      return
   end

   first = true
   if model.objvar
      write(io, "obj")
      first = false
   end

   for i in 1:model.m
      if ! first
         write(io, ", ")
      end
      write(io, "eq$i")
      first = false
   end

   # add sos1 constraints
   for i in 1:length(model.sos1_constraints)
      if ! first
         write(io, ", ")
      end
      write(io, "s1eq$(i)(s1s$(i))")
      first = false
   end

   # add sos2 constraints
   for i in 1:length(model.sos2_constraints)
      if ! first
         write(io, ", ")
      end
      write(io, "s2eq$(i)(s2s$(i))")
      first = false
   end

   # add complementarity constraints
   for (i, comp) in enumerate(model.complementarity_constraints)
      for j in 1:comp.set.dimension
         if ! first
            write(io, ", ")
         end
         write(io, "eq$(i)_$(j)")
      end
      first = false
   end

   writeln(io, ";\n")
end

function translate_coefficient(
   io::GAMSTranslateStream,
   coef::Float64;
   first::Bool=false
)
   if coef < 0.0
      if first && coef == -1.0
         write(io, "-")
      elseif first
         write(io, "-" * dbl2str(-coef) * " * ")
      elseif coef == -1.0
         write(io, " - ")
      else
         write(io, " - " * dbl2str(-coef) * " * ")
      end
   elseif coef > 0.0
      if first && coef == 1.0
      elseif first
         write(io, dbl2str(coef) * " * ")
      elseif coef == 1.0
         write(io, " + ")
      else
         write(io, " + " * dbl2str(coef) * " * ")
      end
   end
   return
end

function translate_variable(
   model::Optimizer,
   idx::Int
)
   if model.variable_info[idx].type == VARTYPE_FREE
      return "x$idx"
   elseif model.variable_info[idx].type == VARTYPE_BINARY
      return "b$idx"
   elseif model.variable_info[idx].type == VARTYPE_INTEGER
      return "i$idx"
   elseif model.variable_info[idx].type == VARTYPE_SEMICONT
      return "sc$idx"
   elseif model.variable_info[idx].type == VARTYPE_SEMIINT
      return "si$idx"
   end
end

function translate_variable(
   io::GAMSTranslateStream,
   model::Optimizer,
   idx::Int
)
   write(io, translate_variable(model, idx))
end

function translate_function(
   io::GAMSTranslateStream,
   model::Optimizer,
   func::MOI.SingleVariable
)
   translate_variable(io, model, func.variable.value)
end

function translate_function(
   io::GAMSTranslateStream,
   model::Optimizer,
   terms::Vector{MOI.ScalarAffineTerm{Float64}}
)
   if length(terms) == 0
      write(io, "0.0")
      return
   end

   first = true
   for (j, term) in enumerate(terms)
      if term.coefficient == 0.0
         continue
      end
      translate_coefficient(io, term.coefficient, first=first)
      translate_variable(io, model, term.variable_index.value)
      first = false
   end
end

function translate_function(
   io::GAMSTranslateStream,
   model::Optimizer,
   func::MOI.ScalarAffineFunction{Float64}
)
   translate_function(io, model, func.terms)
   if func.constant < 0.0
      write(io, " - " * dbl2str(-func.constant))
   elseif func.constant > 0.0
      write(io, " + " * dbl2str(func.constant))
   end
end

function translate_function(
   io::GAMSTranslateStream,
   model::Optimizer,
   terms::Vector{MOI.VectorAffineTerm{Float64}}
)
   sterms = [term.scalar_term for term in terms]
   translate_function(io, model, sterms)
end

function translate_function(
   io::GAMSTranslateStream,
   model::Optimizer,
   func::MOI.ScalarQuadraticFunction{Float64}
)
   naff = length(func.affine_terms)

   if naff + length(func.quadratic_terms) == 0
      write(io, "0.0")
      return
   end

   first = true
   if naff > 0
      translate_function(io, model, func.affine_terms)
      first = false
   end

   for (j, term) in enumerate(func.quadratic_terms)
      if term.coefficient == 0.0
         continue
      end
      idx1 = term.variable_index_1
      idx2 = term.variable_index_2
      if idx1 == idx2
         translate_coefficient(io, term.coefficient / 2.0, first=first)
      else
         translate_coefficient(io, term.coefficient, first=first)
      end
      translate_variable(io, model, idx1.value)
      write(io, " * ")
      translate_variable(io, model, idx2.value)
      first = false
   end

   if func.constant < 0.0
      write(io, " - " * dbl2str(-func.constant))
   elseif func.constant > 0.0
      write(io, " + " * dbl2str(func.constant))
   end
end

function translate_function(
   io::GAMSTranslateStream,
   model::Optimizer,
   func::Expr;
   is_parenthesis::Bool = false
)
   if length(func.args) == 0
      write(io, "0.0")
      return
   end

   op = func.args[1]

   if op in (:+, :-)
      @assert length(func.args) >= 2
      if ! is_parenthesis
         write(io, "(")
      end
      if length(func.args) == 2
         write(io, "(" * string(op))
         translate_function(io, model, func.args[2], is_parenthesis=false)
         write(io, ")")
      else
         translate_function(io, model, func.args[2], is_parenthesis=true)
         for i in 3:length(func.args)
            write(io, " " * string(op) * " ")
            translate_function(io, model, func.args[i], is_parenthesis=(op==:+))
         end
      end
      if ! is_parenthesis
         write(io, ")")
      end

   elseif op in (:*, :/)
      @assert length(func.args) >= 3
      if ! is_parenthesis
         write(io, "(")
      end
      translate_function(io, model, func.args[2], is_parenthesis=false)
      for i in 3:length(func.args)
         write(io, " " * string(op) * " ")
         translate_function(io, model, func.args[i], is_parenthesis=false)
      end
      if ! is_parenthesis
         write(io, ")")
      end

   elseif op == :^
      @assert length(func.args) == 3
      if func.args[3] == 1.0
         translate_function(io, model, func.args[2])
      elseif func.args[3] == 2.0
         write(io, "sqr(")
         translate_function(io, model, func.args[2], is_parenthesis=true)
         write(io, ")")
      elseif func.args[3] isa Real && func.args[3] == round(func.args[3])
         write(io, "power(")
         translate_function(io, model, func.args[2], is_parenthesis=true)
         write(io, ", ")
         translate_function(io, model, func.args[3], is_parenthesis=true)
         write(io, ")")
      else
         translate_function(io, model, func.args[2], is_parenthesis=false)
         write(io, "**")
         translate_function(io, model, func.args[3], is_parenthesis=false)
      end

   elseif op in (:sqrt, :log, :log10, :log2, :exp, :sin, :sinh, :cos, :cosh, :tan, :tanh, :abs, :sign)
      @assert length(func.args) == 2
      write(io, string(op) * "(")
      translate_function(io, model, func.args[2], is_parenthesis=true)
      write(io, ")")

   elseif op in (:acos,)
      @assert length(func.args) == 2
      write(io, "arccos(")
      translate_function(io, model, func.args[2], is_parenthesis=true)
      write(io, ")")

   elseif op in (:asin,)
      @assert length(func.args) == 2
      write(io, "arcsin(")
      translate_function(io, model, func.args[2], is_parenthesis=true)
      write(io, ")")

   elseif op in (:atan,)
      @assert length(func.args) == 2
      write(io, "arctan(")
      translate_function(io, model, func.args[2], is_parenthesis=true)
      write(io, ")")

   elseif op in (:max, :min, :mod)
      @assert length(func.args) >= 2
      write(io, string(op) * "(")
      translate_function(io, model, func.args[2], is_parenthesis=true)
      for i in 3:length(func.args)
         write(io, ", ")
         translate_function(io, model, func.args[i], is_parenthesis=true)
      end
      write(io, ")")

   elseif typeof(op) == Symbol && func.args[2] isa MOI.VariableIndex
      translate_variable(io, model, func.args[2].value)
   else
      error("Unrecognized operation ($op)")
   end
end

function translate_function(
   io::GAMSTranslateStream,
   model::Optimizer,
   func::Float64;
   is_parenthesis::Bool = false
)
   if is_parenthesis && func > 0.0
      write(io, dbl2str(func))
   else
      write(io, "(" * dbl2str(func) * ")")
   end
end

function translate_function(
   io::GAMSTranslateStream,
   model::Optimizer,
   func::Int;
   is_parenthesis::Bool = false
)
   write(io, "$func")
end

function translate_objective(
   io::GAMSTranslateStream,
   model::Optimizer
)
   # do nothing if we don't need objective variable
   if model.model_type == GAMS.MODEL_TYPE_MCP || model.model_type == GAMS.MODEL_TYPE_CNS
      return
   end
   if ! model.objvar
      if ! (typeof(model.objective) == MOI.SingleVariable)
         error("GAMS needs obj variable")
      end
      return
   end

   if model.sense == MOI.MIN_SENSE
      write(io, "obj.. objvar =G= ")
   elseif model.sense == MOI.MAX_SENSE
      write(io, "obj.. objvar =L= ")
   else
      writeln(io, "obj.. objvar =E= 0.0;");
      return
   end

   if ! isnothing(model.nlp_data) && model.nlp_data.has_objective
      obj_expr = MOI.objective_expr(model.nlp_data.evaluator)
      translate_function(io, model, obj_expr, is_parenthesis=true)
   else
      translate_function(io, model, model.objective)
   end

   writeln(io, ";")
   return
end

function translate_equations(
   io::GAMSTranslateStream,
   model::Optimizer
)
   for (i, con) in enumerate(model.linear_le_constraints)
      translate_equations(io, model, i + offset_linear_le(model), con.func, con.set)
   end
   for (i, con) in enumerate(model.linear_ge_constraints)
      translate_equations(io, model, i + offset_linear_ge(model), con.func, con.set)
   end
   for (i, con) in enumerate(model.linear_eq_constraints)
      translate_equations(io, model, i + offset_linear_eq(model), con.func, con.set)
   end
   for (i, con) in enumerate(model.quadratic_le_constraints)
      translate_equations(io, model, i + offset_quadratic_le(model), con.func, con.set)
   end
   for (i, con) in enumerate(model.quadratic_ge_constraints)
      translate_equations(io, model, i + offset_quadratic_ge(model), con.func, con.set)
   end
   for (i, con) in enumerate(model.quadratic_eq_constraints)
      translate_equations(io, model, i + offset_quadratic_eq(model), con.func, con.set)
   end
   for (i, con) in enumerate(model.complementarity_constraints)
      translate_equations(io, model, i + offset_complementarity(model), con.func, con.set)
   end
   for i in 1:model.m_nonlin
      translate_equations(io, model, i + offset_nonlin(model), MOI.constraint_expr(model.nlp_data.evaluator, i))
   end
   writeln(io, "")
   for (i, con) in enumerate(model.sos1_constraints)
      translate_equations(io, model, i, con.func, con.set)
   end
   for (i, con) in enumerate(model.sos2_constraints)
      translate_equations(io, model, i, con.func, con.set)
   end
   writeln(io, "")
   return
end

function translate_equations(
   io::GAMSTranslateStream,
   model::Optimizer,
   idx::Int,
   func::Union{MOI.ScalarAffineFunction{Float64}, MOI.ScalarQuadraticFunction{Float64}},
   set::MOI.LessThan{Float64}
)
   write(io, "eq$idx.. ")
   translate_function(io, model, func)
   writeln(io, " =L= " * dbl2str(set.upper) * ";")
   return
end

function translate_equations(
   io::GAMSTranslateStream,
   model::Optimizer,
   idx::Int,
   func::Union{MOI.ScalarAffineFunction{Float64}, MOI.ScalarQuadraticFunction{Float64}},
   set::MOI.GreaterThan{Float64}
)
   write(io, "eq$idx.. ")
   translate_function(io, model, func)
   writeln(io, " =G= " * dbl2str(set.lower) * ";")
   return
end

function translate_equations(
   io::GAMSTranslateStream,
   model::Optimizer,
   idx::Int,
   func::Union{MOI.ScalarAffineFunction{Float64}, MOI.ScalarQuadraticFunction{Float64}},
   set::MOI.EqualTo{Float64}
)
   write(io, "eq$idx.. ")
   translate_function(io, model, func)
   writeln(io, " =E= " * dbl2str(set.value) * ";")
   return
end

function translate_equations(
   io::GAMSTranslateStream,
   model::Optimizer,
   idx::Int,
   func::Expr
)
   if length(func.args) == 0
      return
   end
   @assert(length(func.args) == 3)

   write(io, "eq$idx.. ")
   translate_function(io, model, func.args[2], is_parenthesis=true)
   if func.args[1] == :(==)
      write(io, " =E= ")
   elseif func.args[1] == :(<=)
      write(io, " =L= ")
   else
      write(io, " =G= ")
   end
   translate_function(io, model, func.args[3], is_parenthesis=true)
   writeln(io, ";")
   return
end

function translate_equations(
   io::GAMSTranslateStream,
   model::Optimizer,
   idx::Int,
   func::MOI.VectorOfVariables,
   set::MOI.SOS1{Float64}
)
   write(io, "s1eq$idx(s1s$idx).. s1x$idx(s1s$idx) =e= ")
   for (i, vi) in enumerate(func.variables)
      if i > 1
         write(io, " + ")
      end
      translate_variable(io, model, vi.value)
      write(io, "\$sameas('$(vi.value)',s1s$idx)")
   end
   writeln(io, ";")
   return
end

function translate_equations(
   io::GAMSTranslateStream,
   model::Optimizer,
   idx::Int,
   func::MOI.VectorOfVariables,
   set::MOI.SOS2{Float64}
)
   write(io, "s2eq$idx(s2s$idx).. s2x$idx(s2s$idx) =e= ")
   for (i, vi) in enumerate(func.variables)
      if i > 1
         write(io, " + ")
      end
      translate_variable(io, model, vi.value)
      write(io, "\$sameas('$(vi.value)',s2s$idx)")
   end
   writeln(io, ";")
   return
end

function translate_equations(
   io::GAMSTranslateStream,
   model::Optimizer,
   idx::Int,
   func::MOI.VectorAffineFunction,
   set::MOI.Complements
)
   for i in 1:set.dimension
      row = filter(term -> term.output_index == i, func.terms)
      write(io, "eq$(idx)_$(i).. ")
      translate_function(io, model, row)
      if func.constants[i] < 0
         write(io, " - " * dbl2str(-func.constants[i]))
      elseif func.constants[i] > 0.0
         write(io, " + " * dbl2str(func.constants[i]))
      end
      writeln(io, " =N= 0;")
   end
end

function translate_vardata(
   io::GAMSTranslateStream,
   model::Optimizer
)
   for (i, var) in enumerate(model.variable_info)
      if is_fixed(var)
         translate_variable(io, model, i)
         writeln(io, ".fx = " * dbl2str(var.lower_bound) * ";")
         continue
      end
      if has_lower_bound(var)
         translate_variable(io, model, i)
         writeln(io, ".lo = " * dbl2str(var.lower_bound) * "; ")
      end
      if has_start(var)
         translate_variable(io, model, i)
         writeln(io, ".l = " * dbl2str(var.start) * "; ")
      end
      if has_upper_bound(var)
         translate_variable(io, model, i)
         writeln(io, ".up = " * dbl2str(var.upper_bound) * ";")
      end
   end
   write(io, "\n")
end

function translate_solve(
   io::GAMSTranslateStream,
   model::Optimizer,
   name::String
)
   # model statement
   if model.model_type == GAMS.MODEL_TYPE_MPEC || model.model_type == GAMS.MODEL_TYPE_MCP
      write(io, "Model $name / ")

      m = model.m + length(model.sos1_constraints) + length(model.sos2_constraints) + length(model.complementarity_constraints)

      first = true
      if model.objvar
         write(io, "obj")
         first = false
      end

      for i in 1:model.m
         if ! first
            write(io, ", ")
         end
         write(io, "eq$i")
         first = false
      end

      # add sos1 constraints
      for i in 1:length(model.sos1_constraints)
         if ! first
            write(io, ", ")
         end
         write(io, "s1eq$(i)(s1s$(i))")
         first = false
      end

      # add sos2 constraints
      for i in 1:length(model.sos2_constraints)
         if ! first
            write(io, ", ")
         end
         write(io, "s2eq$(i)(s2s$(i))")
         first = false
      end

      # add complementarity constraints
      for (i, comp) in enumerate(model.complementarity_constraints)
         d = comp.set.dimension
         for j in 1:d
            if ! first
               write(io, ", ")
            end
            var = filter(term -> term.output_index == j + d, comp.func.terms)
            var_str = translate_variable(model, var[1].scalar_term.variable_index.value)
            write(io, "eq$(i)_$(j).$(var_str)")
            first = false
         end
      end
      writeln(io, " /;")
   else
      writeln(io, "Model $name / all /;")
   end

   # solve statement
   write(io, "Solve $name using ")
   write(io, label(model.model_type))
   if model.model_type != GAMS.MODEL_TYPE_MCP && model.model_type != GAMS.MODEL_TYPE_CNS
      if model.sense == MOI.MAX_SENSE
         write(io, " maximizing ")
      else
         write(io, " minimizing ")
      end
      if model.objvar
         write(io, "objvar")
      elseif typeof(model.objective) == MOI.SingleVariable
         translate_variable(io, model, model.objective.variable.value)
      else
         error("GAMS needs obj variable")
      end
   end
   writeln(io, ";\n")
end


text:
using Flux
using Flux: @epochs
using GeometricFlux
using Serialization
using CUDA
using BSON: @save

DATAFILE="databig.blob"
dat = open(DATAFILE) do io
    deserialize(io)
end

# @show dat
#
model = Flux.Chain(
            CGConv((4, 1)),
            CGConv((4, 1)),
            CGConv((4, 1)),
            CGConv((4, 1)),
            CGConv((4, 1)),
            node_feature,
            X -> sum(X, dims=2),
            Dense(4,10), Dense(10,10),
            Dense(10,4), Dense(4,1)
           )

loss(gf, y) = Flux.Losses.mse(model(gf), y)

NUM_TRAIN = convert(Int, floor(length(dat) * 4/5))

dat_train = [(row["molecule"], parse(Float64, row["R2"])) for row in dat[1:NUM_TRAIN]]
dat_val  = [(row["molecule"], parse(Float64, row["R2"])) for row in dat[NUM_TRAIN+1:end]]

@epochs 10 Flux.train!(loss, Flux.params(model), dat_train, ADAM(0.05))

losses_val = [loss(x,y) for (x,y) in dat_val]
avg_loss_val = 1/length(losses_val) * sum(losses_val)
@show maximum(losses_val)
@show avg_loss_val

@save "model1.bson" model


text: using Double64s, BenchmarkTools


text:
using RDatasets, DataFrames, Gadfly

set_default_plot_size(10inch, 10inch)

barley = dataset("lattice", "barley")
levels!(barley.Year, ["1931", "1932"])

idx = [startswith(String(x),"No.") for x in barley.Variety]
plot(barley[idx,:],
     xgroup="Variety", ygroup="Site", x="Year", y="Yield",
     Geom.subplot_grid(Geom.line, Geom.point))


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

@otestset "GlassCat" begin
    # ensure there aren't any ambiguities
    @test isempty(detect_ambiguities(GlassCat))
    @test isempty(detect_unbound_args(GlassCat))

    include(joinpath("..", "..", "src", "GlassCat", "GlassTypes.jl"))
    include(joinpath(@__DIR__, "..", "..", "deps", "generate.jl"))

    CATALOG_NAME = "TEST_CAT"
    SOURCE_DIR = joinpath(@__DIR__, "..")
    SOURCE_FILE = joinpath(SOURCE_DIR, "$(CATALOG_NAME).agf")
    tmpdir = mktempdir()
    MAIN_FILE = joinpath(tmpdir, "AGF_TEST_CAT.jl")
    JL_DIR = tmpdir

    cat = Dict()

    @testset "Parsing Tests" begin
        cat = sourcefile_to_catalog(SOURCE_FILE)

        # test that all glasses were parsed correctly
        @test "MG463_" ∈ keys(cat)
        @test "MG523" ∈ keys(cat)
        @test "_5MG448" ∈ keys(cat)
        @test "MG452_STAR" ∈ keys(cat)

        # check that parsing all the variables worked for one example
        glass = cat["MG463_"]

        @test glass["raw_name"] == "MG463&"

        @test glass["dispform"] == 2
        @test glass["Nd"] == 1.0
        @test glass["Vd"] == 0.0
        @test glass["exclude_sub"] == 5
        @test glass["status"] == 6
        @test glass["meltfreq"] == -1

        @test glass["TCE"] == 19.8
        @test glass["p"] == 4.63
        @test glass["ΔPgF"] == 0.0
        @test glass["ignore_thermal_exp"] == 3

        @test glass["C1"] == 6.781409960
        @test glass["C2"] == 0.09868538020
        @test glass["C3"] == 0.003469529960
        @test glass["C4"] == 2.277823700
        @test glass["C5"] == 2.126055280
        @test glass["C6"] == 3227.490100
        @test glass["C7"] == 0.000000000
        @test glass["C8"] == 0.000000000
        @test glass["C9"] == 0.000000000
        @test glass["C10"] == 0.000000000

        @test glass["D₀"] == 2.226000000E-05
        @test glass["D₁"] == 5.150100000E-08
        @test glass["D₂"] == 1.309700000E-10
        @test glass["E₀"] == 4.382400000E-05
        @test glass["E₁"] == 4.966000000E-07
        @test glass["λₜₖ"] == 3.728000000E-01
        @test glass["temp"] == 2.002000000E+01

        @test glass["relcost"] == 1.1
        @test glass["CR"] == 1.2
        @test glass["FR"] == 1.3
        @test glass["SR"] == 1.4
        @test glass["AR"] == 1.5
        @test glass["PR"] == 1.6

        @test glass["λmin"] == 2.0
        @test glass["λmax"] == 14.0

        @test glass["transmission"] == [[2.0, 1.0, 3.0], [14.0, 4.0, 5.0]]
    end

    @testset "Module Gen Tests" begin
        generate_jls([CATALOG_NAME], MAIN_FILE, JL_DIR, SOURCE_DIR)
        include(MAIN_FILE)

        @test :MG463_ ∈ names(TEST_CAT, all = true)
        @test :MG523 ∈ names(TEST_CAT, all = true)
        @test :_5MG448 ∈ names(TEST_CAT, all = true)
        @test :MG452_STAR ∈ names(TEST_CAT, all = true)

        # test that the struct matches the dict
        glassd = cat["MG452_STAR"]
        glass = TEST_CAT.MG452_STAR

        @test glass.dispform == glassd["dispform"]

        @test glass.Nd == glassd["Nd"]
        @test glass.Vd == glassd["Vd"]
        @test glass.exclude_sub == glassd["exclude_sub"]
        @test glass.status == glassd["status"]
        @test glass.meltfreq == glassd["meltfreq"]

        @test glass.TCE == glassd["TCE"]
        @test glass.p == glassd["p"]
        @test glass.ΔPgF == glassd["ΔPgF"]
        @test glass.ignore_thermal_exp == glassd["ignore_thermal_exp"]

        @test glass.C1 == glassd["C1"]
        @test glass.C2 == glassd["C2"]
        @test glass.C3 == glassd["C3"]
        @test glass.C4 == glassd["C4"]
        @test glass.C5 == glassd["C5"]
        @test glass.C6 == glassd["C6"]
        @test glass.C7 == glassd["C7"]
        @test glass.C8 == glassd["C8"]
        @test glass.C9 == glassd["C9"]
        @test glass.C10 == glassd["C10"]

        @test glass.D₀ == glassd["D₀"]
        @test glass.D₁ == glassd["D₁"]
        @test glass.D₂ == glassd["D₂"]
        @test glass.E₀ == glassd["E₀"]
        @test glass.E₁ == glassd["E₁"]
        @test glass.λₜₖ == glassd["λₜₖ"]
        @test glass.temp == glassd["temp"]

        @test glass.relcost == glassd["relcost"]
        @test glass.CR == glassd["CR"]
        @test glass.FR == glassd["FR"]
        @test glass.SR == glassd["SR"]
        @test glass.AR == glassd["AR"]
        @test glass.PR == glassd["PR"]

        @test glass.λmin == glassd["λmin"]
        @test glass.λmax == glassd["λmax"]
    end

    @testset "Glass Tests" begin
        include(MAIN_FILE)

        g = TEST_CAT.MG523
        @test isapprox(OpticSim.GlassCat.index(g, ((g.λmin + g.λmax) / 2) * u"μm"), 3.1560980389455593, atol = 1e-14)
        @test_throws ErrorException OpticSim.GlassCat.index(g, (g.λmin - 1) * u"μm")
        @test_throws ErrorException OpticSim.GlassCat.index(g, (g.λmax + 1) * u"μm")
        @test isapprox(OpticSim.GlassCat.absairindex(500 * u"nm"), 1.0002741948670688, atol = 1e-14)
        @test isapprox(OpticSim.GlassCat.absairindex(600 * u"nm", temperature = 35 * u"°C", pressure = 2.0), 1.0005179096900811, atol = 1e-14)
        @test OpticSim.GlassCat.index(OpticSim.GlassCat.Air, 500 * u"nm") == 1.0
        @test OpticSim.GlassCat.index(OpticSim.GlassCat.Air, 600 * u"nm") == 1.0
        # @test isapprox(index(OpticSim.GlassCat.Vacuum, 500 * u"nm"), 0.9997258841958212, atol = 1e-14)

        # test against true values
        g = OpticSim.GlassCat.SCHOTT.N_BK7
        @test isapprox(OpticSim.GlassCat.index(g, 533 * u"nm"), 1.519417351519283, atol = 1e-14)
        @test isapprox(OpticSim.GlassCat.index(g, 533 * u"nm", temperature = 35 * u"°C"), 1.519462486258311, atol = 1e-14)
        @test isapprox(OpticSim.GlassCat.index(g, 533 * u"nm", pressure = 2.0), 1.518994119690216, atol = 1e-14)
        @test isapprox(OpticSim.GlassCat.index(g, 533 * u"nm", temperature = 35 * u"°C", pressure = 2.0), 1.519059871499476, atol = 1e-14)

        # test transmission
        @test OpticSim.GlassCat.absorption(OpticSim.GlassCat.Air, 500 * u"nm") == 0.0
        @test OpticSim.GlassCat.absorption(OpticSim.GlassCat.Air, 600 * u"nm") == 0.0
        # TODO these are currently taken from this package (i.e. regression tests), ideally we would get true values somehow
        @test OpticSim.GlassCat.absorption(g, 500 * u"nm") == 0.0002407228930225164
        @test OpticSim.GlassCat.absorption(g, 600 * u"nm") == 0.00022060722752440445
        @test OpticSim.GlassCat.absorption(g, 3000 * u"nm") == 0.04086604990127926
        @test OpticSim.GlassCat.absorption(g, 600 * u"nm", temperature = 35 * u"°C", pressure = 2.0) == 0.00022075540719494738

        # test that everything is alloc-less
        @test (@allocated OpticSim.GlassCat.absorption(g, 600 * u"nm")) == 0
        @test (@allocated OpticSim.GlassCat.index(g, 533 * u"nm")) == 0
        @test (@allocated OpticSim.GlassCat.SCHOTT.N_BK7) == 0

        @test OpticSim.GlassCat.glassforid(OpticSim.GlassCat.glassid(OpticSim.GlassCat.SCHOTT.N_BK7)) == OpticSim.GlassCat.SCHOTT.N_BK7

        @test OpticSim.GlassCat.isair(OpticSim.GlassCat.Air) == true
        @test OpticSim.GlassCat.isair(OpticSim.GlassCat.SCHOTT.N_BK7) == false

        # test MIL and model glasses
        fit_acc = 0.001
        bk7 = OpticSim.GlassCat.glassfromMIL(517642)
        @test OpticSim.GlassCat.glassforid(OpticSim.GlassCat.glassid(bk7)) == bk7
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(bk7, 0.5875618), 1.517, atol = fit_acc)
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(bk7, 0.533), 1.519417351519283, atol = fit_acc)
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(bk7, 0.743), 1.511997032563557, atol = fit_acc)
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(bk7, 533 * u"nm", temperature = 35 * u"°C", pressure = 2.0), 1.519059871499476, atol = fit_acc)

        t2 = OpticSim.GlassCat.glassfromMIL(1.135635)
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(t2, 0.5875618), 2.135, atol = fit_acc)

        fit_acc = 0.0001
        bk7 = OpticSim.GlassCat.modelglass(1.5168, 64.167336, -0.0009)
        @test OpticSim.GlassCat.glassforid(OpticSim.GlassCat.glassid(bk7)) == bk7
        @test isapprox(OpticSim.GlassCat.index(bk7, 0.5875618), 1.5168, atol = fit_acc)
        @test isapprox(OpticSim.GlassCat.index(bk7, 0.533), 1.519417351519283, atol = fit_acc)
        @test isapprox(OpticSim.GlassCat.index(bk7, 0.743), 1.511997032563557, atol = fit_acc)
        @test isapprox(OpticSim.GlassCat.index(bk7, 533 * u"nm", temperature = 35 * u"°C", pressure = 2.0), 1.519059871499476, atol = fit_acc)

        # test other glass
        @test OpticSim.GlassCat.index(OpticSim.GlassCat.CARGILLE.OG0608, 0.578) == 1.4596475735607324

        # make sure that the other functions work
        OpticSim.GlassCat.plot_indices(OpticSim.GlassCat.SCHOTT.N_BK7; polyfit = true, fiterror = true)
    end
end # testset GlassCat