text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

@testset "adstring" begin
    @test @inferred(adstring((30.4, -1.23), truncate=true)) ==
        @inferred(adstring([30.4, -1.23], truncate=true)) ==
        " 02 01 35.9  -01 13 48"
    @test @inferred(adstring(19.19321, truncate=true)) == "+19 11 35.5"
    @test @inferred(adstring(ten(36,24,15.015), -ten(8,24,36.0428), precision=3)) ==
        " 02 25 37.0010  -08 24 36.043"
    @test adstring.([30.4, -15.63], [-1.23, 48.41], precision=2) ==
        [" 02 01 36.000  -01 13 48.00", " 22 57 28.800  +48 24 36.00"]
    @test adstring.([(-58, 724)]) == [" 20 08 00.0  +724 00 00"]
end

@testset "airtovac" begin
    @test airtovac.([1234 6056.125]) ≈ [1234.0 6057.801930991426]
    @test @inferred(airtovac(2100)) ≈ 2100.666421596007
end

@testset "aitoff" begin
    @test @inferred(aitoff([227.23, 130], [-8.890, -35])) ==
        ([-137.92196683723276, 115.17541338020645],
         [-11.772527357473054, -44.491889962090085])
    @test @inferred(aitoff([375], [2.437])) ==
        ([16.63760711611838],[2.712427279646118])
    @test @inferred(aitoff((227.23, -8.890))) ==
        (-137.92196683723276,-11.772527357473054)
end

@testset "altaz2hadec" begin
    @test @inferred(altaz2hadec(59.086111, 133.30806, 43.07833)) ==
        @inferred(altaz2hadec((59.086111, 133.30806), 43.07833)) ==
        (336.68286017949157, 19.182449588316555)
    @test @inferred(altaz2hadec([15, 25, 35], [25.12, 45.32, -20.3],
                                [-23.44, 45.0, 52.5])) ==
        ([324.9881067314537, 256.7468302330436, 132.4919217875949] ,
         [44.38225395397647, 48.542947077386664, 67.33061196497327])
end

@testset "baryvel" begin
    dvelh_o, dvelb_o = @inferred(baryvel(2000))
    @test dvelh_o ≈ [1.582939967296732e-7, -1.0743272343303577e-7, -4.852410351888098e-8]
    @test dvelb_o ≈ [1.583299045307303e-7, -1.0736465601278539e-7, -4.849238001189245e-8]
    dvelh_o, dvelb_o = @inferred(baryvel(AstroLib.J2000, 1950))
    @test dvelh_o ≈ [-29.85888093436655, -4.684571288755146, -2.0305672776315777]
    @test dvelb_o ≈ [-29.849737241231153, -4.696440255370953, -2.035884790519881]
    dvelh_o, dvelb_o = @inferred(baryvel(jdcnv(1987, 4, 10, 0)))
    @test dvelh_o ≈ [9.514730366039178, -25.83862824722249, -11.202290688985778]
    @test dvelb_o ≈ [9.51331763088569,  -25.84985540772467, -11.207031996297271]
end

@testset "bprecess" begin
    ra, dec = @inferred(bprecess([ten(13, 42, 12.74)*15], [ten(8, 23, 17.69)],
                                 reshape(100*[-15*0.0257, -0.090], 2, 1)))
    @test ra  ≈ [204.93552515632123]
    @test dec ≈ [8.641287183886163]
    ra, dec = @inferred(bprecess(82, 19))
    @test ra  ≈ 81.26467916346334
    @test dec ≈ 18.959495700195394
    ra, dec = @inferred(bprecess([57], [23], 2024))
    @test ra  ≈ [56.26105898810067]
    @test dec ≈ [22.84693298145991]
    ra, dec = @inferred(bprecess([57], [23], reshape([9, 86], 2, 1), parallax=[1],
                                 radvel=[4]))
    @test ra  ≈ [56.25988479854577]
    @test dec ≈ [22.83493370392355]
    ra, dec = @inferred(bprecess((-57, -23), 2024))
    @test ra  ≈  302.2593299643789
    @test dec ≈ -23.150089972802036
    ra, dec = bprecess((-57, -23), [9, 86], parallax=1, radvel=4) # Inferred Type Error
    @test ra  ≈  302.2580376402947
    @test dec ≈ -23.16208183899836
end

@testset "calz_unred" begin
    @test calz_unred.(collect(900:1000:9900), ones(Float64, 10), -0.1) ≈
        [1.0, 0.43189326452379095, 0.5203675483533704, 0.594996469192435,
         0.6569506252451913, 0.7080829505773865, 0.7502392743978797, 0.7861262388745882,
         0.8151258710444882,0.8390325371659836]
end

# The values used for the testset are from running the code. They are slightly
# different from the output of the co_aberration routine of IDL AstroLib, as
# the function here uses an updated method to find mean obliquity
@testset "co_aberration" begin
    d_ra, d_dec = @inferred(co_aberration(jdcnv(1987, 4, 10, 0), ten(2,46,11.331)*15,
                                          ten(49,20,54.54), 1))
    @test d_ra ≈ -18.692441865574867
    @test d_dec ≈ -9.070782150537646
    ao, bo =  @inferred(co_aberration([57555.0, -6.44311e5], [302.282, 69.5667],
                                      [37.1519, 20.6847]))
    @test ao[1] ≈ 21.673056337579048
    @test ao[2] ≈ 18.496516329468466
    @test bo[1] ≈ -6.773070772568567
    @test bo[2] ≈ 2.9205843718089457
end

# The values used for the testset are from running the code. They are slightly
# different from the output of the co_aberration routine of IDL AstroLib, as
# the function here uses an updated method to find mean obliquity
@testset "co_nutate" begin
    an,bn,cn,dn,en = @inferred(co_nutate([jdcnv(2028,11,13,4,56), jdcnv(2013, 4, 16)],
                                         [10, 160],[80,30]))
    @test an ≈ [0.001209279097382776, 0.002465026741191423   ]
    @test bn ≈ [0.0017471459185713911, -0.0018134211486149354]
    @test cn ≈ [0.40904016038217567, 0.4090340058477726      ]
    @test dn ≈ [14.8593894278967, 12.102640377483143         ]
    @test en ≈ [2.7038090372351267, -5.86229256359996        ]
    ra_out, dec_out, eps_out, d_psi_out, d_eps_out = @inferred(co_nutate(2.451545e6,
                                                                         325, 0))
    @test ra_out ≈ -0.0035484441576727477
    @test dec_out ≈ -0.00034017946720967174
    @test eps_out ≈ 0.4090646078966446
    @test d_psi_out ≈ -13.923152677481191
    @test d_eps_out ≈ -5.773909654153591
end

@testset "co_refract" begin
    @test @inferred(co_refract(0.8)) ≈ 0.3714184384944585
    @test co_refract.([5.86,20], 50, 568.967, 273, 0.15, to_observe=true) ≈
        [5.94252628176525 , 20.02627520026167]
    @test @inferred(co_refract(14, 15000)) ≈ 13.990329255193124
end

@testset "ct2lst" begin
    @test ct2lst.(-76.72, -4, [DateTime(2008, 7, 30, 15, 53)]) ≈ [11.356505172312609]
    @test ct2lst.(9, [jdcnv(2015, 11, 24, 12, 21)]) ≈ [17.159574059885927]
end

# Test daycnv with Gregorian Calendar in force.
@testset "daycnv" begin
    @test @inferred(daycnv(2440000.0)) == DateTime(1968, 05, 23, 12)
    # Test daycnv with Julian Calendar in force (same result as IDL AstroLib's
    # daycnv).
    @test @inferred(daycnv(2000000.0)) == DateTime(763, 09, 18, 12)
    @test @inferred(daycnv(0.0)) == DateTime(-4713, 11, 24, 12)
end

@testset "deredd" begin
    by0, m0, c0, ub0 = @inferred(deredd([0.5, -0.5], [0.2, 0.5], [1, 1], [1, 1],
                                        [0.1, 0.3]))
    @test by0 ≈ [-0.3,0.5]
    @test m0  ≈ [1.165,1.0]
    @test c0  ≈ [0.905,1.0]
    @test ub0 ≈ [-0.665,0.3]
end

@testset "eci2geo" begin
    lat, long, alt = @inferred(eci2geo([0], [0], [0], [2452343]))
    @test lat  ≈ [0]
    @test long ≈ [12.992783145436988]
    @test alt  ≈ [-6378.137]
    lat, long, alt = @inferred(eci2geo((6978.137, 0, 0), jdcnv("2015-06-30T14:03:12.857")))
    @test lat  ≈ 0
    @test long ≈ 230.87301833205856
    @test alt  ≈ 600
    # Test `eci2geo' is the inverse of `geo2eci'
    jd = @inferred(get_juldate())
    lat, long, alt = @inferred(eci2geo(geo2eci(10, 10, 10, jd), jd))
    @test lat  ≈ 10
    @test long ≈ 10
    @test alt  ≈ 10
end

# The values used for the testset are from running the code. However they have been
# correlated with the output from eq2hor routine of IDL AstroLib, with
# differences only in the least significant digits.
@testset "eq2hor" begin
    alt_o, az_o, ha_o = eq2hor(259.20238631600944, 49.674907472176095,
                               AstroLib.J2000, "pbo", B1950=true)
    @test alt_o ≈ 43.804935603297004
    @test az_o ≈ 56.74141416977815
    @test ha_o ≈ 291.2750910119525
    alt_o, az_o, ha_o = eq2hor(142.2933457820434, -34.218006262991786, 2e6, 54.435, -34.78,
                               1000.34, ws=true, B1950=true, precession = false,
                               nutate=false, aberration=false,
                               refract=false, pressure = 500.345, temperature = 293.343)
    @test alt_o ≈ 1.3449999999999924
    @test az_o ≈ 359.43
    @test ha_o ≈ 359.3108663499664
    alt_o, az_o, ha_o = eq2hor(3.3222617779538037, 15.190516725395284, 2466879.7083333,
                               "kpno", pressure = 711, temperature = 273)
    @test alt_o ≈ 37.91138916818937
    @test az_o ≈ 264.918333213257
    @test ha_o ≈ 54.61193155973385
    alt_o, az_o, ha_o = @inferred(AstroLib._eq2hor(259.52076321839485, 49.62352289872951,
                                                   Float64(AstroLib.J2000), 43.0783,
                                                   -89.865, 0.0, NaN, NaN, false, false,
                                                   true, true, true, true))
    @test alt_o ≈ 43.687900264047116
    @test az_o ≈ 56.68399934960606
    @test ha_o ≈ 291.0817909922114
    alt_o, az_o = eq2hor(hor2eq(25, 55, 2.05e6, "pbo")[1:2]..., 2.05e6, "pbo")
    @test alt_o ≈ 24.99993224731665
    @test az_o ≈ 54.99993893556545
end

@testset "eqpole" begin
    x, y = @inferred(eqpole([100], [35], southpole=true))
    @test x ≈ [-111.18287262822456]
    @test y ≈ [ -19.604540237028665]
    x, y = @inferred(eqpole([80], [19]))
    @test x ≈ [72.78853915267848]
    @test y ≈ [12.83458333897169]
end

# The values used for the testset are from running the code. However they have been
# correlated with the output from euler routine of IDL AstroLib, with
# differences only in the least significant digits.
@testset "euler" begin
    glong, glat = @inferred(euler(299.590315, 35.201604, 1))
    @test glong ≈ 71.33498957116959
    @test glat ≈ 3.0668335310640984
    ra, dec = @inferred(euler((71.33498957116959, 3.0668335310640984), 2))
    @test ra ≈ 299.590315
    @test dec ≈ 35.201604
    elong, elat = @inferred(euler(3.141592653589793, 0.6143838917832061, 3,
                                  FK4 = true, radians=true))
    @test elong ≈ 2.8679433080257506
    @test elat ≈ 0.557258307291505
    ra, dec = @inferred(euler((2.8679433080257506, 0.557258307291505), 4,
                               FK4 = true, radians=true))
    @test ra ≈ 3.141592653589793
    @test dec ≈ 0.6143838917832061
    ecl, gal = euler.(30.45, 76.54, [5,6])
    @test ecl[1] ≈ 103.50477919192522
    @test ecl[2] ≈ 18.01965967759107
    @test gal[1] ≈ 194.96100731553986
    @test gal[2] ≈ 34.46136801388695
    @test @inferred(euler(183/pi, pi/180, 2, FK4=false, radians=true)) ==
        (5.682517110086799, 0.947078051715398)
    glong, glat = @inferred(euler([0.45, 130], [16.28, 53.65], 5))
    @test glong ≈ [96.9525940157568, 138.09922696730337]
    @test glat ≈ [-43.90672396295434, 46.95527026543361]
    @test_throws ErrorException @inferred(euler((45,45), 7))
end

