text:
#
# Generic version of geotransform, fall through to including the handlers as the final two inputs
#
"""
  geotransform - function to transfom points between coordinate systems

  Usage:

    geotransform(::Type{T}, X)

        Transform geodesy point X to the coordinate system defined by T

        Examples:
            geotransform(ECEF,                   LLA{WGS84}(0,0,0))
            geotransform(ECEF{WGS84},            LLA{WGS84}(0,0,0))
            geotransform(CRS{SRID{:EPSG, 4978}}, LLA{WGS84}(0,0,0))
            geotransform(ENU{LLA{WGS84}(0,0,0)}, LLA{WGS84}(1,1,10))  # ENU frame centered on LLA{WGS84}(0,0,0)

    geotransform(::Type{T}, X, ref)

        Transform geodesy point X to the coordinate system defined by T using the reference data in ref

        Examples:
            geotransform(ECEF{UnknownDatum}, LLA{UnknownDatum}(0,0,0), Geodesy.eWGS84)  # supply the ellipse to use since there's no datum info
            geotransform(ENU, LLA{WGS84}(1,1,10), LLA{WGS84}(0,0,0))                    # ENU frame centered on LLA{WGS84}(0,0,0)

    geotransform(::Type{T}, X, handler_output, handler_input)

        Transform geodesy point X to the coordinate system defined by T using the package indicated by the package handlers

        Example:
            geotransform(ECEF{WGS84}, LLA{WGS84}(0,0,0), Geodesy.Proj4Handler, Geodesy.Proj4Handler)  # Get the Proj4 package to do the transform

    geotransform(::Type{T}, X, ref, handler_output, handler_input)

        Transform geodesy point X to the coordinate system defined by T using the reference data ref and package indicated by the handlers

        Example:

            # Get the Proj4 package to do the transform by suppling the SRIDs seperately
            geotransform(ECEF{UnknownDatum}, LLA{UnknownDatum}(0,0,0), (SRID{:epsg, 4978}, SRID{:epsg, 4326}), Geodesy.Proj4Handler, Geodesy.Proj4Handler)

"""
geotransform(X)         = geotransform(X, get_handler(X))
geotransform(X, Y)      = geotransform(X, Y, get_handler(X), get_handler(Y))
geotransform(X, Y, ref) = geotransform(X, Y, ref, get_handler(X), get_handler(Y))  # with a reference point

# error message for when no handler is defined
function geotransform{Handler1}(X, ::Type{Handler1})
    tOut = isa(X, DataType) ? X : typeof(X)
    tIn = !isa(Y, DataType) ? typeof(Y) : error("geotransform: expected a point type as the second input, got DataType: $(Y)")
    error("No methods defined to transform from $(tIn) to $(tOut)")
end
function geotransform{Handler1, Handler2}(X, Y, ::Type{Handler1}, ::Type{Handler2})
    tOut = isa(X, DataType) ? X : typeof(X)
    tIn = !isa(Y, DataType) ? typeof(Y) : error("geotransform: expected a point type as the second input, got DataType: $(Y)")
    error("No methods defined to transform from $(tIn) to $(tOut)")
end
function geotransform{Handler1, Handler2}(X, Y, ref, ::Type{Handler1}, ::Type{Handler2})
    tOut = isa(X, DataType) ? X : typeof(X)
    tIn = !isa(Y, DataType) ? typeof(Y) : error("geotransform: expected a point type as the second input, got DataType: $(Y)")
    error("No methods defined to transform from $(tIn) to $(tOut) using a reference point")
end

#
# version that return the transform parameters
#
geotransform_params(X)         = geotransform_params(X, get_handler(X))
geotransform_params(X, Y)      = geotransform_params(X, Y, get_handler(X), get_handler(Y))
geotransform_params(X, Y, ref) = geotransform_params(X, Y, ref, get_handler(X), get_handler(Y))  # with a reference point

#
# vectorized version (Proj4 point need special handling)
#
geotransform_vector{T}(X::Vector{T})                              = geotransform_vector(X, get_handler(eltype(T)))
geotransform_vector{T, U <: AbstractVector}(::Type{T}, Y::U)      = geotransform_vector(T, Y, get_handler(T), get_handler(eltype(U)))
geotransform_vector{T, U <: AbstractVector}(::Type{T}, Y::U, ref) = geotransform_vector(T, Y, ref, get_handler(T), get_handler(eltype(U)))  # with a reference point

#
# rules for infering one datum from another
#

# the methods
select_datum(::Type{UnknownDatum}, ::Type{UnknownDatum}) = UnknownDatum
select_datum{T <: KnownDatum}(::Type{T}, ::Type{T}) = T
select_datum{T1 <: KnownDatum, T2 <: KnownDatum}(::Type{T1}, ::Type{T2}) = error("Inputs type and output type have different datums")
select_datum{T <: KnownDatum}(::Type{UnknownDatum}, ::Type{T}) = T
select_datum{T <: KnownDatum}(::Type{T}, ::Type{UnknownDatum}) = T

# from 3 options
select_datum{T,U,V}(::Type{T}, ::Type{U}, ::Type{V}) = select_datum(select_datum(get_datum(T), get_datum(U)), get_datum(V))

# fall through methods
select_datum{T,U}(::Type{T}, ::Type{U}) = select_datum(get_datum(T), get_datum(U))
select_datum{T,U}(::T, ::Type{U})       = select_datum(get_datum(T), get_datum(U))
select_datum{T,U}(::Type{T}, ::U)       = select_datum(get_datum(T), get_datum(U))
select_datum{T,U}(::T, ::U)             = select_datum(get_datum(T), get_datum(U))

##########################################
# Special case, allow LL <-> LLA
##########################################

geotransform{T <: LL, U <: LLA}(::Type{T}, lla::U, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = LL{select_datum(add_param(T,U), U)}(lla.lon, lla.lat)
geotransform{T <: LLA, U <: LL}(::Type{T}, ll::U, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})  = LLA{select_datum(add_param(T,U), U)}(lla.lon, lla.lat, 0.0)

###############################
### LLA to ECEF coordinates ###
###############################

# reference ellipse inferred from the in / output types
function geotransform{T <: ECEF, U <: Union{LL, LLA}}(::Type{T}, X::U, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})
    oT = add_param(T, U)  # insert the datum if needed
    datum = select_datum(get_datum(oT), get_datum(U))
    ellipse = ellipsoid(datum)
    lla_to_ecef(oT, X, ellipse)
end

# reference ellipse specified (N.B. output type is assumed to be what the user wants)
function geotransform{T <: ECEF, U <: Union{LL, LLA}}(::Type{T}, X::U, ellipse::Ellipsoid, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})
    oT = add_param(T, U)  # insert the datum if needed
    lla_to_ecef(oT, X, ellipse)
end

# worker function
function lla_to_ecef{T}(::Type{T}, ll, ellipse::Ellipsoid)
     ϕdeg, λdeg, h = get_lat(ll), get_lon(ll), get_alt(ll)

    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)
    sinλ, cosλ = sind(λdeg), cosd(λdeg)

    N = ellipse.a / sqrt(1 - ellipse.e² * sinϕ^2)  # Radius of curvature (meters)

    x = (N + h) * cosϕ * cosλ
    y = (N + h) * cosϕ * sinλ
    z = (N * (1 - ellipse.e²) + h) * sinϕ

    return T(x, y, z)
end

###############################
### ECEF to LLA coordinates ###
###############################

# reference ellipse inferred from the in / output types (always using LL to parameterise atm)
function geotransform{T <: Union{LL, LLA}, U <: ECEF}(::Type{T}, X::U, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})
    oT = add_param(T, U)  # insert the datum if needed
    datum = select_datum(get_datum(oT), get_datum(U))
    ellipse = ellipsoid(datum)
    ecef_to_lla(oT, X, ellipse)
end

# reference ellipse specified (N.B. output type is assumed to be what the user wants)
function geotransform{T <: Union{LL, LLA}, U <: ECEF}(::Type{T}, X::U, ellipse::Ellipsoid, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})
    oT = add_param(T, U)  # insert the datum if needed
    lla_to_ecef(oT, X, ellipse)
end

# worker function
function ecef_to_lla{T <: Union{LL, LLA}}(::Type{T}, ecef, ellipse::Ellipsoid)
    x, y, z = getX(ecef), getY(ecef), getZ(ecef)

    p = hypot(x, y)
    θ = atan2(z*ellipse.a, p*ellipse.b)
    λ = atan2(y, x)
     ϕ = atan2(z + ellipse.e′² * ellipse.b * sin(θ)^3, p - ellipse.e²*ellipse.a*cos(θ)^3)

    N = ellipse.a / sqrt(1 - ellipse.e² * sin(ϕ)^2)  # Radius of curvature (meters)
    h = p / cos(ϕ) - N

    lla = T(rad2deg(λ), rad2deg(ϕ), h)

    return lla
end

# TESTING - compare accuracy of this method to ecef_to_lla using Geographic lib as a ground truth
function ecef_to_lla_test{T <: Union{LL, LLA}}(::Type{T}, ecef, ellipse::Ellipsoid)

    ellipse = ellipsoid(T)

    # termination tolerances for convergence
    hTol = 1e-6
    latTol = hTol / ellipse.a

    # go back into Geodetic via iterative method
    X, Y, Z = getX(ecef), getY(ecef), getZ(ecef)
    lon = atan2(Y, X)
    R = sqrt(X*X + Y*Y)

    # initialisation for iterative lat and h solution
    h = 0.0
    lat = atan2(Z,(1.0 - ellipse.e²) * R)

    converged = false
    maxIter = 10
    i=1
    while (!converged) && (i <= maxIter)
        Nlat = ellipse.a/(sqrt(1.0 - ellipse.e² * (sin(lat))^2))
        hNew = R/cos(lat) - Nlat
        latNew = atan2(Z, R * (1.0 -  ellipse.e² * Nlat / (Nlat+hNew)))
        converged = abs(latNew - lat) < latTol && abs(h-hNew) < hTol
        lat = latNew
        h = hNew
        i=i+1
    end

    return T(rad2deg(lon), rad2deg(lat), h)

end

###############################
### ECEF to ENU coordinates ###
###############################

