text:
"""
    body2hor(xb, yb, zb, psi, theta, phi)

Transform body coordinates to local horizon.

# Arguments
* `xb, yb, zb`: body coordinates.
* `psi, theta, phi`: Euler angles. Yaw, pitch, roll (rad).
"""
function body2hor(xb, yb, zb, psi, theta, phi)

    s_psi, c_psi = sin(psi), cos(psi)
    s_theta, c_theta = sin(theta), cos(theta)
    s_phi, c_phi = sin(phi), cos(phi)

    xh = (c_theta * c_psi)                         * xb +
         (s_phi * s_theta * c_psi - c_phi * s_psi) * yb +
         (c_phi * s_theta * c_psi + s_phi * s_psi) * zb

    yh = (c_theta * s_psi)                         * xb +
         (s_phi * s_theta * s_psi + c_phi * c_psi) * yb +
         (c_phi * s_theta * s_psi - s_phi * c_psi) * zb

    zh = -s_theta         * xb +
          s_phi * c_theta * yb +
          c_phi * c_theta * zb

    return [xh, yh, zh]
end

"""
    rot_matrix_body2hor(psi, theta, phi)

Rotation matrix from body to local horizon

# Arguments
* `psi, theta, phi`: Euler angles. Yaw, pitch, roll (rad).
"""
function rot_matrix_body2hor(psi, theta, phi)
    sψ, cψ = sin(psi), cos(psi)
    sθ, cθ = sin(theta), cos(theta)
    sϕ, cϕ = sin(phi), cos(phi)

    r = [ cθ*cψ   (sϕ*sθ*cψ - cϕ*sψ)   (cϕ*sθ*cψ + sϕ*sψ);
          cθ*sψ   (sϕ*sθ*sψ + cϕ*cψ)   (cϕ*sθ*sψ - sϕ*cψ);
         -sθ          sϕ*cθ                cϕ*cθ         ]
    return r
end

"""
    body2hor(xb, yb, zb, q0, q1, q2, q3)

Transform body coordinates to local horizon.

# Arguments
* `xb, yb, zb`: body coordinates.
* `q0, q1, q2, q3`: quaternions.
"""
function body2hor(xb, yb, zb, q0, q1, q2, q3)
    q02, q12, q22, q32 = q0*q0, q1*q1, q2*q2, q3*q3

    xh = (q02+q12-q22-q32) * xb + 2*(q1*q2 - q0*q3) * yb + 2*(q1*q3 + q0*q2) * zb
    yh = 2*(q1*q2 + q0*q3) * xb + (q02-q12+q22-q32) * yb + 2*(q2*q3 - q0*q1) * zb
    zh = 2*(q1*q3 - q0*q2) * xb + 2*(q2*q3 + q0*q1) * yb + (q02-q12-q22+q32) * zb
    return [xh, yh, zh]
end

"""
    rot_matrix_body2hor(q0, q1, q2, q3)

Rotation matrix to transform from body to local horizon.

# Arguments
* `q0, q1, q2, q3`: quaternions.

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 41, formula 1.4-23)
"""
function rot_matrix_body2hor(q0, q1, q2, q3)

    q02, q12, q22, q32 = q0*q0, q1*q1, q2*q2, q3*q3

    r = [(q02+q12-q22-q32)    2*(q1*q2 - q0*q3)    2*(q1*q3 + q0*q2);
         2*(q1*q2 + q0*q3)    (q02-q12+q22-q32)    2*(q2*q3 - q0*q1);
         2*(q1*q3 - q0*q2)    2*(q2*q3 + q0*q1)    (q02-q12-q22+q32)]

    return r
end

"""
    hor2body(xb, yb, zb, q0, q1, q2, q3)

Transform local horizon corrdinates to body coordinates.

# Arguments
* `xh, yh, zh`: local horizon coordinates.
* `q0, q1, q2, q3`: quaternions.
"""
function hor2body(xh, yh, zh, q0, q1, q2, q3)

    q02, q12, q22, q32 = q0*q0, q1*q1, q2*q2, q3*q3

    xb = (q02+q12-q22-q32) * xh + 2*(q1*q2 + q0*q3) * yh + 2*(q1*q3 - q0*q2) * zh
    yb = 2*(q1*q2 - q0*q3) * xh + (q02-q12+q22-q32) * yh + 2*(q2*q3 + q0*q1) * zh
    zb = 2*(q1*q3 + q0*q2) * xh + 2*(q2*q3 - q0*q1) * yh + (q02-q12-q22+q32) * zh
    return [xb, yb, zb]
end

"""
    rot_matrix_hor2body(q0, q1, q2, q3)

Rotation matrix to transform from local horizon to body.

# Arguments
* `q0, q1, q2, q3`: quaternions.

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 41, formula 1.4-23)
"""
function rot_matrix_hor2body(q0, q1, q2, q3)
    r = rot_matrix_body2hor(q0, q1, q2, q3)
    return r'
end

"""
    hor2body(xh, yh, zh, psi, theta, phi)

Transform local horizon coordinates to body.

# Arguments
* `xh, yh, zh`: local horizon coordinates.
* `psi, theta, phi`: Euler angles. Yaw, pitch, roll (rad).
"""
function hor2body(xh, yh, zh, psi, theta, phi)

    s_psi, c_psi = sin(psi), cos(psi)
    s_theta, c_theta = sin(theta), cos(theta)
    s_phi, c_phi = sin(phi), cos(phi)

    xb =  c_theta * c_psi * xh +
          c_theta * s_psi * yh -
          s_theta         * zh

    yb = (s_phi * s_theta * c_psi - c_phi * s_psi) * xh +
         (s_phi * s_theta * s_psi + c_phi * c_psi) * yh +
         (s_phi * c_theta)                         * zh

    zb = (c_phi * s_theta * c_psi + s_phi * s_psi) * xh +
         (c_phi * s_theta * s_psi - s_phi * c_psi) * yh +
         (c_phi * c_theta)                         * zh

    return [xb, yb, zb]
end

"""
    rot_matrix_hor2body(psi, theta, phi)

Rotation matrix from body to local horizon

# Arguments
* `psi, theta, phi`: Euler angles. Yaw, pitch, roll (rad).
"""
function rot_matrix_hor2body(psi, theta, phi)
    r = rot_matrix_body2hor(psi, theta, phi)
    return r'
end

"""
    wind2hor(xw, yw, zw, chi, gamma, mu)

Transform wind coordinates to local horizon.

# Arguments
* `xw, yw, zw`: wind coordinates.
* `chi, gamma, mu`: Velocity angles. Yaw (azimuth), pitch(elevation), roll (rad).
"""
function wind2hor(xw, yw, zw, chi, gamma, mu)

    s_chi, c_chi = sin(chi), cos(chi)
    s_gamma, c_gamma = sin(gamma), cos(gamma)
    s_mu, c_mu = sin(mu), cos(mu)

    xh = (c_gamma * c_chi)                       * xw +
         (s_mu * s_gamma * c_chi - c_mu * s_chi) * yw +
         (c_mu * s_gamma * c_chi + s_mu * s_chi) * zw

    yh = (c_gamma * s_chi)                       * xw +
         (s_mu * s_gamma * s_chi + c_mu * c_chi) * yw +
         (c_mu * s_gamma * s_chi - s_mu * c_chi) * zw

    zh = -s_gamma        * xw +
          s_mu * c_gamma * yw +
          c_mu * c_gamma * zw

    return [xh, yh, zh]
end

"""
    hor2wind(xh, yh, zh, chi, gamma, mu)

Transform local horizon coordinates to wind.

# Arguments
* `xh, yh, zh`: local horizon coordinates.
* `chi, gamma, mu`: Velocity angles. Yaw (azimuth), pitch(elevation), roll (rad).
"""
function hor2wind(xh, yh, zh, chi, gamma, mu)

    s_chi, c_chi = sin(chi), cos(chi)
    s_gamma, c_gamma = sin(gamma), cos(gamma)
    s_mu, c_mu = sin(mu), cos(mu)

    xw =  c_gamma * c_chi * xh +
          c_gamma * s_chi * yh -
          s_gamma         * zh

    yw = (s_mu * s_gamma * c_chi - c_chi * s_chi) * xh +
         (s_mu * s_gamma * s_chi + c_mu  * c_chi) * yh +
         (s_mu * c_gamma)                         * zh

    zw = (c_mu * s_gamma * c_chi + s_mu * s_chi) * xh +
         (c_mu * s_gamma * s_chi - s_mu * c_chi) * yh +
         (c_mu * c_gamma)                        * zh

    return [xw, yw, zw]
end

"""
    body2wind(xb, yb, zb, alpha, beta)

Transform body coordinates to wind.

# Arguments
* `xb, yb, zb`: body coordinates.
* `alpha, beta`: Aerodynamic angles. Angle of attack, angle of side-slip (rad).
"""
function body2wind(xb, yb, zb, alpha, beta)

    s_alpha, c_alpha = sin(alpha), cos(alpha)
    s_beta,  c_beta  = sin(beta),  cos(beta)

    xw =  c_alpha*c_beta * xb + s_beta * yb + s_alpha*c_beta * zb
    yw = -c_alpha*s_beta * xb + c_beta * yb - s_alpha*s_beta * zb
    zw = -s_alpha        * xb +                      c_alpha * zb

    return [xw, yw, zw]
end

"""
    wind2body(xw, yw, zw, alpha, beta)

Transform body coordinates to wind.

# Arguments
* `xw, yw, zw`: wind coordinates.
* `alpha, beta`: Aerodynamic angles. Angle of attack, angle of side-slip (rad).
"""
function wind2body(xw, yw, zw, alpha, beta)

    s_alpha, c_alpha = sin(alpha), cos(alpha)
    s_beta,  c_beta  = sin(beta),  cos(beta)

    xb = c_alpha*c_beta * xw - c_alpha*s_beta * yw - s_alpha * zw
    yb = s_beta         * xw + c_beta         * yw
    zb = s_alpha*c_beta * xw - s_alpha*s_beta * yw + c_alpha * zw

    return [xb, yb, zb]
end

"""
    ecef2hor(xecef, yecef, zecef, lat, lon)

Transform ECEF (Earth Fixed Earth Centered) coordinates to local horizon
coordinates using geodetic latitude and longitude.

# Arguments
 * `xecef, yecef, zecef`: ECEF coordinates.
 * `lat`: geodetic latitude (rad).
 * `lon`: longitude (rad).

 # References

- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 36, formula 1.4-9)

 # Notes

 X, Y, Z ECEF are defined in a different way in [1]. In order to reproduce this
 transformation substitute in [1] Y->X, Z->Y, X->Z

 [1] Only defines the opposite transformation. However this must be the
 transpose of hor2ecef.
"""
function ecef2hor(xecef, yecef, zecef, lat, lon)
    clat, slat = cos(lat), sin(lat)
    clon, slon = cos(lon), sin(lon)

    xh = -slat*clon * xecef - slat*slon * yecef + clat * zecef
    yh = -slon      * xecef + clon      * yecef
    zh = -clat*clon * xecef - clat*slon * yecef - slat * zecef

    return [xh, yh, zh]
end

"""
    rot_matrix_ecef2hor(lat, lon)

Rotation matrix to transform from ECEF to local horizon.

# Arguments
* `lat, lon`: geodetic latitude and longitude (rad)

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 36, formula 1.4-9)

# Notes

 X, Y, Z ECEF are defined in a different way in [1]. In order to reproduce this
 transformation substitute in [1] Y->X, Z->Y, X->Z
"""
function rot_matrix_ecef2hor(lat, lon)
    r = rot_matrix_hor2ecef(lat, lon)
    return r'
end

"""
    hor2ecef(xh, yh, zh, lat, lon)

Transform local horizon coordinates to ECEF (Earth Centered Earth Fixed)
coordinates using geodetic latitude and longitude.

# Arguments
* `xh, yh, zh`: local horizon coordinates.
* `lat`: geodetic latitude (rad).
* `lon`: longitude (rad).

# References

- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 36, formula 1.4-9)

# Notes

 X, Y, Z ECEF are defined in a different way in [1]. In order to reproduce this
 transformation substitute in [1] Y->X, Z->Y, X->Z
"""
function hor2ecef(xh, yh, zh, lat, lon)
    clat, slat = cos(lat), sin(lat)
    clon, slon = cos(lon), sin(lon)

    xecef = -slat*clon * xh - slon * yh - clat*clon * zh
    yecef = -slat*slon * xh + clon * yh - clat*slon * zh
    zecef =  clat      * xh +           - slat      * zh

    return [xecef, yecef, zecef]
end

"""
    rot_matrix_hor2ecef(lat, lon)

Rotation matrix to transform from local horizon to ECEF.

# Arguments
* `lat, lon`: geodetic latitude and longitude (rad)

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 36, formula 1.4-9)

# Notes

 X, Y, Z ECEF are defined in a different way in [1]. In order to reproduce this
 transformation substitute in [1] Y->X, Z->Y, X->Z
"""
function rot_matrix_hor2ecef(lat, lon)
    clat, slat = cos(lat), sin(lat)
    clon, slon = cos(lon), sin(lon)

    r = [-slat*clon    -slon    -clat*clon;
         -slat*slon     clon    -clat*slon;
          clat          0.0     -slat     ]
    return r
end

"""
    body2ecef(xb, yb, zb, lat, lon, psi, theta, phi)

Transform body coordinates to ECEF coordinates.

# Arguments
* `xb, yb, zb`: body coordinates.
* `lat`: geodetic latitude (rad).
* `lon`: longitude (rad).
* `psi, theta, phi`: Euler angles. Yaw, pitch, roll (rad).
"""
function body2ecef(xb, yb, zb, lat, lon, psi, theta, phi)
    xh, yh, zh = body2hor(xb, yb, zb, psi, theta, phi)
    xecef, yecef, zecef = hor2ecef(xh, yh, zh, lat, lon)
    return [xecef, yecef, zecef]