@testset "flux2mag" begin
    @test flux2mag.([1.5e-12, 8.7e-15, 4.4e-10]) ≈
        [8.459771852360795, 14.051201868453454, 2.291368308784527]
    @test @inferred(flux2mag(1)) ≈ -21.1
    @test @inferred(flux2mag(5.2e-15)) ≈ 14.609991640913002
    @test @inferred(flux2mag(5.2e-15, 15)) ≈ 20.709991640913003
    @test flux2mag(5.2e-15, ABwave=15) ≈ 27.423535345634598 # Inferred Type Error
end

@testset "gal_uvw" begin
    u, v, w = @inferred(gal_uvw([ten(1,9,42.3)*15], [ten(61,32,49.5)], [627.89], [77.84],
                                [-321.4], [1e3/129], lsr=true))
    @test u ≈ [118.2110474553902]
    @test v ≈ [-466.4828898385057]
    @test w ≈ [88.16573278565097]
    u, v, w = @inferred(gal_uvw(1, 2, 3, 4, 5, 6))
    @test u ≈  4.0228405867158745
    @test v ≈  3.7912174342038227
    @test w ≈ -3.1700191400725464
end

@testset "geo2eci" begin
    x, y, z = @inferred(geo2eci([0], [0], [0], [2452343]))
    @test x ≈ [6214.846433007192]
    @test y ≈ [-1433.9858454345972]
    @test z ≈ [0.0]
    x, y, z = @inferred(geo2eci((0,0,0), jdcnv("2015-06-30T14:03:12.857")))
    @test x ≈ -4024.8671780315185
    @test y ≈ 4947.835465127513
    @test z ≈ 0.0
end

@testset "geo2geodetic" begin
    lat, long, alt = @inferred(geo2geodetic([90], [0], [0], "Jupiter"))
    @test lat  ≈   [90]
    @test long ≈    [0]
    @test alt  ≈ [4638]
    lat, long, alt = @inferred(geo2geodetic((90, 0, 0)))
    @test lat  ≈ 90
    @test long ≈  0
    @test alt  ≈ 21.38499999999931
    lat, long, alt = @inferred(geo2geodetic((43.16, -24.32, 3.87), 8724.32, 8619.19))
    @test lat  ≈  43.849399515234516
    @test long ≈ -24.32
    @test alt  ≈  53.53354478670836
    lat, long, alt = @inferred(geo2geodetic([43.16], [-24.32], [3.87], 8724.32, 8619.19))
    @test lat  ≈ [ 43.849399515234516]
    @test long ≈ [-24.32]
    @test alt  ≈ [ 53.53354478670836]
end

@testset "geo2mag" begin
    lat, long = @inferred(geo2mag(ten(35,0,42), ten(135,46,6), 2016))
    @test lat  ≈  36.86579228937769
    @test long ≈ -60.184060536651614
    lat, long = @inferred(geo2mag([15], [24], 2016))
    @test lat  ≈ [ 11.452100529696096]
    @test long ≈ [-169.86030510727102]
end

@testset "geodetic2geo" begin
    lat, long, alt = @inferred(geodetic2geo([90], [0], [0], "Jupiter"))
    @test lat  ≈ [90]
    @test long == [0]
    @test alt  ≈ [-4638]
    lat, long, alt = @inferred(geodetic2geo((90, 0, 0)))
    @test lat ≈ 90
    @test long == 0
    @test alt ≈ -21.38499999999931
    lat, long, alt = @inferred(geodetic2geo((43.16, -24.32, 3.87), 8724.32, 8619.19))
    @test lat ≈ 42.46772711708433
    @test long == -24.32
    @test alt ≈ -44.52902080669082
    lat, long, alt = @inferred(geodetic2geo([43.16], [-24.32], [3.87], 8724.32, 8619.19))
    @test lat ≈ [42.46772711708433]
    @test long == [-24.32]
    @test alt ≈ [-44.52902080669082]
    # Test geodetic2geo is the inverse of geo2geodetic, within a certain
    # tolerance.
    lat, long, alt = @inferred(geodetic2geo(geo2geodetic(67.2,13.4,1.2)))
    @test lat ≈ 67.2 atol = 1e-8
    @test long == 13.4
    @test alt ≈ 1.2 atol = 1e-9
end

# Test get_date with mixed keywords.
@testset "get_date" begin
    @test @inferred(get_date(DateTime(2001,09,25,14,56,14), old=true,timetag=true)) ==
        @inferred(get_date(2001,09,25,14,56,14, old=true,timetag=true)) ==
        @inferred(get_date("2001-09-25T14:56:14", old=true,timetag=true)) ==
        "25/09/2001T14:56:14"
    @test @inferred(get_date(DateTime(2001,09,25,14,56,14))) ==
        @inferred(get_date(2001,09,25,14,56,14)) ==
        @inferred(get_date("2001-09-25T14:56:14")) ==
        "2001-09-25"
    @test get_date.([DateTime(2024), Date(2016, 3, 14)]) ==
        get_date.([Date(2024), "2016-03-14"]) ==
        get_date.(["2024-01", DateTime(2016, 3, 14)]) == ["2024-01-01", "2016-03-14"]
end

@testset "gcirc" begin
    @test gcirc.(0, [0,1,2], [1,2,3], [2,3,4], [3,4,5]) ≈
        [1.222450611061632, 2.500353926443337, 1.5892569925227757]
    @test @inferred(gcirc(0,  120, -43,   175, +22))     ≈  1.590442261600714
    @test @inferred(gcirc(1, (120, -43),  175, +22))     ≈  415908.56615322345
    @test @inferred(gcirc(2,  120, -43,  (175, +22)))    ≈  296389.3666794745
    @test @inferred(gcirc(0, (120, -43), (175, +22)) )   ≈  1.590442261600714
    @test gcirc.(1, [120], [-43],  175, +22)             ≈ [415908.56615322345]
    @test gcirc.(2,  120, -43,  [175], [+22])            ≈ [296389.3666794745]
    @test_throws ErrorException @inferred(gcirc(3, 0, 0, 0, 0))
end

@testset "hadec2altaz" begin
    alt1, az1 = @inferred(hadec2altaz([0], [11.978165], [ten(43,4,42)]))
    @test alt1 ≈ [58.89983166666667]
    @test az1  ≈ [180.0]
    @test @inferred(hadec2altaz((0, 11.978165), ten(43,4,42), ws=true)[2]) ≈ 0.0
    alt1, az1 = 50, 20
    alt2, az2 = @inferred(hadec2altaz(altaz2hadec(alt1, az1, 40), 40))
    @test alt1 ≈ alt2
    @test az1  ≈ az2
end

@testset "helio_jd" begin
    @test helio_jd(juldate(2016, 6, 15, 11, 40), ten(20, 9, 7.8)*15, ten(37, 9, 7)) ≈
        57554.98808289718
    @test @inferred(helio_jd(1000, 23, 67, B1950=true)) ≈ 999.9997659545342
    @test @inferred(helio_jd(2000, 12, 88, diff=true))  ≈ -167.24845957792076
end

@testset "helio_rv" begin
    @test @inferred(helio_rv(helio_jd(juldate(94, 10, 25, 17, 30), ten(04, 38, 16)*15,
                                      ten(20, 41, 05)), 46487.5303, 2.0563056, -6, 59.3)) ≈
        -62.965570109145034
    @test helio_rv.([0.1, 0.9], 0, 1, 0, 100, 0.6, 45) ≈
        [-45.64994926111004, 89.7820347358485]
end

# The values used for the testset are from running the code. However they have been
# correlated with the output from helio routine of IDL AstroLib, with
# differences only in the least significant digits (except for `hrad`` output of Mars)
@testset "helio" begin
    @test_throws ErrorException @inferred(helio(jdcnv(2005,07,17,2,6,9), 10))
    hrad_out, hlong_out, hlat_out = @inferred(helio(jdcnv(2000,08,23,0), 2, true))
    @test hrad_out ≈ 0.7213758288364316
    @test hlong_out ≈ 3.462574978561256
    @test hlat_out ≈ 0.050393862449261535
    hrad_out, hlong_out, hlat_out = @inferred(helio([AstroLib.J2000], [7]))
    @test hrad_out[1] ≈ 19.921687573575788
    @test hlong_out[1] ≈ 316.4011812518626
    @test hlat_out[1] ≈ -0.6846115653974465
end

# The values used for the testset are from running the code. However they have been
# correlated with the output from hor2eq routine of IDL AstroLib, with
# differences only in the least significant digits.
@testset "hor2eq" begin
    ra_o, dec_o, ha_o = hor2eq(43.6879, 56.684, AstroLib.J2000, "pbo", B1950=true)
    @test ra_o ≈ 259.3943636338339
    @test dec_o ≈ 49.67396411401468
    @test ha_o ≈ 291.0833116221485
    ra_o, dec_o, ha_o = hor2eq(1.345, 359.43, 2e6, 54.435, -34.78, 1000.34, ws=true,
                               B1950=true, precession = false, nutate=false,
                               aberration=false, refract=false, pressure = 500.345,
                               temperature = 293.343)
    @test ra_o ≈ 142.2933457820434
    @test dec_o ≈ -34.218006262991786
    @test ha_o ≈ 359.3108663499664
    ra_o, dec_o, ha_o = hor2eq(ten(37,54,41), ten(264,55,06), 2466879.7083333,
                               "kpno", pressure = 711, temperature = 273)
    @test ra_o ≈ 3.3222617779538037
    @test dec_o ≈ 15.190516725395284
    @test ha_o ≈ 54.61193186104758
    ra_o, dec_o, ha_o = @inferred(AstroLib._hor2eq(43.6879, 56.684, Float64(AstroLib.J2000),
                                                   43.0783, -89.865, 0.0, NaN, NaN, false,
                                                   false, true, true, true, true))
    @test ra_o ≈ 259.52076321839485
    @test dec_o ≈ 49.62352289872951
    @test ha_o ≈ 291.0817908419628
    ra_o, dec_o = hor2eq(eq2hor(45, 60, 2e6, "pbo")[1:2]..., 2e6, "pbo")
    @test ra_o ≈ 45.0001735428487
    @test dec_o ≈ 59.999995143349956
end

# The values used for the testset are from running the code. However they have been
# correlated with the output from imf routine of IDL AstroLib, with
# differences only in the least significant digits.
@testset "imf" begin
    @test_throws ErrorException @inferred(imf([5], [-6.75], [0.9]))
    @test @inferred(imf([0.1, 0.01], [-0.6, -1], [ 0.007, 1.8, 110])) ≈
        [0.49627714725007616, 1.9757149090208912]
    @test imf.([[3],[5]], [[-1.35], [-0.6, -1.7]], [[0.1, 100], [0.007, 1.8, 110]]) ≈
        [[0.038948937298846235], [0.027349915327755464]]
end

# The values used for the testset are from running the code. However they have been
# correlated with the output from ismeuv routine of IDL AstroLib, with
# differences only in the least significant digits.
@testset "ismeuv" begin
    @test @inferred(ismeuv(304, 1e20)) ≈ 58.30508020244554
    @test ismeuv.([50, 75, 343], 1e18) ≈
        [0.004486567212077276, 0.01601219540569213, 0.7811526665834057]
    @test ismeuv.([96, 41, 233], 1e18, 1e17, 1e17) ≈
        [0.04657033979688484, 0.0035293162863089807, 0.30662192095758833]
    @test ismeuv.([96.6, 41.056, 233.19], 1e18, 1e17, 1e17, true) ≈
        [0.04733922036264192, 0.003543553203857853, 0.3103333860793942]
    @test ismeuv.([480, 910], 1e19, 5e17, 5e17) ≈
        [14.273937721143753, 62.68379602891728]
    @test @inferred(ismeuv(4500, 1e18)) ≈ 0
end

@testset "jdcnv" begin
    @test @inferred(jdcnv(-4713, 11, 24, 12)) ≈ 0.0
    @test @inferred(jdcnv(763, 09, 18, 12)) == @inferred(jdcnv("763-09-18T12")) == 2000000
    @test (jd=1234567.89; @inferred(jdcnv(daycnv(jd))) == jd)
    @test jdcnv.([DateTime(2016, 07, 31), "1969-07-20"]) ==
        jdcnv.([Date(2016, 07, 31), DateTime(1969, 07, 20)]) ==
        jdcnv.(["2016-07-31", Date(1969, 07, 20)])
end

