text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Rotations from an Earth-Centered Inertial (ECI) reference frame to another
#   ECI reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export rECItoECI

"""
    rECEFtoECI([T,] ECIo, ECIf, JD_UTC::Number [, eop_data])
    rECEFtoECI([T,] ECIo, JD_UTCo::Number, ECIf, JD_UTCf::Number [, eop_data])

Compute the rotation from an Earth-Centered Inertial (`ECI`) reference frame to
another ECI reference frame. If the origin and destination frame contain only
one *of date* frame, then the first signature is used and `JD_UTC` is the epoch
of this frame. On the other hand, if the origin and destination frame contain
two *of date* frame[^1], e.g. TOD => MOD, then the second signature must be used
in which `JD_UTCo` is the epoch of the origin frame and `JD_UTCf` is the epoch
of the destination frame.

The rotation description that will be used is given by `T`, which can be `DCM`
or `Quaternion`. The origin ECI frame is selected by the input `ECIo` and the
destination ECI frame is selected by the input `ECIf`. The model used to compute
the rotation is specified by the selection of the origin and destination frames.
Currently, there are two models supported: IAU-76/FK5 and IAU-2006 with 2010
conventions (CIO approach only).

[^1]: TEME is an *of date* frame.

# Rotation description

The rotations that aligns the origin ECI frame with the destination ECI frame
can be described by Direction Cosine Matrices or Quaternions. This is selected
by the parameter `T`.

The possible values are:

* `DCM`: The rotation will be described by a Direction Cosine Matrix.
* `Quaternion`: The rotation will be described by a Quaternion.

If no value is specified, then it falls back to `DCM`.

# Conversion model

The model that will be used to compute the rotation is automatically inferred
given the selection of the origin and destination frames. **Notice that mixing
IAU-76/FK5 and IAU-2006/2010 frames is not supported yet.**

# ECI Frame

The supported ECI frames for both origin `ECIo` and destination `ECIf` are:

* `TEME()`: ECI will be selected as the True Equator Mean Equinox (TEME)
            reference frame.
* `TOD()`: ECI will be selected as the True of Date (TOD).
* `MOD()`: ECI will be selected as the Mean of Date (MOD).
* `J2000()`: ECI will be selected as the J2000 reference frame.
* `GCRF()`: ECI will be selected as the Geocentric Celestial Reference Frame
            (GCRF).
* `CIRS()`: ECEF will be selected as the Celestial Intermediate Reference System
            (CIRS).

# EOP Data

The conversion between the frames depends on EOP Data (see `get_iers_eop` and
`read_iers_eop`). If IAU-76/FK5 model is used, then the type of `eop_data` must
be `EOPData_IAU1980`. Otherwise, if IAU-2006/2010 model is used, then the type
of `eop_data` must be `EOPData_IAU2000A`. The following table shows the
requirements for EOP data given the selected frames.

|   Model       |   ECIo  |   ECIf  |    EOP Data   | Function Signature |
|:--------------|:--------|:--------|:--------------|:-------------------|
| IAU-76/FK5    | `GCRF`  | `J2000` | EOP IAU1980   | First              |
| IAU-76/FK5    | `GCRF`  | `MOD`   | EOP IAU1980   | First              |
| IAU-76/FK5    | `GCRF`  | `TOD`   | EOP IAU1980   | First              |
| IAU-76/FK5    | `GCRF`  | `TEME`  | EOP IAU1980   | First              |
| IAU-76/FK5    | `J2000` | `GCRF`  | EOP IAU1980   | First              |
| IAU-76/FK5    | `J2000` | `MOD`   | Not required  | First              |
| IAU-76/FK5    | `J2000` | `TOD`   | Not required  | First              |
| IAU-76/FK5    | `J2000` | `TEME`  | Not required  | First              |
| IAU-76/FK5    | `MOD`   | `GCRF`  | EOP IAU1980   | First              |
| IAU-76/FK5    | `MOD`   | `J2000` | Not required  | First              |
| IAU-76/FK5    | `MOD`   | `TOD`   | Not required  | Second             |
| IAU-76/FK5    | `MOD`   | `TEME`  | Not required  | Second             |
| IAU-76/FK5    | `TOD`   | `GCRF`  | EOP IAU1980   | First              |
| IAU-76/FK5    | `TOD`   | `J2000` | Not required  | First              |
| IAU-76/FK5    | `TOD`   | `MOD`   | Not required  | Second             |
| IAU-76/FK5    | `TOD`   | `TEME`  | Not required  | Second             |
| IAU-76/FK5    | `TEME`  | `GCRF`  | EOP IAU1980   | First              |
| IAU-76/FK5    | `TEME`  | `J2000` | Not required  | First              |
| IAU-76/FK5    | `TEME`  | `MOD`   | Not required  | Second             |
| IAU-76/FK5    | `TEME`  | `TOD`   | Not required  | Second             |
| IAU-2006/2010 | `GCRF`  | `CIRS`  | Not required¹ | First              |
| IAU-2006/2010 | `CIRS`  | `CIRS`  | Not required¹ | Second             |

`¹`: In this case, the terms that account for the free-core nutation and time
dependent effects of the Celestial Intermediate Pole (CIP) position with respect
to the GCRF will not be available, reducing the precision.

## MOD and TOD

In this function, if EOP corrections are not provided, then MOD and TOD frames
will be computed considering the original IAU-76/FK5 theory. Otherwise, the
corrected frame will be used.

# Returns

The rotation description represented by `T` that rotates the origin ECI
reference frame into alignment with the destination ECI reference frame.

# Examples

```julia-repl
julia> eop_IAU1980 = get_iers_eop(:IAU1980);

julia> rECItoECI(DCM, GCRF(), J2000(), DatetoJD(1986, 6, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  1.0          -2.45469e-12   4.56602e-10
  2.45466e-12   1.0          -1.84455e-9
 -4.56602e-10   1.84455e-9    1.0

julia> rECItoECI(Quaternion, TEME(), GCRF(), DatetoJD(1986, 6, 19, 21, 35, 0), eop_IAU1980)
Quaternion{Float64}:
  + 0.9999986335698654 + 1.8300414020900853e-5.i + 0.0006653038276169474.j - 0.0015132396749411375.k

julia> rECItoECI(TOD(), DatetoJD(1986,6,19,21,35,0), TOD(), DatetoJD(1987,5,19,3,0,0), eop_IAU1980)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 1.0          -0.000224087  -9.73784e-5
 0.000224086   1.0          -5.79859e-6
 9.73797e-5    5.77677e-6    1.0

julia> rECItoECI(Quaternion, TOD(), JD_J2000, MOD(), JD_J2000, eop_IAU1980)
Quaternion{Float64}:
  + 0.9999999993282687 - 1.400220690336851e-5.i + 1.3473593746216003e-5.j - 3.107834312843103e-5.k

julia> rECItoECI(J2000(), TEME(), DatetoJD(1986,6,19,21,35,0))
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  0.999995    0.0030265    0.00133055
 -0.00302645  0.999995    -3.86125e-5
 -0.00133066  3.45854e-5   0.999999

julia> eop_IAU2000A = get_iers_eop(:IAU2000A);

julia> rECItoECI(CIRS(), GCRF(), DatetoJD(1986,6,19,21,35,0), eop_IAU2000A)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 0.999999     3.88379e-8  -0.00133066
 7.18735e-9   1.0          3.45882e-5
 0.00133066  -3.45882e-5   0.999999

julia> rECItoECI(Quaternion, CIRS(), GCRF(), DatetoJD(1986,6,19,21,35,0), eop_IAU2000A)
Quaternion{Float64}:
  + 0.9999997785177528 + 1.7294102099105917e-5.i + 0.0006653310148723835.j + 7.912627369563795e-9.k
```
"""
@inline rECItoECI(T_ECIo::T_ECIs, T_ECIf::T_ECIs, JD_UTC::Number,
                  eop_data::EOPData_IAU1980) =
    rECItoECI(DCM, T_ECIo, T_ECIf, JD_UTC, eop_data)

@inline rECItoECI(T_ECIo::T_ECIs_IAU_2006, T_ECIf::T_ECIs_IAU_2006,
                  JD_UTC::Number, eop_data::EOPData_IAU2000A) =
    rECItoECI(DCM, T_ECIo, T_ECIf, JD_UTC, eop_data)

# Specializations for those cases in which we have two *of dates* frames.
@inline rECItoECI(T_ECIo::T_ECIs_of_date, JD_UTCo::Number,
                  T_ECIf::T_ECIs_of_date, JD_UTCf::Number,
                  eop_data::EOPData_IAU1980) =
    rECItoECI(DCM, T_ECIo, JD_UTCo, T_ECIf, JD_UTCf, eop_data)

@inline rECItoECI(T_ECIo::Val{:CIRS}, JD_UTCo::Number,
                  T_ECIf::Val{:CIRS}, JD_UTCf::Number,
                  eop_data::EOPData_IAU2000A) =
    rECItoECI(DCM, T_ECIo, JD_UTCo, T_ECIf, JD_UTCf, eop_data)

# Specializations for those cases that EOP Data is not needed.
@inline rECItoECI(T_ECIo::Val{:J2000},
                  T_ECIf::Union{Val{:MOD}, Val{:TOD}, Val{:TEME}},
                  JD_UTC::Number) =
    rECItoECI(DCM, T_ECIo, T_ECIf, JD_UTC)

@inline rECItoECI(T_ECIo::Union{Val{:MOD}, Val{:TOD}, Val{:TEME}},
                  T_ECIf::Val{:J2000},
                  JD_UTC::Number) =
    rECItoECI(DCM, T_ECIo, T_ECIf, JD_UTC)

@inline rECItoECI(T_ECIo::T_ECIs_of_date, JD_UTCo::Number,
                  T_ECIf::T_ECIs_of_date, JD_UTCf::Number) =
    rECItoECI(DCM, T_ECIo, JD_UTCo, T_ECIf, JD_UTCf)

@inline rECItoECI(T_ECIo::T_ECIs_IAU_2006, T_ECIf::T_ECIs_IAU_2006,
                  JD_UTC::Number) =
    rECItoECI(DCM, T_ECIo, T_ECIf, JD_UTC)

 @inline rECItoECI(T_ECIo::Val{:CIRS}, JD_UTCo::Number,
                   T_ECIf::Val{:CIRS}, JD_UTCf::Number) =
    rECItoECI(DCM, T_ECIo, JD_UTCo, T_ECIf, JD_UTCf)

################################################################################
#                                  IAU-76/FK5
################################################################################

#                                GCRF <=> J2000
# ==============================================================================

function rECItoECI(T::T_ROT, ::Val{:GCRF}, ::Val{:J2000}, JD_UTC::Number,
                   eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    δΔϵ_1980 = eop_data.dEps(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # In this case, we need to convert GCRF back to PEF and then convert to
    # J2000, which is the same conversion from PEF to GCRF **without** the EOP
    # data.
    #
    # TODO: Can I simplify the rotation from TOD with corrections to TOD without
    # corrections?
    r_MOD_GCRF  = rGCRFtoMOD_fk5(T, JD_TT)
    r_PEF_MOD   = rMODtoPEF_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)
    r_MOD_PEF   = rPEFtoMOD_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_J2000_MOD = inv_rotation(r_MOD_GCRF)

    # Compose the full rotation.
    return compose_rotation(r_MOD_GCRF, r_PEF_MOD, r_MOD_PEF, r_J2000_MOD)
end

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:J2000}, T_ECIf::Val{:GCRF},
                  JD_UTC::Number, eop_data::EOPData_IAU1980) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC, eop_data))

#                                 GCRF <=> MOD
# ==============================================================================

function rECItoECI(T::T_ROT, ::Val{:GCRF}, ::Val{:MOD}, JD_UTC::Number,
                   eop_data::EOPData_IAU1980)

    # Get the time in TT.
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Return the rotation.
    return rGCRFtoMOD_fk5(T, JD_TT)
end

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:MOD}, T_ECIf::Val{:GCRF},
                  JD_UTC::Number, eop_data::EOPData_IAU1980) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC, eop_data))

#                                 GCRF <=> TOD
# ==============================================================================

function rECItoECI(T::T_ROT, ::Val{:GCRF}, ::Val{:TOD}, JD_UTC::Number,
                   eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in TT.
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    δΔϵ_1980 = eop_data.dEps(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # Return the rotation.
    r_MOD_GCRF = rGCRFtoMOD_fk5(T, JD_TT)
    r_TOD_MOD  = rMODtoTOD_fk5(T, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compose the full rotation.
    return compose_rotation(r_MOD_GCRF, r_TOD_MOD)
end

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:TOD}, T_ECIf::Val{:GCRF},
                  JD_UTC::Number, eop_data::EOPData_IAU1980) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC, eop_data))

#                                GCRF <=> TEME
# ==============================================================================

function rECItoECI(T::T_ROT, ::Val{:GCRF}, ::Val{:TEME}, JD_UTC::Number,
                   eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in TT.
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    δΔϵ_1980 = eop_data.dEps(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # Return the rotation.
    r_MOD_GCRF = rGCRFtoMOD_fk5(T, JD_TT)
    r_TEME_MOD = rMODtoTEME(T, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compose the full rotation.
    return compose_rotation(r_MOD_GCRF, r_TEME_MOD)
end

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:TEME}, T_ECIf::Val{:GCRF},
                  JD_UTC::Number, eop_data::EOPData_IAU1980) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC, eop_data))

#                                J2000 <=> MOD
# ==============================================================================

function rECItoECI(T::T_ROT, ::Val{:J2000}, ::Val{:MOD}, JD_UTC::Number,
                   eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    δΔϵ_1980 = eop_data.dEps(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # In this case, we need to convert J2000 back to PEF and then convert to
    # MOD. This is necessary because we need to apply EOP corrections to convert
    # to MOD and just a `rGCRFtoMOD_fk5` would yield a frame according to the
    # original IAU-76/FK5 theory.
    #
    # TODO: Can I simplify this rotation?
    r_MOD_J2000 = rGCRFtoMOD_fk5(T, JD_TT)
    r_PEF_MOD   = rMODtoPEF_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_MOD_PEF   = rPEFtoMOD_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compose the full rotation.
    return compose_rotation(r_MOD_J2000, r_PEF_MOD, r_MOD_PEF)
end

function rECItoECI(T::T_ROT, ::Val{:J2000}, ::Val{:MOD}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Compute and return the rotation.
    r_MOD_J2000 = rGCRFtoMOD_fk5(T, JD_TT)

    return r_MOD_J2000
end

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:MOD}, T_ECIf::Val{:J2000},
                  JD_UTC::Number, eop_data::EOPData_IAU1980) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC, eop_data))

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:MOD}, T_ECIf::Val{:J2000},
                  JD_UTC::Number) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC))