end

"""
    rot_matrix_body2ecef(lat, lon, psi, theta, phi)

Rotation matrix to transform from body to ECEF

# Arguments
* `lat`: geodetic latitude (rad).
* `lon`: longitude (rad).
* `psi, theta, phi`: Euler angles. Yaw, pitch, roll (rad).

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 37, formula 1.4-11)
"""
function rot_matrix_body2ecef(lat, lon, psi, theta, phi)
    r1 = rot_matrix_body2hor(psi, theta, phi)
    r2 = rot_matrix_hor2ecef(lat, lon)
    return r2 * r1
end

"""
    ecef2body(xecef, yecef, zecef, lat, lon, psi, theta, phi)

Transform ECEF coordinates to body coordinates.

# Arguments
* `xecef, yecef, zecef`: ECEF (Earth Centered Earth Fixed) coordinates.
* `lat`: geodetic latitude (rad).
* `lon`: longitude (rad).
* `psi, theta, phi`: Euler angles. Yaw, pitch, roll (rad).
"""
function ecef2body(xecef, yecef, zecef, lat, lon, psi, theta, phi)
    xh, yh, zh = ecef2hor(xecef, yecef, zecef, lat, lon)
    xb, yb, zb = hor2body(xh, yh, zh, psi, theta, phi)
    return [xb, yb, zb]
end

"""
    rot_matrix_ecef2body(lat, lon, psi, theta, phi)

Rotation matrix to transform from ECEF to body

# Arguments
* `lat`: geodetic latitude (rad).
* `lon`: longitude (rad).
* `psi, theta, phi`: Euler angles. Yaw, pitch, roll (rad).

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 37, formula 1.4-11)
"""
function rot_matrix_ecef2body(lat, lon, psi, theta, phi)
    r1 = rot_matrix_hor2body(psi, theta, phi)
    r2 = rot_matrix_ecef2hor(lat, lon)
    return r1 * r2
end

"""
    body2ecef(xb, yb, zb, lat, lon, q0, q1, q2, q3)

Transform body coordinates to ECEF coordinates.

# Arguments
* `xb, yb, zb`: body coordinates.
* `lat`: geodetic latitude (rad).
* `lon`: longitude (rad).
* `q0, q1, q2, q3`: quaternions.
"""
function body2ecef(xb, yb, zb, lat, lon, q0, q1, q2, q3)
    xh, yh, zh = body2hor(xb, yb, zb, q0, q1, q2, q3)
    xecef, yecef, zecef = hor2ecef(xh, yh, zh, lat, lon)
    return [xecef, yecef, zecef]
end

"""
    rot_matrix_body2ecef(lat, lon, q0, q1, q2, q3)

Rotation matrix to transform from body to ECEF

# Arguments
* `lat`: geodetic latitude (rad).
* `lon`: longitude (rad).
* `q0, q1, q2, q3`: quaternions.

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 37, formula 1.4-11)
"""
function rot_matrix_body2ecef(lat, lon, q0, q1, q2, q3)
    r1 = rot_matrix_body2hor(q0, q1, q2, q3)
    r2 = rot_matrix_hor2ecef(lat, lon)
    return r2 * r1
end

"""
    ecef2body(xecef, yecef, zecef, lat, lon, q0, q1, q2, q3)

Transform ECEF coordinates to body coordinates.

# Arguments
* `xecef, yecef, zecef`: ECEF (Earth Centered Earth Fixed) coordinates.
* `lat`: geodetic latitude (rad).
* `lon`: longitude (rad).
* `q0, q1, q2, q3`: quaternions.
"""
function ecef2body(xecef, yecef, zecef, lat, lon, q0, q1, q2, q3)
    xh, yh, zh = ecef2hor(xecef, yecef, zecef, lat, lon)
    xb, yb, zb = hor2body(xh, yh, zh, q0, q1, q2, q3)
    return [xb, yb, zb]
end

"""
    rot_matrix_ecef2body(lat, lon, q0, q1, q2, q3)

Rotation matrix to transform from ECEF to body

# Arguments
* `lat`: geodetic latitude (rad).
* `lon`: longitude (rad).
* `q0, q1, q2, q3`: quaternions.

# References
- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 37, formula 1.4-11)
"""
function rot_matrix_ecef2body(lat, lon, q0, q1, q2, q3)
    r1 = rot_matrix_hor2body(q0, q1, q2, q3)
    r2 = rot_matrix_ecef2hor(lat, lon)
    return r1 * r2
end

"""
    euler2quaternion(psi, theta, phi)

Transform Euler angles orientation to quaternion orientation.

# Arguments
* `psi, theta, phi`: Euler angles. Yaw, pitch, roll (rad).

# References
- [1] Zipfel, P. H. (2007). Modeling and simulation of aerospace vehicle
 dynamics. American Institute of Aeronautics and Astronautics.
 (page 126, formula 4.78)
"""
function euler2quaternion(psi, theta, phi)

    s_psi2, c_psi2 = sin(psi/2.0), cos(psi/2.0)
    s_theta2, c_theta2 = sin(theta/2.0), cos(theta/2.0)
    s_phi2, c_phi2 = sin(phi/2.0), cos(phi/2.0)

    q0 = c_psi2*c_theta2*c_phi2 + s_psi2*s_theta2*s_phi2
    q1 = c_psi2*c_theta2*s_phi2 - s_psi2*s_theta2*c_phi2
    q2 = c_psi2*s_theta2*c_phi2 + s_psi2*c_theta2*s_phi2
    q3 = s_psi2*c_theta2*c_phi2 - c_psi2*s_theta2*s_phi2

    return [q0, q1, q2, q3]
end

"""
    quaternion2euler(q0, q1, q2, q3)

Transform quaternion orientation to Euler angles orientation.

# Arguments
* `q0, q1, q2, q3`: quaternions.

# References
- [1] Zipfel, P. H. (2007). Modeling and simulation of aerospace vehicle
 dynamics. American Institute of Aeronautics and Astronautics.
 (page 127, formula 4.82)
"""
function quaternion2euler(q0, q1, q2, q3)
    psi = atan(2 * (q1*q2 + q0*q3), q0*q0 + q1*q1 - q2*q2 - q3*q3)
    theta = asin(-2 * (q1*q3 - q0*q2))
    phi = atan(2 * (q2*q3 + q0*q1), q0*q0 - q1*q1 - q2*q2 + q3*q3)

    return [psi, theta, phi]
end

"""
    llh2ecef(lat, lon, height; ellipsoid=WGS84)

Transform geodetic latitude, longitude and ellipsoidal height to ECEF for the
given ellipsoid (default ellipsoid is WGS84)

# References
- [1] Rogers, R. M. (2007). Applied mathematics in integrated navigation
 systems. American Institute of Aeronautics and Astronautics.
 (Page 75, equations 4.20, 4.21, 4.22)
"""
function llh2ecef(lat, lon, height; ellipsoid=WGS84)
    f = ellipsoid.f
    a = ellipsoid.a
    e2 = ellipsoid.e2

    var = a / sqrt(1.0 - e2 * sin(lat)*sin(lat))

    px = (var + height) * cos(lat)*cos(lon)
    py = (var + height) * cos(lat)*sin(lon)
    pz = (var * (1.0 - e2) + height)* sin(lat)

    return [px, py, pz]
end

"""
    ecef2llh(xecef, yecef, zecef; ellipsoid=WGS84)

Transform ECEF coordinates to geodetic latitude, longitude and ellipsoidal
height for the given ellipsoid (default ellipsoid is WGS84)

# References
- [1] Bowring, B. R. (1976). Transformation from spatial to geographical
 coordinates. Survey review, 23(181), 323-327.
- [2] Bowring, B. R. (1985). The accuracy of geodetic latitude and height
 equations. Survey Review, 28(218), 202-206.

# Notes

* The transformation is direct without iterations as [1] introduced the need to
 iterate for near Earth positions.
* [2] is an updated of increased accuracy of [1]. The former is used in this
 implementation although the latter implementation is commented in the code.
* Model becomes unstable if latitude is close to 90º. An alternative equation
 can be found in [2] equation (16) but has not been implemented.
"""
function ecef2llh(xecef, yecef, zecef; ellipsoid=WGS84)

    x, y, z = xecef, yecef, zecef
    e = sqrt(ellipsoid.e2);
    e2 = ellipsoid.e2
    ϵ2 = ellipsoid.ϵ2
    a = ellipsoid.a
    b = ellipsoid.b

    p = sqrt(x*x + y*y)
    R = sqrt(p*p + z*z)
    θ = atan(z, p)

    # [1] equation (1) does not change in [2]
    lon = atan(y, x)

    # u -> geographical latitude
    # [1] below equation (4) and [2] equation (6)
    # This lead to errors in latitud with maximum value 0.0018"
    #u = atan(a/b * z/x)

    # [2] equation (17) If the latitude is also required to be very accurate
    # for outer-space positions then the value of u for (6) should be obtained
    # from:
    u = atan( b*z / (a*p) * (1 + ϵ2 * b / R) )

    # [1] equation (4)
    #lat = atan((z + ϵ2 * b * sin(u)^3) / (x - e2 * a * cos(u)^3))

    # [2] equation (18)
    lat = atan( (z + ϵ2 * b * sin(u)^3) / (p - e2 * a * cos(u)^3) )

    # [2] after equation (1)
    v = a / sqrt(1.0 - e2*sin(lat)^2)
    # [2] equation (7)
    # height = p*cos(lat) + z*sin(lat) - a*a / v
    # equivalent to [2] equation (8)
    height = R * cos(lat - θ) - a*a / v
    # which is insensitive to the error in latitude calculation (The influence
    # is of order 2 [2] equation (12))
    return [lat, lon, height]
end


text:
module js

import GR

pxwidth = 640
pxheight = 480

@static if VERSION < v"0.7.0-DEV.4762"
    macro cfunction(f, rt, tup)
        :(Base.cfunction($(esc(f)), $(esc(rt)), Tuple{$(esc(tup))...}))
    end
end

id_count = 0
js_running = false

mutable struct JSTermWidget
    identifier::String
    width::Int
    height::Int
    visible::Bool
end