# calling with not enough / too much info on the reference point
geotransform(::Type{ENU}, ecef::ECEF, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = error("You must supply the reference point in the ENU type or as an additional input")
geotransform{T <: Union{LL, LLA}}(::Type{ENU{T}}, ecef::ECEF, ll_ref, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = error("Don't specify both the reference point in the type and provide it as an argument")

# get the reference location from the input type
geotransform{T <: ENU, U <: ECEF}(::Type{T}, X::U, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = ecef_to_enu(T, X, get_refloc(T))

# the reference point is user specified
geotransform(::Type{ENU_NULL}, ecef::ECEF, ll_ref::Union{LLA, LL}, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = ecef_to_enu(ENU_NULL, ecef, ll_ref)
geotransform(::Type{ENU}, ecef::ECEF, ll_ref::Union{LLA, LL}, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = ecef_to_enu(ENU_NULL, ecef, ll_ref)

# if given an arbitrary reference, convert it to LL
geotransform(::Type{ENU}, ecef::ECEF, ref, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = geotransform(T, ecef, LL(ref), GeodesyHandler, GeodesyHandler)
geotransform(::Type{ENU_NULL}, ecef::ECEF, ref, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = geotransform(T, ecef, LL(ref), GeodesyHandler, GeodesyHandler)

# worker function
function ecef_to_enu{T, U <: ECEF, V <: Union{LLA, LL}}(::Type{T}, ecef::U, ll_ref::V)
    ϕdeg, λdeg, h =  get_lat(ll_ref), get_lon(ll_ref), get_alt(ll_ref)

    # get the reference point in ecef
    ecef_ref = V <: LL ? geotransform(U, ll_ref) : geotransform(U, LL(ll_ref))  # omit height here when calculating the reference point in ECEF.  Add it back at the end
    Δx = ecef.x - ecef_ref.x
    Δy = ecef.y - ecef_ref.y
    Δz = ecef.z - ecef_ref.z

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # R = [     -sinλ       cosλ  0.0
    #      -cosλ*sinϕ -sinλ*sinϕ cosϕ
    #       cosλ*cosϕ  sinλ*cosϕ sinϕ]
    #
    # east, north, up = R * [Δx, Δy, Δz]
    east  =      -sinλ * Δx +       cosλ * Δy +  0.0 * Δz
    north = -cosλ*sinϕ * Δx + -sinλ*sinϕ * Δy + cosϕ * Δz
    up    =  cosλ*cosϕ * Δx +  sinλ*cosϕ * Δy + sinϕ * Δz

    return T(east, north, up - h)  # I think the height vector is the same direction as the ellipsoidal height at the ENU's origin
end

###############################################
### Transformation matrices for ECEF -> ENU ###
###############################################

# reference point in the template
geotransform_params{R <: Union{LLA, LL}, T <: ECEF}(::Type{ENU{R}}, ::Type{T}, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = geotransform_params(T, ENU, R, GeodesyHandler, GeodesyHandler)

# reference point as a parameter
function geotransform_params{T <: ECEF, Ref <: Union{LLA, LL}}(::Type{ENU}, ::Type{T}, ll_ref::Ref, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})

    ϕdeg, λdeg, h =  get_lat(ll_ref), get_lon(ll_ref), get_alt(ll_ref)

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # Reference
    ecef_ref = Ref <: LL ? geotransform(T, ll_ref) : geotransform(T, LL(ll_ref))  # omit height here when calculating the reference point in ECEF.  Add it back at the end

    R = @fsa([-sinλ        cosλ        0.0;
              -cosλ*sinϕ   -sinλ*sinϕ    cosϕ;
               cosλ*cosϕ    sinλ*cosϕ    sinϕ])

    t = -R * Vec(ecef_ref.x, ecef_ref.y, ecef_ref.z) - @fsa([0.0, 0.0, h])
    return (R, t)

end

###############################
### ENU to ECEF coordinates ###
###############################

# convert to ECEF when the reference position is included in the ENU template
function geotransform{T <: ECEF, U}(::Type{T}, enu::ENU{U}, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})   # I cant template U into anything because its a value type :-(
    ref = (isa(U, LL) || isa(U, LLA)) ? U : LL(U)
    enu_to_ecef(add_param(T, U), enu, ref)
end

# convert to ECEF when the LL reference position is not included in the ENU template
geotransform{T <: ECEF}(::Type{T}, enu::ENU_NULL, ll_ref::Union{LLA, LL}, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = enu_to_ecef(add_param(T, typeof(ll_ref)), enu, ll_ref)

# if given an arbitrary reference, convert it to LL
geotransform{T <: ECEF}(::Type{T}, enu::ENU_NULL, ref, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = geotransform(T, enu, LL(ref), GeodesyHandler, GeodesyHandler)

# spit an error if provided with a reference point in the template and as a parameter
geotransform{T <: ECEF}(::Type{T}, enu::ENU_NULL, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = error("You must supply the reference point in the ENU type or as an additional input")
geotransform{T <: ECEF}(::Type{T}, enu::ENU, ll_ref::Union{LLA, LL}, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = error("Dont supply both the reference point in thee template and as a parameter")

# worker function
function enu_to_ecef{T <: ECEF, U <: Union{LL, LLA}}(::Type{T}, enu::ENU, ll_ref::U)

    ϕdeg, λdeg, h =  get_lat(ll_ref), get_lon(ll_ref), get_alt(ll_ref)

    # Reference
    ecef_ref = U <: LL ? geotransform(T, ll_ref) : geotransform(T, LL(ll_ref))  # omit height here when calculating the reference point in ECEF.  Add it back later

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # Rᵀ = [-sinλ -cosλ*sinϕ cosλ*cosϕ
    #        cosλ -sinλ*sinϕ sinλ*cosϕ
    #         0.0       cosϕ      sinϕ]
    # Δx, Δy, Δz = Rᵀ * [east, north, up]
    Δx = -sinλ * enu.east + -cosλ*sinϕ * enu.north + cosλ*cosϕ * (enu.up + h)
    Δy =  cosλ * enu.east + -sinλ*sinϕ * enu.north + sinλ*cosϕ * (enu.up + h)
    Δz =   0.0 * enu.east +       cosϕ * enu.north +      sinϕ * (enu.up + h)

    X = ecef_ref.x + Δx
    Y = ecef_ref.y + Δy
    Z = ecef_ref.z + Δz

    return T(X,Y,Z)
end

###############################################
### Transformation matrices for ENU -> ECEF ###
###############################################

geotransform_params{T <: ECEF, R <: Union{LLA, LL}}(::Type{T}, ::Type{ENU{R}}, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) = geotransform_params(T, ENU, R, GeodesyHandler, GeodesyHandler)

# no LLA reference provided in the ENU type so a LLA ref must be supplied
# N.B. this will ignore the LLA refernce in ENU template
function geotransform_params{T <: ECEF, U <: ENU, Ref <: Union{LLA, LL}}(::Type{T}, ::Type{U}, ll_ref::Ref, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})

    ϕdeg, λdeg, h =  get_lat(ll_ref), get_lon(ll_ref), get_alt(ll_ref)

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # Reference
    ecef_ref = Ref <: LL ? geotransform(T, ll_ref) : geotransform(T, LL(ll_ref))  # omit height here when calculating the reference point in ECEF.  Add it back later

    R = @fsa([-sinλ  -cosλ*sinϕ    cosλ*cosϕ;
               cosλ  -sinλ*sinϕ    sinλ*cosϕ;
               0.0    cosϕ            sinϕ])

    t = Vec(ecef_ref.x, ecef_ref.y, ecef_ref.z) + R * @fsa([0.0, 0.0, h])

    return (R, t)

end

##############################
### LLA to ENU coordinates ###
##############################

# no convert the point to ECEF then go from there
geotransform{T <: ENU, U <: Union{LLA, LL}}(::Type{T}, lla::U, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) =
    geotransform(T, ECEF{select_datum(T, U)}(lla), GeodesyHandler, GeodesyHandler)

# convert the point to ECEF then go from there
geotransform{T <: ENU, U <: Union{LLA, LL}, V}(::Type{T}, lla::U, ref::V, ::Type{GeodesyHandler}, ::Type{GeodesyHandler}) =
    geotransform(T, ECEF{select_datum(T, U, V)}(lla), ref, GeodesyHandler, GeodesyHandler)

################################
### ENU to LLA coordinates   ###
################################

geotransform{T <: Union{LLA, LL}}(::Type{T}, enu::ENU_NULL, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})  = error("You must supply the reference point in the ENU type or as an additional input")

# no convert the point to ECEF then go from there
function geotransform{T <: Union{LLA, LL}, U}(::Type{T}, enu::ENU{U}, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})
    eT = ECEF{select_datum(T, U)}  # find a datum somewhere
    T(geotransform(eT, enu, GeodesyHandler, GeodesyHandler))
end

# convert the point to ECEF then go from there
function geotransform{T <: Union{LLA, LL}, U <: ENU, V}(::Type{T}, enu::U, ref::V, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})
    eT = ECEF{select_datum(T, U, V)} # find a datum somewhere
    T(geotransform(eT, enu, ref, GeodesyHandler, GeodesyHandler))
end

#############################################
### Decimal <=> Degrees, Minutes, Seconds ###
#############################################

function decimal2dms(x::Float64)
    d = trunc(x, 0)
    ms = 60 * abs(x - d)
    m = trunc(ms, 0)
    s = 60 * rem(ms, 1)

    return d, m, s
end

function dms2decimal(d::Float64, m::Float64, s::Float64)
    signbit(d) ? d - m/60 - s/3600 :
                 d + m/60 + s/3600
end

#################################################
### A vectorized version of the transforms
#################################################

# vectorizied transform for geodesy types
function geotransform_vector{T, U}(::Type{T}, X::Vector{U}, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})

    # make sure the output parameter is filled
    oT =  add_param(T, U)
    Xout = Vector{oT}(length(X))
    @inbounds for i = 1:length(X)
        Xout[i] = geotransform(oT, X[i])
    end
    return Xout
end

# vectorizied transform for geodesy types
function geotransform_vector{T, U, V}(::Type{T}, X::Vector{U}, ll_ref::V, ::Type{GeodesyHandler}, ::Type{GeodesyHandler})

    # make sure the output parameter is filled
    oT =  add_param(T, U)
    Xout = Vector{oT}(length(X))
    @inbounds for i = 1:length(X)
        Xout[i] = geotransform(oT, X[i], ll_ref)
    end
    return Xout
end


text:
using PortfolioOpt
using Test

@testset "PortfolioOpt.jl" begin
    # Write your tests here.
end

ext: jl


text:
"""
# Get all sources of economic data.

    get_sources()

## Examples

```jldoctest
julia> x = get_sources()
julia> println(x)
```
"""
function get_sources()

    url = "https://api.stlouisfed.org/fred/sources"
    parameters = Dict("api_key" => ENV["API_FRED"],
                      "file_type" => "json")

    response = HTTP.request("GET", url; query = parameters)
    results = JSON.parse(String(response.body))["sources"]

    v_keys = ["link" in keys(results[i]) for i ∈ eachindex(results)]
    link = Any["" for i ∈ eachindex(results)]
    for i ∈ eachindex(results)
        if v_keys[i] == true
            link[i] = results[i]["link"]
        else
            link[i] = missing
        end
    end

    x = DataFrame(
        id = [results[i]["id"] for i ∈ eachindex(results)],
        name = [results[i]["name"] for i ∈ eachindex(results)],
        link = link
    )

    return(x)
end


text:
## --------------------------------------------------------------------
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Human1")

import MAT

import Chemostat_Human1
const H1 = Chemostat_Human1

import Chemostat
const Ch = Chemostat
const ChU = Ch.Utils
const ChLP = Ch.LP

## --------------------------------------------------------------------
# Tools
fba_objval(model) = ChLP.fba(model, H1.BIOMASS_IDER).obj_val

## --------------------------------------------------------------------
# Processing all models
for cell_line in H1.CELL_NAMES
    for (const_name, const_level) in zip(H1.CONSTRAINTS_NAMES,
                                        H1.CONSTRAINTS_LEVELS)
        for model_type in [H1.MODEL_TYPE_KEY, H1.EC_MODEL_TYPE_KEY]

            ChU.println_inmw("\n")
            ChU.tagprintln_inmw("PROCESSING:",
                "\ncellLine:     ", cell_line,
                "\nconst_level:  ", const_level,
                "\nconst_name:   ", const_name,
                "\nmodel_type:   ", model_type,
                "\n"
            )

            model_file = H1.get_human1_model_file(cell_line, const_level, model_type)
            if isfile(model_file)
                msg = string(relpath(model_file), " found!!")
                ChU.println_inmw(msg)
            else
                msg = string("ERROR: ", relpath(model_file), " not found!!")
                ChU.println_inmw(msg)
                continue
            end

            ## --------------------------------------------------------------------
            # Check for existing file
            fva_pp_model_file = H1.get_fva_pp_model_file(cell_line, const_level, model_type)
            if isfile(fva_pp_model_file)
                ChU.println_inmw(relpath(fva_pp_model_file), " exist, skipping!")
                continue
            end

            ## --------------------------------------------------------------------
            # Loading
            model_key, model_mat = first(MAT.matread(model_file));
            model = ChU.MetNet(model_mat; reshape = true);
            ChU.clampfields!(model, [:b, :lb, :ub];
                abs_max = H1.MAX_ABS_BOUND, zeroth = H1.ZEROTH)
            ChU.println_inmw("size: ", size(model))
            println("Orig model objval: ", fba_objval(model))

            ## --------------------------------------------------------------------
            # Redefine all rxns to be a single reaction fwd defined
            function is_exchange(model::ChU.MetNet, ider::ChU.IDER_TYPE)
                idx = ChU.rxnindex(model, ider)
                subsys = model.subSystems[idx]
                return occursin(H1.EXCH_SUBSYS_HINT, string(subsys))
            end
            exchs = filter((rxn) -> is_exchange(model, rxn), model.rxns);
            found_buffer = Set()
            for exch in exchs
                exch in found_buffer && continue
                if endswith(exch, "_REV")
                    rxn_rev = exch
                    rxn = replace(exch, "REV_")
                else
                    rxn = exch
                    rxn_rev = string(exch, "_REV")
                end
                push!(found_buffer, rxn, rxn_rev)

                rxn_exist = rxn in exchs
                rxn_rev_exist = rxn_rev in exchs

                if rxn_exist && rxn_rev_exist
                    # merge bounds
                    ChU.bounds!(model, rxn,
                        # [A <- ]lb == -[ -> A]ub
                        min(ChU.lb(model, rxn), -ChU.ub(model, rxn_rev)),
                        max(ChU.ub(model, rxn), -ChU.lb(model, rxn_rev)),
                    )
                    # Close rev
                    ChU.bounds!(model, rxn_rev, 0.0, 0.0)
                elseif rxn_rev_exist
                    # revert rev reaction
                    ChU.invert_rxn!(model, rxn_rev; rename = rxn)
                end
            end
            ChU.println_inmw("After first preprocessing, objval: ", fba_objval(model))

            # Checking
            for exch in exchs
                if ChU.isbkwd_defined(model, exch) &&
                        ChU.isopen(model, exch)
                    error(string(exch, "is an open backward defined reaction"))
                end
            end

            ## --------------------------------------------------------------------
            # FVA preprocessing
            fva_pp_model = ChLP.fva_preprocess(model, check_obj = H1.BIOMASS_IDER);
            ChU.println_inmw("size: ", size(fva_pp_model))
            ChU.println_inmw("fva preprocessed model, objval: ", fba_objval(model))

            ## --------------------------------------------------------------------
            # Saving
            mkpath(fva_pp_model_file |> dirname)
            ChU.save_data(fva_pp_model_file, ChU.compressed_model(fva_pp_model))

        end

    end
end


text:
@inline Base.zero(a::SA) where {SA<:TupleVector} = zeros(SA)
@inline Base.zero(a::Type{SA}) where {SA<:TupleVector} = zeros(SA)

@inline Base.zeros(::Type{SA}) where {SA<:TupleVector{N}} where N = _zeros(Val(N), SA)
@generated function _zeros(::Val{N}, ::Type{SA}) where SA<:TupleVector{N} where N
    T = eltype(SA)
    if T == Any
        T = Float64
    end
    v = [:(zero($T)) for i = 1:N]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end

@inline Base.ones(::Type{SA}) where {SA<:TupleVector{N}} where N = _ones(Val(N), SA)
@generated function _ones(::Val{N}, ::Type{SA}) where SA<:TupleVector{N} where N
    T = eltype(SA)
    if T == Any
        T = Float64
    end
    v = [:(one($T)) for i = 1:N]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end

@inline Base.fill(val,::SA) where {SA<:TupleVector{N}} where N = _fill(val, Val(N), SA)
@inline Base.fill(val,::Type{SA}) where {SA<:TupleVector{N}} where N = _fill(val, Val(N), SA)
@generated function _fill(val,::Val{N},::Type{SA}) where SA<:TupleVector{N} where N
    v = [:val for i = 1:N]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end

# Also consider randcycle, randperm? Also faster rand!(staticarray, collection)

using Random
import Random: SamplerType, AbstractRNG
@inline Base.rand(rng::AbstractRNG,::Type{SA},dims::Dims) where SA<:TupleVector = rand!(rng, Array{SA}(undef, dims), SA)
@inline Base.rand(rng::AbstractRNG,::SamplerType{SA}) where SA<:TupleVector{N} where N = _rand(rng, Val(N), SA)

@generated function _rand(rng::AbstractRNG,::Val{N},::Type{SA}) where {SA <: TupleVector{N}} where N
    T = eltype(SA)
    if T == Any
        T = Float64
    end
    v = [:(rand(rng, $T)) for i = 1:N]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end

@inline Base.rand(rng::AbstractRNG,range::AbstractArray, ::Type{SA}) where {SA<:TupleVector{N}} where N = _rand(rng, range, Val(N), SA)
@inline Base.rand(range::AbstractArray,::Type{SA}) where {SA<:TupleVector{N}} where N = _rand(Random.GLOBAL_RNG, range, Val(N), SA)
@generated function _rand(rng::AbstractRNG, range::AbstractArray, ::Val{N}, ::Type{SA}) where SA<:TupleVector{N} where N
    v = [:(rand(rng, range)) for i = 1:N]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end

#@inline Base.rand(rng::MersenneTwister, range::AbstractArray, ::Type{SA}) where {SA <: TupleVector} = _rand(rng, range, Val(SA), SA)

@inline Base.randn(rng::AbstractRNG,::Type{SA}) where SA<:TupleVector{N} where N = _randn(rng, Val(N), SA)
@generated function _randn(rng::AbstractRNG,::Val{N},::Type{SA}) where SA<:TupleVector{N} where N
    T = eltype(SA)
    if T == Any
        T = Float64
    end
    v = [:(randn(rng, $T)) for i = 1:N]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end

@inline Random.randexp(rng::AbstractRNG,::Type{SA}) where SA<:TupleVector{N} where N = _randexp(rng, Val(N), SA)
@generated function _randexp(rng::AbstractRNG,::Val{N},::Type{SA}) where SA <: TupleVector{N} where N
    T = eltype(SA)
    if T == Any
        T = Float64
    end
    v = [:(randexp(rng, $T)) for i = 1:N]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end

# Mutable versions

# Why don't these two exist in Base?
# @generated function Base.zeros!{SA <: TupleVector}(a::SA)
# @generated function Base.ones!{SA <: TupleVector}(a::SA)

@inline Base.fill!(a::SA,val) where SA<:TupleVector{N} where N = _fill!(Val(N), a, val)
@generated function _fill!(::Val{N},a::SA,val) where SA<:TupleVector{N} where N
    exprs = [:(a[$i] = valT) for i = 1:N]
    return quote
        @_inline_meta
        valT = convert(eltype(SA), val)
        @inbounds $(Expr(:block, exprs...))
        return a
    end
end

@inline Random.rand!(rng::AbstractRNG,a::SA) where SA<:TupleVector{N} where N = _rand!(rng, Val(N), a)
@generated function _rand!(rng::AbstractRNG,::Val{N},a::SA) where SA <: TupleVector{N} where N
    exprs = [:(a[$i] = rand(rng, eltype(SA))) for i = 1:N]
    return quote
        @_inline_meta
        @inbounds $(Expr(:block, exprs...))
        return a
    end
end

@inline Random.rand!(rng::MersenneTwister,a::SA) where SA<:TupleVector{N,Float64} where N = _rand!(rng, Val(N), a)

@inline Random.randn!(rng::AbstractRNG,a::SA) where SA<:TupleVector{N} where N = _randn!(rng, Val(N), a)
@generated function _randn!(rng::AbstractRNG,::Val{N},a::SA) where SA<:TupleVector{N} where N
    exprs = [:(a[$i] = randn(rng, eltype(SA))) for i = 1:N]
    return quote
        @_inline_meta
        @inbounds $(Expr(:block, exprs...))
        return a
    end
end

@inline Random.randexp!(rng::AbstractRNG,a::SA) where SA<:TupleVector{N} where N = _randexp!(rng, Val(N), a)
@generated function _randexp!(rng::AbstractRNG,::Val{N},a::SA) where SA<:TupleVector{N} where N
    exprs = [:(a[$i] = randexp(rng, eltype(SA))) for i = 1:N]
    return quote
        @_inline_meta
        @inbounds $(Expr(:block, exprs...))
        return a
    end
end


text:
module TowerFields
    # export FQ2, FQ3, FQ4, FQ6, FQ12, AfflinePoint
    # include("Common.jl")

    # module Types
    #     export FQ2, FQ3, FQ4, FQ6, FQ12, AfflinePoint
    #     include("Common.jl")
    # end

    module Common
        const FQ2 = Tuple{T, T} where T
        const FQ3 = Tuple{T, T, T} where T
        const FQ6 = Tuple{T, T, T} where {T <: FQ2}
        const FQ12 = Tuple{T, T} where {T <: FQ6}
        zero(a::FQ2) = (zero(a[1]), zero(a[2]))
        zero(a::FQ3) = (zero(a[1]), zero(a[2]), zero(a[3]))
        zero(a) = Base.zero(a)
        one(a::FQ2) = (one(a[1]), zero(a[2]))
        one(a::FQ3) = (one(a[1]), zero(a[2]), zero(a[3]))
        one(a) = Base.one(a)
    end

    module BigIntField
        using ..Common: FQ2, FQ3, zero, one, iszero, isone
        include("BigIntField.jl")
    end

    module TowerField12
        using ..Common: FQ2, FQ6, FQ12
        using ..BigIntField
        const β = -1
        const ξ = (1, 1)
        const γ = ((0, 0), (1, 0), (0, 0))
        mulroot(x) = -x # β = -1
        mulroot(x::FQ2) = (x[1] - x[2], x[1] + x[2]) # ξ = (1, 1)
        mulroot(x::FQ6) = (mulroot(x[3]), x[1], x[2]) # γ = ((0, 0), (1, 0), (0, 0))
        function create_frobmap(q)
            g2 = BigIntField.powmod(q, mulroot, β, fld(q - 1, 2))
            g2sq = BigIntField.sqmod(q, mulroot, g2)
            g6 = BigIntField.powmod(q, mulroot, ξ, fld(q - 1, 3))
            g6sq = BigIntField.sqmod(q, mulroot, g6)
            g12 = BigIntField.powmod(q, mulroot, γ, fld(q - 1, 2))
            g12sq = BigIntField.sqmod(q, mulroot, g12)
            frobmap(x::FQ2) = (g2, g2sq)
            frobmap(x::FQ6) = (g6, g6sq)
            frobmap(x::FQ12) = (g12, g12sq)
            return frobmap
        end
    end

    module BLS12381Field
        using ..Common: zero, one, iszero, isone, FQ2
        using ..BigIntField
        using ..TowerField12

        # Field Params
        const ABS_X = 15132376222941642752
        const X_IS_NEG = true
        const R = 52435875175126190479447740508185965837690552500527637822603658699938581184513
        const Q = 4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559787

        # Frobenius map function
        const frobmap = TowerField12.create_frobmap(Q)

        # API
        add(a, b) = BigIntField.addmod(Q, a, b)
        double(a) = BigIntField.addmod(Q, a, a)
        sub(a, b) = BigIntField.submod(Q, a, b)
        sub(a) = BigIntField.submod(Q, zero(a), a)
        mul(a, b) = BigIntField.mulmod(Q, TowerField12.mulroot, a, b)
        square(a) = BigIntField.sqmod(Q, TowerField12.mulroot, a)
        inv(a) = BigIntField.invmod(Q, TowerField12.mulroot, a)
        pow(a, b) = BigIntField.powmod(Q, TowerField12.mulroot, a, b)
        frob(a) = BigIntField.frobmod(Q, frobmap, TowerField12.mulroot, a)
        conj(a::FQ2) = (a[1], sub(a[2]))
        powx(a::FQ2) =  X_IS_NEG ? conj(pow(a, ABS_X)) : pow(a, ABS_X)
        div(a, b) = mul(a, inv(b))
        const +, -, *, //, ^, neg = add, sub, mul, div, pow, sub
    end

end


text:
module MathOptInterface

"""
    ModelLike

Abstract supertype for objects that implement the "Model" interface for defining
an optimization problem.
"""
abstract type ModelLike end

# This allows to use `ModelLike`s in broadcast calls without the need to
# embed it in a `Ref`
Base.broadcastable(model::ModelLike) = Ref(model)

Base.show(io::IO, model::ModelLike) = Utilities.print_with_acronym(io, summary(model))

"""
    AbstractOptimizer

Abstract supertype for objects representing an instance of an optimization problem
tied to a particular solver. This is typically a solver's in-memory representation.
In addition to `ModelLike`, `AbstractOptimizer` objects let you solve the
model and query the solution.
"""
abstract type AbstractOptimizer <: ModelLike end

"""
    optimize!(optimizer::AbstractOptimizer)

Start the solution procedure.
"""
function optimize! end

"""
    write_to_file(model::ModelLike, filename::String)

Writes the current model data to the given file.
Supported file types depend on the model type.
"""
function write_to_file end

"""
    read_from_file(model::ModelLike, filename::String)

Read the file `filename` into the model `model`. If `model` is non-empty, this may
throw an error.

Supported file types depend on the model type.

### Note

Once the contents of the file are loaded into the model, users can query the variables via
`get(model, ListOfVariableIndices())`. However, some filetypes, such as LP files, do not
maintain an explicit ordering of the variables. Therefore, the returned list may be in an
arbitrary order. To avoid depending on the order of the indices, users should look up each
variable index by name: `get(model, VariableIndex, "name")`.
"""
function read_from_file end

"""
    is_empty(model::ModelLike)

Returns `false` if the `model` has any model attribute set or has any variables or constraints.
Note that an empty model can have optimizer attributes set.
"""
function is_empty end

"""
    empty!(model::ModelLike)

Empty the model, that is, remove all variables, constraints and model attributes but not optimizer attributes.
"""
function empty! end

"""
    copy_to(dest::ModelLike, src::ModelLike; copy_names=true, warn_attributes=true)

Copy the model from `src` into `dest`. The target `dest` is emptied, and all
previous indices to variables or constraints in `dest` are invalidated. Returns
a dictionary-like object that translates variable and constraint indices from
the `src` model to the corresponding indices in the `dest` model.

If `copy_names` is `false`, the `Name`, `VariableName` and `ConstraintName`
attributes are not copied even if they are set in `src`. If a constraint that
is copied from `src` is not supported by `dest` then an
[`UnsupportedConstraint`](@ref) error is thrown. Similarly, if a model, variable
or constraint attribute that is copied from `src` is not supported by `dest`
then an [`UnsupportedAttribute`](@ref) error is thrown. Unsupported *optimizer*
attributes are treated differently:

* If `warn_attributes` is `true`, a warning is displayed, otherwise,
* the attribute is silently ignored.

### Example

```julia
# Given empty `ModelLike` objects `src` and `dest`.

x = add_variable(src)

is_valid(src, x)   # true
is_valid(dest, x)  # false (`dest` has no variables)

index_map = copy_to(dest, src)
is_valid(dest, x) # false (unless index_map[x] == x)
is_valid(dest, index_map[x]) # true
```
"""
function copy_to end

include("error.jl")
include("indextypes.jl")
include("functions.jl")
include("sets.jl")
include("attributes.jl")
include("constraints.jl")
include("modifications.jl")
include("variables.jl")
include("nlp.jl")

# submodules
include("Utilities/Utilities.jl") # MOI.Utilities
include("Test/Test.jl")           # MOI.Test
include("Bridges/Bridges.jl")     # MOI.Bridges
include("Benchmarks/Benchmarks.jl")

end


text:
function set_bdx(Qbase, rho, u, v, p, ch, nch, icell, j)
    Qbase[icell, j, 1] = rho
    Qbase[icell, j, 2] = u
    Qbase[icell, j, 3] = v
    Qbase[icell, j, 4] = p
    for ns in 1:nch
        Qbase[icell, j, 4+ns] = ch[ns]
    end
    return Qbase
end

function set_bdy(Qbase, rho, u, v, p, ch, nch, jcell, i)
    Qbase[i, jcell, 1] = rho
    Qbase[i, jcell, 2] = u
    Qbase[i, jcell, 3] = v
    Qbase[i, jcell, 4] = p
    for ns in 1:nch
        Qbase[i, jcell, 4+ns] = ch[ns]
    end
    return Qbase
end

function set_boundary(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rhat, nval, nch)
    """bdcon[i][j]
        i : 境界番号(x-,x+ ,y-,y+)
        j=1-6 : "bd1_con":"2",
                "bd1_rho":"1.0",
                "bd1_u"  :"300.0",
                "bd1_v"  :"0.0",
                "bd1_p"  :"1.0",
                "bd1_N2" :"1.0",
                "bd1_N"  :"0.0",
                "bd1_T"  :"300.0",
    """
    rho = 0.0
    u   = 0.0
    v   = 0.0
    p   = 0.0
    ch  = zeros(nch)

    # bd1 = x-
    # bd2 = x+
    bdx         = [1, cellxmax]
    bdx_nei     = [2, cellxmax-1]
    bdx_nei_rev = [cellxmax-1, 2]
    for dim in 1:2
        if Int(bdcon[dim][1]) == 00
            for j in 1:cellymax
                rho = bdcon[dim][2]
                u   = bdcon[dim][3]
                v   = bdcon[dim][4]
                p   = bdcon[dim][5]
                for ns in 1:nch
                    ch[ns] =  bdcon[dim][5+ns]
                end
                set_bdx(Qbase, rho, u, v, p, ch, nch, bdx[dim], j)
            end
        elseif Int(bdcon[dim][1]) == 01
            for j in 1:cellymax
                rho = bdcon[dim][2]
                u   = bdcon[dim][3]
                v   = bdcon[dim][4]
                T   = bdcon[dim][6+nch]
                p   = rho*Rhat[bdx_nei[dim],j]*T
                for ns in 1:nch
                    ch[ns] =  bdcon[dim][5+ns]
                end
                set_bdx(Qbase, rho, u, v, p, ch, nch, bdx[dim], j)
            end
        elseif Int(bdcon[dim][1]) == 20
            for j in 1:cellymax
                rho = Qbase[bdx_nei[dim],j,1]
                u   = Qbase[bdx_nei[dim],j,2]
                v   = Qbase[bdx_nei[dim],j,3]
                p   = Qbase[bdx_nei[dim],j,4]
                for ns in 1:nch
                    ch[ns] =  Qbase[bdx_nei[dim],j,4+ns]
                end
                set_bdx(Qbase, rho, u, v, p, ch, nch, bdx[dim], j)
            end
        elseif Int(bdcon[dim][1]) == 30
            for j in 1:cellymax
                rho = Qbase[bdx_nei[dim],j,1]
                u   = -Qbase[bdx_nei[dim],j,2]
                v   = -Qbase[bdx_nei[dim],j,3]
                p   = Qbase[bdx_nei[dim],j,4]
                for ns in 1:nch
                    ch[ns] =  Qbase[bdx_nei[dim],j,4+ns]
                end
                set_bdx(Qbase, rho, u, v, p, ch, nch, bdx[dim], j)
            end
        elseif Int(bdcon[dim][1]) == 31
            for j in 1:cellymax
                u   = -Qbase[bdx_nei[dim],j,2]
                v   = -Qbase[bdx_nei[dim],j,3]
                p   = Qbase[bdx_nei[dim],j,4]
                T   = bdcon[dim][6+nch]
                rho = p/(Rhat[bdx_nei[dim],j]*T)
                for ns in 1:nch
                    ch[ns] =  Qbase[bdx_nei[dim],j,4+ns]
                end
                set_bdx(Qbase, rho, u, v, p, ch, nch, bdx[dim], j)
            end
        elseif Int(bdcon[dim][1]) == 40
            for j in 1:cellymax
                rho = Qbase[bdx_nei_rev[dim],j,1]
                u   = Qbase[bdx_nei_rev[dim],j,2]
                v   = Qbase[bdx_nei_rev[dim],j,3]
                p   = Qbase[bdx_nei_rev[dim],j,4]
                for ns in 1:nch
                    ch[ns] =  Qbase[bdx_nei_rev[dim],j,4+ns]
                end
                set_bdx(Qbase, rho, u, v, p, ch, nch, bdx[dim], j)
            end
        end
    end

    # bd3 = y-
    # bd4 = y+
    bdy     = [1, cellymax]
    bdy_nei     = [2, cellymax-1]
    bdy_nei_rev = [cellymax-1, 2]
    for dim in 3:4
        if Int(bdcon[dim][1]) == 00
            for i in 1:cellxmax
                rho = bdcon[dim][2]
                u   = bdcon[dim][3]
                v   = bdcon[dim][4]
                p   = bdcon[dim][5]
                for ns in 1:nch
                    ch[ns] =  bdcon[dim][5+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
        elseif Int(bdcon[dim][1]) == 01
            for i in 1:cellxmax
                rho = bdcon[dim][2]
                u   = bdcon[dim][3]
                v   = bdcon[dim][4]
                T   = bdcon[dim][6+nch]
                p   = rho*Rhat[i,bdy_nei[dim-2]]*T
                for ns in 1:nch
                    ch[ns] =  bdcon[dim][5+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
        elseif Int(bdcon[dim][1]) == 02
            temp1 = Int64(round(cellxmax/4)+1) # 1/4流出
            temp2 = Int64((temp1-1)*3+1)       # 1/4~3/4流入
            for i in 1:temp1
                rho = Qbase[i,bdy_nei[dim-2],1]
                u   = Qbase[i,bdy_nei[dim-2],2]
                v   = Qbase[i,bdy_nei[dim-2],3]
                p   = Qbase[i,bdy_nei[dim-2],4]
                for ns in 1:nch
                    ch[ns] =  Qbase[i,bdy_nei[dim-2],4+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
            for i in temp1+1:temp2
                rho = bdcon[dim][2]
                u   = bdcon[dim][3]
                v   = bdcon[dim][4]
                p   = bdcon[dim][5]
                for ns in 1:nch
                    ch[ns] =  bdcon[dim][5+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
            for i in temp2+1:cellxmax
                rho = Qbase[i,bdy_nei[dim-2],1]
                u   = Qbase[i,bdy_nei[dim-2],2]
                v   = Qbase[i,bdy_nei[dim-2],3]
                p   = Qbase[i,bdy_nei[dim-2],4]
                for ns in 1:nch
                    ch[ns] =  Qbase[i,bdy_nei[dim-2],4+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
        elseif Int(bdcon[dim][1]) == 03
            temp1 = Int64(round(cellxmax/4)+1) # 1/4流出
            temp2 = Int64((temp1-1)*3+1)       # 1/4~3/4流入
            for i in 1:temp1
                rho = Qbase[i,bdy_nei[dim-2],1]
                u   = Qbase[i,bdy_nei[dim-2],2]
                v   = Qbase[i,bdy_nei[dim-2],3]
                p   = Qbase[i,bdy_nei[dim-2],4]
                for ns in 1:nch
                    ch[ns] =  Qbase[i,bdy_nei[dim-2],4+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
            for i in temp1+1:temp2
                rho = bdcon[dim][2]
                u   = bdcon[dim][3]
                v   = bdcon[dim][4]
                T   = bdcon[dim][6+nch]
                p   = rho*Rhat[i,bdy_nei[dim-2]]*T
                for ns in 1:nch
                    ch[ns] =  bdcon[dim][5+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
            for i in temp2+1:cellxmax
                rho = Qbase[i,bdy_nei[dim-2],1]
                u   = Qbase[i,bdy_nei[dim-2],2]
                v   = Qbase[i,bdy_nei[dim-2],3]
                p   = Qbase[i,bdy_nei[dim-2],4]
                for ns in 1:nch
                    ch[ns] =  Qbase[i,bdy_nei[dim-2],4+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
        elseif Int(bdcon[dim][1]) == 20
            for i in 1:cellxmax
                rho = Qbase[i,bdy_nei[dim-2],1]
                u   = Qbase[i,bdy_nei[dim-2],2]
                v   = Qbase[i,bdy_nei[dim-2],3]
                p   = Qbase[i,bdy_nei[dim-2],4]
                for ns in 1:nch
                    ch[ns] =  Qbase[i,bdy_nei[dim-2],4+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
        elseif Int(bdcon[dim][1]) == 30
            for i in 1:cellxmax
                rho = Qbase[i,bdy_nei[dim-2],1]
                u   = -Qbase[i,bdy_nei[dim-2],2]
                v   = -Qbase[i,bdy_nei[dim-2],3]
                p   = Qbase[i,bdy_nei[dim-2],4]
                for ns in 1:nch
                    ch[ns] =  Qbase[i,bdy_nei[dim-2],4+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
        elseif Int(bdcon[dim][1]) == 31
            for i in 1:cellxmax
                u   = -Qbase[i,bdy_nei[dim-2],2]
                v   = -Qbase[i,bdy_nei[dim-2],3]
                p   = Qbase[i,bdy_nei[dim-2],4]
                T   = bdcon[dim][6+nch]
                rho = p/(Rhat[i,bdy_nei[dim-2]]*T)
                for ns in 1:nch
                    ch[ns] =  Qbase[i,bdy_nei[dim-2],4+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
        elseif Int(bdcon[dim][1]) == 40
            for i in 1:cellxmax
                rho = Qbase[i,bdy_nei_rev[dim-2],1]
                u   = Qbase[i,bdy_nei_rev[dim-2],2]
                v   = Qbase[i,bdy_nei_rev[dim-2],3]
                p   = Qbase[i,bdy_nei_rev[dim-2],4]
                for ns in 1:nch
                    ch[ns] = Qbase[i,bdy_nei_rev[dim-2],4+ns]
                end
                set_bdy(Qbase, rho, u, v, p, ch, nch, bdy[dim-2],i)
            end
        end
    end

    return Qbase
end

function check_bd(bdcon)
    bd = -1
    for l in 1:4
        if     Int(bdcon[l][1]) == 00
            bd = 00
        elseif Int(bdcon[l][1]) == 01
            bd = 01
        elseif Int(bdcon[l][1]) == 02
            bd = 02
        elseif Int(bdcon[l][1]) == 03
            bd = 03

        elseif Int(bdcon[l][1]) == 20
            bd = 20

        elseif Int(bdcon[l][1]) == 30
            bd = 30
        elseif Int(bdcon[l][1]) == 31
            bd = 31

        elseif Int(bdcon[l][1]) == 40
            bd = 40
        else
            println("\n check boundary condition ! \n")
            throw(UndefVarError(:x))
        end
        println("bd condition")
        println(bd)
    end
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
function readinput(filename)
  input = readlines("input")
  midpoint = findfirst(isequal(""), input)
  p1 = map(n -> parse(Int, n), input[2 : midpoint - 1])
  p2 = map(n -> parse(Int, n), input[midpoint + 2 : end])
  [p1, p2]
end

function combat(decks)
  queues = [copy(deck) for deck=decks]
  while !any(isempty.(queues))
    cards = [popfirst!(q) for q=queues]
    winner = queues[argmax(cards)]
    for card in sort(cards, rev=true)
      push!(winner, card)
    end
  end
  winner = argmax(length.(queues))
  winner, queues[winner]
end

function recursivecombat(decks, game = 1)
  queues = [copy(deck) for deck=decks]
  rounds = Set()
  round = 1
  while !any(isempty.(queues))
    round += 1
    if queues in rounds
      return 1, queues[1]
    end
    push!(rounds, deepcopy(queues))
    cards = [popfirst!(q) for q=queues]
    if all(cards .<= length.(queues))
      recdecks = [ queues[i][1:cards[i]] for i=1:length(queues)]
      windex, _ = recursivecombat(recdecks, game + 1)
    else
      windex = argmax(cards)
    end
    push!(queues[windex], cards[windex])

    for loser in cards[1:end .!=windex]
      push!(queues[windex], loser)
    end
  end
  winner = argmax(length.(queues))
  winner, queues[winner]
end

function deckvalue(deck)
  decklength = length(deck)
  sum([ (decklength - i + 1) * deck[i] for i=1:decklength ])
end

decks = readinput("input")

v1, d1 = combat(decks)
v1value = deckvalue(d1)
println("Part 1: $v1value")

v2, d2 = recursivecombat(decks)
v2value = deckvalue(d2)
println("Part 2: $v2value")


text:
using Documenter, LinearCovarianceModels

makedocs(
    sitename = "LinearCovarianceModels",
    pages = [
        "Introduction" => "index.md",
        ],
    strict=true)

deploydocs(
    repo = "github.com/saschatimme/LinearCovarianceModels.jl.git"
)


text:
module LinearSystems

println("\nLinearSystems: Demonstrates type and size deduction.")

using Modia
using Modia.Blocks: ABCD, Switch
using ModiaMath: plot

@model MySISOABCD begin
  @extends ABCD(x=Float(start=[0.0; 0.0]), A=[-1 -2; 0.5 -2], B=[1; 1], C=[1 1], D=[0], u=1.0)
  @inherits u
@equations begin
  # u = 1.0 # Does not work since size deduction finds the wrong "solution" which is not consistent (inv finds minimum norm solution).
  end
end

result = checkSimulation(MySISOABCD, 10, "x", 0.5000001198147023, logTranslation=true) #logOnFile = true)
plot(result, "x", heading="MySISOABCD", figure=13)

@model MyMIMOABCD begin
  s=Switch(u1=[1; 1], u2=[2; 3])
  lambda=0.5
  omega=10
  @extends ABCD(x=Float(start=[0.0; 0.0]), A=[-2*lambda omega^2; -1 0], B=[1 0; 0 1], C=[1 0; 0 1], D=zeros(2,2))
  @inherits u
@equations begin
  s.sw = time>5
  u = s.y
  end
end

result = checkSimulation(MyMIMOABCD, 10, "x", 0.004778288940817275, storeEliminated=false) # storeEliminated=false needed. Investigate
plot(result, "x", heading="MyMIMOABCD with generic switch", figure=14)

end


text:
using HypothesisTests, Test
using HypothesisTests: default_tail

@testset "Mann-Whitney" begin
@testset "Basic exact test" begin
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([1:10;], [2.1:2:21;]))) - 0.0232) <= 1e-4
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([2.1:2:21;], [1:10;]))) - 0.0232) <= 1e-4
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([1.5:10:100;], [2.1:2:21;]))) - 0.0068) <= 1e-4
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([2.1:2:21;], [1.5:10:100;]))) - 0.0068) <= 1e-4
	@test default_tail(ExactMannWhitneyUTest([1:10;], [2.1:2:21;])) == :both
	show(IOBuffer(), ExactMannWhitneyUTest([1:10;], [2.1:2:21;]))
end

@testset "Exact with ties" begin
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([1:10;], [1:10;]))) - 1) <= 1e-4
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([1:10;], [2:11;]))) - 0.5096) <= 1e-4
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([2:11;], [1:10;]))) - 0.5096) <= 1e-4
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([1:10;], [1:5; ones(5)]))) - 0.0057) <= 1e-4
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([1:5; ones(5)], [1:10;]))) - 0.0057) <= 1e-4
	show(IOBuffer(), ExactMannWhitneyUTest([1:10;], [1:10;]))
