text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related with the equinox-based model IAU-2006 with IAU-2010
#   conventions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] Capitaine, N., Wallace, P. T (2006). High precision methods for locating
#       the celestial intermediate pole and origin. Astronomy & Astrophysics.
#
#   [3] IERS (2010). Transformation between the International Terrestrial
#       Reference System and the Geocentric Celestial Reference System. IERS
#       Technical Note No. 36, Chapter 5.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export r_tirs_to_ers_iau2006, r_ers_to_tirs_iau2006,
       r_ers_to_mod_iau2006, r_mod_to_ers_iau2006,
       r_mod_to_mj2000_iau2006, r_mj2000_to_mod_iau2006,
       r_mj2000_to_gcrf_iau2006, r_gcrf_to_mj2000_iau2006,
       r_tirs_to_mod_iau2006, r_mod_to_tirs_iau2006

################################################################################
#                      IAU-2006 equinox-based reductions
################################################################################
#
# The conversion between the Geocentric Celestial Reference Frame (GCRF) to the
# International Terrestrial Reference Frame (ITRF) is done by means of:
#
#              GCRF <=> MJ2000 <=> MOD <=> ERS <=> TIRS <=> ITRF
#
# in which:
#   - TIRS: Terrestrial Intermediate Reference System.
#   - ERS: Earth Reference System.
#   - MOD: Mean of Date reference frame.
#   - MJ2000: J2000 mean equatorial frame (mean Equator, mean equinox dynamical
#             system).
#   - GCRF: Geocentric Celestial Reference Frame.
#
# Every rotation will be coded as a function using the IAU-2006 theory with the
# IAU-2010 conventions. The API is:
#
#   function r_<Origin Frame>_to_<Destination Frame>_iau2006
#
# The arguments vary depending on the origin and destination frame and should be
# verified using the function documentation.
#
################################################################################
#
################################################################################
#                               Single Rotations
################################################################################

#                                 TIRS <=> ERS
# ==============================================================================

"""
    r_tirs_to_ers_iau2006([T::Type,] JD_UT1::Number, JD_TT::Number, δΔΨ_2000::Number = 0)

Compute the rotation that aligns the Terrestrial Intermediate Reference System
(TIRS) with the Earth Reference System (ERS) at the Julian Day `JD_UT1` [UT1]
and `JD_TT` [Terrestrial Time]. This algorithm uses the IAU-2006 theory.

Notice that one can provide corrections for the nutation in longitude
(`δΔψ_2000`) [rad] that are usually obtained from IERS EOP Data (see
[`get_iers_eop`](@ref) and [`deps_dpsi`](@ref)). This corrections are related to
Free Core Nutation (FCN) that models the effect of a liquid Earth core.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TIRS frame with the ERS frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The reference frames TIRS and ERS are separated by a rotation about the Z-axis
of the Greenwhich apparent sidereal angle (GAST). This angle is computed using
the IAU-2006 theory, which consist of obtaining the Earth Rotation Angle (ERA)
and subtracting the result of the Equation of Origins (EO).
"""
function r_tirs_to_ers_iau2006(
    JD_UT1::Number,
    JD_TT::Number,
    δΔΨ_2000::Number = 0
)
    return r_tirs_to_ers_iau2006(DCM, JD_UT1, JD_TT)
end

function r_tirs_to_ers_iau2006(
    T::Type,
    JD_UT1::Number,
    JD_TT::Number,
    δΔΨ_2000::Number = 0
)
    # In this theory, the rotation of Earth is taken into account by the Earth
    # Rotation Angle, which is the angle between the Conventional International
    # Origin (CIO) and the Terrestrial Intermediate Origin (TIO) [1]. The latter
    # is a reference meridian on Earth that is located about 100m away from
    # Greenwich meridian along the equator of the Celestial Intermediate Pole
    # (CIP) [1].
    θ_era = 2π * (0.7790572732640 + 1.00273781191135448 * (JD_UT1 - JD_J2000))
    θ_era = mod(θ_era, 2π)

    # Compute the Equation of the Origins (EO).
    ~, ~, ~, EO = nutation_eo_iau2006(JD_TT, 0, δΔΨ_2000)

    # Compute the Greenwich apparent sidereal angle (GAST).
    θ_gast2000 = θ_era - EO

    # Compute the rotation between the TIRS and ERS.
    return angle_to_rot(T, -θ_gast2000, 0, 0, :ZYX)
end

"""
    r_ers_to_tirs_iau2006(JD_UT1::Number, JD_TT::Number, δΔΨ_2000::Number = 0)

Compute the rotation that aligns the Earth Reference System (ERS) with the
Terrestrial Intermediate Reference System (TIRS) at the Julian Day `JD_UT1`
[UT1] and `JD_TT` [Terrestrial Time]. This algorithm uses the IAU-2006 theory.

Notice that one can provide corrections for the nutation in longitude
(`δΔψ_2000`) [rad] that are usually obtained from IERS EOP Data (see
[`get_iers_eop`](@ref) and [`deps_dpsi`](@ref)). This corrections are related to
Free Core Nutation (FCN) that models the effect of a liquid Earth core.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the ERS frame with the TIRS frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The reference frames TIRS and ERS are separated by a rotation about the Z-axis
of the Greenwhich apparent sidereal angle (GAST). This angle is computed using
the IAU-2006 theory, which consist of obtaining the Earth Rotation Angle (ERA)
and subtracting the result of the Equation of Origins (EO).
"""
function r_ers_to_tirs_iau2006(
    JD_UT1::Number,
    JD_TT::Number,
    δΔΨ_2000::Number = 0
)
    return r_ers_to_tirs_iau2006(DCM, JD_UT1, JD_TT)
end

function r_ers_to_tirs_iau2006(
    T::Type,
    JD_UT1::Number,
    JD_TT::Number,
    δΔΨ_2000::Number = 0
)
    return inv_rotation(r_tirs_to_ers_iau2006(T, JD_UT1, JD_TT, δΔΨ_2000))
end

#                                ERS <=> MOD
# ==============================================================================

"""
    r_ers_to_mod_iau2006([T::Type,] JD_TT::Number, δΔϵ_2000::Number = 0, δΔΨ_2000::Number = 0)

Compute the rotation that aligns the Earth Reference System (ERS) with the
Mean of Date (MOD) reference frame at Julian day `JD_TT` [Terrestrial Time].
This algorithm uses the IAU-2006 theory.

Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_2000`) and in longitude (`δΔψ_2000`) [rad] that are usually obtained from
IERS EOP Data (see [`get_iers_eop`](@ref) and [`deps_dpsi`](@ref)). This
corrections are related to Free Core Nutation (FCN) that models the effect of a
liquid Earth core.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

!!! info
    The reference systems ERS and MOD are separated by the nutation of the pole.

# Returns

The rotation that aligns the ERS frame with the MOD frame. The rotation
representation is selected by the optional parameter `T`.
"""
function r_ers_to_mod_iau2006(
    JD_TT::Number,
    δΔϵ_2000::Number = 0,
    δΔΨ_2000::Number = 0
)
    return r_ers_to_mod_iau2006(DCM, JD_TT, δΔϵ_2000, δΔΨ_2000)
end

function r_ers_to_mod_iau2006(
    T::Type,
    JD_TT::Number,
    δΔϵ_2000::Number = 0,
    δΔΨ_2000::Number = 0
)
    # Compute the angles used to compute the nutation.
    mϵ_2000, Δϵ_2000, ΔΨ_2000, ~ = nutation_eo_iau2006(JD_TT, δΔϵ_2000, δΔΨ_2000)
    return angle_to_rot(T, mϵ_2000 + Δϵ_2000, ΔΨ_2000, -mϵ_2000, :XZX)
end

"""
    r_mod_to_ers_iau2006([T::Type,] JD_TT::Number, δΔϵ_2000::Number = 0, δΔΨ_2000::Number = 0)

Compute the rotation that aligns the Mean of Date (MOD) reference frame with the
Earth Reference System (ERS) at Julian day `JD_TT` [Terrestrial Time]. This
algorithm uses the IAU-2006 theory.

Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_2000`) and in longitude (`δΔψ_2000`) [rad] that are usually obtained from
IERS EOP Data (see [`get_iers_eop`](@ref) and [`deps_dpsi`](@ref)). This
corrections are related to Free Core Nutation (FCN) that models the effect of a
liquid Earth core.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the MOD frame with the ERS frame. The rotation
representation is selected by the optional parameter `T`.
"""
function r_mod_to_ers_iau2006(
    JD_TT::Number,
    δΔϵ_2000::Number = 0,
    δΔΨ_2000::Number = 0
)
    return r_mod_to_ers_iau2006(DCM, JD_TT, δΔϵ_2000, δΔΨ_2000)
end

function r_mod_to_ers_iau2006(
    T::Type,
    JD_TT::Number,
    δΔϵ_2000::Number = 0,
    δΔΨ_2000::Number = 0
)
    return inv_rotation(r_ers_to_mod_iau2006(T, JD_TT, δΔϵ_2000, δΔΨ_2000))
end

#                                MOD <=> MJ2000
# ==============================================================================

"""
    r_mod_to_mj2000_iau2006([T::Type,] JD_TT::Number)

Compute the rotation that aligns the Mean of Date (MOD) reference frame with the
J2000 mean equatorial frame at Julian day `JD_TT` [Terrestrial Time]. This
algorithm uses the IAU-2006 theory.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the MOD frame with the MJ2000 frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The J2000 reference frame here is not equal to the previous definition in FK5
theory. It is the reason why it is internally called `MJ2000`. According to
- **[3]**:

> The mean equinox of J2000.0 to be considered is not the “rotational dynamical
> mean equinox of J2000.0” as used in the past, but the “inertial dynamical mean
> equinox of J2000.0” to which the recent numerical or analytical solutions
> refer.  The latter is associated with the ecliptic in the inertial sense,
> which is the plane perpendicular to the angular momentum vector of the orbital
> motion of the Earth-Moon barycenter as computed from the velocity of the
> barycenter relative to an inertial system. The rotational equinox is
> associated with the ecliptic in the rotational sense, which is perpendicular
> to the angular momentum vector computed from the velocity referred to the
> rotating orbital plane of the Earth-Moon barycenter. (The difference between
> the two angular momenta is the angular momentum associated with the rotation
> of the orbital plane.)

# References

- **[3]**: IERS (2010). Transformation between the International Terrestrial
    Reference System and the Geocentric Celestial Reference System. IERS
    Technical Note No. 36, Chapter 5.
"""
r_mod_to_mj2000_iau2006(JD_TT::Number) = r_mod_to_mj2000_iau2006(DCM, JD_TT)

function r_mod_to_mj2000_iau2006(T::Type, JD_TT::Number)
    # Compute the angles used in the precession model.
    Ψ_a, ω_a, χ_a = precession_iau2006(JD_TT)
    ϵ_0 = 84381.406 * π / 648_000

    # NOTE: According to [2], the matrix as written in [1](p. 218) rotates the
    # MJ2000 to the MOD. Hence, we need the inverse matrix. Furthermore, the
    # equation in [1](p. 218, eq. 3-73) uses mϵ_2000 instead of ϵ_0 as in
    # [2](eq. 12).
    return compose_rotation(
        angle_to_rot(T, -χ_a, ω_a, Ψ_a, :ZXZ),
        angle_to_rot(T, -ϵ_0, 0, 0, :XYZ)
    )
end

