text:
# Ported to Julia by Andy Ferris, 2016, and re-released under an MIT license.
#/**
# * Copyright (c) Charles Karney (2008-2015) <charles@karney.com> and licensed
# * under the MIT/X11 License.  For more information, see
# * http://geographiclib.sourceforge.net/
# **********************************************************************/

"""
eatanhe(x, e2) returns e*atanh(e*x) where e = sqrt(e2)
e2 is a scalar; x can be any shape.
"""
function eatanhe(x, e2)
    e = sqrt(abs(e2))
    if (e2 >= 0)
        y = e * atanh(e * x)
    else
        y = -e * atan(e * x)
    end
end

"""
%TAUF   tan(phi)
%
%   TAUF(taup, e2) returns tangent of phi in terms of taup the tangent of
%   chi.  e2, the square of the eccentricity, is a scalar; taup can be any
%   shape.
"""
function tauf(taup, e2)
    numit = 5
    e2m = 1 - e2     # This possibly disagrees with the C++ code in the prolate case
    tau = taup / e2m
    stol = 0.1 * sqrt(eps(Float64)) * max(1.0, abs(taup))
    g = isfinite(tau)
    for i = 1 : numit
        if !g
            break
        end
        tau1 = hypot(1, tau)
        sig = sinh( eatanhe( tau / tau1, e2 ) )
        taupa = hypot(1, sig) * tau - sig * tau1
        dtau = (taup - taupa) * (1 + e2m * tau*tau) / (e2m * tau1 * hypot(1, taupa))
        tau = tau + dtau
        g = g & (abs(dtau) >= stol)
    end
    return tau
end

"""
taupf(tau, e2) returns tangent of chi in terms of tau the tangent of
phi.  e2, the square of the eccentricity, is a scalar; tau can be any
shape.
"""
function taupf(tau, e2)
  tau1 = hypot(1, tau)
  sig = sinh( eatanhe( tau / tau1, e2 ) )
  return hypot(1, sig) * tau - sig * tau1
end

# Code adapted from MATLAB port of geograhiclib, which was licensed under
# a permissive MIT-like license and the following copyright notice
#     Copyright (c) 2016, Charles Karney
#     All rights reserved.

"""
Evaluate polynomial with coefficients p[i] at x
"""
function polyval(p, x)
    y = 1
    out = p[end]
    for i = length(p)-1:-1:1
        y *= x
        out += p[i] * y
    end
    return out
end

"""
LatFix(x) returns x is it is in the range [-90, 90]; otherwise it
returns NaN.
"""
function LatFix(x::Number)
    y = x
    if abs(x) > 90
        y = NaN
    end
    return y
end

function AngDiff(x,y)
    AngNormalize(AngNormalize(y) - AngNormalize(x))
    #(d, t) = sumx(AngNormalize(x), AngNormalize(-y))
    #d = - AngNormalize(d)
    #if d == 180 && t < 0
    #    d = -180
    #end
    #return sumx(d, -t)[1]
end

"""
%AngNormalize:  Reduce angle to range [-180, 180)
%
%   x = AngNormalize(x) reduces angles to the range [-180, 180).
"""
function AngNormalize(x::Number)
    x = rem(x, 360)
    if x >= 180
        x = x - 360
    elseif x < -180
        x = x + 360
    end

    return x
end

function atan2d(y, x)
    # In order to minimize round-off errors, this function rearranges the
    # arguments so that result of atan is in the range [-pi/4, pi/4] before
    # converting it to degrees and mapping the result to the correct
    # quadrant.
    q = 0
    if (abs(y) > abs(x))
        tmp = x
        x = y
        y = tmp
        q = 2
    end
    if (x < 0)
        x = -x
        q = q + 1
    end
    # here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4]
    ang = atan(y, x) * 180 / pi
    # Note that atan2d(-0.0, 1.0) will return -0.  However, we expect that
    # atan2d will not be called with y = -0.  If need be, include
    #
    #   case 0: ang = 0 + ang; break;
    #
    # and handle mpfr as in AngRound.
    if q == 1
        ang = (y > 0 ? 180 : -180) - ang
    elseif q == 2
        ang =  90 - ang
    elseif q == 3
        ang = -90 + ang
    end
    return ang
end

"""
    immutable TransverseMercator{MaxPow}
    TransverseMercator(datum)
    TransverseMercator(ellipsoid)
    TransverseMercator(a, f, [::Type{Val{MaxPow}} = Val{6}])

Cache of ellipsoidal calculations necessary for transverse-Mercator and
polar-stereographic transformations. Series expansion coefficients up to order
`MaxPow` (between 4 and 8, default 6) are calculated and stored for fast
transverse-Mercator and UTM calculations.
"""
struct TransverseMercator{MaxPow}
    a::Float64
    f::Float64
    e2::Float64
    es::Float64
    e2m::Float64
    c::Float64
    n::Float64
    a1::Float64
    b1::Float64
    alp::NTuple{MaxPow, Float64}
    bet::NTuple{MaxPow, Float64}
end

TransverseMercator(datum) = TransverseMercator(ellipsoid(datum))
TransverseMercator(el::Ellipsoid) = TransverseMercator(el.a, 1-el.b/el.a)
TransverseMercator(a,f) = TransverseMercator(a,f,Val{6}) # Default to sixth-order expansion
function TransverseMercator(a::Float64, f::Float64, ::Type{Val{MaxPow}}) where MaxPow
    e2 = f * (2 - f)
    es = (f < 0 ? -1 : 1) * sqrt(abs(e2))
    e2m = (1 - e2)
    # c = sqrt( pow(1 + e, 1 + e) * pow(1 - e, 1 - e) ) )
    # See, for example, Lee (1976), p 100.
    c = sqrt(e2m) * exp(eatanhe(1.0, e2)) # Maybe revert to C++ version with es
    n = f / (2 - f)

    (1 - f) * exp(eatanhe(1.0, es))

    if !(isa(MaxPow,Int)) || MaxPow < 4 || MaxPow > 8
        error("MaxPow must be 4, 5, 6, 7 or 8")
    end

    # Stack-allocated constants
    if MaxPow == 4
        b1coeff = (1, 16, 64, 64)

        alpcoeff = (# alp[1]/n^1, polynomial in n of order 3
                    164, 225, -480, 360, 720,
                    # alp[2]/n^2, polynomial in n of order 2
                    557, -864, 390, 1440,
                    # alp[3]/n^3, polynomial in n of order 1
                    -1236, 427, 1680,
                    # alp[4]/n^4, polynomial in n of order 0
                    49561, 161280)

        betcoeff = (# bet[1]/n^1, polynomial in n of order 3
                    -4, 555, -960, 720, 1440,
                    # bet[2]/n^2, polynomial in n of order 2
                    -437, 96, 30, 1440,
                    # bet[3]/n^3, polynomial in n of order 1
                    -148, 119, 3360,
                    # bet[4]/n^4, polynomial in n of order 0
                    4397, 161280)
    elseif MaxPow == 5
        b1coeff = (1, 16, 64, 64)

        alpcoeff = (# alp[1]/n^1, polynomial in n of order 4
                    -635, 328, 450, -960, 720, 1440,
                    # alp[2]/n^2, polynomial in n of order 3
                    4496, 3899, -6048, 2730, 10080,
                    # alp[3]/n^3, polynomial in n of order 2
                    15061, -19776, 6832, 26880,
                    # alp[4]/n^4, polynomial in n of order 1
                    -171840, 49561, 161280,
                    # alp[5]/n^5, polynomial in n of order 0
                    34729, 80640)

        betcoeff = (# bet[1]/n^1, polynomial in n of order 4
                    -3645, -64, 8880, -15360, 11520, 23040,
                    # bet[2]/n^2, polynomial in n of order 3
                    4416, -3059, 672, 210, 10080,
                    # bet[3]/n^3, polynomial in n of order 2
                    -627, -592, 476, 13440,
                    # bet[4]/n^4, polynomial in n of order 1
                    -3520, 4397, 161280,
                    # bet[5]/n^5, polynomial in n of order 0
                    4583, 161280)
    elseif MaxPow == 6
        b1coeff = (1, 4, 64, 256, 256)

        alpcoeff = (# alp[1]/n^1, polynomial in n of order 5
                    31564, -66675, 34440, 47250, -100800, 75600, 151200,
                    # alp[2]/n^2, polynomial in n of order 4
                    -1983433, 863232, 748608, -1161216, 524160, 1935360,
                    # alp[3]/n^3, polynomial in n of order 3
                    670412, 406647, -533952, 184464, 725760,
                    # alp[4]/n^4, polynomial in n of order 2
                    6601661, -7732800, 2230245, 7257600,
                    # alp[5]/n^5, polynomial in n of order 1
                    -13675556, 3438171, 7983360,
                    # alp[6]/n^6, polynomial in n of order 0
                    212378941, 319334400)

        betcoeff = (# bet[1]/n^1, polynomial in n of order 5
                    384796, -382725, -6720, 932400, -1612800, 1209600, 2419200,
                    # bet[2]/n^2, polynomial in n of order 4
                    -1118711, 1695744, -1174656, 258048, 80640, 3870720,
                    # bet[3]/n^3, polynomial in n of order 3
                    22276, -16929, -15984, 12852, 362880,
                    # bet[4]/n^4, polynomial in n of order 2
                    -830251, -158400, 197865, 7257600,
                    # bet[5]/n^5, polynomial in n of order 1
                    -435388, 453717, 15966720,
                    # bet[6]/n^6, polynomial in n of order 0
                    20648693, 638668800)
    elseif MaxPow == 7
        b1coeff = (1, 4, 64, 256, 256)

        alpcoeff = (# alp[1]/n^1, polynomial in n of order 6
                    1804025, 2020096, -4267200, 2204160, 3024000, -6451200, 4838400, 9676800,
                    # alp[2]/n^2, polynomial in n of order 5
                    4626384, -9917165, 4316160, 3743040, -5806080, 2620800, 9676800,
                    # alp[3]/n^3, polynomial in n of order 4
                    -67102379, 26816480, 16265880, -21358080, 7378560, 29030400,
                    # alp[4]/n^4, polynomial in n of order 3
                    155912000, 72618271, -85060800, 24532695, 79833600,
                    # alp[5]/n^5, polynomial in n of order 2
                    102508609, -109404448, 27505368, 63866880,
                    # alp[6]/n^6, polynomial in n of order 1
                    -12282192400, 2760926233, 4151347200,
                    # alp[7]/n^7, polynomial in n of order 0
                    1522256789, 1383782400)

        betcoeff = (# bet[1]/n^1, polynomial in n of order 6
                    -5406467, 6156736, -6123600, -107520, 14918400, -25804800, 19353600,
                    38707200,
                    # bet[2]/n^2, polynomial in n of order 5
                    829456, -5593555, 8478720, -5873280, 1290240, 403200, 19353600,
                    # bet[3]/n^3, polynomial in n of order 4
                    9261899, 3564160, -2708640, -2557440, 2056320, 58060800,
                    # bet[4]/n^4, polynomial in n of order 3
                    14928352, -9132761, -1742400, 2176515, 79833600,
                    # bet[5]/n^5, polynomial in n of order 2
                    -8005831, -1741552, 1814868, 63866880,
                    # bet[6]/n^6, polynomial in n of order 1
                    -261810608, 268433009, 8302694400,
                    # bet[7]/n^7, polynomial in n of order 0
                    219941297, 5535129600)
    elseif MaxPow == 8
        b1coeff = (25, 64, 256, 4096, 16384, 16384)

        alpcoeff = (# alp[1]/n^1, polynomial in n of order 7
                    -75900428, 37884525, 42422016, -89611200, 46287360, 63504000, -135475200,
                    101606400, 203212800,
                    # alp[2]/n^2, polynomial in n of order 6
                    148003883, 83274912, -178508970, 77690880, 67374720, -104509440,
                    47174400, 174182400,
                    # alp[3]/n^3, polynomial in n of order 5
                    318729724, -738126169, 294981280, 178924680, -234938880, 81164160,
                    319334400,
                    # alp[4]/n^4, polynomial in n of order 4
                    -40176129013, 14967552000, 6971354016, -8165836800, 2355138720,
                    7664025600,
                    # alp[5]/n^5, polynomial in n of order 3
                    10421654396, 3997835751, -4266773472, 1072709352, 2490808320,
                    # alp[6]/n^6, polynomial in n of order 2
                    175214326799, -171950693600, 38652967262, 58118860800,
                    # alp[7]/n^7, polynomial in n of order 1
                    -67039739596, 13700311101, 12454041600,
                    # alp[8]/n^8, polynomial in n of order 0
                    1424729850961, 743921418240)

        betcoeff = (# bet[1]/n^1, polynomial in n of order 7
                    31777436, -37845269, 43097152, -42865200, -752640, 104428800, -180633600,
                    135475200, 270950400,
                    # bet[2]/n^2, polynomial in n of order 6
                    24749483, 14930208, -100683990, 152616960, -105719040, 23224320, 7257600,
                    348364800,
                    # bet[3]/n^3, polynomial in n of order 5
                    -232468668, 101880889, 39205760, -29795040, -28131840, 22619520,
                    638668800,
                    # bet[4]/n^4, polynomial in n of order 4
                    324154477, 1433121792, -876745056, -167270400, 208945440, 7664025600,
                    # bet[5]/n^5, polynomial in n of order 3
                    457888660, -312227409, -67920528, 70779852, 2490808320,
                    # bet[6]/n^6, polynomial in n of order 2
                    -19841813847, -3665348512, 3758062126, 116237721600,
                    # bet[7]/n^7, polynomial in n of order 1
                    -1989295244, 1979471673, 49816166400,
                    # bet[8]/n^8, polynomial in n of order 0
                    191773887257, 3719607091200)
    end

    m = div(MaxPow, 2)
    b1 = polyval(b1coeff[1:m+1], n*n) / (b1coeff[m+2] * (1 + n))

    # a1 is the equivalent radius for computing the circumference of ellipse.
    a1 = b1 * a
    o = 1
    d = n
    _alp = Vector{Float64}(undef, MaxPow)
    _bet = Vector{Float64}(undef, MaxPow)
    for l = 1:MaxPow
        m = MaxPow - l
        _alp[l] = d * polyval(alpcoeff[o:o+m], n) / alpcoeff[o + m + 1]
        _bet[l] = d * polyval(betcoeff[o:o+m], n) / betcoeff[o + m + 1]
        o += m + 2
        d *= n
    end
    # Post condition: o == sizeof(alpcoeff) / sizeof(real) &&
    # o == sizeof(betcoeff) / sizeof(real)
    alp = (_alp...,)
    bet = (_bet...,)

    return TransverseMercator{MaxPow}(a,f,e2,es,e2m,c,n,a1,b1,alp,bet)