function inject_js()
  global comm
  comm = nothing
  _js_fallback = "https://gr-framework.org/downloads/gr-latest.js"
  _gr_js = if isfile(joinpath(ENV["GRDIR"], "lib", "gr.js"))
    _gr_js = try
      _gr_js = open(joinpath(ENV["GRDIR"], "lib", "gr.js")) do f
        _gr_js = read(f, String)
        _gr_js
      end
    catch e
      nothing
    end
    _gr_js
  end
  _jsterm = """
    if (typeof grJSTermRunning === 'undefined' || !grJstermReady) {
      BOXZOOM_THRESHOLD = 5;
      BOXZOOM_TRIGGER_THRESHHOLD = 1000;
      MAX_KERNEL_CONNECTION_ATTEMPTS = 25;
      KERNEL_CONNECT_WAIT_TIME = 100;
      REFRESH_PLOT_TIMEOUT = 100;

      var comm = undefined;
      var idcount = 0;
      var widgets = [];
      var jupyterRunning = false;

      function saveLoad(url, callback, maxtime) {
        let script = document.createElement('script');
        script.onload = function() {
          callback();
        };
        script.onerror = function() {
          console.error(url + ' can not be loaded.');
        };
        script.src = url;
        document.head.appendChild(script);
        setTimeout(function() {
          if (!grJstermReady) {
            console.error(url + ' can not be loaded.');
          }
        }, maxtime);
      }

      function sendEvt(data) {
        if (jupyterRunning) {
          comm.send({
            "type": "evt",
            "content": data
          });
        }
      }

      function jsLoaded() {
        grJstermReady = true;
        for (let or in onready) {
          or();
          onready = [];
        }
      }

      function canvasRemoved(id) {
        if (jupyterRunning) {
          comm.send({
            "type": "removed",
            "content": id
          });
        }
      }

      function saveData(data, id) {
        if (jupyterRunning) {
          comm.send({
            "type": "save",
            "content": {
              "id": id,
              "data": JSON.stringify(data)
            }
          });
        }
      }

      function registerComm(kernel) {
        kernel.comm_manager.register_target('jsterm_comm', function(c) {
          c.on_msg(function(msg) {
            if(typeof msg.content.data.type !== 'undefined' && msg.content.data.type == 'evt') {
              if (typeof widgets[msg.content.data.json.plot] !== 'undefined') {
                widgets[msg.content.data.json.plot].msgHandleEvent(msg.content.data.json);
              }
            } else {
              draw(msg);
            }
          });
          c.on_close(function() {});
          window.addEventListener('beforeunload', function(e) {
            c.close();
          });
          comm = c;
        });
      }

      function onLoad() {
        if (typeof Jupyter !== 'undefined') {
          jupyterRunning = true;
          initKernel(1);
        } else {
          drawSavedData();
        }
      };

      function initKernel(attempt) {
        let kernel = Jupyter.notebook.kernel;
        if (typeof kernel === 'undefined' || kernel == null) {
          if (attempt < MAX_KERNEL_CONNECTION_ATTEMPTS) {
            setTimeout(function() {
              initKernel(attempt + 1);
            }, KERNEL_CONNECT_WAIT_TIME);
          }
        } else {
          registerComm(kernel);
          Jupyter.notebook.events.on('kernel_ready.Kernel', function() {
            registerComm(kernel);
            for (let key in widgets) {
              widgets[key].init();
            }
          });
          drawSavedData();
        }
      }

      function draw(msg) {
        if (!grJstermReady) {
          onready.push(function() {
            return draw(msg);
          });
        } else if (!GR.is_ready) {
          GR.ready(function() {
            return draw(msg);
          });
        } else {
          if (typeof widgets[msg.content.data.canvasid] === 'undefined') {
            widgets[msg.content.data.canvasid] = new JSTermWidget(idcount, msg.content.data.canvasid);
            idcount += 1;
          }
          widgets[msg.content.data.canvasid].draw(msg);
        }
      }

      function drawSavedData() {
        let data = document.getElementsByClassName("jsterm-data");
        for (let i = 0; i < data.length; i++) {
          let msg = data[i].innerText;
          draw(JSON.parse(msg));
        }
      }

      if (document.readyState!='loading') {
        onLoad();
      } else if (document.addEventListener) {
        document.addEventListener('DOMContentLoaded', onLoad);
      } else document.attachEvent('onreadystatechange', function() {
        if (document.readyState=='complete') {
          onLoad();
        }
      });

      function JSTermWidget(id, htmlId) {

        this.init = function() {
          this.canvas = undefined;
          this.overlayCanvas = undefined;
          this.args = undefined;
          this.id = id;  // context id for meta.c (switchmeta)
          this.gr = undefined;
          this.htmlId = htmlId;

          this.waiting = false;
          this.oncanvas = function() {};

          // event handling
          this.pinching = false;
          this.panning = false;
          this.prevMousePos = undefined;
          this.boxzoom = false;
          this.keepAspectRatio = true;
          this.boxzoomTriggerTimeout = undefined;
          this.boxzoomPoint = [0, 0];
        }
        this.init();

        this.getCoords = function(event) {
          let rect = this.canvas.getBoundingClientRect();
          //TODO mind the canvas-padding if necessary!
          return [Math.floor(event.clientX - rect.left), Math.floor(event.clientY - rect.top)];
        };

        this.grEventinput = function(mouseargs) {
          this.gr.switchmeta(this.id);
          this.gr.inputmeta(mouseargs);
          this.gr.current_canvas = this.canvas;
          this.gr.current_context = this.gr.current_canvas.getContext('2d');
          this.gr.select_canvas();
          this.gr.plotmeta();
        };

        this.handleWheel = function(x, y, angle_delta) {
          if (typeof this.boxzoomTriggerTimeout !== 'undefined') {
            clearTimeout(this.boxzoomTriggerTimeout);
          }
          let mouseargs = this.gr.newmeta();
          this.gr.meta_args_push(mouseargs, "x", "i", [x]);
          this.gr.meta_args_push(mouseargs, "y", "i", [y]);
          this.gr.meta_args_push(mouseargs, "angle_delta", "d", [angle_delta]);
          this.grEventinput(mouseargs);
        };

        this.mouseHandleWheel = function (event) {
          let coords = this.getCoords(event);
          sendEvt({
            "x": coords[0],
            "y": coords[1],
            "angle_delta": event.deltaY,
            "type": "mousewheel",
            "plot": this.htmlId  // HTML canvas` string id
          });
          this.handleWheel(coords[0], coords[1], event.deltaY);
          event.preventDefault();
        };

        this.handleMousedown = function(x, y, button, ctrlKey) {
          if (typeof this.boxzoomTriggerTimeout !== 'undefined') {
            clearTimeout(this.boxzoomTriggerTimeout);
          }
          if (button == 0) {
            this.overlayCanvas.style.cursor = 'move';
            this.panning = true;
            this.boxzoom = false;
            this.prevMousePos = [x, y];
            this.boxzoomTriggerTimeout = setTimeout(function() {this.startBoxzoom(x, y, ctrlKey);}.bind(this), BOXZOOM_TRIGGER_THRESHHOLD);
          } else if (button == 2) {
            this.startBoxzoom(x, y, ctrlKey);
          }
        };

        this.mouseHandleMousedown = function (event) {
          let coords = this.getCoords(event);
          sendEvt({
            "x": coords[0],
            "y": coords[1],
            "button": event.button,
            "ctrlKey": event.ctrlKey,
            "type": "mousedown",
            "plot": this.htmlId  // HTML canvas` string id
          });
          this.handleMousedown(coords[0], coords[1], event.button, event.ctrlKey);
          event.preventDefault();
        };

        this.startBoxzoom = function(x, y, ctrlKey) {
          this.panning = false;
          this.boxzoom = true;
          if (ctrlKey) {
            this.keepAspectRatio = false;
          }
          this.boxzoomPoint = [x, y];
          this.overlayCanvas.style.cursor = 'nwse-resize';
        };

        this.handleMouseup = function(x, y, button) {
          if (typeof this.boxzoomTriggerTimeout !== 'undefined') {
            clearTimeout(this.boxzoomTriggerTimeout);
          }
          if (this.boxzoom) {
            if ((Math.abs(this.boxzoomPoint[0] - x) >= BOXZOOM_THRESHOLD) && (Math.abs(this.boxzoomPoint[1] - y) >= BOXZOOM_THRESHOLD)) {
              let mouseargs = this.gr.newmeta();
              if (this.boxzoomPoint[0] < x) {
                this.gr.meta_args_push(mouseargs, "left", "i", [this.boxzoomPoint[0]]);
                this.gr.meta_args_push(mouseargs, "right", "i", [x]);
              } else {
                this.gr.meta_args_push(mouseargs, "right", "i", [this.boxzoomPoint[0]]);
                this.gr.meta_args_push(mouseargs, "left", "i", [x]);
              }
              if (this.boxzoomPoint[1] < y) {
                this.gr.meta_args_push(mouseargs, "top", "i", [this.boxzoomPoint[1]]);
                this.gr.meta_args_push(mouseargs, "bottom", "i", [y]);
              } else {
                this.gr.meta_args_push(mouseargs, "bottom", "i", [this.boxzoomPoint[1]]);
                this.gr.meta_args_push(mouseargs, "top", "i", [y]);
              }
              if (this.keepAspectRatio) {
                this.gr.meta_args_push(mouseargs, "keepAspectRatio", "i", [1]);
              } else {
                this.gr.meta_args_push(mouseargs, "keepAspectRatio", "i", [0]);
              }
              this.grEventinput(mouseargs);
            }
          }
          this.prevMousePos = undefined;
          this.overlayCanvas.style.cursor = 'auto';
          this.panning = false;
          this.boxzoom = false;
          this.keepAspectRatio = true;
          let context = this.overlayCanvas.getContext('2d');
          context.clearRect(0, 0, this.overlayCanvas.width, this.overlayCanvas.height);
        };

        this.mouseHandleMouseup = function (event) {
          let coords = this.getCoords(event);
          sendEvt({
            "x": coords[0],
            "y": coords[1],
            "button": event.button,
            "type": "mouseup",
            "plot": this.htmlId  // HTML canvas` string id
          });
          this.handleMouseup(coords[0], coords[1], event.button);
          event.preventDefault();
        };

        this.handleLeave = function() {
          if (typeof this.boxzoomTriggerTimeout !== 'undefined') {
            clearTimeout(this.boxzoomTriggerTimeout);
          }
          this.overlayCanvas.style.cursor = 'auto';
          this.panning = false;
          this.prevMousePos = undefined;
          if (this.boxzoom) {
            let context = this.overlayCanvas.getContext('2d');
            context.clearRect(0, 0, this.overlayCanvas.width, this.overlayCanvas.height);
          }
          this.boxzoom = false;
          this.keepAspectRatio = true;
        };

        this.mouseHandleMouseleave = function(event) {
          this.handleLeave();
        };

        this.handleMousemove = function(x, y) {
          if (this.panning) {
            if (typeof this.boxzoomTriggerTimeout !== 'undefined') {
              clearTimeout(this.boxzoomTriggerTimeout);
            }
            let mouseargs = this.gr.newmeta();
            this.gr.meta_args_push(mouseargs, "x", "i", [this.prevMousePos[0]]);
            this.gr.meta_args_push(mouseargs, "y", "i", [this.prevMousePos[1]]);
            this.gr.meta_args_push(mouseargs, "xshift", "i", [x - this.prevMousePos[0]]);
            this.gr.meta_args_push(mouseargs, "yshift", "i", [y - this.prevMousePos[1]]);
            this.grEventinput(mouseargs);
            this.prevMousePos = [x, y];
          } else if (this.boxzoom) {
            let context = this.overlayCanvas.getContext('2d');
            let diff = [x - this.boxzoomPoint[0], y - this.boxzoomPoint[1]];
            if (this.keepAspectRatio) {
              if (Math.abs(diff[0]) / this.overlayCanvas.width > Math.abs(diff[1]) / this.overlayCanvas.height) {
                let factor = Math.abs(x - this.boxzoomPoint[0]) / this.overlayCanvas.width;
                diff[1] = Math.sign(diff[1]) * factor * this.overlayCanvas.height;
              } else {
                let factor = Math.abs(y - this.boxzoomPoint[1]) / this.overlayCanvas.height;
                diff[0] = Math.sign(diff[0]) * factor * this.overlayCanvas.width;
              }
            }
            context.clearRect(0, 0, this.overlayCanvas.width, this.overlayCanvas.height);
            if (diff[0] * diff[1] >= 0) {
              this.overlayCanvas.style.cursor = 'nwse-resize';
            } else {
              this.overlayCanvas.style.cursor = 'nesw-resize';
            }
            context.fillStyle = '#FFAAAA';
            context.strokeStyle = '#FF0000';
            context.beginPath();
            context.rect(this.boxzoomPoint[0], this.boxzoomPoint[1], diff[0], diff[1]);
            context.globalAlpha = 0.2;
            context.fill();
            context.globalAlpha = 1.0;
            context.stroke();
            context.closePath();
          }
        };

        this.mouseHandleMousemove = function (event) {
          let coords = this.getCoords(event);
          sendEvt({
            "x": coords[0],
            "y": coords[1],
            "type": "mousemove",
            "plot": this.htmlId  // HTML canvas` string id
          });
          this.handleMousemove(coords[0], coords[1]);
          event.preventDefault();
        };

        this.handleDoubleclick = function(x, y) {
          let mouseargs = this.gr.newmeta();
          this.gr.meta_args_push(mouseargs, "x", "i", [x]);
          this.gr.meta_args_push(mouseargs, "y", "i", [y]);
          this.gr.meta_args_push(mouseargs, "key", "s", "r");
          this.grEventinput(mouseargs);
        };

        this.mouseHandleDoubleclick = function(event) {
          let coords = this.getCoords(event);
          sendEvt({
            "x": coords[0],
            "y": coords[1],
            "type": "doubleclick",
            "plot": this.htmlId  // HTML canvas` string id
          });
          this.handleDoubleclick(coords[0], coords[1]);
          event.preventDefault();
        };

        this.msgHandleEvent = function(msg) {
          switch (msg.type) {
            case "mousewheel":
              this.handleWheel(msg.x, msg.y, msg.angle_delta);
              break;
            case "mouseup":
              this.handleMouseup(msg.x, msg.y, msg.button);
              break;
            case "mousedown":
              this.handleMousedown(msg.x, msg.y, msg.button, msg.ctrlKey);
              break;
            case "mousemove":
              this.handleMousemove(msg.x, msg.y);
              break;
            case "doubleclick":
              this.handleDoubleclick(msg.x, msg.y);
              break;
            default:
              break;
          }
        };

        this.draw = function(msg) {
          if (this.waiting) {
            this.oncanvas = function() {
              return this.draw(msg);
            };
          } else {
            if (document.getElementById('jsterm-' + msg.content.data.canvasid) == null) {
              canvasRemoved(msg.content.data.canvasid);
              this.canvas = undefined;
              this.waiting = true;
              this.oncanvas = function() {
                return draw(msg);
              };
              setTimeout(function() {
                this.refreshPlot(msg, 0);
              }.bind(this), REFRESH_PLOT_TIMEOUT);
            } else {
              if (document.getElementById('jsterm-data-' + this.htmlId) == null) {
                saveData(msg, msg.content.data.canvasid);
              }
              if (typeof this.canvas === 'undefined' || typeof this.overlayCanvas === 'undefined') {
                this.canvas = document.getElementById('jsterm-' + this.htmlId);
                this.overlayCanvas = document.getElementById('jsterm-overlay-' + this.htmlId);
                this.overlayCanvas.addEventListener('DOMNodeRemoved', function() {
                  canvasRemoved(msg.content.data.canvasid);
                  this.canvas = undefined;
                  this.waiting = true;
                  this.oncanvas = function() {};
                });
                this.overlayCanvas.style.cursor = 'auto';

                //registering event handler
                this.overlayCanvas.addEventListener('wheel', function(evt) { this.mouseHandleWheel(evt); }.bind(this));
                this.overlayCanvas.addEventListener('mousedown', function(evt) { this.mouseHandleMousedown(evt); }.bind(this));
                this.overlayCanvas.addEventListener('mousemove', function(evt) { this.mouseHandleMousemove(evt); }.bind(this));
                this.overlayCanvas.addEventListener('mouseup', function(evt) { this.mouseHandleMouseup(evt); }.bind(this));
                this.overlayCanvas.addEventListener('mouseleave', function(evt) { this.mouseHandleMouseleave(evt); }.bind(this));
                this.overlayCanvas.addEventListener('dblclick', function(evt) { this.mouseHandleDoubleclick(evt); }.bind(this));
                this.overlayCanvas.addEventListener('contextmenu', function(event) {
                  event.preventDefault();
                  return false;
                });
              }
              if (typeof this.gr === 'undefined') {
                this.gr = new GR('jsterm-' + this.htmlId);
              }
              if (typeof this.args === 'undefined') {
                this.args = this.gr.newmeta();
              }
              this.waiting = false;
              this.gr.switchmeta(this.id);
              this.gr.current_canvas = this.canvas; //TODO is this always set? (check)
              this.gr.current_context = this.gr.current_canvas.getContext('2d');
              this.gr.select_canvas();
              this.gr.meta_args_push(this.args, "size", "dd", [this.canvas.width, this.canvas.height]);
              this.gr.readmeta(this.args, msg.content.data.json);
              this.gr.plotmeta(this.args);
            }
          }
        };

        this.refreshPlot = function(msg, count) {
          if (document.getElementById('jsterm-' + this.htmlId) == null) {
            setTimeout(function() {
              this.refreshPlot(msg, count + 1);
            }.bind(this), REFRESH_PLOT_TIMEOUT);
          } else {
            this.waiting = false;
            if (typeof this.oncanvas !== 'undefined') {
              this.oncanvas();
            }
          }
        };
      }
    }
    var grJSTermRunning = true;"""
  if _gr_js === nothing
      _gr_js = string("""
        JSTerm.saveLoad('""", _js_fallback, """', jsLoaded, 10000);
        var grJstermReady = false;
      """)
  else
    _gr_js = string(_gr_js, "var grJstermReady = true;")
  end
  display(HTML(string("""
    <script type="text/javascript">
      """, _gr_js, """
      """, _jsterm, """
    </script>
  """)))