"""
    r_mj2000_to_mod_iau2006([T::Type,] JD_TT::Number)

Compute the rotation that aligns the J2000 mean equatorial frame with the Mean
of Date (MOD) reference frame with the at Julian day `JD_TT` [Terrestrial Time].
This algorithm uses the IAU-2006 theory.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the MJ2000 frame with the MOD frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The J2000 reference frame here is not equal to the previous definition in FK5
theory. It is the reason why it is internally called `MJ2000`. According to
**[3]**:

> The mean equinox of J2000.0 to be considered is not the “rotational dynamical
> mean equinox of J2000.0” as used in the past, but the “inertial dynamical mean
> equinox of J2000.0” to which the recent numerical or analytical solutions
> refer.  The latter is associated with the ecliptic in the inertial sense,
> which is the plane perpendicular to the angular momentum vector of the orbital
> motion of the Earth-Moon barycenter as computed from the velocity of the
> barycenter relative to an inertial system. The rotational equinox is
> associated with the ecliptic in the rotational sense, which is perpendicular
> to the angular momentum vector computed from the velocity referred to the
> rotating orbital plane of the Earth-Moon barycenter. (The difference between
> the two angular momenta is the angular momentum associated with the rotation
> of the orbital plane.)

# References

- **[3]**: IERS (2010). Transformation between the International Terrestrial
    Reference System and the Geocentric Celestial Reference System. IERS
    Technical Note No. 36, Chapter 5.
"""
r_mj2000_to_mod_iau2006(JD_TT::Number) = r_mj2000_to_mod_iau2006(DCM, JD_TT)

function r_mj2000_to_mod_iau2006(T::Type, JD_TT::Number)
    return inv_rotation(r_mod_to_mj2000_iau2006(T, JD_TT))
end

################################################################################
#                               MJ2000 <=> GCRF
################################################################################

"""
    r_mj2000_to_gcrf_iau2006([T::Type,] JD_TT::Number = 0)

Compute the rotation that aligns the J2000 mean equatorial frame with the
Geocentric Celestial Reference Frame (GCRF). This algorithm uses the IAU-2006
theory. Notice that this rotation is just a bias matrix that does not depend on
the date. However, this function receives the argument `JD_TT` just to keep the
API compatibility.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

!!! info
    According to **[1]**, the frame bias that converts MJ2000 <=> GCRF is not a
    precise transformation for all the times.

# Returns

The rotation that aligns the MJ2000 frame with the MOD frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
r_mj2000_to_gcrf_iau2006(JD_TT::Number = 0) = r_mj2000_to_gcrf_iau2006(DCM, JD_TT)

function r_mj2000_to_gcrf_iau2006(T::Type, JD_TT::Number = 0)
    # Auxiliary variables.
    d2r = π / 180
    a2d = 1 / 3600
    a2r = a2d * d2r

    δα₀ = -0.0146 * a2r
    ξ₀  = -0.041775 * sin(84381.448 * a2r) * a2r
    η₀  = -0.0068192 * a2r

    return angle_to_rot(T, η₀, -ξ₀, -δα₀, :XYZ)
end

"""
    r_gcrf_to_mj2000_iau2006([T::Type,] JD_TT::Number = 0)

Compute the rotation that aligns the Geocentric Celestial Reference Frame (GCRF)
with the J2000 mean equatorial frame. This algorithm uses the IAU-2006 theory.
Notice that this rotation is just a bias matrix that does not depend on the
date. However, this function receives the argument `JD_TT` just to keep the API
compatibility.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

!!! info
    According to **[1]**, the frame bias that converts MJ2000 <=> GCRF is not a
    precise transformation for all the times.

# Returns

The rotation that aligns the MJ2000 frame with the MOD frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
r_gcrf_to_mj2000_iau2006(JD_TT::Number = 0) = r_gcrf_to_mj2000_iau2006(DCM, JD_TT)

function r_gcrf_to_mj2000_iau2006(T::Type, JD_TT::Number = 0)
    return inv_rotation(r_mj2000_to_gcrf_iau2006(T, JD_TT))
end

################################################################################
#                              Multiple Rotations
################################################################################

# The functions with multiple rotations must be added here only when the it will
# decrease the computational burden compared to calling the functions with the
# single rotations.

"""
    r_tirs_to_mod_iau2006([T::Type,] JD_UT1::Number, JD_TT::Number, δΔϵ_2000::Number = 0, δΔΨ_2000::Number = 0)

Compute the rotation that aligns the Terrestrial Intermediate Reference System
(TIRS) with the Mean of Date (MOD) reference frame at the Julian Day
`JD_UT1` [UT1] and `JD_TT` [Terrestrial Time]. This algorithm uses the IAU-2006
theory.

Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_2000`) and in longitude (`δΔψ_2000`) [rad] that are usually obtained
from IERS EOP Data (see [`get_iers_eop`](@ref) and [`deps_dpsi`](@ref)). This
corrections are related to Free Core Nutation (FCN) that models the effect of a
liquid Earth core.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

!!! info
    This composed rotation TIRS <=> ERS <=> MOD is implemented as a new function
    because the single rotations TIRS <=> ERS and ERS <=> MOD call the function
    `nutation_eo`, which has a high computational burden. In this case, the
    composed algorithm is about 2x faster than calling those function
    separately.

# Returns

The rotation that aligns the TIRS frame with the ERS frame. The rotation
representation is selected by the optional parameter `T`.
"""
function r_tirs_to_mod_iau2006(
    JD_UT1::Number,
    JD_TT::Number,
    δΔϵ_2000::Number = 0,
    δΔΨ_2000::Number = 0
)
    return r_tirs_to_mod_iau2006(DCM, JD_UT1, JD_TT, δΔϵ_2000, δΔΨ_2000)
end

function r_tirs_to_mod_iau2006(
    T::Type,
    JD_UT1::Number,
    JD_TT::Number,
    δΔϵ_2000::Number = 0,
    δΔΨ_2000::Number = 0
)
    # In this theory, the rotation of Earth is taken into account by the Earth
    # Rotation Angle, which is the angle between the Conventional International
    # Origin (CIO) and the Terrestrial Intermediate Origin (TIO) [1]. The latter
    # is a reference meridian on Earth that is located about 100m away from
    # Greenwich meridian along the equator of the Celestial Intermediate Pole
    # (CIP) [1].
    θ_era = 2π * (0.7790572732640 + 1.00273781191135448 * (JD_UT1 - JD_J2000))
    θ_era = mod(θ_era, 2π)

    # Compute the Equation of the Origins (EO).
    mϵ_2000, Δϵ_2000, ΔΨ_2000, EO = nutation_eo_iau2006(JD_TT, δΔϵ_2000, δΔΨ_2000)

    # Compute the Greenwich apparent sidereal angle (GAST).
    θ_gast2000 = θ_era - EO

    # Compute the rotation between the TIRS and ERS.
    r_ERS_TIRS = angle_to_rot(T, -θ_gast2000, 0, 0, :ZYX)

    # Compute the rotation between ERS and MOD.
    r_MOD_ERS = angle_to_rot(T, mϵ_2000 + Δϵ_2000, ΔΨ_2000, -mϵ_2000, :XZX)

    # Return the composed rotation.
    return compose_rotation(r_ERS_TIRS, r_MOD_ERS)
end

"""
    r_mod_to_tirs_iau2006([T::Type,] JD_UT1::Number, JD_TT::Number, δΔϵ_2000::Number = 0, δΔΨ_2000::Number = 0)

Compute the rotation that aligns the Mean of Date (MOD) reference frame with the
Terrestrial Intermediate Reference System (TIRS) at the Julian Day `JD_UT1`
[UT1] and `JD_TT` [Terrestrial Time]. This algorithm uses the IAU-2006
theory.

Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_2000`) and in longitude (`δΔψ_2000`) \\[rad] that are usually obtained
from IERS EOP Data (see [`get_iers_eop`](@ref) and [`deps_dpsi`](@ref)). This
corrections are related to Free Core Nutation (FCN) that models the effect of a
liquid Earth core.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

!!! info
    This composed rotation TIRS <=> ERS <=> MOD is implemented as a new function
    because the single rotations TIRS <=> ERS and ERS <=> MOD call the function
    `nutation_eo`, which has a high computational burden. In this case, the
    composed algorithm is about 2x faster than calling those function
    separately.

# Returns

The rotation that aligns the TIRS frame with the ERS frame. The rotation
representation is selected by the optional parameter `T`.
"""
function r_mod_to_tirs_iau2006(
    JD_UT1::Number,
    JD_TT::Number,
    δΔϵ_2000::Number = 0,
    δΔΨ_2000::Number = 0
)
    return r_mod_to_tirs_iau2006(DCM, JD_UT1, JD_TT, δΔϵ_2000, δΔΨ_2000)
end

function r_mod_to_tirs_iau2006(
    T::Type,
    JD_UT1::Number,
    JD_TT::Number,
    δΔϵ_2000::Number = 0,
    δΔΨ_2000::Number = 0
)
    return inv_rotation(r_tirs_to_mod_iau2006(T, JD_UT1, JD_TT, δΔϵ_2000, δΔΨ_2000))
end


text:
module TestCopulas

using Copulas
using Base.Test

##############
## asym_fgm ##
##############

uLow = 0.001
uHigh = 0.999
testGrid = [uLow, 0.1:0.1:0.9, uHigh]

epsTol = 1e-8

cop = Copulas.AsymFGM(0.4)
        Copulas.IndepPC_Cpp([0.3]),
        Copulas.IndepPC_Cpp([0.3]),
        ## Copulas.ClaytonPC_Cpp([4.3]),
        ## Copulas.ClaytonPC_Cpp([2]),

Copulas.testHandV(cop, epsTol, testGrid)
Copulas.testCdf(cop, epsTol, testGrid)
Copulas.testPdf(cop, testGrid)

end

ext: jl