end

const wgs84_tm = TransverseMercator(wgs84_ellipsoid)
const airy1830_tm = TransverseMercator(airy1830)
const clarke1866_tm = TransverseMercator(clarke1866)
const grs80_tm = TransverseMercator(grs80)

TransverseMercator(::WGS84) = wgs84_tm
TransverseMercator(::OSGB36) = airy1830_tm
TransverseMercator(::NAD27) = clarke1866_tm
TransverseMercator(::GDA94) = grs80_tm

# // Engsager and Poder (2007) use trigonometric series to convert between phi
# // and phip.  Here are the series...
# //
# // Conversion from phi to phip:
# //
# //     phip = phi + sum(c[j] * sin(2*j*phi), j, 1, 6)
# //
# //       c[1] = - 2 * n
# //              + 2/3 * n^2
# //              + 4/3 * n^3
# //              - 82/45 * n^4
# //              + 32/45 * n^5
# //              + 4642/4725 * n^6;
# //       c[2] =   5/3 * n^2
# //              - 16/15 * n^3
# //              - 13/9 * n^4
# //              + 904/315 * n^5
# //              - 1522/945 * n^6;
# //       c[3] = - 26/15 * n^3
# //              + 34/21 * n^4
# //              + 8/5 * n^5
# //              - 12686/2835 * n^6;
# //       c[4] =   1237/630 * n^4
# //              - 12/5 * n^5
# //              - 24832/14175 * n^6;
# //       c[5] = - 734/315 * n^5
# //              + 109598/31185 * n^6;
# //       c[6] =   444337/155925 * n^6;
# //
# // Conversion from phip to phi:
# //
# //     phi = phip + sum(d[j] * sin(2*j*phip), j, 1, 6)
# //
# //       d[1] =   2 * n
# //              - 2/3 * n^2
# //              - 2 * n^3
# //              + 116/45 * n^4
# //              + 26/45 * n^5
# //              - 2854/675 * n^6;
# //       d[2] =   7/3 * n^2
# //              - 8/5 * n^3
# //              - 227/45 * n^4
# //              + 2704/315 * n^5
# //              + 2323/945 * n^6;
# //       d[3] =   56/15 * n^3
# //              - 136/35 * n^4
# //              - 1262/105 * n^5
# //              + 73814/2835 * n^6;
# //       d[4] =   4279/630 * n^4
# //              - 332/35 * n^5
# //              - 399572/14175 * n^6;
# //       d[5] =   4174/315 * n^5
# //              - 144838/6237 * n^6;
# //       d[6] =   601676/22275 * n^6;
# //
# // In order to maintain sufficient relative accuracy close to the pole use
# //
# //     S = sum(c[i]*sin(2*i*phi),i,1,6)
# //     taup = (tau + tan(S)) / (1 - tau * tan(S))
#
# // In Math::taupf and Math::tauf we evaluate the forward transform explicitly
# // and solve the reverse one by Newton's method.
# //
# // There are adapted from TransverseMercatorExact (taup and taupinv).  tau =
# // tan(phi), taup = sinh(psi)

"""
    (x, y, γ, k) = transverse_mercator_forward(lon0, lat, lon, k0, tm::TransverseMercator)

Perform transverse-Mercator projection of `lat` and `lon` with respect to reference
meridian `lat0` and horizontal scaling `k0` (`= 0.9996` for UTM) using a
series expansion approach (see `TransverseMercator`). `γ` and `k` are the local
convergence and scaling factors, respectively.
"""
function transverse_mercator_forward(lon0, lat, lon, k0, tm::TransverseMercator{MaxPow}) where MaxPow
    lat = LatFix(lat)
    lon = AngDiff(lon0, lon)
    if abs(lon) > 50
        error("Relative longitude is too large for transverse-Mercator expansion, and must be less than 50° (received $(lon)°). Possibly used wrong UTM zone.")
    end
    # Explicitly enforce the parity
    latsign = lat < 0 ? -1 : 1
    lonsign = lon < 0 ? -1 : 1
    lon *= lonsign
    lat *= latsign
    backside = lon > 90
    if backside
        if (lat == 0)
            latsign = -1
        end
        lon = 180 - lon
    end

    sphi = sind(lat)
    cphi = cosd(lat)
    slam = sind(lon)
    clam = cosd(lon)

    # phi = latitude
    # phi' = conformal latitude
    # psi = isometric latitude
    # tau = tan(phi)
    # tau' = tan(phi')
    # [xi', eta'] = Gauss-Schreiber TM coordinates
    # [xi, eta] = Gauss-Krueger TM coordinates
    #
    # We use
    #   tan(phi') = sinh(psi)
    #   sin(phi') = tanh(psi)
    #   cos(phi') = sech(psi)
    #   denom^2    = 1-cos(phi')^2*sin(lam)^2 = 1-sech(psi)^2*sin(lam)^2
    #   sin(xip)   = sin(phi')/denom          = tanh(psi)/denom
    #   cos(xip)   = cos(phi')*cos(lam)/denom = sech(psi)*cos(lam)/denom
    #   cosh(etap) = 1/denom                  = 1/denom
    #   sinh(etap) = cos(phi')*sin(lam)/denom = sech(psi)*sin(lam)/denom
    if (lat != 90)
        tau = sphi / cphi
        taup = taupf(tau, tm.e2) # TODO maybe switch to the C++ version of taupf using tm.es?
        xip = atan(taup, clam)
        # Used to be
        #   etap = Math::atanh(sin(lam) / cosh(psi));
        etap = asinh(slam / hypot(taup, clam))
        # convergence and scale for Gauss-Schreiber TM (xip, etap) -- gamma0 =
        # atan(tan(xip) * tanh(etap)) = atan(tan(lam) * sin(phi'));
        # sin(phi') = tau'/sqrt(1 + tau'^2)
        # Krueger p 22 (44)
        gamma = atan2d(slam * taup, clam * hypot(1.0, taup))
        # k0 = sqrt(1 - _e2 * sin(phi)^2) * (cos(phi') / cos(phi)) * cosh(etap)
        # Note 1/cos(phi) = cosh(psip);
        # and cos(phi') * cosh(etap) = 1/hypot(sinh(psi), cos(lam))
        #
        # This form has cancelling errors.  This property is lost if cosh(psip)
        # is replaced by 1/cos(phi), even though it's using "primary" data (phi
        # instead of psip).
        k = sqrt(tm.e2m + tm.e2 * cphi*cphi) * hypot(1.0, tau) / hypot(taup, clam)
    else
        xip = pi/2
        etap = 0.0
        gamma = lon
        k = tm.c
    end

    # {xi',eta'} is {northing,easting} for Gauss-Schreiber transverse Mercator
    # (for eta' = 0, xi' = bet). {xi,eta} is {northing,easting} for transverse
    # Mercator with constant scale on the central meridian (for eta = 0, xip =
    # rectifying latitude).  Define
    #
    #   zeta = xi + i*eta
    #   zeta' = xi' + i*eta'
    #
    # The conversion from conformal to rectifying latitude can be expressed as
    # a series in _n:
    #
    #   zeta = zeta' + sum(h[j-1]' * sin(2 * j * zeta'), j = 1..maxpow_)
    #
    # where h[j]' = O(_n^j).  The reversion of this series gives
    #
    #   zeta' = zeta - sum(h[j-1] * sin(2 * j * zeta), j = 1..maxpow_)
    #
    # which is used in Reverse.
    #
    # Evaluate sums via Clenshaw method.  See
    #    http:#mathworld.wolfram.com/ClenshawRecurrenceFormula.html
    #
    # Let
    #
    #    S = sum(c[k] * F[k](x), k = 0..N)
    #    F[n+1](x) = alpha(n,x) * F[n](x) + beta(n,x) * F[n-1](x)
    #
    # Evaluate S with
    #
    #    y[N+2] = y[N+1] = 0
    #    y[k] = alpha(k,x) * y[k+1] + beta(k+1,x) * y[k+2] + c[k]
    #    S = c[0] * F[0](x) + y[1] * F[1](x) + beta(1,x) * F[0](x) * y[2]
    #
    # Here we have
    #
    #    x = 2 * zeta'
    #    F[n](x) = sin(n * x)
    #    a(n, x) = 2 * cos(x)
    #    b(n, x) = -1
    #    [ sin(A+B) - 2*cos(B)*sin(A) + sin(A-B) = 0, A = n*x, B = x ]
    #    N = maxpow_
    #    c[k] = _alp[k]
    #    S = y[1] * sin(x)
    #
    # For the derivative we have
    #
    #    x = 2 * zeta'
    #    F[n](x) = cos(n * x)
    #    a(n, x) = 2 * cos(x)
    #    b(n, x) = -1
    #    [ cos(A+B) - 2*cos(B)*cos(A) + cos(A-B) = 0, A = n*x, B = x ]
    #    c[0] = 1; c[k] = 2*k*_alp[k]
    #    S = (c[0] - y[2]) + y[1] * cos(x)
    c0 = cos(2 * xip)
    ch0 = cosh(2 * etap)
    s0 = sin(2 * xip)
    sh0 = sinh(2 * etap)
    ar = 2 * c0 * ch0
    ai = -2 * s0 * sh0 # 2 * cos(2*zeta')

    n = MaxPow

    xi0 = (isodd(n) ? tm.alp[n] : 0.0)
    eta0 = 0.0
    xi1 = 0.0
    eta1 = 0.0

    # Accumulators for dzeta/dzeta'
    if isodd(n)
        yr0 = 2 * MaxPow * tm.alp[n]
        n = n - 1
    else
        yr0 = 0.0
    end
    yi0 = 0.0
    yr1 = 0.0
    yi1 = 0.0
    while n > 0
        xi1  = ar * xi0 - ai * eta0 - xi1 + tm.alp[n]
        eta1 = ai * xi0 + ar * eta0 - eta1
        yr1 = ar * yr0 - ai * yi0 - yr1 + 2 * n * tm.alp[n]
        yi1 = ai * yr0 + ar * yi0 - yi1

        n = n - 1

        xi0  = ar * xi1 - ai * eta1 - xi0 + tm.alp[n]
        eta0 = ai * xi1 + ar * eta1 - eta0
        yr0 = ar * yr1 - ai * yi1 - yr0 + 2 * n * tm.alp[n]
        yi0 = ai * yr1 + ar * yi1 - yi0

        n = n - 1
    end

    ar *= 0.5
    ai *= 0.5      # cos(2*zeta')
    yr1 = 1 - yr1 + ar * yr0 - ai * yi0
    yi1 =   - yi1 + ai * yr0 + ar * yi0
    ar = s0 * ch0
    ai = c0 * sh0  # sin(2*zeta')

    xi  = xip  + ar * xi0 - ai * eta0
    eta = etap + ai * xi0 + ar * eta0

    # Fold in change in convergence and scale for Gauss-Schreiber TM to
    # Gauss-Krueger TM.

    gamma -= atan2d(yi1, yr1)
    k *= tm.b1 * hypot(yr1, yi1)
    y = tm.a1 * k0 * (backside ? pi - xi : xi) * latsign
    x = tm.a1 * k0 * eta * lonsign
    if backside
        gamma = 180 - gamma
    end
    gamma *= latsign * lonsign
    gamma = AngNormalize(gamma)
    k *= k0

    return (x, y, gamma, k)