end

function JSTermWidget(name::String, id::Int64)
  global id_count, js_running
  if GR.isijulia()
    id_count += 1
    if !js_running
      inject_js()
      js_running = true
    end
    JSTermWidget(string(name, id), pxwidth, pxheight, false)
  else
    error("JSTermWidget is only available in IJulia environments")
  end
end

function jsterm_display(widget::JSTermWidget)
  global pxwidth, pxheight
  if GR.isijulia()
    display(HTML(string("<div style=\"position: relative; width: ", widget.width, "px; height: ", widget.height, "px;\"><canvas id=\"jsterm-overlay-", widget.identifier, "\" style=\"position:absolute; top: 0; right: 0; z-index: 1;\" width=\"", widget.width, "\" height=\"", widget.height, "\"></canvas>
        <canvas id=\"jsterm-", widget.identifier, "\" style=\"position: absolute; top: 0; right: 0; z-index: 0;\"width=\"", widget.width, "\" height=\"", widget.height, "\"></canvas>")))
    widget.visible = true
  else
    error("jsterm_display is only available in IJulia environments")
  end
end

comm = nothing

evthandler = function(data)
  nothing
end

function register_evthandler(f)
  global evthandler
  evthandler = f
end

function send_evt(f)
  global comm
  if comm === nothing
    error("JSTerm comm not initialized.")
  else
    Main.IJulia.send_comm(comm, Dict("json" => f, "type" => "evt"))
  end
end

function jsterm_send(widget::JSTermWidget, data::String)
  global js_running, draw_condition, comm, PXWIDTH, PXHEIGHT, counter
  if GR.isijulia()
    if comm === nothing
      comm = Main.IJulia.Comm("jsterm_comm")
      comm.on_close = function comm_close_callback(msg)
        global js_running
        js_running = false
      end
      comm.on_msg = function comm_msg_callback(msg)
        if haskey(msg.content["data"], "type")
          if msg.content["data"]["type"] == "removed"
            jswidgets[msg.content["data"]["content"]].visible = false
            jsterm_display(jswidgets[msg.content["data"]["content"]])
          elseif msg.content["data"]["type"] == "save"
            display(HTML(string("<div style=\"display:none;\" id=\"jsterm-data-", msg.content["data"]["content"]["id"], "\" class=\"jsterm-data\">", msg.content["data"]["content"]["data"], "</div>")))
          elseif msg.content["data"]["type"] == "evt"
            evthandler(msg.content["data"]["content"])
          end
        end
      end
    end
    Main.IJulia.send_comm(comm, Dict("json" => data, "canvasid" => widget.identifier))
  else
    error("jsterm_send is only available in IJulia environments")
  end
end

function recv(name::Cstring, id::Int64, msg::Cstring)
    # receives string from C and sends it to JS via Comm
    global js_running
    if !js_running
      inject_js()
      js_running = true
    end
    _name = unsafe_string(name)
    if haskey(jswidgets, string(_name, id))
        widget = jswidgets[string(_name, id)]
    else
        widget = JSTermWidget(_name, id)
        jswidgets[string(_name, id)] = widget
    end
    if !widget.visible
      jsterm_display(widget)
    end
    jsterm_send(widget, unsafe_string(msg))
    return convert(Int32, 1)
end

function send(name::Cstring, id::Int64)
    # Dummy function, not in use
    return convert(Cstring, "String")
end

jswidgets = nothing
send_c = nothing
recv_c = nothing

function initjs()
    global jswidgets, send_c, recv_c

    jswidgets = Dict{String, JSTermWidget}()
    send_c = @cfunction(send, Cstring, (Cstring, Int64))
    recv_c = @cfunction(recv, Int32, (Cstring, Int64, Cstring))
    send_c, recv_c
end

end # module

ext: jl


text:
module Registries

using ..FastAI
using ..FastAI.Datasets
using ..FastAI.Datasets: DatasetLoader, DataDepLoader, isavailable, loaddata, typify

import Markdown
using DataDeps
using FeatureRegistries
using FeatureRegistries: Registry, Field
using InlineTest

_formatblock(t::Type{<:Tuple}) = _formatblock(Tuple(t.types))
_formatblock(t::Tuple) = map(_formatblock, t)
_formatblock(T::Type) = T
_formatblock(::T) where T = T

function blocktypesmatch(
        BSupported::Type,
        BWanted::Type)
    # true if both types are part of the same type tree
    return BSupported <: BWanted || BWanted <: BSupported
end
function blocktypesmatch(B1::Type{<:Tuple}, B2::Type{<:Tuple})
    all(blocktypesmatch(b1, b2) for (b1, b2) in zip(B1.types, B2.types))
end

blocktypesmatch(BSupported::Type, ::Type{Any}) = true
blocktypesmatch(BSupported::Type{Any}, ::Type) = true
blocktypesmatch(BSupported::Type{Any}, ::Type{Any}) = true
blocktypesmatch(bsupported, bwanted) = blocktypesmatch(typify(bsupported), typify(bwanted))

@testset "`blocktypesmatch`" begin
    @test blocktypesmatch(FastAI.Image, FastAI.Image{2})
    @test blocktypesmatch(FastAI.Image, FastAI.Image)
    @test blocktypesmatch(FastAI.Image{2}, FastAI.Image)
    @test blocktypesmatch(Tuple{FastAI.Image}, Tuple{FastAI.Image})
    @test blocktypesmatch(Tuple{FastAI.Image{2}}, Tuple{FastAI.Image})
    @test blocktypesmatch(Tuple{FastAI.Image{2}}, Any)
    @test blocktypesmatch(FastAI.Image{2}(), FastAI.Image{2})
    @test blocktypesmatch(FastAI.Image, FastAI.Image{2}())
    @test blocktypesmatch((FastAI.Image{2}(), FastAI.Label(1:10)), (FastAI.Image, FastAI.Label))
    @test blocktypesmatch((FastAI.Image{2}(), FastAI.AbstractBlock), (FastAI.Image, FastAI.Label))
    @test !blocktypesmatch(FastAI.Image{2}(), Label(1:10))
end

include("datasets.jl")
include("tasks.jl")
include("recipes.jl")

export datasets,
    learningtasks,
    datarecipes,
    find,
    info,
    load

function __init__()
    _registerdatasets(DATASETS)
end

end


text:
module TableView

using Tables
using WebIO, JSExpr, JSON, Dates, UUIDs
using Observables: @map

export showtable

const ag_grid_imports = []

function __init__()
    version = readchomp(joinpath(@__DIR__, "..", "ag-grid.version"))
    empty!(ag_grid_imports)
    for f in ["ag-grid.js", "ag-grid.css", "ag-grid-light.css", "ag-grid-dark.css"]
        push!(ag_grid_imports, normpath(joinpath(@__DIR__, "..", "deps", "ag-grid-$(version)", f)))
    end
    pushfirst!(ag_grid_imports, normpath(joinpath(@__DIR__, "rowNumberRenderer.js")))
end

to_css_size(s::AbstractString) = s
to_css_size(s::Real) = "$(s)px"

struct IteratorAndFirst{F, T}
    first::F
    source::T
    len::Int
    function IteratorAndFirst(x)
        len = Base.haslength(x) ? length(x) : 0
        first = iterate(x)
        return new{typeof(first), typeof(x)}(first, x, len)
    end
    function IteratorAndFirst(first, x)
        len = Base.haslength(x) ? length(x) + 1 : 1
        return new{typeof(first), typeof(x)}(first, x, len)
    end
end
Base.IteratorSize(::Type{IteratorAndFirst{F, T}}) where {F, T} = Base.IteratorSize(T)
Base.length(x::IteratorAndFirst) = x.len
Base.IteratorEltype(::Type{IteratorAndFirst{F, T}}) where {F, T} = Base.IteratorEltype(T)
Base.eltype(x::IteratorAndFirst) = eltype(x.source)
Base.iterate(x::IteratorAndFirst) = x.first
function Base.iterate(x::IteratorAndFirst, st)
    st === nothing && return nothing
    return iterate(x.source, st)
end

showtable(table::AbstractMatrix; kwargs...) = showtable(Tables.table(table); kwargs...)

"""
    showtable(table; options::Dict{Symbol, Any} = Dict{Symbol, Any}(), option_mutator! = identity,
              dark = false, height = :auto, width = "100%", cell_changed = nothing)

Return a `WebIO.Scope` that displays the provided `table`.

Optional arguments:
    - `options`: Directly passed to agGrid's `Grid` constructor. Refer to the
                 [documentation](https://www.ag-grid.com/documentation/) for more info.
    - `options_mutator!`: Runs on the `options` dictionary populated by TableView and allows for
                          customizing the grid (at your own risk -- you can break the package by
                          supplying invalid options).
    - `dark`: Switch to a dark theme.
    - `title`: Displayed above the table if non-empty;
    - `height`/`width`: CSS attributes specifying the output height and with.
    - `cell_changed`: Either `nothing` or a function that takes a single argument with the fields
                      `"new"`, `"old"`, `"row"`, and `"col"`. This function is called whenever the
                      user edits a table field. Note that all values will be strings, so you need to
                      do the necessary conversions yourself.
"""
function showtable(table;
        options::Dict{Symbol, Any} = Dict{Symbol, Any}(),
        option_mutator! = nothing,
        dark::Bool = false,
        title::String = "",
        height = :auto,
        width = "100%",
        cell_changed = nothing
    )
    rows = Tables.rows(table)
    it_sz = Base.IteratorSize(rows)
    has_len = it_sz isa Base.HasLength || it_sz isa Base.HasShape
    tablelength = has_len ? length(rows) : nothing

    if height === :auto
        height = 500
        if tablelength !== nothing
            # header + footer height ≈ 40px, 28px per row
            height = min(50 + tablelength*28, height)
        end
    end

    schema = Tables.schema(rows)
    if schema === nothing
        st = iterate(rows)
        rows = IteratorAndFirst(st, rows)
        names = Symbol[]
        types = []
        if st !== nothing
            row = st[1]
            for nm in propertynames(row)
                push!(names, nm)
                push!(types, typeof(getproperty(row, nm)))
            end
            schema = Tables.Schema(names, types)
        else
            # no schema and no rows
        end
    else
        names = schema.names
        types = schema.types
    end

    async = tablelength === nothing || tablelength > 10_000

    w = Scope(imports = ag_grid_imports)

    onCellValueChanged = @js function () end
    if cell_changed != nothing
        onedit = Observable(w, "onedit", Dict{Any,Any}(
                                                "row" => 0,
                                                "col" => 0,
                                                "old" => 0,
                                                "new" => 0
                                              ))
        on(onedit) do x
            cell_changed(x)
        end

        onCellValueChanged = @js function (ev)
            $onedit[] = Dict(
                                "row" => ev.rowIndex,
                                "col" => ev.colDef.headerName,
                                "old" => ev.oldValue,
                                "new" => ev.newValue
                            )
        end
    end

    coldefs = [Dict(
                :headerName => string(n),
                :editable => cell_changed !== nothing,
                :headerTooltip => string(types[i]),
                :field => string(n),
                :sortable => !async,
                :resizable => true,
                :type => types[i] <: Union{Missing, T where T <: Number} ? "numericColumn" : nothing,
                :filter => async ? false : types[i] <: Union{Missing, T where T <: Dates.Date} ? "agDateColumnFilter" :
                         types[i] <: Union{Missing, T where T <: Number} ? "agNumberColumnFilter" : true,
               ) for (i, n) in enumerate(names)]

    pushfirst!(coldefs, Dict(
        :headerName => "Row",
        :editable => false,
        :headerTooltip => "",
        :field => "__row__",
        :sortable => !async,
        :resizable => true,
        :type => "numericColumn",
        :cellRenderer => "rowNumberRenderer",
        :filter => false
    ))

    options[:onCellValueChanged] = onCellValueChanged
    options[:columnDefs] = coldefs
    options[:multiSortKey] = "ctrl"
    options[:rowSelection] = "multiple"

    for e in ["onCellClicked", "onCellDoubleClicked", "onRowClicked", "onCellFocused", "onCellKeyDown"]
        o = Observable{Any}(w, e, nothing)
        handler = @js function (ev)
            @var x = Dict()
            if ev.rowIndex !== undefined
                x["rowIndex"] = ev.rowIndex + 1
            end
            if ev.colDef !== undefined
                x["column"] = ev.colDef.headerName
            end
            $o[] = x
        end
        options[Symbol(e)] = handler
    end

    id = string("grid-", string(uuid1())[1:8])
    w.dom = dom"div"(
                dom"div"(
                    title,
                    style = Dict(
                        "background-color" => dark ? "#1c1f20" : "#F5F7F7",
                        "color" => dark ? "#F5F7F7" : "#1c1f20",
                        "height" => isempty(title) ? "0" : "18px",
                        "padding" => isempty(title) ? "0" : "5px",
                        "font-family" => """-apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, Oxygen-Sans, Ubuntu, Cantarell, "Helvetica Neue", sans-serif"""
                    )
                ),
                dom"div"(className = "ag-theme-balham$(dark ? "-dark" : "")",
                     style = Dict("width" => to_css_size(width),
                                  "height" => to_css_size(height)),
                     id = id
                )
            )

    showfun = async ? _showtable_async! : _showtable_sync!

    showfun(w, schema, types, rows, tablelength, id, options, option_mutator!)

    return w