end

@testset "Exact with ties and unequal lengths" begin
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([1:10;], [2:2:24;]))) - 0.0118) <= 1e-4
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest([2:2:24;], [1:10;]))) - 0.0118) <= 1e-4
	show(IOBuffer(), ExactMannWhitneyUTest([1:10;], [2:2:24;]))
end

@testset "Approximate test" begin
    @test abs(@inferred(pvalue(ApproximateMannWhitneyUTest([1:10;], [2.1:2:21;]))) - 0.0257) <= 1e-4
    @test abs(@inferred(pvalue(ApproximateMannWhitneyUTest([2.1:2:21;], [1:10;]))) - 0.0257) <= 1e-4
    @test abs(@inferred(pvalue(ApproximateMannWhitneyUTest([1.5:10:100;], [2.1:2:21;]))) - 0.0091) <= 1e-4
    @test abs(@inferred(pvalue(ApproximateMannWhitneyUTest([2.1:2:21;], [1.5:10:100;]))) - 0.0091) <= 1e-4
    @test default_tail(ApproximateMannWhitneyUTest([1:10;], [2.1:2:21;])) == :both
	show(IOBuffer(), ApproximateMannWhitneyUTest([1:10;], [2.1:2:21;]))
end

@testset "Approximate with ties" begin
    @test abs(@inferred(pvalue(ApproximateMannWhitneyUTest([1:10;], [1:10;]))) - 1) <= 1e-4
    @test abs(@inferred(pvalue(ApproximateMannWhitneyUTest([1:10;], [2:11;]))) - 0.4948) <= 1e-4
    @test abs(@inferred(pvalue(ApproximateMannWhitneyUTest([2:11;], [1:10;]))) - 0.4948) <= 1e-4
    @test abs(@inferred(pvalue(ApproximateMannWhitneyUTest([1:10;], [1:5; ones(5)]))) - 0.0076) <= 1e-4
    @test abs(@inferred(pvalue(ApproximateMannWhitneyUTest([1:5; ones(5)], [1:10;]))) - 0.0076) <= 1e-4
	show(IOBuffer(), ApproximateMannWhitneyUTest([1:10;], [1:10;]))