end

"""
    (lat, lon, γ, k) = transverse_mercator_reverse(lon0, x, y, k0, tm::TransverseMercator)

Invert transverse-Mercator projection of `x` and `y` with respect to reference
meridian `lat0` and horizontal scaling `k0` (`= 0.9996` for UTM) using a
series expansion approach (see `TransverseMercator`). `γ` and `k` are the local
convergence and scaling factors, respectively.
"""
function transverse_mercator_reverse(lon0, x, y, k0, tm::TransverseMercator{MaxPow}) where MaxPow
    # This undoes the steps in transverse_mercator_forward.  The wrinkles are: (1) Use of the
    # reverted series to express zeta' in terms of zeta. (2) Newton's method
    # to solve for phi in terms of tan(phi).
    xi = y / (tm.a1 * k0)
    eta = x / (tm.a1 * k0)

    # Explicitly enforce the parity
    xisign = xi < 0 ? -1 : 1
    etasign = eta < 0 ? -1 : 1
    xi *= xisign
    eta *= etasign
    backside = xi > pi/2
    if backside
        xi = pi - xi
    end

    c0 = cos(2 * xi)
    ch0 = cosh(2 * eta)
    s0 = sin(2 * xi)
    sh0 = sinh(2 * eta)

    ar = 2 * c0 * ch0
    ai = -2 * s0 * sh0 # 2 * cos(2*zeta)
    n = MaxPow

    # Accumulators for zeta'
    xip0 = (isodd(n) ? -tm.bet[n] : 0.0)
    etap0 = 0.0
    xip1 = 0.0
    etap1 = 0.0

    # Accumulators for dzeta'/dzeta
    if isodd(n)
        yr0 =  2 * MaxPow * tm.bet[n]
        n -= 1
    else
        yr0 = 0.0
    end
    yi0 = 0.0
    yr1 = 0.0
    yi1 = 0.0

    while (n > 0)
        xip1  = ar * xip0 - ai * etap0 - xip1 - tm.bet[n]
        etap1 = ai * xip0 + ar * etap0 - etap1
        yr1 = ar * yr0 - ai * yi0 - yr1 - 2 * n * tm.bet[n]
        yi1 = ai * yr0 + ar * yi0 - yi1

        n = n - 1

        xip0  = ar * xip1 - ai * etap1 - xip0 - tm.bet[n]
        etap0 = ai * xip1 + ar * etap1 - etap0
        yr0 = ar * yr1 - ai * yi1 - yr0 - 2 * n * tm.bet[n]
        yi0 = ai * yr1 + ar * yi1 - yi0

        n = n - 1
    end

    ar *= 0.5
    ai *= 0.5  # cos(2*zeta')
    yr1 = 1 - yr1 + ar * yr0 - ai * yi0
    yi1 =   - yi1 + ai * yr0 + ar * yi0
    ar = s0 * ch0
    ai = c0 * sh0 # sin(2*zeta)

    xip  = xi  + ar * xip0 - ai * etap0
    etap = eta + ai * xip0 + ar * etap0

    # Convergence and scale for Gauss-Schreiber TM to Gauss-Krueger TM.
    gamma = atand(yi1, yr1)
    k = tm.b1 / hypot(yr1, yi1)
    # JHS 154 has
    #
    #   phi' = asin(sin(xi') / cosh(eta')) (Krueger p 17 (25))
    #   lam = asin(tanh(eta') / cos(phi')
    #   psi = asinh(tan(phi'))

    s = sinh(etap)
    c = max(0.0, cos(xip))  # cos(pi/2) might be negative
    r = hypot(s, c)

    if (r != 0)
        lon = atand(s, c)# Krueger p 17 (25)

        # Use Newton's method to solve for tau
        sxip = sin(xip)
        tau = tauf(sxip/r, tm.e2) # TODO maybe change to C++ es version
        gamma += atand(sxip * tanh(etap), c) # Krueger p 19 (31)
        lat = atand(tau)
        # Note cos(phi') * cosh(eta') = r
        k *= sqrt(tm.e2m + tm.e2 / (1 + tau*tau)) * hypot(1.0, tau) * r
    else
        lat = 90.0
        lon = 0.0
        k *= tm.c
    end

    lat *= xisign
    if backside
        lon = 180 - lon
    end
    lon *= etasign
    if abs(lon) > 50
        error("Relative longitude is too large for transverse-Mercator expansion, and must be less than 50° (calculated $(lon)°). Possibly used wrong UTM zone.")
    end
    lon = AngNormalize(lon + lon0)
    if backside
        gamma = 180 - gamma
    end
    gamma *= xisign * etasign
    gamma = AngNormalize(gamma)
    k *= k0

    return (lat, lon, gamma, k)
end


text:
for f in (:erf, :erfc)

    @eval function($f)(x::BigFloat, r::RoundingMode)
        setrounding(BigFloat, r) do
            ($f)(x)
        end
    end

    @eval ($f)(a::Interval{Float64}) = convert(Interval{Float64}, ($f)(big53(a)))

end

function erf(a::Interval{T}) where T
    isempty(a) && return a
    @round( erf(a.lo), erf(a.hi) )
end

function erfc(a::Interval{T}) where T
    isempty(a) && return a
    @round( erfc(a.hi), erfc(a.lo) )
end

erfinv(x, r::RoundingMode{:Down}) = _erfinv(x).lo
erfinv(x, r::RoundingMode{:Up}) = _erfinv(x).hi
erfcinv(x, r::RoundingMode{:Down}) = _erfcinv(x).hi
erfcinv(x, r::RoundingMode{:Up}) = _erfcinv(x).lo

erfinv(a::BigFloat) = mid(_erfinv(a))
erfcinv(a::BigFloat) = mid(_erfcinv(a))

function _erfinv(a::T) where T
    domain = Interval{T}(-1, 1)
    a ∉ domain && return DomainError("$a is not in [-1, 1]")
    f = x -> erf(x) - a
    fp = x->2/sqrt(pi_interval(T)) * exp(-x^2)
    rts = roots(f, Interval{T}(-Inf, Inf), 1e-20, deriv=fp)
    @assert length(rts) == 1 # && rts[1].status == :unique

    rts[1].interval
end

function erfinv(a::Interval{T}) where T
    domain = Interval{T}(-1, 1)
    a = a ∩ domain

    isempty(a) && return a
    @round(erfinv(a.lo), erfinv(a.hi))
end

function _erfcinv(a::T) where T
    domain = Interval{T}(0, 2)
    a ∉ domain && return DomainError("$a is not in [0, 2]")
    f = x -> erfc(x) - a
    fp = x -> -2/sqrt(pi_interval(T)) * exp(-x^2)
    rts = roots(f, Interval{T}(-Inf, Inf), 1e-20, deriv=fp)
    @assert length(rts) == 1 # && rts[1].status == :unique

    rts[1].interval
end

function erfcinv(a::Interval{T}) where T
    domain = Interval{T}(0, 2)
    a = a ∩ domain

    isempty(a) && return a
    @round(erfcinv(a.hi), erfcinv(a.lo))
end

ext: jl


text:
# Autogenerated wrapper script for PicoSAT_jll for i686-linux-musl
export libpicosat, picogcnf, picomcs, picomus, picosat