end

function _showtable_sync!(w, schema, types, rows, tablelength, id, options, option_mutator!)
    options[:rowData] = JSONText(table2json(schema, rows, types))

    (option_mutator! !== nothing) && option_mutator!(options, w)

    license = get(ENV, "AG_GRID_LICENSE_KEY", nothing)
    handler = @js function (RowNumberRenderer, agGrid)
        @var gridOptions = $options
        @var el = document.getElementById($id)
        gridOptions.components = Dict(
            "rowNumberRenderer" => RowNumberRenderer
        )
        if $(license !== nothing)
            agGrid.LicenseManager.setLicenseKey($license)
        end
        this.table = @new agGrid.Grid(el, gridOptions)
        gridOptions.columnApi.autoSizeAllColumns()
    end
    onimport(w, handler)
end

function _showtable_async!(w, schema, types, rows, tablelength, id, options, option_mutator!)
    rowparams = Observable(w, "rowparams", Dict("startRow" => 1,
                                                "endRow" => 100,
                                                "successCallback" => @js v -> nothing))
    requestedrows = Observable(w, "requestedrows", JSONText("{}"))
    on(rowparams) do x
        requestedrows[] = JSONText(table2json(schema, rows, types, requested = [x["startRow"] + 1, x["endRow"] + 1]))
    end

    onjs(requestedrows, @js function (val)
        ($rowparams[]).successCallback(val, $(tablelength))
    end)

    options[:maxConcurrentDatasourceRequests] = 1
    options[:cacheBlockSize] = 1000
    options[:maxBlocksInCache] = 100
    options[:rowModelType] = "infinite"
    options[:datasource] = Dict(
        "getRows" =>
            @js function (rowParams)
                $rowparams[] = rowParams
            end
        ,
        "rowCount" => tablelength
    )
    license = get(ENV, "AG_GRID_LICENSE_KEY", nothing)

    (option_mutator! !== nothing) && option_mutator!(options, w)

    handler = @js function (RowNumberRenderer, agGrid)
        @var gridOptions = $options
        @var el = document.getElementById($id)

        gridOptions.components = Dict(
            "rowNumberRenderer" => RowNumberRenderer
        )
        if $(license !== nothing)
            agGrid.LicenseManager.setLicenseKey($license)
        end
        this.table = @new agGrid.Grid(el, gridOptions)
        gridOptions.columnApi.autoSizeAllColumns()
    end
    onimport(w, handler)
end

function _is_javascript_safe(x::Integer)
    min_safe_int = -(Int64(2)^53-1)
    max_safe_int = Int64(2)^53-1
    min_safe_int < x < max_safe_int
end

function _is_javascript_safe(x::AbstractFloat)
    min_safe_float = -(Float64(2)^53-1)
    max_safe_float = Float64(2)^53-1
    min_safe_float < x < max_safe_float
end

# directly write JSON instead of allocating temporary dicts etc
function table2json(schema, rows, types; requested = nothing)
    io = IOBuffer()
    rowwriter = JSON.Writer.CompactContext(io)
    JSON.begin_array(rowwriter)
    ser = JSON.StandardSerialization()
    for (i, row) in enumerate(rows)
        if requested != nothing && (i < first(requested) || i > last(requested))
            continue
        end
        JSON.delimit(rowwriter)
        columnwriter = JSON.Writer.CompactContext(io)
        JSON.begin_object(columnwriter)
        Tables.eachcolumn(schema, row) do val, ind, name
            if val isa Real && isfinite(val) && _is_javascript_safe(val)
                JSON.show_pair(columnwriter, ser, name, val)
            elseif val === nothing || val === missing
                JSON.show_pair(columnwriter, ser, name, repr(val))
            else
                JSON.show_pair(columnwriter, ser, name, sprint(print, val))
            end
        end
        JSON.end_object(columnwriter)
    end
    JSON.end_array(rowwriter)
    String(take!(io))
end

include("TableViewAdvanced.jl")
export TableViewAdvanced

end


text:
# sensemap-sim.jl

using MIRT: ir_mri_sensemap_sim
using MIRT: jim, prompt, image_geom, ndgrid
import MIRT: ir_mri_smap_r
import MIRT: ir_mri_smap1

using SpecialFunctions: ellipk, ellipe
using Test: @test, @testset, @test_throws, @inferred
using Plots: Plot, plot, plot!, contour!, scatter!, quiver, quiver!, Arrow

# many plots
function sensemap_test()
	nx,ny = 32,32
	smap, info = ir_mri_sensemap_sim(:all ;
		dims=(nx,ny), scale=:ssos_center,
	)
	smap = reshape(smap, nx, ny, 1, info.ncoilpr, info.nring)

	x = info.x
	y = info.y
	data = info.data
	down2 = a -> vec(a[1:3:end,1:3:end])
	down3 = a -> vec(a[1:3:end,1:3:end,1:3:end])

	for ir = 1:info.nring
		for ic = 1:info.ncoilpr
			pc = data[ir,ic]
			xr = pc.xr
			zr = pc.zr

			p = Array{Any}(undef, 7)

			# coordinates
			p[1] = jim(x, y, xr, xlabel="x", ylabel="y", "xr")
			p[2] = jim(x, y, zr, "zr")

			sx = pc.sx
			sy = pc.sy
			sz = pc.sz
			tmp = sqrt.(sx.^2 + sz.^2)
			(xx,yy) = ndgrid(x,y)

			p[3] = jim(x, y, sx, "sx")
			p[4] = jim(x, y, sy, "sy")
			p[5] = jim(x, y, sz, "sz")
			p[6] = quiver(down2(xx), down2(yy), title="(sx,sy)",
					quiver=(down2(sx./tmp), down2(sz./tmp)))

			# see final field components vs phase (for nz=0)
			(xx,yy) = ndgrid(x,y)
			bx = pc.bx
			by = pc.by
			bb = sqrt.(bx.^2 + by.^2)
			p[7] = jim(x, y, angle.(smap[:,:,1,ic,1]), "phase", color=:hsv)
			quiver!(down2(xx), down2(yy), title="(bx,by)",
					quiver=(down3((bx./bb)), down3((by./bb))))
			plot(p...)
			prompt()
		end
	end

end

# test code to explore when r is near 0
function ir_mri_smap_r_test()
	r0 = LinRange(0,5e-7,101)
	z0 = 0.4
	t0 = ir_mri_smap_r.(r0, z0)
	slope = 3*pi * z0 / ((1+z0^2)^2.5)
	plot(r0, t0, label="fun", xlabel="r")
	plot!(r0, slope * r0, label="line")
end

"""
ellipke_plot()
"""
function ellipke_plot()
	m = LinRange(0,1,101)
	(k,e) = (ellipk.(m), ellipe.(m))
	plot(xaxis=[0,1], yaxis=[0,3π/2])
	plot!(m, k, label="k", xlabel="m")
	plot!(m, e, label="e")
	plot!(ytick=((0:3)*pi/2, ["0", "π/2", "π", "3π/2"]))
end

"""
    ir_mri_sensemap_sim_test1()
test ir_mri_smap1 routine, cf Fig. 4 of grivich:00:tmf
"""
function ir_mri_sensemap_sim_test1()
	a = 1
	x = LinRange{Float32}(-2,2,99)
	y = LinRange{Float32}(-2,2,97)
#	zlist = Float32[0.001, 0.1, 0.2, 0.5, 1.0]
	zlist = Float32[0.1, 0.2, 0.5, 1.0]
#	(xx,yy,zz) = ndgrid(x, y, zlist)
	tmp = Iterators.product(x, y, zlist)
	tmp = [ir_mri_smap1(i[1], i[2], i[3], a) for i in tmp]
	smap_x = [p[1] for p in tmp]
	smap_y = [p[2] for p in tmp]
	smap_z = [p[3] for p in tmp]
#	(smap_x, smap_y, smap_z) = ir_mri_smap1.(xx, yy, zz, a)
	smap_b = @. sqrt(smap_x^2 + smap_y^2)
#	return plot(ir_mri_sensemap_sim_test1_show(smap_x, x, y, zlist, "x")...)

	return plot(layout=(4,length(zlist)),
		ir_mri_sensemap_sim_test1_show(smap_x, x, y, zlist, "x")...,
		ir_mri_sensemap_sim_test1_show(smap_y, x, y, zlist, "y")...,
		ir_mri_sensemap_sim_test1_show(smap_z, x, y, zlist, "z")...,
		ir_mri_sensemap_sim_test1_show(smap_b, x, y, zlist, "b")...,
	)
end

"""
    ir_mri_sensemap_sim_test1_show()
"""
function ir_mri_sensemap_sim_test1_show(smap, x, y, zlist, title)
	clim = (-20,20)
	nz = length(zlist)
	pl = Array{Any}(undef, nz)
	for iz = 1:nz
		pl[iz] = jim(x, y, smap[:,:,iz],
			title=title, clim=clim, xtick=[-2,2], ytick=[-2,2])
		blim = (zlist[iz] < 0.5) ? [7,12,19] : [1,3,5]
		contour!(pl[iz], x, y, abs.(smap[:,:,iz])', levels=blim,
			color=:blue, colorbar_entry=false)
		contour!(pl[iz], x, y, abs.(smap[:,:,iz])', levels=[0,0].+0.001,
			color=:green, colorbar_entry=false)
		title == "b" && plot!(xlabel="z = $(zlist[iz])")
	end
	return pl
end

"""
    ir_mri_sensemap_sim_show2()
returns plot
"""
function ir_mri_sensemap_sim_show2(smap, x, y, dx, dy, nlist, plist, rlist)
	if ndims(smap) == 3
		(nx,ny,ncoil) = size(smap)
	elseif ndims(smap) == 4
		(nx,ny,nz,ncoil) = size(smap)
	else
		throw("unknown ndims(smap) = $(ndims(smap))")
	end

	pl = Array{Any}(undef, ncoil, 3)
	clim = (0, maximum(abs.(smap)))
	xmax = maximum(abs.([vec(x); vec(y); vec(plist[:,:,[1,2]])]))
	ymax = xmax
	for ic=1:ncoil
		tmp = smap[:,:,ic]
		p = jim(x, y, abs.(tmp), clim=clim, "Magn. $ic")
		plot!(p, xlim=[-1,1]*1.1*xmax, xtick=(-1:1) * nx/2 * dx)
		plot!(p, ylim=[-1,1]*1.1*xmax, ytick=(-1:1) * nx/2 * dy)

		scatter!(p, [0], [0], marker=:o, label="", color=:green) # center
		scatter!(p, vec(plist[:,:,1]), vec(plist[:,:,2]),
			marker=:o, label="", color=:blue) # coil location
		xdir = nlist[ic,1,2]
		ydir = nlist[ic,1,1]
		r = rlist[ic,1]
		plot!(p, plist[ic,1,1].+r*xdir*[-1,1], plist[ic,1,2].+r*ydir*[1,-1],
			label="", color=:blue, linewidth=3) # coil
		pl[ic,1] = p

		ph = angle.(tmp) # show raw phase (understandable with hsv colormap)
		p = jim(x, y, ph, "Phase", clim=(-pi,pi), color=:hsv)
		plot!(p, xlim=[-1,1]*1.1*xmax, xtick=(-1:1) * nx/2 * dx)
		plot!(p, ylim=[-1,1]*1.1*xmax, ytick=(-1:1) * nx/2 * dy)
		pl[ic,2] = p
	end

	ssos = sqrt.(sum(abs2.(smap), dims=ndims(smap)))
	ssos = ssos / ssos[Int(end/2),Int(end/2)]

	p = jim(x, y, ssos, "SSoS (norm.)",
		#	xlim=[-1,1]*1.1*xmax,
			xtick=(-1:1) * nx/2 * dx,
			ytick=(-1:1) * ny/2 * dy)
	pl[1,3] = p

	if true # quiver plot for 1st coil
		bx = real(smap[:,:,1])
		by = imag(smap[:,:,1])
		(xx,yy) = ndgrid(x,y)
		down = a -> vec(a[1:3:end,1:3:end])
		p = quiver(down(xx), down(yy), quiver=(down(bx), down(by)),
				aspect_ratio = 1,
				arrow = Arrow(:simple, :head, 0.01, 0.01), # no effect!?
				title = "Field in x-y plane")
		plot!(p, xlim=[-1,1]*1.1*xmax, xtick=(-1:1) * nx/2 * dx)
		plot!(p, ylim=[-1,1]*1.1*xmax, ytick=(-1:1) * nx/2 * dy)
		pl[2,3] = p
	end

	for ic=3:ncoil
		pl[ic,3] = plot(xaxis=:off,yaxis=:off,grid=:off) # kludge
	end

	return plot(pl...)
end