text:
function adapt_stage_one(design, data, solver; p0 = nothing, p1 = nothing, lambda1 = 10, lambda2 = 10, grid_resolution = 0.05, mincpr = 0.5)

  params = parameters(design)

  n1         = length(data)
  data1      = data[1:n1]
  n1old      = interimsamplesize(design)
  n1star     = min(n1, n1old)
  tmp        = data1[1:n1star]
  x1star     = tmp[.!DataFrames.isna.(tmp)] |> sum # number of observed resoponses on joint range
  nna1sstar  = tmp[DataFrames.isna.(tmp)] |> length # number of NA on joint range
  x1         = data[.!DataFrames.isna.(data1)] |> sum # number of observed responses on entire range
  nna1       = data[DataFrames.isna.(data1[1:n1])] |> length
  na_inds    = find(DataFrames.isna.(data))
  nna        = length(na_inds)

  if (n1 == n1old) & (nna1 == 0) # nothing to do!
    return(design)
  end

  # construct integer programming model
  params.samplespace.specialcvalues =
    [ params.samplespace.specialcvalues;
      0:(maxsamplesize(params.samplespace) - 1)
    ]
  ipm = IPModel(params.samplespace, n1)
  m   = ipm.m
  y   = ipm.y

  # overwrite null and alternative if necessary and
  # construct grid for conditional error control
  p0        = p0 == nothing ? null(params) : p0
  p1        = p1 == nothing ? mcrv(params) : p1
  null_grid = collect(grid_resolution:grid_resolution:p0)
  alt_grid  = collect(p1:grid_resolution:(1 - grid_resolution))

  # create array of possible stage one outcomes (all possible imputations, 2^(# missing))
  stage_one = repmat(transpose(data), 2^nna)
  if nna > 0
    for i in 0:(2^nna - 1)
      stage_one[i + 1, na_inds] = [parse(Int, bin(i, nna)[j]) for j in 1:nna]
    end
  end
  n_imp = size(stage_one, 1)

  # enforce conditional properties for all possible stage-one outcomes
  JuMP.@variable(ipm.m, max_cpr_violation_1 >= 0)
  for i in 1:n_imp
    data_imp = stage_one[i, :] # imputed stage one data
    # stage two must be invariant under imputation
    if i > 1 # i == 1 corresponds to only 0s imputed -> minimal x1
      x1 = Integer(sum(data_imp))
      for n in ipm.nvals, c in ipm.cvals
        if !((c - 1) in ipm.cvals) # c - 1 not available, cannot produce locally constant second stage
          JuMP.@constraint(m,
            y[x1, n, c] == 0
          )
        else # c must be increasing with x1 to guarantee invariant stage two!
          JuMP.@constraint(m,
            y[x1, n, c] - y[x1 - 1, n, c - 1] == 0
          )
        end
      end
    end

    # unconditional type one error rate control, sufficient condition:
    for p in null_grid
      JuMP.@constraint(m,
        cpr_new(ipm, data_imp, p) <= cpr_old(design, data_imp, p)
      )
    end

    # make sure that trial continues if possible
    JuMP.@constraint(m,
      max_cpr_violation_1 >= cpr_old(design, data_imp, p1) - cpr_new(ipm, data_imp, p1) # max_cpr_violation_1 is larger than the biggest undershooting
    )
  end

  p_grid = collect(grid_resolution:grid_resolution:(1 - grid_resolution))

  # ensure monotonicity of cpr in x1
  for x1 in 1:n1, p in [p0; p1]
    JuMP.@constraint(m,
      sum(_cpr(x1, n1, n, c, p) * y[x1, n, c] - _cpr(x1 - 1, n1, n, c, p) * y[x1 - 1, n, c] for
        n in ipm.nvals, c in ipm.cvals
      ) >= 0
    )
  end

  # ensure new design maintains type one error as well
  JuMP.@constraint(ipm.m,
    pwr_new(ipm, n1, p0) <= power(design, p0)
  )

  # define auxiliary slack variable greater than maximal undershooting of minimal cpr
  mincpr_old = minimum(power.(design, 0:n1old, p1))
  JuMP.@variable(ipm.m, max_cpr_violation_2 >= 0)
  for x1 in 0:n1 # for any possible outcome, max_cpr_violation is larger than max(mincpr, mincpr_old) - cpr
    cvalsfinite   = ipm.cvals[isfinite.(ipm.cvals)]
    cvalsinfinite = ipm.cvals[.!isfinite.(ipm.cvals)]
    JuMP.@constraint(m,
      max_cpr_violation_2 >=
      max(mincpr, mincpr_old) -
      sum(_cpr(x1, n1, n, c, p1) * y[x1, n, c] for
        n  in ipm.nvals,
        c  in cvalsfinite
      ) -
      sum(y[x1, n, c] for # deactivate constraint for early stopping
        n  in ipm.nvals,
        c  in cvalsinfinite
      )
    )
  end

  # compute expected absolute power difference old/new
  # JuMP.@variable(ipm.m, abs_pwr_diff[p in p_grid] >= 0)
  # for p in p_grid
  #   JuMP.@constraint(ipm.m,
  #     pwr_new(ipm, n1, p) - power(design, p) <= abs_pwr_diff[p]
  #   )
  #   JuMP.@constraint(ipm.m,
  #     -(pwr_new(ipm, n1, p) - power(design, p)) <= abs_pwr_diff[p]
  #   )
  # end
  # p_RV = Distributions.Beta(1, 1)
  # JuMP.@expression(ipm.m, avg_abs_pwr_diff,
  #   sum(abs_pwr_diff[p] * Distributions.pdf(p_RV, p) for p in p_grid)
  # )

  JuMP.@variable(ipm.m, max_abs_pwr_diff >= 0)
  JuMP.@variable(ipm.m, max_abs_ess_diff >= 0)
  for p in p_grid
    JuMP.@constraint(ipm.m,
      pwr_new(ipm, n1, p) - power(design, p) <= max_abs_pwr_diff
    )
    JuMP.@constraint(ipm.m,
      -(pwr_new(ipm, n1, p) - power(design, p)) <= max_abs_pwr_diff
    )
    JuMP.@constraint(ipm.m,
      ess_new(ipm, n1, p) - mean(SampleSize(design, p)) <= max_abs_ess_diff
    )
    JuMP.@constraint(ipm.m,
      -(ess_new(ipm, n1, p) - mean(SampleSize(design, p))) <= max_abs_ess_diff
    )
  end

  # objective: expected quadratic difference of sample size under the two designs
  # JuMP.@objective(m, :Min, diff_sample_size(design, ipm, n1) + lambda1 * diff_pwr(design, ipm, n1) + lambda2 * max_cpr_violation_2 + lambda2 * max_cpr_violation_1)
  JuMP.@objective(m, :Min,
    max_abs_ess_diff + lambda1 * max_abs_pwr_diff + lambda2 * max_cpr_violation_1 + lambda2 * max_cpr_violation_2
  )

  JuMP.setsolver(ipm.m, solver)
  JuMP.solve(ipm.m) in (:Optimal, :UserLimit) ? nothing : error("no feasible solution reached")
  design = extractsolution(ipm, params)
  return design

end

function adapt_stage_two(design, data; p0 = nothing, p1 = nothing)

  params   = parameters(design)
  n1       = interimsamplesize(design)
  data1    = data[1:n1]
  x1       = Integer(data1[.!DataFrames.isna.(data1)] |> sum) # number of observed responses stage one
  nna1     = DataFrames.isna.(data1) |> sum
  na1_inds = find(DataFrames.isna.(data1))
  data2    = data[(n1 + 1):end]
  x2       = Integer(data2[.!DataFrames.isna.(data2)] |> sum) # number of observed responses stage one
  nna2     = DataFrames.isna.(data2) |> sum
  na2_inds = find(DataFrames.isna.(data2))
  nna      = nna2 + nna1
  n        = length(data)
  n2       = n - n1
  nold     = samplesize(design, x1)
  n2old    = nold - interimsamplesize(design) # invariant under imputation assumed
  cold     = criticalvalue(design, x1)
  c2old    = cold - x1 # c2old is invariant under missing in stage one

  p0        = p0 == nothing ? null(params) : p0
  p1        = p1 == nothing ? mcrv(params) : p1

  if n == nold # nothing to do!
    return(design)
  end

  new_design = deepcopy(design)

  # adjust to new sample size
  new_design.n[(x1:(x1 + nna1)) + 1] = n

  function conditional_probability_to_reject(n2, c2, data2, p)
    if n2 <= length(data2)
      x2star = sum(data2[1:n2])
      return x2star > c2 ? 1.0 : 0.0
    else
      x2star = sum(data2)
      if x2star > c2
        return 1.0
      else
        dn2 = n2 - length(data2) # remaining patients
        dx2 = c2 - x2star # need dx2 + 1 responses to have x2star + dx2 + 1 > c2
        return sum(Distributions.pdf.(Distributions.Binomial(dn2, p), (dx2 + 1):dn2))
      end
    end
  end

  # create array of possible stage one outcomes (all possible imputations, 2^(# missing))
  stage_two = repmat(transpose(data2), 2^nna2)
  if nna2 > 0
    for i in 0:(2^nna2 - 1)
      stage_two[i + 1, na2_inds] = [parse(Int, bin(i, nna2)[j]) for j in 1:nna2]
    end
  end
  n_imp = size(stage_two, 1)

  function possible_c(c)
    c2 = c - x1
    if c2 < 0
      return false
    end
    res = falses(n_imp)
    for i in 1:n_imp
      data2   = stage_two[i, :]
      cpr_old = conditional_probability_to_reject(n2old, c2old, data2, p0)
      cpr_new = conditional_probability_to_reject(n2, c2, data2, p0)
      # check if conditional probability of rejection is smaller
      # but also check conditional on stage one outcome only
      if !isfinite(c2old)
        println(c2old)
      end
      # (sum(dbinom.((c2 + 1):n2, n2, p0)) <= sum(dbinom.((c2old + 1):n2old, n2old, p0)))
      if (cpr_new <= cpr_old) & (sum(dbinom.((c2 + 1):n2, n2, p0)) <= power(design, x1, p0))
        res[i] = true
      else
        res[i] = false
      end
    end
    return all(res)
  end

  c = n
  while possible_c(c - 1)
    c -= 1
    if c == 0
      break
    end
  end
  if c == n
    c = Inf
  end

  new_design.c[(x1:(x1 + nna1)) + 1] = c + (0:nna1)

  return(new_design)

end

function diff_sample_size(design, ipm, n1; alpha = 1, beta = 1)
  # compute expected relative absolute deviation from sample size of original design
  # for old design n() and new n'()
  # default uses uniform beta binomial distributions, alpha and beta of
  # beta mixture can be adjusted though
  dn     = JuMP.AffExpr(0.0)
  n1_old = interimsamplesize(design)
  X1     = Distributions.BetaBinomial(n1_old, alpha, beta)
  X1_bar = Distributions.BetaBinomial(n1, alpha, beta)
  for x1 in 0:n1_old, x1_bar in 0:n1, n in ipm.nvals, c in ipm.cvals
    JuMP.push!(dn,
      abs(n - samplesize(design, x1)) / samplesize(design, x1) * Distributions.pdf(X1, x1) * Distributions.pdf(X1_bar, x1_bar),
      ipm.y[x1_bar, n, c]
    )
  end
  return dn
end

function diff_pwr(design, ipm, n1; alpha = 1, beta = 1)
  # compute expected relative absolute deviation from sample size of original design
  # for old design n() and new n'()
  # default uses uniform beta binomial distributions, alpha and beta of
  # beta mixture can be adjusted though
  dn     = JuMP.AffExpr(0.0)
  n1_old = interimsamplesize(design)
  pdf_X1(x1, p) = Distributions.pdf(Distributions.Binomial(n1_old, p), x1)
  pdf_X1_bar(x1_bar, p) = Distributions.pdf(Distributions.Binomial(n1, p), x1_bar)
  for x1 in 0:n1_old, x1_bar in 0:n1, n in ipm.nvals, c in ipm.cvals
    JuMP.push!(dn,
      QuadGK.quadgk(
        p -> abs(
          _cpr(x1, n1_old, samplesize(design, x1), criticalvalue(design, x1), p)
          - _cpr(x1_bar, n1, n, c, p)
        ) * pdf_X1(x1, p) * pdf_X1_bar(x1_bar, p) * Distributions.pdf(Distributions.Beta(alpha, beta), p),
        0, 1, abstol = 0.0001
      )[1],
      ipm.y[x1_bar, n, c]
    )
  end
  return dn
end

function cpr_new(ipm::IPModel, data, p)
  # conditional probability to reject of new design

  x1 = Integer(sum(data))
  n1 = length(data)
  cpr = JuMP.AffExpr(0.0)
  for n in ipm.nvals, c in ipm.cvals
    JuMP.push!(cpr,
      _cpr(x1, n1, n, c, p),
      ipm.y[x1, n, c]
    )
  end
  return cpr
end

function pwr_new(ipm::IPModel, n1, p)
  # probability to reject of new design

  X1  = Distributions.Binomial(n1, p)
  x1_range = collect(0:n1)
  pwr = JuMP.AffExpr(0.0)
  for x1 in x1_range, n in ipm.nvals, c in ipm.cvals
    JuMP.push!(pwr,
      Distributions.pdf(X1, x1) * _cpr(x1, n1, n, c, p),
      ipm.y[x1, n, c]
    )
  end
  return pwr
end

function ess_new(ipm::IPModel, n1, p)

  X1  = Distributions.Binomial(n1, p)
  x1_range = collect(0:n1)
  ess = JuMP.AffExpr(0.0)
  for x1 in x1_range, n in ipm.nvals, c in ipm.cvals
    JuMP.push!(ess,
      Distributions.pdf(X1, x1) * n,
      ipm.y[x1, n, c]
    )
  end

  return ess

end

function cpr_old(design::Design, data, p)
  # conditional probability to reject of old design

  n1old = interimsamplesize(design)
  n1new = length(data)

  if n1old == n1new
    x1 = Integer(sum(data))
    return power(design, x1, p)
  end
  if n1old <= n1new
    x1star  = Integer(sum(data[1:n1old])) # number of responses observed in
                                          # first stage under original design
    if samplesize(design, x1star) == n1old # trial stops early
      return power(design, x1star, p)
    else # trial continues to stage two
      if n1new >= samplesize(design, x1star) # outcome under original trial
                                              # known
        x = Integer(sum(data[1:samplesize(design, x1star)]))
        return x > criticalvalue(design, x1star) ? 1.0 : 0.0
      else # outcome still not necessarily decided
        n2star = samplesize(design, x1star) - n1new # number of samples to
                                                    # complete original stage two
        xstar   = Integer(sum(data)) # total responses observed so far
        c2star  = criticalvalue(design, x1star) - xstar # c2star + 1 responses
                                                        # needed in remaining
                                                        # n2star patients
        return sum(dbinom.((c2star + 1):n2star, n2star, p))
      end
    end
  end
  if n1old > n1new # premature interim analysis
    dn1 = n1old - n1new
    x1star = Integer(sum(data))
    return sum(dbinom.(0:dn1, dn1, p) .* power.(design, x1star + (0:dn1), p))
  end