end

@testset "Tests for automatic selection" begin
    @test abs(@inferred(pvalue(MannWhitneyUTest([1:10;], [2.1:2:21;]))) - 0.0232) <= 1e-4
    @test abs(@inferred(pvalue(MannWhitneyUTest([1:10;], [2:11;]))) - 0.4948) <= 1e-4
	show(IOBuffer(), MannWhitneyUTest([1:10;], [2.1:2:21;]))
end

@testset "Issue #39" begin
    @test abs(@inferred(pvalue(ExactMannWhitneyUTest(Float32[1:10;], Float32[2:11;]))) - 0.5096) <= 1e-4
end

@testset "Issue #126 pvalues above 1" begin
	A = [1.34937, 1.75722,0.276514, 1.04546, 1.69085, 0.738085, 2.36313]
	B = [2.62325, 1.16533, 1.1327, 0.728714]

	m = ExactMannWhitneyUTest(A,B)
	p = @inferred(pvalue(m; tail = :both))
	@test p == 1

	A = [12,10,7,6,3,1]
	B = [11,9,8,5,4,2]

	m = MannWhitneyUTest(A,B)
	p = @inferred(pvalue(m; tail = :both))

	@test p == 1

	m = ApproximateMannWhitneyUTest(A, B)
	p = @inferred(pvalue(m; tail = :both))