JLLWrappers.@generate_wrapper_header("PicoSAT")
JLLWrappers.@declare_library_product(libpicosat, "libpicosat.so")
JLLWrappers.@declare_executable_product(picogcnf)
JLLWrappers.@declare_executable_product(picomcs)
JLLWrappers.@declare_executable_product(picomus)
JLLWrappers.@declare_executable_product(picosat)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpicosat,
        "lib/libpicosat.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        picogcnf,
        "bin/picogcnf",
    )

    JLLWrappers.@init_executable_product(
        picomcs,
        "bin/picomcs",
    )

    JLLWrappers.@init_executable_product(
        picomus,
        "bin/picomus",
    )

    JLLWrappers.@init_executable_product(
        picosat,
        "bin/picosat",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    base dispatch for formal subset on taxonomies.

    a is formally a subset of A, even if they are identical,
    e.g. a = [1,2,3] and A = [1,2,3].

    returns true if a ⊆ A, false otherwise.
"""
function Base.:⊆(a::AbstractTaxonomy, A::AbstractTaxonomy)
    # how can this evaluate false?
    # 1. a is lower in rank than A, but a is not a clade of A
    # 2. a is higher rank than A, so function call doesn't make sense. throw warning
    # 3. a and A are not the same type of taxonomy, so function call doesn't make sense. throw error

    # how does this evaluate true?
    # 1. All of the species in a are also in A
    # 2. a is exactly the same as A
end

"""
    base dispatch for strict subset on taxonomies.

    a = [1,2,3] and A = [1,2,3] yields a ⊊ A = false,
    where a = [1,2] with the same A would yield a ⊊ A = true.

    returns true if a ⊊ A, false otherwise.
"""
function Base.:⊊(a::AbstractTaxonomy, A::AbstractTaxonomy)
    # how can this evaluate false?
    # 1. a is lower or equal to the rank than A, but a is not a clade of A
    # 2. a is higher rank or equal in rank to A, so function call doesn't make sense. throw warning
    # 3. a and A are not the same type of taxonomy, so function call doesn't make sense. throw error

    # how does this evaluate true?
    # 1. All of the species in a are also in A and a is a lower rank than A
end

"""
    Base.:∈(a::AbstractTaxonomy, A::AbstractTaxonomy)

    Evaluates same as ⊆(a,A) .
"""
Base.:∈(a::AbstractTaxonomy, A::AbstractTaxonomy) = Base.:⊆(a,A)
⊂(a::AbstractTaxonomy, A::AbstractTaxonomy) = Base.:⊆(a,A)

"""
    Base.:∩(A::AbstractTaxonomy, B::AbstractTaxonomy)

    Base dispatch for intersection of two taxonomies.

    Should return a taxonomy with the set of taxa
    that exist within both `A` and `B`.
"""
function Base.:∩(A::AbstractTaxonomy, B::AbstractTaxonomy)
end

"""
    Base.∪:(A::AbstractTaxonomy, B::AbstractTaxonomy)

    Base dispatch for union of two taxnomies.

    Should return a taxonomy with all of the species
    combined between `A` and `B`.

"""
function Base.:∪(A::AbstractTaxonomy, B::AbstractTaxonomy)

    # Should throw an error if :
    #   A and B are not of the same taxonomy type
    # Should throw a warning if :
    #   No common taxa found, so A and B merged via polytomy at the root
    # Otherwise:
    #   should find the highest order rank within A and B and root there.
end


text:
# Define symbol inputs for different variable types
const Infinite = :Infinite
const Point = :Point
const Hold = :Hold

## Extend Base.copy for new variable types
# GeneralVariableRef
Base.copy(v::GeneralVariableRef) = v

# InfiniteVariableRef
function Base.copy(v::InfiniteVariableRef,
                   new_model::InfiniteModel)::InfiniteVariableRef
    return InfiniteVariableRef(new_model, v.index)
end

# HoldVariableRef
function Base.copy(v::HoldVariableRef,
                   new_model::InfiniteModel)::HoldVariableRef
    return HoldVariableRef(new_model, v.index)
end

# PointVariableRef
function Base.copy(v::PointVariableRef,
                   new_model::InfiniteModel)::PointVariableRef
     return PointVariableRef(new_model, v.index)
 end

# Extend other Base functions
function Base.:(==)(v::T, w::U)::Bool where {T <: GeneralVariableRef,
                                             U <: GeneralVariableRef}
    return v.model === w.model && v.index == w.index && T == U
end
Base.broadcastable(v::GeneralVariableRef) = Ref(v)

# Extend JuMP functions
JuMP.isequal_canonical(v::GeneralVariableRef, w::GeneralVariableRef) = v == w
JuMP.variable_type(model::InfiniteModel) = GeneralVariableRef
function JuMP.variable_type(model::InfiniteModel, type::Symbol)
    if type == Infinite
        return InfiniteVariableRef
    elseif type == Point
        return PointVariableRef
    elseif type == Hold
        return HoldVariableRef
    elseif type == Parameter
        return ParameterRef
    else
        error("Invalid variable type.")
    end
end

# Check parameter tuple, ensure all elements contain parameter references
function _check_parameter_tuple(_error::Function, prefs::Tuple)
    types = [typeof(pref) for pref in prefs]
    num_params = length(types)
    valid_types = zeros(Bool, num_params)
    for i in eachindex(types)
        if types[i] == ParameterRef || types[i] <: AbstractArray{<:ParameterRef}
            valid_types[i] = true
        end
    end
    if !all(valid_types)
        _error("Invalid parameter type(s) given.")
    end
    return
end

# Convert parameter tuple s.t. array elements are SparseAxisArrays
function _make_formatted_tuple(prefs::Tuple)::Tuple
    converted_prefs = ()
    for pref in prefs
        if isa(pref, ParameterRef) || isa(pref, Number)
            converted_prefs = (converted_prefs..., pref)
        else
            converted_prefs = (converted_prefs...,
                               convert(JuMPC.SparseAxisArray, pref))
        end
    end
    return converted_prefs
end

# Ensure each element onyl contains parameters with same group ID
function _check_tuple_groups(_error::Function, prefs::Tuple)
    valid_elements = _only_one_group.(prefs)
    if sum(valid_elements) != length(prefs)
        _error("Each parameter tuple element must have contain only infinite " *
               "parameters with the same group ID.")
    end
    groups = _group.(prefs)
    if length(unique(groups)) != length(groups)
        _error("Cannot double specify infinite parameter references.")
    end
    return
end

# Ensure parameter values match shape of parameter reference tuple stored in the
# infinite variable reference
function _check_tuple_shape(_error::Function,
                            infinite_variable_ref::InfiniteVariableRef,
                            values::Tuple)
    prefs = parameter_refs(infinite_variable_ref)
    container = JuMPC.SparseAxisArray
    if length(prefs) != length(values)
        _error("The dimensions of the infinite parameter values must match " *
               "those defined for the infinite variable.")
    end
    for i in eachindex(values)
        if isa(prefs[i], ParameterRef) && !(isa(values[i], Number))
            _error("The dimensions and array type of the infinite parameter " *
                   "values must match those defined for the infinite variable.")
        elseif isa(prefs[i], container) && !isa(values[i], container)
            _error("The dimensions and array type of the infinite parameter " *
                   "values must match those defined for the infinite variable.")
        elseif isa(prefs[i], container)
            if keys(prefs[i].data) != keys(values[i].data)
                _error("Index keys of infinite parameter values don't match " *
                       "those defined for the infinite variable.")
            end
        end
    end
    return
end

# Used to ensure values don't violate parameter bounds
function _check_tuple_values(_error::Function, inf_vref::InfiniteVariableRef,
                             param_values::Tuple)
    prefs = parameter_refs(inf_vref)
    for i in eachindex(prefs)
        if isa(prefs[i], ParameterRef)
            if !supports_in_set(param_values[i], infinite_set(prefs[i]))
                _error("Parameter values violate parameter bounds.")
            end
        else
            for (k, pref) in prefs[i].data
                if !supports_in_set(param_values[i].data[k], infinite_set(pref))
                    _error("Parameter values violate parameter bounds.")
                end
            end
        end
    end
    return
end

# Update point variable info to consider the infinite variable
function _update_point_info(info::JuMP.VariableInfo, ivref::InfiniteVariableRef)
    if JuMP.has_lower_bound(ivref) && !info.has_fix && !info.has_lb
        info = JuMP.VariableInfo(true, JuMP.lower_bound(ivref),
                                 info.has_ub, info.upper_bound,
                                 info.has_fix, info.fixed_value,
                                 info.has_start, info.start,
                                 info.binary, info.integer)
    end
    if JuMP.has_upper_bound(ivref) && !info.has_fix && !info.has_ub
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 true, JuMP.upper_bound(ivref),
                                 info.has_fix, info.fixed_value,
                                 info.has_start, info.start,
                                 info.binary, info.integer)
    end
    if JuMP.is_fixed(ivref) && !info.has_fix  && !info.has_lb  && !info.has_ub
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 info.has_ub, info.upper_bound,
                                 true, JuMP.fix_value(ivref),
                                 info.has_start, info.start,
                                 info.binary, info.integer)
    end
    if !(JuMP.start_value(ivref) === NaN) && !info.has_start
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 info.has_ub, info.upper_bound,
                                 info.has_fix, info.fixed_value,
                                 true, JuMP.start_value(ivref),
                                 info.binary, info.integer)
    end
    if JuMP.is_binary(ivref) && !info.integer
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 info.has_ub, info.upper_bound,
                                 info.has_fix, info.fixed_value,
                                 info.has_start, info.start,
                                 true, info.integer)
    end
    if JuMP.is_integer(ivref) && !info.binary
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 info.has_ub, info.upper_bound,
                                 info.has_fix, info.fixed_value,
                                 info.has_start, info.start,
                                 info.binary, true)
    end
    return info
end

# Check that parameter_bounds argument is valid
function _check_bounds(bounds::ParameterBounds; _error = error)
    for (pref, set) in bounds.intervals
        # check that respects lower bound
        if JuMP.has_lower_bound(pref) && (set.lower_bound < JuMP.lower_bound(pref))
                _error("Specified parameter lower bound exceeds that defined " *
                       "for $pref.")
        end
        # check that respects upper bound
        if JuMP.has_upper_bound(pref) && (set.upper_bound > JuMP.upper_bound(pref))
                _error("Specified parameter upper bound exceeds that defined " *
                       "for $pref.")
        end
    end
    return
end

## Check to ensure correct inputs and build variables and return
# InfiniteVariable
function _make_variable(_error::Function, info::JuMP.VariableInfo, ::Val{Infinite};
                        parameter_refs::Union{ParameterRef,
                                              AbstractArray{<:ParameterRef},
                                              Tuple, Nothing} = nothing,
                        extra_kw_args...)::InfiniteVariable
    # check for unneeded keywords
    for (kwarg, _) in extra_kw_args
        _error("Keyword argument $kwarg is not for use with infinite variables.")
    end
    # check that we have been given parameter references
    if parameter_refs == nothing
        _error("Parameter references not specified, use the var(params...) " *
               "syntax or the parameter_refs keyword argument.")
    end
    # make sure the parameters are in a tuple
    if !isa(parameter_refs, Tuple)
        parameter_refs = (parameter_refs, )
    end
    # check tuple for validity and format
    _check_parameter_tuple(_error, parameter_refs)
    parameter_refs = _make_formatted_tuple(parameter_refs)
    _check_tuple_groups(_error, parameter_refs)
    # make the variable and return
    return InfiniteVariable(info, parameter_refs)
end

# PointVariable
function _make_variable(_error::Function, info::JuMP.VariableInfo, ::Val{Point};
                        infinite_variable_ref::Union{InfiniteVariableRef,
                                                     Nothing} = nothing,
                        parameter_values::Union{Number,
                                                AbstractArray{<:Number},
                                                Tuple, Nothing} = nothing,
                        extra_kw_args...)::PointVariable
    # check for unneeded keywords
    for (kwarg, _) in extra_kw_args
        _error("Keyword argument $kwarg is not for use with point variables.")
    end
    # ensure the needed arguments are given
    if parameter_values == nothing || infinite_variable_ref == nothing
        _error("Must specify the infinite variable and the values of its " *
               "infinite parameters")
    end
    # format as tuple if needed
    if !isa(parameter_values, Tuple)
        parameter_values = (parameter_values, )
    end
    # check information and prepare format
    parameter_values = _make_formatted_tuple(parameter_values)
    _check_tuple_shape(_error, infinite_variable_ref, parameter_values)
    _check_tuple_values(_error, infinite_variable_ref, parameter_values)
    info = _update_point_info(info, infinite_variable_ref)
    # make variable and return
    return PointVariable(info, infinite_variable_ref, parameter_values)
end

# HoldVariable
function _make_variable(_error::Function, info::JuMP.VariableInfo, ::Val{Hold};
                        parameter_bounds::ParameterBounds = ParameterBounds(),
                        extra_kw_args...)::HoldVariable
    # check for unneeded keywords
    for (kwarg, _) in extra_kw_args
        _error("Keyword argument $kwarg is not for use with hold variables.")
    end
    # check that the bounds don't violate parameter domains
    _check_bounds(parameter_bounds)
    # make variable and return
    return HoldVariable(info, parameter_bounds)
end

# Fallback method
function _make_variable(_error::Function, info::JuMP.VariableInfo, type;
                        extra_kw_args...)
    _error("Unrecognized variable type $type, should be Infinite, " *
           "Point, or Hold.")
end

"""
    JuMP.build_variable(_error::Function, info::JuMP.VariableInfo,
                        var_type::Symbol;
                        [parameter_refs::Union{ParameterRef,
                                              AbstractArray{<:ParameterRef},
                                              Tuple, Nothing} = nothing,
                        infinite_variable_ref::Union{InfiniteVariableRef,
                                                     Nothing} = nothing,
                        parameter_values::Union{Number, AbstractArray{<:Number},
                                                Tuple, Nothing} = nothing,
                        parameter_bounds::Union{Dict{ParameterRef, IntervalSet},
                                                Nothing} = nothing])

Extend the `JuMP.build_variable` function to accomodate `InfiniteOpt`
variable types. Returns the appropriate variable Datatype (i.e.,
[`InfiniteVariable`](@ref), [`PointVariable`](@ref), and
[`HoldVariable`](@ref)). Primarily this method is to be used internally by the
appropriate constructor macros [`@infinite_variable`](@ref),
[`@point_variable`](@ref), and [`@hold_variable`](@ref). However, it can be
called manually to build `InfiniteOpt` variables. Errors if an unneeded keyword
argument is given or if the keywoard arguments are formatted incorrectly (e.g.,
`parameter_refs` contains repeated parameter references when an infinite variable
is defined). Also errors if needed kewword arguments are negated.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel())
julia> @infinite_parameter(m, 0 <= t <= 1)
t