end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
@testset "tags" begin
    @testset "exact fit tags" begin
        tf = TiffImages.TiffFile(UInt32)
        t1 = TiffImages.Tag(UInt16(TiffImages.IMAGEWIDTH), 0x00000200)
        @test write(tf, t1) # should fit so true

        seekstart(tf.io)
        t2 = read(tf, TiffImages.Tag)
        @test t1 == t2
    end

    @testset "small tags" begin
        tf = TiffImages.TiffFile(UInt32)
        t1 = TiffImages.Tag(UInt16(TiffImages.COMPRESSION), 0x0001)
        @test write(tf, t1) # should fit but needs padding
        @test position(tf.io) == 12 # should be full length

        seekstart(tf.io)
        t2 = read(tf, TiffImages.Tag)
        @test t2.data == 0x0001
    end

    @testset "large tags" begin
        tf = TiffImages.TiffFile(UInt32)
        offsets = UInt32[8, 129848, 259688, 389528]
        t1 = TiffImages.Tag(UInt16(TiffImages.STRIPOFFSETS), offsets)

        @test !write(tf, t1) # should fail to write into tf since it's too large
        @test position(tf.io) == 12

        data = Array{UInt8}(undef, 12)
        seekstart(tf.io)
        read!(tf.io, data)

        @test all(data .== 0x00)

        seekstart(tf.io)
        @test write(tf, t1, 12)

        seekstart(tf.io)
        t2 = read(tf, TiffImages.Tag)

        # make sure the data field contains our single offset from above
        @test typeof(t2.data) <: TiffImages.RemoteData
        @test t2.data.position == 12
    end

    @testset "String length equal to offset size" begin
        tf = TiffImages.TiffFile(UInt32)
        t1 = TiffImages.Tag(TiffImages.SOFTWARE, "test")

        # should fail because it's too large to fit
        @test !write(tf, t1)
        @test position(tf.io) == 12

        t3 = TiffImages.Tag(TiffImages.SOFTWARE, "tes\0")
        seekstart(tf.io)
        @test write(tf, t3)
        @test position(tf.io) == 12

        seekstart(tf.io)
        t4 = read(tf, TiffImages.Tag)

        # Since this is a NUL-terminated string of length 4 it should fit
        @test getfield(t4, :data) == "tes\0"
        @test t4.data == "tes"
    end
end

@testset "ifds" begin
    tf = TiffImages.TiffFile(UInt32)
    ifd = TiffImages.IFD(UInt32)
    ifd[TiffImages.IMAGEDESCRIPTION] = "Testing IFD read/write"
    ifd[TiffImages.IMAGEWIDTH] = UInt32(512)
    ifd[TiffImages.IMAGELENGTH] = UInt32(512)
    ifd[TiffImages.XRESOLUTION] = Rational{UInt32}(72, 1)
    ifd[TiffImages.COMPRESSION] = TiffImages.COMPRESSION_NONE
    ifd[TiffImages.SAMPLESPERPIXEL] = 1
    ifd[TiffImages.PHOTOMETRIC] = TiffImages.PHOTOMETRIC_MINISBLACK
    ifd[TiffImages.STRIPBYTECOUNTS] = UInt32[512*512]
    ifd[TiffImages.BITSPERSAMPLE] = 8
    ifd[TiffImages.ICCPROFILE] = Any[0x00, 0x00, 0x0b, 0xe8, 0x00, 0x00, 0x00, 0x00, 0x02, 0x00, 0x00, 0x00, 0x6d]

    write(tf, ifd)

    seekstart(tf.io)
    read_ifd, next_ifd = read(tf, TiffImages.IFD)
    TiffImages.load!(tf, read_ifd)

    @test all(ifd .== read_ifd)

    expected = TiffImages.IFDLayout(1, 512, 512, 262144,
                                    UInt8, UInt8, FixedPointNumbers.Normed{UInt8,8},
                                    TiffImages.COMPRESSION_NONE,
                                    TiffImages.PHOTOMETRIC_MINISBLACK)
    @test TiffImages.output(ifd) == expected

    delete!(ifd, TiffImages.COMPRESSION)

    @test TiffImages.output(ifd) == expected
end

@testset "Simple 2D image" begin
    filepath = get_example("house.tif")
    img = TiffImages.load(filepath)

    img2 = TiffImages.DenseTaggedImage(Gray{Float64}.(img.data));
    path, io = mktemp()
    write(io, img2)
    img3 = TiffImages.load(path)
    @test eltype(img3) == Gray{Float64}
end

@testset "3D image" begin
    filepath = get_example("mri.tif")
    img = TiffImages.load(filepath)

    img2 = TiffImages.DenseTaggedImage(RGB{Float64}.(img.data));
    path, io = mktemp()
    write(io, img2)
    img3 = TiffImages.load(path)
    @test eltype(img3) == RGB{Float64}
end

@testset "BigTIFF saving" begin
    filepath = get_example("house.tif")
    img = TiffImages.load(filepath)

    ifd = TiffImages._constructifd(img.data, UInt64);
    img2 = TiffImages.DenseTaggedImage(img.data, ifd)
    path, io = mktemp()
    write(io, img2)
    img3 = TiffImages.load(path)
    @test eltype(img3) == GrayA{N0f8}
end

@testset "save function" begin
    data = rand(RGB{N0f8}, 128, 128)
    filepath = tempname()

    TiffImages.save(filepath, data)
    @test all(data .== TiffImages.load(filepath))
end


text:
using Distributed

@time @everywhere include("../BioClim/src/required.jl")

## Get & prepare data
@time @everywhere begin
    # Load data from CSV files
    df = CSV.read("/Users/kiristern/Documents/GitHub/SDM/data/pred_prey/Puma_concolor.csv", header=true)
    # Prepare data (select columns, arrange values)
    df = prepare_csvdata(df)
    # Separate species
    taxa_occ = [df[df.species .== u,:] for u in unique(df.species)]

    # Define coordinates range
    lon_range = (-165.0, -30.0)
    lat_range = (75.0, -60.0)
end

cd("/Users/kiristern/Documents/GitHub/SDM/")
## Get the worldclim data
@time wc_vars = pmap(x -> worldclim(x)[lon_range, lat_range], 1:19);

# Make predictions for all species
@time predictions = pmap(x -> species_bclim(x, wc_vars), taxa_occ);

begin
    # Create Y -> site-by-species community data table
    Y = zeros(Int64, (prod(size(predictions[1])),length(taxa_occ)))
    # Fill Y with community predictions
    @progress for gc in eachindex(predictions[1].grid) # loop for all sites
        # Group predictions for all species in site
        R = map(x -> x.grid[gc], predictions)
        # Fill Y with binary values -> 1 if species prediction for site != NaN, 0 if == NaN
        global Y[gc,:] = .!isnan.(R)
    end

    ## Arrange prediction values as grid
    # Create empty grid
    predict_grid = zeros(Float64, size(predictions[1]))

    # Create SDMLayer with prediction values
    predict_SDM = SDMLayer(predict_grid, predictions[1].left, predictions[1].right, predictions[1].bottom, predictions[1].top)
end

## Plot results
sdm_plot = plotSDM(predict_SDM)

## custom heatmap colour gradient
heatmap!(
    sdm_plot,
    longitudes(prediction), latitudes(prediction), prediction.grid,
    aspectratio = 92.60/60.75, c = cgrad([:lightblue, :blue, :darkblue]),
    clim = (0.0, maximum(filter(!isnan, prediction.grid)))
    )

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
module ShaleDrillingLikelihood_Drilling_Test

const DOBTIME = false

using ShaleDrillingLikelihood

using Test
using StatsFuns
using Random
using BenchmarkTools
using Base.Threads
using InteractiveUtils
using Distributions

using Calculus
using Optim
using LinearAlgebra

using LinearAlgebra: checksquare
using Optim: minimizer

using ShaleDrillingLikelihood: SimulationDraws,
    DataDrill,
    TestDrillModel,
    DrillingTmpVars,
    DrillUnit,
    AbstractDrillRegime,
    DrillLease,
    ObservationDrill,
    InitialDrilling,
    DevelopmentDrilling,
    theta_drill_ρ,
    _model,
    update!,
    loglik_drill_lease!,
    loglik_drill_unit!,
    simloglik_drill_unit!,
    simloglik_drill_data!,
    _y

println("testing drilling likelihood")

@testset "Drilling Likelihood" begin

    Random.seed!(911)

    theta = [-0.3, 2.0, -2.0, -0.75, 0.5]
    @test 0.5 == theta_drill_ρ(TestDrillModel(), theta)

    data = DataDrill(
        TestDrillModel(), theta;
        minmaxleases=1:1,
        num_i=100, nper_initial=10:40, nper_development=10:40,
        num_zt=200, tstart=1:50
    )

    sim = SimulationDraws(100, data)
    println("number of periods is $(length(_y(data)))")

    grad = zeros(length(theta))
    hess = zeros(length(theta),length(theta))
    dtv = DrillingTmpVars(data)

    # fill!(grad, 0)
    # update!(sim, theta_drill_ρ(_model(data), theta))
    # @code_warntype simloglik_drill_unit!(grad, data[1], theta, view(sim, 1), true)
    # simloglik_drill_unit!(grad, data[3], theta, view(sim, 3), false)
    LL1 = simloglik_drill_data!(grad, hess, data, theta, sim, true)
    LL2 = simloglik_drill_data!(grad, hess, data, theta, sim, false)
    @test isfinite(LL1)
    @test isfinite(LL2)
    @test LL1 ≈ LL2

    # check that we don't update gradient here
    simloglik_drill_data!(grad, hess, data, theta, sim, true)
    gradcopy = copy(grad)
    simloglik_drill_data!(grad, hess, data, theta.*2, sim, false)
    @test all(grad .== gradcopy)

    # @testset "warntypes for drilling" begin
    #     unit = data[1]
    #     lease = unit[InitialDrilling()][1]
    #     simi = view(sim, 1)
    #     @code_warntype loglik_drill_lease!(  grad, lease, theta, simi[1], dtv[threadid()], true)
    #     @code_warntype loglik_drill_unit!(   grad, data[1], theta, simi[1], dtv[threadid()], true)
    #     @code_warntype simloglik_drill_unit!(grad, unit, theta, simi, true)
    # end

    # check finite difference
    @testset "check finite difference" begin
        let theta = theta.*0.5
            fdgrad = Calculus.gradient(x -> simloglik_drill_data!(grad, hess, data, x, sim, false), theta)
            fill!(grad, 0)
            simloglik_drill_data!(grad, hess, data, theta, sim, true)
            @test grad ≈ fdgrad
            @test !(grad ≈ zeros(length(grad)))
        end
    end

    let k = length(theta), hess = zeros(k,k)
        # @show @benchmark simloglik_drill_data!($grad, $hess, $data, $theta, $sim, false, false)
        # @show @benchmark simloglik_drill_data!($grad, $hess, $data, $theta, $sim, true, false)
        if DOBTIME
            print("")
            @show @benchmark simloglik_drill_data!($grad, $hess, $data, $theta, $sim, true, true)
            print("")
        end
    end

    println("initial tests done")

    @testset "drilling likelihood optimization" begin

        k = length(theta)

        function f(x)
            update!(sim, theta_drill_ρ(_model(data), x))
            LL = simloglik_drill_data!(zeros(k), zeros(k,k), data, x, sim, false, false)
            return -LL
        end

        function fg!(grad, x)
            tmphess = zeros(k,k)
            fill!(grad, 0)
            update!(sim, theta_drill_ρ(_model(data), x))
            LL = simloglik_drill_data!(grad, tmphess, data, x, sim, true, false)
            grad .*= -1
            return -LL
        end

        function h!(hess, x)
            grad = zeros(k)
            checksquare(hess) == length(x) || throw(DimensionMismatch())
            update!(sim, theta_drill_ρ(_model(data), x))
            LL = simloglik_drill_data!(grad, hess, data, x, sim, true, true)
            grad .*= -1
            return -LL
        end

        function invH0(x::AbstractVector)
            grad = zeros(k)
            hess = zeros(k,k)
            simloglik_drill_data!(grad, hess, data, x, sim, true, true)
            return inv(hess)
        end

        println("getting ready to optmize")
        odfg  = OnceDifferentiable(f, fg!, fg!, theta)
        tdfgh = TwiceDifferentiable(f, fg!, fg!, h!, theta)
        res = optimize(tdfgh, theta*0.5, BFGS(;initial_invH = invH0), Optim.Options(allow_f_increases=true, show_trace=true))
        @show res
        @show minimizer(res), theta

        vcovinv = invH0(res.minimizer)
        se = sqrt.(diag(vcovinv))
        tstats = minimizer(res) ./ se
        pvals = cdf.(Normal(), -2*abs.(tstats))
        coef_and_se = hcat(theta, minimizer(res), se, tstats, pvals)

        err = theta .- res.minimizer
        waldtest = err'*vcovinv*err
        @test ccdf(Chisq(length(theta)), waldtest) > 0.05
        @show coef_and_se

        rhohat = last(minimizer(res))
        rho_se = last(se)
        rho_t = (rhohat - last(theta)) / rho_se
        rho_p = ccdf(Normal(), -2*abs(rho_t))
        @show rhohat, last(theta), rho_t, rho_p

        # Base.showarray(stdout, coef_and_se)
    end