"""
    ir_mri_sensemap_sim_test2( ; chat)
return plot with 2D example
"""
function ir_mri_sensemap_sim_test2()

	@test_throws String ir_mri_sensemap_sim(scale=:bug)
	@test_throws String ir_mri_sensemap_sim(nring=2, orbit_start = [1,2,3])

	(smap,t) = ir_mri_sensemap_sim(
		Vector{Tuple{Int,Int}}(undef, 0) ;
		dims=(32,32), scale=:ssos_center,
	)

	if true # check rotational symmetry in 4-coil case
		for ic=2:4
			tmp = rotl90(smap[:,:,1], ic-1)
			@test isapprox(abs.(tmp), abs.(smap[:,:,ic]))
			p1 = angle.(tmp) .+ (ic-1) * pi/2 # add pi/2 to rotated
			p2 = angle.(smap[:,:,ic])
			tmp = ComplexF32.(cis.(p2 - p1))
			@test isapprox(tmp, ones(size(tmp))) # trick: equivs mod 2*pi
		end
	end

	return ir_mri_sensemap_sim_show2(smap,
		t.x, t.y, t.dx, t.dy, t.nlist, t.plist, t.rlist)
end

"""
    ir_plot3_cube(x,y,z)
"""
function ir_plot3_cube(x,y,z)
	x1 = x[1]
	x2 = x[end]
	y1 = y[1]
	y2 = y[end]
	z1 = z[1]
	z2 = z[end]
	x = [x1,x2,x2,x1,x1,x1,x2,x2,x1,x1]
	y = [y1,y1,y2,y2,y1,y1,y1,y2,y2,y1]
	z = [z1,z1,z1,z1,z1,z2,z2,z2,z2,z2]
	plot(x,y,z, label="", xlabel="x", ylabel="y", zlabel="z")
end

"""
    ir_mri_sensemap_sim_show3()
shows coil geometry but not the 3D smap
"""
function ir_mri_sensemap_sim_show3(smap, x, y, z, dx, dy, dz,
	nlist, plist, rlist, olist, ulist, nring, ncoilpr, rcoil)

	pcolor = ['c', 'g', 'r']
	pcolor = i -> pcolor[1+rem(i,3)]

	ir_plot3_cube(x,y,z)
	scatter!(plist[:,:,1], plist[:,:,2], plist[:,:,3], label="") # coil centers

#=
	# coil normals - todo: 3d quiver not yet working
	tmp1 = reshape(plist, :, 3)
	tmp2 = reshape(nlist, :, 3)
	quiver!(tmp1[:,1], tmp1[:,2], tmp1[:,3],
		quiver=(tmp2[:,1], tmp2[:,2], tmp2[:,3]), label="3d")
	prompt()