julia> info = VariableInfo(false, 0, false, 0, false, 0, false, 0, false, false);

julia> inf_var = build_variable(error, info, Infinite, parameter_refs = t)
InfiniteVariable{Int64,Int64,Int64,Int64}(VariableInfo{Int64,Int64,Int64,Int64}(false, 0, false, 0, false, 0, false, 0, false, false), (t,))

julia> ivref = add_variable(m, inf_var, "var_name")
var_name(t)

julia> pt_var = build_variable(error, info, Point, infinite_variable_ref = ivref,
                               parameter_values = 0.5)
PointVariable{Int64,Int64,Int64,Float64}(VariableInfo{Int64,Int64,Int64,Float64}(false, 0, false, 0, false, 0, true, 0.0, false, false), var_name(t), (0.5,))

julia> hd_var = build_variable(error, info, Hold)
HoldVariable{Int64,Int64,Int64,Int64}(VariableInfo{Int64,Int64,Int64,Int64}(false, 0, false, 0, false, 0, false, 0, false, false), Subdomain bounds (0): )
```
"""
function JuMP.build_variable(_error::Function, info::JuMP.VariableInfo,
                             var_type::Symbol;
                             macro_error::Union{Function, Nothing} = nothing,
                             kw_args...)
    if macro_error != nothing
        _error = macro_error # replace with macro error function
    end
    # make the variable and conduct necessary checks
    return _make_variable(_error, info, Val(var_type); kw_args...)
end

# Used to update the model.param_to_vars field
function _update_param_var_mapping(vref::InfiniteVariableRef, prefs::Tuple)
    model = JuMP.owner_model(vref)
    pref_list = _list_parameter_refs(prefs)
    for pref in pref_list
        if haskey(model.param_to_vars, JuMP.index(pref))
            push!(model.param_to_vars[JuMP.index(pref)], JuMP.index(vref))
        else
            model.param_to_vars[JuMP.index(pref)] = [JuMP.index(vref)]
        end
    end
    return
end

# check the pref tuple contains only valid parameters
function _check_parameters_valid(model::InfiniteModel, prefs::Tuple)
    pref_list = _list_parameter_refs(prefs)
    for pref in pref_list
        !JuMP.is_valid(model, pref) && error("Invalid Parameter reference " *
                                             "provided.")
    end
    return
end

# Used to add point variable support to parameter supports if necessary
function _update_param_supports(inf_vref::InfiniteVariableRef,
                                param_values::Tuple)
    prefs = parameter_refs(inf_vref)
    for i in eachindex(prefs)
        if isa(prefs[i], ParameterRef)
            add_supports(prefs[i], param_values[i])
        else
            for (k, v) in prefs[i].data
                add_supports(v, param_values[i].data[k])
            end
        end
    end
    return
end

# Used to update mapping infinite_to_points
function _update_infinite_point_mapping(pvref::PointVariableRef,
                                        ivref::InfiniteVariableRef)
    model = JuMP.owner_model(pvref)
    if haskey(model.infinite_to_points, JuMP.index(ivref))
        push!(model.infinite_to_points[JuMP.index(ivref)], JuMP.index(pvref))
    else
        model.infinite_to_points[JuMP.index(ivref)] = [JuMP.index(pvref)]
    end
    return
end

# Validate parameter bounds and add support(s) if needed
function _validate_bounds(model::InfiniteModel, bounds::ParameterBounds; _error = error)
    for (pref, set) in bounds.intervals
        # check validity
        !JuMP.is_valid(model, pref) && _error("Parameter bound reference " *
                                              "is invalid.")
        # ensure has a support if a point constraint was given
        if set.lower_bound == set.upper_bound
            add_supports(pref, set.lower_bound)
        end
    end
    return
end

## Make the variable reference and do checks/mapping updates
# InfiniteVariable
function _check_make_variable_ref(model::InfiniteModel,
                                  v::InfiniteVariable)::InfiniteVariableRef
    _check_parameters_valid(model, v.parameter_refs)
    vref = InfiniteVariableRef(model, model.next_var_index)
    _update_param_var_mapping(vref, v.parameter_refs)
    return vref
end

# PointVariable
function _check_make_variable_ref(model::InfiniteModel,
                                  v::PointVariable)::PointVariableRef
    ivref = v.infinite_variable_ref
    !JuMP.is_valid(model, ivref) && error("Invalid infinite variable " *
                                          "reference.")
    vref = PointVariableRef(model, model.next_var_index)
    _update_param_supports(ivref, v.parameter_values)
    _update_infinite_point_mapping(vref, ivref)
    return vref
end

# HoldVariable
function _check_make_variable_ref(model::InfiniteModel,
                                  v::HoldVariable)::HoldVariableRef
    _validate_bounds(model, v.parameter_bounds)
    vref = HoldVariableRef(model, model.next_var_index)
    if length(v.parameter_bounds.intervals) != 0
        model.has_hold_bounds = true
    end
    return vref
end

# Fallback
function _check_make_variable_ref(model::InfiniteModel, v)
    error("Invalid variable object type.")
end

"""
    JuMP.add_variable(model::InfiniteModel, var::InfOptVariable, [name::String = ""])

Extend the [`JuMP.add_variable`](@ref JuMP.add_variable(::JuMP.Model, ::JuMP.ScalarVariable, ::String))
function to accomodate `InfiniteOpt` variable types. Adds a variable to an
infinite model `model` and returns an appropriate variable reference (i.e.,
[`InfiniteVariableRef`](@ref), [`PointVariableRef`](@ref), or
[`HoldVariableRef`](@ref)). Primarily intended to be an internal function of the
constructor macros [`@infinite_variable`](@ref), [`@point_variable`](@ref), and
[`@hold_variable`](@ref). However, it can be used in combination with
[`JuMP.build_variable`](@ref) to add variables to an infinite model object.
Errors if invalid parameters reference(s) or an invalid infinite variable
reference is included in `var`.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel())
julia> @infinite_parameter(m, t in [0, 10]);

julia> info = VariableInfo(false, 0, false, 0, false, 0, false, 0, false, false);

julia> inf_var = build_variable(error, info, Infinite, parameter_refs = t);

julia> ivref = add_variable(m, inf_var, "var_name")
var_name(t)

julia> pt_var = build_variable(error, info, Point, infinite_variable_ref = ivref,
                               parameter_values = 0.5);

julia> pvref = add_variable(m, pt_var, "var_alias")
var_alias

julia> hd_var = build_variable(error, info, Hold);

julia> hvref = add_variable(m, hd_var, "var_name")
var_name
```
"""
function JuMP.add_variable(model::InfiniteModel, var::InfOptVariable,
                           name::String = "")
    model.next_var_index += 1
    vref = _check_make_variable_ref(model, var)
    model.vars[JuMP.index(vref)] = var
    JuMP.set_name(vref, name)
    if var.info.has_lb
        newset = MOI.GreaterThan(convert(Float64, var.info.lower_bound))
        cref = JuMP.add_constraint(JuMP.owner_model(vref),
                                   JuMP.ScalarConstraint(vref, newset))
        _set_lower_bound_index(vref, JuMP.index(cref))
        model.constr_in_var_info[JuMP.index(cref)] = true
    end
    if var.info.has_ub
        newset = MOI.LessThan(convert(Float64, var.info.upper_bound))
        cref = JuMP.add_constraint(JuMP.owner_model(vref),
                                   JuMP.ScalarConstraint(vref, newset))
        _set_upper_bound_index(vref, JuMP.index(cref))
        model.constr_in_var_info[JuMP.index(cref)] = true
    end
    if var.info.has_fix
        newset = MOI.EqualTo(convert(Float64, var.info.fixed_value))
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(vref, newset))
        _set_fix_index(vref, JuMP.index(cref))
        model.constr_in_var_info[JuMP.index(cref)] = true
    end
    if var.info.binary
        cref = JuMP.add_constraint(JuMP.owner_model(vref),
                                   JuMP.ScalarConstraint(vref, MOI.ZeroOne()))
        _set_binary_index(vref, JuMP.index(cref))
        model.constr_in_var_info[JuMP.index(cref)] = true
    elseif var.info.integer
        cref = JuMP.add_constraint(JuMP.owner_model(vref),
                                   JuMP.ScalarConstraint(vref, MOI.Integer()))
        _set_integer_index(vref, JuMP.index(cref))
        model.constr_in_var_info[JuMP.index(cref)] = true
    end
    model.var_in_objective[JuMP.index(vref)] = false
    return vref
end

"""
    JuMP.owner_model(vref::GeneralVariableRef)::InfiniteModel

Extend [`JuMP.owner_model`](@ref JuMP.owner_model(::JuMP.AbstractVariableRef)) function
for `InfiniteOpt` variables. Returns the infinite model associated with `vref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 <= vref <= 1))
julia> owner_model(vref)
An InfiniteOpt Model
Feasibility problem with:
Variable: 1
`HoldVariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 1 constraint
`HoldVariableRef`-in-`MathOptInterface.LessThan{Float64}`: 1 constraint
Names registered in the model: vref
Optimizer model backend information:
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
```
"""
JuMP.owner_model(vref::GeneralVariableRef)::InfiniteModel = vref.model

"""
    JuMP.index(v::GeneralVariableRef)::Int

Extent [`JuMP.index`](@ref JuMP.index(::JuMP.VariableRef)) to return the index of a
`InfiniteOpt` variable.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref))
julia> index(vref)
1
```
"""
JuMP.index(v::GeneralVariableRef)::Int = v.index

"""
    used_by_constraint(vref::InfOptVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a constraint.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref))
julia> used_by_constraint(vref)
false
```
"""
function used_by_constraint(vref::InfOptVariableRef)::Bool
    return haskey(JuMP.owner_model(vref).var_to_constrs, JuMP.index(vref))
end

"""
    used_by_measure(vref::InfOptVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a measure.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref); m.var_to_meas[1] = [1])
julia> used_by_measure(vref)
true
```
"""
function used_by_measure(vref::InfOptVariableRef)::Bool
    return haskey(JuMP.owner_model(vref).var_to_meas, JuMP.index(vref))
end

"""
    used_by_objective(vref::InfOptVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by the objective.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref); m.var_in_objective[1] = true)