@testset "jprecess" begin
    ra, dec = @inferred(jprecess([ten(13, 39, 44.526)*15], [ten(8, 38, 28.63)],
                                 reshape(100*[-15*0.0259, -0.093], 2, 1)))
    @test ra  ≈ [205.5530845731372]
    @test dec ≈ [8.388247441904628]
    ra, dec = @inferred(jprecess(82, 19))
    @test ra  ≈ 82.73568745151148
    @test dec ≈ 19.036972917272056
    ra, dec = @inferred(jprecess([57], [23], 2024))
    @test ra  ≈ [57.74049975335702]
    @test dec ≈ [23.150053754297726]
    ra, dec = @inferred(jprecess([57], [23], reshape([9, 86], 2, 1), parallax=[1],
                        radvel=[4]))
    @test ra  ≈ [57.74180294549785]
    @test dec ≈ [23.16200582079095]
    ra, dec = @inferred(jprecess((-57, -23), 2024))
    @test ra  ≈ 303.73910971499015
    @test dec ≈ -22.846895476784482
    ra, dec = jprecess((-57, -23), [9, 86], parallax=1, radvel=4) # Inferred Type Error
    @test ra  ≈ 303.7402950607101
    @test dec ≈ -22.834931625610313
end

# Test juldate with Gregorian Calendar in force.  This also makes sure precision
# of the result is high enough.  Note that "juldate(dt::DateTime) =
# Dates.datetime2julian(dt)-2.4e6" would not be precise.
@testset "juldate" begin
    @test @inferred(juldate(DateTime(2016, 1, 1, 8))) ≈ (57388.5 + 1//3)
    # Test juldate with Julian Calendar in force, for different centuries.  This
    # also makes sure precision of the result is high enough.
    @test @inferred(juldate(1582, 10, 1, 20))    ≈ (-100843 + 1//3)
    @test @inferred(juldate("1000-01-01T20"))    ≈ (-313692 + 1//3)
    @test @inferred(juldate("100-10-25T20"))     ≈ (-642119 + 1//3)
    @test @inferred(juldate(-4713, 1, 1, 12))    ≈ -2.4e6
    @test @inferred(juldate(2016, 06, 30, 00, 05, 53, 120)) ≈
        jdcnv(2016, 06, 30, 00, 05, 53, 120) - 2.4e6
    # Test daycnv and juldate together, with Gregorian Calendar in force.  Note that
    # they are not expected to be one the inverse of the other during Julian
    # Calendar.
    @test (dt=DateTime(2016, 1, 1, 20, 45, 33, 457);
           @inferred(daycnv(juldate(dt) + 2.4e6)) == dt)
end

@testset "kepler_solver" begin
    for e in 0:0.1:1
        M = collect(0:1e-3:2pi)
        E = mod2pi.(kepler_solver.(M, e))
        # Make sure E is the solution of the Kepler's Equation with high precision.
        @test M ≈ E .- e .* sin.(E) rtol = 1e-15
    end
    @test kepler_solver.([pi/4, pi/6, 8pi/3], 0) ≈ [pi/4, pi/6, 2pi/3]
    @test @inferred(kepler_solver(0, 1)) == 0.0
    @test_throws AssertionError @inferred(kepler_solver(pi, -0.5))
    @test_throws AssertionError @inferred(kepler_solver(pi,  1.5))
end

@testset "lsf_rotate" begin
    vel, lsf = @inferred(lsf_rotate(3, 90))
    @test @inferred(length(vel)) == @inferred(length(lsf)) == 61
    vel, lsf = @inferred(lsf_rotate(5, 10))
    @test vel ≈ collect(-10.0:5.0:10.0)
    @test lsf ≈ [0.0, 0.556914447710896, 0.6933098861837907, 0.556914447710896, 0.0]
end

@testset "mag2flux" begin
    @test @inferred(mag2flux(4.83, 21.12))               ≈ 4.1686938347033296e-11
    @test flux2mag(mag2flux(15, ABwave=12.), ABwave=12)  ≈ 15.0 # Inferred Type Error
    @test @inferred(mag2flux(8.3))                       ≈ 1.7378008287493692e-12
    @test @inferred(mag2flux(8.3, 12))                   ≈ 7.58577575029182e-9
    @test mag2flux(8.3, ABwave=12)                       ≈ 3.6244115683017193e-7 # Inferred Type Error
end

@testset "mag2geo" begin
    lat, long = @inferred(mag2geo(90, 0, 2016))
    @test lat  ≈ 86.395
    @test long ≈ -166.29000000000002
    lat, long = @inferred(mag2geo([15], [24], 2016))
    @test lat  ≈ [11.702066965890157]
    @test long ≈ [-142.6357492442842]
    # Test geo2mag is approximately the inverse of mag2geo
    lat, long = @inferred(geo2mag(mag2geo(12.34, 56.78, 2016)..., 2016))
    @test lat  ≈ 12.34
    @test long ≈ 56.78
end

@testset "mean_obliquity" begin
    @test @inferred(mean_obliquity(AstroLib.J2000)) ≈ 0.4090926006005829
    @test mean_obliquity.(jdcnv.([DateTime(1916, 09, 22, 03, 39),
                             DateTime(2063, 10, 13, 09)])) ≈
        [0.4092816887615259, 0.40894777540460037]
end

@testset "month_cnv" begin
    @test month_cnv.([" januavv  ", "SEPPES ", " aUgUsT", "la"]) == [1, 9, 8, -1]
    @test month_cnv.([2, 12, 6], short=true, low=true) == ["feb", "dec", "jun"]
    @test @inferred(month_cnv(5, up=true)) == "MAY"
    @test (list=[1, 2, 3]; month_cnv.(month_cnv.(list)) == list)
    @test (list=["July", "March", "November"]; month_cnv.(month_cnv.(list)) == list)
end

@testset "moonpos" begin
    ra, dec, dis, lng, lat = @inferred(moonpos(jdcnv(1992, 4, 12)))
    @test ra  ≈ 134.68846854844108
    @test dec ≈ 13.768366630560255
    @test dis ≈ 368409.68481612665
    @test lng ≈ 133.16726428105378
    @test lat ≈ -3.2291264192144356
    ra, dec, dis, lng, lat = @inferred(moonpos([2457521], radians=true))
    @test ra  ≈ [2.2587950290926178]
    @test dec ≈ [0.26183388011392217]
    @test dis ≈ [385634.68772395694]
    @test lng ≈ [2.232459255739553]
    @test lat ≈ [-0.059294466326164315]
end

@testset "mphase" begin
    @test mphase.([2457520, 2457530, 2457650]) ≈
        [0.2781695910737857, 0.9969808583803166, 0.9580708477591693]
end

@testset "nutate" begin
    long, obl = @inferred(nutate(jdcnv(1987, 4, 10)))
    @test long ≈ -3.787931077110755
    @test obl  ≈  9.442520698644401
    long, obl = @inferred(nutate(2457521))
    @test long ≈ -4.401443629818089
    @test obl  ≈ -9.26823431959121
    long, obl = @inferred(nutate([2457000, 2458000]))
    @test long ≈ [ 4.327189321653877, -9.686089990639474]
    @test obl  ≈ [-9.507794266102866, -6.970768250588256]
end

@testset "paczynski" begin
    @test @inferred(paczynski(-1e-10)) ≈  -1e10
    @test @inferred(paczynski(1e-1))   ≈  10.037461005722337
    @test @inferred(paczynski(-1))     ≈  -1.3416407864998738
    @test @inferred(paczynski(10))     ≈   1.0001922892047386
    @test @inferred(paczynski(-1e10))  ≈  -1
end

@testset "planck_freq" begin
    @test @inferred(planck_freq(2000, 5000)) ≈ 6.1447146126144004e-30
end

@testset "planck_wave" begin
    @test @inferred(planck_wave(2000, 5000)) ≈ 8.127064833530511e-24
end

# The values used for the testset are from running the code. However they have been
# correlated with the output from planet_coord routine of IDL AstroLib, with
# differences only in the least significant digits
@testset "planet_coords" begin
    @test_throws ErrorException @inferred(planet_coords(DateTime(2013, 07, 22,
                                                                 03, 19, 06),0))
    ra_out, dec_out = @inferred(planet_coords([AstroLib.J2000, 2.45e6], [2,8]))
    @test ra_out[1] ≈ 239.8965221579066
    @test ra_out[2] ≈ 294.55483837772476
    @test dec_out[1] ≈ -18.450868549676304
    @test dec_out[2] ≈ -20.992319312642262
    ra_out, dec_out = @inferred(planet_coords(2.45e6, 9))
    @test ra_out ≈ 238.3131013864547
    @test dec_out ≈ -6.964788781133789
    @test @inferred(planet_coords(juldate(1), 3) == (0, 0))
end

@testset "polrec" begin
    x, y = @inferred(polrec([1, 2, 3], [pi, pi/2.0, pi/4.0]))
    @test x ≈ [-1.0, 0.0, 1.5*sqrt(2.0)]
    @test y ≈ [ 0.0, 2.0, 1.5*sqrt(2.0)]
    x, y = @inferred(polrec((2, 135), degrees=true))
    @test x ≈ -sqrt(2)
    @test y ≈  sqrt(2)
end

@testset "posang" begin
    @test @inferred(posang(1, ten(13, 25, 13.5), ten(54, 59, 17),
                           ten(13, 23, 55.5), ten(54, 55, 31))) ≈ -108.46011246802047
    @test posang.(0, [0,1,2], [1,2,3], [2,3,4], [3,4,5]) ≈
        [1.27896824717634, 1.6840484573313608, 0.2609280020139511]
    @test @inferred(posang(0,  120, -43,   175, +22))     ≈ -1.5842896165356724
    @test @inferred(posang(1, (120, -43),  175, +22))     ≈ 82.97831348792039
    @test @inferred(posang(2,  120, -43,  (175, +22)))    ≈ 50.02816530382374
    @test @inferred(posang(0, (120, -43), (175, +22)))    ≈ -1.5842896165356724
    @test posang.(1, [120], [-43],  175, +22)  ≈ [82.97831348792039]
    @test posang.(2,  120, -43,  [175], [+22]) ≈ [50.02816530382374]
    @test_throws ErrorException @inferred(posang(3, 0, 0, 0, 0))
end

@testset "precess" begin
    ra1, dec1 = @inferred(precess((ten(2,31,46.3)*15, ten(89,15,50.6)), 2000, 1985))
    @test ra1  ≈ 34.09470328718033
    @test dec1 ≈ 89.19647174928589
    ra2, dec2 = @inferred(precess([ten(21, 59, 33.053)*15], [ten(-56, 59, 33.053)],
                                  1950, 1975, FK4=true))
    @test ra2  ≈ [330.3144305418865]
    @test dec2 ≈ [-56.87186126487889]
end

# The values used for the testset are from running the code. However they have been
# correlated with the output from precess_cd routine of IDL AstroLib, with
# differences only in the least significant digits.
@testset "precess_cd" begin
    @test @inferred(precess_cd([30 60; 60 90], 1950, 2000, [13, 8], [43, 23])) ≈
        [30.919029003435927 62.343060521017435;
         61.93905850970097 93.56509103294071  ]
    @test @inferred(precess_cd([30 60; 60 90], 2000, 1950, [13, 8], [43, 23])) ≈
        [30.919029003435927 62.343060521017435;
         61.93905850970097 93.56509103294071  ]
    @test @inferred(precess_cd([12.45 56.7; 66 89], 2000, 1985, [67.4589455, 0.345345],
                               [37.94291666666666, 89.26405555555556])) ≈
        [963.4252080520984 4387.890452343343  ;
         5107.504395433958 6887.55936949333   ]
    @test @inferred(precess_cd([30.0 28.967; 60.45 90.65], 2000, 1975, [13, 10.658],
                               [35.54, 67], true)) ≈
        [64.78429186351575 62.637156996728194 ;
         130.49379143419722 195.9699513801844 ]
end

@testset "precess_xyz" begin
    x1, y1, z1 = @inferred(precess_xyz((1.2, 2.3, 1.7), 2000, 2050))
    @test x1 ≈ 1.165933061423247
    @test y1 ≈ 2.313228746401996
    @test z1 ≈ 1.7057470102860104
    x2, y2, z2 = @inferred(precess_xyz([0.7, -2.4], [3.3, 6.6], [0, 4], 2000, 2016))
    @test x2 ≈ [0.688187142071843,   -2.429815562246262]
    @test y2 ≈ [3.3024835038223532,   6.591359330834213]
    @test z2 ≈ [0.001079105285993004, 3.9962455511755794]
end

@testset "premat" begin
    @test @inferred(premat(1967, 1982, FK4=true)) ≈
        [0.9999933170034135    -0.0033529069683496567 -0.0014573823699636742;
         0.00335290696825777    0.9999943789886484    -2.443304965138481e-6 ;
         0.0014573823701750721 -2.4431788671274868e-6  0.9999989380147651   ]
    @test @inferred(premat(1995, 2003)) ≈
        [ 0.9999980977132219    -0.0017889257711428855 -0.0007773766929507687;
          0.0017889257711354528  0.9999983998707707    -6.953448226403318e-7 ;
          0.0007773766929678732 -6.953257000046125e-7   0.9999996978424512   ]
end

@testset "radec" begin
    @test @inferred(radec(15.90, -0.85)) == (1.0, 3.0, 36.0, -0.0, 51.0, 0.0)
    @test @inferred(radec(-0.85,15.9)) == (23.0,56.0,36.0,15.0,54.0,0.0)
    @test @inferred(radec(-20,4,hours=true)) == (4.0,0.0,0.0,4.0,0.0,0.0)
    @test @inferred(radec([15.90, -0.85], [-0.85,15.9])) ==
        ([1.0, 23.0], [3.0, 56.0], [36.0, 36.0],
         [-0.0, 15.0], [51.0, 54.0], [0.0, 0.0])
end

@testset "recpol" begin
    r = a = zeros(Float64, 3)
    r, a = @inferred(recpol([0, sqrt(2.0), 2.0*sqrt(3.0)], [0, sqrt(2.0), 2.0]))
    @test r ≈ [0.0,  2.0,  4.0]
    @test a ≈ [0.0, pi/4.0, pi/6.0]
    r, a = @inferred(recpol(1, 1))
    @test r ≈ sqrt(2.0)
    @test a ≈ pi/4.0
    # Test polrec is the inverse of recpol
    xi, yi, x, y, = 6.3, -2.7, 0.0, 0.0
    x, y = @inferred(polrec(recpol((xi, yi), degrees=true), degrees=true))
    @test x ≈ xi
    @test y ≈ yi
end

@testset "rhotheta" begin
    ρ, θ = @inferred(rhotheta(41.623, 1934.008, 0.2763, 0.907, 59.025, 23.717, 219.907,
                              1980))
    @test ρ ≈ 0.41101776646245836
    @test θ ≈ 318.4242564860495
end

# Test also it's the inverse of ten.
@testset "sixty" begin
    @test @inferred(sixty(-51.36))                    ≈ [-51.0, 21.0, 36.0]
    @test @inferred(ten(sixty(-0.10934835545824395))) ≈ -0.10934835545824395
    @test @inferred(sixty(1))                         ≈ [1.0, 0.0, 0.0]
end

@testset "sphdist" begin
    @test sphdist.([0,1,2], [1,2,3], [2,3,4], [3,4,5]) ≈
        [1.222450611061632, 2.500353926443337, 1.5892569925227762]
    @test @inferred(sphdist(120, -43, 175, +22))      ≈  1.5904422616007134
    @test sphdist.([120], [-43], 175, +22) ≈ [1.5904422616007134]
    @test sphdist.(120, -43, [175], [+22]) ≈ [1.5904422616007134]
end

@testset "sunpos" begin
    ra, dec, lon, obl = @inferred(sunpos(jdcnv(1982, 5, 1)))
    @test ra  ≈ 37.88589057369026
    @test dec ≈ 14.909699471099517
    @test lon ≈ 40.31067053890748
    @test obl ≈ 23.440840980112657
    ra, dec, lon, obl = @inferred(sunpos(jdcnv.([DateTime(2016, 5, 10)]), radians=true))
    @test ra  ≈ [0.8259691339090751]
    @test dec ≈ [0.3085047454107549]
    @test lon ≈ [0.8687853454154388]
    @test obl ≈ [0.40901175207670365]
    ra, dec, lon, obl = @inferred(sunpos([2457531]))
    @test ra  ≈ [59.71655864208797]
    @test dec ≈ [20.52127006818727]
    @test lon ≈ [61.824436793991545]
    @test obl ≈ [23.434648653514724]
end

# Make sure string and numerical inputs are consistent (IDL implementation of "ten" is not).
@testset "ten" begin
    @test @inferred(ten(0, -23, 34)) == @inferred(ten((0, -23, 34))) ==
        @inferred(ten([0, -23, 34])) == ten(" : 0 -23 :: 34") == -0.37388888888888894
    @test @inferred(ten(-0.0, 60)) == @inferred(ten((-0.0, 60))) ==
        @inferred(ten([-0.0, 60])) == @inferred(ten("-0.0 60")) == -1.0
    @test @inferred(ten(-5, -60, -3600)) == @inferred(ten((-5, -60, -3600))) ==
        @inferred(ten([-5, -60, -3600])) == @inferred(ten("  -5: :-60: -3600")) == -3.0
    @test ten("") == 0.0
    @test ten.([0, -0.0, -5], [-23, 60, -60], [34, 0, -3600]) ==
        ten.([(0, -23,34), ":-0.0:60", (-5, -60, -3600)]) ==
        ten.(["0   -23 :: 34", (-0.0, 60), " -5:-60: -3600"]) ==
        [-0.37388888888888894, -1.0, -3.0]
    @test ten.([12.0, -0.0], [24, 30]) == ten.([" 12::24", " -0:30: "]) == [12.4, -0.5]
end

@testset "tic_one" begin
    min2, tic1 = @inferred(tic_one(30.2345, 12.74, 10))
    @test min2 ≈ 30.333333333333332
    @test tic1 ≈ 7.554820000000081
    min2, tic1 = @inferred(tic_one(45, 50, 4, true))
    @test min2 ≈ 46.0
    @test tic1 ≈ 50.0
    min2, tic1 = @inferred(tic_one(pi\8, tics(90, 45, 1000, 10)...))
    @test min2 ≈ 2.5
    @test tic1 ≈ 1.0318357862412286
end

@testset "ticpos" begin
    @test ticpos.([16,8,4],[1024,512,256], [150,75,37.5]) ==
                  [(256.0, 4, "Degrees"), (128.0, 2, "Degrees"), (64.0, 1, "Degrees")]
    @test @inferred(ticpos(2, 512, 75)) == (128.0, 30, "Arc Minutes")
    @test @inferred(ticpos(1.5, 512, 75)) == (85.33333333333333, 15, "Arc Minutes")
    @test @inferred(ticpos(1.5, 512, 50)) == (56.888888888888886, 10, "Arc Minutes")
    @test @inferred(ticpos(1.6, 1024, 50)) == (53.333333333333336, 5, "Arc Minutes")
    @test @inferred(ticpos(0.2, 512, 75)) == (85.33333333333333, 2, "Arc Minutes")
    @test @inferred(ticpos(0.5, 512, 10)) == (17.066666666666666, 1, "Arc Minutes")
    @test @inferred(ticpos(0.1, 1024, 50)) == (85.33333333333333, 30, "Arc Seconds")
    @test @inferred(ticpos(0.08, 1024, 40)) == (53.333333333333336, 15, "Arc Seconds")
    @test @inferred(ticpos(0.025, 512, 50)) == (56.888888888888886, 10, "Arc Seconds")
    @test @inferred(ticpos(pi/100, 1024, 40)) == (45.27073936836133, 5, "Arc Seconds")
    @test @inferred(ticpos(0.06, 2048, 20)) == (18.962962962962965, 2, "Arc Seconds")
    @test @inferred(ticpos(0.016, 1024, 20)) == (17.77777777777778, 1, "Arc Seconds")
end

@testset "tics" begin
    @test @inferred(tics(30, 90, 30, 1)) == (3.8666666666666667, 480)
    @test @inferred(tics(30, 90, 3, 3, true)) == (4.0, 240)
    @test @inferred(tics(30, 70, 3, 1, true)) == (0.75, 60)
    @test tics.([30,50],[70,60], [6,12], [3,0.5], [true, false]) ==
                [(3.75,120), (0.55,30)]
    @test @inferred(tics(45, 55, 30, 0.5)) == (0.725, 15)
    @test @inferred(tics(45, 60, 10, 0.1)) == (0.1, 10)
    @test @inferred(tics(55, 60, 100.0, 1/2)) == (0.66, 2)
    @test @inferred(tics(25, 30, 50, 2, true)) == (2.45, 1)
    @test @inferred(tics(20, 80, 600, 0.03)) == (0.04159722222222222, 0.25)
    @test @inferred(tics(25, 75, 500, 0.02)) == (0.02772222222222222, 0.16666666666666666)
    @test @inferred(tics(10, 12, 25, 0.01)) == (0.016666666666666666, 0.08333333333333333)
    @test @inferred(tics(20, 80, 6000, 0.03)) ==
        (0.055546296296296295, 0.03333333333333333)
    @test @inferred(tics(30, 60, 200, 0.02, true)) ==
        (0.02763888888888889, 0.016666666666666666)
    @test @inferred(tics(60, 70, 125, 0.001)) ==
        (0.0017222222222222222, 0.008333333333333333)
    @test @inferred(tics(10, 12, 25, 0.01, true)) == (0.01, 0.0033333333333333335)
    @test @inferred(tics(130, 180, 1000, 0.0004)) == (0.000555, 0.0016666666666666668)
    @test @inferred(tics(60, 70, 5500//2, 0.003)) ==
        (0.003818055555555556, 0.0008333333333333334)
    @test @inferred(tics(30, 150, 4000, 0.002, true)) ==
        (0.0027770833333333337, 0.0003333333333333333)
    @test @inferred(tics(9.5, 14.5, 5000, 0.002)) ==
        (0.002777222222222222, 0.00016666666666666666)
    @test @inferred(tics(90, 45, 1000, 10)) == (11.1, -30)
    ticsize, incr = @inferred(tics(30, 70, 50, 0.1))
    @test ticsize ≈ 0.10208333333333333
    @test incr ≈ 5
    ticsize, incr = @inferred(tics(pi/3, pi/2, 60.0, 7.5, true))
    @test ticsize ≈ 14.085212463632736
    @test incr ≈ 0.5
end

@testset "true_obliquity" begin
    @test @inferred(true_obliquity(AstroLib.J2000)) ≈ 0.4090646078966446
    @test true_obliquity.(jdcnv.([DateTime(2016, 08, 23, 03, 39, 06),
                             DateTime(763, 09, 18, 12)])) ≈
        [0.4090133706884892, 0.41188965892279295]
end

@testset "kepler_solver" begin
    @test @inferred(trueanom(8pi/3, 0.7))              ≈ 2.6657104039293764
    @test trueanom.([pi/4, pi/6, 8pi/3], 0)            ≈ [pi/4, pi/6, 2pi/3]
    @test @inferred(trueanom(3pi/2, 0.8))              ≈ -2.498091544796509
    @test @inferred(trueanom(0.1, 1))                  ≈ pi
    @test_throws AssertionError @inferred(trueanom(pi, -0.5))
    @test_throws DomainError @inferred(trueanom(pi,  1.5))
end

# The values used for the testset are from running the code. However they have been
# correlated with the output from uvbybeta routine of IDL AstroLib, with
# differences only in the least significant digits.
@testset "uvbybeta" begin
    @test_throws ErrorException @inferred uvbybeta(NaN, NaN, NaN, 9, NaN)
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(-0.009, 0.11, 0.68, 1)
    @test te_o ≈ 12719.770257555097
    @test mv_o ≈ -0.08913388245565224
    @test eby_o ≈ 0.038295081967213124
    @test delm_o ≈ -0.007787863429965028
    @test radius_o ≈ 2.549253655861009
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.415, -0.069, 0.308, 2)
    @test te_o ≈ 12280.215047866226
    @test mv_o ≈ -7.112880191374548
    @test eby_o ≈ 0.4040160642570281
    @test isnan(delm_o)
    @test radius_o ≈ 90.29892932880738
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.273, -0.051, 0.051, 3)
    @test te_o ≈ 18733.204098765906
    @test mv_o ≈ -5.453550292618578
    @test eby_o ≈ 0.33158225857187795
    @test isnan(delm_o)
    @test radius_o ≈ 27.825306168356892
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.29343, 0.09, 0.001, 4)
    @test te_o ≈ 25713.1878000067
    @test mv_o ≈ -4.816603735163776
    @test eby_o ≈ 0.6430112570356474
    @test isnan(delm_o)
    @test radius_o ≈ 0.16611917129253512
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.03, 0.21345, 0.3976, 5)
    @test te_o ≈ 10363.815662133818
    @test mv_o ≈ 4.411076600000009
    @test eby_o ≈ -0.038394585955939545
    @test delm_o ≈ -0.03113954500000002
    @test radius_o ≈ 0.4863555386334208
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.18547, -0.04356, 0.2576, 6)
    @test te_o ≈ 7308.315628104092
    @test mv_o ≈ 7.562860949999999
    @test eby_o ≈ -0.00736941600000007
    @test delm_o ≈ 0.2449243827527609
    @test radius_o ≈ 0.17062904371865295
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.18547, -0.04356, 0.2576, 6,
                                                             2.72001)
    @test te_o ≈ 7308.315628104092
    @test mv_o ≈ 6.1817999999999955
    @test eby_o ≈ -0.07476300000000027
    @test delm_o ≈ 0.21601800458420034
    @test radius_o ≈ 0.3223035928541916
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.18547, -0.04356, 0.2576, 6,
                                                             2.82)
    @test te_o ≈ 7334.8182052629045
    @test mv_o ≈ 7.767123299999994
    @test eby_o ≈ 0.003229999999999733
    @test delm_o ≈ 0.2454149000000012
    @test radius_o ≈ 0.154376398171302
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.18547, -0.04356, 0.2576, 6,
                                                             2.8201)
    @test te_o ≈ 7343.722680158638
    @test mv_o ≈ 7.682471615999958
    @test eby_o ≈ 0.0043099999999998695
    @test delm_o ≈ 0.24506736642000115
    @test radius_o ≈ 0.16018990575428574
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.1923, 0.2186, 0.2783, 6)
    @test te_o ≈ 7252.900460401587
    @test mv_o ≈ 7.219177499999998
    @test eby_o ≈ -0.011614630064603726
    @test delm_o ≈ -0.019005720258414266
    @test radius_o ≈ 0.2024586440816956
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.216, 0.167, 0.785, 7)
    @test te_o ≈ 7108.688707930596
    @test mv_o ≈ 1.538642957897997
    @test eby_o ≈ 0.005426551873147473
    @test delm_o ≈ 0.010495724945951107
    @test radius_o ≈ 2.8675336941552185
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.206, 0.162, 0.786, 7,
                                                             2.601)
    @test te_o ≈ 7144.241056060125
    @test mv_o ≈ -0.8920873665838025
    @test eby_o ≈ -0.12713573819611615
    @test delm_o ≈ 0.05272866901000031
    @test radius_o ≈ 8.706345749177931
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.394, 0.184, 0.382, 8)
    @test te_o ≈ 5755.671513413262
    @test mv_o ≈ 3.7737408311021916
    @test eby_o ≈ -0.028753800373693867
    @test delm_o ≈ 0.034414237187698926
    @test radius_o ≈ 1.6597544196694192
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.6501, 0.184, 0.382, 8)
    @test te_o ≈ 5452.352129194577
    @test mv_o ≈ 5.192225306337874
    @test eby_o ≈ 0.19958051948051947
    @test delm_o ≈ 0.032731030075137835
    @test radius_o ≈ 0.9924525316819005
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.7901, 0.184, 0.0145, 8)
    @test te_o ≈ 3977.185981725903
    @test mv_o ≈ 9.893069966897801
    @test eby_o ≈ -0.022297953964194334
    @test delm_o ≈ 0.5456850156962505
    @test radius_o ≈ 0.26255047396110764
    te_o, mv_o, eby_o, delm_o, radius_o = @inferred uvbybeta(0.7901, 0.037, 0.039, 8)
    @test te_o ≈ 3985.833020767027
    @test mv_o ≈ 10.356215131980953
    @test eby_o ≈ 0.0027659846547315762
    @test delm_o ≈ 0.6944746189113468
    @test radius_o ≈ 0.2106831009861503