end

end


text:
#
# Confidence Interval Tools

abstract type ConfidenceInterval end

Base.broadcastable(ci::ConfidenceInterval) = Ref(ci)

# LowerUpperConfidenceInterval
Base.@kwdef struct LowerUpperConfidenceInterval <: ConfidenceInterval
    lower::Float64
    upper::Float64
    α::Float64 = 0.05
    estimate::Float64 = (lower + upper)/2
    target = nothing
    method = nothing
end

function Base.show(io::IO, ci::LowerUpperConfidenceInterval)
    print(io, "lower = ", round(ci.lower,sigdigits=4))
    print(io, ", upper = ", round(ci.upper,sigdigits=4))
    print(io, ", α = ", ci.α)
    print(io, "  (", ci.target,")")
end

@recipe function f(bands::AbstractVector{<:ConfidenceInterval})
    x = [Float64(location(band.target)) for band in bands]
    x, bands
end

@recipe function f(x, bands::AbstractVector{<:ConfidenceInterval};  show_ribbon=true)
    lower = getproperty.(bands, :lower)
    upper = getproperty.(bands, :upper)
    estimate = getproperty.(bands, :estimate)

	background_color_legend --> :transparent
	foreground_color_legend --> :transparent
    grid --> nothing
    framestyle --> :box
    legend --> :topleft

    if !show_ribbon
        _label = get(plotattributes, :label, "CI")
        _linestyle = get(plotattributes, :linestyle, :dash)
        _linecolor = get(plotattributes, :linecolor, :black)
        label := [_label nothing]
        linecolor := [_linecolor _linecolor]
        linestyle := [_linestyle _linestyle]
        return x, [lower upper]
    else
	    cis_ribbon  = estimate .- lower, upper .- estimate
	    fillalpha --> 0.36
	    seriescolor --> "#018AC4"
	    ribbon --> cis_ribbon
        linealpha --> 0
        return x, estimate
    end
end

function gaussian_ci(se; maxbias=0.0, α=0.05)
    if iszero(se)
        return maxbias
    end
    maxbias = abs(maxbias)
    rel_bias = maxbias/se
    if abs(rel_bias) > 6
        pm = quantile(Normal(), 1-α) + abs(rel_bias)
    else
        pm = sqrt(quantile(NoncentralChisq(1, abs2(rel_bias)), 1-α))
    end
    se*pm
end

Base.@kwdef struct BiasVarianceConfidenceInterval <: ConfidenceInterval
    target        = nothing
    method        = nothing
    α::Float64       = 0.05
    estimate::Float64
    se::Float64
    maxbias::Float64 = 0.0
    halflength::Float64 = gaussian_ci(se; maxbias=maxbias, α=α)
    lower::Float64 = estimate - halflength
    upper::Float64 = estimate + halflength
end

function Base.show(io::IO, ci::BiasVarianceConfidenceInterval)
    print(io, "lower = ", round(ci.lower,sigdigits=4))
    print(io, ", upper = ", round(ci.upper,sigdigits=4))
    print(io, ", est. = ", round(ci.estimate,sigdigits=4))
    print(io, ", se = ", round(ci.se,sigdigits=4))
    print(io, ", maxbias = ", round(ci.maxbias,sigdigits=4))
    print(io, ", α = ", ci.α)
    print(io, "  (", ci.target,")")
end

Base.@kwdef struct BisectionPair
    var1::Float64
    max_bias1::Float64
    estimate1::Float64
    var2::Float64
    max_bias2::Float64
    estimate2::Float64
    cov::Float64
end

Base.broadcastable(pair::BisectionPair) = Ref(pair)

function confint(pair::BisectionPair, λ ; α=0.05)
    _se = sqrt( abs2(1-λ)*pair.var1 + abs2(λ)*pair.var2 + 2*λ*(1-λ)*pair.cov)
    _maxbias = (1-λ)*pair.max_bias1 + λ*pair.max_bias2
    _estimate = (1-λ)*pair.estimate1 + λ*pair.estimate2
    bw = gaussian_ci(_se; maxbias=_maxbias, α=α)
    _estimate -bw , _estimate +bw
end


text:
function run_list_spark_glvprk(idae)
    run_list(idae, :TableauSPARKGLVPRK, (
            ( TableauSPARKGLVPRK(1),     "spark_glvprk1" ),
            ( TableauSPARKGLVPRK(2),     "spark_glvprk2" ),
            ( TableauSPARKGLVPRK(3),     "spark_glvprk3" ),
            ( TableauSPARKGLVPRK(4),     "spark_glvprk4" ),
        ))
end

function run_list_spark_glrk(idae)
    run_list(idae, :TableauSPARKGLRK, (
            ( TableauSPARKGLRK(1),       "spark_glrk1" ),
            ( TableauSPARKGLRK(2),       "spark_glrk2" ),
            ( TableauSPARKGLRK(3),       "spark_glrk3" ),
            ( TableauSPARKGLRK(4),       "spark_glrk4" ),
        ))
end

function run_list_spark_lobatto(idae)
    run_list(idae, :TableauSPARKLobatto, (
            ( TableauSPARKLobABC(2),     "spark_lobIIIABC2"    ),
            ( TableauSPARKLobABC(3),     "spark_lobIIIABC3"    ),
            ( TableauSPARKLobABC(4),     "spark_lobIIIABC4"    ),
            ( TableauSPARKLobABC(5),     "spark_lobIIIABC5"    ),

            ( TableauSPARKLobABD(2),     "spark_lobIIIABD2"    ),
            ( TableauSPARKLobABD(3),     "spark_lobIIIABD3"    ),
            ( TableauSPARKLobABD(4),     "spark_lobIIIABD4"    ),
            ( TableauSPARKLobABD(5),     "spark_lobIIIABD5"    ),
        ))
end

function run_list_vspark_internal_projection(idae)
    run_list(idae, :TableauVSPARKInternalProjection, (
            ( TableauVSPARKGLRKpInternal(1),                    "vspark_pInternal_glrk1" ),
            ( TableauVSPARKGLRKpInternal(2),                    "vspark_pInternal_glrk2" ),
            ( TableauVSPARKGLRKpInternal(3),                    "vspark_pInternal_glrk3" ),
            ( TableauVSPARKGLRKpInternal(4),                    "vspark_pInternal_glrk4" ),
        ); phi_average = Φ -> sum(Φ))
end

function run_list_vspark_modified_internal_projection(idae)
    run_list(idae, :TableauVSPARKModifiedInternalProjection, (
            ( TableauVSPARKGLRKpModifiedInternal(1),            "vspark_pModifiedInternal_glrk1" ),
            ( TableauVSPARKGLRKpModifiedInternal(2),            "vspark_pModifiedInternal_glrk2" ),
            ( TableauVSPARKGLRKpModifiedInternal(3),            "vspark_pModifiedInternal_glrk3" ),
            ( TableauVSPARKGLRKpModifiedInternal(4),            "vspark_pModifiedInternal_glrk4" ),
        ); phi_average = Φ -> sum(Φ))
end