julia> used_by_objective(vref)
true
```
"""
function used_by_objective(vref::InfOptVariableRef)::Bool
    return JuMP.owner_model(vref).var_in_objective[JuMP.index(vref)]
end

"""
    is_used(vref::InfOptVariableRef)::Bool

Return a `Bool` indicating if `vref` is used in the model.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, 0 <= vref))
julia> is_used(vref)
true
```
"""
function is_used(vref::InfOptVariableRef)::Bool
    return used_by_measure(vref) || used_by_constraint(vref) || used_by_objective(vref)
end

"""
    used_by_point_variable(vref::InfiniteVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a point variable.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @infinite_variable(m, vref(@infinite_parameter(m, t in [0, 1]))))
julia> used_by_point_variable(vref)
false
```
"""
function used_by_point_variable(vref::InfiniteVariableRef)::Bool
    return haskey(JuMP.owner_model(vref).infinite_to_points, JuMP.index(vref))
end

"""
    used_by_reduced_variable(vref::InfiniteVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a reduced infinite variable.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @infinite_variable(m, vref(@infinite_parameter(m, t in [0, 1]))))
julia> used_by_reduced_variable(vref)
false
```
"""
function used_by_reduced_variable(vref::InfiniteVariableRef)::Bool
    return haskey(JuMP.owner_model(vref).infinite_to_reduced, JuMP.index(vref))
end

"""
    is_used(vref::InfiniteVariableRef)::Bool

Return a `Bool` indicating if `vref` is used in the model.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @infinite_variable(m, vref(@infinite_parameter(m, t in [0, 1]))))
julia> is_used(vref)
false
```
"""
function is_used(vref::InfiniteVariableRef)::Bool
    if used_by_measure(vref) || used_by_constraint(vref)
        return true
    end
    if used_by_point_variable(vref)
        for vindex in JuMP.owner_model(vref).infinite_to_points[JuMP.index(vref)]
            if is_used(PointVariableRef(JuMP.owner_model(vref), vindex))
                return true
            end
        end
    end
    if used_by_reduced_variable(vref)
        for rindex in JuMP.owner_model(vref).infinite_to_reduced[JuMP.index(vref)]
            rvref = ReducedInfiniteVariableRef(JuMP.owner_model(vref), rindex)
            if used_by_constraint(rvref) || used_by_measure(rvref)
                return true
            end
        end
    end
    return false
end

"""
    JuMP.delete(model::InfiniteModel, vref::InfOptVariableRef)

Extend [`JuMP.delete`](@ref JuMP.delete(::JuMP.Model, ::JuMP.VariableRef)) to delete
`InfiniteOpt` variables and their dependencies. Errors if variable is invalid,
meaning it has already been deleted or it belongs to another model.

**Example**
```julia-repl
julia> print(model)
Min measure(g(t)*t) + z
Subject to
 z ≥ 0.0
 g(t) + z ≥ 42.0
 g(0.5) = 0
 t ∈ [0, 6]

julia> delete(model, g)

julia> print(model)
Min measure(t) + z
Subject to
 z ≥ 0.0
 z ≥ 42.0
 t ∈ [0, 6]
```
"""
function JuMP.delete(model::InfiniteModel, vref::InfOptVariableRef)
    @assert JuMP.is_valid(model, vref) "Variable is invalid."
    # update the optimizer model status
    if is_used(vref)
        set_optimizer_model_ready(model, false)
    end
    # remove variable info constraints associated with vref
    if JuMP.has_lower_bound(vref)
        JuMP.delete_lower_bound(vref)
    end
    if JuMP.has_upper_bound(vref)
        JuMP.delete_upper_bound(vref)
    end
    if JuMP.is_fixed(vref)
        JuMP.unfix(vref)
    end
    if JuMP.is_binary(vref)
        JuMP.unset_binary(vref)
    elseif JuMP.is_integer(vref)
        JuMP.unset_integer(vref)
    end
    # remove dependencies from measures and update them
    if used_by_measure(vref)
        for mindex in model.var_to_meas[JuMP.index(vref)]
            if isa(model.measures[mindex].func, InfOptVariableRef)
                model.measures[mindex] = Measure(zero(JuMP.AffExpr),
                                                 model.measures[mindex].data)
            else
                _remove_variable(model.measures[mindex].func, vref)
            end
            JuMP.set_name(MeasureRef(model, mindex),
                           _make_meas_name(model.measures[mindex]))
        end
        # delete mapping
        delete!(model.var_to_meas, JuMP.index(vref))
    end
    # remove dependencies from measures and update them
    if used_by_constraint(vref)
        for cindex in model.var_to_constrs[JuMP.index(vref)]
            if isa(model.constrs[cindex].func, InfOptVariableRef)
                model.constrs[cindex] = JuMP.ScalarConstraint(zero(JuMP.AffExpr),
                                                      model.constrs[cindex].set)
            else
                _remove_variable(model.constrs[cindex].func, vref)
            end
        end
        # delete mapping
        delete!(model.var_to_constrs, JuMP.index(vref))
    end
    # remove from objective if vref is in it
    if used_by_objective(vref)
        if isa(model.objective_function, InfOptVariableRef)
            model.objective_function = zero(JuMP.AffExpr)
        else
            _remove_variable(model.objective_function, vref)
        end
    end
    # do specific updates if vref is infinite
    if isa(vref, InfiniteVariableRef)
        # update parameter mapping
        all_prefs = _list_parameter_refs(parameter_refs(vref))
        for pref in all_prefs
            filter!(e -> e != JuMP.index(vref),
                    model.param_to_vars[JuMP.index(pref)])
            if length(model.param_to_vars[JuMP.index(pref)]) == 0
                delete!(model.param_to_vars, JuMP.index(pref))
            end
        end
        # delete associated point variables and mapping
        if used_by_point_variable(vref)
            for index in model.infinite_to_points[JuMP.index(vref)]
                JuMP.delete(model, PointVariableRef(model, index))
            end
            delete!(model.infinite_to_points, JuMP.index(vref))
        end
        # delete associated reduced variables and mapping
        if used_by_reduced_variable(vref)
            for index in model.infinite_to_reduced[JuMP.index(vref)]
                JuMP.delete(model, ReducedInfiniteVariableRef(model, index))
            end
            delete!(model.infinite_to_reduced, JuMP.index(vref))
        end
    end
    # update mappings if is point variable
    if isa(vref, PointVariableRef)
        ivref = infinite_variable_ref(vref)
        filter!(e -> e != JuMP.index(vref),
                model.infinite_to_points[JuMP.index(ivref)])
        if length(model.infinite_to_points[JuMP.index(ivref)]) == 0
            delete!(model.infinite_to_points, JuMP.index(ivref))
        end
    end
    # delete the variable information
    delete!(model.var_in_objective, JuMP.index(vref))
    delete!(model.vars, JuMP.index(vref))
    delete!(model.var_to_name, JuMP.index(vref))
    return
end

"""
    JuMP.is_valid(model::InfiniteModel, vref::InfOptVariableRef)::Bool

Extend [`JuMP.is_valid`](@ref JuMP.is_valid(::JuMP.Model, ::JuMP.VariableRef))
to accomodate `InfiniteOpt` variables.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @hold_variable(model, vref))
julia> is_valid(model, vref)
true
```
"""
function JuMP.is_valid(model::InfiniteModel, vref::InfOptVariableRef)::Bool
    return (model === JuMP.owner_model(vref) && JuMP.index(vref) in keys(model.vars))
end

"""
    JuMP.num_variables(model::InfiniteModel)::Int

Extend [`JuMP.num_variables`](@ref JuMP.num_variables(::JuMP.Model)) to return the
number of `InfiniteOpt` variables assigned to `model`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @hold_variable(model, vref[1:3]))
julia> num_variables(model)
3
```
"""
JuMP.num_variables(model::InfiniteModel)::Int = length(model.vars)

# Include all the extension functions for manipulating the properties associated
# with VariableInfo
include("variable_info.jl")

"""
    JuMP.name(vref::InfOptVariableRef)::String

Extend [`JuMP.name`](@ref JuMP.name(::JuMP.VariableRef)) to return the names of
`InfiniteOpt` variables.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref, base_name = "var_name"))
julia> name(vref)
"var_name"
```
"""
function JuMP.name(vref::InfOptVariableRef)::String
    return JuMP.owner_model(vref).var_to_name[JuMP.index(vref)]
end

"""
    JuMP.set_name(vref::HoldVariableRef, name::String)

Extend [`JuMP.set_name`](@ref JuMP.set_name(::JuMP.VariableRef, ::String)) to set
names of hold variables.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, vref))
julia> set_name(vref, "var_name")

julia> name(vref)
"var_name"
```
"""
function JuMP.set_name(vref::HoldVariableRef, name::String)
    JuMP.owner_model(vref).var_to_name[JuMP.index(vref)] = name
    JuMP.owner_model(vref).name_to_var = nothing
    return
end

"""
    infinite_variable_ref(vref::PointVariableRef)::InfiniteVariableRef

Return the `InfiniteVariableRef` associated with the point variable `vref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> @infinite_variable(model, T(t))
T(t)

julia> vref = @point_variable(model, T(0))
T(0)

julia> infinite_variable_ref(vref)
T(t)
```
"""
function infinite_variable_ref(vref::PointVariableRef)::InfiniteVariableRef
    return JuMP.owner_model(vref).vars[JuMP.index(vref)].infinite_variable_ref
end

"""
    parameter_values(vref::PointVariableRef)::Tuple

Return the support point associated with the point variable `vref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> @infinite_variable(model, T(t))
T(t)

julia> vref = @point_variable(model, T(0))
T(0)

julia> parameter_values(vref)
(0,)
```
"""
function parameter_values(vref::PointVariableRef)::Tuple
    return JuMP.owner_model(vref).vars[JuMP.index(vref)].parameter_values
end

# Internal function used to change the parameter value tuple of a point variable
function _update_variable_param_values(vref::PointVariableRef, pref_vals::Tuple)
    info = JuMP.owner_model(vref).vars[JuMP.index(vref)].info
    ivref = JuMP.owner_model(vref).vars[JuMP.index(vref)].infinite_variable_ref
    JuMP.owner_model(vref).vars[JuMP.index(vref)] = PointVariable(info, ivref,
                                                                  pref_vals)
    return
end

# Get root name of infinite variable
function _root_name(vref::InfiniteVariableRef)
    name = JuMP.name(vref)
    return name[1:findfirst(isequal('('), name)-1]
end

## Return the parameter value as an appropriate string
# Number
function _make_str_value(value::Number)::String
    return string(JuMP._string_round(value))
end

# AbstractArray
function _make_str_value(value::AbstractArray)::String
    if length(keys(value)) <= 4
        str_value = "["
        counter = 1
        for key in sort(collect(keys(value)))
            if counter != length(keys(value))
                str_value *= JuMP._string_round(value[key]) * ", "
            else
                str_value *= JuMP._string_round(value[key]) * "]"
            end
            counter += 1
        end
        return string(str_value)
    else
        value = [value[key] for key in sort(collect(keys(value)))]
        return string("[", JuMP._string_round(first(value)), ", ..., ",
                      JuMP._string_round(last(value)), "]")
    end
end

"""
    JuMP.set_name(vref::PointVariableRef, name::String)

Extend [`JuMP.set_name`](@ref JuMP.set_name(::JuMP.VariableRef, ::String)) to set
the names of point variables.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> @infinite_variable(model, T(t))
T(t)

julia> vref = @point_variable(model, T(0))
T(0)

julia> set_name(vref, "new_name")

julia> name(vref)
"new_name"
```
"""
function JuMP.set_name(vref::PointVariableRef, name::String)
    if length(name) == 0
        inf_var_ref = infinite_variable_ref(vref::PointVariableRef)
        name = _root_name(inf_var_ref)
        values = JuMP.owner_model(vref).vars[JuMP.index(vref)].parameter_values
        name = string(name, "(")
        for i in eachindex(values)
            if i != length(values)
                name *= _make_str_value(values[i]) * ", "
            else
                name *= _make_str_value(values[i]) * ")"
            end
        end
    end
    JuMP.owner_model(vref).var_to_name[JuMP.index(vref)] = name
    JuMP.owner_model(vref).name_to_var = nothing
    return
end

"""
    parameter_refs(vref::InfiniteVariableRef)::Tuple