=#

	if true # coils
		for ir = 1:nring
			for ic = 1:ncoilpr
				tmp = LinRange(0, 2*pi, 50)
				tmp = cos.(tmp) * olist[ic,ir,:]' + sin.(tmp) * ulist[ic,ir,:]'
				tmp = repeat(plist[ic,ir,:]', size(tmp,1), 1) + rcoil * tmp
				plot!(tmp[:,1], tmp[:,2], tmp[:,3], label="")
			#	patch(tmp[:,1], tmp[:,2], tmp[:,3], pcolor(ir),
			#		edgecolor=:none, facealpha=0.5)
			end
		end
	end
	plot!()
end

"""
    ir_mri_sensemap_sim_test3( ; chat)
return plot that illustrates 3D sense maps
"""
function ir_mri_sensemap_sim_test3( ; chat::Bool=false)
	nring = 3
	ncoil = 4 * nring
	ig = image_geom(nx=16, ny=14, nz=10, fov=200, dz=20, mask=:circ) # 20cm fov
#	ig = image_geom(nx=72, ny=48, nz=12, fov=22, zfov=10) % michelle

	(smap,t) = #@inferred # todo-i fails
		ir_mri_sensemap_sim(
			Vector{Tuple{Int,Int}}(undef, 0) ;
			dims=(ig.nx, ig.ny, ig.nz),
			dx=ig.dx, dz=ig.dz,
			orbit_start = 1*[0,45,0],
			rcoil=70, nring=nring, ncoil=ncoil,
		)

	if true
		tmp = smap .* repeat(ig.mask, 1,1,1,ncoil)
		jim(ncol=ig.nz, tmp, abswarn=false)
		chat && prompt()
		tmp = permutedims(tmp, [1,3,2,4]) # [nx nz ny ncoil] z cuts are smooth
		jim(ncol=ig.ny, tmp, abswarn=false)
		chat && prompt()
		jim(ncol=1, tmp[:,:,round(Int,end/2),:], abswarn=false)
		chat && prompt()
	end

	ir_mri_sensemap_sim_show3(smap, t.x, t.y, t.z, t.dx, t.dy, t.dz,
			t.nlist, t.plist, t.rlist, t.olist, t.ulist,
			t.nring, t.ncoilpr, t.rcoil)
end

@testset "1D" begin
	@inferred ir_mri_smap_r(5e-7, 0.4)
	plot(
		ir_mri_smap_r_test(),
		ellipke_plot(),
	)
end
prompt()

@testset "smap1" begin
	@inferred ir_mri_smap1(1, 2, 3, 4)
	ir_mri_sensemap_sim_test1() # basic test
end
prompt()

#@inferred ir_mri_sensemap_sim() # fails: Union of 3D, 4D

ir_mri_sensemap_sim_test2() # 2d test
prompt()

ir_mri_sensemap_sim_test3() # 3d test
prompt()

sensemap_test()

#=
@test typeof(ir_mri_sensemap_sim_test0()) <: Plots.Plot # ellipk
@test typeof(ir_mri_sensemap_sim_test1()) <: Plots.Plot # basic test
@test typeof(ir_mri_sensemap_sim_test2()) <: Plots.Plot # 2d test
@test typeof(ir_mri_sensemap_sim_test3()) <: Plots.Plot # 3d test
=#

#= for plotting
ir_mri_sensemap_sim_test3(chat=true)
ir_mri_sensemap_sim_test3(chat=true)
=#


text: include("freqtable.jl")


text:
using ShaleDrillingLikelihood.SDLParameters
using DataFrames
using Formatting
using TexTables
using JLD2, FileIO
using Dates
using LinearAlgebra
using Statistics
using DataStructures: OrderedDict

using ShaleDrillingLikelihood: idx_ρ, _ρ, _dρdθρ, theta_0, vw_revenue,
    split_theta, revenue,
    theta_produce_σ2η, theta_produce_σ2u, theta_produce_β,
    average_cost_df, theta_indexes

global MGL_CORP_INC_TAX
global CAPITAL_MGL_TAX_RATE

function sigma_epsilon_and_capital_tau(_alpha0, _gamma0, sig_u, sig_eta)
    NOMINAL_ANNUAL_DISCOUNT_RATE = 0.125
    BETA_ANNUAL_NOMINAL = 1/(1+NOMINAL_ANNUAL_DISCOUNT_RATE)
    MONTHLY_INFLATION = 0x1.006b55c832502p+0  #  1.00163780343638
    BETA_ANNUAL_REAL  = (MONTHLY_INFLATION)^12 / (1+NOMINAL_ANNUAL_DISCOUNT_RATE)

    MEDIAN_PERF              = 4428
    MEDIAN_LOGPERF           = log(MEDIAN_PERF)
    GAMMA_Q_T                = 0x1.32694dfd7a02ep+2  #  4.78767728570533    cumulative pdxn FE for pdxn month 240
    GATH_COMP_TRTMT_PER_MCF  = 0.42 + 0.07           # midstream charges
    global MGL_CORP_INC_TAX  = 0.402                 # mgl corp inc tax rate
    IDC_FRACTION             = 0.8                   # fraction of drilling costs that are intangible and expensed
    DEPRECIATION_YEARS       = 7                     # straight-line depreciation for remaining asset over these years
    WELL_COST                = 10                    # in millions of dollars
    MEDIAN_ROYALTY_RATE      = 0.2
    RELATIVE_CAPITAL_MGL_TAX = IDC_FRACTION + (1-IDC_FRACTION) / DEPRECIATION_YEARS * sum( BETA_ANNUAL_REAL.^(1:DEPRECIATION_YEARS) )
    global CAPITAL_MGL_TAX_RATE  = MGL_CORP_INC_TAX * RELATIVE_CAPITAL_MGL_TAX     # 0.3769469 vs 36.8% in HKL (Jan 2019 WP)

    eulers_gamma = Base.MathConstants.eulergamma

    _sigma_epsilon = exp(
        _gamma0 + (sig_u^2 + sig_eta^2)/2 + MEDIAN_LOGPERF + GAMMA_Q_T + log(1-MGL_CORP_INC_TAX) - _alpha0
    )

    return _sigma_epsilon, CAPITAL_MGL_TAX_RATE
end

texify(x::String) = x[1] == '\\' ? ("\$ " * x * " \$") : x

function coef_translate(nm)
    coefdict = Dict(
        "\\theta_\\rho" => "\\rho",
        "\\beta_1"    =>  "Log median house value",
        "\\beta_2"    =>  "Out-of-state owners (share)",
        "\\beta_3"    =>  "Pct impervious",
        "\\beta_4"    =>  "Log OGIP",
        "\\kappa_1"   =>  "0.125 | 0.1667",
        "\\kappa_2"   =>  "0.1667 | 0.1875",
        "\\kappa_3"   =>  "0.1875 | 0.2",
        "\\kappa_4"   =>  "0.2 | 0.225",
        "\\kappa_5"   =>  "0.225 | 0.25",
        "\\gamma_{3}" =>  "Intercept",
        "\\gamma_{10}" =>  "Intercept",
    )
    if nm in keys(coefdict)
        nm2 = coefdict[nm]
    else
        nm2 = nm
    end

    nm3 = replace.(nm2, r"(\d{4})\.\d+" => s"\1")

    return texify(nm3)
end

function nms_coef_se_sumstat(jld2file,
    DATADIR = "E:/projects/haynesville/intermediate_data";
    first_cost_year=FIRST_COST_YEAR, last_cost_year = LAST_COST_YEAR
)

    # retrieve estimates
    println_time_flush("Loading results from $jld2file")
    file = jldopen(jld2file, "r")
        DATAPATH = file["DATAPATH"]
        M        = file["M"]
        ddmnovf  = file["ddm_novf"]
        theta    = file["theta1"]
        LL       = file["LL"]
        g     = file["grad"]
        vcov  = file["invhess"]
    close(file)
    println_time_flush("loaded jld2")

    # info from DDM
    # -------------------

    REWARD = reward(ddmnovf)
    npsi = length(psispace(ddmnovf))
    nz = first(length.(zspace(ddmnovf)))

    # load in data from disk
    # -------------------
    rdatapath = joinpath(DATADIR, DATAPATH)
    println_time_flush("loading $rdatapath")

    data_royalty    = DataRoyalty(REWARD, rdatapath)
    data_produce    = DataProduce(REWARD, rdatapath)
    data_drill_prim = DataDrillPrimitive(REWARD, rdatapath)
    data_drill      = DataDrill(data_drill_prim, ddmnovf)

    # full dataset
    dataset_full = DataSetofSets(data_drill, data_royalty, data_produce, CoefLinks(data_drill))

    idx_d, idx_r, idx_p = theta_indexes(dataset_full)
    theta_d, theta_r, theta_p = split_thetas(dataset_full, theta)

    # make coefs
    # -------------------

    coef = theta
    se = sqrt.(diag(vcov))
    nms = coef_translate.(coefnames(dataset_full))

    # update rho
    # -------------------

    idx_rho = idx_ρ(REWARD)
    thetarho = coef[idx_rho]
    se[idx_rho] *= _dρdθρ(thetarho)
    coef[idx_rho] = _ρ(thetarho)

    # scale
    # -------------------

    theta_rev = vw_revenue(REWARD, theta_d)
    _alpha0 = theta_0(revenue(REWARD), theta_rev)
    _gamma0 = last(theta_produce_β(data_produce, theta_p))
    sig_eta = theta_produce_σ2η(data_produce, theta_p)
    sig_u   = theta_produce_σ2u(data_produce, theta_p)

    sigeps, captau = sigma_epsilon_and_capital_tau(_alpha0, _gamma0, sig_u, sig_eta)

    # costs
    # --------

    cost_df = average_cost_df(dataset_full, theta)
    date_range = first_cost_year .<= year.(cost_df[!,:date]) .<= last_cost_year
    avgcost = mean( cost_df[date_range,:cost2] ) * sigeps / (1-captau)

    sumstats = OrderedDict(
        "\$\\sigma_\\epsilon\$" => sigeps,
        "Avg cost" => -avgcost,
        "Log lik" => LL,
        "Num \$z\$ gridpoints" => nz,
        "Num \$\\psi\$ gridpoints" => npsi,
        "Num simulations" => M,
    )

    idxnew = vcat(idx_r[2:end], idx_d, idx_p[end-2:end])

    return (nms[idxnew], coef[idxnew], se[idxnew], sumstats)

end

function regcol(title, nms, coefs, ses, sumstats)

    rc = RegCol(title)

    for x in zip(nms, coefs, ses)
        setcoef!(rc, x...)
    end

    setstats!(rc, sumstats)
    rc["Avg cost"].format = "{:.2f}"
    rc["Log lik"].format = "{:.2f}"

    return rc
end

addlinespace(textbl)             = replace(textbl, r"(\\\\)(\n)(?!\s*(Log [Ll]ik|\$ \\sigma\\_\\epsilon \$|Avg cost|Num|\\bottom|&))" => s"\1 \\addlinespace[1pt] \2")
hline_to_midrule(textbl)         = replace(textbl, r"hline" => s"cmidrule(lr)") # cmidrule(lr){2-5}\\cmidrule(lr){6-9}")
header(capt,lab)                 = "\\begin{table}\n\\centering\n\\begin{adjustbox}{max totalheight={\\textheight}, max width={\\textwidth}}\n\\begin{threeparttable}\n\t\\caption{$capt}\\label{$lab}\n\n"
tablenotes(note)                 = "\\begin{tablenotes}[flushleft]\n\\item\\scriptsize{\n\t" * note * "\n}\n\\end{tablenotes}\n"
right_to_left(textbl)            = replace(textbl, r"\\begin\{tabular\}\{r\|([^\}]+)\}" => s"\\begin{tabular}{l\1}")
footer()                         = "\\end{threeparttable}\n\\end{adjustbox}\n\\end{table}"
add_leasing_heading(textable)    = replace(textable, r"(\n +\$ *\\psi\^0 *\$)"   => s" & \\multicolumn{4}{c}{\\emph{Leasing}   } \\\\ \\addlinespace[1pt]\1")
add_drilling_heading(textable)   = replace(textable, r"(\n +\$ *\\alpha\\_\{c,2008} *\$)" => s" & \\multicolumn{4}{c}{\\emph{Drilling}  } \\\\ \\addlinespace[1pt]\1")
add_production_heading(textable) = replace(textable, r"(\n +Intercept)"          => s" & \\multicolumn{4}{c}{\\emph{Production}} \\\\ \\addlinespace[1pt]\1")
fix_underscores(textable)        = replace(textable, r"\\_" => s"_")

function make_table(rt::IndexedTable, caption::String, label::String, note::String)
    out = header(caption, label) *
    (rt |> to_tex |> addlinespace |> hline_to_midrule |> right_to_left |> add_leasing_heading |> add_drilling_heading |> add_production_heading |> fix_underscores) *
    tablenotes(note) *footer()
    return out
end


text:
widen_eltype(::Type{T}) where {T<:Number} = widen(T)
widen_eltype(::Type{SVector{N,T}}) where {N,T<:Number} = SVector{N,widen(T)}
widen_eltype(::Type{Vector{T}}) where {T<:Number} = Vector{widen(T)}

# We test the generic functionality of a domain.
# These tests check whether the given domain correctly implements the
# interface of a domain.
function test_generic_domain(d::Domain)
    @test isreal(d) == isreal(eltype(d))
    @test isreal(d) == isreal(numtype(d))

    @test convert(Domain{eltype(d)}, d) == d
    @test convert(Domain{widen_eltype(eltype(d))}, d) == d
    @test prectype(convert_prectype(d, BigFloat)) == BigFloat

    if !isempty(d)
        x = point_in_domain(d)
        @test x ∈ d
        @test approx_in(x, d, 0.01)
        @test_throws ErrorException approx_in(x, d, -1)
        @test in.([x,x], d) == in.([x,x], Ref(d))
        @test approx_in.([x,x], d, 0.01) == approx_in.([x,x], Ref(d), 0.01)
    else
        try
            x = point_in_domain(d)
            @test false
        catch
        end
    end
    @test canonicaldomain(DomainSets.Equal(), d) == d
    if canonicaldomain(d) == d
        @test mapto_canonical(d) == IdentityMap{eltype(d)}(dimension(d))
        @test mapfrom_canonical(d) == IdentityMap{eltype(d)}(dimension(d))
    else
        cd = canonicaldomain(d)
        @test mapfrom_canonical(d) == mapto(cd, d)
        @test mapto_canonical(d) == mapto(d, cd)
    end
    if iscomposite(d)
        @test ncomponents(d) == length(components(d))
        els = components(d)
        @test all([component(d,i) == els[i] for i in 1:ncomponents(d)])
    end
end

@testset "generic domain interface" begin
    domains = [
        0..1,
        UnitInterval(),
        ChebyshevInterval(),
        HalfLine(),
        NegativeHalfLine(),
        UnitInterval()^3,
        (0..1) × (2.0..3) × (3..4),
        UnitBall(),
        VectorUnitBall(),
        VectorUnitBall(8),
        UnitDisk(),
        VectorUnitDisk(),
        UnitCircle(),
        VectorUnitCircle(),
        UnitSphere(),
        VectorUnitSphere(),
        UnitSimplex(Val(2)),
        UnitSimplex(2),
        DomainSets.WrappedDomain(0..2.0)
    ]

    @testset "generic domains" begin
        for domain in domains
            test_generic_domain(domain)
        end
    end

    @testset "generic functionality" begin
        struct SomeDomain <: Domain{Float64}
        end
        @test_throws MethodError 0.5 ∈ SomeDomain()
        @test_throws MethodError approx_in(0.5, SomeDomain())

        if VERSION < v"1.6-"
            @test_logs (:warn, "`in`: incompatible combination of point: Tuple{Float64,Float64} and domain eltype: SArray{Tuple{2},Float64,1,2}. Returning false.") (0.5,0.2) ∈ UnitCircle()
            @test_logs (:warn, "`in`: incompatible combination of point: Tuple{Float64,Float64} and domain eltype: SArray{Tuple{2},Float64,1,2}. Returning false.") approx_in((0.5,0.2), UnitCircle())
        else
            @test_logs (:warn, "`in`: incompatible combination of point: Tuple{Float64, Float64} and domain eltype: SVector{2, Float64}. Returning false.") (0.5,0.2) ∈ UnitCircle()
            @test_logs (:warn, "`in`: incompatible combination of point: Tuple{Float64, Float64} and domain eltype: SVector{2, Float64}. Returning false.") approx_in((0.5,0.2), UnitCircle())
        end

        # some functionality in broadcast
        @test 2 * (1..2) == 2 .* (1..2)
        @test (1..2) * 2 == (1..2) .* 2
        @test (1..2) / 2 ≈ (0.5..1)
        @test 2 \ (1..2) ≈ (0.5..1)
        @test_throws MethodError (0..1) + 0.4

        # promotion
        @test DomainSets.promote_domains() == ()
        s1 = Set([0..1,2..3,3..4.0])
        @test s1 isa Set{<:Domain{Float64}}
        @test DomainSets.promote_domains(s1) == s1
        s2 = Set([0..1,2..3,3..4.0, Point(2)])
        @test s2 isa Set{Domain}
        @test DomainSets.promote_domains(s2) isa Set{<:Domain{Float64}}
        @test DomainSets.promote(0..1.0, [1,2,3]) isa Tuple{Interval,Vector{Float64}}
        @test DomainSets.promote([1,2,3], 0..1.0) isa Tuple{Vector{Float64},Interval}

        # compatible point-domain pairs
        @test DomainSets.iscompatiblepair(0.5, 0..1)
        @test DomainSets.iscompatiblepair(0.5, 0..1.0)
        @test !DomainSets.iscompatiblepair(0.5, UnitBall())
        @test DomainSets.iscompatiblepair(0.5, [0.3])
        @test DomainSets.iscompatiblepair(0, [0.3])
        @test DomainSets.iscompatiblepair(0.5, [1, 2, 3])
        @test DomainSets.iscompatiblepair(0, Set([1, 2, 3]))
        @test DomainSets.iscompatiblepair(0.5, Set([1, 2, 3]))
        @test DomainSets.iscompatiblepair(0, Set([1.0, 2, 3]))

        @test DomainSets.convert_eltype(Float64, Set([1,2])) isa Set{Float64}
        @test_throws ErrorException DomainSets.convert_eltype(Float64, (1,2)) isa NTuple{2,Float64}
        @test DomainSets.convert_eltype(Int, (1,2)) == (1,2)
    end
end


text:
#__precompile__()
module NumericalBasis

#   This file is part of NumericalBasis.jl. It is licensed under the MIT license
#   Copyright (C) 2018 Michael Reed

using SyntaxTree, ForceImport, Reduce
export floatset, optimal, roots, polyhorner, polyfactors, polyexpand

function floatset(T::DataType,N;scale=x->x)
    l = scale(eps(T))
    u = scale(prevfloat(T(Inf)))
    return l:(u-l)/(N-1):u
end

polyhorner(x,a) = polyhorner(x,a,1)
polyhorner(x,a::Array{<:Any,1},k) = k==length(a) ? a[k] : Algebra.:+(a[k],Algebra.:*(x,polyhorner(x,a,k+1)))

polyfactors(x,a) = polyfactors(x,a,1)
polyfactors(x,a::Array{<:Any,1},k) = k==length(a) ? Algebra.:-(x,a[k]) : Algebra.:*(Algebra.:-(x,a[k]),polyfactors(x,a,k+1))

polyexpand(x,a) = polyexpand(x,a,length(a))
polyexpand(x,a::Array{<:Any,1},k) = k==1 ? a[k] : Algebra.:+(Algebra.:*(a[k],Algebra.:^(x,k-1)),polyexpand(x,a,k-1))

roots = factor

function optimal(expr)
    h = horner(expr)
    f = factor(h)
    eh = SyntaxTree.exprval(h)[1]
    ef = SyntaxTree.exprval(f)[1]
    eo = SyntaxTree.exprval(expr)[1]
    if eh ≤ ef
        return eh ≤ eo ? h : expr
    else
        return ef ≤ eo ? f : expr
    end
end

__init__() = nothing

end # module


text:
function mdl_line_ode!(
    device_states::AbstractArray{T},
    output_ode::AbstractArray{T},
    voltage_r_from::T,
    voltage_i_from::T,
    voltage_r_to::T,
    voltage_i_to::T,
    current_r_from::AbstractArray{T},
    current_i_from::AbstractArray{T},
    current_r_to::AbstractArray{T},
    current_i_to::AbstractArray{T},
    branch::BranchWrapper,
) where {T <: ACCEPTED_REAL_TYPES}
    L = PSY.get_x(branch)
    R = PSY.get_r(branch)
    ω_b = get_system_base_frequency(branch) * 2 * π

    Il_r = device_states[1]
    Il_i = device_states[2]
    output_ode[1] = (ω_b / L) * ((voltage_r_from - voltage_r_to) - (R * Il_r - L * Il_i))
    output_ode[2] = (ω_b / L) * ((voltage_i_from - voltage_i_to) - (R * Il_i + L * Il_r))

    current_r_from[1] -= Il_r
    current_i_from[1] -= Il_i
    current_r_to[1] += Il_r
    current_i_to[1] += Il_i
end


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

export sv_ecef_to_ecef

"""
    sv_ecef_to_ecef(sv::OrbitStateVector, args...)

Convert the orbit state vector `sv` from an ECEF frame to another ECEF frame.
The arguments `args...` must match those of the function
[`r_ecef_to_ecef`](@ref) **wihtout** the rotation representation.
"""
function sv_ecef_to_ecef(sv::OrbitStateVector, args...)
    D = r_ecef_to_ecef(DCM, args...)

    # Since both frames does not have a significant angular velocity between
    # them, then we just need to convert the representations.
    return OrbitStateVector(t = sv.t, r = D * sv.r, v = D * sv.v, a = D * sv.a)
end


text:
struct AffineMean{T<:Real,V<:AbstractVector{T},O} <: PriorMean{T}
    w::V
    b::Vector{T}
    nDim::Int
    opt::O
end

"""
    AffineMean(w::Vector, b::Real; opt = ADAM(0.01))
    AffineMean(dims::Int; opt=ADAM(0.01))

## Arguments
- `w::Vector` : Weight vector
- `b::Real` : Bias
- `dims::Int` : Number of features per vector

Construct an affine operation on `X` : `μ₀(X) = X * w + b` where `w` is a vector and `b` a scalar
Optionally give an optimiser `opt` (`Adam(α=0.01)` by default)
"""
function AffineMean(w::V, b::T; opt=ADAM(0.01)) where {T<:Real,V<:AbstractVector{T}}
    return AffineMean{T,V,typeof(opt)}(copy(w), [b], length(w), opt)
end

function AffineMean(dims::Int; opt=ADAM(0.01))
    return AffineMean{Float64,Vector{Float64},typeof(opt)}(
        randn(Float64, dims), [0.0], dims, opt
    )
end

function Base.show(io::IO, ::MIME"text/plain", μ₀::AffineMean)
    return print(
        io, "Affine Mean Prior (size(w) = ", length(μ₀.w), ", b = ", first(μ₀.b), ")"
    )
end

function (μ₀::AffineMean{T})(x::AbstractVector) where {T<:Real}
    # μ₀.nDim == size(x, 1) || error(
    # "Number of dimensions of prior weight W (",
    # size(μ₀.w),
    # ") and X (",
    # size(x),
    # ") do not match",
    # )
    return dot.(x, Ref(μ₀.w)) .+ first(μ₀.b)
end

function init_priormean_state(hyperopt_state, μ₀::AffineMean)
    μ₀_state = (; w=Optimisers.init(μ₀.opt, μ₀.w), b=Optimisers.init(μ₀.opt, μ₀.b))
    return merge(hyperopt_state, (; μ₀_state))
end

function update!(μ₀::AffineMean{T}, hyperopt_state, grad) where {T<:Real}
    μ₀_state = hyperopt_state.μ₀_state
    w, Δw = Optimisers.apply(μ₀.opt, μ₀_state.w, μ₀.w, grad.w)
    b, Δb = Optimisers.apply(μ₀.opt, μ₀_state.b, μ₀.b, grad.b)
    μ₀.w .+= Δw
    μ₀.b .+= Δb
    return merge(hyperopt_state, (; μ₀_state=(; w, b)))
end


text:
using Distributions

include("confidence_interval.jl")

function one_prop_z_int(p, n, lvl)
    distribution = Normal()

    standard_error = sqrt((p * (1-p)) / n)
    critical_value = quantile(distribution, lvl + (1 - lvl) / 2)

    lowest_value = round(p - critical_value * standard_error, 4)
    highest_value = round(p + critical_value * standard_error, 4)

    return ConfidenceInterval(lowest_value, highest_value)
end

function one_prop_z_test(p, p0, n, alpha, alt)
    distribution = Normal(p0, sqrt((p0 * (1-p0)) / n))

    if alt == "<"
        return round(cdf(distribution, p), 4)
    elseif alt == ">"
        return round(1 - cdf(distribution, p), 4)
    elseif alt == "!="
        if p > p0
            return round(2 * (1 - cdf(distribution, p)), 4)
        else
            return round(2 * cdf(distribution, p), 4)
        end
    else
        return -1
    end
end

function two_prop_z_test(p1, n1, p2, n2, alpha, alt)
    p_pooled = (p1*n1 + p2*n2) / (n1 + n2)
    p_diff = p1 - p2
    distribution = Normal(0, sqrt((p_pooled*(1-p_pooled)/n1) + (p_pooled*(1-p_pooled)/n2)))

    if alt == "<"
        return round(cdf(distribution, p_diff), 4)
    elseif alt == ">"
        return round(1 - cdf(distribution, p_diff), 4)
    elseif alt == "!="
        if p_diff > 0
            return round(2 * (1 - cdf(distribution, p_diff)), 4)
        else
            return round(2 * cdf(distribution, p_diff), 4)
        end
    else
        return -1
    end
end

function two_prop_z_int(p1, n1, p2, n2, lvl)
    standard_error = sqrt(((p1 * (1-p1)) / n1) + ((p2 * (1-p2)) / n2))

    standard_normal = Normal()
    critical_value = quantile(standard_normal, lvl + (1 - lvl) / 2)

    lowest_value = round((p1-p2) - critical_value * standard_error, 4)
    highest_value = round((p1-p2) + critical_value * standard_error, 4)

    return ConfidenceInterval(lowest_value, highest_value)
end

function z_interval(x̄, sx, n, lvl)
    distribution = Normal()

    if n < 30
        println("Conditions weren't met, using t-distribution...\n")
        distribution = TDist(n - 1)
    end

    standard_error = sx / sqrt(n)
    critical_value = quantile(distribution, lvl + (1 - lvl) / 2)

    lowest_value = round(x̄ - critical_value * standard_error, 4)
    highest_value = round(x̄ + critical_value * standard_error, 4)

    return ConfidenceInterval(lowest_value, highest_value)
end

function z_test(x̄, μ0, sx, n, alpha, alt)
    distribution = Normal(μ0, sx / sqrt(n))

    if alt == "<"
        return round(cdf(distribution, x̄), 4)
    elseif alt == ">"
        return round(1 - cdf(distribution, x̄), 4)
    elseif alt == "!="
        if x̄ > μ0
            return round(2 * (1 - cdf(distribution, x̄)), 4)
        else
            return round(2 * cdf(distribution, x̄), 4)
        end
    else
        return -1
    end
end

function two_sample_z_interval(x̄1, sx1, n1, x̄2, sx2, n2, lvl)
    distribution = Normal()

    if n1 < 30 || n2 < 30
        println("Conditions weren't met, using t-distribution...\n")
        distribution = TDist(n1 + n2 - 2)
    end

    standard_error = sqrt((((n1 - 1) * sx1 ^ 2) + ((n2 - 1) * sx2 ^ 2)) / (n1 + n2 - 2)) * sqrt(1 / n1 + 1 / n2)
    critical_value = quantile(distribution, lvl + (1 - lvl) / 2)

    lowest_value = round(x̄1 - x̄2 - critical_value * standard_error, 4)
    highest_value = round(x̄1 - x̄2 + critical_value * standard_error, 4)

    return ConfidenceInterval(lowest_value, highest_value)
end

function two_sample_z_test(x̄1, sx1, n1, x̄2, sx2, n2, lvl, alt)
    z = (x̄1 - x̄2) / sqrt((sx1^2 / n1 + sx2^2 / n2))

    standard_normal = Normal()

    if alt == "<"
        return round(cdf(standard_normal, abs(z)), 4)
    elseif alt == ">"
        return round(1 - cdf(standard_normal, abs(z)), 4)
    elseif alt == "!="
        return round(2 * (1 - cdf(standard_normal, abs(z))), 4)
    else
        return -1
    end
end

max_stars_repo_name: jackcook/statulator

max_issues_repo_name: jackcook/statulator

max_forks_repo_name: jackcook/statulator


text:
using OMEinsum

function contract(gb::EinGraph)
    EinCode((gb.labels...,), ())(gb.tensors...)[]
end


text:
include("mod_todo.jl")
using .ToDoApp, Dates

todo1 = ToDo(1, "Getting groceries", false, Date("2022-04-01", "yyyy-mm-dd"), 5)
# ToDoApp.ToDo(...) is also correct

print_todo(todo1)
# I still have to do: Getting groceries
# A todo created at: 2022-04-01

# helper(todo1)
# ERROR: UndefVarError: helper not defined


text:
#=

  A Round of Golf puzzle in Julia ConstraintSolver.jl

  From http://brownbuffalo.sourceforge.net/RoundOfGolfClues.html
  """
  Title: A Round of Golf
  Author: Ellen K. Rodehorst
  Publication: Dell Favorite Logic Problems
  Issue: Summer, 2000
  Puzzle #: 9
  Stars: 1

  When the Sunny Hills Country Club golf course isn't in use by club members,
  of course, it's open to the club's employees. Recently, jack and three other
  workers at the golf course got together on their day off to play a round of
  eighteen holes of golf.
  Afterward, all four, including Mr. green, went to the clubhouse to total
  their scorecards. Each man works at a different job (one is a short-order
  cook), and each shot a different score in the game. No one scored below
  70 or above 85 strokes. From the clues below, can you discover each man's
  full name, job and golf score?

  1. bill, who is not the maintenance man, plays golf often and had the lowest
  score of the foursome.
  2. Mr. clubb, who isn't paul, hit several balls into the woods and scored ten
  strokes more than the pro-shop clerk.
  3. In some order, frank and the caddy scored four and seven more strokes than
  Mr. sands.
  4. Mr. carter thought his score of 78 was one of his better games, even
     though frank's score  was lower.
  5. None of the four scored exactly 81 strokes.

  Determine: First Name - Last Name - Job - score
  """

  Compare with the F1 model:
  http://www.f1compiler.com/samples/A 20Round 20of 20Golf.f1.html

  Solution:
             jack, bill, paul, frank
             clubb sands carter green
             maint cook  caddy clerk
             85    71    78    75
  first_name: [1, 2, 3, 4]
  last_name : [4, 1, 2, 3]
  job       : [2, 1, 4, 3]
  score     : [85, 71, 78, 75]

  Note: This model is quite messy since ConstraintSolver.jl v0.5.3 don't
  have support for a good element (x[y] == z) and a little limited form
  of indicator/reification. Ole is working on all of this as we speak... :-)

  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 a_round_of_golf(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"=>true,
                                                            "logging"=>[],

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

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

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

                                                            "time_limit"=>3,

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

    n = 4
    jack,bill,paul,frank = 1:n
    first_name = [jack,bill,paul,frank]

    @variable(model, 1 <= last_name[1:n] <= n, Int)
    green,clubb,sands,carter = last_name

    @variable(model, 1 <= job[1:n] <= n, Int)
    cook,maintenance_man,clerk,caddy = job

    @variable(model, 70 <= score[1:n] <= 85, Int)
    score_jack,score_bill,score_paul,score_frank = score

    @constraint(model, last_name in CS.AllDifferentSet())
    @constraint(model, job in CS.AllDifferentSet())
    @constraint(model, score in CS.AllDifferentSet())

    #  1. bill, who is not the maintenance man, plays golf often and had
    #     the lowest score of the foursome.
    @constraint(model, bill != maintenance_man)

    @constraint(model, score_bill <= score_jack) # should be '<'
    @constraint(model, score_bill <= score_paul)
    @constraint(model, score_bill <= score_frank)

    #  2. Mr. clubb, who isn"t paul, hit several balls into the woods and
    #     scored ten strokes more than the pro-shop clerk.
    @constraint(model, clubb != paul)

    @variable(model, 70 <= score_clubb <= 85, Int)
    @variable(model, 70 <= score_clerk <= 85, Int)
    my_element(model,clubb,score,score_clubb)
    my_element(model,clerk,score,score_clerk)
    @constraint(model, score_clubb == score_clerk + 10)

    #  3. In some order, frank and the caddy scored four and seven more
    #     strokes than Mr. sands.
    @constraint(model, frank != caddy)
    @constraint(model, frank != sands)
    @constraint(model, caddy != sands)

    @variable(model, 70 <= score_sands <= 85, Int)
    @variable(model, 70 <= score_caddy <= 85, Int)
    @variable(model, 70 <= score_carter <= 85, Int)
    my_element(model, sands,score,score_sands)
    my_element(model, caddy,score,score_caddy)
    my_element(model, carter,score,score_carter)
    #=
    # Picat code:
    (
       (scorefrank = scoresands + 4 #/\
       scorecaddy = scoresands + 7)
       #\/
       (scorefrank = scoresands + 7 #/\
       scorecaddy = scoresands + 4)
    ),
    =#
    c1 = @variable(model, [1:2], Bin )
    @constraint(model, c1[1] := {score_frank == score_sands + 4} )
    @constraint(model, c1[2] := {score_caddy == score_sands + 7} )

    c2 = @variable(model, [1:2], Bin )
    @constraint(model, c2[1] := {score_frank == score_sands + 7} )
    @constraint(model, c2[2] := {score_caddy == score_sands + 4} )

    c3 = @variable(model, [1:2], Bin )
    @constraint(model, c3[1] := {sum(c1) == 2})
    @constraint(model, c3[2] := {sum(c2) == 2})

    c = @variable(model, [1:1], Bin )
    @constraint(model, c[1] := {sum(c3) == 1})
    @constraint(model, c[1] == 1)

    #  4. Mr. carter thought his score of 78 was one of his better games, even
    #  though frank"s score was lower.
    @constraint(model,frank != carter)

    #  score[carter] = 78,
    @constraint(model, score_carter == 78)
    @constraint(model, score_frank <= score_carter) # should be '<'

    #  5. None of the four scored exactly 81 strokes.
    for s in score
        @constraint(model, s != 81)
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    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")
                last_namex = convert.(Integer,JuMP.value.(last_name; result=sol))
                jobx = convert.(Integer,JuMP.value.(job; result=sol))
                scorex = convert.(Integer,JuMP.value.(score; result=sol))
                c1x = convert.(Integer,JuMP.value.(c1; result=sol))
                c2x = convert.(Integer,JuMP.value.(c2; result=sol))
                c3x = convert.(Integer,JuMP.value.(c3; result=sol))
                cx = convert.(Integer,JuMP.value.(c; result=sol))
                println("last_name:$last_namex")
                println("job      :$jobx")
                println("score    :$scorex")
                println("c1:$c1x c2:$c2x c3:$c3x c:$cx")
                println()
            end
        end
    else
        println("status:$status")
    end

    return status
end

@time a_round_of_golf()

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using ExtremeVertexDesigns
using Test

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


text:
@testset "OpenCL.Context" begin
    @testset "OpenCL.Context constructor" begin
        @test_throws MethodError (cl.Context([]))
        for platform in cl.platforms()
            for device in cl.devices(platform)
                ctx = cl.Context(device)
                @test ctx != nothing
                finalize(ctx)
            end
        end
    end

    @testset "OpenCL.Context platform properties" begin
        for platform in cl.platforms()
            try
                cl.Context(cl.CL_DEVICE_TYPE_CPU)
            catch err
                @test typeof(err) == cl.CLError
                # CL_DEVICE_NOT_FOUND could be throw for GPU only drivers
                @test err.desc in (:CL_INVALID_PLATFORM,
                                         :CL_DEVICE_NOT_FOUND)
            end

            if platform[:name] == "Portable Computing Language"
                warn("Skipping OpenCL.Context platform properties for " *
                     "Portable Computing Language Platform")
                continue
            end

            properties = [(cl.CL_CONTEXT_PLATFORM, platform)]
            for (cl_dev_type, sym_dev_type) in [(cl.CL_DEVICE_TYPE_CPU, :cpu),
                                                (cl.CL_DEVICE_TYPE_GPU, :gpu)]
                if !cl.has_device_type(platform, sym_dev_type)
                    continue
                end
                @test cl.Context(sym_dev_type, properties=properties) != nothing
                @test cl.Context(cl_dev_type, properties=properties) != nothing
                ctx = cl.Context(cl_dev_type, properties=properties)
                @test isempty(cl.properties(ctx)) == false
                test_properties = cl.properties(ctx)

                @test test_properties == properties

                platform_in_properties = false
                for (t, v) in test_properties
                    if t == cl.CL_CONTEXT_PLATFORM
                        @test v[:name] == platform[:name]
                        @test v == platform
                        platform_in_properties = true
                        break
                    end
                end
                @test platform_in_properties
            end
            try
                ctx2 = cl.Context(cl.CL_DEVICE_TYPE_ACCELERATOR,
                                  properties=properties)
            catch err
                @test typeof(err) == cl.CLError
                @test err.desc == :CL_DEVICE_NOT_FOUND
            end
        end
    end

    @testset "OpenCL.Context create_some_context" begin
        @test cl.create_some_context() != nothing
        @test typeof(cl.create_some_context()) == cl.Context
    end

   @testset "OpenCL.Context parsing" begin
        for platform in cl.platforms()
            properties = [(cl.CL_CONTEXT_PLATFORM, platform)]
            parsed_properties = cl._parse_properties(properties)

            @test isodd(length(parsed_properties))
            @test parsed_properties[end] == 0
            @test parsed_properties[1] == cl.cl_context_properties(cl.CL_CONTEXT_PLATFORM)
            @test parsed_properties[2] == cl.cl_context_properties(platform.id)
        end
    end
end