end

# Test airtovac is its inverse (it isn't true only around 2000, just avoid those values)
@testset "vactoair" begin
    @test @inferred(vactoair(2000)) ≈ 1999.3526230448367
    @test airtovac.(vactoair.(collect(1000:300:4000))) ≈ collect(1000:300:4000)
end

@testset "xyz" begin
    x, y, z, vx, vy, vz = @inferred(xyz([51200.5 + 64 / 86400], 2000))
    @test x  ≈ [0.5145687092402946]
    @test y  ≈ [-0.7696326261820777]
    @test z  ≈ [-0.33376880143026394]
    @test vx ≈ [0.014947267514081075]
    @test vy ≈ [0.008314838205475709]
    @test vz ≈ [0.003606857607574784]
end

@testset "ydn2md" begin
    @test ydn2md.(2016, [60, 234]) == [Date(2016, 02, 29), Date(2016, 08, 21)]
    @test @inferred(ymd2dn(ydn2md(2016, 60))) == 60
end

@testset "ymd2dn" begin
    @test ymd2dn.([Date(2015,3,5), Date(2016,3,5)]) == [64, 65]
    @test @inferred(ydn2md(2016, ymd2dn(Date(2016, 09, 16)))) == Date(2016, 09, 16)
end

@testset "zenpos" begin
    ra, dec = @inferred(zenpos(2.457514070138889e6, 45, 45))
    @test ra ≈ 1.9915758420649625
    @test dec ≈ 0.7853981633974483
    ra, dec = @inferred(zenpos(DateTime(2015, 11, 24, 13, 21), 43.16, -24.32, 4))
    @test ra ≈ 3.1232737646297757
    @test dec ≈ 0.7532841051607526
    ra, dec = @inferred(zenpos(jdcnv(2017, 01, 30, 04, 30), ten(35,0,42), ten(135,46,6)))
    @test ra ≈ 5.809762417608341
    @test dec ≈ 0.6110688599440813