Return the `ParameterRef`(s) associated with the infinite variable `vref`. This
is formatted as a Tuple of containing the parameter references as they inputted
to define `vref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> @infinite_variable(model, T(t))
T(t)

julia> parameter_refs(T)
(t,)
```
"""
function parameter_refs(vref::InfiniteVariableRef)
    return JuMP.owner_model(vref).vars[JuMP.index(vref)].parameter_refs
end

# Internal function used to change the parameter reference tuple of an infinite
# variable
function _update_variable_param_refs(vref::InfiniteVariableRef, prefs::Tuple)
    info = JuMP.owner_model(vref).vars[JuMP.index(vref)].info
    JuMP.owner_model(vref).vars[JuMP.index(vref)] = InfiniteVariable(info, prefs)
    return
end

"""
    set_parameter_refs(vref::InfiniteVariableRef, prefs::Tuple)

Specify a new parameter reference tuple `prefs` for the infinite variable `vref`.
Note each element must contain a single parameter reference or an array of
parameter references. Errors if a parameter is double specified or if an element
contains parameters with different group IDs.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> @infinite_variable(model, T(t))
T(t)

julia> @infinite_parameter(model, x[1:2] in [-1, 1])
2-element Array{ParameterRef,1}:
 x[1]
 x[2]

julia> set_parameter_refs(T, (t, x))

julia> parameter_refs(T)
(t,   [2]  =  x[2]
  [1]  =  x[1])
```
"""
function set_parameter_refs(vref::InfiniteVariableRef, prefs::Tuple)
    _check_parameter_tuple(error, prefs)
    prefs = _make_formatted_tuple(prefs)
    _check_tuple_groups(error, prefs)
    _update_variable_param_refs(vref, prefs)
    JuMP.set_name(vref, _root_name(vref))
    if is_used(vref)
        set_optimizer_model_ready(JuMP.owner_model(vref), false)
    end
    return
end

"""
    add_parameter_ref(vref::InfiniteVariableRef,
                      pref::Union{ParameterRef, AbstractArray{<:ParameterRef}})

Add additional parameter reference or group of parameter references to be
associated with the infinite variable `vref`. Errors if the parameter references
are already added to the variable or if the added parameters have different
group IDs.

```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> @infinite_variable(model, T(t))
T(t)

julia> @infinite_parameter(model, x[1:2] in [-1, 1])
2-element Array{ParameterRef,1}:
 x[1]
 x[2]

julia> add_parameter_ref(T, x)

julia> name(T)
"T(t, x)"
```
"""
function add_parameter_ref(vref::InfiniteVariableRef,
                       pref::Union{ParameterRef, AbstractArray{<:ParameterRef}})
    set_parameter_refs(vref, (parameter_refs(vref)..., pref))
    return
end

"""
    parameter_bounds(vref::HoldVariableRef)::ParameterBounds

Return the [`ParameterBounds`](@ref) object associated with the hold variable
`vref`. It contains a dictionary where each key is a `ParameterRef` which points
to an `IntervalSet` that that defines a sub-domain for `vref` relative to that
parameter reference.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10])
t

julia> @hold_variable(model, vref, parameter_bounds = (t in [0, 2]))
vref

julia> parameter_bounds(vref)
Subdomain bounds (1): t ∈ [0, 2]
```
"""
function parameter_bounds(vref::HoldVariableRef)::ParameterBounds
    return JuMP.owner_model(vref).vars[JuMP.index(vref)].parameter_bounds
end

"""
    has_parameter_bounds(vref::HoldVariableRef)::Bool

Return a `Bool` indicating if `vref` is limited to a sub-domain as defined
by parameter bound.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10])
t

julia> @hold_variable(model, vref, parameter_bounds = (t in [0, 2]))
vref

julia> has_parameter_bounds(vref)
true
```
"""
function has_parameter_bounds(vref::HoldVariableRef)::Bool
    return length(parameter_bounds(vref)) != 0
end

# Other variable types
function has_parameter_bounds(vref::GeneralVariableRef)::Bool
    return false
end

# Internal function used to change the parameter bounds of a hold variable
function _update_variable_param_bounds(vref::HoldVariableRef,
                                       bounds::ParameterBounds)
    info = JuMP.owner_model(vref).vars[JuMP.index(vref)].info
    JuMP.owner_model(vref).vars[JuMP.index(vref)] = HoldVariable(info, bounds)
    return
end

## Check that the bounds dictionary is compadable with existing dependent measures
function _check_meas_bounds(bounds::ParameterBounds, data::AbstractMeasureData;
                            _error = error)
    if !measure_data_in_hold_bounds(data, bounds)
        _error("New bounds don't span existing dependent measure bounds.")
    end
    return
end

# Update the current bounds to overlap with the new bounds if possible
function _update_bounds(bounds1::Dict, bounds2::Dict; _error = error)
    # check each new bound
    for (pref, set) in bounds2
        # we have a new bound
        if !haskey(bounds1, pref)
            bounds1[pref] = set
        # the previous set and the new one do not overlap
        elseif set.lower_bound > bounds1[pref].upper_bound || set.upper_bound < bounds1[pref].lower_bound
            _error("Sub-domains of constraint and/or hold variable(s) do not" *
                   " overlap. Consider changing the parameter bounds of the" *
                   " constraint and/or hold variable(s).")
        # we have an existing bound
        else
            # we have a new stricter lower bound to update with
            if set.lower_bound > bounds1[pref].lower_bound
                bounds1[pref] = IntervalSet(set.lower_bound, bounds1[pref].upper_bound)
            end
            # we have a new stricter upper bound to update with
            if set.upper_bound < bounds1[pref].upper_bound
                bounds1[pref] = IntervalSet(bounds1[pref].lower_bound, set.upper_bound)
            end
        end
    end
    return
end

## Update the variable bounds if it has any
# GeneralVariableRef
function _update_var_bounds(vref::GeneralVariableRef,
                            constr_bounds::ParameterBounds)
    return
end

# HoldVariableRef
function _update_var_bounds(vref::HoldVariableRef,
                            constr_bounds::ParameterBounds)
    if has_parameter_bounds(vref)
        _update_bounds(constr_bounds.intervals, parameter_bounds(vref).intervals)
    end
    return
end

# MeasureRef
function _update_var_bounds(mref::MeasureRef,
                            constr_bounds::ParameterBounds)
    vrefs = _all_function_variables(measure_function(mref))
    for vref in vrefs
        _update_var_bounds(vref, constr_bounds)
    end
    return
end

## Rebuild the constraint bounds (don't change in case of error)
# BoundedScalarConstraint
function _rebuild_constr_bounds(c::BoundedScalarConstraint,
                                var_bounds::ParameterBounds; _error = error)
    # prepare new constraint
    vrefs = _all_function_variables(c.func)
    c_new = BoundedScalarConstraint(c.func, c.set, copy(c.orig_bounds), c.orig_bounds)
    # look for bounded hold variables and update bounds
    for vref in vrefs
        _update_var_bounds(vref, c_new.bounds)
    end
    # check if the constraint bounds have and update if doesn't
    if length(c_new.bounds) == 0
        c_new = JuMP.ScalarConstraint(c.func, c.set)
    end
    return c_new
end

# ScalarConstraint
function _rebuild_constr_bounds(c::JuMP.ScalarConstraint,
                                var_bounds::ParameterBounds; _error = error)
    return BoundedScalarConstraint(c.func, c.set, var_bounds, ParameterBounds())
end

"""
    set_parameter_bounds(vref::HoldVariableRef, bounds::ParameterBounds;
                         [force = false])

Specify a new dictionary of parameter bounds `bounds` for the hold variable `vref`.
These are stored in a [`ParameterBounds`](@ref) object which contains a dictionary.
Note the dictionary keys must be `ParameterRef`s and the values must be
`IntervalSet`s that indicate a particular sub-domain for which `vref` is defined.
This is meant to be primarily used by [`@set_parameter_bounds`](@ref) which
provides a more intuitive syntax.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10])
t

julia> @hold_variable(model, vref)
vref

julia> set_parameter_bounds(vref, ParameterBounds(Dict(t => IntervalSet(0, 2))))

julia> parameter_bounds(vref)
Subdomain bounds (1): t ∈ [0, 2]
```
"""
function set_parameter_bounds(vref::HoldVariableRef, bounds::ParameterBounds;
                              force = false, _error = error)
    if has_parameter_bounds(vref) && !force
        _error("$vref already has parameter bounds. Consider adding more using " *
               "`add_parameter_bounds` or overwriting them by setting " *
               "the keyword argument `force = true`")
    else
        # check that bounds are valid and add support(s) if necessary
        _check_bounds(bounds, _error = _error)
        # check dependent measures
        meas_cindices = []
        if used_by_measure(vref)
            for mindex in JuMP.owner_model(vref).var_to_meas[JuMP.index(vref)]
                meas = JuMP.owner_model(vref).measures[mindex]
                _check_meas_bounds(bounds, meas.data, _error = _error)
                if used_by_constraint(MeasureRef(JuMP.owner_model(vref), mindex))
                    indices = JuMP.owner_model(vref).meas_to_constrs[mindex]
                    meas_cindices = [meas_cindices; indices]
                end
            end
        end
        # set the new bounds
        _validate_bounds(JuMP.owner_model(vref), bounds, _error = _error)
        _update_variable_param_bounds(vref, bounds)
        # check and update dependent constraints
        if used_by_constraint(vref) || length(meas_cindices) != 0
            for cindex in unique([meas_cindices; JuMP.owner_model(vref).var_to_constrs[JuMP.index(vref)]])
                constr = JuMP.owner_model(vref).constrs[cindex]
                new_constr = _rebuild_constr_bounds(constr, bounds,
                                                    _error = _error)
                JuMP.owner_model(vref).constrs[cindex] = new_constr
            end
        end
        # update status
        JuMP.owner_model(vref).has_hold_bounds = true
        if is_used(vref)
            set_optimizer_model_ready(JuMP.owner_model(vref), false)
        end
    end
    return
end

## Check and update the constraint bounds (don't change in case of error)
# BoundedScalarConstraint
function _update_constr_bounds(bounds::ParameterBounds, c::BoundedScalarConstraint;
                               _error = error)
    new_bounds_dict = copy(c.bounds.intervals)
    _update_bounds(new_bounds_dict, bounds.intervals, _error = _error)
    return BoundedScalarConstraint(c.func, c.set, ParameterBounds(new_bounds_dict),
                                   c.orig_bounds)
end

# ScalarConstraint
function _update_constr_bounds(bounds::ParameterBounds, c::JuMP.ScalarConstraint;
                               _error = error)
    return BoundedScalarConstraint(c.func, c.set, bounds, ParameterBounds())
end

"""
    add_parameter_bound(vref::HoldVariableRef, pref::ParameterRef,
                        lower::Number, upper::Number)

Add an additional parameter bound to `vref` such that it is defined over the
sub-domain based on `pref` from `lower` to `upper`. This is primarily meant to be
used by [`@add_parameter_bounds`](@ref).