end # testset
end # module


text:
using Random

const Series = AbstractVector{Int}
const Ensemble = AbstractArray{Int,2}

mutable struct TransferDist
    b::Int
    k::Int
    ncond::Int
    counts::Int
    states::Vector{Int}
    histories::Vector{Int}
    sources::Vector{Int}
    predicates::Vector{Int}
    function TransferDist(b, k, n)
        if b < 2
            throw(DomainError(b, "base must be at least 2"))
        end
        if k < 1
            throw(DomainError(k, "history length must be at least 1"))
        end
        if n < 0
            throw(DomainError(n, "number of conditions must be non-negative"))
        end
        q = b^k
        r = b^n
        states = zeros(Int, b*b*q*r)
        histories = zeros(Int, q*r)
        sources = zeros(Int, b*q*r)
        predicates = zeros(Int, b*q*r)
        new(b, k, n, 0, states, histories, sources, predicates)
    end
end

isvalid(t::TransferDist) = t.counts != 0

function observe!(t::TransferDist, src::Series, dst::Series, cond::Union{Nothing,Ensemble} = nothing)
    if size(src) != size(dst)
        throw(ArgumentError("length(src) must equal length(dst)"))
    end
    if length(dst) <= t.k
        throw(ArgumentError("length(dst) must be greater than the history length"))
    end
    if isnothing(cond) && t.ncond != 0
        throw(ArgumentError("expected $(t.ncond) conditions, got 0"))
    end
    if !isnothing(cond) && size(cond, 1) != t.ncond
        throw(ArgumentError("expected $(t.ncond) conditions, got $(size(cond, 1))"))
    end
    if !isnothing(cond) && size(cond, 2) != length(dst)
        throw(ArgumentError("size(cond, 2) must be equal to length(dist)"))
    end

    t.counts += length(dst) - t.k
    src_state, future, state, source, predicate = 0, 0, 0, 0, 0
    history, q = 0, 1
    for i in 1:t.k
        q *= t.b
        history *= t.b
        history += dst[i]
    end
    for i in t.k+1:length(src)
        back = 0
        for c in 1:t.ncond
            back = back * t.b + cond[c, i-1]
        end
        history += back * q
        src_state = src[i-1]
        future = dst[i]
        source = history * t.b + src_state
        predicate = history * t.b + future
        state = predicate * t.b + src_state

        t.states[state + 1] += 1
        t.histories[history + 1] += 1
        t.sources[source + 1] += 1
        t.predicates[predicate + 1] += 1

        history = predicate - q * (dst[i - t.k] + back * t.b)
    end

    t
end