end


text:
#Matrices related to stochastic processes

import Base: iterate
export AiryProcess, BrownianProcess, WhiteNoiseProcess, next!

abstract type StochasticProcess{T<:Real} end

struct WhiteNoiseProcess{T<:Real} <: StochasticProcess{T}
    dt::T
end

struct BrownianProcess{T<:Real} <: StochasticProcess{T}
    dt::T
end

struct AiryProcess{S<:Real, T<:Real} <: StochasticProcess{T}
    dt::T
    beta::S
end

done(p::StochasticProcess{T}, x...) where {T} = false #Processes can always go on forever

#######################
# White noise process #
#######################

iterate(p::WhiteNoiseProcess, state=()) = (randn()*sqrt(p.dt), ())

####################
# Brownian process #
####################

function iterate(p::BrownianProcess{T}, state=zero(T)) where {T}
    newx = state + randn()*sqrt(p.dt)
    return newx, newx
end

################
# Airy process #
################

"""
Like next, but update only the state of the AiryProcess

Skip the eigenvalue computation, which gets expensive
"""
function next!(p::AiryProcess{T}, S::SymTridiagonal{T}=SymTridiagonal(T[-(2/p.dt^2)], T[])) where {T} # TODO maybe change this to iterate!
    t = (size(S, 1)-1)*p.dt

    #Discredited Airy operator plus diagonal noise
    x = inv(p.dt^2)
    push!(S.dv, -2x - t + 2/sqrt(p.dt*p.beta)*randn())
    push!(S.ev, x)

    return S
end
function iterate(p::AiryProcess{T}, S::SymTridiagonal{T}= SymTridiagonal(T[-(2/p.dt^2)], T[])) where {T}
    t = (size(S, 1)-1)*p.dt

    #Discredited Airy operator plus diagonal noise
    x = inv(p.dt^2)
    push!(S.dv, -2x - t + 2/sqrt(p.dt*p.beta)*randn())
    push!(S.ev, x)

    return (eigmax(S), S)
end

ext: jl


text:
module Renderer

export respond, json, redirect_to, html, flax, include_asset, has_requested, css_asset, js_asset
export respond_with_json, respond_with_html

using Genie, Util, JSON, Genie.Configuration, HttpServer, App, Router, Logger, Macros

if IS_IN_APP
  eval(:(using $(App.config.html_template_engine), $(App.config.json_template_engine)))
  eval(:(const HTMLTemplateEngine = $(App.config.html_template_engine)))
  eval(:(const JSONTemplateEngine = $(App.config.json_template_engine)))

  export HTMLTemplateEngine, JSONTemplateEngine

  const DEFAULT_LAYOUT_FILE = App.config.renderer_default_layout_file
else
  const DEFAULT_LAYOUT_FILE = :app
end

const CONTENT_TYPES = Dict{Symbol,String}(
  :html   => "text/html",
  :plain  => "text/plain",
  :json   => "application/json",
  :js     => "text/javascript",
  :xml    => "text/xml",
)

const VIEWS_FOLDER = "views"
const LAYOUTS_FOLDER = "layouts"