#                                J2000 <=> TOD
# ==============================================================================

function rECItoECI(T::T_ROT, ::Val{:J2000}, ::Val{:TOD}, JD_UTC::Number,
                   eop_data::EOPData_IAU1980)

    arcsec2rad = π/648000

    # Get the time in UT1 and TT.
    JD_UT1 = JD_UTCtoUT1(JD_UTC, eop_data)
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    #
    # TODO: The difference is small, but should it be `JD_TT` or `JD_UTC`?
    δΔϵ_1980 = eop_data.dEps(JD_UTC)*arcsec2rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC)*arcsec2rad

    # In this case, we need to convert J2000 back to PEF and then convert to
    # TOD. This is necessary because we need to apply EOP corrections to convert
    # to TOD.
    #
    # TODO: Can I simplify this rotation?
    r_MOD_J2000 = rGCRFtoMOD_fk5(T, JD_TT)
    r_PEF_MOD   = rMODtoPEF_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_MOD_PEF   = rPEFtoMOD_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)
    r_TOD_MOD   = rMODtoTOD_fk5(T, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compose the full rotation.
    return compose_rotation(r_MOD_J2000, r_PEF_MOD, r_MOD_PEF, r_TOD_MOD)
end

function rECItoECI(T::T_ROT, ::Val{:J2000}, ::Val{:TOD}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Compute and return the composed rotation.
    r_MOD_J2000 = rGCRFtoMOD_fk5(T, JD_TT)
    r_TOD_MOD   = rMODtoTOD_fk5(T, JD_TT, 0, 0)

    return compose_rotation(r_MOD_J2000, r_TOD_MOD)
end

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:TOD}, T_ECIf::Val{:J2000},
                  JD_UTC::Number, eop_data::EOPData_IAU1980) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC, eop_data))

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:TOD}, T_ECIf::Val{:J2000},
                  JD_UTC::Number) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC))

#                                J2000 <=> TEME
# ==============================================================================

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:J2000}, T_ECEFf::Val{:TEME},
                  JD_UTC::Number, eop_data::EOPData_IAU1980) =
    rECItoECI(T, T_ECIo, T_ECEFf, JD_UTC)

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:TEME}, T_ECEFf::Val{:J2000},
                  JD_UTC::Number, eop_data::EOPData_IAU1980) =
    rECItoECI(T, T_ECIo, T_ECEFf, JD_UTC)

function rECItoECI(T::T_ROT, ::Val{:J2000}, ::Val{:TEME}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT  = JD_UTCtoTT(JD_UTC)

    # Return the rotation.
    return rGCRFtoTEME(T, JD_TT, 0, 0)
end

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:TEME}, T_ECIf::Val{:J2000},
                  JD_UTC::Number) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC))

#                          Between MOD, TOD, and TEME
# ==============================================================================

function rECItoECI(T::T_ROT, T_ECIo::T_ECIs_of_date, JD_UTCo::Number,
                   T_ECIf::T_ECIs_of_date, JD_UTCf::Number,
                   eop_data::EOPData_IAU1980)

    # In this case, we convert origin to GCRF and then convert back to the
    # destination. This is necessary because the user may want to change the
    # epoch.
    r_GCRF_ECIo = rECItoECI(T, T_ECIo,     Val(:GCRF), JD_UTCo, eop_data)
    r_ECIf_GCRF = rECItoECI(T, Val(:GCRF), T_ECIf,     JD_UTCf, eop_data)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end

function rECItoECI(T::T_ROT, T_ECIo::T_ECIs_of_date, JD_UTCo::Number,
                   T_ECIf::T_ECIs_of_date, JD_UTCf::Number)

    # In this case, in which we do not have EOP data, we convert origin to J2000
    # and then convert back to the destination. This is necessary because the
    # user may want to change the epoch.
    r_GCRF_ECIo = rECItoECI(T, T_ECIo,      Val(:J2000), JD_UTCo)
    r_ECIf_GCRF = rECItoECI(T, Val(:J2000), T_ECIf,      JD_UTCf)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end

################################################################################
#                                IAU-2006/2010
################################################################################

#                                GCRF <=> CIRS
# ==============================================================================

function rECItoECI(T::T_ROT, ::Val{:GCRF}, ::Val{:CIRS}, JD_UTC::Number,
                   eop_data::EOPData_IAU2000A)

    arcsec2rad = π/648000

    # Get the time in TT.
    JD_TT = JD_UTCtoTT(JD_UTC)

    # Get the EOP data related to the desired epoch.
    dX = eop_data.dX(JD_UTC)*arcsec2rad
    dY = eop_data.dY(JD_UTC)*arcsec2rad

    # Compute and return the rotation.
    return rGCRFtoCIRS_iau2006(T, JD_TT, dX, dY)
end

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:CIRS}, T_ECIf::Val{:GCRF},
                  JD_UTC::Number, eop_data::EOPData_IAU2000A) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC, eop_data))