function run_list_vspark_lobatto_IIIAIIIB_projection(idae)
    run_list(idae, :TableauVSPARKLobattoIIIAIIIBProjection, (
            ( TableauVSPARKGLRKpLobattoIIIAIIIB(1),                 "vspark_pLobattoIIIAIIIB_glrk1"    ),
            ( TableauVSPARKGLRKpLobattoIIIAIIIB(2),                 "vspark_pLobattoIIIAIIIB_glrk2"    ),
            ( TableauVSPARKGLRKpLobattoIIIAIIIB(3),                 "vspark_pLobattoIIIAIIIB_glrk3"    ),
            ( TableauVSPARKGLRKpLobattoIIIAIIIB(4),                 "vspark_pLobattoIIIAIIIB_glrk4"    ),

            ( TableauVSPARKLobattoIIIAIIIBpLobattoIIIAIIIB(2),      "vspark_pLobattoIIIAIIIB_lobIIIAB2" ),
            ( TableauVSPARKLobattoIIIAIIIBpLobattoIIIAIIIB(3),      "vspark_pLobattoIIIAIIIB_lobIIIAB3" ),
            ( TableauVSPARKLobattoIIIAIIIBpLobattoIIIAIIIB(4),      "vspark_pLobattoIIIAIIIB_lobIIIAB4" ),
            ( TableauVSPARKLobattoIIIAIIIBpLobattoIIIAIIIB(5),      "vspark_pLobattoIIIAIIIB_lobIIIAB5" ),

            ( TableauVSPARKLobattoIIIBIIIApLobattoIIIAIIIB(2),      "vspark_pLobattoIIIAIIIB_lobIIIBA2" ),
            ( TableauVSPARKLobattoIIIBIIIApLobattoIIIAIIIB(3),      "vspark_pLobattoIIIAIIIB_lobIIIBA3" ),
            ( TableauVSPARKLobattoIIIBIIIApLobattoIIIAIIIB(4),      "vspark_pLobattoIIIAIIIB_lobIIIBA4" ),
            ( TableauVSPARKLobattoIIIBIIIApLobattoIIIAIIIB(5),      "vspark_pLobattoIIIAIIIB_lobIIIBA5" ),
        ))
end

function run_list_vspark_lobatto_IIIBIIIA_projection(idae)
    run_list(idae, :TableauVSPARKLobattoIIIBIIIAProjection, (
            ( TableauVSPARKGLRKpLobattoIIIBIIIA(1),                 "vspark_pLobattoIIIBIIIA_glrk1"     ),
            ( TableauVSPARKGLRKpLobattoIIIBIIIA(2),                 "vspark_pLobattoIIIBIIIA_glrk2"     ),
            ( TableauVSPARKGLRKpLobattoIIIBIIIA(3),                 "vspark_pLobattoIIIBIIIA_glrk3"     ),
            ( TableauVSPARKGLRKpLobattoIIIBIIIA(4),                 "vspark_pLobattoIIIBIIIA_glrk4"     ),

            ( TableauVSPARKLobattoIIIAIIIBpLobattoIIIBIIIA(2),      "vspark_pLobattoIIIBIIIA_lobIIIAB2" ),
            ( TableauVSPARKLobattoIIIAIIIBpLobattoIIIBIIIA(3),      "vspark_pLobattoIIIBIIIA_lobIIIAB3" ),
            ( TableauVSPARKLobattoIIIAIIIBpLobattoIIIBIIIA(4),      "vspark_pLobattoIIIBIIIA_lobIIIAB4" ),
            ( TableauVSPARKLobattoIIIAIIIBpLobattoIIIBIIIA(5),      "vspark_pLobattoIIIBIIIA_lobIIIAB5" ),

            ( TableauVSPARKLobattoIIIBIIIApLobattoIIIBIIIA(2),      "vspark_pLobattoIIIBIIIA_lobIIIBA2" ),
            ( TableauVSPARKLobattoIIIBIIIApLobattoIIIBIIIA(3),      "vspark_pLobattoIIIBIIIA_lobIIIBA3" ),
            ( TableauVSPARKLobattoIIIBIIIApLobattoIIIBIIIA(4),      "vspark_pLobattoIIIBIIIA_lobIIIBA4" ),
            ( TableauVSPARKLobattoIIIBIIIApLobattoIIIBIIIA(5),      "vspark_pLobattoIIIBIIIA_lobIIIBA5" ),
        ))
end

function run_list_vspark_modified_lobatto_IIIAIIIB_projection(idae)
    run_list(idae, :TableauVSPARKModifiedLobattoIIIAIIIBProjection, (
            ( TableauVSPARKGLRKpModifiedLobattoIIIAIIIB(1),             "vspark_pModifiedLobattoIIIAIIIB_glrk1"     ),
            ( TableauVSPARKGLRKpModifiedLobattoIIIAIIIB(2),             "vspark_pModifiedLobattoIIIAIIIB_glrk2"     ),
            ( TableauVSPARKGLRKpModifiedLobattoIIIAIIIB(3),             "vspark_pModifiedLobattoIIIAIIIB_glrk3"     ),
            ( TableauVSPARKGLRKpModifiedLobattoIIIAIIIB(4),             "vspark_pModifiedLobattoIIIAIIIB_glrk4"     ),

            ( TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIAIIIB(2),  "vspark_pModifiedLobattoIIIAIIIB_lobIIIAB2" ),
            ( TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIAIIIB(3),  "vspark_pModifiedLobattoIIIAIIIB_lobIIIAB3" ),
            ( TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIAIIIB(4),  "vspark_pModifiedLobattoIIIAIIIB_lobIIIAB4" ),
            ( TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIAIIIB(5),  "vspark_pModifiedLobattoIIIAIIIB_lobIIIAB5" ),

            ( TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIAIIIB(2),  "vspark_pModifiedLobattoIIIAIIIB_lobIIIBA2" ),
            ( TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIAIIIB(3),  "vspark_pModifiedLobattoIIIAIIIB_lobIIIBA3" ),
            ( TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIAIIIB(4),  "vspark_pModifiedLobattoIIIAIIIB_lobIIIBA4" ),
            ( TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIAIIIB(5),  "vspark_pModifiedLobattoIIIAIIIB_lobIIIBA5" ),
        ))
end

function run_list_vspark_modified_lobatto_IIIBIIIA_projection(idae)
    run_list(idae, :TableauVSPARKModifiedLobattoIIIBIIIAProjection, (
            ( TableauVSPARKGLRKpModifiedLobattoIIIBIIIA(1),             "vspark_pModifiedLobattoIIIBIIIA_glrk1"     ),
            ( TableauVSPARKGLRKpModifiedLobattoIIIBIIIA(2),             "vspark_pModifiedLobattoIIIBIIIA_glrk2"     ),
            ( TableauVSPARKGLRKpModifiedLobattoIIIBIIIA(3),             "vspark_pModifiedLobattoIIIBIIIA_glrk3"     ),
            ( TableauVSPARKGLRKpModifiedLobattoIIIBIIIA(4),             "vspark_pModifiedLobattoIIIBIIIA_glrk4"     ),

            ( TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIBIIIA(2),  "vspark_pModifiedLobattoIIIBIIIA_lobIIIAB2" ),
            ( TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIBIIIA(3),  "vspark_pModifiedLobattoIIIBIIIA_lobIIIAB3" ),
            ( TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIBIIIA(4),  "vspark_pModifiedLobattoIIIBIIIA_lobIIIAB4" ),
            ( TableauVSPARKLobattoIIIAIIIBpModifiedLobattoIIIBIIIA(5),  "vspark_pModifiedLobattoIIIBIIIA_lobIIIAB5" ),

            ( TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIBIIIA(2),  "vspark_pModifiedLobattoIIIBIIIA_lobIIIBA2" ),
            ( TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIBIIIA(3),  "vspark_pModifiedLobattoIIIBIIIA_lobIIIBA3" ),
            ( TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIBIIIA(4),  "vspark_pModifiedLobattoIIIBIIIA_lobIIIBA4" ),
            ( TableauVSPARKLobattoIIIBIIIApModifiedLobattoIIIBIIIA(5),  "vspark_pModifiedLobattoIIIBIIIA_lobIIIBA5" ),
        ))
end

function run_list_vspark_midpoint_projection(idae)
    run_list(idae, :TableauVSPARKMidpointProjection, (
        ( TableauVSPARKGLRKpMidpoint(1),                            "vspark_pMidpoint_glrk1"     ),
        ( TableauVSPARKGLRKpMidpoint(2),                            "vspark_pMidpoint_glrk2"     ),
        ( TableauVSPARKGLRKpMidpoint(3),                            "vspark_pMidpoint_glrk3"     ),
        ( TableauVSPARKGLRKpMidpoint(4),                            "vspark_pMidpoint_glrk4"     ),

        ( TableauVSPARKLobattoIIIAIIIBpMidpoint(2),                 "vspark_pMidpoint_lobIIIAB2" ),
        ( TableauVSPARKLobattoIIIAIIIBpMidpoint(3),                 "vspark_pMidpoint_lobIIIAB3" ),
        ( TableauVSPARKLobattoIIIAIIIBpMidpoint(4),                 "vspark_pMidpoint_lobIIIAB4" ),
        ( TableauVSPARKLobattoIIIAIIIBpMidpoint(5),                 "vspark_pMidpoint_lobIIIAB5" ),

        ( TableauVSPARKLobattoIIIBIIIApMidpoint(2),                 "vspark_pMidpoint_lobIIIBA2" ),
        ( TableauVSPARKLobattoIIIBIIIApMidpoint(3),                 "vspark_pMidpoint_lobIIIBA3" ),
        ( TableauVSPARKLobattoIIIBIIIApMidpoint(4),                 "vspark_pMidpoint_lobIIIBA4" ),
        ( TableauVSPARKLobattoIIIBIIIApMidpoint(5),                 "vspark_pMidpoint_lobIIIBA5" ),
        ))
end

function run_list_vspark_modified_midpoint_projection(idae)
    run_list(idae, :TableauVSPARKModifiedMidpointProjection, (
        ( TableauVSPARKGLRKpModifiedMidpoint(1),                    "vspark_pModifiedMidpoint_glrk1"     ),
        ( TableauVSPARKGLRKpModifiedMidpoint(2),                    "vspark_pModifiedMidpoint_glrk2"     ),
        ( TableauVSPARKGLRKpModifiedMidpoint(3),                    "vspark_pModifiedMidpoint_glrk3"     ),
        ( TableauVSPARKGLRKpModifiedMidpoint(4),                    "vspark_pModifiedMidpoint_glrk4"     ),

        ( TableauVSPARKLobattoIIIAIIIBpModifiedMidpoint(2),         "vspark_pModifiedMidpoint_lobIIIAB2" ),
        ( TableauVSPARKLobattoIIIAIIIBpModifiedMidpoint(3),         "vspark_pModifiedMidpoint_lobIIIAB3" ),
        ( TableauVSPARKLobattoIIIAIIIBpModifiedMidpoint(4),         "vspark_pModifiedMidpoint_lobIIIAB4" ),
        ( TableauVSPARKLobattoIIIAIIIBpModifiedMidpoint(5),         "vspark_pModifiedMidpoint_lobIIIAB5" ),

        ( TableauVSPARKLobattoIIIBIIIApModifiedMidpoint(2),         "vspark_pModifiedMidpoint_lobIIIBA2" ),
        ( TableauVSPARKLobattoIIIBIIIApModifiedMidpoint(3),         "vspark_pModifiedMidpoint_lobIIIBA3" ),
        ( TableauVSPARKLobattoIIIBIIIApModifiedMidpoint(4),         "vspark_pModifiedMidpoint_lobIIIBA4" ),
        ( TableauVSPARKLobattoIIIBIIIApModifiedMidpoint(5),         "vspark_pModifiedMidpoint_lobIIIBA5" ),
        ))
end

function run_list_vspark_symmetric_projection(idae)
run_list(idae, :TableauVSPARKSymmetricProjection, (
        ( TableauVSPARKGLRKpSymmetric(1),                           "vspark_pSymmetric_glrk1"     ),
        ( TableauVSPARKGLRKpSymmetric(2),                           "vspark_pSymmetric_glrk2"     ),
        ( TableauVSPARKGLRKpSymmetric(3),                           "vspark_pSymmetric_glrk3"     ),
        ( TableauVSPARKGLRKpSymmetric(4),                           "vspark_pSymmetric_glrk4"     ),

        ( TableauVSPARKLobattoIIIAIIIBpSymmetric(2),                "vspark_pSymmetric_lobIIIAB2" ),
        ( TableauVSPARKLobattoIIIAIIIBpSymmetric(3),                "vspark_pSymmetric_lobIIIAB3" ),
        ( TableauVSPARKLobattoIIIAIIIBpSymmetric(4),                "vspark_pSymmetric_lobIIIAB4" ),
        ( TableauVSPARKLobattoIIIAIIIBpSymmetric(5),                "vspark_pSymmetric_lobIIIAB5" ),

        ( TableauVSPARKLobattoIIIBIIIApSymmetric(2),                "vspark_pSymmetric_lobIIIBA2" ),
        ( TableauVSPARKLobattoIIIBIIIApSymmetric(3),                "vspark_pSymmetric_lobIIIBA3" ),
        ( TableauVSPARKLobattoIIIBIIIApSymmetric(4),                "vspark_pSymmetric_lobIIIBA4" ),
        ( TableauVSPARKLobattoIIIBIIIApSymmetric(5),                "vspark_pSymmetric_lobIIIBA5" ),
    ))
end


text:
#testing the sbit state of the Gnum.
set_onesided!{ESS,FSS}(x::Gnum{ESS,FSS}) = (x.scratchpad.flags |= GNUM_ONESIDED_MASK; x)
set_twosided!{ESS,FSS}(x::Gnum{ESS,FSS}) = (x.scratchpad.flags &= ~GNUM_ONESIDED_MASK; x)
is_onesided{ESS,FSS}(x::Gnum{ESS,FSS}) = ((x.scratchpad.flags & GNUM_ONESIDED_MASK != 0) || (is_nan(x)))
is_twosided{ESS,FSS}(x::Gnum{ESS,FSS}) = ((x.scratchpad.flags & GNUM_ONESIDED_MASK == 0) && (!is_nan(x)))