observe!(t::TransferDist, src::Series, dst::Series, cond::Series) = observe!(t, src, dst, cond')

function (t::TransferDist)()
    if !isvalid(t)
        error("TransferDist has no observed data")
    end
    te = 0.0
    for (h, history) in enumerate(t.histories)
        if history == 0
            continue
        end
        for f in 1:t.b
            p = (h - 1) * t.b + f
            predicate = t.predicates[p]
            if predicate == 0
                continue
            end
            for src in 1:t.b
                ss = (h - 1) * t.b + src
                source = t.sources[ss]
                if source == 0
                    continue
                end
                s = (p - 1) * t.b + src
                state = t.states[s]
                if state == 0
                    continue
                end
                te += state * log2((state * history) / (source * predicate))
            end
        end
    end
    te / t.counts
end

function transferentropy(src::Series, dst::Series, cond::Ensemble, k::Int)
    xs = src .- minimum(src)
    ys = dst .- minimum(dst)
    cs = cond .- minimum(cond; dims=2)
    b = max(1, maximum(xs), maximum(ys), maximum(cs)) + 1
    tedist = TransferDist(b, k, size(cond,1))
    observe!(tedist, xs, ys, cond)
    tedist()
end

function transferentropy(src::Series, dst::Series, k::Int)
    xs = src .- minimum(src)
    ys = dst .- minimum(dst)
    b = max(1, maximum(xs), maximum(ys)) + 1
    tedist = TransferDist(b, k, 0)
    observe!(tedist, xs, ys)
    tedist()
end

transferentropy(src::Series, dst::Series, cond::Series, k::Int) = transferentropy(src, dst, cond', k)

struct Significance
    value::Float64
    p::Float64
    se::Float64
end

function transferentropy(rng::AbstractRNG, src::Series, dst::Series, cond::Union{Series,Ensemble}, k::Int, nperm::Int)
    if nperm < 10
        throw(DomainError(nperm, "expected at least 10 permutation, got $(nperm)"))
    end

    te = transferentropy(src, dst, cond, k)
    count = 1
    for _ in 1:nperm
        permutedsrc = src[randperm(rng, length(src))]
        count += transferentropy(permutedsrc, dst, cond, k) >= te
    end
    p = count / (nperm + 1)
    se = sqrt((p * (1 - p)) / (nperm + 1))

    Significance(te, p, se)
end
transferentropy(src, dst, cond, k, nperm) = transferentropy(Random.GLOBAL_RNG, src, dst, cond, k, nperm)

function transferentropy(rng::AbstractRNG, src::Series, dst::Series, k::Int, nperm::Int)
    if nperm < 10
        throw(DomainError(nperm, "expected at least 10 permutation, got $(nperm)"))
    end

    te = transferentropy(src, dst, k)
    count = 1
    for _ in 1:nperm
        permutedsrc = src[randperm(rng, length(src))]
        count += transferentropy(permutedsrc, dst, k) >= te
    end
    p = count / (nperm + 1)
    se = sqrt((p * (1 - p)) / (nperm + 1))

    Significance(te, p, se)
end
transferentropy(src, dst, k, nperm) = transferentropy(Random.GLOBAL_RNG, src, dst, k, nperm)

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
const tJSite = TagType"tJ"

function siteinds(::tJSite,
                  N::Int; kwargs...)
  conserve_qns = get(kwargs,:conserve_qns,false)
  conserve_sz = get(kwargs,:conserve_sz,conserve_qns)
  conserve_nf = get(kwargs,:conserve_nf,conserve_qns)
  conserve_parity = get(kwargs,:conserve_parity,conserve_qns)
  if conserve_sz && conserve_nf
    em = QN(("Nf",0,-1),("Sz", 0)) => 1
    up = QN(("Nf",1,-1),("Sz",+1)) => 1
    dn = QN(("Nf",1,-1),("Sz",-1)) => 1
    return [Index(em,up,dn;tags="Site,tJ,n=$n") for n=1:N]
  elseif conserve_nf
    zer = QN("Nf",0,-1) => 1
    one = QN("Nf",1,-1) => 2
    return [Index(zer,one;tags="Site,tJ,n=$n") for n=1:N]
  elseif conserve_sz
    em = QN(("Sz", 0),("Pf",0,-2)) => 1
    up = QN(("Sz",+1),("Pf",1,-2)) => 1
    dn = QN(("Sz",-1),("Pf",1,-2)) => 1
    return [Index(em,up,dn;tags="Site,tJ,n=$n") for n=1:N]
  elseif conserve_parity
    zer = QN("Pf",0,-2) => 1
    one = QN("Pf",1,-2) => 2
    return [Index(zer,one;tags="Site,tJ,n=$n") for n=1:N]
  end
  return [Index(3,"Site,tJ,n=$n") for n=1:N]
end

function state(::tJSite,
               st::AbstractString)
  if st == "0" || st == "Emp"
    return 1
  elseif st == "Up" || st == "↑"
    return 2
  elseif st == "Dn" || st == "↓"
    return 3
  end
  throw(ArgumentError("State string \"$st\" not recognized for tJ site"))
  return 0
end

function op(::tJSite,
            s::Index,
            opname::AbstractString)::ITensor
  Emp = s(1)
  EmpP = s'(1)
  Up = s(2)
  UpP = s'(2)
  Dn = s(3)
  DnP = s'(3)

  Op = emptyITensor(s',dag(s))
  if opname == "Nup"
    Op[UpP, Up] = 1.
  elseif opname == "Ndn"
    Op[DnP, Dn] = 1.
  elseif opname == "Ntot"
    Op[UpP, Up] = 1.
    Op[DnP, Dn] = 1.
  elseif opname == "Cup" || opname == "Aup"
    Op[EmpP, Up] = 1.
  elseif opname == "Cdagup" || opname == "Adagup"
    Op[UpP, Emp] = 1.
  elseif opname == "Cdn" || opname == "Adn"
    Op[EmpP, Dn] = 1.
  elseif opname == "Cdagdn" || opname == "Adagdn"
    Op[DnP, Emp] = 1.
  elseif opname == "FermiPhase" || opname == "FP"
    Op[UpP, Up] = -1.
    Op[EmpP, Emp] = 1.
    Op[DnP, Dn] = -1.
  elseif opname == "Fup"
    Op[UpP, Up] = -1.
    Op[EmpP, Emp] = 1.
    Op[DnP, Dn] = 1.
  elseif opname == "Fdn"
    Op[UpP, Up] = 1.
    Op[EmpP, Emp] = 1.
    Op[DnP, Dn] = -1.
  elseif opname == "Sᶻ" || opname == "Sz"
    Op[UpP, Up] = 0.5
    Op[DnP, Dn] = -0.5
  elseif opname == "Sˣ" || opname == "Sx"
    Op[UpP, Dn] = 1.0
    Op[DnP, Up] = 1.0
  elseif opname == "S⁺" || opname == "Splus"
    Op[UpP, Dn] = 1.
  elseif opname == "S⁻" || opname == "Sminus"
    Op[DnP, Up] = 1.
  elseif opname == "Emp" || opname == "0"
    pEmp = emptyITensor(s)
    pEmp[Emp] = 1.
    return pEmp
  elseif opname == "Up" || opname == "↑"
    pU = emptyITensor(s)
    pU[Up] = 1.
    return pU
  elseif opname == "Dn" || opname == "↓"
    pD = emptyITensor(s)
    pD[Dn] = 1.
    return pD
  else
    throw(ArgumentError("Operator name '$opname' not recognized for tJSite"))
  end
  return Op
end

function has_fermion_string(::tJSite,
            s::Index,
            opname::AbstractString)::Bool
  if opname=="Cup" || opname=="Cdagup" || opname=="Cdn" || opname=="Cdagdn"
    return true
  end
  return false
end


text:
using DataDeps

register(DataDep(
  "yeast",
  "http://archive.ics.uci.edu/ml/datasets/Yeast",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/yeast/yeast.data",
  "7cf61776fc04f527f93bf57a327b863893a1225d82df02d457e8950173218258",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=10,
        feature_indices=2:9,
        id_indices=1,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# Register functions and handle literals
"""
$(SIGNATURES)

Registers a function call as a primitive for the `Operation` graph of the
ModelingToolkit IR. Example:

```julia
@register f(x,y)
```

registers `f` as a possible two-argument function.
"""
macro register(sig)
    splitsig = splitdef(:($sig = nothing))
    name = splitsig[:name]
    args = splitsig[:args]
    typargs = typed_args(args)
    defs = :()
    for typarg in typargs
        splitsig[:args] = typarg
        splitsig[:body] = :(Operation($name, Expression[$(args...)]))
        defs = :($defs; $(combinedef(splitsig)))
    end
    esc(defs)
end
# Create all valid combinations of Expression,Number for function signature
function typed_args(args)
    cases = Iterators.product(Tuple((:Expression, :Number) for i = 1:length(args))...)
    typargs = Vector{Any}[]
    for case in cases
        typarg = [Expr(:(::), args[i], case[i]) for i = 1:length(args)]
        push!(typargs, typarg)
    end
    typargs[1:end-1] # Last case is all Numbers but that should be specialized
end

# Binary & unary operators and functions
import DiffRules, SpecialFunctions, NaNMath
for (M, f, arity) in DiffRules.diffrules()
    fun = :($M.$f)
    sig = arity == 1 ? :($fun(x)) :
          arity == 2 && :($fun(x,y))
    @eval @register $sig
end

for fun ∈ [:!]
    basefun = Expr(:., Base, QuoteNode(fun))
    sig = :($basefun(x))
    @eval @register $sig
end

for fun ∈ [:<, :>, :(==), :&, :|, :div]
    basefun = Expr(:., Base, QuoteNode(fun))
    sig = :($basefun(x,y))
    @eval @register $sig
end

# ifelse
#@register Base.ifelse(cond,t,f)

# special cases
Base.:^(x::Expression,y::T) where T <: Integer = Operation(Base.:^, Expression[x, y])
Base.:^(x::Expression,y::T) where T <: Rational = Operation(Base.:^, Expression[x, y])

@register Base.conj(x)
@register Base.getindex(x,i)
Base.getindex(x::Operation,i::Int64) = Operation(getindex,[x,i])
Base.one(::Operation) = 1


text:
# this file is only used for testing - do not remove
Pair{String,Any}[]


text:
module VisualizationTests

using Test

@testset "visualization_data" begin include("VisualizationDataTests.jl") end

@testset "PrintOpTrees" begin include("PrintOpTreesTests.jl") end

@testset "Vtk" begin include("VtkTests.jl") end

end # module


text:
s = """
73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450
"""
s2 = replace(s, "\n" => "")
a = split(s2, "")
a2 = map(x -> parse(Int, x), a)
len = length(a2)

window = 13
maxsofar = 0
for i = 1:(len-window+1)
  global maxsofar
  p = prod(a2[i:i+window-1])
  maxsofar = max(maxsofar, p)
end
println(maxsofar)

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
y = [1, 2, 3, 4, 5, 6]
plot(y, mcol="red")
y = 0


text:
using BenchmarkTools
include("forest_fire.jl")

agent_properties = [:model]

height=100
d=0.6
nsteps = 100
steps_to_collect_data = collect(1:nsteps);

function counter(model::AbstractModel)
  on_fire = 0
  green = 0
  burned = 0
  for tree in model.agents
    if tree.status == 1
      green += 1
    elseif tree.status == 2
      on_fire += 1
    else
      burned += 1
    end
  end
  return green, on_fire, burned
end
aggregators = [counter]
size_range = 100:100:1000

results = Float64[]
for width in size_range
  b = @benchmarkable data=step!(dummystep, forest_step!, forest, nsteps, agent_properties, aggregators, steps_to_collect_data) setup=(forest=model_initiation(d=d, griddims=($width, height), seed=2))

  j = run(b)
  push!(results, minimum(j.times)/1e9)  # convert them to seconds
end

results

# 0.124292399
# 0.261890796
# 0.382272875
# 0.526156528
# 0.643057459
# 0.795138978
# 1.012378349
# 1.180964479
# 1.248486157
# 1.369958317


text:
"""
Configuration structure for Synchrony, it defines the Synchrony Web API endpoint and port.
"""
mutable struct GraffConfig
    apiEndpoint::String
    accessKey::String
    secretKey::String
    region::String
    userId::String
    robotId::String
    sessionId::String
    GraffConfig(apiEndpoint::String, accessKey::String, secretKey::String, region::String, userId::String, robotId::String, sessionId::String) = new(apiEndpoint, accessKey, secretKey, region, userId, robotId, sessionId)
    GraffConfig(apiEndpoint::String, accessKey::String, secretKey::String, region::String, userId::String) = new(apiEndpoint, accessKey, secretKey, region, userId, "", "")
end

function show(io::IO, c::GraffConfig)
    println(io, "GraffSDK configuration:")
    println(io, " - API endpoint: $(c.apiEndpoint)")
    println(io, " - API Access key: $(c.accessKey)")
    println(io, " - Region: $(c.region)")
    println(io, " - User ID: $(c.userId)")
    println(io, " - Robot ID: $(c.robotId)")
    println(io, " - Session ID: $(c.sessionId)")
end

"""
Standardized error response for any request. This is returned if any requests fail server-side.
"""
struct ErrorResponse
    message::String
    returnCode::Int
end

# The current Graff config.
global __graffConfig = nothing


text:
using Test

@testset "ObjectDetectors" begin
    include("static_object_detector_tests.jl")
    include("bounding_box_detection_tests.jl")
    include("oriented_bounding_box_detection_tests.jl")
    include("polygon_detection_tests.jl")
    include("iou_tests.jl")
end


text:
using SymbolicUtils: FnType, Sym
using Setfield

const IndexMap = Dict{Char,Char}(
            '0' => '₀',
            '1' => '₁',
            '2' => '₂',
            '3' => '₃',
            '4' => '₄',
            '5' => '₅',
            '6' => '₆',
            '7' => '₇',
            '8' => '₈',
            '9' => '₉')

struct VariableDefaultValue end

"""
$(TYPEDEF)

A named variable which represents a numerical value. The variable is uniquely
identified by its `name`, and all variables with the same `name` are treated
as equal.

# Fields
$(FIELDS)

For example, the following code defines an independent variable `t`, a parameter
`α`, a function parameter `σ`, a variable `x`, which depends on `t`, a variable
`y` with no dependents, a variable `z`, which depends on `t`, `α`, and `x(t)`
and parameters `β₁` and `β₂`.

```julia
σ = Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(:σ)) # left uncalled, since it is used as a function
w = Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(:w)) # unknown, left uncalled
x = Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(:x))(t)  # unknown, depends on `t`
y = Num(Variable(:y))   # unknown, no dependents
z = Num(Variable{Symbolics.FnType{NTuple{3,Any},Real}}(:z))(t, α, x)  # unknown, multiple arguments
β₁ = Num(Variable(:β, 1)) # with index 1
β₂ = Num(Variable(:β, 2)) # with index 2

expr = β₁ * x + y^α + σ(3) * (z - t) - β₂ * w(t - 1)
```
"""
struct Variable{T} <: Function # backward compat
    """The variable's unique name."""
    name::Symbol
    Variable(name) = Sym{Real}(name)
    Variable{T}(name) where T = Sym{T}(name)
    function Variable{T}(name, indices...) where T
        var_name = Symbol("$(name)$(join(map_subscripts.(indices), "ˏ"))")
        Sym{T}(var_name)
    end
end

function Variable(name, indices...)
    var_name = Symbol("$(name)$(join(map_subscripts.(indices), "ˏ"))")
    Variable(var_name)
end

# TODO: move this to Symutils
function Sym{T}(name, i, indices...) where T
    var_name = Symbol("$(name)$(join(map_subscripts.((i, indices...,)), "ˏ"))")
    Sym{T}(var_name)
end

function map_subscripts(indices)
    str = string(indices)
    join(IndexMap[c] for c in str)
end

rename(x::Sym,name) = @set! x.name = name
function rename(x::Symbolic, name)
    if operation(x) isa Sym
        @assert x isa Term
        @set! x.f = rename(operation(x), name)
        @set! x.hash = Ref{UInt}(0)
        return x
    else
        error("can't rename $x to $name")
    end
end

function unwrap_runtime_var(v)
    isruntime = Meta.isexpr(v, :$) && length(v.args) == 1
    isruntime && (v = v.args[1])
    return isruntime, v
end

# Build variables more easily
function _parse_vars(macroname, type, x, transform=identity)
    ex = Expr(:block)
    var_names = Symbol[]
    # if parsing things in the form of
    # begin
    #     x
    #     y
    #     z
    # end
    x = x isa Tuple && first(x) isa Expr && first(x).head == :tuple ? first(x).args : x # tuple handling
    x = flatten_expr!(x)
    cursor = 0
    isoption(ex) = Meta.isexpr(ex, [:vect, :vcat, :hcat])
    while cursor < length(x)
        cursor += 1
        v = x[cursor]

        # We need lookahead to the next `v` to parse
        # `@variables x [connect=Flow,unit=u]`
        nv = cursor < length(x) ? x[cursor+1] : nothing
        val = unit = connect = options = nothing

        # x = 1, [connect = flow; unit = u"m^3/s"]
        if Meta.isexpr(v, :(=))
            v, val = v.args
            if Meta.isexpr(val, :tuple) && length(val.args) == 2 && isoption(val.args[2])
                options = val.args[2].args
                val = val.args[1]
            end
        end

        # x [connect = flow; unit = u"m^3/s"]
        if isoption(nv)
            options = nv.args
            cursor += 1
        end

        isruntime, v = unwrap_runtime_var(v)
        iscall = Meta.isexpr(v, :call)
        isarray = Meta.isexpr(v, :ref)
        issym  = v isa Symbol
        @assert iscall || isarray || issym "@$macroname expects a tuple of expressions or an expression of a tuple (`@$macroname x y z(t) v[1:3] w[1:2,1:4]` or `@$macroname x y z(t) v[1:3] w[1:2,1:4] k=1.0`)"

        if iscall
            isruntime, fname = unwrap_runtime_var(v.args[1])
            call_args = map(last∘unwrap_runtime_var, @view v.args[2:end])
            var_name, expr = construct_vars(fname, type, call_args, val, options, transform, isruntime)
        else
            var_name, expr = construct_vars(v, type, nothing, val, options, transform, isruntime)
        end

        push!(var_names, var_name)
        push!(ex.args, expr)
    end
    rhs = build_expr(:vect, var_names)
    push!(ex.args, rhs)
    return ex
end

function construct_vars(v, type, call_args, val, prop, transform, isruntime)
    issym  = v isa Symbol
    isarray = isa(v, Expr) && v.head == :ref
    if isarray
        var_name = v.args[1]
        isruntime, var_name = unwrap_runtime_var(var_name)
        indices = v.args[2:end]
        expr = _construct_array_vars(isruntime ? var_name : Meta.quot(var_name), type, call_args, val, prop, indices...)
    else
        var_name = v
        expr = construct_var(isruntime ? var_name : Meta.quot(var_name), type, call_args, val, prop)
    end
    lhs = isruntime ? gensym(var_name) : var_name
    rhs = :($transform($expr))
    lhs, :($lhs = $rhs)
end

function option_to_metadata_type(::Val{opt}) where {opt}
    throw(Base.Meta.ParseError("unknown property type $opt"))
end

function setprops_expr(expr, props)
    isnothing(props) && return expr
    for opt in props
        if !Meta.isexpr(opt, :(=))
            throw(Base.Meta.ParseError(
                "Variable properties must be in " *
                "the form of `a = b`. Got $opt."))
        end

        lhs, rhs = opt.args

        @assert lhs isa Symbol "the lhs of an option must be a symbol"
        expr = :($setmetadata($expr,
                              $(option_to_metadata_type(Val{lhs}())),
                       $rhs))
    end
    expr
end

function construct_var(var_name, type, call_args, val, prop)
    expr = if call_args === nothing
        :($Num($Sym{$type}($var_name)))
    elseif !isempty(call_args) && call_args[end] == :..
        :($Num($Sym{$FnType{Tuple, $type}}($var_name))) # XXX: using Num as output
    else
        :($Num($Sym{$FnType{NTuple{$(length(call_args)), Any}, $type}}($var_name)($(map(x->:($value($x)), call_args)...))))
    end

    if val !== nothing
        expr = :($setmetadata($expr, $VariableDefaultValue, $val))
    end

    return setprops_expr(expr, prop)
end

function construct_var(var_name, type, call_args, val, prop, ind)
    # TODO: just use Sym here
    expr = if call_args === nothing
        :($Num($Sym{$type}($var_name, $ind...)))
    elseif !isempty(call_args) && call_args[end] == :..
        :($Num($Sym{$FnType{Tuple{Any}, $type}}($var_name, $ind...))) # XXX: using Num as output
    else
        :($Num($Sym{$FnType{NTuple{$(length(call_args)), Any}, $type}}($var_name, $ind...)($(map(x->:($value($x)), call_args)...))))
    end
    if val !== nothing
        expr = :($setmetadata($expr, $VariableDefaultValue, $val isa AbstractArray ? $val[$ind...] : $val))
    end

    return setprops_expr(expr, prop)
end

function _construct_array_vars(var_name, type, call_args, val, prop, indices...)
    :(map(Iterators.product($(indices...))) do ind
        $(construct_var(var_name, type, call_args, val, prop, :ind))
    end)
end

"""

Define one or more unknown variables.

```julia
@variables t α σ(..) β[1:2]
@variables w(..) x(t) y z(t, α, x)

expr = β[1]* x + y^α + σ(3) * (z - t) - β[2] * w(t - 1)
```

`(..)` signifies that the value should be left uncalled.

Sometimes it is convenient to define arrays of variables to model things like `x₁,…,x₃`.
The `@variables` macro supports this with the following syntax:

```julia
julia> @variables x[1:3]
1-element Vector{Vector{Num}}:
 [x₁, x₂, x₃]

julia> @variables y[2:3, 1:5:6] # support for arbitrary ranges and tensors
1-element Vector{Matrix{Num}}:
 [y₂ˏ₁ y₂ˏ₆; y₃ˏ₁ y₃ˏ₆]

julia> @variables t z[1:3](t) # also works for dependent variables
2-element Vector{Any}:
 t
  Num[z₁(t), z₂(t), z₃(t)]
```

Note that `@variables` returns a vector of all the defined variables.

`@variables` can also take runtime symbol values by the `\$` interpolation
operator, and in this case, `@variables` doesn't automatically assign the value,
instead, it only returns a vector of symbolic variables. All the rest of the
syntax also applies here.

```julia
julia> a, b, c = :runtime_symbol_value, :value_b, :value_c
:runtime_symbol_value

julia> vars = @variables t \$a \$b(t) \$c[1:3](t)
3-element Vector{Num}:
      t
 runtime_symbol_value
   value_b(t)
       Num[value_c₁(t), value_c₂(t), value_c₃(t)]

julia> (t, a, b, c)
(t, :runtime_symbol_value, :value_b, :value_c)
```
"""
macro variables(xs...)
    esc(_parse_vars(:variables, Real, xs))
end

TreeViews.hastreeview(x::Sym) = true

function TreeViews.treelabel(io::IO,x::Sym,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Text(x.name))
end


text:
# # DUMMY MODELS FOR TESTING

# ## PARTICLE TRACKER (without methods lifted)

# Consider an object that tracks a particle in one dimension, moving,
# in discrete time, at a speed proportional to the distance away from
# some moving target. The particle is initially at rest at the origin,
# which is where the target also begins.

# Calling `train!` on the object moves it along for the specified
# number of time steps, while calling `ingest!` updates the target
# position.

mutable struct Particle
    position::Float64
    target::Float64
    velocity::Float64
    η::Float64 # learning rate
    training_losses::Vector{Float64}
    Particle(η) = new(0.0, 0.0, 0.0, η, Float64[])
end

Particle(; η=0.1) = Particle(η)

# native train!
function train!(model::Particle)
    model.velocity = model.η*(model.target - model.position)
    model.position = model.position + model.velocity
end

loss(model::Particle) =  abs(model.target - model.position)

function train!(model, n)
    training_losses = map(1:n) do _
        train!(model)
        loss(model)
    end
    model.training_losses = training_losses
    return nothing
end

training_losses(model::Particle) = model.training_losses

function ingest!(model::Particle, target)
    model.target = target
    return nothing
end

# temporary testing area

# model = Particle()
# ingest!(model, 1)
# train!(model, 1)
# @assert loss(model) ≈ 0.9
# ingest!(model, -0.9)
# train!(model, 1)
# @assert loss(model) ≈ 0.9

# model = Particle()
# ingest!(model, 1)
# train!(model, 2)
# @assert training_losses(model) ≈ [0.9, 0.81]


text:
module QuantumLattices
    using Reexport: @reexport

    include("Interfaces.jl")
    include("Prerequisites/Prerequisites.jl")
    include("Mathematics/Mathematics.jl")
    include("Essentials/Essentials.jl")

    @reexport using .Interfaces
    @reexport using .Essentials
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

@testset "TestModel" begin

#---
##

using Test, ACE, JuLIP
using ACE.Testing: ToyModel

#---

@info("Finite difference test of toymodel with different setups")
for species in ( [:Fe, ], [:Fe, :Al], [:Fe, :Ti, :Al] )
   @info("    species = $species")
   V = ToyModel(species)
   at = rand(V, nrepeat=3)
   energy(V, at)
   F = forces(V, at)
   println(@test JuLIP.Testing.fdtest(V, at))
end

#---
##

end


text:
function add!(n::NetworkP,extra::Extracellular)

    extra.v=zeros(Float64,length(n.t))
    coeffs=Array(Extra_coeffs,0)

    for j=1:length(fieldnames(n.neur))-4
        mycoeffs=add_extra(getfield(n.neur,j),extra)
        append!(coeffs,mycoeffs)
    end

    add_extra_(n,extra,coeffs)
end

function add!(n::NetworkP,intra::Intracellular)
    add_intra_(n,intra)
    (ind,newx)=neuron_ind(intra.neur,n.neur.c[intra.mtype])
    l=add_intra(ind,newx,intra.node,length(n.t),getfield(n.neur,intra.mtype))
    push!(n.neur.i,l)
    nothing
end

function add!(n::NetworkP,s::Stim)
    add_stim_(n,s)
    (ind,newx)=neuron_ind(s.neur,n.neur.c[s.mtype])
    l=add_stim(ind,newx,s.node,s.Is,getfield(n.neur,s.mtype))
    push!(n.neur.s,l)
end

function run!(n::NetworkP,init=false)

    #get initial conditions if uninitialized
    if init==true
        initialcon!(n.neur)
    end

    @inbounds main(n.neur)

    for j=1:length(n.extra)
        for k=1:length(fieldnames(n.neur))-4
            for l=1:length(n.neur.c[k])
                n.extra[j].v+=fetch_extra(getfield(n.neur,k),l,j)
            end
        end
    end

    for j=1:length(n.intra)
        @inbounds n.intra[j].v=fetch_intra(getfield(n.neur,n.intra[j].mtype),n.intra[j].mtype,n.neur.i[j])
    end

    reset_t!(n.neur)

    nothing
end

main{T<:Neuron}(n::Array{T,1})=(@inbounds for i=1:length(n);main(n[i]);end;nothing)

main{T<:Puddle}(p::Array{T,1})=(@inbounds main(p[1]))

function main(p::Puddle)

    id=myid()-1

    for p.i in 1:p.imax

        take!(p.r)

        for j=1:length(p.h.s)
            @inbounds p.n[p.h.s[j].neur].rhs[p.h.s[j].node]+=p.h.s[j].Is[p.i]
        end

        main(p.n)

        for j=1:length(p.h.e)
            for k=1:length(p.n)
                @inbounds p.h.e[j].v[p.i]+=a_mult_b(p.h.e[j].coeffs[k],p.n[k])
            end
        end

        for j=1:length(p.h.i)
            @inbounds p.h.i[j].v[p.i]=p.n[p.h.i[j].neur].v[p.h.i[j].node]
        end

        put!(p.r,id)

    end

    nothing
end

initialcon!{T<:Neuron}(n::Array{T,1})=(@inbounds for i=1:length(n);initialcon!(n[i]);end;nothing)

initialcon!{T<:Puddle}(p::Array{T,1})=(@inbounds initialcon!(p[1]); nothing)

initialcon!(p::Puddle)=(initialcon!(p.n); p.i=1)

function neuron_ind(x::Int64,y::Array{UnitRange{Int64},1})

    ind=0
    newx=0
    mysum=0

    for i in y

        if (x<=i[end])&&(x>=i[1])
            newx=x-mysum
            ind+=1
            break
        end

        mysum+=length(i[1])
        ind+=1
    end

    (ind,newx)

end

function get_dims(dims)

    y=Array(UnitRange{Int64},0)

    for i in dims
        push!(y,i[1])
    end
    y
end

function add_intra(p::Int64,neur::Int64,node::Int64,t::Int64,n::DArray)

    remotecall_fetch(((x,neur,node,t)->(push!(localpart(x)[1].h.i,(Intracellular(1,neur,node,zeros(Float64,t))));length(localpart(x)[1].h.i))),p+1,n,neur,node,t)

end

function add_stim(p::Int64,neur::Int64,node::Int64,s::Array{Float64,1},n::DArray)

    remotecall_fetch(((x,neur,node,t)->(push!(localpart(x)[1].h.s,(Stim(t,1,neur,node,0.0,0.0)));length(localpart(x)[1].h.s))),p+1,n,neur,node,s)

end

function fetch_intra(nd::DArray,p::Int64,ind::Int64)
    remotecall_fetch(((n,j)->localpart(n)[1].h.i[j].v),p+1,nd,ind)
end

function add_extra(n::DArray,e::Extracellular)

    mycoeffs=Array(Extra_coeffs,0)

    for p in workers()
        append!(mycoeffs,add_extra(n,e,p))
    end

    mycoeffs
end

function add_extra(n::DArray,e::Extracellular,p::Int64)
    remotecall_fetch(((x,y)->(c=extracellular(localpart(x)[1].n,y);push!(localpart(x)[1].h.e,typeof(e)(e.xyz,c,zeros(Float64,length(e.v))));c)),p,n,e)
end

function fetch_extra(nd::DArray,p::Int64,ind::Int64)
    remotecall_fetch(((n,j)->localpart(n)[1].h.e[j].v),p+1,nd,ind)
end

reset_t!{T<:Puddle}(p::Array{T,1})=(@inbounds for i=1:length(p);reset_t!(p[i]);end;nothing)

reset_t!(p::Puddle)=(p.i=1; nothing)

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
"""
Test costs for paths in an nxn grid lattice. Compare analysis with theory
"""

using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using JLD
using Parameters

# Params
num_trials = 3::Int64
size = 10
mpaths = 2

# performance data
perf_data = []
# Generate network:
net = GridNetwork(size, size)
# Collect performance statistics for 1 userpair and 1 maxpath
performance, dummy, dummy, dummy = net_performance(net, num_trials, 2, max_paths=mpaths)
push!(perf_data, performance)

# Extract data from performance data
loss_arr = collect(map(x->x["loss"], perf_data))
z_arr = collect(map(x->x["Z"], perf_data))

# TODO: This isn't strictly necessary, benchmarking does the averaging for us already.
ave_E = mean(loss_arr)
ave_Z = mean(z_arr)

# Analytic costs:
function ave_pathlength(n)
    return 2/3 * n
end

theory_E = dB_to_P(ave_pathlength(size))
theory_Z = dB_to_Z(ave_pathlength(size))

theory_E2 = dB_to_P(ave_pathlength(size))
theory_Z2 = dB_to_Z(ave_pathlength(size))

theory_2E, theory_2Z = QuNet.purify(theory_E, theory_E2, theory_Z, theory_Z2)

# # Test purification on sample data set
# raw_path_statistics = (12., 11., 9., 3., 3., 10., 8., 8., 5., 4., 3., 12., 4., 9., 10., 1., 13.,
# 3., 8., 6., 3., 5., 10., 7., 8., 4., 10., 8., 3., 5., 8., 9., 8., 1., 5., 6., 10., 6., 6., 12., 12., 2., 10., 9.,
# 5., 3., 10.)
#
# pur_statistics = []
# for pathlength in raw_path_statistics
#     P = dB_to_P(pathlength)
#     Z = dB_to_Z(pathlength)
#     push!(pur_statistics, QuNet.purify_PBS(P, P, Z, Z))
# end
# # println(pur_statistics)
# pur_E_statistics = mean(collect(pur_statistics[i][1] for i in 1:length(pur_statistics)))
# ave_Z_statistics = mean(collect(pur_statistics[i][2] for i in 1:length(pur_statistics)))

println("")
println("ave_E: $ave_E")
println("theory_2E: $theory_2E")
println("ave_Z: $ave_Z")
println("theory_2Z: $theory_2Z")
# println("pur_E_statistics: $ave_E_statistics")
# println("pur_Z_statistics: $ave_Z_statistics")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
module GridapMHDTests

using Test

@time @testset "main" begin include("MainTests.jl") end

@time @testset "Hunt" begin include("HuntTests.jl") end

end # module


text:
#testing for pressure network
cd(@__DIR__)
using Pkg
Pkg.activate(".")

using Flux
using Plots
using BSON: @load
using Random
using Statistics

#load network to test
@load "Models/Pressure_overall_best.bson" Pressure_overall_best
Pressure_Model = Pressure_overall_best
#function for parsing data
function get_data(filename)
    f1 = open(filename)
    tempdata = readlines(f1)
    tempdata[1] = tempdata[1][4:end]  #cutting of beggining bc some weird chars
    close(f1)
    data = []
    regvec = zeros(length(tempdata),
                   length(split(tempdata[1], ",")))
    #parsing datavector to float array... not regularized
    for i in 1:length(tempdata)
        regvec[i,:] =  parse.(Float64,split(tempdata[i], ","))
    end
    return regvec
end

#normalizing data
function norm_data(data, scalingmatrix=nothing)
    #getting regularization values and normalizing data
    regvec = data
    normvec = zeros(length(regvec[:,1]),
                    length(regvec[1,:]))
    #checking if optional scalingmatrix was passed, if not will calc mean and std
    if scalingmatrix == nothing
        scalingmatrix = zeros(length(regvec[1,:]),2)
        for i in 1:length(regvec[1,:])
            scalingmatrix[i,1] = Statistics.mean(regvec[:,i]) #first index is mean
            scalingmatrix[i,2] = Statistics.std(regvec[:,i])  #second index is std
        end
        for i in 1:length(regvec[1,:])
            #normalizing regvec
            normvec[:,i] = (regvec[:,i] .- scalingmatrix[i,1]) ./ scalingmatrix[i,2]
        end
    else
        #for case of normalizing test data by training std and mean
        for i in 1:length(regvec[1,:])
            #normalizing regvec
            normvec[:,i] = (regvec[:,i] .- scalingmatrix[i,1]) ./ scalingmatrix[i,2]
        end
    end
    return normvec, scalingmatrix
end

#function to convert vector of vector to matrix
function vecvec_to_matrix(vecvec)
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = zeros(dim1, dim2)
    for i in 1:dim1
        for j in 1:dim2
            my_array[i,j] = vecvec[i][j]
        end
    end
    return my_array
end

#getting training data
fname = "Dataset/FINAL_DATA.txt"
data_tmp = get_data(fname)
#seperating test and train data
train_data_tmp = []; test_data_tmp = [];
for i in 1:length(data_tmp[:,1])
    if data_tmp[i,2] in [0.3,0.38,0.46,0.53,0.61,0.65] #training case
        push!(train_data_tmp, data_tmp[i,:])
    else #testing case
        push!(test_data_tmp, data_tmp[i,:])
    end
end
train_data = vecvec_to_matrix(train_data_tmp)
test_data = vecvec_to_matrix(test_data_tmp)
#normalizing data
train_norm, scalingmatrix = norm_data(train_data)
test_norm, scalingmatrix = norm_data(test_data,scalingmatrix)
data_norm, scalingmatrix = norm_data(data_tmp,scalingmatrix)

alt_mean = scalingmatrix[1,1]; alt_std = scalingmatrix[1,2];
mach_mean = scalingmatrix[2,1]; mach_std = scalingmatrix[2,2];
fault_mean = scalingmatrix[3,1];fault_std = scalingmatrix[3,2]
pressure_mean = scalingmatrix[4,1]; pressure_std = scalingmatrix[4,2]
#converting data to format which works with network
data = []
for i in 1:length(train_norm[:,1])
    #adjust x and y for fault detection here
    push!(data,(train_norm[i,[1,2,3]],train_norm[i,4]))
end

#Normalized Root Mean Squared Error
function RMSE(data_set,model)
    #calculate root mean squared error
    N = length(data_set[:,1])
    print(N,'\n')
    MSE_sum = 0.0
    for i in 1:N
        MSE_sum += 1/N * (data_set[i,4] - model(data_set[i,1:3])[1])^2.0
    end
    return sqrt(MSE_sum)  # no norm needed bc data is already normalized
end

print("Error Anlysis")
print("\nOverall root mean squared error: ")
print(RMSE(data_norm,Pressure_Model))
print("\nTest root mean squared error: ")
print(RMSE(test_norm,Pressure_Model))
print("\nTrain root mean squared error: ")
print(RMSE(train_norm,Pressure_Model))

#grab all 0 and 0.99 mach values from training set
#TODO: this is like interval of prediction (goes +/- 0.23ish)
#"interval of prediction" for this model is like 10x smaller than fault model
zero_fault = [];  ninetynine_fault = [];

for i in 1:length(train_norm[:,1])
    if train_data[i,3] == 0.0
        push!(zero_fault, train_norm[i,4])
    elseif train_data[i,3] == 0.99
        push!(ninetynine_fault, train_norm[i,4])
    end
end
N = length(zero_fault)
#calculating average percent change along const. mach/alt line
MSE_change = 0.0
for i in 1:N
    global MSE_change += 1/N * (ninetynine_fault[i] - zero_fault[i])^2.0
end
RMSE_change = sqrt(MSE_change)
print("\nAverage RMSE change along constant mach/alt lines is: $RMSE_change")

#Plotting Stuff
plotline = true;altitudeslice = true; #set these to true/false to adjust what is plotted
#lists to be plotted
mach_SCAT = [];pressure_SCAT = [];fault_SCAT = [];model_SCAT = [];
machline = []; pressureline = []; faultline = [];
#specified_alt/mach for altitude slice and plotline respectively

specified_alt_tmp = 1525
specified_alt = (specified_alt_tmp - alt_mean)/alt_std
specified_mach = (0.49 - mach_mean)/mach_std
#organizing plotting data
for i in 1:length(test_norm[:,1])
    mach = test_norm[i,2]
    fault = test_norm[i,3]
    pressure = test_norm[i,4]
    alt = test_norm[i,1]
    if altitudeslice == true
        if alt == specified_alt
            push!(mach_SCAT,mach)
            push!(fault_SCAT,fault)
            push!(pressure_SCAT,pressure)
            push!(model_SCAT,Pressure_Model(test_norm[i,1:3])[1])
        end
    else
        push!(mach_SCAT,mach)
        push!(fault_SCAT,fault)
        push!(pressure_SCAT,pressure)
        push!(model_SCAT,Pressure_Model(test_norm[i,1:3])[1])
    end
    #plot line stuff
    if mach == specified_mach && alt == specified_alt
        push!(machline, mach)
        push!(pressureline, pressure)
        push!(faultline, fault)
    end
end
if altitudeslice == true
    title = "$specified_alt_tmp"
else
    title = "All"
end
#creating 3d scatter for showing network results
scatter(fault_SCAT,mach_SCAT,pressure_SCAT,label="Actual",
        title="Predictions with Testing Data [$title]",
        xlabel="Fault Parameter",
        ylabel="Mach",
        zlabel="Pressure")
scatter!(fault_SCAT,mach_SCAT,model_SCAT,label="Predicted")
#plotline
if plotline == true
    plot!(faultline,machline,pressureline, label="Constant Mach/Alt")
end
current()

max_stars_repo_name: NoahLopezGit/PitotTubeNeuralNetwork

max_issues_repo_name: NoahLopezGit/PitotTubeNeuralNetwork

max_forks_repo_name: NoahLopezGit/PitotTubeNeuralNetwork


text:
module FlexFloat

             #  more is imported than is overloaded
import Base: -, +, *, /, %, ^,
             convert, promote_rule, promote_type, r_promote,
             show, showcompact, string, parse,
             frexp, ldexp,
             isinteger, isfinite, isreal, isinf, isnan,
             sign, signbit, copysign, flipsign, abs,
             decompose,
             sizeof, reinterpret,
             typemin, typemax, realmin, realmax,
             zero, one,  eps,
             sizeof, reinterpret,
             isequal,isless, (!=),(>=),(>),(<=),(<),(==),
             trunc, round, floor, ceil,
             fld, rem, mod, mod1, rem1, fld1,
             min, max, minmax,
             sqrt,exp,expm1,log,log1p,
             sin,cos,tan,csc,sec,cot,
             asin,acos,atan,acsc,asec,acot,
             sinh,cosh,tanh,csch,sech,coth,
             asinh,acosh,atanh,acsch,asech,acoth,
             gamma,lgamma,erf,erfinv,
             start, next, done

export Flex64, Flux64,
       ClCl, ClOp, OpCl, OpOp,
       clcl, clop, opcl, opop,
       widened, narrowed, opened, closed,
       glb, lub, supremum, infimum,
       dia, rad, mid, mig, mag, dev, abs, dist,
       value, isexact, isinexact,
       isclosedclosed, isclosedopen, isopenclosed, isopenopen,
       showmidpointradius, showmidrad

# provide these modular services
include("module/AdjacentFloat/AdjacentFloat.jl")
using .AdjacentFloat
include("module/RoundFast/RoundFast.jl")
using .RoundFast
include("module/Reinterpret/Reinterpret.jl")
using .Reinterpret

typealias Float AbstractFloat

include("type/flexible.jl")
include("type/primitive.jl")
include("type/convert.jl")
include("type/interval.jl")
include("type/compare.jl")
include("type/distill.jl")
include("type/cvtqualia.jl")
include("type/io.jl")

include("arith/valueops.jl")
include("arith/add.jl")
include("arith/sub.jl")
include("arith/mul.jl")
include("arith/div.jl")
include("arith/sqrt.jl")

if isdefined(Main, :CRlibm)
  include("math/moreCRlibm.jl")
  include("math/withCRlibm.jl")
else
  include("math/approxlibm.jl")
end

if isdefined(Main, :Distributions)
    include("math/distributions.jl")
end

end # module


text:
# BEGIN 1.2.1  deprecations
@deprecate get_energy_value(val::BatteryEMS) get_energy_surplus_cost(
    get_operation_cost(val),
)
@deprecate get_penalty_cost(val::BatteryEMS) get_energy_shortage_cost(
    get_operation_cost(val),
)

@deprecate set_energy_value!(val::BatteryEMS, v) set_energy_shortage_cost!(
    get_operation_cost(val),
    v,
)
@deprecate set_penalty_cost!(val::BatteryEMS, v) set_energy_surplus_cost!(
    get_operation_cost(val),
    v,
)

export get_energy_value
export get_penalty_cost
export set_energy_value!
export set_penalty_cost!

# END 1.2.1 deprecations

@deprecate remove_components!(t, sys::System) remove_components!(sys, t)

@deprecate CurrentControl(kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad) VoltageModeControl(
    kpv,
    kiv,
    kffv,
    rv,
    lv,
    kpc,
    kic,
    kffi,
    ωad,
    kad,
)
@deprecate CurrentControl(; kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad) VoltageModeControl(;
    kpv,
    kiv,
    kffv,
    rv,
    lv,
    kpc,
    kic,
    kffi,
    ωad,
    kad,
)