"""
    html(resource::Symbol, action::Symbol, layout::Symbol = DEFAULT_LAYOUT_FILE, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Dict{Symbol,String}

Invokes the HTML renderer of the underlying configured templating library.
"""
function html(resource::Union{Symbol,String}, action::Union{Symbol,String}, layout::Union{Symbol,String} = DEFAULT_LAYOUT_FILE, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Dict{Symbol,String}
  HTMLTemplateEngine.html(resource, action, layout; parse_vars(vars)...)
end

"""
    respond_with_html(resource::Symbol, action::Symbol, layout::Symbol = DEFAULT_LAYOUT_FILE, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Response

Invokes the HTML renderer of the underlying configured templating library and wraps it into a `HttpServer.Response`.
"""
function respond_with_html(resource::Union{Symbol,String}, action::Union{Symbol,String}, layout::Union{Symbol,String} = DEFAULT_LAYOUT_FILE, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Response
  html(resource, action, layout, check_nulls; vars...) |> respond
end

function flax(resource::Union{Symbol,String}, action::Union{Symbol,String}, layout::Union{Symbol,String} = DEFAULT_LAYOUT_FILE, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Dict{Symbol,String}
  HTMLTemplateEngine.flax(resource, action, layout; parse_vars(vars)...)
end

"""
    json(resource::Symbol, action::Symbol, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Dict{Symbol,String}

Invokes the JSON renderer of the underlying configured templating library.
"""
function json(resource::Union{Symbol,String}, action::Union{Symbol,String}, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Dict{Symbol,String}
  JSONTemplateEngine.json(resource, action; parse_vars(vars)...)
end

"""
    respond_with_json(resource::Symbol, action::Symbol, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Response

Invokes the JSON renderer of the underlying configured templating library and wraps it into a `HttpServer.Response`.
"""
function respond_with_json(resource::Union{Symbol,String}, action::Union{Symbol,String}, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Response
  json(resource, action, check_nulls; vars...) |> respond
end

"""
    redirect_to(location::String, code::Int = 302, headers = Dict{AbstractString,AbstractString}()) :: Response

Sets redirect headers and prepares the `Response`.
"""
function redirect_to(location::String, code = 302, headers = Dict{AbstractString,AbstractString}()) :: Response
  headers["Location"] = location
  respond(Dict{Symbol,AbstractString}(:plain => "Redirecting you to $location"), code, headers)
end
function redirect_to(named_route::Symbol, code = 302, headers = Dict{AbstractString,AbstractString}()) :: Response
  redirect_to(link_to(named_route), code, headers)
end

"""
    has_requested(content_type::Symbol) :: Bool

Checks wheter or not the requested content type matches `content_type`.
"""
function has_requested(content_type::Symbol) :: Bool
  task_local_storage(:__params)[:response_type] == content_type
end

"""
    respond{T}(body::Dict{Symbol,T}, code::Int = 200, headers = Dict{AbstractString,AbstractString}()) :: Response

Constructs a `Response` corresponding to the content-type of the request.
"""
function respond(body::Dict{Symbol,T}, code::Int = 200, headers = Dict{AbstractString,AbstractString}())::Response where {T}
  sbody::String =   if haskey(body, :json)
                      headers["Content-Type"] = CONTENT_TYPES[:json]
                      body[:json]
                    elseif haskey(body, :html)
                      headers["Content-Type"] = CONTENT_TYPES[:html]
                      body[:html]
                    elseif haskey(body, :js)
                      headers["Content-Type"] = CONTENT_TYPES[:js]
                      body[:js]
                    elseif haskey(body, :plain)
                      headers["Content-Type"] = CONTENT_TYPES[:plain]
                      body[:plain]
                    else
                      Logger.log("Unsupported Content-Type", :err)
                      Logger.log(body)
                      Logger.@location

                      error("Unsupported Content-Type")
                    end

  Response(code, headers, sbody)
end
function respond(response::Tuple, headers = Dict{AbstractString,AbstractString}()) :: Response
  respond(response[1], response[2], headers)
end
function respond(response::Response) :: Response
  response
end
function respond{T}(body::String, params::Dict{Symbol,T}) :: Response
  r = params[:RESPONSE]
  r.data = body

  r |> respond
end
function respond(body::String) :: Response
  respond(Response(body))
end

"""
    http_error(status_code; id = "resource_not_found", code = "404-0001", title = "Not found", detail = "The requested resource was not found")

Constructs an error `Response`.
"""
function http_error(status_code; id = "resource_not_found", code = "404-0001", title = "Not found", detail = "The requested resource was not found")
  respond(detail, status_code, Dict{AbstractString,AbstractString}())
end

"""
    include_asset(asset_type::Symbol, file_name::String; fingerprinted = App.config.assets_fingerprinted) :: String

Returns the path to an asset. `asset_type` can be one of `:js`, `:css`. `file_name` should not include the extension.
`fingerprinted` is a `Bool` indicated wheter or not fingerprinted (unique hash) should be added to the asset's filename (used in production to invalidate caches).
"""
function include_asset(asset_type::Symbol, file_name::String; fingerprinted::Bool = App.config.assets_fingerprinted) :: String
  suffix = fingerprinted ? "-" * App.ASSET_FINGERPRINT * ".$(asset_type)" : ".$(asset_type)"
  "/$asset_type/$(file_name)$(suffix)"
end
function include_asset(asset_type::Symbol, file_name::Symbol; fingerprinted::Bool = App.config.assets_fingerprinted) :: String
  include_asset(asset_type, string(file_name), fingerprinted = fingerprinted)
end

"""
    css_asset(file_name::String; fingerprinted::Bool = App.config.assets_fingerprinted) :: String

Path to a css asset. `file_name` should not include the extension.
`fingerprinted` is a `Bool` indicated wheter or not fingerprinted (unique hash) should be added to the asset's filename (used in production to invalidate caches).
"""
function css_asset(file_name::String; fingerprinted::Bool = App.config.assets_fingerprinted) :: String
  include_asset(:css, file_name, fingerprinted = fingerprinted)
end

"""
    js_asset(file_name::String; fingerprinted::Bool = App.config.assets_fingerprinted) :: String

Path to a js asset. `file_name` should not include the extension.
`fingerprinted` is a `Bool` indicated wheter or not fingerprinted (unique hash) should be added to the asset's filename (used in production to invalidate caches).
"""
function js_asset(file_name::String; fingerprinted::Bool = App.config.assets_fingerprinted) :: String
  include_asset(:js, file_name, fingerprinted = fingerprinted)
end

function parse_vars(vars)
  pos_counter = 1
  for pair in vars
    if pair[1] != :check_nulls
      pos_counter += 1
      continue
    end

    for p in pair[2]
      if ! isa(p[2], Nullable)
        push!(vars, p[1] => p[2])
        continue
      end

      if isnull(p[2])
        return error_404()
      else
        push!(vars, p[1] => Base.get(p[2]))
      end
    end
  end

  vars
end

end


text:
include("emptygridworld.jl")
include("fourrooms.jl")
include("gotodoor.jl")
include("doorkey.jl")
include("collectgems.jl")
include("dynamicobstacles.jl")


text:
try
    using Revise
catch e
    @warn "Error initializing Revise" exception=(e, catch_backtrace())
end

atreplinit() do repl
    try
        @eval using OhMyREPL
    catch e
        @warn "Error initializing OhMyREPL" e
    end
end

max_stars_repo_name: garborg/mydots

max_issues_repo_name: garborg/mydots

max_forks_repo_name: garborg/mydots


text:
function Base.getindex(r::Union{AbstractStepRangeLen,AbstractLinRange}, i::Integer)
    Base.@_inline_meta
    @boundscheck checkbounds(r, i)
    unsafe_getindex(r, i)
end

###
### AbstractStepRangeLen
###
function _getindex_hiprec(r::AbstractStepRangeLen, i::Integer)  # without rounding by T
    u = i - _offset(r)
    return _ref(r) + u * step(r)
end

function _getindex_hiprec(
    r::AbstractStepRangeLen{<:Any,<:TwicePrecision,<:TwicePrecision},
    i::Integer
   )
    u = i - _offset(r)
    shift_hi, shift_lo = u * gethi(step_hp(r)), u * getlo(step_hp(r))
    x_hi, x_lo = add12(refhi(r), shift_hi)
    x_hi, x_lo = add12(x_hi, x_lo + (shift_lo + reflo(r)))
    return TwicePrecision(x_hi, x_lo)
end

for RT in (:StepMRangeLen,:StepSRangeLen)
    @eval begin
        function Base.getindex(
            r::$(RT){T,TwicePrecision{T},TwicePrecision{T}},
            s::OrdinalRange{<:Integer}
           ) where {T}
            @boundscheck checkbounds(r, s)
            soffset = 1 + round(Int, (_offset(r) - first(s))/step(s))
            soffset = clamp(soffset, 1, length(s))
            ioffset = first(s) + (soffset-1)*step(s)
            if step(s) == 1 || length(s) < 2
                newstep = step_hp(r)
            else
                newstep = Base.twiceprecision(step_hp(r)*step(s), Base.nbitslen(T, length(s), soffset))
            end
            if ioffset == _offset(r)
                return similar_type(r)(_ref(r), newstep, length(s), max(1,soffset))
            else
                return similar_type(r)(_ref(r) + (ioffset-_offset(r))*step_hp(r), newstep, length(s), max(1,soffset))
            end
        end
    end
end

# although these should technically not need to be completely typed for
# each, dispatch ignores TwicePrecision on the static version and only
# uses the first otherwise
function Base.unsafe_getindex(
    r::StepSRangeLen{T,TwicePrecision{T},TwicePrecision{T}},
    i::Integer
   ) where {T}
    # Very similar to _getindex_hiprec, but optimized to avoid a 2nd call to add12
    Base.@_inline_meta
    u = i - _offset(r)
    shift_hi, shift_lo = u * stephi(r), u * steplo(r)
    x_hi, x_lo = add12(refhi(r), shift_hi)
    return T(x_hi + (x_lo + (shift_lo + reflo(r))))
end

function Base.unsafe_getindex(
    r::StepMRangeLen{T,TwicePrecision{T},TwicePrecision{T}},
    i::Integer
   ) where {T}
    # Very similar to _getindex_hiprec, but optimized to avoid a 2nd call to add12
    Base.@_inline_meta
    u = i - _offset(r)
    shift_hi, shift_lo = u * stephi(r), u * steplo(r)
    x_hi, x_lo = add12(refhi(r), shift_hi)
    return T(x_hi + (x_lo + (shift_lo + reflo(r))))
end

function Base.unsafe_getindex(r::StepSRangeLen{T,R,S}, i::Integer) where {T,R,S}
    return T(_ref(r) + (i - _offset(r)) * step_hp(r))
end
function Base.unsafe_getindex(r::StepMRangeLen{T,R,S}, i::Integer) where {T,R,S}
    return T(_ref(r) + (i - _offset(r)) * step_hp(r))
end

###
### AbstractLinRange
###
function Base.unsafe_getindex(r::AbstractLinRange, i::Integer)
    return Base.lerpi(i-1, r.lendiv, r.start, r.stop)
end

function Base.getindex(r::AbstractLinRange, s::OrdinalRange{<:Integer})
    Base.@_inline_meta
    @boundscheck checkbounds(r, s)
    vfirst = unsafe_getindex(r, first(s))
    vlast  = unsafe_getindex(r, last(s))
    return LinMRange(vfirst, vlast, length(s))
end

function Base.getindex(r::LinSRange, s::Union{OneToSRange{T},UnitSRange{T},StepSRange{T}}) where {T<:Integer}
    Base.@_inline_meta
    @boundscheck checkbounds(r, s)
    vfirst = unsafe_getindex(r, first(s))
    vlast  = unsafe_getindex(r, last(s))
    return LinMRange(vfirst, vlast, length(s))
end

###
### StaticUnitRange
###
_in_unit_range(v::StaticUnitRange, val, i::Integer) = i > 0 && val <= last(v) && val >= first(v)

function Base.getindex(v::StaticUnitRange{T}, i::Integer) where T
    Base.@_inline_meta
    val = convert(T, first(v) + (i - 1))
    @boundscheck _in_unit_range(v, val, i) || throw(BoundsError(v, i))
    val
end

function Base.getindex(v::StaticUnitRange{T}, i::Integer) where {T<:Base.OverflowSafe}
    Base.@_inline_meta
    val = v.start + (i - 1)
    @boundscheck _in_unit_range(v, val, i) || throw(BoundsError(v, i))
    val % T
end


text:
export FaradayStep!,AmpereStep!,_update_dual_constitutive!,_update_primal_constitutive!

include("FaradayStep!.jl")
include("AmpereStep!.jl")
include("_update_dual_constitutive!.jl")
include("_update_primal_constitutive!.jl")


text:
#
#	text_manipulate.jl
#
#					Jul/05/2018
# --------------------------------------------------------------------
function text_to_dict_proc(fname_in)
	delimiter = "\t"
	dict_aa = text_to_dict_proc_exec(delimiter,fname_in)
#
	return	dict_aa
end
# --------------------------------------------------------------------
function text_to_dict_proc_exec(delimiter,fname_in)
	string = open(readstring, fname_in, "r")
	lines = split(string,"\n")
#
	dict_aa = Dict()
	llx = length(lines)
	println(llx)
	for it=1:llx
		line = lines[it]
		if (1 < length(line))
			cols = split(line,delimiter)
			key = cols[1]
			K = ["name","population","date_mod"]
			V = [cols[2],cols[3],cols[4]]
#			unit_aa = Dict(K,V)
			unit_aa = Dict("name" => cols[2],"population" => cols[3], "date_mod" => cols[4])
			dict_aa[key] = unit_aa
		end
	end
#
	return	dict_aa
end
# --------------------------------------------------------------------
function dict_display_proc(dict_aa)
#
	for key in sort(collect(keys(dict_aa)))
		print(key, '\t', dict_aa[key]["name"])
		print('\t', dict_aa[key]["population"])
		println('\t', dict_aa[key]["date_mod"])
	end

end
# --------------------------------------------------------------------
function dict_append_proc(dict_aa,key,name,population,date_mod)
	K = ["name","population","date_mod"]
	V = [name,population,date_mod]
#	unit_aa = Dict(K,V)
	unit_aa = Dict("name" => name,"population" => population, "date_mod" => date_mod)
	dict_aa[key] = unit_aa
#
	return	dict_aa
end
# --------------------------------------------------------------------
function dict_to_text_proc(dict_aa,fname_out)
	delimiter = "\t"
	dict_to_text_proc_exec(delimiter,dict_aa,fname_out)
#
end
#
# --------------------------------------------------------------------
function dict_to_text_proc_exec(delimiter,dict_aa,fname_out)
	fp_out = open(fname_out,"w")
#
	for key in sort(collect(keys(dict_aa)))
		str_out = string(key,delimiter,dict_aa[key]["name"])
		str_out = string(str_out,delimiter,dict_aa[key]["population"])
		str_out = string(str_out,delimiter,dict_aa[key]["date_mod"])
		str_out = string(str_out,'\n')
		write(fp_out,str_out)
	end
#
	close(fp_out)
end
#
# --------------------------------------------------------------------
function dict_update_proc(dict_aa,key,population_in)
#
	today = Libc.strftime("%F", time())
#
	unit_aa = Dict()
	unit_aa["name"] = dict_aa[key]["name"]
	unit_aa["population"] = population_in
	unit_aa["date_mod"] = today
	dict_aa[key] = unit_aa
#
	return	dict_aa
end
#
# --------------------------------------------------------------------
function csv_to_dict_proc(fname_in)
	delimiter = ","
	dict_aa = text_to_dict_proc_exec(delimiter,fname_in)
#
	return	dict_aa
end
# --------------------------------------------------------------------
function dict_to_csv_proc(dict_aa,fname_out)
	delimiter = ","
	dict_to_text_proc_exec(delimiter,dict_aa,fname_out)
#
end
#
# --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using GeoMakie, Makie

lons = LinRange(-179.5, 179.5, 360)
lats = LinRange(-89.5, 89.5, 180)

field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

source = LonLat()
dest = WinkelTripel()

xs, ys = xygrid(lons, lats)
Proj4.transform!(source, dest, vec(xs), vec(ys))

scene = surface(xs, ys; color = field, shading = false, show_axis = false, scale_plot = false)

geoaxis!(scene, -180, 180, -90, 90; crs = (src = source, dest = dest,))

coastlines!(scene, 1; crs = (src = source, dest = dest,))


text:
#############
# Promotion #
#############

# Promotion between two tensors promote the eltype and promotes
# symmetric tensors to tensors

@inline function Base.promote_rule(::Type{SymmetricTensor{order, dim, A, M}},
                                   ::Type{SymmetricTensor{order, dim, B, M}}) where {dim, A, B, order, M}
    SymmetricTensor{order, dim, promote_type(A, B), M}
end

@inline function Base.promote_rule(::Type{Tensor{order, dim, A, M}},
                                   ::Type{Tensor{order, dim, B, M}}) where {dim, A, B, order, M}
    Tensor{order, dim, promote_type(A, B), M}
end

@inline function Base.promote_rule(::Type{SymmetricTensor{order, dim, A, M1}},
                                   ::Type{Tensor{order, dim, B, M2}}) where {dim, A, B, order, M1, M2}
    Tensor{order, dim, promote_type(A, B), M2}
end

@inline function Base.promote_rule(::Type{Tensor{order, dim, A, M1}},
                                   ::Type{SymmetricTensor{order, dim, B, M2}}) where {dim, A, B, order, M1, M2}
    Tensor{order, dim, promote_type(A, B), M1}
end

# inlined promote (promote in Base is not inlined)
@inline function Base.promote(S1::T, S2::S) where {T <: AbstractTensor, S <: AbstractTensor}
    return convert(promote_type(T, S), S1), convert(promote_type(T, S), S2)
end
@inline Base.promote(S1::AbstractTensor{order, dim, T}) where {order, dim, T} = convert(Tensor{order, dim, T}, S1)

# base promotion that only promotes SymmetricTensor to Tensor but leaves eltype
@inline function promote_base(S1::Tensor{order, dim}, S2::SymmetricTensor{order, dim}) where {order, dim}
    return S1, convert(Tensor{order, dim}, S2)
end
@inline function promote_base(S1::SymmetricTensor{order, dim}, S2::Tensor{order, dim}) where {order, dim}
    return convert(Tensor{order, dim}, S1), S2
end

###############
# Conversions #
###############

# Identity conversions
@inline Base.convert(::Type{Tensor{order, dim, T}}, t::Tensor{order, dim, T}) where {order, dim, T} = t
@inline Base.convert(::Type{Tensor{order, dim, T, M}}, t::Tensor{order, dim, T, M}) where {order, dim, T, M} = t
@inline Base.convert(::Type{SymmetricTensor{order, dim, T}}, t::SymmetricTensor{order, dim, T}) where {order, dim, T} = t
@inline Base.convert(::Type{SymmetricTensor{order, dim, T, M}}, t::SymmetricTensor{order, dim, T, M}) where {order, dim, T, M} = t

# Change element type
@inline function Base.convert(::Type{Tensor{order, dim, T1}}, t::Tensor{order, dim, T2}) where {order, dim, T1, T2}
    apply_all(Tensor{order, dim}, @inline function(i) @inbounds T1(t.data[i]); end)
end

@inline function Base.convert(::Type{SymmetricTensor{order, dim, T1}}, t::SymmetricTensor{order, dim, T2}) where {order, dim, T1, T2}
    apply_all(SymmetricTensor{order, dim}, @inline function(i) @inbounds T1(t.data[i]); end)
end

# Peel off the M but define these so that convert(typeof(...), ...) works
@inline Base.convert(::Type{Tensor{order, dim, T1, M}}, t::Tensor{order, dim})                   where {order, dim, T1, M} = convert(Tensor{order, dim, T1}, t)
@inline Base.convert(::Type{SymmetricTensor{order, dim, T1, M}}, t::SymmetricTensor{order, dim}) where {order, dim, T1, M} = convert(SymmetricTensor{order, dim, T1}, t)
@inline Base.convert(::Type{Tensor{order, dim, T1, M}}, t::SymmetricTensor{order, dim})          where {order, dim, T1, M} = convert(Tensor{order, dim, T1}, t)
@inline Base.convert(::Type{SymmetricTensor{order, dim, T1, M}}, t::Tensor{order, dim})          where {order, dim, T1, M} = convert(SymmetricTensor{order, dim, T1}, t)

@inline Base.convert(::Type{Tensor{order, dim}}, t::SymmetricTensor{order, dim, T}) where {order, dim, T} = convert(Tensor{order, dim, T}, t)
@inline Base.convert(::Type{SymmetricTensor{order, dim}}, t::Tensor{order, dim, T}) where {order, dim, T} = convert(SymmetricTensor{order, dim, T}, t)
@inline Base.convert(::Type{Tensor}, t::SymmetricTensor{order, dim, T})             where {order, dim, T} = convert(Tensor{order, dim, T}, t)
@inline Base.convert(::Type{SymmetricTensor}, t::Tensor{order, dim, T})             where {order, dim, T} = convert(SymmetricTensor{order, dim, T}, t)

# SymmetricTensor -> Tensor
@inline function Base.convert(::Type{Tensor{2, dim, T1}}, t::SymmetricTensor{2, dim, T2}) where {dim, T1, T2}
    Tensor{2, dim}(@inline function(i,j) @inbounds T1(t[i,j]); end)
end

@inline function Base.convert(::Type{Tensor{4, dim, T1}}, t::SymmetricTensor{4, dim, T2}) where {dim, T1, T2}
    Tensor{4, dim}(@inline function(i,j,k,l) @inbounds T1(t[i,j,k,l]); end)
end

# Tensor -> SymmetricTensor
@inline function Base.convert(::Type{SymmetricTensor{order, dim, T1}}, t::Tensor{order, dim}) where {dim, order, T1}
    if issymmetric(t)
        return convert(SymmetricTensor{order, dim, T1}, symmetric(t))
    else
        throw(InexactError(:convert, SymmetricTensor{order, dim, T1}, t))
    end
end


text:
println("\n\n\nRunning tests ...")

using RandomFerns
using Base.Test

using FunctionalData

D = 2
N = 100
C = 10
S = 10
data = @p map (1:C) (x->randn(D,N).+5x.+randn(D,1)) | unstack | flatten
labels = @p map (1:C) (x->x*ones(1,N)) | unstack | flatten

a = RandFerns(data,labels);
# testdata = @p meshgrid linspace(-15,15,200) linspace(-15,15,200)
r = predict(a, data)
@test sum(vec(r).!=vec(labels)) < 100

a = RandFerns(data,labels; regression = true);
r = predict(a, data)
@test mean(abs(vec(r)-vec(labels))) < 0.5

println("   done running tests!")


text:
x = 10
println("x --> ", x)

x = x + 155
println("x --> ", x)

x = "Hello World!"
println("x --> ", x)

δ = 0.00001
println("δ --> ", δ)

print(π)


text:
@reset
@testset "RHF" begin
    @testset "Conventional" begin
        wf1 = @energy rhf;
        @test isapprox(wf1.energy, -74.9650028737304410, rtol=tol)
    end

    @testset "DF" begin
        @set {
              scf_alg df
              basis cc-pvdz
              jkfit cc-pvdz-jkfit
             }
        wf2 = @energy rhf;
        @test isapprox(wf2.energy, -76.02317515357174, rtol=tol)
    end
end


text:
module SymbolicRegression

# Types
export Population,
    PopMember,
    HallOfFame,
    Options,
    Node,

    #Functions:
    EquationSearch,
    SRCycle,
    calculateParetoFrontier,
    countNodes,
    copyNode,
    printTree,
    stringTree,
    evalTreeArray,
    differentiableEvalTreeArray,
    node_to_symbolic,
    symbolic_to_node,
    custom_simplify,
    simplifyWithSymbolicUtils,
    combineOperators,
    genRandomTree,

    #Operators:
    plus,
    sub,
    mult,
    square,
    cube,
    pow,
    div,
    log_abs,
    log2_abs,
    log10_abs,
    log1p_abs,
    acosh_abs,
    sqrt_abs,
    neg,
    greater,
    relu,
    logical_or,
    logical_and,

    # special operators
    gamma,
    erf,
    erfc,
    atanh_clip

using Distributed
using Printf: @printf
using Pkg
using Random: seed!
using FromFile
using Reexport
@reexport using LossFunctions

@from "Core.jl" import CONST_TYPE, maxdegree, Dataset, Node, copyNode, Options, plus, sub, mult, square, cube, pow, div, log_abs, log2_abs, log10_abs, log1p_abs, sqrt_abs, acosh_abs, neg, greater, greater, relu, logical_or, logical_and, gamma, erf, erfc, atanh_clip
@from "Utils.jl" import debug, debug_inline, is_anonymous_function
@from "EquationUtils.jl" import countNodes, printTree, stringTree
@from "EvaluateEquation.jl" import evalTreeArray, differentiableEvalTreeArray
@from "CheckConstraints.jl" import check_constraints
@from "MutationFunctions.jl" import genRandomTree
@from "LossFunctions.jl" import EvalLoss, Loss, scoreFunc
@from "PopMember.jl" import PopMember, copyPopMember
@from "Population.jl" import Population, bestSubPop
@from "HallOfFame.jl" import HallOfFame, calculateParetoFrontier, string_dominating_pareto_curve
@from "SingleIteration.jl" import SRCycle, OptimizeAndSimplifyPopulation
@from "InterfaceSymbolicUtils.jl" import node_to_symbolic, symbolic_to_node
@from "CustomSymbolicUtilsSimplification.jl" import custom_simplify
@from "SimplifyEquation.jl" import simplifyWithSymbolicUtils, combineOperators, simplifyTree
@from "ProgressBars.jl" import ProgressBar, set_multiline_postfix

include("Configure.jl")
include("Deprecates.jl")

"""
    EquationSearch(X, y[; kws...])

Perform a distributed equation search for functions which
describe the mapping f(X[:, j]) ≈ y[j]. Options are
configured using SymbolicRegression.Options(...),
which should be passed as a keyword argument to options.
One can turn off parallelism with `numprocs=0`,
which is useful for debugging and profiling.

# Arguments
- `X::AbstractMatrix{T}`:  The input dataset to predict `y` from.
    The first dimension is features, the second dimension is rows.
- `y::AbstractVector{T}`: The values to predict. Only a single feature
    is allowed, so `y` is a 1D array.
- `niterations::Int=10`: The number of iterations to perform the search.
    More iterations will improve the results.
- `weights::Union{AbstractVector{T}, Nothing}=nothing`: Optionally
    weight the loss for each `y` by this value (same shape as `y`).
- `varMap::Union{Array{String, 1}, Nothing}=nothing`: The names
    of each feature in `X`, which will be used during printing of equations.
- `options::Options=Options()`: The options for the search, such as
    which operators to use, evolution hyperparameters, etc.
- `numprocs::Union{Int, Nothing}=nothing`:  The number of processes to use,
    if you want `EquationSearch` to set this up automatically. By default
    this will be `4`, but can be any number (you should pick a number <=
    the number of cores available).
- `procs::Union{Array{Int, 1}, Nothing}=nothing`: If you have set up
    a distributed run manually with `procs = addprocs()` and `@everywhere`,
    pass the `procs` to this keyword argument.
- `runtests::Bool=true`: Whether to run (quick) tests before starting the
    search, to see if there will be any problems during the equation search
    related to the host environment.

# Returns
- `hallOfFame::HallOfFame`: The best equations seen during the search.
    hallOfFame.members gives an array of `PopMember` objects, which
    have their tree (equation) stored in `.tree`. Their score (loss)
    is given in `.score`. The array of `PopMember` objects
    is enumerated by size from `1` to `options.maxsize`.
"""
function EquationSearch(X::AbstractMatrix{T}, y::AbstractVector{T};
        niterations::Int=10,
        weights::Union{AbstractVector{T}, Nothing}=nothing,
        varMap::Union{Array{String, 1}, Nothing}=nothing,
        options::Options=Options(),
        numprocs::Union{Int, Nothing}=nothing,
        procs::Union{Array{Int, 1}, Nothing}=nothing,
        runtests::Bool=true
       ) where {T<:Real}

    dataset = Dataset(X, y,
                     weights=weights,
                     varMap=varMap)
    serial = (procs == nothing && numprocs == 0)
    parallel = !serial

    if runtests
        testOptionConfiguration(T, options)
        testDatasetConfiguration(dataset, options)
    end

    if dataset.weighted
        avgy = sum(dataset.y .* dataset.weights)/sum(dataset.weights)
        baselineMSE = Loss(dataset.y, ones(T, dataset.n) .* avgy, dataset.weights, options)
    else
        avgy = sum(dataset.y)/dataset.n
        baselineMSE = Loss(dataset.y, ones(T, dataset.n) .* avgy, options)
    end

    if options.seed !== nothing
        seed!(options.seed)
    end
    # Start a population on every process
    allPopsType = parallel ? Future : Tuple{Population,HallOfFame}
    allPops = allPopsType[]
    # Set up a channel to send finished populations back to head node
    channels = [RemoteChannel(1) for j=1:options.npopulations]
    bestSubPops = [Population(dataset, baselineMSE, npop=1, options=options, nfeatures=dataset.nfeatures) for j=1:options.npopulations]
    hallOfFame = HallOfFame(options)
    actualMaxsize = options.maxsize + maxdegree
    frequencyComplexity = ones(T, actualMaxsize)
    curmaxsize = 3
    if options.warmupMaxsizeBy == 0f0
        curmaxsize = options.maxsize
    end

    we_created_procs = false
    ##########################################################################
    ### Distributed code:
    ##########################################################################
    if parallel
        if numprocs == nothing && procs == nothing
            numprocs = 4
            procs = addprocs(4)
            we_created_procs = true
        elseif numprocs == nothing
            numprocs = length(procs)
        elseif procs == nothing
            procs = addprocs(numprocs)
            we_created_procs = true
        end
        if we_created_procs
            project_path = splitdir(Pkg.project().path)[1]
            activate_env_on_workers(procs, project_path, options)
            import_module_on_workers(procs, @__FILE__, options)
        end
        move_functions_to_workers(procs, options, dataset)
        if runtests
            test_module_on_workers(procs, options)
        end

        if runtests
            test_entire_pipeline(procs, dataset, options)
        end
    end
    cur_proc_idx = 1
    # Get the next worker process to give a job:
    function next_worker()::Int
        if parallel
            idx = ((cur_proc_idx-1) % numprocs) + 1
            cur_proc_idx += 1
            return procs[idx]
        else
            return 0
        end
    end

    for i=1:options.npopulations
        worker_idx = next_worker()
        new_pop = if parallel
            (@spawnat worker_idx Population(dataset, baselineMSE,
                                            npop=options.npop, nlength=3,
                                            options=options,
                                            nfeatures=dataset.nfeatures),
            HallOfFame(options))
        else
            (Population(dataset, baselineMSE,
                        npop=options.npop, nlength=3,
                        options=options,
                        nfeatures=dataset.nfeatures), HallOfFame(options))
        end
        push!(allPops, new_pop)
    end
    # 2. Start the cycle on every process:
    for i=1:options.npopulations
        worker_idx = next_worker()
        allPops[i] = if parallel
            @spawnat worker_idx let
                tmp_pop, tmp_best_seen = SRCycle(dataset, baselineMSE, fetch(allPops[i])[1], options.ncyclesperiteration, curmaxsize, copy(frequencyComplexity)/sum(frequencyComplexity), verbosity=options.verbosity, options=options)
                tmp_pop = OptimizeAndSimplifyPopulation(dataset, baselineMSE, tmp_pop, options, curmaxsize)
                if options.batching
                    for i_member=1:(options.maxsize + maxdegree)
                        tmp_best_seen.members[i_member].score = scoreFunc(dataset, baselineMSE, tmp_best_seen.members[i_member].tree, options)
                    end
                end
                (tmp_pop, tmp_best_seen)
            end
        else
            tmp_pop, tmp_best_seen = SRCycle(dataset, baselineMSE, allPops[i][1], options.ncyclesperiteration, curmaxsize, copy(frequencyComplexity)/sum(frequencyComplexity), verbosity=options.verbosity, options=options)
            tmp_pop = OptimizeAndSimplifyPopulation(dataset, baselineMSE, tmp_pop, options, curmaxsize)
            if options.batching
                for i_member=1:(options.maxsize + maxdegree)
                    tmp_best_seen.members[i_member].score = scoreFunc(dataset, baselineMSE, tmp_best_seen.members[i_member].tree, options)
                end
            end
            (tmp_pop, tmp_best_seen)
        end
    end

    debug(options.verbosity > 0 || options.progress, "Started!")
    total_cycles = options.npopulations * niterations
    cycles_remaining = total_cycles
    if options.progress
        progress_bar = ProgressBar(1:cycles_remaining;
                                   width=options.terminal_width)
        cur_cycle = nothing
        cur_state = nothing
    end

    last_print_time = time()
    num_equations = 0.0
    print_every_n_seconds = 5
    equation_speed = Float32[]

    if parallel
        for i=1:options.npopulations
            # Start listening for each population to finish:
            @async put!(channels[i], fetch(allPops[i]))
        end
    end

    while cycles_remaining > 0
        for i=1:options.npopulations
            # Non-blocking check if a population is ready:
            population_ready = parallel ? isready(channels[i]) : true
            if population_ready
                # Take the fetch operation from the channel since its ready
                (cur_pop, best_seen) = parallel ? take!(channels[i]) : allPops[i]
                cur_pop::Population
                best_seen::HallOfFame
                bestSubPops[i] = bestSubPop(cur_pop, topn=options.topn)

                #Try normal copy...
                bestPops = Population([member for pop in bestSubPops for member in pop.members])

                for (i_member, member) in enumerate(Iterators.flatten((cur_pop.members, best_seen.members[best_seen.exists])))
                    part_of_cur_pop = i_member <= length(cur_pop.members)
                    size = countNodes(member.tree)
                    if part_of_cur_pop
                        frequencyComplexity[size] += 1
                    end
                    actualMaxsize = options.maxsize + maxdegree
                    if size < actualMaxsize && all([member.score < hallOfFame.members[size2].score for size2=1:size])
                        hallOfFame.members[size] = copyPopMember(member)
                        hallOfFame.exists[size] = true
                    end
                end

                # Dominating pareto curve - must be better than all simpler equations
                dominating = calculateParetoFrontier(dataset, hallOfFame, options)
                open(options.hofFile, "w") do io
                    println(io,"Complexity|MSE|Equation")
                    for member in dominating
                        println(io, "$(countNodes(member.tree))|$(member.score)|$(stringTree(member.tree, options, varMap=dataset.varMap))")
                    end
                end
                cp(options.hofFile, options.hofFile*".bkup", force=true)

                # Try normal copy otherwise.
                if options.migration
                    for k in rand(1:options.npop, round(Int, options.npop*options.fractionReplaced))
                        to_copy = rand(1:size(bestPops.members, 1))
                        cur_pop.members[k] = PopMember(
                            copyNode(bestPops.members[to_copy].tree),
                            bestPops.members[to_copy].score)
                    end
                end

                if options.hofMigration && size(dominating, 1) > 0
                    for k in rand(1:options.npop, round(Int, options.npop*options.fractionReplacedHof))
                        # Copy in case one gets used twice
                        to_copy = rand(1:size(dominating, 1))
                        cur_pop.members[k] = PopMember(
                           copyNode(dominating[to_copy].tree), dominating[to_copy].score
                        )
                    end
                end

                cycles_remaining -= 1
                if cycles_remaining == 0
                    break
                end
                worker_idx = next_worker()
                allPops[i] = if parallel
                    @spawnat worker_idx let
                        tmp_pop, tmp_best_seen = SRCycle(
                            dataset, baselineMSE, cur_pop, options.ncyclesperiteration,
                            curmaxsize, copy(frequencyComplexity)/sum(frequencyComplexity),
                            verbosity=options.verbosity, options=options)
                        tmp_pop = OptimizeAndSimplifyPopulation(dataset, baselineMSE, tmp_pop, options, curmaxsize)
                        if options.batching
                            for i_member=1:(options.maxsize + maxdegree)
                                tmp_best_seen.members[i_member].score = scoreFunc(dataset, baselineMSE, tmp_best_seen.members[i_member].tree, options)
                            end
                        end
                        (tmp_pop, tmp_best_seen)
                    end
                else
                    tmp_pop, tmp_best_seen = SRCycle(
                        dataset, baselineMSE, cur_pop, options.ncyclesperiteration,
                        curmaxsize, copy(frequencyComplexity)/sum(frequencyComplexity),
                        verbosity=options.verbosity, options=options)
                    tmp_pop = OptimizeAndSimplifyPopulation(dataset, baselineMSE, tmp_pop, options, curmaxsize)
                    if options.batching
                        for i_member=1:(options.maxsize + maxdegree)
                            tmp_best_seen.members[i_member].score = scoreFunc(dataset, baselineMSE, tmp_best_seen.members[i_member].tree, options)
                        end
                    end
                    (tmp_pop, tmp_best_seen)
                end
                if parallel
                    @async put!(channels[i], fetch(allPops[i]))
                end

                cycles_elapsed = total_cycles - cycles_remaining
                if options.warmupMaxsizeBy > 0
                    fraction_elapsed = 1f0 * cycles_elapsed / total_cycles
                    if fraction_elapsed > options.warmupMaxsizeBy
                        curmaxsize = options.maxsize
                    else
                        curmaxsize = 3 + floor(Int, (options.maxsize - 3) * fraction_elapsed / options.warmupMaxsizeBy)
                    end
                end
                num_equations += options.ncyclesperiteration * options.npop / 10.0

                if options.progress
                    # set_postfix(iter, Equations=)
                    equation_strings = string_dominating_pareto_curve(hallOfFame, baselineMSE,
                                                                      dataset, options,
                                                                      avgy)
                    set_multiline_postfix(progress_bar, equation_strings)
                    if cur_cycle == nothing
                        (cur_cycle, cur_state) = iterate(progress_bar)
                    else
                        (cur_cycle, cur_state) = iterate(progress_bar, cur_state)
                    end
                end
            end
            elapsed = time() - last_print_time
            #Update if time has passed, and some new equations generated.
            if elapsed > print_every_n_seconds && num_equations > 0.0
                # Dominating pareto curve - must be better than all simpler equations
                current_speed = num_equations/elapsed
                average_over_m_measurements = 10 #for print_every...=5, this gives 50 second running average
                push!(equation_speed, current_speed)
                if length(equation_speed) > average_over_m_measurements
                    deleteat!(equation_speed, 1)
                end
                if options.verbosity > 0
                    @printf("\n")
                    average_speed = sum(equation_speed)/length(equation_speed)
                    @printf("Cycles per second: %.3e\n", round(average_speed, sigdigits=3))
                    cycles_elapsed = total_cycles - cycles_remaining
                    @printf("Progress: %d / %d total iterations (%.3f%%)\n",
                            cycles_elapsed, total_cycles,
                            100.0*cycles_elapsed/total_cycles)
                    equation_strings = string_dominating_pareto_curve(hallOfFame, baselineMSE,
                                                                      dataset, options,
                                                                      avgy)
                    print(equation_strings)
                end
                last_print_time = time()
                num_equations = 0.0
            end
        end
        sleep(1e-3)
    end
    if we_created_procs
        rmprocs(procs)
    end
    ##########################################################################
    ### Distributed code^
    ##########################################################################
    return hallOfFame
end

function EquationSearch(X::AbstractMatrix{T1}, y::AbstractVector{T2}; kw...) where {T1<:Real,T2<:Real}
    U = promote_type(T1, T2)
    EquationSearch(convert(AbstractMatrix{U}, X), convert(AbstractVector{U}, y); kw...)
end

end #module SR


text:
@adjoint Base.reduce(::typeof(hcat), V::AbstractVector{<:AbstractVector}) =
    reduce(hcat, V), dV -> (nothing, collect(eachcol(dV)),)


text:
 function j(n,k)
    p, i, seq=[0:n-1], 0, Int[]
    while !isempty(p)
        i=(i+k-1)%length(p)
        push!(seq,splice!(p,i+1))
    end
    @sprintf("Prisoner killing order: %s.\nSurvivor: %i",replace(chomp(string(seq[1:end-1])),"\n",", "),seq[end])
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Stores the OperationProblem definitions to be used in the simulation. When creating the
SimulationProblems, the order in which the problems are created determines the order on which
the simulation is executed.
"""
mutable struct SimulationProblems
    op_problems::OrderedDict{Symbol, OperationsProblem}
    names::Vector{Symbol}
    function SimulationProblems(; kwargs...)
        prob_dict = OrderedDict(kwargs...)
        new(prob_dict, collect(keys(prob_dict)))
    end
end

Base.getindex(problems::SimulationProblems, key::Symbol) = problems.op_problems[key]
Base.getindex(problems::SimulationProblems, key) = getindex(problems, Symbol(key))

Base.length(problems::SimulationProblems) = length(problems.op_problems)
Base.first(problems::SimulationProblems) = first(problems.op_problems)
Base.iterate(problems::SimulationProblems, args...) = iterate(problems.op_problems, args...)

get_problem_names(problems::SimulationProblems) = problems.names

function get_problem_number(problems::SimulationProblems, name)
    return findfirst(x -> x == Symbol(name), get_problem_names(problems))
end

function determine_horizons!(problems::SimulationProblems)
    horizons = OrderedDict{Symbol, Int}()
    for (name, problem) in problems.op_problems
        optimization_container = get_optimization_container(problem)
        settings = get_settings(optimization_container)
        horizon = get_horizon(settings)
        if horizon == UNSET_HORIZON
            sys = get_system(problem)
            horizon = PSY.get_forecast_horizon(sys)
            set_horizon!(settings, horizon)
        end
        horizons[name] = horizon
    end
    return horizons
end

function determine_step_resolution(intervals)
    return first(intervals)[2][1]
end

function initialize_simulation_internals!(problems::SimulationProblems, uuid::Base.UUID)
    for (ix, (name, problem)) in enumerate(problems.op_problems)
        info = SimulationInfo(
            ix,
            # JDNOTE: Making conversion to avoid breaking things
            String(name),
            0,
            0,
            Set{CacheKey}(),
            0,
            Dict{Int, FeedForwardChronology}(),
            false,
            uuid,
        )
        set_simulation_info!(problem, info)
        settings = get_settings(problem)
        set_use_parameters!(settings, true)
    end
end


text:
using CSyntax
using Test

include("for.jl")
include("switch.jl")
include("ref.jl")
include("static.jl")


text:
function get_date(obj::PosVelChebyshev)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_position_velocity_acceleration(obj::PosVelChebyshev, arg0::AbsoluteDate)
    return jcall(obj, "getPositionVelocityAcceleration", PVCoordinates, (AbsoluteDate,), arg0)
end

function get_position_velocity_acceleration(obj::PosVelChebyshev, arg0::FieldAbsoluteDate)
    return jcall(obj, "getPositionVelocityAcceleration", FieldPVCoordinates, (FieldAbsoluteDate,), arg0)
end

function in_range(obj::PosVelChebyshev, arg0::AbsoluteDate)
    return jcall(obj, "inRange", jboolean, (AbsoluteDate,), arg0)
end