function rECItoECI(T::T_ROT, ::Val{:GCRF}, ::Val{:CIRS}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT = JD_UTCtoTT(JD_UTC)

    # Compute and return the rotation.
    return rGCRFtoCIRS_iau2006(T, JD_TT)
end

@inline rECItoECI(T::T_ROT, T_ECIo::Val{:CIRS}, T_ECIf::Val{:GCRF},
                  JD_UTC::Number) =
    inv_rotation(rECItoECI(T, T_ECIf, T_ECIo, JD_UTC))

#                                 Between CIRS
# ==============================================================================

function rECItoECI(T::T_ROT, T_ECIo::Val{:CIRS}, JD_UTCo::Number,
                   T_ECIf::Val{:CIRS}, JD_UTCf::Number,
                   eop_data::EOPData_IAU2000A)

    # In this case, we convert origin to GCRF and then convert back to the
    # destination. This is necessary because the user may want to change the
    # epoch.
    r_GCRF_ECIo = rECItoECI(T, T_ECIo,     Val(:GCRF), JD_UTCo, eop_data)
    r_ECIf_GCRF = rECItoECI(T, Val(:GCRF), T_ECIf,     JD_UTCf, eop_data)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end

function rECItoECI(T::T_ROT, T_ECIo::Val{:CIRS}, JD_UTCo::Number,
                   T_ECIf::Val{:CIRS}, JD_UTCf::Number)

    # In this case, we convert origin to GCRF and then convert back to the
    # destination. This is necessary because the user may want to change the
    # epoch.
    r_GCRF_ECIo = rECItoECI(T, T_ECIo,     Val(:GCRF), JD_UTCo)
    r_ECIf_GCRF = rECItoECI(T, Val(:GCRF), T_ECIf,     JD_UTCf)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end


text:
#=

  Crypto problem (alphametic) in Julia ConstraintSolver.jl

  This is a standard alphametic problem in mathematical recreations,
  constraint programming etc.

  From GLPK:s model cryto.mod.

  """
  This problem comes from the newsgroup rec.puzzle.
  The numbers from 1 to 26 are assigned to the letters of the alphabet.
  The numbers beside each word are the total of the values assigned to
  the letters in the word (e.g. for LYRE: L, Y, R, E might be to equal
  5, 9, 20 and 13, or any other combination that add up to 47).
  Find the value of each letter under the equations:

  BALLET  45     GLEE  66     POLKA      59     SONG     61
  CELLO   43     JAZZ  58     QUARTET    50     SOPRANO  82
  CONCERT 74     LYRE  47     SAXOPHONE 134     THEME    72
  FLUTE   30     OBOE  53     SCALE      51     VIOLIN  100
  FUGUE   50     OPERA 65     SOLO       37     WALTZ    34

  Solution:
  A, B,C, D, E,F, G, H, I, J, K,L,M, N, O, P,Q, R, S,T,U, V,W, X, Y, Z
  5,13,9,16,20,4,24,21,25,17,23,2,8,12,10,19,7,11,15,3,1,26,6,22,14,18

  Reference:
  Koalog Constraint Solver <http://www.koalog.com/php/jcs.php>,
  Simple problems, the crypto-arithmetic puzzle ALPHACIPHER.
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function crypto(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    BALLET     =  45
    CELLO      =  43
    CONCERT    =  74
    FLUTE      =  30
    FUGUE      =  50
    GLEE       =  66
    JAZZ       =  58
    LYRE       =  47
    OBOE       =  53
    OPERA      =  65
    POLKA      =  59
    QUARTET    =  50
    SAXOPHONE  = 134
    SCALE      =  51
    SOLO       =  37
    SONG       =  61
    SOPRANO    =  82
    THEME      =  72
    VIOLIN     = 100
    WALTZ      =  34

    n = 26
    @variable(model, 1 <= x[1:n] <= n, Int)
    A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z = x

    @constraint(model, x in CS.AllDifferent())

    @constraint(model,B + A + L + L + E + T == BALLET)
    @constraint(model,C + E + L + L + O == CELLO)
    @constraint(model,C + O + N + C + E + R + T == CONCERT)
    @constraint(model,F + L + U + T + E == FLUTE)
    @constraint(model,F + U + G + U + E == FUGUE)
    @constraint(model,G + L + E + E == GLEE)
    @constraint(model,J + A + Z + Z == JAZZ)
    @constraint(model,L + Y + R + E == LYRE)
    @constraint(model,O + B + O + E == OBOE)
    @constraint(model,O + P + E + R + A == OPERA)
    @constraint(model,P + O + L + K + A == POLKA)
    @constraint(model,Q + U + A + R + T + E + T == QUARTET)
    @constraint(model,S + A + X + O + P + H + O + N + E == SAXOPHONE)
    @constraint(model,S + C + A + L + E == SCALE)
    @constraint(model,S + O + L + O == SOLO)
    @constraint(model,S + O + N + G == SONG)
    @constraint(model,S + O + P + R + A + N + O == SOPRANO)
    @constraint(model,T + H + E + M + E == THEME)
    @constraint(model,V + I + O + L + I + N == VIOLIN)
    @constraint(model,W + A + L + T + Z == WALTZ)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

@time crypto()

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.7.0
#     language: julia
#     name: julia-1.7
# ---

# %%
using LinearAlgebra
using Distributions
using Plots
using Random
using Optim
using QuadGK

"""最小二乗法による多項式フィッティング"""
function polynomialfit(x, y, d)
    n = length(y)   # データサイズ
    X = x .^ (0:d)' # design matrix
    β̂ = X \ y       # Juliaでは X \ y の一発で最小二乗法が可能
    f(x) = evalpoly(x, β̂)
    ŷ = f.(x)
    r = ŷ - y
    σ̂ = norm(r)/√n  # モデルの標準偏差パラメータの推定
    dist = product_distribution([Normal(μ, σ̂) for μ in ŷ])
    loglik = logpdf(dist, y)  # 最大化された対数尤度
    AIC = -2loglik + 2(d + 2) # 所謂赤池情報量規準
    (; x, y, d, β̂, σ̂, f, loglik, AIC)
end

"""optimize函数を使った最尤法による多項式フィッティング"""
function polynomialfit_optim(x, y, d)
    negloglik(w) = begin
        β = @view w[1:end-1]
        σ = exp(w[end])
        n = length(y)
        dist = product_distribution([Normal(evalpoly(x[i], β), σ) for i in 1:n])
        -loglikelihood(dist, y)
    end
    o = optimize(negloglik, zeros(d+2), Optim.Options(store_trace=true, extended_trace = true, iterations = 10^4))
    ŵ = o.minimizer
    β̂ = ŵ[1:end-1]
    σ̂ = exp(ŵ[end])
    loglik = -o.minimum
    AIC = 2o.minimum + 2(d + 2)
    f(x) = evalpoly(x, β̂)
    (; x, y, d, β̂, σ̂, f, loglik, AIC, negloglik, o)
end

function plotpolyfit(result; legend=:topleft, titlefontsize=10, kwargs...)
    (; x, y, d, f, loglik, AIC) = result
    xs = range(extrema(x)...; length=1000)
    plot(; legend, titlefontsize)
    scatter!(x, y; label="data")
    plot!(xs, f.(xs); label="degree-$d fitting")
    title!("degree = $d,  log likelihood = $(round(loglik; digits=2)),  AIC = $(round(AIC; digits=2))")
    plot!(; kwargs...)
end

# %%
# データの生成

rng = MersenneTwister(1234567890)

n = 21
x = range(-1, 1; length=n)
σ₀ = 0.3
f₀(x) = sinpi(x)
noise = σ₀*randn(rng, length(x))
y = f₀.(x) + noise

xs = range(extrema(x)...; length=1000)
P = plot(; legend=:topleft)
scatter!(x, y; label="data")
plot!(xs, sinpi.(xs); label="sin(πx)")
title!("data and true law"; titlefontsize=10)
plot!(; size=(400, 300))

# %%
β₀ = [isodd(k) ? (-1)^(k÷2)*π^k/factorial(k) : 0.0 for k in 0:3]
@show β₀ σ₀;

# %%
(; β̂, σ̂, loglik, AIC) = polynomialfit(x, y, 3)
@show β̂ σ̂ loglik AIC;

# %%
(; β̂, σ̂, loglik, AIC) = polynomialfit_optim(x, y, 3)
@show β̂ σ̂ loglik AIC;

# %%
Q = plotpolyfit(polynomialfit(x, y, 3))
R = plotpolyfit(polynomialfit(x, y, 9))
S = plotpolyfit(polynomialfit(x, y, 15))

plot(P, Q, R, S; size=(800, 600), layout=(2, 2), titlefontsize=9)

# %% [markdown]
# 対数尤度(log likelihood)の大きさは「もっともらしさ」の度合いではない！
#
# 対数尤度の大きさはモデルのデータへのフィッティングの良さの指標の1つに過ぎない。
#
# 上の結果を見れば分かるように、フィッティングで使う多項式の次数を上げた方がフィッティングの度合いが上がって、対数尤度も大きくなるが、次数が大きな場合の結果は全然「もっともらしい」とは言えない。
#
# 尤度を最大化することによってモデルのパラメータ値を調節する最尤法は「最も尤もらしいパラメータ値」を探すのではなく、単に「モデルをデータに最もフィットさせるようなパラメータ値」を探しているだけである。
#
# 最尤法は「最も尤もらしいパラメータ値」を探す方法ではない！

# %%
function plot_loglikAIC(x, y)
    result = polynomialfit.(Ref(x), Ref(y), 0:18)
    d = getproperty.(result, :d)
    loglik = getproperty.(result, :loglik)
    AIC = getproperty.(result, :AIC)
    dbest = d[argmin(AIC)]

    T = plot(d, loglik; label="log likelihood", legend=:topleft)
    plot!(; xtick=d, xlabel="degree", ylabel="log likelihood")
    title!("log likelihoods of polynomial fitting"; titlefontsize=10)

    U = plot(d, AIC; label="AIC", legend=:topright)
    plot!(; xtick=d, xlabel="degree", ylabel="AIC")
    title!("AICs of polynomial fitting (minimized at degree $dbest)"; titlefontsize=10)

    plot(T, U; size=(800, 300), titlefontsize=9, leftmargin=3Plots.mm, bottommargin=3Plots.mm)
end

# %%
plot_loglikAIC(x, y)

# %%
function animate_polyfit(x, y, d, f₀, σ₀; gifname="polyfit$d.gif", fps=20, ylim=(-10, 45), ytick=-100:5:100)
    (; x, y, d, negloglik, o) = polynomialfit_optim(x, y, d)
    a, b = extrema(x)
    n = length(y)
    value = -getproperty.(o.trace, :value)
    centroid = getindex.(getproperty.(o.trace, :metadata), "centroid")

    ngenerr = Vector{Float64}(undef, length(centroid))
    Threads.@threads for t in keys(centroid)
        β = centroid[t][1:end-1]
        f(x) = evalpoly(x, β)
        σ = exp(centroid[t][end])
        GE = (1/2)*(log(2π*σ^2) + σ₀^2/σ^2 + (1/σ^2)*1/(b-a)*quadgk(x -> (f(x) - f₀(x))^2, a, b)[1])
        ngenerr[t] = n*GE
    end

    xs = range(extrema(x)...; length=1000)
    ylim = isnothing(ylim) ? extrema([-value; ngenerr]) .+ (-1, 1) : ylim

    L = length(value)
    tstep = max(1, L ÷ 400)

    anim = @animate for t in [fill(1, fps); 1:tstep:L; fill(L, 3fps)]
        β = centroid[t][1:end-1]
        f(x) = evalpoly(x, β)
        P = plot(; legend=false)
        scatter!(x, y)
        plot!(xs, f.(xs); ylim=(-1.5, 1.5))
        title!("degree-$d fitting of size-$n data"; titlefontsize=10)

        Q = plot(; legend=:topleft)
        plot!(-value[1:t]; label="(-1)×(log likelihood)")
        plot!(ngenerr[1:t]; label="n×(generalization error)")
        plot!(;xlim=extrema(axes(value, 1)).+(-0.05L, 0.05L), ylim, ytick)
        title!("t = $t"; titlefontsize=10)

        plot(P, Q; size=(800, 400))    end
    gif(anim, gifname; fps)
end

# %%
animate_polyfit(x, y, 1, f₀, σ₀)

# %%
animate_polyfit(x, y, 3, f₀, σ₀)

# %%
animate_polyfit(x, y, 9, f₀, σ₀)

# %%
animate_polyfit(x, y, 15, f₀, σ₀)

# %%
GE₀ = (1/2)*(log(2π*σ₀^2) + σ₀^2/σ₀^2)
n * GE₀

# %%
rng = MersenneTwister(1234567891)
y = f₀.(x) + σ₀*randn(rng, length(x))
plot_loglikAIC(x, y)

# %%
animate_polyfit(x, y, 3, f₀, σ₀; gifname="polyfit3a.gif")

# %%
animate_polyfit(x, y, 5, f₀, σ₀; gifname="polyfit5a.gif")

# %%
rng = MersenneTwister(1234567892)
y = f₀.(x) + σ₀*randn(rng, length(x))
plot_loglikAIC(x, y)

# %%
animate_polyfit(x, y, 3, f₀, σ₀; gifname="polyfit3b.gif")

# %%
animate_polyfit(x, y, 5, f₀, σ₀; gifname="polyfit5b.gif")

# %%
rng = MersenneTwister(1234567893)
y = f₀.(x) + σ₀*randn(rng, length(x))
plot_loglikAIC(x, y)

# %%
animate_polyfit(x, y, 3, f₀, σ₀; gifname="polyfit3c.gif")

# %%
animate_polyfit(x, y, 5, f₀, σ₀; gifname="polyfit5c.gif")

# %%
animate_polyfit(x, y, 9, f₀, σ₀; gifname="polyfit9c.gif")

# %%
rng = MersenneTwister(1234567894)
y = f₀.(x) + σ₀*randn(rng, length(x))
plot_loglikAIC(x, y)

# %%
animate_polyfit(x, y, 3, f₀, σ₀; gifname="polyfit3d.gif")

# %%
animate_polyfit(x, y, 5, f₀, σ₀; gifname="polyfit5d.gif")

# %%
rng = MersenneTwister(1234567895)
y = f₀.(x) + σ₀*randn(rng, length(x))
plot_loglikAIC(x, y)

# %%
animate_polyfit(x, y, 3, f₀, σ₀; gifname="polyfit3e.gif")

# %%
animate_polyfit(x, y, 5, f₀, σ₀; gifname="polyfit5e.gif")

# %%
rng = MersenneTwister(1234567896)
y = f₀.(x) + σ₀*randn(rng, length(x))
plot_loglikAIC(x, y)

# %%
animate_polyfit(x, y, 3, f₀, σ₀; gifname="polyfit3f.gif")

# %%
animate_polyfit(x, y, 5, f₀, σ₀; gifname="polyfit5f.gif")

# %%
rng = MersenneTwister(1234567897)
y = f₀.(x) + σ₀*randn(rng, length(x))
plot_loglikAIC(x, y)

# %%
animate_polyfit(x, y, 3, f₀, σ₀; gifname="polyfit3g.gif")

# %%
animate_polyfit(x, y, 5, f₀, σ₀; gifname="polyfit5g.gif")

# %%
rng = MersenneTwister(1234567898)
y = f₀.(x) + σ₀*randn(rng, length(x))
plot_loglikAIC(x, y)

# %%
animate_polyfit(x, y, 3, f₀, σ₀; gifname="polyfit3h.gif")

# %%
animate_polyfit(x, y, 5, f₀, σ₀; gifname="polyfit5h.gif")

# %%
rng = MersenneTwister(1234567899)
y = f₀.(x) + σ₀*randn(rng, length(x))
plot_loglikAIC(x, y)

# %%
animate_polyfit(x, y, 3, f₀, σ₀; gifname="polyfit3i.gif")

# %%
animate_polyfit(x, y, 5, f₀, σ₀; gifname="polyfit5i.gif")

# %%

max_stars_repo_path: 0025/simple example of overfitting.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0025/simple example of overfitting.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0025/simple example of overfitting.jl

max_forks_repo_name: genkuroki/public


text:
besselh0(z) = besselj0(z) + im*bessely0(z)
besselh1(z) = besselj1(z) + im*bessely1(z)

"Evaluate the Helmholtz 2D single layer potential kernel"
function hh_slp(x, y, wavenumber)
    T = eltype(x)
    if norm(x-y) < eps(T)
        return zero(Complex{T})
    else
        return im * besselh0(wavenumber*norm(x-y)) / 4
    end
end

"Evaluate the Helmholtz 2D double layer potential kernel"
hh_dlp(x, y, wavenumber, normal_y, z = norm(x-y)) =
    -(normal_y' * (x-y)) * im * wavenumber/4 * besselh1(wavenumber*z) / z

"Diagonal limiting value of the Helmholtz 2D double layer potential kernel"
function hh_dlp_limit(x, y, wavenumber, normal_y, grad_y, grad_d_y)
    -(grad_d_y[1] * grad_y[2] - grad_y[1] * grad_d_y[2]) / (4*norm(grad_y)^3*pi) +0im
end

"Evaluate the Helmholtz 2D double layer potential kernel"
function hh_dlp_kernel(x, y, tau, wavenumber, param)
    z = norm(x-y)
    if abs(z) > 100*eps(eltype(x))
        hh_dlp(x, y, wavenumber, normal(param, tau), z)
    else
        # Evaluate the limiting value
        hh_dlp_limit(x, y, wavenumber, normal(param, tau),
            gradient(param, tau), gradient_derivative(param, tau))
    end
end

hh_slpt(x, y, wavenumber, normal_x) = error("to do")
hh_hypersingular(x, y, wavenumber, normal_x, normal_y) = error("to do")

"Supertype of all Helmholtz kernels."
abstract type HelmholtzKernel <: BoundaryKernel end

"Wavenumber of the Helmholtz equation"
wavenumber(op::HelmholtzKernel) = op.wavenumber

export Helmholtz_SLP_2D
"The single layer potential kernel of the Helmholtz equation."
struct Helmholtz_SLP_2D{T} <: HelmholtzKernel
    wavenumber  ::  T
end

Helmholtz_SLP_2D(wavenumber::Integer) = Helmholtz_SLP_2D(float(wavenumber))

(kernel::Helmholtz_SLP_2D)(t, tau, param, x, y) = hh_slp(x, y, wavenumber(kernel))
(kernel::Helmholtz_SLP_2D)(x, tau, param, y) = hh_slp(x, y, wavenumber(kernel))

is_symmetric(::Helmholtz_SLP_2D) = true

singularity(::Helmholtz_SLP_2D) = LogSingularDiagonal()

export Helmholtz_DLP_2D
"The double layer potential kernel of the Helmholtz equation in 2D."
struct Helmholtz_DLP_2D{T} <: HelmholtzKernel
    wavenumber  ::  T
end

Helmholtz_DLP_2D(wavenumber::Integer) = Helmholtz_DLP_2D(float(wavenumber))

BasisFunctions.name(kernel::Helmholtz_DLP_2D) = "2D Helmholtz double layer potential kernel"

(kernel::Helmholtz_DLP_2D)(t, tau, param, x, y) =
    hh_dlp_kernel(x, y, tau, wavenumber(kernel), param)
(kernel::Helmholtz_DLP_2D)(x, tau, param, y) =
    hh_dlp_kernel(x, y, tau, wavenumber(kernel), param)

is_symmetric(::Helmholtz_DLP_2D) = false

# The kernel is continuous but its derivative is singular
singularity(::Helmholtz_DLP_2D) = LogSingularDiagonal()

## A few standard boundary conditions

export make_boundary_condition_planewave,
    make_parboundary_condition_planewave,
    make_boundary_condition_pointsource,
    make_parboundary_condition_pointsource

boundary_condition_wave(x, wavenumber, direction, amplitude) =
    amplitude * exp(im*wavenumber*dot(direction, x))

"Return a field function that evaluates to a plane wave."
make_boundary_condition_planewave(wavenumber, direction, amplitude = 1) =
    (x,y) -> boundary_condition_wave(SVector(x,y), wavenumber, direction, amplitude)

"Return a parametric boundary function that evaluates to a plane wave."
make_parboundary_condition_planewave(param, wavenumber, direction, amplitude = 1) =
    t -> boundary_condition_wave(applymap(param, t), wavenumber, direction, amplitude)

boundary_condition_pointsource(x, center, wavenumber, amplitude) =
    amplitude * besselh(0, 1, wavenumber*norm(x-center))

"Return a field function that evaluates to a point source."
make_boundary_condition_pointsource(center, wavenumber, amplitude = 1) =
    (x,y) -> boundary_condition_pointsource(SVector(x,y), center, wavenumber, amplitude)

"Return a parametric boundary function that evaluates to a point source."
make_parboundary_condition_pointsource(param, center, wavenumber, amplitude = 1) =
    t -> boundary_condition_pointsource(applymap(param, t), center, wavenumber, amplitude)


text:
#####################################
#     Define Auxillary functions
#####################################

function readin(n::Int64)
	println("Begin read-in and formatting for n = ",n,"...")
	@load "data"*string(n)*".jld2" K M negator dof_map
	println("Finsihes read-in and formatting for n = ",n)
	println("--------------------------------------------")
	return K::SparseMatrixCSC{Float64,Int64}, M, negator, dof_map
end

function interpolate_mesh(sol_coarse::Array{Float64,1},n_coarse::Int64,n_fine::Int64)
	# Input: coarse solution: sol_coarse = (y,p,lambda)
	#	 coarse DoF-map: dof_coarse ordered by x-coordinate first, y-coordinate second with stepsize h
	#	 Fine   DoF-map: dof_fine   ordered by x-coordinate first, y-coordinate second with stepsize h/2
	#
	# Assumptions: both discretizations are uniform square grids in the unit square with seq. numbering
	#
	# Output: Linear interpolation of the coarse solution on the fine mesh.
	dof_fine = load("data"*string(n_fine)*".jld2", "dof_map")
	dof_coarse = load("data"*string(n_coarse)*".jld2", "dof_map")
	#Find stepsizes:
	h_coarse = dof_coarse[2,:][1] - dof_coarse[1,:][1]
	h_fine = h_coarse/2

	# Find number of nodes per size (assuming unit square)
	side_coarse = Int64(1/h_coarse + 1)
	side_fine =  Int64(1/h_fine + 1)

	# With which we can calculate the number of points
	n_block_coarse = side_coarse^2
	n_block_fine = side_fine^2

	# This allows us to allocate memory for the interpolated solution
	y_fine = zeros(side_fine^2) # State
	p_fine = zeros(side_fine^2) # adjoint
	l_fine = zeros(side_fine^2) # multiplier

	# Now loop over the coarse DoF-handler
	# This loop will transfer the known values from the coarse
	# grid to the fine grid and interpolate the coarse points between them

	# to map the coarse grid to the fine we need to keep track of the indices
	indexadjust=0 # This variable will ensure dof_coarse[i]=dof_fine[i+indexadjust]
	for i = 1:n_block_coarse
		#println("i = ",i)
		# Current x,y values:
		x_coarse = dof_coarse[i,:][1]
		y_coarse = dof_coarse[i,:][2]
		# Extract current values from coarse solution
		#statetemp = sol_coarse[i]
		#adjointtemp = sol_coarse[i+n_block_coarse]
		#lagrangetemp = sol_coarse[i+2*n_block_coarse]

		# In each step we will transfer the known value and
		# construct at most 3 points in the interpolated solution
		# If possible we interpolate to the right, up and right/up diagonally
		# ----------
		# Thus, in each step we need to check wether we are in the right conrner or top of the coarse grid
		# if neither: we send 4 values into the coarse grid and update indexajust+=1
		# if at the right: we only interpolate up and update indexajust += side_fine as we move to the next row
		# If we are at the top and not at the right: we interpolate right and update indexadjust +=1

		# Check if we are at the top (in course grid), i.e. y=1, we use h_fine
		if( y_coarse + h_fine > 1)
			# Now we are at the top, check if we are also on the right i.e. x=1
			if( x_coarse + h_fine > 1)
				# This is the final point we need only transfer the final value.
				y_fine[i+indexadjust] = sol_coarse[i]
				p_fine[i+indexadjust] = sol_coarse[i+n_block_coarse]
				l_fine[i+indexadjust] = sol_coarse[i+2*n_block_coarse]
			else
				# We are at the top but not at the far right so we need oo interpolate
				# one step to the right and set indexadjust += 1

				# Transfer known value:
				y_fine[i+indexadjust] = sol_coarse[i]
				p_fine[i+indexadjust] = sol_coarse[i+n_block_coarse]
				l_fine[i+indexadjust] = sol_coarse[i+2*n_block_coarse]

				# Interpolate one step to the right:
				y_fine[i+indexadjust+1] = (sol_coarse[i] + sol_coarse[i+1])/2
				p_fine[i+indexadjust+1] = (sol_coarse[i+n_block_coarse] + sol_coarse[i+1+n_block_coarse])/2
				l_fine[i+indexadjust+1] = (sol_coarse[i+2*n_block_coarse] + sol_coarse[i+1+2*n_block_coarse])/2

				#Update indexadjust +=1
				indexadjust = 1 + indexadjust
			end

		# If not at the top check if we are at the right:
		elseif( x_coarse + h_fine > 1 )
			# Transfer known value:
			y_fine[i+indexadjust] = sol_coarse[i]
			p_fine[i+indexadjust] = sol_coarse[i+n_block_coarse]
			l_fine[i+indexadjust] = sol_coarse[i+2*n_block_coarse]

			# Interpolate up
			y_fine[i+indexadjust+side_fine] = (sol_coarse[i]+sol_coarse[i+side_coarse])/2
			p_fine[i+indexadjust+side_fine] = (sol_coarse[i+n_block_coarse]+sol_coarse[i+side_coarse+n_block_coarse])/2
			l_fine[i+indexadjust+side_fine] = (sol_coarse[i+2*n_block_coarse]+sol_coarse[i+side_coarse+2*n_block_coarse])/2

			#Update indexadjust += n_block_fine
			indexadjust = indexadjust + side_fine

		# We are neither at the nor the right, i.e. this is a "normal" point
		#and we should transfer 4 values to the fine grid
		else
			# Transfer known value:
			y_fine[i+indexadjust] = sol_coarse[i]
			p_fine[i+indexadjust] = sol_coarse[i+n_block_coarse]
			l_fine[i+indexadjust] = sol_coarse[i+2*n_block_coarse]

			# Interpolate up:
			y_fine[i+indexadjust+side_fine] = (sol_coarse[i]+sol_coarse[i+side_coarse])/2
			p_fine[i+indexadjust+side_fine] = (sol_coarse[i+n_block_coarse]+sol_coarse[i+side_coarse+n_block_coarse])/2
			l_fine[i+indexadjust+side_fine] = (sol_coarse[i+2*n_block_coarse]+sol_coarse[i+side_coarse+2*n_block_coarse])/2

			# Interpolate one step to the right:
			y_fine[i+indexadjust+1] = (sol_coarse[i] + sol_coarse[i+1])/2
			p_fine[i+indexadjust+1] = (sol_coarse[i+n_block_coarse] + sol_coarse[i+1+n_block_coarse])/2
			l_fine[i+indexadjust+1] = (sol_coarse[i+2*n_block_coarse] + sol_coarse[i+1+2*n_block_coarse])/2

			# Interpolate diagonally up-right:
			y_fine[i+indexadjust+side_fine+1] = ( sol_coarse[i] + sol_coarse[i+1] + sol_coarse[i+side_coarse] + sol_coarse[i+side_coarse+1] )/4
			p_fine[i+indexadjust+side_fine+1] = ( sol_coarse[i+n_block_coarse] + sol_coarse[i+1+n_block_coarse] + sol_coarse[i+side_coarse+n_block_coarse] + sol_coarse[i+side_coarse+1+n_block_coarse] )/4
			l_fine[i+indexadjust+side_fine+1] = ( sol_coarse[i+2*n_block_coarse] + sol_coarse[i+1+2*n_block_coarse] + sol_coarse[i+side_coarse+2*n_block_coarse] + sol_coarse[i+side_coarse+1+2*n_block_coarse] )/4
			#Update indexadjust +=1
			indexadjust = 1 + indexadjust

		end

	end
	return [y_fine; p_fine ; l_fine]
end

function diag_sqrt_inv(A)
	# Returns inv(sqrt(A)) for a diagonal matrix A
	n_block,m_block = size(A)
	A_sqrt_inv = zeros(n_block)
	for i=1:n_block
		A_sqrt_inv[i]=1/sqrt(A[i,i])
	end
	return spdiagm(0 => A_sqrt_inv)
end

function diag_sqrt(A)
	# Returns (sqrt(A)) for a diagonal matrix A
	n_block,m_block = size(A)
	A_sqrt = zeros(n_block)
	for i=1:n_block
		A_sqrt[i]=sqrt(A[i,i])
	end
	return spdiagm(0 => A_sqrt)
end

function desired_state(func::Function,dof_map::Array{Float64,2})
	# Evaluates desired state on Dofs and returns vector
	n_block::Int64 = size(dof_map)[1]
	yd::Array{Float64,1} = zeros(n_block)
	for i=1:n_block
		xtemp,ytemp=dof_map[i,:]
		yd[i] = func(xtemp,ytemp)
	end
	return yd
end

#########################################
########           END           ########
#########################################

max_stars_repo_name: dravinsi/demoPDEOPT

max_issues_repo_name: dravinsi/demoPDEOPT

max_forks_repo_name: dravinsi/demoPDEOPT


text:
include(joinpath(pwd(), "examples/direct_policy_optimization/simulate.jl"))

using Random

# Unpack trajectories
x̄, ū = unpack(z̄, prob)
x, u = unpack(z, prob)
θ = get_policy(z, prob_dpo)

# Simulation setup
model_sim = Quadrotor{RK3, FixedTime}(n, m, d,
                  mass,
                  J,
                  Jinv,
                  g,
                  L,
                  kf,
                  km)

x1_sim = copy(x1)
T_sim = 10 * T

using Random
Random.seed!(1)

W = Distributions.MvNormal(zeros(model_sim.n),
	Diagonal(1.0e-5 * ones(model_sim.n)))
w = rand(W, T_sim)

tf_nom = sum([ū[t][end] for t = 1:T-1])
t_nom = range(0, stop = tf_nom, length = T)
t_sim_nom = range(0, stop = tf_nom, length = T_sim)

tf_dpo = sum([u[t][end] for t = 1:T-1])
t_nom_dpo = range(0, stop = tf_dpo, length = T)
t_sim_nom_dpo = range(0, stop = tf_dpo, length = T_sim)

dt_sim_nom = tf_nom / (T_sim - 1)
dt_sim_dpo = tf_dpo / (T_sim - 1)

# Simulate
N_sim = 100

Random.seed!(1)
W0 = Distributions.MvNormal(zeros(model_sim.n),
	Diagonal(1.0e-3 * ones(model_sim.n)))
w0 = rand(W0, N_sim)

z_lqr1 = []
u_lqr1 = []
J_lqr1 = []
Jx_lqr1 = []
Ju_lqr1 = []

for i = 1:N_sim
	z0_sim = vec(copy(x1_sim) + w0[:, i])
	println("sim: $i")
	@show z0_sim[1:3]
	_z_lqr1, _u_lqr1, _J_lqr1, _Jx_lqr1, _Ju_lqr1 = simulate(
		model_sim,
		policy, K,
	    x̄, ū,
		Q, R,
		T_sim, ū[1][end],
		z0_sim,
		w,
		_norm = 2,
		ul = ul_nom, uu = uu1,
		u_idx = (1:model.m - 1))

	push!(z_lqr1, _z_lqr1)
	push!(u_lqr1, _u_lqr1)
	push!(J_lqr1, _J_lqr1)
	push!(Jx_lqr1, _Jx_lqr1)
	push!(Ju_lqr1, _Ju_lqr1)
end

@show mean(J_lqr1)
@show std(J_lqr1)
@show mean(Jx_lqr1)
@show std(Jx_lqr1)
@show mean(Ju_lqr1)
@show std(Ju_lqr1)

z_lqr2 = []
u_lqr2 = []
J_lqr2 = []
Jx_lqr2 = []
Ju_lqr2 = []

for i = 1:N_sim
	z0_sim = vec(copy(x1_sim) + w0[:, i])
	println("sim: $i")
	@show z0_sim[1:3]

	_z_lqr2, _u_lqr2, _J_lqr2, _Jx_lqr2, _Ju_lqr2 = simulate(
		model_sim,
		policy,
		K,
	    x̄, ū,
		Q, R,
		T_sim, ū[1][end],
		z0_sim,
		w,
		_norm = 2,
		ul = ul_nom,
		uu = uu2,
		u_idx = (1:model.m - 1))

	push!(z_lqr2, _z_lqr2)
	push!(u_lqr2, _u_lqr2)
	push!(J_lqr2, _J_lqr2)
	push!(Jx_lqr2, _Jx_lqr2)
	push!(Ju_lqr2, _Ju_lqr2)
end

@show mean(J_lqr2)
@show std(J_lqr2)
@show mean(Jx_lqr2)
@show std(Jx_lqr2)
@show mean(Ju_lqr2)
@show std(Ju_lqr2)

z_lqr3 = []
u_lqr3 = []
J_lqr3 = []
Jx_lqr3 = []
Ju_lqr3 = []

for i = 1:N_sim
	z0_sim = vec(copy(x1_sim) + w0[:, i])
	println("sim: $i")
	@show z0_sim[1:3]

	_z_lqr3, _u_lqr3, _J_lqr3, _Jx_lqr3, _Ju_lqr3 = simulate(
		model_sim,
		policy, K,
	    x̄, ū,
		Q, R,
		T_sim, ū[1][end],
		z0_sim,
		w,
		_norm = 2,
		ul = ul_nom, uu = uu3,
		u_idx = (1:model.m - 1))

	push!(z_lqr3, _z_lqr3)
	push!(u_lqr3, _u_lqr3)
	push!(J_lqr3, _J_lqr3)
	push!(Jx_lqr3, _Jx_lqr3)
	push!(Ju_lqr3, _Ju_lqr3)
end
idx_finite = isfinite.(J_lqr3)
@show mean(J_lqr3[idx_finite])
@show std(J_lqr3[idx_finite])
@show mean(Jx_lqr3[idx_finite])
@show std(Jx_lqr3[idx_finite])
@show mean(Ju_lqr3[idx_finite])
@show std(Ju_lqr3[idx_finite])

z_lqr4 = []
u_lqr4 = []
J_lqr4 = []
Jx_lqr4 = []
Ju_lqr4 = []

for i = 1:N_sim
	z0_sim = vec(copy(x1_sim) + w0[:, i])
	println("sim: $i")
	@show z0_sim[1:3]

	_z_lqr4, _u_lqr4, _J_lqr4, _Jx_lqr4, _Ju_lqr4 = simulate(
		model_sim,
		policy, K,
	    x̄, ū,
		Q, R,
		T_sim, ū[1][end],
		z0_sim, w,
		_norm = 2,
		ul = ul_nom, uu = uu4,
		u_idx = (1:model.m - 1))

	push!(z_lqr4, _z_lqr4)
	push!(u_lqr4, _u_lqr4)
	push!(J_lqr4, _J_lqr4)
	push!(Jx_lqr4, _Jx_lqr4)
	push!(Ju_lqr4, _Ju_lqr4)
end
idx_finite = isfinite.(J_lqr4)
convert.(Int, range(1, stop = 100, length = 100)[.!isfinite.(J_lqr4)])
@show mean(J_lqr4[idx_finite])
@show std(J_lqr4[idx_finite])
@show mean(Jx_lqr4[idx_finite])
@show std(Jx_lqr4[idx_finite])
@show mean(Ju_lqr4[idx_finite])
@show std(Ju_lqr4[idx_finite])

z_dpo1 = []
u_dpo1 = []
J_dpo1 = []
Jx_dpo1 = []
Ju_dpo1 = []

for i = 1:N_sim
	z0_sim = vec(copy(x1_sim) + w0[:, i])
	println("sim: $i")
	@show z0_sim[1:3]

	_z_dpo1, _u_dpo1, _J_dpo1, _Jx_dpo1, _Ju_dpo1 = simulate(
		model_sim,
		policy, θ,
	    x, u,
		Q, R,
		T_sim, u[1][end],
		z0_sim, w,
		_norm = 2,
		ul = ul_nom, uu = uu1,
		u_idx = (1:model.m - 1))

	push!(z_dpo1, _z_dpo1)
	push!(u_dpo1, _u_dpo1)
	push!(J_dpo1, _J_dpo1)
	push!(Jx_dpo1, _Jx_dpo1)
	push!(Ju_dpo1, _Ju_dpo1)
end

@show mean(J_dpo1)
@show std(J_dpo1)
@show mean(Jx_dpo1)
@show std(Jx_dpo1)
@show mean(Ju_dpo1)
@show std(Ju_dpo1)

z_dpo2 = []
u_dpo2 = []
J_dpo2 = []
Jx_dpo2 = []
Ju_dpo2 = []

for i = 1:N_sim
	z0_sim = vec(copy(x1_sim) + w0[:, i])
	println("sim: $i")
	@show z0_sim[1:3]
	_z_dpo2, _u_dpo2, _J_dpo2, _Jx_dpo2, _Ju_dpo2 = simulate(
		model_sim,
		policy, θ,
	    x, u,
		Q, R,
		T_sim, u[1][end],
		z0_sim, w,
		_norm = 2,
		ul = ul_nom, uu = uu2,
		u_idx = (1:model.m - 1))

	push!(z_dpo2, _z_dpo2)
	push!(u_dpo2, _u_dpo2)
	push!(J_dpo2, _J_dpo2)
	push!(Jx_dpo2, _Jx_dpo2)
	push!(Ju_dpo2, _Ju_dpo2)
end

@show mean(J_dpo2)
@show std(J_dpo2)
@show mean(Jx_dpo2)
@show std(Jx_dpo2)
@show mean(Ju_dpo2)
@show std(Ju_dpo2)

z_dpo3 = []
u_dpo3 = []
J_dpo3 = []
Jx_dpo3 = []
Ju_dpo3 = []

for i = 1:N_sim
	z0_sim = vec(copy(x1_sim) + w0[:, i])
	println("sim: $i")
	@show z0_sim[1:3]

	_z_dpo3, _u_dpo3, _J_dpo3, _Jx_dpo3, _Ju_dpo3 = simulate(
		model_sim,
		policy, θ,
	    x, u,
		Q, R,
		T_sim, u[1][end],
		z0_sim, w,
		_norm = 2,
		ul = ul_nom, uu = uu3,
		u_idx = (1:model.m - 1))
	push!(z_dpo3, _z_dpo3)
	push!(u_dpo3, _u_dpo3)
	push!(J_dpo3, _J_dpo3)
	push!(Jx_dpo3, _Jx_dpo3)
	push!(Ju_dpo3, _Ju_dpo3)
end

@show mean(J_dpo3)
@show std(J_dpo3)
@show mean(Jx_dpo3)
@show std(Jx_dpo3)
@show mean(Ju_dpo3)
@show std(Ju_dpo3)

z_dpo4 = []
u_dpo4 = []
J_dpo4 = []
Jx_dpo4 = []
Ju_dpo4 = []

for i = 1:N_sim
	z0_sim = vec(copy(x1_sim) + w0[:, i])
	println("sim: $i")
	@show z0_sim[1:3]

	_z_dpo4, _u_dpo4, _J_dpo4, _Jx_dpo4, _Ju_dpo4 = simulate(
		model_sim,
		policy, θ,
	    x, u,
		Q, R,
		T_sim, u[1][end],
		z0_sim, w,
		_norm = 2,
		ul = ul_nom, uu = uu4,
		u_idx = (1:model.m - 1))
	push!(z_dpo4, _z_dpo4)
	push!(u_dpo4, _u_dpo4)
	push!(J_dpo4, _J_dpo4)
	push!(Jx_dpo4, _Jx_dpo4)
	push!(Ju_dpo4, _Ju_dpo4)
end

@show mean(J_dpo4)
@show std(J_dpo4)
@show mean(Jx_dpo4)
@show std(Jx_dpo4)
@show mean(Ju_dpo4)
@show std(Ju_dpo4)


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Qualification Round 2014
# Problem B: Cookie Clicker Alpha

const initial_rate = 2.0

# Calculates the amount of time (in seconds) needed to get x cookies when farms
# cost c cookies and increase production by f, if using the optimal strategy.
function besttime(c, f, x)
	t = 0
	r = initial_rate
	while x/r > c/r + x/(r+f)
		t += c/r
		r += f
	end
	t + x/r
end

function main()
	n = int(readline())
	for i = 1:n
		c, f, x = float(split(readline()))
		t = besttime(c, f, x)
		println("Case #$i: ", t)
	end
end

main()


text:
# phantomrestore index, phantomgetarray

function AE_phantomrestore(line::Ref{UInt}, ns::Namespace, ph::AetherlangObject{AetherlangPhantomCollection}, inx::AetherlangObject)::AetherlangObject
    if inx.dataref < 1 || inx.dataref > length(ph.dataref.ref)
        throw(AetherlangError("`phantomrestore` cannot restore index $(inx.dataref)"))
    end
    nph::AetherlangPhantomCollection = copy(ph.dataref)
    try
        pop!(nph.addindx, inx.dataref)
    catch
    end
    AetherlangObject(nph)
end

function AE_phantomrestore(line::Ref{UInt}, ns::Namespace, ph::AetherlangObject{AetherlangPhantomCollection})::AetherlangObject
    AetherlangObject(ph.dataref.ref)
end

betterphantoms_namespace_modify = Namespace(
    "phantomrestore" => AetherlangObject(AE_phantomrestore)
)

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
# CAVEAT: do not forget this!
#
# -> TestProblem does not belong to the NLS_Solver module.
#    But one must overload these functions, hence we import them.
#
import NLS_Solver: parameter_size, residue_size, eval_r, eval_r_J

include("Rosenbrock.jl")
include("PowellSingular.jl")


text:
module HDF5Plugins
    include("TestUtils.jl")

    using .TestUtils

    function register()
    end

    function __init__()
        if get(ENV, "JULIA_HDF5PLUGINS_REGISTER", "1") != "0"
            register()
        end
    end
end


text:
using MRIReco

##### simple example ####

N = 256
I = shepp_logan(N)
I = circularShutterFreq!(I,1)
cmap = 1im*quadraticFieldmap(N,N,125*2pi)

# simulation parameters
params = Dict{Symbol, Any}()
params[:simulation] = "fast"
params[:trajName] = "Spiral"
params[:numProfiles] = 1
params[:numSamplingPerProfile] = N*N
params[:windings] = 128
params[:AQ] = 3.0e-2
params[:correctionMap] = cmap[:,:,1]

# do simulation
acqData = simulation(I, params)

# reco parameters
params = Dict{Symbol, Any}()
params[:reco] = "direct"
params[:reconSize] = (N,N)
Ireco = reconstruction(acqData, params)

# export image
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/fieldmapReco1.png")
exportImage(filename, abs.(Ireco[:,:,1,1,1]) )

# now with fieldmap correction
params[:reco] = "standard"
params[:regularization] = "L2"
params[:iterations] = 3
params[:solver] = "admm"
params[:correctionMap] = cmap
params[:alpha] = 1.75
params[:m] = 4.0
params[:K] = 28

Ireco = reconstruction(acqData, params)

# export images
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/fieldmapReco2.png")
exportImage(filename, abs.(Ireco[:,:,1,1,1]) )

# export fieldmap image
cmap_ = imag(cmap[:,:,1]) .+ 125*2pi
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/fieldmap.png")
exportImage(filename, abs.(cmap_) )


text:
using ExpmV
using Base.Test

for i = 1:100, d in 10:10:60, herm in [true, false]

    d = 60

    r = sprandn(d,d,.1)+1im*sprandn(d,d,.1)
    if herm
        r = (r-r')/2
    end

    rv = randn(d)+1im*randn(d)
    rv = rv / norm(rv,2)

    rt = randn()

    x = expmv(rt,r,rv)

    @test_approx_eq_eps norm(x-expm(full(rt*r))*rv,2) 0.0 1e-9
end


text:
module Markowitz
using LinearAlgebra

export markowitz, unit_sum, add_constraint, frontier, smooth, optimal

function simplex(mu,A,b,L,U,epsilon=1e-10)
    n = size(mu,1)
    m = size(A,1)
    Ai = zeros(m,m)
    L = [ L; zeros(m) ]
    U = [ U; Inf*ones(m) ]
    z = zeros(n+m)
    x = zeros(n+m)
    state = Vector{Symbol}(undef,n+m)
    nartificial = m
    adjRate = Vector{Float64}(undef,n+m)
    phase = 1
    OUT = 1:n
    state[OUT] .= :LOW
    x[OUT] = L[1:n]
    z[OUT] .= 0
    IN = n .+ (1:m)
    tmp = b - A * L[1:n]
    for i in 1:m
        Ai[i,i] = (tmp[i]<0 ? -1 : 1)
    end
    A = hcat(A, Ai)
    state[IN] .= :IN
    x[IN] = abs.(tmp)
    z[IN] .= -1
    loop = 0
    while loop < 10*(n+m)
        loop = loop + 1
        price = -z[IN]' * Ai[1:length(IN),:]
        profit = ( z[OUT]' + price * A[:,OUT]) .* map(x-> x==:HIGH ? -1 : 1, state[OUT])'
        if maximum(profit) < epsilon
            if phase == 1
                error("unfeasible or degenerate")
            else
                muChanged = false
                for i in 1:length(OUT)
                    if profit[i] > -epsilon
                        muChanged = true
                        @warn "tweaking mu[$i] to ensure the solution is unique"
                        mu[OUT[i]] = mu[OUT[i]] + epsilon * (state[OUT[i]] == :LOW ? -1 : 1)
                    end
                end
                return(Dict(:IN => IN, :OUT => OUT, :state => state[1:n],
                            :x => x[1:n], :Ai => Ai, :mu => muChanged ? mu : nothing))
            end
        end
        enters = OUT[argmax(profit)[2]]
        directionIn = (state[enters]==:HIGH ? (:DOWN) : (:UP))
        adjRate = -Ai * A[:,enters] * (directionIn==:DOWN ? -1 : 1)
        exits = enters
        exitsIdx = NaN
        directionOut = directionIn
        if ( enters > n || U[enters] == Inf )
            theta = Inf
        else
            theta = U[enters]-L[enters]
        end
        for i in 1:m
            if adjRate[i] < -epsilon
                tmp = ( L[IN[i]] - x[IN[i]] ) /adjRate[i]
                if tmp < theta
                    theta = tmp
                    exits = IN[i]
                    exitsIdx = i
                    directionOut = :DOWN
                end
            elseif (adjRate[i] > epsilon) && (U[IN[i]] < Inf)
                tmp = ( U[IN[i]] - x[IN[i]] ) / adjRate[i]
                if tmp < theta
                    theta = tmp
                    exits = IN[i]
                    exitsIdx = i
                    directionOut = :UP
                end
            end
        end
        if theta == Inf
            error("simplex unbounded")
        end
        x[IN] = x[IN] + theta * adjRate
        x[enters] = x[enters] + theta * (directionIn==:UP ? 1 : -1)
        OUT = OUT[OUT.!=enters]
        IN = sort([ enters; IN ])
        state[enters] = :IN
        IN = IN[IN.!=exits]
        if exits <= n
            OUT = sort([ exits; OUT ])
        end
        state[exits] = (directionOut==:UP ? (:HIGH) : (:LOW))
        if !isnan(exitsIdx)
            for i in 1:m
                if i != exitsIdx
                    Ai[i,:] = Ai[i,:] - adjRate[i] / adjRate[exitsIdx] * Ai[exitsIdx,:]
                end
            end
            Ai[exitsIdx,:] = Ai[exitsIdx,:] / adjRate[exitsIdx] * (directionIn==:UP ? -1 : 1)
            entersIdx = findall(IN .== enters)[1]
            if entersIdx > exitsIdx
                tmp = Ai[exitsIdx,:]
                for i in exitsIdx:(entersIdx-1)
                    Ai[i,:] = Ai[i+1,:]
                end
                Ai[entersIdx,:]= tmp
            elseif entersIdx < exitsIdx
                tmp = Ai[exitsIdx,:]
                for i in exitsIdx:-1:(entersIdx+1)
                    Ai[i,:] = Ai[i-1,:]
                end
                Ai[entersIdx,:] = tmp
            end
        end
        if phase == 1 && exits > n
            nartificial = nartificial-1
            if nartificial == 0
                phase = 2
                z = mu
            end
        end
    end
end

struct Change
    inout::Symbol
    direction::Symbol
    variable::Int
    lambda::Float64
end

function cla(C,mu,A,b,L,U,epsilon=1e-10,lambdaTarget=1e-6)
    @assert det(C)>=0
    scale = maximum(C)
    C = C/scale
    n = size(mu,1)
    m = size(A,1)
    alpha = zeros(n+m)
    beta = zeros(n+m)
    xi = zeros(n+m)
    bbar = zeros(n+m)
    extAi = zeros(n,m)
    M = [C     A'    ;
         A zeros(m,m)]
    lambdaE = Inf
    lambdaEold = NaN
    weights = Array{Float64}(undef,0,n)
    ret = Vector{Float64}(undef,0)
    vol = Vector{Float64}(undef,0)
    lambda = Vector{Float64}(undef,0)
    simp = simplex(mu,A,b,L,U,epsilon)
    if simp[:mu] != nothing; mu = simp[:mu]; end
    x = simp[:x]
    IN = simp[:IN]
    OUT = simp[:OUT]
    state = simp[:state]
    Ai = simp[:Ai]
    alpha[[ state.!=:IN; falses(m) ]] = x[state.!=:IN]
    IN = [IN; n.+(1:m)]
    for in in IN
        bbar[in] = (in <= n ? 0 : b[in-n]) - dot(M[in,OUT], x[OUT])
    end
    extAi[IN[1:m],1:m] = Ai[1:m,1:m]
    Mi = [ zeros(n,n)           extAi            ;
           extAi'     -Ai'*C[IN[1:m],IN[1:m]]*Ai ]
    loop = 0
    E = undef
    V = undef
    a0 = undef
    a1 = undef
    a2 = undef
    while (loop < 10*(n+m))
        loop = loop + 1
        changes = []
        for in in IN
            alpha[in] = dot(Mi[in,IN], bbar[IN])
            beta[in] = dot(Mi[in,IN[IN.<=n]], mu[IN[IN.<=n]])
            if in <= n
                if beta[in] > epsilon
                    push!(changes,Change(:OUT,:DOWN,in,(L[in]-alpha[in])/beta[in]))
                elseif (U[in] < Inf) && (beta[in] < -epsilon)
                    push!(changes,Change(:OUT,:UP,in,(U[in]-alpha[in])/beta[in]))
                end
            end
        end
        for out in OUT
            gamma = dot(M[out,:], alpha)
            delta = -mu[out] + dot(M[out,:], beta)
            if (state[out] == :LOW) && (delta > epsilon)
                push!(changes,Change(:IN,:UP,out,-gamma/delta))
            elseif (state[out] == :HIGH) && (delta < -epsilon)
                push!(changes,Change(:IN,:DOWN,out,-gamma/delta))
            end
        end
        changes = sort(changes,by=x->-x.lambda)
        if length(changes) > 1 && changes[1].lambda == changes[2].lambda
            @warn "more than one variable change could be done, will continue but note that the degenerate case is not properly handled"
        end
        lambdaEold = lambdaE
        if length(changes) > 0
            change = changes[1]
            changes = []
            lambdaE = max(0,change.lambda)
        else
            lambdaE = 0
        end
        x[IN[1:(length(IN)-m)]] = alpha[IN[1:(length(IN)-m)]] + lambdaE * beta[IN[1:(length(IN)-m)]]
        dEdLambdaE = beta[IN[1:(length(IN)-m)]]' * mu[IN[1:(length(IN)-m)]]
        if dEdLambdaE < epsilon
            E = x[1:n]' * mu[1:n]
            V = x[1:n]' * C * x[1:n]
        else
            a2 = 1/dEdLambdaE
            a1 = 2*(lambdaEold-a2*E)
            a0 = V-a1*E-a2*E*E
            E = E + (lambdaE-lambdaEold)*dEdLambdaE
            V = a0 + a1*E + a2*E*E
        end
        weights = vcat(weights,x')
        push!(lambda,lambdaE)
        push!(ret,E)
        push!(vol,sqrt(V*scale))
        if lambdaE < lambdaTarget
            return([weights,ret,vol,lambda])
        end
        if change.inout == :OUT
            out = change.variable
            alpha[out] = x[out]
            beta[out] = 0
            IN = IN[IN.!=out]
            if out <= n
                OUT = sort([ out; OUT ])
            end
            state[out] = (change.direction==:UP ? (:HIGH) : (:LOW))
            Mi = Mi - Mi[:,[out]] * Mi[[out],:] / Mi[out,out]
            bbar[IN] = bbar[IN] - M[IN,out] * x[out]
        else
            in = change.variable
            xi[IN] = Mi[IN,IN] * M[IN,[in]]
            xij = M[in,in] - dot(M[in,IN], xi[IN])
            Mi[IN,IN] = Mi[IN,IN] + xi[IN] * xi[IN]' / xij
            Mi[IN,in] = -xi[IN]/xij
            Mi[in,IN] = -xi[IN]/xij
            Mi[in,in] = 1/xij
            bbar[IN] = bbar[IN] + M[IN,in] * x[in]
            OUT = OUT[OUT.!=in]
            IN = sort([ in; IN ])
            state[in] = :IN
            bbar[in] = -dot(M[in,OUT], x[OUT])
        end
    end
end

mutable struct Form0
    V::Matrix{Float64}
    E::Vector{Float64}
    lower::Vector{Float64}
    upper::Vector{Float64}
    names::Vector{String}
    lhs::Matrix{Float64}
    op::Vector{Char}
    rhs::Vector{Float64}
end

function Base.display(m::Form0)
    if length(m.names) > 0
        println("Assets")
        Base.display(m.names)
    end
    println("Covariance matrix")
    Base.display(m.V)
    println("Expected returns")
    Base.display(m.E')
    println("Asset bounds")
    Base.display(vcat(m.lower', m.upper'))
    println("Additional constraints")
    Base.display(hcat(m.lhs, m.op, m.rhs))
end

function markowitz(E::Union{Vector,Matrix},V::Matrix;lower=0,upper=Inf,names=[])
    n = length(E)
    @assert size(V) == (n,n)
    if length(lower) == 1
        lower = lower*ones(n)
    end
    if length(upper) == 1
        upper = upper*ones(n)
    end
    Form0(V,vec(E),vec(lower),vec(upper),vec(names),
          Matrix{Float64}(undef,0,n),Vector{Symbol}(undef,0),Vector{Float64}(undef,0))
end

function add_constraint(m::Form0,lhs::Vector,op::Char,rhs::Number)
    add_constraint(m,reshape(lhs,1,length(lhs)),op,rhs)
end

function add_constraint(m::Form0,lhs::Matrix,op::Char,rhs::Number)
    m.lhs = vcat(m.lhs, lhs)
    push!(m.op, op)
    push!(m.rhs, rhs)
    m
end

function unit_sum(m::Form0)
    add_constraint(m, ones(length(m.E)), '=', 1)
end

struct Form2
    V::Matrix{Float64}
    E::Vector{Float64}
    lower::Vector{Float64}
    upper::Vector{Float64}
    A::Matrix{Float64}
    b::Vector{Float64}
end

function Base.display(m::Form2)
    println("Covariance matrix")
    Base.display(m.V)
    println("Expected returns")
    Base.display(m.E')
    println("Asset bounds")
    Base.display(vcat(m.lower', m.upper'))
    println("Equality constraints")
    Base.display(hcat(m.A, NaN*ones(length(m.b)), m.b))
end

function form2(m::Form0)
    n = length(m.E)
    nc = length(m.rhs)
    slack = sum(m.op.!='=')
    E = zeros(n+slack)
    E[1:n] = m.E
    V = zeros(n+slack,n+slack)
    V[1:n,1:n] = m.V
    lower = zeros(n+slack)
    lower[1:n] = m.lower
    upper = Inf*ones(n+slack)
    upper[1:n] = m.upper
    A = zeros(nc,n+slack)
    col = n+1
    for row in 1:nc
        A[row,1:n] = m.lhs[row,:]
        if m.op[row] != '='
            if m.op[row] == '<'
                A[row,col] = 1
            elseif m.op[row] == '>'
                A[row,col] = -1
            else
                error("invalid constraint type ",m.op[row])
            end
            col = col+1
        end
    end
    b = m.rhs
    Form2(V,E,lower,upper,A,b)
end

struct Frontier
    weights::Matrix{Float64}
    ret::Vector{Float64}
    vol::Vector{Float64}
    lambda::Vector{Float64}
    problem::Union{Form0,Form2}
end

function simplex(m::Form2)
    simplex(m.E,m.A,m.b,m.lower,m.upper)
end

function frontier(m::Form2,form0=nothing)
    weights,ret,vol,lambda = cla(m.V,m.E,m.A,m.b,m.lower,m.upper)
    if form0 == nothing
        Frontier(weights,ret,vol,lambda,m)
    else
        Frontier(weights[:,1:length(form0.E)],ret,vol,lambda,form0)
    end
end

function simplex(m::Form0)
    simplex(form2(m))
end

function frontier(m::Form0)
    frontier(form2(m),m)
end

function smooth(f,N=100)
    frontier = [ f.vol[1] f.ret[1] ]
    lengths = sqrt.((diff(f.ret)/(f.ret[1]-f.ret[end])).^2+(diff(f.vol)/(f.vol[1]-f.vol[end])).^2)
    nsegments = ceil.(N*lengths/sum(lengths))
    for i in 1:length(nsegments)
        w0 = f.weights[i,:]
        w1 = f.weights[i+1,:]
        if nsegments[i] > 0
            for p in (1:nsegments[i])/nsegments[i]
                w = w0 + p*(w1-w0)
                frontier = vcat(frontier, [ sqrt(w'*f.problem.V*w) f.problem.E'*w ])
            end
        end
    end
    frontier
end

struct Portfolio
    vol::Float64
    ret::Float64
    weights::Vector{Float64}
end

function optimal(f::Frontier,ret=nothing)
    if ret == nothing
        Portfolio(f.vol[end],f.ret[end],f.weights[end,:])
    else
        idx = findfirst(x-> x<ret, f.ret)
        if idx < 2
            error("the target return is not inside the frontier! min: ", f.ret[end], " max: ", f.ret[1])
        else
            p = (f.ret[idx-1]-ret)/(f.ret[idx-1]-f.ret[idx])
            w = f.weights[idx-1,:] + p*(f.weights[idx,:]-f.weights[idx-1,:])
            Portfolio(sqrt(w'*f.problem.V*w), f.problem.E'*w, w)
        end
    end
end

end


text:
using DataFrames

"""
    pic_diagnostics_transfer( particle_group, kernel_smoother, efield_dofs)

Compute ``\\sum_{particles} w  p e(x_p) ``

- `particle_group`
- `kernel_smoother` : Kernel smoother
- `efield_dofs` : coefficients of efield

"""
function pic_diagnostics_transfer(
    particle_group::ParticleGroup{1,1}, kernel_smoother, efield_dofs::Vector{Float64}
)
    transfer = 0.0

    @inbounds for i_part in 1:(particle_group.n_particles)
        xi = particle_group.array[1, i_part]
        vi = particle_group.array[2, i_part]
        wi =
            particle_group.array[3, i_part] *
            particle_group.charge *
            particle_group.common_weight
        efield = evaluate(kernel_smoother, xi, efield_dofs)

        transfer += (vi * efield) * wi
    end

    return transfer
end

export TimeHistoryDiagnosticsSpin

"""
    TimeHistoryDiagnosticsSpin( particle_group, maxwell_solver,
                            kernel_smoother_0, kernel_smoother_1 )

Context to save and plot diagnostics

- `particle_group` : Particles data
- `maxwell_solver` : Maxwell solver
- `kernel_smoother_0` : Mesh coupling operator
- `kernel_smoother_1` : Mesh coupling operator
- `data` : DataFrame containing time history values
"""
struct TimeHistoryDiagnosticsSpin
    particle_group::ParticleGroup
    maxwell_solver::Maxwell1DFEM
    kernel_smoother_0::ParticleMeshCoupling1D
    kernel_smoother_1::ParticleMeshCoupling1D
    data::DataFrame
    diagnostics::Vector{Float64}
    potential_energy::Vector{Float64}

    function TimeHistoryDiagnosticsSpin(
        particle_group::ParticleGroup,
        maxwell_solver::Maxwell1DFEM,
        kernel_smoother_0::ParticleMeshCoupling1D,
        kernel_smoother_1::ParticleMeshCoupling1D,
    )
        data = DataFrame(;
            Time=Float64[],
            KineticEnergy=Float64[],
            Kineticspin=Float64[],
            Momentum1=Float64[],
            Momentum2=Float64[],
            Momentum3=Float64[],
            Momentum4=Float64[],
            Momentum5=Float64[],
            Momentum6=Float64[],
            Momentum7=Float64[],
            Momentum8=Float64[],
            Momentum9=Float64[],
            Momentum10=Float64[],
            PotentialEnergyE1=Float64[],
            PotentialEnergyE2=Float64[],
            PotentialEnergyE3=Float64[],
            PotentialEnergyB2=Float64[],
            PotentialEnergyB3=Float64[],
            Transfer=Float64[],
            ErrorPoisson=Float64[],
        )

        diagnostics = zeros(Float64, 12)
        potential_energy = zeros(Float64, 5)

        return new(
            particle_group,
            maxwell_solver,
            kernel_smoother_0,
            kernel_smoother_1,
            data,
            diagnostics,
            potential_energy,
        )
    end
end

export write_step!

"""
    write_step!( thdiag, time, degree, efield_dofs, bfield_dofs,
                 efield_dofs_n, efield_poisson)

write diagnostics for PIC
- `time` : Time
- `afield_dofs[1]` : Magnetic Potential Ay
- `afield_dofs[2]` : Magnetic Potential Az
- `efield_dofs[1]` : Longitudinal Electric field Ex
- `efield_poisson` : Electric field compute from Poisson equation
- `efield_dofs[2]` : Ey
- `efield_dofs[3]` : Ez
- `degree` : Spline degree
"""
function write_step!(
    thdiag::TimeHistoryDiagnosticsSpin,
    time,
    degree,
    efield_dofs,
    afield_dofs,
    efield_dofs_n,
    efield_poisson,
    propagator,
)
    nn = thdiag.kernel_smoother_0.n_dofs
    tmp = zeros(Float64, nn)

    fill!(thdiag.diagnostics, 0.0)
    fill!(thdiag.potential_energy, 0.0)

    for i_part in 1:(thdiag.particle_group.n_particles)
        fill!(propagator.j_dofs[1], 0.0)
        fill!(propagator.j_dofs[2], 0.0)

        xi = thdiag.particle_group.array[1, i_part]
        vi = thdiag.particle_group.array[2, i_part]
        wi =
            thdiag.particle_group.array[3, i_part] *
            thdiag.particle_group.mass *
            thdiag.particle_group.common_weight
        s1 = thdiag.particle_group.array[4, i_part]
        s2 = thdiag.particle_group.array[5, i_part]
        s3 = thdiag.particle_group.array[6, i_part]
        # Kinetic energy
        v2 = evaluate(thdiag.kernel_smoother_0, xi[1], afield_dofs[1])
        v3 = evaluate(thdiag.kernel_smoother_0, xi[1], afield_dofs[2])
        thdiag.diagnostics[1] += 0.5 * (vi[1]^2 + v2[1]^2 + v3[1]^2) * wi[1]
        add_charge!(propagator.j_dofs[2], propagator.kernel_smoother_1, xi, 1.0)
        compute_rderivatives_from_basis!(
            propagator.j_dofs[1], propagator.maxwell_solver, propagator.j_dofs[2]
        )
        thdiag.diagnostics[2] +=
            propagator.HH * (
                afield_dofs[2]' * propagator.j_dofs[1] * wi * s2 -
                afield_dofs[1]' * propagator.j_dofs[1] * wi * s3
            )

        # \int (x f) dx dp ds
        thdiag.diagnostics[3] += xi[1] * wi[1]
        # \int (x s f) dx dp ds
        thdiag.diagnostics[4] += xi[1] * wi[1] * s1
        thdiag.diagnostics[5] += xi[1] * wi[1] * s2
        thdiag.diagnostics[6] += xi[1] * wi[1] * s3
        # \int (s f) dx dp ds
        thdiag.diagnostics[7] += wi[1] * s1
        thdiag.diagnostics[8] += wi[1] * s2
        thdiag.diagnostics[9] += wi[1] * s3
        # \int (s x B) dx dp ds
        thdiag.diagnostics[10] += (
            afield_dofs[1]' * propagator.j_dofs[1] * wi * s2 +
            afield_dofs[2]' * propagator.j_dofs[1] * wi * s3
        )
        thdiag.diagnostics[11] += afield_dofs[1]' * propagator.j_dofs[1] * wi * s1 * (-1.0)
        thdiag.diagnostics[12] += afield_dofs[2]' * propagator.j_dofs[1] * wi * s1 * (-1.0)
    end

    transfer = pic_diagnostics_transfer(
        thdiag.particle_group, thdiag.kernel_smoother_0, efield_dofs[1]
    )

    thdiag.potential_energy[1] =
        0.5 *
        inner_product(thdiag.maxwell_solver, efield_dofs[1], efield_dofs[1], degree - 1)
    thdiag.potential_energy[2] =
        0.5 * inner_product(thdiag.maxwell_solver, efield_dofs[2], efield_dofs[2], degree)
    thdiag.potential_energy[3] =
        0.5 * inner_product(thdiag.maxwell_solver, efield_dofs[3], efield_dofs[3], degree)

    compute_lderivatives_from_basis!(tmp, thdiag.maxwell_solver, afield_dofs[1])

    thdiag.potential_energy[4] =
        0.5 * l2norm_squared(thdiag.maxwell_solver, tmp, degree - 1)

    compute_lderivatives_from_basis!(tmp, thdiag.maxwell_solver, afield_dofs[2])

    thdiag.potential_energy[5] =
        0.5 * l2norm_squared(thdiag.maxwell_solver, tmp, degree - 1)

    return push!(
        thdiag.data,
        (
            time,
            thdiag.diagnostics...,
            thdiag.potential_energy...,
            transfer,
            maximum(abs.(efield_dofs[1] .- efield_poisson)),
        ),
    )
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "perl"
version = v"5.30.0"

# Collection of sources required to build perl
# with a few extra modules for polymake
sources = [
    "https://www.cpan.org/src/5.0/perl-$version.tar.gz" =>
    "851213c754d98ccff042caa40ba7a796b2cee88c5325f121be5cbb61bbf975f2",

    "https://cpan.metacpan.org/authors/id/I/IS/ISHIGAKI/JSON-4.01.tar.gz" =>
    "631593a939d4510e6ed76402556f38a34b20007237828670282e975712e0b1ed",

    "https://cpan.metacpan.org/authors/id/J/JO/JOSEPHW/XML-Writer-0.625.tar.gz" =>
    "e080522c6ce050397af482665f3965a93c5d16f5e81d93f6e2fe98084ed15fbe",

    "https://cpan.metacpan.org/authors/id/S/SH/SHLOMIF/XML-LibXML-2.0132.tar.gz" =>
    "721452e3103ca188f5968ab06d5ba29fe8e00e49f4767790882095050312d476",

    "https://cpan.metacpan.org/authors/id/S/SH/SHLOMIF/XML-LibXSLT-1.96.tar.gz" =>
    "2a5e374edaa2e9f9d26b432265bfea9b4bb7a94c9fbfef9047b298fce844d473",

    "https://cpan.metacpan.org/authors/id/J/JS/JSTOWE/TermReadKey-2.38.tar.gz" =>
    "5a645878dc570ac33661581fbb090ff24ebce17d43ea53fd22e105a856a47290",

    "https://cpan.metacpan.org/authors/id/H/HA/HAYASHI/Term-ReadLine-Gnu-1.36.tar.gz" =>
    "9a08f7a4013c9b865541c10dbba1210779eb9128b961250b746d26702bab6925",

    "https://cpan.metacpan.org/authors/id/G/GR/GRANTM/XML-SAX-1.02.tar.gz" =>
    "4506c387043aa6a77b455f00f57409f3720aa7e553495ab2535263b4ed1ea12a",

    "https://cpan.metacpan.org/authors/id/P/PE/PERIGRIN/XML-NamespaceSupport-1.12.tar.gz" =>
    "47e995859f8dd0413aa3f22d350c4a62da652e854267aa0586ae544ae2bae5ef",

    "https://cpan.metacpan.org/authors/id/G/GR/GRANTM/XML-SAX-Base-1.09.tar.gz" =>
    "66cb355ba4ef47c10ca738bd35999723644386ac853abbeb5132841f5e8a2ad0",

    "https://cpan.metacpan.org/authors/id/M/MA/MANWAR/SVG-2.84.tar.gz" =>
    "ec3d6ddde7a46fa507eaa616b94d217296fdc0d8fbf88741367a9821206f28af",

    "./bundled"
]

# Bash recipe for building
script = raw"""
perldir=`ls -1d perl-*`
cd $WORKSPACE/srcdir/
for dir in *;
do
   [[ "$dir" == "perl-"* ]] && continue;
   [[ "$dir" == "patches" ]] && continue;
   sed -i '1s/^/$ENV{PERL_CORE}=0;/' $dir/Makefile.PL
   mv $dir $perldir/cpan/${dir%-*};
done

cd $perldir/
if [[ $target != x86_64-linux* ]] && [[ $target != i686-linux* ]]; then
   # build host miniperl
   src=`pwd`
   mkdir $src/host
   cd $src/host
   $src/Configure -des -Dusedevel -Dmksymlinks -Dosname=linux -Dcc=$CC_FOR_BUILD -Dld=$LD_FOR_BUILD -Dlibs=-lm
   make -j${nproc} miniperl
   make -j${nproc} generate_uudmap
   cp -p miniperl $prefix/bin/miniperl-for-build
   cd ..
   cp ../patches/config-$target.sh config.sh
   atomic_patch -p1 ../patches/cross-nolibchecks.patch
   ./Configure -K -S
else
   atomic_patch -p1 ../patches/cross-nolibchecks.patch
   ./Configure -des -Dcc="$CC" -Dprefix=$prefix -Duseshrplib -Dsysroot=/opt/$target/$target/sys-root -Dccflags="-I${prefix}/include -I${prefix}/include/libxml2" -Dldflags="-L${prefix}/lib -Wl,-rpath,${prefix}/lib" -Dlddlflags="-shared -L${prefix}/lib -Wl,-rpath,${prefix}/lib"
fi

make -j${nproc} depend
make -j${nproc}

make install

pushd $prefix/lib
ln -s perl5/*/*/CORE/libperl.${dlext} libperl.${dlext}
popd

if [[ $target == *linux* ]]; then
patchelf --set-rpath $(patchelf --print-rpath ${prefix}/bin/perl | sed -e "s#${prefix}#\$ORIGIN/..#g") ${prefix}/bin/perl
for lib in ${prefix}/lib/perl5/*/*/auto/XML/LibXML/LibXML.${dlext} \
           ${prefix}/lib/perl5/*/*/auto/XML/LibXSLT/LibXSLT.${dlext};
do
   patchelf --set-rpath $(patchelf --print-rpath ${lib} | sed -e "s#${prefix}/lib#\$ORIGIN/../../../../../..#g") ${lib};
done
fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
    Linux(:i686, libc=:glibc)
    MacOS(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "perl", :perl)
    LibraryProduct(prefix, "libperl", :libperl)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/XML2-v2.9.9+0/build_XML2.v2.9.9.jl",
    "https://github.com/benlorenz/XSLTBuilder/releases/download/v1.1.33/build_XSLTBuilder.v1.1.33.jl",
    "https://github.com/benlorenz/readlineBuilder/releases/download/v8.0/build_readline.v8.0.0.jl",
    "https://github.com/benlorenz/ncursesBuilder/releases/download/v6.1/build_ncurses.v6.1.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# PA = LU factorization, Julia 0.6
# L, U, P <- mylup( A )
# L, U, p <- lu( A )

# I don't know how to print a pretty matrix like that in REPL.
function pp(A, name)
	println("--- $(name) ---")
	for i in 1:size(A, 1)
		for j in 1:size(A,2)
			print(A[i,j], "  ")
		end
		println("")
	end
	println("---  ---  ---")
end
# TODO: use writedlm(STDOUT, matrix)

function mylup(A, debug=true)
	# FIXME: NO ARGUMENT CHECK
	P = eye(size(A, 1))
	L = eye(size(A,1))
	U = copy(A)
	if debug
		pp(P, "P orig")
		pp(L, "L orig")
		pp(U, "U orig")
	end

	# select the best pivot for the first row
	cursor = indmax(abs.(U[1:end,1]))
	P[[1,cursor],:] = P[[cursor,1],:]
	U[[1,cursor],:] = U[[cursor,1],:]
	if U[1,1] == 0
		error("Pivot is 0 at step 1, the matrix is singular")
	end
	if debug
		println("> round $(1), indmax $(cursor), exchange $(1),$(cursor)")
		pp(P, "P exchange")
		pp(U, "U exhange")
	end

	# process row 2:end
	for i in 2:size(A,1)

		# select pivot for row i, i > 1
		cursor = indmax(abs.(U[i:end,i])) +i-1
		P[[i,cursor],:] = P[[cursor,i],:]
		U[[i,cursor],:] = U[[cursor,i],:]
		if U[i,i] == 0
			error("Pivot is 0 at step $(i), the matrix is singular")
		end
		if debug
			println("> round $(i), indmax $(cursor), exchange $(i),$(cursor)")
			pp(P, "P exchange")
			pp(U, "U exhange")
		end

		for j in 1:i-1
			multiplier = U[i,j]/U[j,j]
			U[i,:] -= multiplier .* U[j,:]
			L[i,j] = multiplier
			if debug
				println("> multiplier $(multiplier), row$(i) - $(multiplier) * row$(j)")
				pp(U, "U mult")
				pp(L, "L mult")
			end
		end
	end
	if debug
		pp(P*A-L*U, "P*A-L*U")
		println(maximum(P*A-L*U), " XXX")
	end
	return L, U, P
end

function compare(a, debug=false)
	L, U, P = mylup(a, debug)
	println("size $(size(a)), norm (LU-PA) = $(norm(L*U - P*a)), max element diff $(maximum(L*U-P*a))")
end

a = [ 1.0 1 1; 1 1 3; 2 5 8] # FIXME: weirdly fails if this is Array{Int64,2} with bloody InexactError
compare(a, false)

b = [ 1 -1 0 0 ; -1 2 -1 0 ; 0 -1 2 -1; 0 0 -1 2 ]
compare(b, false)

x = 20
c = rand(x,x)
compare(c, false)

for s in 1:99
	A = rand(s, s)
	compare(A)
end

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
module Config

using Logging
using Pkg

const config = Dict{String, Any}()

function load(path::String)
    if isfile(path)
        @info "*** Using config file: $path"
        merge!(config, Pkg.TOML.parsefile(path))
        return
    end
    @warn "config file $path not found."
end

function load()
    path = get(ENV, "TBB_CONFIG", nothing)
    if path == nothing
        env = get(ENV, "ENV", "priv")
        path = pwd() * "/config/app.$(env).toml"
        if get(ENV, "TRAVIS", "false") == "true"
            env = "travis"
            path = ((@__DIR__) |> dirname) * "/config/app.$(env).toml"
        end
    end
    load(path)
end

function get_config(key::String, default=nothing)
    length(config) == 0 && load()
    keys = split(key, ".")
    ret = config
    for kitem in keys
        ret = get(ret, kitem, Dict())
    end
    isa(ret, Dict) && length(ret) == 0 ? default : ret
end

function set_config(key::String, value)
    length(config) == 0 && load()
    keys = split(key, ".")
    node = config
    for kitem in keys[1:end-1]
        if haskey(node, kitem)
            node = node[kitem]
        else
            node[kitem] = Dict()
            node = node[kitem]
        end
    end
    node[keys[end]] = value
end

function log_level()
    log_level_str = lowercase(get_config("server.log_level", "debug"))

    (log_level_str == "debug") ? Logging.Debug :
    (log_level_str == "info")  ? Logging.Info  :
    (log_level_str == "warn")  ? Logging.Warn  : Logging.Error
end

end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/N2Oforcing.jl")

    add_comp!(m, n2oforcing)

    set_param!(m, :n2oforcing, :c_N2Oconcentration, readpagedata(m,"test/validationdata/$valdir/c_n2oconcentration.csv"))
    set_param!(m, :n2oforcing, :c_CH4concentration, readpagedata(m,"test/validationdata/$valdir/c_ch4concentration.csv"))

    ##running Model
    run(m)

    forcing=m[:n2oforcing,:f_N2Oforcing]
    forcing_compare=readpagedata(m,"test/validationdata/$valdir/f_n2oforcing.csv")

    @test forcing ≈ forcing_compare rtol=1e-3
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
@testset "Gtilde" begin
  Random.seed!(0)

  # NOTE: With Ni of 10000, should have around 10 it/s without loglike computation.
  # NOTE: With Ni of 1000, should have around > 100 it/s without loglike computation.
  NC = 1100
  NT = 1150
  QC = 100
  QT = 150
  yC = rand(SkewT(1, .6, 9, -5), NC - QC)
  yT = rand(SkewT(3, 1, 9, -10), NT - QT)
  K = 5
  onm_prior = OrderedNormalMeanPrior(K, Normal(0, 6), truncated(Normal(1, .3), 0, Inf))

  for arg in [(skew=true, tdist=true, mu=onm_prior), (skew=false, tdist=false, mu=nothing)]
    for beta in (nothing, 0, 1)

      model = beta === nothing ?
              Gtilde(yC, yT, K, 1, skew=arg.skew, tdist=arg.tdist, mu=arg.mu) :
              Gtilde(yC, yT, K, beta, skew=arg.skew, tdist=arg.tdist, mu=arg.mu)

      init = MCMC.make_init_state(model)
      spl = make_sampler(model, init=init)

      print_model_info(model); @test true

      nburn = beta === nothing ? 2 : 50
      nsamps = beta === nothing ? 2 : 100
      thin = 4
      callback = make_callback(model, nburn=nburn, nsamps=nsamps, thin=thin)
      chain, metrics = mcmc(spl, nsamps, init=init, nburn=nburn, thin=thin, callback=callback)
      @test length(metrics[:loglike]) == nsamps
      @assert length(unique(chain)) == nsamps
      cdd.printsummary(chain, metrics)

      # Hellinger
      if beta == 1 || beta == 0
        @test let
          pzero = Pzero(NC=NC, NT=NT, QC=QC, QT=QT)
          H = cdd.hellinger(chain, pzero, 5000)
          H_mean = mean(H)
          println("H: ", H_mean)
          0 < H_mean < 1
        end
      end

      # NOTE: This is just to check that the samples are updated every iteration.
      # Note that as discrete parameters (`lambda_i`) are present, the number of
      # unique samples may be < `nsamps`.
      if beta !== nothing
        syms = keys(init)
        for s in syms
          # println(s)
          samps = getindex.(chain, s)
          if s == :nu  # Metropolis
            if cdd.uses_tdist(model)
              @test length(unique(samps)) > 2
            else
              @test length(unique(samps)) == 1
            end
          elseif s in (:lambdaC, :lambdaT)  # discrete
            @test length(unique(samps)) > nsamps * 0.1
          elseif !cdd.uses_skew(model) && s in (:zetaC, :zetaT, :psi, :phi)
            @test length(unique(samps)) == 1
          elseif !cdd.uses_tdist(model) && s in (:vC, :vT)
            @test length(unique(samps)) == 1
          else
            @test length(unique(samps)) == nsamps
          end
        end
      end
    end
  end

  # TODO: Test that Ordered prior works.
end


text:
using Dates

struct Qualification
    Name::AbstractString
    Granted::Date
    Expires::Date
    CanExpire::Bool
    Qualification(Name::AbstractString,Granted::Date,Expires::Date,CanExpire::Bool) = new(Name,Granted,Expires,CanExpire)
    Qualification(Name::AbstractString,Granted::Date,Expires::Date) = Qualification(Name,Granted,Expires,true)
    Qualification(Name::AbstractString,Granted::Date) = Qualification(Name,Granted,Date(0),false)
    Qualification(Name::AbstractString) = Qualification(Name,Date(0),Date(0),false)
end

import Base.==
==(a::AbstractString,b::Qualification) = a==b.Name
==(a::Qualification,b::AbstractString) = a.Name==b

abstract type Person end
mutable struct Aircrew <: Person
    Name::AbstractString
    Qualifications::AbstractArray{Qualification}
    Aircrew(name::AbstractString,qualifications::AbstractArray{Qualification}) = new(name,qualifications)
    Aircrew(name::AbstractString) = Aircrew(name,Qualification[])
end

mutable struct NonAircrew <: Person
    Name::AbstractString
    NonAircrew(name::AbstractString) = new(name)
end


text:
"""
    componentgraph(V, copEV, bicon_comps)

Topological Gift Wrapping algorithm on 2D skeletons.

This is the offline part of the TGW algorithm. It takes in input a model and its
biconnected components mapping and evaluates usefull informations:
 1. Number of biconnected components.
 2. Component Graph of the biconnected structure.
 3. The 1-cells structure (UNMODIFIED).                                         # Could be removed?
 4. Association between non-dangling 2-cells and their orientation (∀ component)
 5. Association between 3-cells and 2-cells (with orientation, ∀ component).
 6. Association between 3-cells and their orientation (∀ component).
 7. Shell bounding boxes of the components.

See also: [`CAGD.planar_arrangement_2`](@ref).
It uses:
 - [`CAGD.get_external_cycle`](@ref)
 - [`CAGD.pre_containment_test`](@ref)
 - [`CAGD.prune_containment_graph`](@ref)
 - [`CAGD.transitive_reduction!`](@ref)
---
# Examples
```jldoctest
# Papillon
julia> V = [0.0 0.0; 0.0 3.0; 2.0 1.5; 4.0 0.0; 4.0 3.0];

julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
		[1 1 0 0 0] #1 -> 1,2
		[0 1 1 0 0] #2 -> 2,3
		[1 0 1 0 0] #3 -> 3,1
		[0 0 1 1 0] #4 -> 3,4
		[0 0 0 1 1] #5 -> 4,5
		[0 0 1 0 1] #6 -> 3,5
	]));

julia> copFE = SparseArrays.sparse(Array{Int8, 2}([
		[1 1 1 0 0 0] #1 -> 1,2,3
		[0 0 0 1 1 1] #2 -> 4,5,6
		[1 1 1 1 1 1] #3 -> 1,2,3,4,5,6 External
	]));

julia> bicon_comps = [[1, 2, 3], [4, 5, 6]];

julia> CAGD.componentgraph(V, copEV, bicon_comps)[1]
2

julia> CAGD.componentgraph(V, copEV, bicon_comps)[2]
2×2 SparseMatrixCSC{Int8,Int64} with 0 stored entries

julia> CAGD.componentgraph(V, copEV, bicon_comps)[4]
2-element Array{SparseMatrixCSC{Int8,Int64},1}:

  [1, 1]  =  -1
  [3, 1]  =  -1
  [1, 2]  =  1
  [2, 2]  =  -1
  [2, 3]  =  1
  [3, 3]  =  1

  [1, 3]  =  -1
  [3, 3]  =  -1
  [1, 4]  =  1
  [2, 4]  =  -1
  [2, 5]  =  1
  [3, 5]  =  1

julia> CAGD.componentgraph(V, copEV, bicon_comps)[5]
2-element Array{SparseMatrixCSC{Int8,Int64},1}:

  [1, 1]  =  -1
  [1, 2]  =  -1
  [1, 3]  =  1

  [1, 1]  =  1
  [1, 2]  =  1
  [1, 3]  =  -1

julia> CAGD.componentgraph(V, copEV, bicon_comps)[6]
2-element Array{SparseVector{Int8,Int64},1}:
   [1]  =  1
  [2]  =  1
  [3]  =  -1
   [1]  =  -1
  [2]  =  -1
  [3]  =  1

julia> CAGD.componentgraph(V, copEV, bicon_comps)[7]
2-element Array{Any,1}:
 ([0.0 0.0], [2.0 3.0])
 ([2.0 0.0], [4.0 3.0])
```
"""
function componentgraph(V, copEV, bicon_comps)
    # arrangement of isolated components
	n = size(bicon_comps, 1)
   	shells = Array{Lar.Chain, 1}(undef, n)
	boundaries = Array{Lar.ChainOp, 1}(undef, n)
	EVs = Array{Lar.ChainOp, 1}(undef, n)
    # for each component
	for p in 1 : n
		ev = copEV[sort(bicon_comps[p]), :]
        # computation of 2-cells
		fe = Lar.Arrangement.minimal_2cycles(V, ev)
        # exterior cycle
		shell_num = CAGD.get_external_cycle(V, ev, fe)
        # decompose each fe (co-boundary local to component)
		EVs[p] = ev
		tokeep = setdiff(1:fe.m, shell_num)
		boundaries[p] = fe[tokeep, :]
		shells[p] = fe[shell_num, :]
    end
    # computation of bounding boxes of isolated components
	shell_bboxes = []
	for i in 1 : n
    	vs_indexes = (abs.(EVs[i]')*abs.(shells[i])).nzind
   		push!(shell_bboxes, Lar.bbox(V[vs_indexes, :]))
	end
    # computation and reduction of containment graph
	containment_graph = CAGD.pre_containment_test(shell_bboxes)
	containment_graph = CAGD.prune_containment_graph(V,EVs,shells,containment_graph)
	CAGD.transitive_reduction!(containment_graph)
	return n, containment_graph, V, EVs, boundaries, shells, shell_bboxes
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 2 - COMPONENT GRAPH (TGW) METHODS
#    get_external_cycle
#    pre_containment_test
#    prune_containment_graph
#    transitive_reduction!
#-------------------------------------------------------------------------------

"""
    get_external_cycle(model::CAGD.Model)::Union{Int64, Nothing}
    get_external_cycle(V::Lar.Points, EV::Lar.ChainOp, FE::Lar.ChainOp)

Evaluates the index of the external 3D-cell of a Topological Model.

This method looks for and retrieve the external cycle of the `model`.
If the cell does not exist then it returns `nothing`.
The method is also callable by the signature `(V, EV, FE)`.

See also: [`CAGD.componentgraph`](@ref).
---
# Examples
```jldoctest
# ``K^4`` graph
julia> model = CAGD.Model(hcat([[0.0,0.0], [4.0,0.0], [2.0,3.0], [2.0,1.5]]...));
julia> CAGD.addModelCells!(model, 1, [[1,2],[2,3],[3,1],[1,4],[2,4],[3,4]]);
julia> CAGD.addModelCells!(model, 2, [[1,4,5],[2,5,6],[3,4,6],[1,2,3]]);

julia> CAGD.get_external_cycle(model)
4

julia> CAGD.deleteModelCell!(model, 2, 4);
julia> typeof(CAGD.get_external_cycle(model))
Nothing
```
"""
function get_external_cycle(model::CAGD.Model)::Union{Int64, Nothing}
    CAGD.get_external_cycle(
        convert(Lar.Points, model.G'),
        model.T[1],
        model.T[2]
    )
end

function get_external_cycle(
        V::Lar.Points,
        EV::Lar.ChainOp,
        FE::Lar.ChainOp
    )::Union{Int64, Nothing}

    FV = abs.(FE)*EV

    vs = sparsevec(mapslices(sum, abs.(EV), dims=1)').nzind
    minv_x1 = maxv_x1 = minv_x2 = maxv_x2 = pop!(vs)
    for i in vs
        if V[i, 1] > V[maxv_x1, 1]
            maxv_x1 = i
        elseif V[i, 1] < V[minv_x1, 1]
            minv_x1 = i
        end
        if V[i, 2] > V[maxv_x2, 2]
            maxv_x2 = i
        elseif V[i, 2] < V[minv_x2, 2]
            minv_x2 = i
        end
    end
    cells = intersect(
        FV[:, minv_x1].nzind,
        FV[:, maxv_x1].nzind,
        FV[:, minv_x2].nzind,
        FV[:, maxv_x2].nzind
    )
    if length(cells) == 1
        return cells[1]
    else
        for c in cells
            if Lar.face_area(V, EV, FE[c, :]) < 0
                return c
            end
        end
    end
end

"""
    pre_containment_test(
        bboxes::Array{Tuple{Array{Float64,2},Array{Float64,2}},1}
    )::SparseMatrixCSC{Int8,Int64}

Generate the containment graph associated to `bboxes`.

The function evaluate a `SparseArrays{Int8, n, n}` PRE-containmnet graph meaning
it stores `1` when the bbox associated to the `i`-th component contains entirely
the bbox associated to the `j`-th one. (where `n = length(bboxes)`)

See also: [`CAGD.componentgraph`](@ref).

---

# Examples
```jldoctest
# Planar Tiles
julia> bboxes = [
		([0.0 0.0], [1.0 1.0])
		([0.0 0.0], [0.4 0.4])
		([0.6 0.0], [1.0 0.4])
		([0.0 0.6], [0.4 1.0])
		([0.6 0.6], [1.0 1.0])
		([2.0 3.0], [2.0 3.0]) # no intersection
	];

julia> CAGD.pre_containment_test(bboxes)
6x6 SparseMatrixCSC{Int8,Int64} with 4 stored entries:
  [2, 1]  =  1
  [3, 1]  =  1
  [4, 1]  =  1
  [5, 1]  =  1
```
"""
function pre_containment_test(
        bboxes#::Array{Tuple{Array{Float64,2},Array{Float64,2}},1}              #
    )::SparseMatrixCSC{Int8,Int64}
    n = length(bboxes)
    containment_graph = spzeros(Int8, n, n)

    for i in 1:n
        for j in 1:n
            if i != j && Lar.bbox_contains(bboxes[j], bboxes[i])
                containment_graph[i, j] = 1
            end
        end
    end

    return containment_graph
end

"""
    prune_containment_graph(
            V::Lar.Points,
            EVs::Array{Lar.ChainOp, 1},
            shells::Array{Lar.Chain,1},
            graph::SparseMatrixCSC{Int8, Int64}
        )::SparseMatrixCSC{Int8, Int64}

Prunes the containment `graph` from the non-included faces.

This method prunes the containment graph eliminating the included bouning boxes
that do not corresponds to included faces.

Do note that this method expects to work on biconnectet clusters of faces.

See also: [`CAGD.componentgraph`](@ref).
```
"""
function prune_containment_graph(
        V::Lar.Points,
        EVs::Array{Lar.ChainOp, 1},
        shells::Array{Lar.Chain,1},
        graph::SparseMatrixCSC{Int8, Int64}
    )::SparseMatrixCSC{Int8, Int64}

    graph = deepcopy(graph) # copy needed
    n = length(EVs)
    length(shells) == n || throw(ArgumentError("EVs and shells not coherent"))
    size(graph) == (n, n) || throw(ArgumentError("graph dim not coherent"))

    for i in 1 : n
        # Take a point of the i-th shell by taking a vertex of its first edge
        an_edge = shells[i].nzind[1]
        origin_index = EVs[i][an_edge, :].nzind[1]
        origin = V[origin_index, :]

        # Check the containment for every other cell by testing if such a point
        #   is inner to the face defined byt that cell
        for j in 1 : n
            if i != j && graph[i, j] == 1
                shell_edge_indexes = shells[j].nzind
                ev = EVs[j][shell_edge_indexes, :]

                if !CAGD.point_in_face(origin, V, ev)
                    graph[i, j] = 0
                end
            end
         end

     end
     return dropzeros(graph)
end

"""
    transitive_reduction!(graph)

Prune the redundancies in containment `graph`.

The ``m``-layer containment `graph` is reduced to the corresponding 1-layer
containment `graph` where only the first level nested components are considered.

See also: [`CAGD.componentgraph`](@ref).

---

# Examples
```jldoctest
julia> containment = [
		0 0 0 0 0 0
		0 0 0 0 0 0
		1 0 0 0 0 0
		1 0 0 0 0 0
		1 0 1 0 0 0
		1 0 1 0 1 0
	];

julia> CAGD.transitive_reduction!(containment)

julia> containment
6×6 Array{Int64,2}:
 0  0  0  0  0  0
 0  0  0  0  0  0
 1  0  0  0  0  0
 1  0  0  0  0  0
 0  0  1  0  0  0
 0  0  0  0  1  0
```
"""
function transitive_reduction!(graph)
    n = size(graph, 1)
    for j in 1:n
        for i in 1:n
            if graph[i, j] > 0
                for k in 1:n
                    if graph[j, k] > 0
                        graph[i, k] = 0
                    end
                end
            end
        end
    end
end