```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10])
t

julia> @hold_variable(model, vref)
vref

julia> add_parameter_bound(vref, t, 0, 2)

julia> parameter_bounds(vref)
Subdomain bounds (1): t ∈ [0, 2]
```
"""
function add_parameter_bound(vref::HoldVariableRef, pref::ParameterRef,
                             lower::Number, upper::Number; _error = error)
    # check the new bounds
    new_bounds = ParameterBounds(Dict(pref => IntervalSet(lower, upper)))
    _check_bounds(new_bounds, _error = _error)
    # check dependent measures
    meas_cindices = []
    if used_by_measure(vref)
        for mindex in JuMP.owner_model(vref).var_to_meas[JuMP.index(vref)]
            meas = JuMP.owner_model(vref).measures[mindex]
            _check_meas_bounds(new_bounds, meas.data, _error = _error)
            if used_by_constraint(MeasureRef(JuMP.owner_model(vref), mindex))
                indices = JuMP.owner_model(vref).meas_to_constrs[mindex]
                meas_cindices = [meas_cindices; indices]
            end
        end
    end
    # check and update dependent constraints
    if used_by_constraint(vref) || length(meas_cindices) != 0
        for cindex in unique([meas_cindices; JuMP.owner_model(vref).var_to_constrs[JuMP.index(vref)]])
            constr = JuMP.owner_model(vref).constrs[cindex]
            new_constr = _update_constr_bounds(new_bounds, constr, _error = _error)
            JuMP.owner_model(vref).constrs[cindex] = new_constr
        end
    end
    _validate_bounds(JuMP.owner_model(vref), new_bounds, _error = _error)
    # add the bounds
    parameter_bounds(vref).intervals[pref] = IntervalSet(lower, upper)
    # update status
    JuMP.owner_model(vref).has_hold_bounds = true
    if is_used(vref)
        set_optimizer_model_ready(JuMP.owner_model(vref), false)
    end
    return
end

"""
    delete_parameter_bound(vref::HoldVariableRef, pref::ParameterRef)

Delete the parameter bound of the hold variable `vref` associated with the
infinite parameter `pref` if `vref` has such a bound. Note that any other
parameter bounds will be unaffected. Any constraints that employ `vref` will
be updated accordingly.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, x[1:2] in [0, 10])
2-element Array{ParameterRef,1}:
 x[1]
 x[2]

julia> @hold_variable(model, z, parameter_bounds = (x in [0, 1]))
z

julia> delete_parameter_bound(z, x[2])

julia> parameter_bounds(z)
Subdomain bounds (1): x[1] ∈ [0, 1]
```
"""
function delete_parameter_bound(vref::HoldVariableRef, pref::ParameterRef)
    # get the current bounds
    bounds = parameter_bounds(vref)
    # check if there are bounds for pref and act accordingly
    if haskey(bounds.intervals, pref)
        delete!(bounds.intervals, pref)
        # check for dependent measures that are used by constraints
        meas_cindices = []
        if used_by_measure(vref)
            for mindex in JuMP.owner_model(vref).var_to_meas[JuMP.index(vref)]
                if used_by_constraint(MeasureRef(JuMP.owner_model(vref), mindex))
                    indices = JuMP.owner_model(vref).meas_to_constrs[mindex]
                    meas_cindices = [meas_cindices; indices]
                end
            end
        end
        # check and update dependent constraints
        if used_by_constraint(vref) || length(meas_cindices) != 0
            for cindex in unique([meas_cindices; JuMP.owner_model(vref).var_to_constrs[JuMP.index(vref)]])
                constr = JuMP.owner_model(vref).constrs[cindex]
                new_constr = _rebuild_constr_bounds(constr, bounds)
                JuMP.owner_model(vref).constrs[cindex] = new_constr
            end
        end
    end
    return
end

"""
    delete_parameter_bounds(vref::HoldVariableRef)

Delete all the parameter bounds of the hold variable `vref`. Any constraints
that employ `vref` will be updated accordingly.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, x[1:2] in [0, 10])
2-element Array{ParameterRef,1}:
 x[1]
 x[2]

julia> @hold_variable(model, z, parameter_bounds = (x in [0, 1]))
z

julia> delete_parameter_bounds(z)

julia> parameter_bounds(z)
Subdomain bounds (0):
```
"""
function delete_parameter_bounds(vref::HoldVariableRef)
    # get the current bounds
    bounds = parameter_bounds(vref)
    # check if there are bounds and act accordingly
    if length(bounds) > 0
        _update_variable_param_bounds(vref, ParameterBounds())
        # check for dependent measures that are used by constraints
        meas_cindices = []
        if used_by_measure(vref)
            for mindex in JuMP.owner_model(vref).var_to_meas[JuMP.index(vref)]
                if used_by_constraint(MeasureRef(JuMP.owner_model(vref), mindex))
                    indices = JuMP.owner_model(vref).meas_to_constrs[mindex]
                    meas_cindices = [meas_cindices; indices]
                end
            end
        end
        # check and update dependent constraints
        if used_by_constraint(vref) || length(meas_cindices) != 0
            for cindex in unique([meas_cindices; JuMP.owner_model(vref).var_to_constrs[JuMP.index(vref)]])
                constr = JuMP.owner_model(vref).constrs[cindex]
                new_constr = _rebuild_constr_bounds(constr, bounds)
                JuMP.owner_model(vref).constrs[cindex] = new_constr
            end
        end
    end
    return
end

"""
    JuMP.set_name(vref::InfiniteVariableRef, root_name::String)

Extend [`JuMP.set_name`](@ref JuMP.set_name(::JuMP.VariableRef, ::String)) to set
names of infinite variables. Adds on to `root_name` the ending `(prefs...)`
where the parameter reference names are listed in the same format as input in
the parameter reference tuple.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @infinite_variable(m, vref(@infinite_parameter(m, t in [0, 1]))))
julia> name(vref)
"vref(t)"

julia> set_name(vref, "new_name")

julia> name(vref)
"new_name(t)"
```
"""
function JuMP.set_name(vref::InfiniteVariableRef, root_name::String)
    if length(root_name) == 0
        root_name = "noname"
    end
    # TODO do something about SparseAxisArrays (report array of values in order)
    # TODO list as vector and use ... like REPL if there are a lot
    prefs = parameter_refs(vref)
    param_names = _root_names(prefs)
    param_name_tuple = "("
    for i = 1:length(param_names)
        if i != length(param_names)
            param_name_tuple *= string(param_names[i], ", ")
        else
            param_name_tuple *= string(param_names[i])
        end
    end
    param_name_tuple *= ")"
    var_name = string(root_name, param_name_tuple)
    JuMP.owner_model(vref).var_to_name[JuMP.index(vref)] = var_name
    JuMP.owner_model(vref).name_to_var = nothing
    return
end

# Make a variable reference
function _make_variable_ref(model::InfiniteModel, index::Int)::GeneralVariableRef
    if isa(model.vars[index], InfiniteVariable)
        return InfiniteVariableRef(model, index)
    elseif isa(model.vars[index], PointVariable)
        return PointVariableRef(model, index)
    else
        return HoldVariableRef(model, index)
    end
end

"""
    JuMP.variable_by_name(model::InfiniteModel,
                          name::String)::Union{GeneralVariableRef, Nothing}

Extend [`JuMP.variable_by_name`](@ref JuMP.variable_by_name(::JuMP.Model, ::String))
for `InfiniteModel` objects. Return the variable reference assoociated with a
variable name. Errors if multiple variables have the same name. Returns nothing
if no such name exists.

**Examples**
```jldoctest; setup = :(using InfiniteOpt, JuMP; m = InfiniteModel(); @hold_variable(m, base_name = "var_name"))
julia> variable_by_name(m, "var_name")
var_name

julia> variable_by_name(m, "fake_name")

```
"""
function JuMP.variable_by_name(model::InfiniteModel,
                               name::String)::Union{GeneralVariableRef, Nothing}
    if model.name_to_var === nothing
        # Inspired from MOI/src/Utilities/model.jl
        model.name_to_var = Dict{String, Int}()
        for (var, var_name) in model.var_to_name
            if haskey(model.name_to_var, var_name)
                # -1 is a special value that means this string does not map to
                # a unique variable name.
                model.name_to_var[var_name] = -1
            else
                model.name_to_var[var_name] = var
            end
        end
    end
    index = get(model.name_to_var, name, nothing)
    if index isa Nothing
        return nothing
    elseif index == -1
        error("Multiple variables have the name $name.")
    else
        return _make_variable_ref(model, index)
    end
end

"""
    JuMP.all_variables(model::InfiniteModel)::Vector{GeneralVariableRef}

Extend [`JuMP.all_variables`](@ref JuMP.all_variables(::JuMP.Model)) to return a
list of all the variable references associated with `model`.

**Examples**
```julia-repl
julia> all_variables(model)
4-element Array{GeneralVariableRef,1}:
 y(t)
 w(t, x)
 y(0)
 z
```
"""
function JuMP.all_variables(model::InfiniteModel)::Vector{GeneralVariableRef}
    vrefs_list = Vector{GeneralVariableRef}(undef, JuMP.num_variables(model))
    indexes = sort([index for index in keys(model.vars)])
    counter = 1
    for index in indexes
        vrefs_list[counter] = _make_variable_ref(model, index)
        counter += 1
    end
    return vrefs_list
end


text:
## Here is the code from unit 1

# This treats everything as a dense matrix :(
A = [0 1 1 1 1 0 1 0 0 0
     1 0 0 0 0 0 0 0 0 0
     0 0 0 0 1 1 1 0 0 0
     1 1 0 0 1 1 0 1 0 0
     1 1 1 1 0 0 0 0 1 0
     0 0 1 0 0 0 1 0 1 1
     0 0 0 0 0 1 0 0 0 1
     0 0 0 0 0 0 0 0 1 0
     0 0 0 0 0 0 0 1 0 0
     0 0 0 0 0 0 0 0 0 0]
d = sum(A,dims=2)
##
using Printf
for i=1:10
  for j=1:10
    if i==j
      @printf("1 & ")
    elseif A[j,i] != 0
      @printf("-\\alpha/%i & ", d[j])
    else
      @printf("0 & ")
    end
  end
  println("\\\\")
end
## Build the I - alpha D^+ A'
alpha = 0.85
M = zeros(size(A)...)
for i=1:10
  for j=1:10
    if i==j
      M[i,j] = 1
    elseif A[j,i] != 0
      M[i,j] = -alpha*A[j,i]/d[j]
    else
      M[i,j] = 0
    end
  end
end
M[:,end] .= -alpha/10
M[end,end] += 1
M[:,1]
##
b = ones(size(A,1))/size(A,1)*(1-alpha)
x = M\b

##
using DelimitedFiles
using SparseArrays
ei,ej = copy.(eachcol(Int.(readdlm("wiki-simple.edges"))))
xy = readdlm("wiki-simple.xy")
pages = readlines("wiki-simple.nodes")
n = size(xy,1)
A = sparse(ei,ej,1,n,n)

## Simple PageRank
function simplepagerank(A,α,v)
  @assert(0 ≤ α < 1, "needs probably α")
  @assert(all(vi -> vi ≥ 0, v), "needs non-negative v")
  v = v ./ sum(v) # we can normalize for them.
  d = vec(sum(A,dims=2)) # compute the degrees
  x = copy(v) # start of with v
  nsteps = 2*ceil(Int,log(eps(1.0))/log(α)) # upper bound on steps
  for i=1:nsteps
    x = α*(A'*(x./d)) .+ (1-α).*v
  end
  return x/sum(x) # renormalize to probability
end
pr = simplepagerank(A,0.85,ones(n)./n)