#the ignore_side utility is used for operations that do identity checks before
#proceeding with calculations, and flags that a side has already been checked
#and should not be altered.  the `should_calculate` directive runs this downstream
@generated function ignore_side!{ESS,FSS,side}(x::Gnum{ESS,FSS}, ::Type{Val{side}})
  :(x.$side.flags |= GNUM_IGNORE_SIDE_MASK; nothing)
end
function ignore_both_sides!{ESS,FSS}(x::Gnum{ESS,FSS})
  x.lower.flags |= GNUM_IGNORE_SIDE_MASK
  x.upper.flags |= GNUM_IGNORE_SIDE_MASK
  nothing
end
function clear_ignore_sides!{ESS,FSS}(x::Gnum{ESS,FSS})
  x.lower.flags &= ~GNUM_IGNORE_SIDE_MASK
  x.upper.flags &= ~GNUM_IGNORE_SIDE_MASK
  nothing
end
#reports on whether or not the side referred to should be calculated.
@generated function should_calculate{ESS,FSS,side}(x::Gnum{ESS,FSS}, ::Type{Val{side}})
  if (side == :lower)
    :((x.lower.flags & GNUM_IGNORE_SIDE_MASK) == 0)
  else
    :(is_twosided(x) && ((x.upper.flags & GNUM_IGNORE_SIDE_MASK) == 0))
  end
end

function put_unum!{ESS,FSS}(src::Unum{ESS,FSS}, dest::Gnum{ESS,FSS})
  #puts a unum into the gnum, assuming it is going to be a single-sided unum.
  copy_unum!(src, dest.lower)
  set_flags!(dest, LOWER_UNUM)
  set_onesided!(dest)
  nothing
end

function get_unum!{ESS,FSS}(src::Gnum{ESS,FSS}, dest::Unum{ESS,FSS})
  #puts a unum into the gnum, assuming it already is a single-sided unum.
  (is_twosided(src) && !is_nan(src)) && throw(ArgumentError("Error: Gnum represents a Ubound"))
  is_nan(src) && return nan(Unum{ESS,FSS})
  force_from_flags!(src, dest, LOWER_UNUM) || copy_unum!(src.lower, dest)
  nothing
end

@generated function put_unum!{ESS,FSS,side}(src::Unum{ESS,FSS}, dest::Gnum{ESS,FSS}, ::Type{Val{side}})
  #sets the unum value on either side of the Gnum (or possibly the scratchpad.)
  quote
    copy_unum!(src, dest.$side)
    set_g_flags!(dest.$side)
    nothing
  end
end

@generated function get_unum!{ESS,FSS,side}(src::Gnum{ESS,FSS}, dest::Unum{ESS,FSS}, ::Type{Val{side}})
  #retrieves the unum value from either side of the Gnum (or possibly the scratchpad.)
  quote
    force_from_flags!(src, dest, Val{$side}) || copy_unum!(src.$side, dest);
    nothing
  end
end

function put_ubound!{ESS,FSS}(src::Ubound{ESS,FSS}, dest::Gnum{ESS,FSS})
  #fills the Gnum data from a source Ubound.
  copy_unum!(src.lower, dest.lower)
  copy_unum!(src.upper, dest.upper)
  set_flags!(dest)
  nothing
end

function get_ubound!{ESS,FSS}(src::Gnum{ESS,FSS}, dest::Ubound{ESS,FSS})
  #fills the Ubound data from a source Gnum.
  (is_nan(src) || is_onesided(src)) && throw(ArgumentError("Error:  Gnum represents a Unum"))
  #be sure to check if one of the flags is thrown before copying, otherwise
  #undefined results may occur.
  force_from_flags!(src, dest.lower, LOWER_UNUM) || copy_unum!(src.lower, dest.lower)
  force_from_flags!(src, dest.upper, UPPER_UNUM) || copy_unum!(src.upper, dest.upper)
  nothing
end

doc"""
`set_g_flags!(v::Unum{ESS,FSS})` sets the flags on any given unum to the appropriate
g-layer flag scheme.
"""
function set_g_flags!{ESS,FSS}(v::Unum{ESS,FSS})
  #clear the mask, keeping only flags and scratchpad values.
  is_inf(v) &&  (v.flags |= GNUM_INF_MASK; return)
  is_mmr(v) &&  (v.flags |= GNUM_MMR_MASK; return)
  is_sss(v) &&  (v.flags |= GNUM_SSS_MASK; return)
  is_zero(v) && (v.flags |= GNUM_ZERO_MASK; return)
  v.flags &= GNUM_SFLAGS_MASK | UNUM_FLAG_MASK
end

doc"""
  `set_flags!(::Gnum{ESS,FSS}, ::Type{Val{side}})` sets flags on one side of the
  gnum by examining the value.
"""
@generated function set_flags!{ESS,FSS,side}(v::Gnum{ESS,FSS}, ::Type{Val{side}})
  quote
    is_nan(v.$side) && (v.scratchpad.flags |= GNUM_NAN_MASK; return)
    #clear the mask, keeping only flags and scratchpad values.
    v.$side.flags &= GNUM_SFLAGS_MASK | UNUM_FLAG_MASK
    set_g_flags!(v.$side)
  end
end

function clear_gflags!{ESS,FSS}(target::Unum{ESS,FSS})
  target.flags &= ~(GNUM_FLAG_MASK)
end

function copy_unum_with_gflags!{ESS,FSS}(src::Unum{ESS,FSS}, dest::Unum{ESS,FSS})
  gflags = src.flags & GNUM_FLAG_MASK
  copy_unum!(src, dest)
  dest.flags |= gflags
end

function additive_inverse!{ESS,FSS}(target::Gnum{ESS,FSS})
  #lazy eval on this process.
  is_nan(target) && return

  if is_twosided(target)
    #swap the order of lower/upper via a buffer intermediate.
    copy_unum_with_gflags!(target.lower, target.buffer)
    copy_unum_with_gflags!(target.upper, target.lower)
    copy_unum_with_gflags!(target.buffer, target.upper)
    #next swap the parities on lower and upper
    target.upper.flags $= UNUM_SIGN_MASK
  end
  target.lower.flags $= UNUM_SIGN_MASK
end

doc"""
  `emit_data(::Gnum{ESS,FSS})` takes the contents of a gnum and decides if it's
  represents a solo unum or a ubound.  It then allocates the appropriate type and
  emits that as a result.
"""
function emit_data{ESS,FSS}(src::Gnum{ESS,FSS})
  #be ready to release a utype as a result.
  res::Utype
  #check to see if we're a NaN
  (is_nan(src)) && return nan(Unum{ESS,FSS})
  #check to see if we're a single unum
  if (is_onesided(src))
    #prepare the result by allocating.
    res = zero(Unum{ESS,FSS})
    #put the value in the allocated space.
    get_unum!(src, res)
    res
  else
    #this time, we know it's a ubound.
    #prepare the result by allocating.
    res = Ubound{ESS,FSS}()
    #put the value in the allocated space.
    get_ubound!(src, res)
    (res.lower == res.upper) ? res.lower : res
  end
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using SAASGP4, Compat

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

tmpdir = mktempdir()

@testset "DllMain" begin
    @testset "getLastErrMsg" begin
        @test SAASGP4.getLastErrMsg() isa AbstractString
        @test SAASGP4.getLastErrMsg() == ""
    end

    @testset "getLastInfoMsg" begin
        @test SAASGP4.getLastInfoMsg() isa AbstractString
        @test SAASGP4.getLastInfoMsg() == ""
    end

    @testset "mainGetInfo" begin
        @test SAASGP4.mainGetInfo() isa AbstractString
        @test length(SAASGP4.mainGetInfo()) > 0
    end

    @testset "getInitDllNames" begin
        @test SAASGP4.getInitDllNames() isa AbstractString
        @test SAASGP4.getInitDllNames() == "Initialized Dlls: DllMain,TimeFunc,EnvConst,AstroFunc,Tle,Sgp4Prop,"
    end

    @testset "openLogFile" begin
        @test_nowarn SAASGP4.openLogFile(joinpath(tmpdir, "log.txt"))
        @test isfile(joinpath(tmpdir, "log.txt"))
    end

    @testset "closeLogFile" begin
        @test_nowarn SAASGP4.closeLogFile()
    end
end

@testset "EnvConst" begin
    @testset "envGetInfo" begin
        @test SAASGP4.envGetInfo() isa AbstractString
        @test length(SAASGP4.envGetInfo()) > 0
    end

    @testset "envGetGeoStr" begin
        @test SAASGP4.envGetGeoStr() isa AbstractString
        @test SAASGP4.envGetGeoStr() == "WGS-72"
    end

    @testset "envSetGeoStr" begin
        @testset "value: $value" for value in ["WGS-84", "EGM-96", "WGS-72", "JGM2", "SEM68R", "GEM5", "GEM9"]
            SAASGP4.envSetGeoStr(value)
            @test SAASGP4.envGetGeoStr() == value
        end

        # Reset to WGS-72
        SAASGP4.envSetGeoStr("WGS-72")
        @test SAASGP4.envGetGeoStr() == "WGS-72"

        # Set to invalid value is no-op
        SAASGP4.envSetGeoStr("WGS-52")
        @test SAASGP4.envGetGeoStr() == "WGS-72"
    end

    @testset "envGetGeoIdx" begin
        @test SAASGP4.envGetGeoIdx() isa Integer
        @test SAASGP4.envGetGeoIdx() == 72
    end

    @testset "envSetGeoIdx" begin
        @testset "id: $id" for id in [84, 96, 72, 2, 68, 5, 9]
            SAASGP4.envSetGeoIdx(id)
            @test SAASGP4.envGetGeoIdx() == id
        end

        # Reset to 72
        SAASGP4.envSetGeoIdx(72)
        @test SAASGP4.envGetGeoIdx() == 72

        # Set to invalid number is no-op
        SAASGP4.envSetGeoIdx(1)
        @test SAASGP4.envGetGeoIdx() == 72
    end

    @testset "envGetGeoConst id: $id" for id in 1:11
        @test SAASGP4.envGetGeoConst(id) isa AbstractFloat
    end

    @testset "envGetFkIdx" begin
        @test SAASGP4.envGetFkIdx() isa Integer
        @test SAASGP4.envGetFkIdx() == 5
    end

    @testset "envSetFkIdx" begin
        SAASGP4.envSetFkIdx(4)
        @test SAASGP4.envGetFkIdx() == 4

        # Set to invalid number is no-op
        SAASGP4.envSetFkIdx(8)
        @test SAASGP4.envGetFkIdx() == 4

        # Reset to 5
        SAASGP4.envSetFkIdx(5)
        @test SAASGP4.envGetFkIdx() == 5
    end

    @testset "envGetFkConst id: $id" for id in 1:11
        @test SAASGP4.envGetFkConst(id) isa AbstractFloat
    end

    @testset "envGetFkPtr" begin
        @test SAASGP4.envGetFkPtr() isa Integer
        # TODO: Check if usable in ThetaGrnwch
    end

    @testset "envSaveFile" begin
        path = joinpath(tmpdir, "envFile.txt")
        @test_nowarn SAASGP4.envSaveFile(path)
        @test isfile(path)

        # Give empty file path
        @test_throws ErrorException SAASGP4.envSaveFile("")

        # Give impossible file path
        @test_throws ErrorException SAASGP4.envSaveFile(joinpath(tmpdir, "bogus", "envFile.txt"))
    end

    @testset "envLoadFile" begin
        # Check for actuall op to env constants
        @test_nowarn SAASGP4.envLoadFile(joinpath(@__DIR__, "wgs84fk4.txt"))
        @test SAASGP4.envGetFkIdx() == 4
        @test SAASGP4.envGetGeoIdx() == 84

        # Reset constants and check for file saved in `envSaveFile`
        @test_nowarn SAASGP4.envLoadFile(joinpath(tmpdir, "envFile.txt"))
        @test SAASGP4.envGetFkIdx() == 5
        @test SAASGP4.envGetGeoIdx() == 72

        # Give non-existant file path
        @test_throws ErrorException SAASGP4.envLoadFile(joinpath(tmpdir, "bogus.txt"))
    end
end

@testset "AstroFunc" begin
    @testset "astroFuncGetInfo" begin
        @test SAASGP4.astroFuncGetInfo() isa AbstractString
        @test length(SAASGP4.astroFuncGetInfo()) > 0
    end
end

@testset "TimeFunc" begin
    @testset "timeFuncGetInfo" begin
        @test SAASGP4.timeFuncGetInfo() isa AbstractString
        @test length(SAASGP4.timeFuncGetInfo()) > 0
    end

    @testset "DTGToUTC" begin
        dtg = "00051.47568104"
        ds50UTC = 18313.47568104
        @test SAASGP4.DTGToUTC(dtg) isa AbstractFloat
        @test SAASGP4.DTGToUTC(dtg) == ds50UTC
    end
end

@testset "Tle" begin
    @testset "tleGetInfo" begin
        @test SAASGP4.tleGetInfo() isa AbstractString
        @test length(SAASGP4.tleGetInfo()) > 0
    end

    line1 = "1 90021U RELEAS14 00051.47568104 +.00000184 +00000+0 +00000-4 0 0814"
    line2 = "2 90021   0.0222 182.4923 0000720  45.6036 131.8822  1.00271328 1199"
    @testset "tleAddSatFrLines" begin
        satkey = SAASGP4.tleAddSatFrLines(line1, line2)
        @test satkey isa Integer
        SAASGP4.tleRemoveSat(satkey)
        # TODO: Check if `satkey` can actually be initialized, or if count is augmented
    end

    @testset "tleRemoveSat" begin
        satkey = SAASGP4.tleAddSatFrLines(line1, line2)
        @test_nowarn SAASGP4.tleRemoveSat(satkey)
        # TODO: Check if `satkey` is actually removed

        @test_throws ErrorException SAASGP4.tleRemoveSat(42)
    end

    @testset "tleRemoveAllSats" begin
        SAASGP4.tleAddSatFrLines(line1, line2)
        @test_nowarn SAASGP4.tleRemoveAllSats()
        # TODO: Check if all sats are actually removed
    end

    @testset "tleGetCount" begin
        count = SAASGP4.tleGetCount()
        @test count isa Integer
        @test count == 0

        SAASGP4.tleAddSatFrLines(line1, line2)
        @test SAASGP4.tleGetCount() == 1
        SAASGP4.tleRemoveAllSats()
    end

    @testset "tleLoadFile" begin
        @test_nowarn SAASGP4.tleLoadFile(joinpath(@__DIR__, "verify/rel14.inp"))
        @test SAASGP4.tleGetCount() == 14
        SAASGP4.tleRemoveAllSats()
    end
end

@testset "Sgp4Prop" begin
    @testset "sgp4GetInfo" begin
        @test SAASGP4.sgp4GetInfo() isa AbstractString
        @test length(SAASGP4.sgp4GetInfo()) > 0
    end

    @testset "sgp4SetLicFilePath" begin
        @test_nowarn SAASGP4.sgp4SetLicFilePath(normpath(joinpath(@__DIR__, "../deps/usr/lib/")))
        # TODO: Test for op. This only affects the library if it's called before `sgp4Init`.
        #       But as the path is set during `__init__`, it is effectivly tested before this point.
    end

    line1 = "1 90021U RELEAS14 00051.47568104 +.00000184 +00000+0 +00000-4 0 0814"
    line2 = "2 90021   0.0222 182.4923 0000720  45.6036 131.8822  1.00271328 1199"
    satkey = SAASGP4.tleAddSatFrLines(line1, line2)

    @testset "sgp4InitSat" begin
        @test_nowarn SAASGP4.sgp4InitSat(satkey)
        # TODO: Check if `satkey` can actually be used to propagate

        @test_throws ErrorException SAASGP4.sgp4InitSat(42)
    end

    @testset "sgp4RemoveSat" begin
        @test_nowarn SAASGP4.sgp4RemoveSat(satkey)
        # TODO: Check if `satkey` is actually removed

        @test_throws ErrorException SAASGP4.sgp4RemoveSat(42)
    end

    @testset "sgp4RemoveAllSats" begin
        SAASGP4.sgp4InitSat(satkey)
        @test_nowarn SAASGP4.sgp4RemoveAllSats()
        # TODO: Check if all sats are actually removed
    end

    SAASGP4.sgp4InitSat(satkey)

    @testset "sgp4PropDs50UTC" begin
        @testset begin
            ds50UTC = 18313.47568104 + 1.0 / 1440

            pos, vel, llh, mse = SAASGP4.sgp4PropDs50UTC(satkey, ds50UTC)

            @test pos isa Array{T, 1} where T <: AbstractFloat
            @test length(pos) === 3
            @test vel isa Array{T, 1} where T <: AbstractFloat
            @test length(vel) === 3
            @test llh isa Array{T, 1} where T <: AbstractFloat
            @test length(llh) === 3
            # TODO: check values against references?

            @test mse isa AbstractFloat
            @test mse ≈ 1.0
        end

        @test_throws ErrorException SAASGP4.sgp4PropDs50UTC(42, 12345.678)
    end

    @testset "sgp4PropDs50UTC!" begin
        @testset begin
            ds50UTC = 18313.47568104 + 1.0 / 1440

            pos = Array{Cdouble, 1}(undef, 3)
            vel = Array{Cdouble, 1}(undef, 3)
            llh = Array{Cdouble, 1}(undef, 3)
            mse = Ref{Cdouble}(0.0)

            SAASGP4.sgp4PropDs50UTC!(pos, vel, llh, mse, satkey, ds50UTC)

            @test pos isa Array{T, 1} where T <: AbstractFloat
            @test length(pos) === 3
            @test vel isa Array{T, 1} where T <: AbstractFloat
            @test length(vel) === 3
            @test llh isa Array{T, 1} where T <: AbstractFloat
            @test length(llh) === 3
            # TODO: check values against references?

            @test mse[] isa AbstractFloat
            @test mse[] ≈ 1.0
        end

        @testset begin
            pos = Array{Cdouble, 1}(undef, 3)
            vel = Array{Cdouble, 1}(undef, 3)
            llh = Array{Cdouble, 1}(undef, 3)
            mse = Ref{Cdouble}(0.0)
            @test_throws ErrorException SAASGP4.sgp4PropDs50UTC!(pos, vel, llh, mse, 42, 12345.678)
        end
    end

    @testset "sgp4PropMse" begin
        @testset begin
            mse = 1.0

            pos, vel, llh, ds50UTC = SAASGP4.sgp4PropMse(satkey, mse)

            @test pos isa Array{T, 1} where T <: AbstractFloat
            @test length(pos) === 3
            @test vel isa Array{T, 1} where T <: AbstractFloat
            @test length(vel) === 3
            @test llh isa Array{T, 1} where T <: AbstractFloat
            @test length(llh) === 3
            # TODO: check values against references?

            @test ds50UTC isa AbstractFloat
            @test ds50UTC ≈ 18313.47568104 + 1.0 / 1440
        end

        @test_throws ErrorException SAASGP4.sgp4PropMse(42, 123.4567)
    end

    @testset "sgp4PropMse!" begin
        @testset begin
            mse = 1.0

            pos = Array{Cdouble, 1}(undef, 3)
            vel = Array{Cdouble, 1}(undef, 3)
            llh = Array{Cdouble, 1}(undef, 3)
            ds50UTC = Ref{Cdouble}(0.0)

            SAASGP4.sgp4PropMse!(pos, vel, llh, ds50UTC, satkey, mse)

            @test pos isa Array{T, 1} where T <: AbstractFloat
            @test length(pos) === 3
            @test vel isa Array{T, 1} where T <: AbstractFloat
            @test length(vel) === 3
            @test llh isa Array{T, 1} where T <: AbstractFloat
            @test length(llh) === 3
            # TODO: check values against references?

            @test ds50UTC[] isa AbstractFloat
            @test ds50UTC[] ≈ 18313.47568104 + 1.0 / 1440
        end

        @testset begin
            pos = Array{Cdouble, 1}(undef, 3)
            vel = Array{Cdouble, 1}(undef, 3)
            llh = Array{Cdouble, 1}(undef, 3)
            ds50UTC = Ref{Cdouble}(0.0)
            @test_throws ErrorException SAASGP4.sgp4PropMse!(pos, vel, llh, ds50UTC, 42, 123.4567)
        end
    end

    @testset "sgp4GetPropOut" begin
        tle1line1 = "1 90001U SGP4-VAL 93051.47568104 +.00000184  00000 0  00000-4 0 0814"
        tle1line2 = "2 90001   0.0221 182.4922 0000720  45.6036 131.8822  1.0027132801199"

        satkey1 = SAASGP4.tleAddSatFrLines(tle1line1, tle1line2)
        SAASGP4.sgp4InitSat(satkey1)

        @test_throws ErrorException SAASGP4.sgp4GetPropOut(satkey1, 0)

        SAASGP4.sgp4PropMse(satkey1, 1.0)

        data = SAASGP4.sgp4GetPropOut(satkey1, 1)
        @test data isa Array{T, 1} where T <: AbstractFloat
        @test length(data) === 1

        data = SAASGP4.sgp4GetPropOut(satkey1, 2)
        @test data isa Array{T, 1} where T <: AbstractFloat
        @test length(data) === 3

        data = SAASGP4.sgp4GetPropOut(satkey1, 3)
        @test data isa Array{T, 1} where T <: AbstractFloat
        @test length(data) === 6

        data = SAASGP4.sgp4GetPropOut(satkey1, 4)
        @test data isa Array{T, 1} where T <: AbstractFloat
        @test length(data) === 6

        @test_throws ErrorException SAASGP4.sgp4GetPropOut(satkey1, 0)
        @test_throws ErrorException SAASGP4.sgp4GetPropOut(satkey1, 5)

        # TODO: check values against references?
    end
end

rm(tmpdir, recursive = true)


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

JLLWrappers.@generate_main_file_header("Ogg")
JLLWrappers.@generate_main_file("Ogg", UUID("e7412a2a-1a6e-54c0-be00-318e2571c051"))
end  # module Ogg_jll


text:
# Making initiaizations for the environment and packages

println("Loading environment...")
sleep(1.5)
using Pkg

try
    Pkg.activate("fc2finalproject/")
    Pkg.status()
catch
    error("Environment was either not found or failed loading")
end

println("Loading packages...")

begin
    using Flux
    using Flux: @epochs
    using Plots, Statistics
    using ProgressMeter, Alert
    using BSON: @save
    using BSON: @load
    using DelimitedFiles
    gr()
end

println("Packages loaded successfully")
println("Loading auxiliary functions...")
sleep(1.5)
function tuple2vec(test::Tuple)
    return [test[1] , test[2]]
end

function generatetraindata(tuplemesh)
    row, col = size(tuplemesh)
    data = Matrix{Float64}(undef, 2, row*col)

    current = 0 # initialize
    for k in 1:col, m in 1:row
        current += 1
        data[:,current] = tuple2vec( tuplemesh[m,k] )

    end
    return data
end

println("Auxiliary functions loaded successfully")
println("Loading data...")
sleep(1.5)

intensities = readdlm("data.csv",',')

ydims, xdims = size(intensities)

sleep(1.5)
println("Your picture is of size $xdims px by $ydims px.")

sleep(1.5)
println("Generating training data...")

z_orthogonal_complement = [ (x,y) for y in 1:xdims, x in 1:ydims ]
Xtrain = generatetraindata(z_orthogonal_complement)
Ytrain = Float32.( Matrix( vec(intensities)' ) )
data = Flux.Data.DataLoader((Xtrain, Ytrain), batchsize=100,shuffle=true);

println("Data loaded successfully")
sleep(1.5)

goal = heatmap( unique(Xtrain[1,:]), unique(Xtrain[2,:]),
reverse(reshape(Ytrain, ydims,xdims)), c=:grays, clims = (0,1));
savefig("goal.pdf")

println("Finding last save state model.bson")
search = findall(x-> x.== "model.bson", readdir() )

if length(search) < 1

    println("Could not find a .bson file with a model")
    sleep(1.5)
    println("Creating model from scrach")
    sleep(1.5)

    function NeuralNetwork()
        return Chain(
                Dense(2, 200,relu),
                Dense(200, 100,relu),
                Dense(100, 50,relu),
                Dense(50, 25,relu),
                Dense(25, 10,relu),
                Dense(10,1,x->σ.(x))
                )
    end

    global model = NeuralNetwork()
    println("New model created successfully")
    sleep(1.5)

    im0 = heatmap( unique(Xtrain[1,:]), unique(Xtrain[2,:]),
    reverse(reshape(model(Xtrain), ydims,xdims)), c=:grays, clims = (0,1));
    savefig("History/initialguess.pdf")
else
    sleep(1.5)
    @load "model.bson"  model
    println("model.bson loaded successfully")
    sleep(1.5)

end

opt = Descent(0.05)
loss(x, y) = sum(Flux.Losses.binarycrossentropy(model(x), y))

save_state = 0; #initialize
total_iterations = 0; # initialize
ps = Flux.params(model);
epochs = 100000;

@showprogress 1 "Training..." for i in 1:epochs
    global save_state
    global total_iterations
    Flux.train!(loss, ps, data, opt)

    if (i % 5000) == 0
        # println("Saving... $i iterations have been completed")
        save_state += 1
        total_iterations += 5000
        @save "model.bson" model
        @save "Checkpoints/model_c$(save_state).bson" model
        image = heatmap( unique(Xtrain[1,:]), unique(Xtrain[2,:]),
        reverse(reshape(model(Xtrain), ydims,xdims)), c=:grays, clims = (0,1) );
        savefig("History/result$(total_iterations).pdf")
    end

end


text:
"""
Sieve of Eratosthenes is an algorithm for finding all the primes upto a limit `n`.

Reference:
-https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
"""
function eratosthenes(n)
    primes = fill(true,n)
    primes[1] = false
    for p = 2:n
        primes[p] || continue
        for i = 2:div(n,p)
            primes[p*i] = false
        end
    end
    findall(primes)
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia