text:
@testset "a2af" begin
    @test a2af(4, 2.345) == ('+', 134, 21, 30, 9706)
end

@testset "a2tf" begin
    @test a2tf(4, -3.01234) == ('-', 11, 30, 22, 6484)
end

@testset "ab" begin
    pnat = [-0.76321968546737951,-0.60869453983060384,-0.21676408580639883]
    v = [2.1044018893653786e-5,-8.9108923304429319e-5,-3.8633714797716569e-5]
    s = 0.99980921395708788
    bm1 = 0.99999999506209258
    ppr = ab(pnat, v, s, bm1)
    @test ppr[1] ≈ -0.7631631094219556269 atol=1e-12
    @test ppr[2] ≈ -0.6087553082505590832 atol=1e-12
    @test ppr[3] ≈ -0.2167926269368471279 atol=1e-12
    @test_throws ArgumentError ab(pnat[1:2], v, s, bm1)
    @test_throws ArgumentError ab(pnat, v[1:2], s, bm1)
end

@testset "ae2hd" begin
   a = 5.5
   e = 1.1
   p = 0.7

   h, d = ERFA.ae2hd(a, e, p)

   @test h ≈ 0.5933291115507309663 atol=1e-14
   @test d ≈ 0.9613934761647817620 atol=1e-14
end

@testset "af2a" begin
    r = af2a('-', 45, 13, 27.2)
    @test r ≈ -0.7893115794313644842 atol=1e-15
    r = af2a('+', 45, 13, 27.2)
    @test r ≈ 0.7893115794313644842 atol=1e-15
    @test_throws ERFAException af2a('+', 361, 13, 27.2)
    @test_throws ERFAException af2a('+', 45, 60, 27.2)
    @test_throws ERFAException af2a('+', 45, 13, 60.2)
end

@testset "anp" begin
    r = ERFA._anp(-0.1)
    @test r ≈ 6.183185307179586477 atol=1e-15
    @test r ≈ mod2pi(-0.1) atol=1-15
end

@testset "anpm" begin
    r = anpm(-4.0)
    @test r ≈ 2.283185307179586477 atol=1e-15
end

@testset "apcg" begin
    date1 = 2456165.5
    date2 = 0.401182685
    ebpv = [[0.901310875,-0.417402664,-0.180982288],
            [0.00742727954,0.0140507459,0.00609045792]]
    ehp = [0.903358544,-0.415395237,-0.180084014]
    astrom = apcg(date1, date2, ebpv, ehp)
    @test astrom.pmt ≈ 12.65133794027378508 atol=1e-11
    @test astrom.eb[1] ≈ 0.901310875 atol=1e-12
    @test astrom.eb[2] ≈ -0.417402664 atol=1e-12
    @test astrom.eb[3] ≈ -0.180982288 atol=1e-12
    @test astrom.eh[1] ≈ 0.8940025429324143045 atol=1e-12
    @test astrom.eh[2] ≈ -0.4110930268679817955 atol=1e-12
    @test astrom.eh[3] ≈ -0.1782189004872870264 atol=1e-12
    @test astrom.em ≈ 1.010465295811013146 atol=1e-12
    @test astrom.v[1] ≈ 0.4289638913597693554e-4 atol=1e-16
    @test astrom.v[2] ≈ 0.8115034051581320575e-4 atol=1e-16
    @test astrom.v[3] ≈ 0.3517555136380563427e-4 atol=1e-16
    @test astrom.bm1 ≈ 0.9999999951686012981 atol=1e-12
    @test astrom.bpn[1,1] ≈ 1.0 atol=1e-10
    @test astrom.bpn[2,1] ≈ 0.0 atol=1e-10
    @test astrom.bpn[3,1] ≈ 0.0 atol=1e-10
    @test astrom.bpn[1,2] ≈ 0.0 atol=1e-10
    @test astrom.bpn[2,2] ≈ 1.0 atol=1e-10
    @test astrom.bpn[3,2] ≈ 0.0 atol=1e-10
    @test astrom.bpn[1,3] ≈ 0.0 atol=1e-10
    @test astrom.bpn[2,3] ≈ 0.0 atol=1e-10
    @test astrom.bpn[3,3] ≈ 1.0 atol=1e-10
    ebpve = [[-0.417402664,-0.180982288],
            [0.00742727954,0.0140507459,0.00609045792]]
    @test_throws ArgumentError apcg(date1, date2, ebpve, ehp)
    @test_throws ArgumentError apcg(date1, date2, ebpv, ehp[1:2])
end

@testset "apcg13" begin
    date1 = 2456165.5
    date2 = 0.401182685
    astrom = apcg13(date1, date2)
    @test astrom.pmt ≈ 12.65133794027378508 atol=1e-12
    @test astrom.eb[1] ≈ 0.9013108747340644755 atol=1e-12
    @test astrom.eb[2] ≈ -0.4174026640406119957 atol=1e-12
    @test astrom.eb[3] ≈ -0.1809822877867817771 atol=1e-12
    @test astrom.eh[1] ≈ 0.8940025429255499549 atol=1e-12
    @test astrom.eh[2] ≈ -0.4110930268331896318 atol=1e-12
    @test astrom.eh[3] ≈ -0.1782189006019749850 atol=1e-12
    @test astrom.em ≈ 1.010465295964664178 atol=1e-12
    @test astrom.v[1] ≈ 0.4289638912941341125e-4 atol=1e-16
    @test astrom.v[2] ≈ 0.8115034032405042132e-4 atol=1e-16
    @test astrom.v[3] ≈ 0.3517555135536470279e-4 atol=1e-16
    @test astrom.bm1 ≈ 0.9999999951686013142 atol=1e-12
    @test astrom.bpn[1,1] ≈ 1.0 atol=1e-10
    @test astrom.bpn[2,1] ≈ 0.0 atol=1e-10
    @test astrom.bpn[3,1] ≈ 0.0 atol=1e-10
    @test astrom.bpn[1,2] ≈ 0.0 atol=1e-10
    @test astrom.bpn[2,2] ≈ 1.0 atol=1e-10
    @test astrom.bpn[3,2] ≈ 0.0 atol=1e-10
    @test astrom.bpn[1,3] ≈ 0.0 atol=1e-10
    @test astrom.bpn[2,3] ≈ 0.0 atol=1e-10
    @test astrom.bpn[3,3] ≈ 1.0 atol=1e-10
end

@testset "apci" begin
    date1 = 2456165.5
    date2 = 0.401182685
    ebpv = [[0.901310875,-0.417402664,-0.180982288],
            [0.00742727954,0.0140507459,0.00609045792]]
    ehp = [0.903358544,-0.415395237,-0.180084014]
    x =  0.0013122272
    y = -2.92808623e-5
    s =  3.05749468e-8
    astrom = apci(date1, date2, ebpv, ehp, x, y, s)
    @test astrom.pmt ≈ 12.65133794027378508 atol=1e-11
    @test astrom.eb[1] ≈ 0.901310875 atol=1e-12
    @test astrom.eb[2] ≈ -0.417402664 atol=1e-12
    @test astrom.eb[3] ≈ -0.180982288 atol=1e-12
    @test astrom.eh[1] ≈ 0.8940025429324143045 atol=1e-12
    @test astrom.eh[2] ≈ -0.4110930268679817955 atol=1e-12
    @test astrom.eh[3] ≈ -0.1782189004872870264 atol=1e-12
    @test astrom.em ≈ 1.010465295811013146 atol=1e-12
    @test astrom.v[1] ≈ 0.4289638913597693554e-4 atol=1e-16
    @test astrom.v[2] ≈ 0.8115034051581320575e-4 atol=1e-16
    @test astrom.v[3] ≈ 0.3517555136380563427e-4 atol=1e-16
    @test astrom.bm1 ≈ 0.9999999951686012981 atol=1e-12
    @test astrom.bpn[1,1] ≈ 0.9999991390295159156 atol=1e-12
    @test astrom.bpn[2,1] ≈ 0.4978650072505016932e-7 atol=1e-12
    @test astrom.bpn[3,1] ≈ 0.1312227200000000000e-2 atol=1e-12
    @test astrom.bpn[1,2] ≈ -0.1136336653771609630e-7 atol=1e-12
    @test astrom.bpn[2,2] ≈ 0.9999999995713154868 atol=1e-12
    @test astrom.bpn[3,2] ≈ -0.2928086230000000000e-4 atol=1e-12
    @test astrom.bpn[1,3] ≈ -0.1312227200895260194e-2 atol=1e-12
    @test astrom.bpn[2,3] ≈ 0.2928082217872315680e-4 atol=1e-12
    @test astrom.bpn[3,3] ≈ 0.9999991386008323373 atol=1e-12
    ebpve = [[-0.417402664,-0.180982288],
            [0.00742727954,0.0140507459,0.00609045792]]
    @test_throws ArgumentError apci(date1, date2, ebpve, ehp, x, y, s)
    @test_throws ArgumentError apci(date1, date2, ebpv, ehp[1:2], x, y, s)
end

@testset "apci13" begin
    date1 = 2456165.5
    date2 = 0.401182685
    astrom, eo = apci13(date1, date2)
    @test astrom.pmt ≈ 12.65133794027378508 atol=1e-11
    @test astrom.eb[1] ≈ 0.9013108747340644755 atol=1e-12
    @test astrom.eb[2] ≈ -0.4174026640406119957 atol=1e-12
    @test astrom.eb[3] ≈ -0.1809822877867817771 atol=1e-12
    @test astrom.eh[1] ≈ 0.8940025429255499549 atol=1e-12
    @test astrom.eh[2] ≈ -0.4110930268331896318 atol=1e-12
    @test astrom.eh[3] ≈ -0.1782189006019749850 atol=1e-12
    @test astrom.em ≈ 1.010465295964664178 atol=1e-12
    @test astrom.v[1] ≈ 0.4289638912941341125e-4 atol=1e-16
    @test astrom.v[2] ≈ 0.8115034032405042132e-4 atol=1e-16
    @test astrom.v[3] ≈ 0.3517555135536470279e-4 atol=1e-16
    @test astrom.bm1 ≈ 0.9999999951686013142 atol=1e-12
    @test astrom.bpn[1,1] ≈ 0.9999992060376761710 atol=1e-12
    @test astrom.bpn[2,1] ≈ 0.4124244860106037157e-7 atol=1e-12
    @test astrom.bpn[3,1] ≈ 0.1260128571051709670e-2 atol=1e-12
    @test astrom.bpn[1,2] ≈ -0.1282291987222130690e-7 atol=1e-12
    @test astrom.bpn[2,2] ≈ 0.9999999997456835325 atol=1e-12
    @test astrom.bpn[3,2] ≈ -0.2255288829420524935e-4 atol=1e-12
    @test astrom.bpn[1,3] ≈ -0.1260128571661374559e-2 atol=1e-12
    @test astrom.bpn[2,3] ≈ 0.2255285422953395494e-4 atol=1e-12
    @test astrom.bpn[3,3] ≈ 0.9999992057833604343 atol=1e-12
    @test eo ≈ -0.2900618712657375647e-2 atol=1e-12
end

@testset "apco" begin
    date1 = 2456384.5
    date2 = 0.970031644
    ebpv = [[-0.974170438,-0.211520082,-0.0917583024],
            [0.00364365824,-0.0154287319,-0.00668922024]]
    ehp = [-0.973458265,-0.209215307,-0.0906996477]
    x = 0.0013122272
    y = -2.92808623e-5
    s = 3.05749468e-8
    theta = 3.14540971
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    sp = -3.01974337e-11
    refa = 0.000201418779
    refb = -2.36140831e-7
    astrom = apco(date1, date2, ebpv, ehp, x, y, s, theta, elong, phi, hm, xp, yp, sp, refa, refb)
    @test astrom.pmt ≈ 13.25248468622587269 atol=1e-11
    @test astrom.eb[1] ≈ -0.9741827110630897003 atol=1e-12
    @test astrom.eb[2] ≈ -0.2115130190135014340 atol=1e-12
    @test astrom.eb[3] ≈ -0.09179840186968295686 atol=1e-12
    @test astrom.eh[1] ≈ -0.9736425571689670428 atol=1e-12
    @test astrom.eh[2] ≈ -0.2092452125848862201 atol=1e-12
    @test astrom.eh[3] ≈ -0.09075578152261439954 atol=1e-12
    @test astrom.em ≈ 0.9998233241710617934 atol=1e-12
    @test astrom.v[1] ≈ 0.2078704992916728762e-4 atol=1e-16
    @test astrom.v[2] ≈ -0.8955360107151952319e-4 atol=1e-16
    @test astrom.v[3] ≈ -0.3863338994288951082e-4 atol=1e-16
    @test astrom.bm1 ≈ 0.9999999950277561236 atol=1e-12
    @test astrom.bpn[1,1] ≈ 0.9999991390295159156 atol=1e-12
    @test astrom.bpn[2,1] ≈ 0.4978650072505016932e-7 atol=1e-12
    @test astrom.bpn[3,1] ≈ 0.1312227200000000000e-2 atol=1e-12
    @test astrom.bpn[1,2] ≈ -0.1136336653771609630e-7 atol=1e-12
    @test astrom.bpn[2,2] ≈ 0.9999999995713154868 atol=1e-12
    @test astrom.bpn[3,2] ≈ -0.2928086230000000000e-4 atol=1e-12
    @test astrom.bpn[1,3] ≈ -0.1312227200895260194e-2 atol=1e-12
    @test astrom.bpn[2,3] ≈ 0.2928082217872315680e-4 atol=1e-12
    @test astrom.bpn[3,3] ≈ 0.9999991386008323373 atol=1e-12
    @test astrom.along ≈ -0.5278008060295995734 atol=1e-12
    @test astrom.xpl ≈ 0.1133427418130752958e-5 atol=1e-17
    @test astrom.ypl ≈ 0.1453347595780646207e-5 atol=1e-17
    @test astrom.sphi ≈ -0.9440115679003211329 atol=1e-12
    @test astrom.cphi ≈ 0.3299123514971474711 atol=1e-12
    @test astrom.diurab ≈ 0 atol=1e-10
    @test astrom.eral ≈ 2.617608903970400427 atol=1e-12
    @test astrom.refa ≈ 0.2014187790000000000e-3 atol=1e-15
    @test astrom.refb ≈ -0.2361408310000000000e-6 atol=1e-18
    ebpve = [[-0.211520082,-0.0917583024],
            [0.00364365824,-0.0154287319,-0.00668922024]]
    @test_throws ArgumentError apco(date1, date2, ebpve, ehp, x, y, s, theta, elong,
                                         phi, hm, xp, yp, sp, refa, refb)
    @test_throws ArgumentError apco(date1, date2, ebpve, ehp[1:2], x, y, s, theta,
                                         elong, phi, hm, xp, yp, sp, refa, refb)
end

@testset "apco13" begin
    utc1 = 2456384.5
    utc2 = 0.969254051
    dut1 = 0.1550675
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    p = 2.47230737e-7
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    phpa = 731.0
    tc = 12.8
    rh = 0.59
    wl = 0.55
    astrom, eo = apco13(utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test astrom.pmt ≈ 13.25248468622475727 atol=1e-11
    @test astrom.eb[1] ≈ -0.9741827107321449445 atol=1e-12
    @test astrom.eb[2] ≈ -0.2115130190489386190 atol=1e-12
    @test astrom.eb[3] ≈ -0.09179840189515518726 atol=1e-12
    @test astrom.eh[1] ≈ -0.9736425572586866640 atol=1e-12
    @test astrom.eh[2] ≈ -0.2092452121602867431 atol=1e-12
    @test astrom.eh[3] ≈ -0.09075578153903832650 atol=1e-12
    @test astrom.em ≈ 0.9998233240914558422 atol=1e-12
    @test astrom.v[1] ≈ 0.2078704994520489246e-4 atol=1e-16
    @test astrom.v[2] ≈ -0.8955360133238868938e-4 atol=1e-16
    @test astrom.v[3] ≈ -0.3863338993055887398e-4 atol=1e-16
    @test astrom.bm1 ≈ 0.9999999950277561004 atol=1e-12
    @test astrom.bpn[1,1] ≈ 0.9999991390295147999 atol=1e-12
    @test astrom.bpn[2,1] ≈ 0.4978650075315529277e-7 atol=1e-12
    @test astrom.bpn[3,1] ≈ 0.001312227200850293372 atol=1e-12
    @test astrom.bpn[1,2] ≈ -0.1136336652812486604e-7 atol=1e-12
    @test astrom.bpn[2,2] ≈ 0.9999999995713154865 atol=1e-12
    @test astrom.bpn[3,2] ≈ -0.2928086230975367296e-4 atol=1e-12
    @test astrom.bpn[1,3] ≈ -0.001312227201745553566 atol=1e-12
    @test astrom.bpn[2,3] ≈ 0.2928082218847679162e-4 atol=1e-12
    @test astrom.bpn[3,3] ≈ 0.9999991386008312212 atol=1e-12
    @test astrom.along ≈ -0.5278008060295995733 atol=1e-12
    @test astrom.xpl ≈ 0.1133427418130752958e-5 atol=1e-17
    @test astrom.ypl ≈ 0.1453347595780646207e-5 atol=1e-17
    @test astrom.sphi ≈ -0.9440115679003211329 atol=1e-12
    @test astrom.cphi ≈ 0.3299123514971474711 atol=1e-12
    @test astrom.diurab ≈ 0 atol=1e-10
    @test astrom.eral ≈ 2.617608909189664000 atol=1e-12
    @test astrom.refa ≈ 0.2014187785940396921e-3 atol=1e-15
    @test astrom.refb ≈ -0.2361408314943696227e-6 atol=1e-18
    @test eo ≈ -0.003020548354802412839 atol=1e-14
    @test_logs (:warn,) apco13(2.5245935e6, 0.0, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test_throws ERFAException apco13(-1e9, 0.0, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl)
end

@testset "apcs" begin
    date1 = 2456384.5
    date2 = 0.970031644
    pv = [[-1836024.09,1056607.72,-5998795.26],
          [-77.0361767,-133.310856,0.0971855934]]
    ebpv = [[-0.974170438,-0.211520082,-0.0917583024],
            [0.00364365824,-0.0154287319,-0.00668922024]]
    ehp = [-0.973458265,-0.209215307,-0.0906996477]
    astrom = apcs(date1, date2, pv, ebpv, ehp)
    @test astrom.pmt ≈ 13.25248468622587269 atol=1e-11
    @test astrom.eb[1] ≈ -0.9741827110630456169 atol=1e-12
    @test astrom.eb[2] ≈ -0.2115130190136085494 atol=1e-12
    @test astrom.eb[3] ≈ -0.09179840186973175487 atol=1e-12
    @test astrom.eh[1] ≈ -0.9736425571689386099 atol=1e-12
    @test astrom.eh[2] ≈ -0.2092452125849967195 atol=1e-12
    @test astrom.eh[3] ≈ -0.09075578152266466572 atol=1e-12
    @test astrom.em ≈ 0.9998233241710457140 atol=1e-12
    @test astrom.v[1] ≈ 0.2078704993282685510e-4 atol=1e-16
    @test astrom.v[2] ≈ -0.8955360106989405683e-4 atol=1e-16
    @test astrom.v[3] ≈ -0.3863338994289409097e-4 atol=1e-16
    @test astrom.bm1 ≈ 0.9999999950277561237 atol=1e-12
    @test astrom.bpn[1,1] ≈ 1 atol=1e-10
    @test astrom.bpn[2,1] ≈ 0 atol=1e-10
    @test astrom.bpn[3,1] ≈ 0 atol=1e-10
    @test astrom.bpn[1,2] ≈ 0 atol=1e-10
    @test astrom.bpn[2,2] ≈ 1 atol=1e-10
    @test astrom.bpn[3,2] ≈ 0 atol=1e-10
    @test astrom.bpn[1,3] ≈ 0 atol=1e-10
    @test astrom.bpn[2,3] ≈ 0 atol=1e-10
    @test astrom.bpn[3,3] ≈ 1 atol=1e-10
    ebpve = [[-0.211520082,-0.0917583024],
            [0.00364365824,-0.0154287319,-0.00668922024]]
    @test_throws ArgumentError apcs(date1, date2, ebpve, ebpv, ehp)
    @test_throws ArgumentError apcs(date1, date2, pv, ebpve, ehp)
    @test_throws ArgumentError apcs(date1, date2, pv, ebpv, ehp[1:2])
end

@testset "apcs13" begin
    date1 = 2456165.5
    date2 = 0.401182685
    pv = [[-6241497.16,401346.896,-1251136.04],
          [-29.264597,-455.021831,0.0266151194]]
    astrom = apcs13(date1, date2, pv)
    @test astrom.pmt ≈ 12.65133794027378508 atol=1e-11
    @test astrom.eb[1] ≈ 0.9012691529023298391 atol=1e-12
    @test astrom.eb[2] ≈ -0.4173999812023068781 atol=1e-12
    @test astrom.eb[3] ≈ -0.1809906511146821008 atol=1e-12
    @test astrom.eh[1] ≈ 0.8939939101759726824 atol=1e-12
    @test astrom.eh[2] ≈ -0.4111053891734599955 atol=1e-12
    @test astrom.eh[3] ≈ -0.1782336880637689334 atol=1e-12
    @test astrom.em ≈ 1.010428384373318379 atol=1e-12
    @test astrom.v[1] ≈ 0.4279877294121697570e-4 atol=1e-16
    @test astrom.v[2] ≈ 0.7963255087052120678e-4 atol=1e-16
    @test astrom.v[3] ≈ 0.3517564013384691531e-4 atol=1e-16
    @test astrom.bm1 ≈ 0.9999999952947980978 atol=1e-12
    @test astrom.bpn[1,1] ≈ 1 atol=1e-10
    @test astrom.bpn[2,1] ≈ 0 atol=1e-10
    @test astrom.bpn[3,1] ≈ 0 atol=1e-10
    @test astrom.bpn[1,2] ≈ 0 atol=1e-10
    @test astrom.bpn[2,2] ≈ 1 atol=1e-10
    @test astrom.bpn[3,2] ≈ 0 atol=1e-10
    @test astrom.bpn[1,3] ≈ 0 atol=1e-10
    @test astrom.bpn[2,3] ≈ 0 atol=1e-10
    @test astrom.bpn[3,3] ≈ 1 atol=1e-10
    pve = [[401346.896,-1251136.04],
          [-29.264597,-455.021831,0.0266151194]]
    @test_throws ArgumentError apcs13(date1, date2, pve)
end

@testset "aper" begin
    theta = 5.678
    pmt = 0.
    eb = zeros(Cdouble, 3)
    eh = zeros(Cdouble, 3)
    em = 0.
    v = zeros(Cdouble, 3)
    bm1 = 0.
    bpn = zeros(Cdouble, 9)
    along = 1.234
    phi, xpl, ypl, sphi, cphi, diurab, eral, refa, refb = 0., 0., 0., 0., 0., 0., 0., 0., 0.
    astrom = ERFA.ASTROM(pmt, eb, eh, em, v, bm1, bpn, along,
                         phi, xpl, ypl, sphi, cphi, diurab, eral, refa, refb)
    astrom = aper(theta, astrom)
    @test astrom.eral ≈ 6.912000000000000000 atol=1e-12
end

@testset "aper13" begin
    ut11 = 2456165.5
    ut12 = 0.401182685
    pmt = 0.
    eb = zeros(Cdouble, 3)
    eh = zeros(Cdouble, 3)
    em = 0.
    v = zeros(Cdouble, 3)
    bm1 = 0.
    bpn = zeros(Cdouble, 9)
    along = 1.234
    phi, xpl, ypl, sphi, cphi, diurab, eral, refa, refb = 0., 0., 0., 0., 0., 0., 0., 0., 0.
    astrom = ERFA.ASTROM(pmt, eb, eh, em, v, bm1, bpn, along,
                         phi, xpl, ypl, sphi, cphi, diurab, eral, refa, refb)
    astrom = aper13(ut11, ut12, astrom)
    @test astrom.eral ≈ 3.316236661789694933 atol=1e-12
end

@testset "apio" begin
    sp = -3.01974337e-11
    theta = 3.14540971
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    refa = 0.000201418779
    refb = -2.36140831e-7
    astrom = apio(sp, theta, elong, phi, hm, xp, yp, refa, refb)
    @test astrom.along ≈ -0.5278008060295995734 atol=1e-12
    @test astrom.xpl ≈ 0.1133427418130752958e-5 atol=1e-17
    @test astrom.ypl ≈ 0.1453347595780646207e-5 atol=1e-17
    @test astrom.sphi ≈ -0.9440115679003211329 atol=1e-12
    @test astrom.cphi ≈ 0.3299123514971474711 atol=1e-12
    @test astrom.diurab ≈ 0.5135843661699913529e-6 atol=1e-12
    @test astrom.eral ≈ 2.617608903970400427 atol=1e-12
    @test astrom.refa ≈ 0.2014187790000000000e-3 atol=1e-15
    @test astrom.refb ≈ -0.2361408310000000000e-6 atol=1e-18
end

@testset "apio13" begin
    utc1 = 2456384.5
    utc2 = 0.969254051
    dut1 = 0.1550675
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    phpa = 731.0
    tc = 12.8
    rh = 0.59
    wl = 0.55
    astrom = apio13(utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test astrom.along ≈ -0.5278008060295995733 atol=1e-12
    @test astrom.xpl ≈ 0.1133427418130752958e-5 atol=1e-17
    @test astrom.ypl ≈ 0.1453347595780646207e-5 atol=1e-17
    @test astrom.sphi ≈ -0.9440115679003211329 atol=1e-12
    @test astrom.cphi ≈ 0.3299123514971474711 atol=1e-12
    @test astrom.diurab ≈ 0.5135843661699913529e-6 atol=1e-12
    @test astrom.eral ≈ 2.617608909189664000 atol=1e-12
    @test astrom.refa ≈ 0.2014187785940396921e-3 atol=1e-15
    @test astrom.refb ≈ -0.2361408314943696227e-6 atol=1e-18
    @test_logs (:warn,) apio13(2.5245935e6, 0.0, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test_throws ERFAException apio13(-1e9, 0.0, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl)
end

@testset "atcc13" begin
   rc = 2.71
   dc = 0.174
   pr = 1e-5
   pd = 5e-6
   px = 0.1
   rv = 55.0
   date1 = 2456165.5
   date2 = 0.401182685

   ra, da = atcc13(rc, dc, pr, pd, px, rv, date1, date2)

   @test ra ≈ 2.710126504531372384 atol=1e-12
   @test da ≈ 0.1740632537628350152 atol=1e-12
end

@testset "atccq" begin
   date1 = 2456165.5
   date2 = 0.401182685
   astrom, eo = apci13(date1, date2)
   rc = 2.71
   dc = 0.174
   pr = 1e-5
   pd = 5e-6
   px = 0.1
   rv = 55.0

   ra, da = atccq(rc, dc, pr, pd, px, rv, astrom)

   @test ra ≈ 2.710126504531372384 atol=1e-12
   @test da ≈ 0.1740632537628350152 atol=1e-12
end

@testset "atci13" begin
    rc = 2.71
    dc = 0.174
    pr = 1e-5
    pd = 5e-6
    px = 0.1
    rv = 55.0
    date1 = 2456165.5
    date2 = 0.401182685
    ri, di, eo = atci13(rc, dc, pr, pd, px, rv, date1, date2)
    @test ri ≈ 2.710121572969038991 atol=1e-12
    @test di ≈ 0.1729371367218230438 atol=1e-12
    @test eo ≈ -0.002900618712657375647 atol=1e-14
end

@testset "atciq" begin
    date1 = 2456165.5
    date2 = 0.401182685
    astrom, eo = apci13(date1, date2)
    rc = 2.71
    dc = 0.174
    pr = 1e-5
    pd = 5e-6
    px = 0.1
    rv = 55.0
    ri, di = atciq(rc, dc, pr, pd, px, rv, astrom)
    @test ri ≈ 2.710121572969038991 atol=1e-12
    @test di ≈ 0.1729371367218230438 atol=1e-12
end

@testset "atciqn" begin
    date1 = 2456165.5
    date2 = 0.401182685
    astrom, eo = apci13(date1, date2)
    rc = 2.71
    dc = 0.174
    pr = 1e-5
    pd = 5e-6
    px = 0.1
    rv = 55.0
    b1 = ERFA.LDBODY(0.00028574, 3e-10,
                     [[-7.81014427,-5.60956681,-1.98079819];
                      [0.0030723249,-0.00406995477,-0.00181335842]])
    b2 = ERFA.LDBODY(0.00095435, 3e-9,
                     [[0.738098796, 4.63658692,1.9693136];
                      [-0.00755816922, 0.00126913722, 0.000727999001]])
    b3 = ERFA.LDBODY(1.0, 6e-6,
                     [[-0.000712174377, -0.00230478303, -0.00105865966];
                      [6.29235213e-6, -3.30888387e-7, -2.96486623e-7]])
    b = [b1; b2; b3]
    ri, di = atciqn(rc, dc, pr, pd, px, rv, astrom, b)
    @test ri ≈ 2.710122008105325582 atol=1e-12
    @test di ≈ 0.1729371916491459122 atol=1e-12
end

@testset "atciqz" begin
    date1 = 2456165.5
    date2 = 0.401182685
    astrom, eo = apci13(date1, date2)
    rc = 2.71
    dc = 0.174
    ri, di = atciqz(rc, dc, astrom)
    @test ri ≈ 2.709994899247599271 atol=1e-12
    @test di ≈ 0.1728740720983623469 atol=1e-12
end

@testset "atco13" begin
    rc = 2.71
    dc = 0.174
    pr = 1e-5
    pd = 5e-6
    px = 0.1
    rv = 55.0
    utc1 = 2456384.5
    utc2 = 0.969254051
    dut1 = 0.1550675
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    phpa = 731.0
    tc = 12.8
    rh = 0.59
    wl = 0.55
    aob, zob, hob, dob, rob, eo = atco13(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong,
                                         phi, hm, xp, yp, phpa, tc, rh, wl)
    @test aob ≈ 0.9251774485485515207e-1 atol=1e-12
    @test zob ≈ 1.407661405256499357 atol=1e-12
    @test hob ≈ -0.9265154431529724692e-1 atol=1e-12
    @test dob ≈ 0.1716626560072526200 atol=1e-12
    @test rob ≈ 2.710260453504961012 atol=1e-12
    @test eo ≈ -0.003020548354802412839 atol=1e-14
    @test_logs (:warn,) atco13(rc, dc, pr, pd, px, rv, 2.5245935e6, 0.0, dut1, elong, phi,
                               hm, xp, yp, phpa, tc, rh, wl)
    @test_throws ERFAException atco13(rc, dc, pr, pd, px, rv, -1e9, 0.0, dut1, elong, phi,
                                      hm, xp, yp, phpa, tc, rh, wl)
end

@testset "atic13" begin
    ri = 2.710121572969038991
    di = 0.1729371367218230438
    date1 = 2456165.5
    date2 = 0.401182685
    rc, dc, eo = atic13(ri, di, date1, date2)
    @test rc ≈ 2.710126504531374930 atol=1e-12
    @test dc ≈ 0.1740632537628342320 atol=1e-12
    @test eo ≈ -0.002900618712657375647 atol=1e-14
end

@testset "aticq" begin
    ri = 2.710121572969038991
    di = 0.1729371367218230438
    date1 = 2456165.5
    date2 = 0.401182685
    astrom, eo = apci13(date1, date2)
    rc, dc = aticq(ri, di, astrom)
    @test rc ≈ 2.710126504531374930 atol=1e-12
    @test dc ≈ 0.1740632537628342320 atol=1e-12
end

@testset "aticqn" begin
    date1 = 2456165.5
    date2 = 0.401182685
    astrom, eo = apci13(date1, date2)
    ri = 2.709994899247599271
    di = 0.1728740720983623469
    b1 = ERFA.LDBODY(0.00028574, 3e-10,
                     [[-7.81014427,-5.60956681,-1.98079819];
                      [0.0030723249,-0.00406995477,-0.00181335842]])
    b2 = ERFA.LDBODY(0.00095435, 3e-9,
                     [[0.738098796, 4.63658692,1.9693136];
                      [-0.00755816922, 0.00126913722, 0.000727999001]])
    b3 = ERFA.LDBODY(1.0, 6e-6,
                     [[-0.000712174377, -0.00230478303, -0.00105865966];
                      [6.29235213e-6, -3.30888387e-7, -2.96486623e-7]])
    b = [b1; b2; b3]
    rc, dc = aticqn(ri, di, astrom, b)
    @test rc ≈ 2.709999575032685412 atol=1e-12
    @test dc ≈ 0.1739999656317778034 atol=1e-12
end

@testset "atio13" begin
    ri = 2.710121572969038991
    di = 0.1729371367218230438
    utc1 = 2456384.5
    utc2 = 0.969254051
    dut1 = 0.1550675
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    phpa = 731.0
    tc = 12.8
    rh = 0.59
    wl = 0.55
    aob, zob, hob, dob, rob = atio13(ri, di, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test aob ≈ 0.9233952224895122499e-1 atol=1e-12
    @test zob ≈ 1.407758704513549991 atol=1e-12
    @test hob ≈ -0.9247619879881698140e-1 atol=1e-12
    @test dob ≈ 0.1717653435756234676 atol=1e-12
    @test rob ≈ 2.710085107988480746 atol=1e-12
    @test_logs (:warn,) atio13(ri, di, 2.5245935e6, 0.0, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test_throws ERFAException atio13(ri, di, -1e9, 0.0, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl)
end

@testset "atioq" begin
    utc1 = 2456384.5
    utc2 = 0.969254051
    dut1 = 0.1550675
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    phpa = 731.0
    tc = 12.8
    rh = 0.59
    wl = 0.55
    astrom = apio13(utc1, utc2, dut1, elong, phi, hm, xp, yp,
                         phpa, tc, rh, wl)
    ri = 2.710121572969038991
    di = 0.1729371367218230438
    aob, zob, hob, dob, rob = atioq(ri, di, astrom)
    @test aob ≈ 0.9233952224895122499e-1 atol=1e-12
    @test zob ≈ 1.407758704513549991 atol=1e-12
    @test hob ≈ -0.9247619879881698140e-1 atol=1e-12
    @test dob ≈ 0.1717653435756234676 atol=1e-12
    @test rob ≈ 2.710085107988480746 atol=1e-12
end

@testset "atoc13" begin
    utc1 = 2456384.5
    utc2 = 0.969254051
    dut1 = 0.1550675
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    phpa = 731.0
    tc = 12.8
    rh = 0.59
    wl = 0.55
    ob1 = 2.710085107986886201
    ob2 = 0.1717653435758265198
    rc, dc = atoc13("r", ob1, ob2, utc1, utc2, dut1,
                         elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test rc ≈ 2.709956744659136129 atol=1e-12
    @test dc ≈ 0.1741696500898471362 atol=1e-12
    ob1 = -0.09247619879782006106
    ob2 = 0.1717653435758265198
    rc, dc = atoc13("h", ob1, ob2, utc1, utc2, dut1,
                         elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test rc ≈ 2.709956744659734086 atol=1e-12
    @test dc ≈ 0.1741696500898471362 atol=1e-12
    ob1 = 0.09233952224794989993
    ob2 = 1.407758704513722461
    rc, dc = atoc13("a", ob1, ob2, utc1, utc2, dut1,
                         elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test rc ≈ 2.709956744659734086 atol=1e-12
    @test dc ≈ 0.1741696500898471366 atol=1e-12
    rc1, dc1 = @test_logs (:warn,) atoc13("foo", ob1, ob2, utc1, utc2, dut1,
                                          elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test rc == rc1
    @test dc == dc1
    @test_logs (:warn,) atoc13("a", ob1, ob2, 2.5245935e6, 0.0, dut1, elong, phi, hm,
                               xp, yp, phpa, tc, rh, wl)
    @test_throws ERFAException atoc13("a", ob1, ob2, -1e9, 0.0, dut1, elong, phi, hm,
                                      xp, yp, phpa, tc, rh, wl)
end

@testset "atoi13" begin
    utc1 = 2456384.5
    utc2 = 0.969254051
    dut1 = 0.1550675
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    phpa = 731.0
    tc = 12.8
    rh = 0.59
    wl = 0.55
    ob1 = 2.710085107986886201
    ob2 = 0.1717653435758265198
    ri, di = atoi13("r", ob1, ob2, utc1, utc2, dut1,
                         elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test ri ≈ 2.710121574447540810 atol=1e-12
    @test di ≈ 0.1729371839116608778 atol=1e-12
    ob1 = -0.09247619879782006106
    ob2 = 0.1717653435758265198
    ri, di = atoi13("h", ob1, ob2, utc1, utc2, dut1,
                         elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test ri ≈ 2.710121574448138676 atol=1e-12
    @test di ≈ 0.1729371839116608778 atol=1e-12
    ob1 = 0.09233952224794989993
    ob2 = 1.407758704513722461
    ri, di = atoi13("a", ob1, ob2, utc1, utc2, dut1,
                         elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test ri ≈ 2.710121574448138676 atol=1e-12
    @test di ≈ 0.1729371839116608781 atol=1e-12
    @test_logs (:warn,) atoi13("a", ob1, ob2, 2.5245935e6, 0.0, dut1,
                               elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    @test_throws ERFAException atoi13("a", ob1, ob2, -1e9, 0.0, dut1,
                                      elong, phi, hm, xp, yp, phpa, tc, rh, wl)
end

@testset "atoiq" begin
    utc1 = 2456384.5
    utc2 = 0.969254051
    dut1 = 0.1550675
    elong = -0.527800806
    phi = -1.2345856
    hm = 2738.0
    xp = 2.47230737e-7
    yp = 1.82640464e-6
    phpa = 731.0
    tc = 12.8
    rh = 0.59
    wl = 0.55
    astrom = apio13(utc1, utc2, dut1,
                         elong, phi, hm, xp, yp, phpa, tc, rh, wl)
    ob1 = 2.710085107986886201
    ob2 = 0.1717653435758265198
    ri, di = atoiq("r", ob1, ob2, astrom)
    @test ri ≈ 2.710121574447540810 atol=1e-12
    @test di ≈ 0.17293718391166087785 atol=1e-12
    ob1 = -0.09247619879782006106
    ob2 = 0.1717653435758265198
    ri, di = atoiq("h", ob1, ob2, astrom)
    @test ri ≈ 2.710121574448138676 atol=1e-12
    @test di ≈ 0.1729371839116608778 atol=1e-12
    ob1 = 0.09233952224794989993
    ob2 = 1.407758704513722461
    ri, di = atoiq("a", ob1, ob2, astrom)
    @test ri ≈ 2.710121574448138676 atol=1e-12
    @test di ≈ 0.1729371839116608781 atol=1e-12
end


text:
using ExternalCmds
using Test

@testset "ExternalCmds.jl" begin

    if Sys.isunix()

        println()

        # ------------------------------------
        # run_cmd
        token = "TESTING_RUN_CMD"
        pid, out = run_cmd(`echo $(token)`)
        @test pid isa Integer
        @test contains(out, token)
        println()

        # ------------------------------------
        # tee_run
        stdout_log = tempname()
        rm(stdout_log; force = true)
        stderr_log = stdout_log
        for (token, long) in [
                ("TESTING_SHORT_RUN", false),
                ("TESTING_LONG_RUN", true),
            ]

            pid = tee_run(`echo $(token)`;
                long,
                stdout_log, stderr_log,
                stdout_tee_ios = [stdout], stderr_tee_ios = [],
                savetime = 5.0,
            )
            @test pid isa Integer

            # wait for answer
            while true
                if isfile(stdout_log)
                    out = read(stdout_log, String)
                    if contains(out, token)
                        @test true
                        break
                    end
                end
                sleep(1.0)
            end
            println()
        end # for (token, long)

        # ------------------------------------
        token = "TESTING_RUN_BASH SINGLELINE"
        pid, out = run_bash("echo $(token)"; run_fun = run_cmd)
        @test pid isa Integer
        @test contains(out, token)
        println()

        token = "TESTING_RUN_BASH MULTILINE"
        buff_file = tempname()
        pid, out = run_bash(["echo $(token)"]; buff_file, run_fun = run_cmd, rm_buff = true)
        @test pid isa Integer
        @test contains(out, token)
        println()

        # ------------------------------------
        # clear
        rm(stdout_log; force = true)
        rm(buff_file; force = true)

    end # if Sys.isunix()

end

ext: jl


text:
### Derivative routines

export derivative_linop

## Wrapper for derivative operator

function derivative_linop(m_geom::ModelGeom)

	return LinFunctionMod2Mod(
		m -> Model(m_geom, Dz(    m.m, m_geom.dz)),
		m -> Model(m_geom, Dz_adj(m.m, m_geom.dz)))

end

## First order

function Dz(m::Array{Float32, 2}, dz::Float32)

		Dz_m = Array{Float32, 2}(undef, size(m))
		Dz_m[1:end-1, :]  = (m[2:end, :]-m[1:end-1, :])/dz
		Dz_m[end,     :] .= 0f0

		return Dz_m

end

function Dz_adj(m::Array{Float32, 2}, dz::Float32)

		Dz_m = Array{Float32, 2}(undef, size(m))
		Dz_m[1,       :] .= (m[1,       :]              )/dz
		Dz_m[2:end-1, :] .= (m[2:end-1, :]-m[1:end-2, :])/dz
		Dz_m[end,     :] .= (             -m[end-1,   :])/dz

		return -Dz_m

end

## Second order

function Dzz(u::Array{Float32, 3}, dz::Float32; acc::Int64 = 2) # Array size: np, nz, nb

		if acc == 2
			return Dzz2(u, dz)
		elseif acc == 4
			return Dzz4(u, dz)
		end

end

function Dzz2(u::Array{Float32, 3}, dz::Float32)

	Dzz_u = Array{Float32, 3}(undef, size(u))
	Dzz_u[:, 1,       :] = (                -2f0*u[:, 1,       :]+u[:, 2,     :])/dz^2f0
	Dzz_u[:, 2:end-1, :] = (u[:, 1:end-2, :]-2f0*u[:, 2:end-1, :]+u[:, 3:end, :])/dz^2f0
	Dzz_u[:, end,     :] = (u[:, end-1,   :]-2f0*u[:, end,     :]               )/dz^2f0

	return Dzz_u

end

function Dzz4(u::Array{Float32, 3}, dz::Float32)

	Dzz_u = Array{Float32, 3}(undef, size(u))
	Dzz_u[:, 1,       :] = (                                                 -5f0/2f0*u[:, 1,       :]+4f0/3f0*u[:, 2,       :]-1f0/12f0*u[:, 3,     :])/dz^2f0
	Dzz_u[:, 2,       :] = (                          4f0/3f0*u[:, 1,       :]-5f0/2f0*u[:, 2,       :]+4f0/3f0*u[:, 3,       :]-1f0/12f0*u[:, 4,     :])/dz^2f0
	Dzz_u[:, 3:end-2, :] = (-1f0/12f0*u[:, 1:end-4, :]+4f0/3f0*u[:, 2:end-3, :]-5f0/2f0*u[:, 3:end-2, :]+4f0/3f0*u[:, 4:end-1, :]-1f0/12f0*u[:, 5:end, :])/dz^2f0
	Dzz_u[:, end-1,   :] = (-1f0/12f0*u[:, end-3,   :]+4f0/3f0*u[:, end-2,   :]-5f0/2f0*u[:, end-1,   :]+4f0/3f0*u[:, end,     :]                       )/dz^2f0
	Dzz_u[:, end,     :] = (-1f0/12f0*u[:, end-2,   :]+4f0/3f0*u[:, end-1,   :]-5f0/2f0*u[:, end,     :]                                               )/dz^2f0

	return Dzz_u

end


text:
function try_sat_pure(model,V0,f!,T,result,error_val,converged)
    converged[] =  false
    try
        res = sat_pure(model,V0,f!,T)
        result[] = res
    catch e #normally, failures occur  near the critical point
        error_val[] = e
    end

    (P_sat,V_l,V_v) = result[]
    ε = abs(V_l-V_v)/(eps(typeof(V_l-V_v)))
    if ε > 8
        #if ΔV > ε then Vl and Vv are different values
        converged[] =  true
    end

    return nothing
end

"""
    sat_pure(model::EoSModel, T, V0 = x0_sat_pure(model,T))

Performs a single component saturation equilibrium calculation, at the specified temperature `T`, of one mol of pure sustance specified by `model`

Returns `(p₀, Vₗ, Vᵥ)` where `p₀` is the saturation pressure (in Pa), `Vₗ` is the liquid saturation volume (in m³) and `Vᵥ` is the vapour saturation volume (in m³).

If the calculation fails, returns  `(NaN, NaN, NaN)`

`V0` is `[log10(Vₗ₀),log10(Vᵥ₀)]` , where `Vₗ₀`  and `Vᵥ₀` are initial guesses for the liquid and vapour volumes.
"""
function sat_pure(model::EoSModel, T, V0 = x0_sat_pure(model,T))
    V_lb = lb_volume(model,SA[1.0])
    TYPE = promote_type(typeof(T),typeof(V_lb))
    nan = zero(TYPE)/zero(TYPE)
    f! = (F,x) -> Obj_Sat(model, F, T, exp10(x[1]), exp10(x[2]),V_lb)
    res0 = (nan,nan,nan)
    result = Ref(res0)
    error_val = Ref{Any}(nothing)
    converged = Ref{Bool}(false)

    try_sat_pure(model,V0,f!,T,result,error_val,converged)
    if converged[]
        return result[]
        #V0 = log10.([V_l[i],V_v[i]])
    end
    return res0
end

#=
based on:
DOI: 10.1007/s10910-007-9272-4
Journal of Mathematical Chemistry, Vol. 43, No. 4, May 2008 (© 2007)
The van der Waals equation: analytical and approximate solutions
=#
function x0_sat_pure_crit(model,T,T_c,P_c,V_c)
    _1 = one(T)
    Tr = T/T_c
    Trm1 = _1-Tr
    Trmid = sqrt(Trm1)
    P_sat0 = (_1 -4.0*(Trm1)+4.8(Trm1*Trm1))*P_c
    #c = 1/Vr
    c_l = 1.0+2.0*Trmid + 0.4*Trm1 - 0.52*Trmid*Trm1 +0.115*Trm1*Trm1
    c_v = 1.0-2.0*Trmid + 0.4*Trm1 + 0.52*Trmid*Trm1 +0.207*Trm1*Trm1
    Vl0 = (1/c_l)*V_c
    Vv0 = (1/c_v)*V_c
    #Vl = volume_compress(model,P_sat0,T,V0=Vl0)
    #Vv = volume_compress(model,P_sat0,T,V0=Vv0)
    return [log10(Vl0),log10(Vv0)]
end

function sat_pure(model::EoSModel,V0,f!,T)
    r = Solvers.nlsolve(f!, V0,LineSearch(Newton()))
    #@show typeof(r)
    #@show f!(rand(2),r.info.zero)
    Vsol = Solvers.x_sol(r)
    V_l = exp10(Vsol[1])
    V_v = exp10(Vsol[2])
    P_sat = pressure(model,V_v,T)
    return (P_sat,V_l,V_v)
end

function Obj_Sat(model::EoSModel, F, T, V_l, V_v,V_lb)
    #components = model.components
    fun(x) = eos(model, x[2], T,SA[x[1]])
    df(x)  = ForwardDiff.gradient(fun,x)
    df_l = df(SA[one(V_l*T),V_l*one(T)])
    df_v = df(SA[one(V_v),V_v*one(T)])
    (p_scale,μ_scale) = scale_sat_pure(model)
    #T̄ = T/T_scale(model)
    F[1] = (df_l[2]-df_v[2])*p_scale
    F[2] = (df_l[1]-df_v[1])*μ_scale
    return F
end
#=
function Obj_Sat(model::ABCubicModel, F, T, V_l, V_v,V_lb)
    #components = model.components
    ar(V) = a_res(model,V,T)
    function dar(V)
        a,da = Solvers.f∂f(ar,V)
        @show Vval = ForwardDiff.value(V)
        @show aval = ForwardDiff.value(da)
        @show rp = pressure(model,Vval,T)

        p = rp
        @show pval = ForwardDiff.value(p)

         Z = p*V/(R̄*T)
        ϕ = a + (Z-1) - log(Z)
        @show phival = ForwardDiff.value(ϕ)

        return p,exp(ϕ)
    end
    pl,ϕl = dar(V_l)
    pv,ϕv = dar(V_v)
    (p_scale,μ_scale) = scale_sat_pure(model)
    #T̄ = T/T_scale(model)
    F[1] = (pl-pv)*p_scale#*exp(5e-10*(V_l-V_v)^-2)*exp(1e-7*(V_l-V_lb)^-2)
    F[2] = (ϕl-ϕv)*μ_scale#*exp(5e-10*(V_l-V_v)^-2)*exp(1e-7*(V_l-V_lb)^-2)
    return F
end
=#

"""
    crit_pure(model::EoSModel,x0=nothing)

Calculates the critical point of a single component modelled by `model`.

Returns `(Tc, pc, Vc)` where `Tc` is the critical temperature (in K), `pc` is the critical pressure (in Pa) and `Vc` is the critical volume (in  m³)
"""
function crit_pure(model::EoSModel,x0=nothing)
    T̄  = T_scale(model)
    f! = (F,x) -> obj_crit(model, F, x[1]*T̄, exp10(x[2]))
    if x0 === nothing
        x0 = x0_crit_pure(model)
    end
    r  = Solvers.x_sol(Solvers.nlsolve(f!, x0))
    T_c = r[1]*T̄
    V_c = exp10(r[2])
    p_c = pressure(model, V_c, T_c)
    return (T_c, p_c, V_c)
end

function obj_crit(model::EoSModel, F, T_c, V_c)
    ∂²A∂V², ∂³A∂V³ = ∂²³f(model, V_c, T_c, SA[1.0])
    F[1] = -∂²A∂V²
    F[2] = -∂³A∂V³
    return F
end
"""
    enthalpy_vap(model::EoSModel, T)

Calculates `ΔH`, the difference between saturated vapour and liquid enthalpies at temperature `T`, in J
"""
function enthalpy_vap(model::EoSModel, T)
    (P_sat,V_l,V_v) = sat_pure(model,T)
   #= _dfl,fl =  ∂f(model,V_l,T,SA[1.0])
    _dfv,fv =  ∂f(model,V_v,T,SA[1.0])
    dVl,dTl = _dfl
    dVv,dTv = _dfv
    H_l = fl  - dVl*V_l - dTl*T
    H_v = fv  - dVv*V_v - dTv*T =#
    H_v = VT_enthalpy(model,V_v,T)
    H_l = VT_enthalpy(model,V_l,T)
    H_vap=H_v -H_l
    return H_vap
end

"""
    acentric_factor(model::EoSModel)

calculates the acentric factor using its definition:

    ω = -log10(psatᵣ) -1, at Tᵣ = 0.7
To do so, it calculates the critical temperature (using `crit_pure`) and performs a saturation calculation (with `sat_pure`)

"""
function acentric_factor(model)
    T_c,p_c,_ = crit_pure(model)
    T = 0.7*T_c
    p = first(sat_pure(model,T))
    p_r = p/p_c
    return -log10(p_r) - 1.0
end

function p_sat_t0(model,p)
    #ln(pr) = h*(1-1/tr)
    #1/(1-ln(pr)/h)*tc = t
    tc = T_scale(model)
    pc = p_scale(model)
    scale = 0.7
    t7 = scale*tc
    p7 = first(sat_pure(model,t7))
    h = 2.3333333333333335*log(pc/p7)
    #T/Tc = 1(1-log(p/pc)/h)
    return 1/(1-log(p/pc)/h)*tc
end

function sat_pure_p(model,p)
    #=
    model_pred = single_sat_aprox(model)
    t_pred = p_sat_t0(model,p) #prediction temperature
    Vv = V_zeros(_pt,model,p,t_pred)

    V1 = first(Vv)
    #(V1 == V2) && return V1
    V2 = last(Vv)

    =#
    T0 = p_sat_t0(model,p)
    px = p

    _A(V, t) = eos(model, V, t,SA[1.0])
    V1old = zero(p)
    V2old = zero(p)
    Told = zero(p)
    V1,V2 = x0_sat_pure(model,T0)
    Tx = T0
    for i = 1:20
        if i > 1
            V1old = V1
            V2old = V2
        end
        A = z -> _A(z, Tx)
        _V1 = Threads.@spawn volume(model,$p,$Tx;phase=:liquid)
        _V2 = Threads.@spawn volume(model,$p,$Tx;phase=:gas)
        V1 = fetch(_V1)
        V2 = fetch(_V2)

        if abs(V1 - V1old) / V1 < 1e-15 && i > 1
            #println("V1 condition")
            break
        elseif abs(V2 - V2old) / V2 < 1e-15 && i > 1
            #println("V2 condition")
            break
        end

        _px(T) = (_A(V1, T) - _A(V2, T)) / (V2 - V1) - p
        #_px(T) = exp(_A(V1,T)-_A(V2,T)) - one(Tx)
        Told = Tx
        Tx = Solvers.ad_newton(_px, Tx)
        if abs(Told - Tx) < 1e-10 * Tx
            #println("P condition")
            break
        end
    end

    return (Tx,V1,V2)
end

function naive_sat_pure_p(model,p)
    T0 = p_sat_t0(model,p)
    ft0(t) = log(first(sat_pure(model,t))/p)
    tt0 =T0
    res = Solvers.ad_newton(ft0,tt0)
    T = res
    Vl = volume(model,p,T,phase=:l)
    Vv = volume(model,p,T,phase=:v)
    return (T,Vl,Vv)
end

export sat_pure, crit_pure, enthalpy_vap

function spinodals(model,T,Vx = nothing)
    T7 = (0.7)*one(T)*T_scale(model)
    #@show T/T_scale(model)
    #f0(_V) = 1/last(p∂p∂V(model,_V,T))

    lb_v = lb_volume(model)
    f0(k) = last(p∂p∂V(model,lb_v + exp(k),T))

    if Vx === nothing
    if T > T7
        sp_l7 = volume_compress(model,zero(T),T7)
        Tc,pc,Vc = crit_pure(model)
        V00 = log(sp_l7 - lb_v)
        sp_l = Roots.find_zero(f0,V00)

        sp_l = lb_v + exp(sp_l)
    else
        #zero pressure approach, the volume_compress function tends to stop at a spinodal at the liquid side
        sp_l = volume_compress(model,zero(T),T)
    end
    else
        V00 = log(Vx - lb_v)
        sp_l = Roots.find_zero(f0,V00)
        sp_l = lb_v + exp(sp_l)

    end
    #now calculate gas spinodal
    #@show Pmin = pressure(model,sp_l,T)
    #p could be negative, and a volume search fails here.
    #@show volume_virial(model,Pmin,T)
    return sp_l#sp_v)
end
#=
function vsa(model,T)
    k = -0.5*one(T)#*(p_scale(model)/lb_volume(model))
    B = second_virial_coefficient(model,T)
    #=basis
    obtain the value at which ∂p/∂V = -0.5
    0.5 is a randomly chosen value.
    the idea is to choose the lowest value possible such
    as not pass the gas spinodal.
    with the value calculated, a pressure and then an aproximate liquid volume is calculated.

    The main equation is:
    z = 1+B/V = pV/RT
    p = RT/V + RTB/V2
    ∂p/∂V = -RT/V2 -2RTB/V3 = 0.5
    -RTV -2RTB = 0.5V3
    0 = 0.5V3 + RTV + 2RTB
    0.5V3 + 0V2 + RTV + 2RTB = 0
    (0.5,0,RT,2RTB)

    on B = B(V)
    p = RT/V + RTB/V2
    ∂p/∂V = -RT/V2 -2RT/V3*B + ∂B*RT/V2 = 0.5
    ∂p/∂V = -RT(1+∂B)/V2 -2RT/V3*B = 0.5

    -RTV -2RTB = 0.5V3
    0 = 0.5V3 + RTV + 2RTB
    0.5V3 + 0V2 + RTV + 2RTB = 0
    (0.5,0,RT,2RTB)
    =#
    #k*one(T),zero(T),R̄*T,2*B
    RT = R̄*T
    Vv = -2*B

    #Vv = Vv_vsa(B,zero(T),T,k)
    p,dpdV = p∂p∂V(model,Vv,T)
    #now, to calculate a new better aproximate of the function at that volume:
    #p = p0 + dpdV0(V-V0) / multiplying by V/RT and adding 1-1
    #pV/RT = 1-1+p0V/RT + dpdV0(V-V0)*V/RT
    #z = 1 + (-1+p0V/RT + dpdV0(V-V0)*V/RT)*V/V
    #z = 1 + B/V
    Bi(V) = (-1+p*V/RT + dpdV*(V-Vv)*V/RT)*V

    @show dpdV
    @show Bi(Vv)
    @show B
    return Vv
end
=#

#=

=#


text:
using Manifolds, Manopt, Test, ManifoldsBase
using LinearAlgebra: I

using Random
Random.seed!(42)
# Test the additional manifold functions
#
@testset "Additional Manifold functions" begin
    @testset "mid point & reflect" begin
        M = Sphere(2)
        p = [1.0, 0.0, 0.0]
        q = [0.0, 1.0, 0.0]

        r = mid_point(M, p, q)
        r2 = similar(r)
        mid_point!(M, r2, p, q)
        r3 = shortest_geodesic(M, p, q, 0.5)
        r4 = mid_point(M, p, q, q)
        @test isapprox(M, r, r2)
        @test isapprox(M, r2, r3)
        @test isapprox(M, r3, r4)
        r5 = similar(r4)
        mid_point!(M, r5, p, q, q)
        @test isapprox(M, r4, r5)

        r4 = mid_point(M, p, -p, q)
        r5 = similar(r4)
        mid_point!(M, r5, p, -p, q)
        @test isapprox(M, r4, q)
        @test isapprox(M, r4, r5)

        @test isapprox(M, reflect(M, p, q), -q)
        qA = similar(q)
        reflect!(M, qA, p, q)
        @test isapprox(M, qA, -q)
        f = x -> x
        @test reflect(M, f, q) == q
        qA = similar(q)
        reflect!(M, qA, f, q)
        @test qA == q

        M2 = Euclidean(2)
        p2 = [1.0, 0.0]
        q2 = [0.0, 1.0]
        s = mid_point(M2, p2, q2)
        s2 = similar(s)
        mid_point!(M2, s2, p2, q2)
        @test s == s2
        @test s == (p2 + q2) / 2
        s = mid_point(M2, p2, q2, s)
        @test s == s2
        s2 = similar(s)
        mid_point!(M2, s2, p2, q2, s)
        @test s == s2

        M = Circle()
        p = 0
        q = π
        @test mid_point(M, p, q, 1.0) ≈ π / 2
        @test mid_point(M, p, q, -1.0) ≈ -π / 2
        @test mid_point(M, 0, π / 2) ≈ π / 4
    end
    @testset "random" begin
        Mc = Circle()
        pc = random_point(Mc)
        @test is_point(Mc, pc, true)
        Xc = random_tangent(Mc, pc)
        @test is_vector(Mc, pc, Xc, true)

        Me = Euclidean(3)
        pe = random_point(Me)
        @test is_point(Me, pe, true)
        pe2 = random_point(Me, :Gaussian)
        @test is_point(Me, pe2, true)
        pe3 = random_point(Me, :Gaussian, 1.0)
        @test is_point(Me, pe3, true)
        Xe = random_tangent(Me, pe)
        @test is_vector(Me, pe, Xe, true)

        Mg = Grassmann(3, 2)
        pg = random_point(Mg)
        @test is_point(Mg, pg, true)
        Xg = random_tangent(Mg, pg)
        @test is_vector(Mg, pg, Xg, true; atol=1e-14)

        Mp = ProductManifold(Mg, Me)
        pp = random_point(Mp)
        @test is_point(Mp, pp, true)
        Xp = random_tangent(Mp, pp)
        @test is_vector(Mp, pp, Xp, true; atol=1e-14)

        Mp2 = PowerManifold(Me, NestedPowerRepresentation(), 4)
        pp2 = random_point(Mp2)
        @test is_point(Mp2, pp2)
        Xp2 = random_tangent(Mp2, pp2)
        @test is_vector(Mp2, pp2, Xp2)

        Mp3 = PowerManifold(Me, ArrayPowerRepresentation(), 4)
        pp3 = random_point(Mp3)
        @test is_point(Mp3, pp3)
        Xp3 = random_tangent(Mp3, pp3)
        @test is_vector(Mp3, pp3, Xp3)

        Mr = Rotations(3)
        pr = random_point(Mr)
        @test is_point(Mr, pr)
        Xr = random_tangent(Mr, pr)
        @test is_vector(Mr, pr, Xr)

        Mr2 = Rotations(1) # only one point, so there is no randomness, but well
        pr2 = random_point(Mr2)
        @test is_point(Mr2, pr2)
        @test is_vector(Mr2, pr2, random_tangent(Mr2, pr2))

        Mspd = SymmetricPositiveDefinite(3)
        pspd = random_point(Mspd)
        @test is_point(Mspd, pspd; atol=10^(-14))
        Xspd = random_tangent(Mspd, pspd)
        @test is_vector(Mspd, pspd, Xspd, true; atol=10^(-15))
        Xspd2 = random_tangent(Mspd, pspd, Val(:Rician))
        @test is_vector(Mspd, pspd, Xspd2, true; atol=10^(-15))

        Mst = Stiefel(3, 2)
        pst = random_point(Mst)
        @test is_point(Mst, pst, true)
        Xst = random_tangent(Mst, pst)
        @test is_vector(Mst, pst, Xst, true)

        Msp = Sphere(2)
        psp = random_point(Msp)
        @test is_point(Msp, psp, true)
        Xsp = random_tangent(Msp, psp)
        @test is_vector(Msp, psp, Xsp, true; atol=10^(-15))

        Mh = Hyperbolic(2)
        ph = [0.0, 0.0, 1.0]
        Xh = random_tangent(Mh, ph)
        @test is_vector(Mh, ph, Xh, true)

        MT = TangentBundle(Msp)
        pT = random_point(MT)
        @test is_point(MT, pT, true)
        XT = random_tangent(MT, pT)
        @test is_vector(MT, pT, XT, true; atol=10^(-15))

        Mfr = FixedRankMatrices(4, 6, 2)
        pfr = random_point(Mfr)
        @test is_point(Mfr, pfr, true)
        Xfr = random_tangent(Mfr, pfr)
        @test is_vector(Mfr, pfr, Xfr, true; atol=10^(-14))

        Mse1 = Euclidean(3)
        Mse2 = Rotations(3)
        Mse = SpecialEuclidean(3)
        pse = random_point(Mse)
        @test is_point(Mse1, pse.parts[1], true)
        @test is_point(Mse2, pse.parts[2], true)
        Xse = random_tangent(Mse, pse)
        @test is_vector(Mse1, pse.parts[1], Xse.parts[1], true)
        @test is_vector(Mse2, pse.parts[2], Xse.parts[2], true)
        pse2 = random_point(Mse, :Gaussian)
        @test is_point(Mse1, pse2.parts[1], true)
        @test is_point(Mse2, pse2.parts[2], true)
        Xse2 = random_tangent(Mse, pse2, :Gaussian)
        @test is_vector(Mse1, pse2.parts[1], Xse2.parts[1], true)
        @test is_vector(Mse2, pse2.parts[2], Xse2.parts[2], true)
    end
end


text:
#=
Bitpair operations
==================

Copyright (c) 2017 Ben J. Ward & Luis Yanes

This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.
=#

@inline function count_nonzero_bitpairs(x::UInt64)
    out = x & 0x5555555555555555
    out |= (x & 0xAAAAAAAAAAAAAAAA) >> 1
    return count_ones(out)
end

@inline function count_zero_bitpairs(x::UInt64)
    return 32 - count_nonzero_bitpairs(x)
end

@inline function count_one_bitpairs(x::UInt64)
    out = x & 0x5555555555555555
    out &= (x & 0xAAAAAAAAAAAAAAAA) >> 1
    return count_ones(out)
end


text:
"""
    bisect_interval(a::Arf, b::Arf; log_midpoint::Bool = false)

Returns two tuples, `(a, midpoint)` and `(midpoint, b)`, which
corresponds to splitting the interval in half.

If `log_midpoint = true` then the midpoint is computed in logarithmic
scale. If `a, b > 0` this sets the midpoint to `exp((log(a) + log(b))
/ 2)`. If `a, b < 0` the midpoint is `-exp((log(-a) + log(-b)) / 2)`.
If zero is contained in the interval then currently a normal bisection
is performed (i.e. the midpoint is `(a + b) / 2`).

Logarithmic bisection can be useful if the function in consideration
has a logarithmic behaviour and the interval has numbers exponentially
close to zero, for example `[1e-10000, 1e-1]`. In that case normal
bisection could give very slow convergence.

TODO: Think about how to handle a logarithmic bisection when the
interval contains zero.

The value of `midpoint` is aliased in the two tuples and care should
therefore be taken if doing inplace operations on it.
"""
function bisect_interval(a::Arf, b::Arf; log_midpoint::Bool = false)
    if log_midpoint && a > 0
        c = Arb(a)
        d = Arb(b)
        Arblib.log!(c, c)
        Arblib.log!(d, d)
        Arblib.add!(c, c, d)
        Arblib.mul_2exp!(c, c, -1)
        Arblib.exp!(c, c)
        mid = midpoint(c)
    elseif log_midpoint && b < 0
        c = Arb(a)
        d = Arb(b)
        Arblib.neg!(c, c)
        Arblib.neg!(d, d)
        Arblib.log!(c, c)
        Arblib.log!(d, d)
        Arblib.add!(c, c, d)
        Arblib.mul_2exp!(c, c, -1)
        Arblib.exp!(c, c)
        Arblib.neg!(c, c)
        mid = midpoint(c)
    else
        mid = a + b
        Arblib.mul_2exp!(mid, mid, -1)
    end

    return (a, mid), (mid, b)
end

"""
    bisect_interval(a::Arb, b::Arb; log_midpoint::Bool = false)

Returns two tuples, `(a, midpoint)` and `(midpoint, b)`, which
corresponds to splitting the interval in half.

The method currently doesn't support `log_midpoint = true`, it is
there to make the signature the same as for the `Arf` version.

The value of `midpoint` is aliased in the two tuples and care should
therefore be taken if doing inplace operations on it.
"""
function bisect_interval(a::Arb, b::Arb; log_midpoint::Bool = false)
    log_midpoint && throw(ArgumentError("log_midpoint currently not supported for Arb"))

    mid = a + b
    Arblib.mul_2exp!(mid, mid, -1)

    return (a, mid), (mid, b)
end

"""
    bisect_interval_recursive(a::T, b::T, depth::Integer; log_midpoint::Bool = false) where {T<:Union{Arf,Arb}}

Recursively bisect the interval `[a, b]` a number of `depth` times.
The bisection is done using [`bisect_interval`](@ref) and it returns a
vector with `2^depth` intervals.

It accepts the same arguments as [`bisect_interval`](@ref).
"""
function bisect_interval_recursive(
    a::T,
    b::T,
    depth::Integer;
    log_midpoint::Bool = false,
) where {T<:Union{Arf,Arb}}
    depth >= 0 || Throw(ArgumentError("depth needs to be non-negative, got $depth"))

    res = Vector{NTuple{2,T}}(undef, 2^depth)
    res[1] = (a, b)
    @inbounds for i = 1:depth
        for j in reverse(1:2^(i-1))
            res[2j-1], res[2j] = bisect_interval(res[j]...; log_midpoint)
        end
    end

    return res
end

"""
    check_tolerance(x::Arb; atol = nothing, rtol = nothing)

Return `true` if `x` satisfies the absolute tolerance `atol` or the
relative tolerance `rtol`.

The absolute tolerance is satisfied if the diameter of `x` is less
than or equal to `atol`.

The relative tolerance is satisfied if the diameter of `x` divided by
the absolute value of `x` is less than or equal to `rtol`. If `x`
contains zero this is only satisfied if `x` is exactly zero.

Both `atol` and `rtol` can be given as `nothing`, in which case that
check is skipped. If both of them are `nothing` then it returns `true`
if `x` is finite.
"""
function check_tolerance(x::Arb; atol = nothing, rtol = nothing)
    isfinite(x) || return false
    isnothing(atol) && isnothing(rtol) && return true

    error = radius(Arb, x)
    Arblib.mul_2exp!(error, error, 1)

    !isnothing(atol) && error <= atol && return true

    if !isnothing(rtol)
        Arblib.contains_zero(x) && return iszero(error)

        return error <= rtol * abs(x)
    else
        return false
    end
end

"""
    check_interval(a::Arf, b::Arf)

Throw an error if `a` and `b` does not define a finite interval `[a,
b]`. It checks that `a` and `b` are both finite and satisfy `a <= b`.
"""
function check_interval(a::Arf, b::Arf)
    isfinite(a) && isfinite(b) ||
        throw(ArgumentError("a and b must be finite, got a = $a and b = $b"))
    a <= b || throw(ArgumentError("must have a <= b, got a = $a and b = $b"))
    return nothing
end

"""
    check_interval(Bool, a::Arf, b::Arf)

Return `true` if `a` and `b` defines a finite interval `[a, b]`. It
checks that `a` and `b` are both finite and satisfy `a <= b`.
"""
check_interval(::Type{Bool}, a::Arf, b::Arf) = isfinite(a) && isfinite(b) && a <= b

"""
    format_interval(a, b; digits = 5)

Return a nicely formatted string representing the interval `[a, b]`.

It prints the interval either on the form `[a, b]` or on the form `[c
+/- d]` depending on the width of the width of the interval.

If `a` or `b` is non-finite it always prints an interval of the form
`[a, b]`. Otherwise it checks if the relative error is less than
`10^-digits`, if it is then it prints it as a ball, otherwise as an
interval. When printed in interval form it prints at most `digits`
digits.
"""
function format_interval(a::Arf, b::Arf; digits = 5)
    if isfinite(a) && isfinite(b)
        ball = Arb((a, b))
        if check_tolerance(ball, rtol = Arb(10)^(-digits))
            return string(ball)
        end
    end

    return "[" *
           Base.MPFR._string(convert(BigFloat, a), digits) *
           ", " *
           Base.MPFR._string(convert(BigFloat, b), digits) *
           "]"
end

"""
    taylor_remainder(p::ArbSeries, x::Arb)

Compute `p[degree] * (x - midpoint(x))^degree`, where `degree` is the
degree of `p`

This corresponds to the Lagrange form of the remainder term in a
Taylor expansion around the point `midpoint(x)` of degree `degree -
1` valid on the full interval `x`.

For example, for a function `f`, a ball `x`, any `degree >= 0` and
```
p = f(ArbSeries((x, 1), degree = degree + 1))
q = f(ArbSeries((midpoint(x), 1), degree = degree))
```
we have for any `y ∈ x` that
```
f(y) ∈ q(y) + taylor_remainder(p, x)
```
"""
function taylor_remainder(p::ArbSeries, x::Arb)
    restterm = zero(x)
    Arblib.set!(Arblib.radref(restterm), Arblib.radref(x))
    Arblib.pow!(restterm, restterm, unsigned(Arblib.degree(p)))
    Arblib.mul!(restterm, restterm, Arblib.ref(p, Arblib.degree(p)))
    return restterm
end


text:
"""
Adapted from 2D Types:
Geometric predicates for 3D Vectors, cartesian points (px,py,pz), homogeneuos points (pw,px,py,pz), and planes
	Julia can inline these.
"""

Vector3 = Tuple{Real,Real,Real}
Point3C = Tuple{Real,Real,Real}
Plane3 = Tuple{Real,Real,Real,Real}
Point3H = Tuple{Real,Real,Real,Real}
Point3 = Union{Point3C,Point3H}
Point3L = Tuple{Real,Real,Real,Real,Real} # last cord is lift
Sphere3 = Tuple{Real,Real,Real,Real,Real} # x,y,z,q

toPoint3H((x,y,z)::Point3C)::Point3H = (one(x),x,y,z)
lift(x::Float32,y::Float32,z::Float32)::Point3L = (x,y,z,muladd(widen(z),widen(z), muladd(widen(y),widen(y), widemul(x,x))) ) # promote to double
lift((x,y,z)::Point3C)::Point3L = (x,y,z, muladd(z,z, muladd(y,y, x*x)) ) # promote to double
lift(a::Array{<:Real,1}) = (a[1], a[2], a[3], sum(a.*widen.(a)))
dot(u,v) = sum(u.*v)
# cartesian and homogeneous orientation determinants

# cartesian plane equation
function plane(p::Point3C, q::Point3C, r::Point3C)::Plane3
	pq = q.-p
	pr = r.-p
	return (det3(p,pq,pr,1,2,3), -det2(pq,pr,2,3), det2(pq,pr,1,3), -det2(pq,pr,1,2))
end

# computing homogeneous plane equations
plane(p::Point3H,q::Point3H,r::Point3H)::Plane3 = (det3(p,q,r,2,3,4),
													-det3(p,q,r,1,3,4),
													det3(p,q,r,1,2,4),
													-det3(p,q,r,1,2,3))

orientptpl((x,y,z)::Point3C,(pw,px,py,pz)::Plane3) = orientptpl(toPoint3H(x,y,z),(pw,px,py,pz))
orientptpl((w,x,y,z)::Point3H,(pw,px,py,pz)::Plane3) = muladd(px,x,muladd(py,y,muladd(pz,z,pw)))
abovePl(q::Point3,p::Plane3)::Bool = orientptpl(q,p)>zero(p[1])
belowPl(q::Point3,p::Plane3)::Bool = orientptpl(q,p)<zero(p[1])

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
module ExactD

using SparseArrays, LinearAlgebra

include("./states.jl")
include("./many_body_operators.jl")
include("./entropy.jl")
include("./auxiliary.jl")
include("./expected.jl")
include("./exponential_ops.jl")

end  # module ExactD


text:
module GLMMetrics

using Statistics, LinearAlgebra

using ..RelayUtils

export PerformanceMetric, ROCArea, JSDivergence, PRArea, BinomialLikelihood,
    RelativeBinomialLikelihood, RRI

export eval_and_store!, key_name, better_than, worst_value, findbest

abstract type PerformanceMetric end

@inline better_than(::Type{<:PerformanceMetric}, x::Real, y::Real) = x > y
@inline worst_value(::Type{<:PerformanceMetric}) = -Inf
@inline nanmean(x::AbstractVector{<:Real}) = mean(filter(!isnan, x))
# ============================================================================ #
function findbest(p::T) where T<:PerformanceMetric
    best = worst_value(T)
    kb = 1
    for k in eachindex(p.x)
        if better_than(T, p.x[k], best)
            best = p.x[k]
            kb = k
        end
    end
    return best, kb
end
# ============================================================================ #
struct ROCArea <: PerformanceMetric
    x::Vector{Float64}
end

ROCArea(nfold::Integer) = ROCArea(fill(NaN, nfold))

function eval_and_store!(r::ROCArea, y::AbstractVector{<:Real}, yp::AbstractVector{<:Real}, k::Integer)
    r.x[k], _ = RelayUtils.roca(y, yp)
    return r
end

Statistics.mean(r::ROCArea) = nanmean(r.x)
function Base.copy!(dst::ROCArea, src::ROCArea)
    dst.x .= src.x
    return dst
end
key_name(::Type{ROCArea}) = "roca"

# ============================================================================ #
struct PRArea <: PerformanceMetric
    x::Vector{Float64}
end

PRArea(nfold::Integer) = PRArea(fill(NaN, nfold))

function eval_and_store!(r::PRArea, y::AbstractVector{<:Real}, yp::AbstractVector{<:Real}, k::Integer)
    r.x[k], _ = RelayUtils.precision_recall(y, yp)
    return r
end

Statistics.mean(r::PRArea) = nanmean(r.x)
function Base.copy!(dst::PRArea, src::PRArea)
    dst.x .= src.x
    return dst
end
key_name(::Type{PRArea}) = "pra"

# ============================================================================ #
struct JSDivergence <: PerformanceMetric
    x::Vector{Float64}
    p::Vector{Float64}
end

JSDivergence(nfold::Integer) = JSDivergence(fill(NaN, nfold), fill(NaN, nfold))

function eval_and_store!(j::JSDivergence, y::AbstractVector{<:Real}, yp::AbstractVector{<:Real}, k::Integer)
    #                                                        niter, nbin
    j.x[k], j.p[k] = RelayUtils.permute_jensen_shannon(y, yp, 1000, 100)
    return j
end

Statistics.mean(j::JSDivergence) = nanmean(j.x), nanmean(j.p)
function Base.copy!(dst::JSDivergence, src::JSDivergence)
    dst.x .= src.x
    dst.p .= src.p
    return dst
end
key_name(::Type{JSDivergence}) = "jsd"
# ============================================================================ #
struct BinomialLikelihood <: PerformanceMetric
    x::Vector{Float64}
end
BinomialLikelihood(nfold::Integer) = BinomialLikelihood(fill(NaN, nfold))

# NOTE: smaller is better for negative log likelihood of Binomial PDF
@inline better_than(::Type{BinomialLikelihood}, x::Real, y::Real) = x < y
@inline worst_value(::Type{BinomialLikelihood}) = +Inf

function eval_and_store!(r::BinomialLikelihood, y::AbstractVector{<:Real}, yp::AbstractVector{<:Real}, k::Integer)
    # r.x[k] = -(dot(log.(yp), y) + dot(1.0 .- y, log.(1.0 .- yp)))
    r.x[k] = -RelayUtils.binomial_lli_turbo(y, yp) / length(y)
    return r
end

Statistics.mean(r::BinomialLikelihood) = nanmean(r.x)
function Base.copy!(dst::BinomialLikelihood, src::BinomialLikelihood)
    dst.x .= src.x
    return dst
end
key_name(::Type{BinomialLikelihood}) = "bili"
# ============================================================================ #
struct RelativeBinomialLikelihood <: PerformanceMetric
    x::Vector{Float64}
end
RelativeBinomialLikelihood(nfold::Integer) = RelativeBinomialLikelihood(fill(NaN, nfold))

# for relative likelihood, larger values are better
function eval_and_store!(r::RelativeBinomialLikelihood, y::Vector{Bool}, yp::Vector{Float64}, k::Integer)
    r.x[k] = (RelayUtils.binomial_lli_turbo(yp, y) - RelayUtils.binomial_lli(y)) / length(y) #null
    return r
end

Statistics.mean(r::RelativeBinomialLikelihood) = nanmean(r.x)
function Base.copy!(dst::RelativeBinomialLikelihood, src::RelativeBinomialLikelihood)
    dst.x .= src.x
    return dst
end
key_name(::Type{RelativeBinomialLikelihood}) = "rbili"
# ============================================================================ #
# Bernoulli variant of cross-validated (CV) single-spike information, which in
# this context becomes CV single-spike Bernoulli information *OR* the
# single-spike relative relay information (RRI)
# c.f. Williamson et al. (2015) PLOS Comp Bio, eq. 17, 20, and 59
struct RRI <: PerformanceMetric
    x::Vector{Float64}
    pred::Vector{Float64}
end
RRI(nfold::Integer, nret::Integer=0) = RRI(fill(NaN, nfold), fill(NaN, nret))

# for RRI, larger values are better
function eval_and_store!(r::RRI, y::Vector{Bool}, yp::Vector{Float64}, k::Integer, kspk::Vector{<:Integer}=Int[])
    # divide by log(2) so our units are bits-per-event
    r.x[k] = (RelayUtils.binomial_lli_turbo(yp, y) - RelayUtils.binomial_lli(y)) / length(y) / log(2)
    if !isempty(kspk)
        r.pred[kspk] .= yp
    end
    return r
end

Statistics.mean(r::RRI) = nanmean(r.x)
function Base.copy!(dst::RRI, src::RRI)
    dst.x .= src.x
    return dst
end
key_name(::Type{RRI}) = "rri"
# ============================================================================ #
end


text:
using Plots

f(t) = exp(im*t)
t = 0:0.01:2pi
z = f.(t)

x = real.(z)
y = imag.(z)

plot(x, y)
savefig("eulers_formula.png")


text:
struct SAGA_basic_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg}
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                  # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    γ::Maybe{R}             # stepsize
    SAG::Bool               # to activate SAG
end

mutable struct SAGA_basic_state{R<:Real,Tx}
    s::Array{Tx}            # table of x_j- γ_j/N nabla f_j(x_j)
    γ::R                    # stepsize
    av::Tx                  # the running average
    z::Tx
    # some extra placeholders
    ind::Int         # running idx set
    ∇f_temp::Tx             # placeholder for gradients
    w::Tx
end

function SAGA_basic_state(s, γ::R, av::Tx, z::Tx) where {R,Tx}
    return SAGA_basic_state{R,Tx}(s, γ, av, z, 1, copy(av), copy(av))
end

function Base.iterate(iter::SAGA_basic_iterable{R,C,Tx}) where {R,C,Tx}
    N = iter.N
    # updating the stepsize
    if iter.γ === nothing
        if iter.L === nothing
            @warn "smoothness parameter absent"
            return nothing
        else
            L_M = maximum(iter.L)
            γ = iter.SAG ? 1 / (16 * L_M) : γ = 1 / (3 * L_M)
        end
    else
        γ = iter.γ # provided γ
    end
    # computing the gradients and updating the table
    s = Vector{Tx}(undef, 0)
    for i = 1:N
        ∇f, ~ = gradient(iter.F[i], iter.x0)
        push!(s, ∇f) # table of x_i
    end
    #initializing the vectors
    av = sum(s) ./ N # the running average
    z, ~ = prox(iter.g, (1 - γ) .* iter.x0, γ)
    state = SAGA_basic_state(s, γ, av, z)
    return state, state
end

function Base.iterate(iter::SAGA_basic_iterable{R}, state::SAGA_basic_state{R}) where {R}

    state.ind = rand(1:iter.N)
    gradient!(state.∇f_temp, iter.F[state.ind], state.z)
    if iter.SAG
        @. state.av += (state.∇f_temp - state.s[state.ind]) / iter.N
        @. state.w = state.z - state.γ * state.av
    else
        @. state.w = state.z - state.γ * (state.∇f_temp - state.s[state.ind] + state.av)
        @. state.av += (state.∇f_temp - state.s[state.ind]) / iter.N
    end
    prox!(state.z, iter.g, state.w, state.γ)
    state.s[state.ind] .= state.∇f_temp  #update x_i

    return state, state
end

solution(state::SAGA_basic_state) = state.z

#TODO: minibatch

max_issues_repo_name: kul-forbes/CIAOAlgorithms

max_forks_repo_name: kul-forbes/CIAOAlgorithms


text:
module Stepwise

using Statistics, LinearAlgebra
using Parameters
using Base: ==

export AbstractDirection, Bidirection, Unidirection, ForwardBackward, BackwardForward, Forward, Backward,
    AbstractCriterion, AIC, BIC,
    AbstractAlgorithm, StepwiseFeatureSelection,
    AbstractAlgResult, SFSResult,
    Model, Null, Full, null_or_full, ==, update!,
    stepwise

abstract type AbstractDirection end
abstract type Bidirection <: AbstractDirection end
abstract type Unidirection <: AbstractDirection end
struct ForwardBackward <: Bidirection end
struct BackwardForward <: Bidirection end
struct Forward <: Unidirection end
struct Backward <: Unidirection end

abstract type AbstractCriterion end
struct AIC <: AbstractCriterion end
struct BIC <: AbstractCriterion end
(::AIC)(𝐗, 𝐲, 𝚺⁻¹) = 𝐲' * (𝚺⁻¹ - 𝚺⁻¹ * 𝐗 * pinv(𝐗' * 𝚺⁻¹ * 𝐗) * 𝐗' * 𝚺⁻¹) * 𝐲 + 2 * size(𝐗, 2)
(::BIC)(𝐗, 𝐲, 𝚺⁻¹) = 𝐲' * (𝚺⁻¹ - 𝚺⁻¹ * 𝐗 * pinv(𝐗' * 𝚺⁻¹ * 𝐗) * 𝐗' * 𝚺⁻¹) * 𝐲 + log(length(𝐲)) * size(𝐗, 2)

abstract type AbstractAlgorithm end
@with_kw struct StepwiseFeatureSelection{S <: AbstractDirection,T <: AbstractCriterion,U <: Integer} <: AbstractAlgorithm
    direction::S = ForwardBackward()
    criterion::T = AIC()
    max_step::U = typemax(Int)
end

# define "model"
struct Null end; struct Full end
null_or_full(::Forward) = Null(); null_or_full(::ForwardBackward) = Null()
null_or_full(::Backward) = Full(); null_or_full(::BackwardForward) = Full()
struct Model{T <: Integer}
    M::Vector{T} # the set of features in this model
    Mᶜ::Vector{T} # the set of features not included in this model
    Model{T}(::Null, dims::T) where {T} = new{T}(T[], collect(T, 1:dims))
    Model{T}(::Full, dims::T) where {T} = new{T}(collect(T, 1:dims), T[])
end
Base.:(==)(m::Model, m′::Model) = issetequal(m.M, m′.M) && issetequal(m.Mᶜ, m′.Mᶜ)
Model(direction::AbstractDirection, dims::T) where {T} = Model{T}(null_or_full(direction), dims)
## update model
function update!(model::Model, idx)
    @unpack M, Mᶜ = model
    if idx > 0
        @assert idx ∈ Mᶜ
        push!(M, idx); setdiff!(Mᶜ, idx)
    else
        @assert -idx ∈ M
        setdiff!(M, -idx); push!(Mᶜ, -idx)
    end
end

abstract type AbstractAlgResult end
struct SFSResult{S,T} <: AbstractAlgResult
    model::Model{S}
    history::T
end
function Base.show(io::IO, ::MIME"text/plain", result::SFSResult{S,T}) where {S,T}
    println(io, "SFSResult{$S,$T}")
    println(io, " selected features: ", result.model.M)
    println(io, " history of this algorithm: ", result.history)
end

# stepwise feature selection
function stepwise(𝐗, 𝐲, 𝚺, direction::Bidirection, criterion::AbstractCriterion; max_step=typemax(Int))
    𝚺⁻¹ = I / 𝚺
    d = size(𝐗, 2)
    model = Model(direction, d)
    history = typeof(d)[]
    L_prev = Inf

    for i = 1:max_step
        @unpack M, Mᶜ = model
        minL = Inf
        idx = 0
        # delete
        for j ∈ M
            L = criterion(𝐗[:, setdiff(M, j)], 𝐲, 𝚺⁻¹)
            L < minL && (minL = L; idx = -j)
        end
        # add
        for j′ ∈ Mᶜ
            L = criterion(𝐗[:, union(M, j′)], 𝐲, 𝚺⁻¹)
            L < minL && (minL = L; idx = +j′)
        end

        # if the criterion is not improved, then finish this algorithm
        L_prev ≤ minL && (return SFSResult(model, history))

        # update
        L_prev = minL
        push!(history, idx)
        update!(model, idx)
    end
    return SFSResult(model, history)
end
stepwise(𝐗, 𝐲_obs, 𝚺) = stepwise(𝐗, 𝐲_obs, 𝚺, ForwardBackward(), AIC()) # default
function stepwise(𝐗, 𝐲, 𝚺, direction::Forward, criterion::AbstractCriterion; max_step=typemax(Int))
    𝚺⁻¹ = I / 𝚺
    d = size(𝐗, 2)
    model = Model(direction, d)
    history = typeof(d)[]
    L_prev = Inf

    for i = 1:max_step
        @unpack M, Mᶜ = model
        minL = Inf
        idx = 0
        # add
        for j′ ∈ Mᶜ
            L = criterion(𝐗[:, union(M, j′)], 𝐲, 𝚺⁻¹)
            L < minL && (minL = L; idx = +j′)
        end

        # if the criterion is not improved, then finish this algorithm
        L_prev ≤ minL && (return SFSResult(model, history))

        # update
        L_prev = minL
        push!(history, idx)
        update!(model, idx)
    end
    return SFSResult(model, history)
end
function stepwise(𝐗, 𝐲, 𝚺, direction::Backward, criterion::AbstractCriterion; max_step=typemax(Int))
    𝚺⁻¹ = I / 𝚺
    d = size(𝐗, 2)
    model = Model(direction, d)
    history = typeof(d)[]
    L_prev = Inf

    for i = 1:max_step
        M = model.M
        minL = Inf
        idx = 0
        # delete
        for j ∈ M
            L = criterion(𝐗[:, setdiff(M, j)], 𝐲, 𝚺⁻¹)
            L < minL && (minL = L; idx = -j)
        end

        # if the criterion is not improved, then finish this algorithm
        L_prev ≤ minL && (return SFSResult(model, history))

        # update
        L_prev = minL
        push!(history, idx)
        update!(model, idx)
    end
    return SFSResult(model, history)
end
function stepwise(𝐗, 𝐲, 𝚺, alg::StepwiseFeatureSelection)
    @unpack direction, criterion, max_step = alg
    stepwise(𝐗, 𝐲, 𝚺, direction, criterion; max_step=max_step)
end # wrapper of stepwise

end # module Stepwise


text:
export Synset, word_count, words

struct Synset
    offset::Int
    lex_filenum::Int
    word_counts::Dict{String, Int}
    synset_type::Char
    pos::Char
    pointers::Vector{Pointer}
    gloss::String
end

const ∅ = Synset(
    -1,
    -1,
    Dict{String, Int}(),
    '-',
    '-',
    Vector{Pointer}(),
    ""
)

is_nothing(synset::Synset) = synset.offset == -1

parse_int_hex(s) = parse(Int, string("0x", s))

function Synset(raw_line::AbstractString, pos::Char)
    dl_parts = split(strip(raw_line), " | ")
    line = split(dl_parts[1], SPACE)
    gloss = join(dl_parts[2:end], " | ")

    offset = parse(Int, popfirst!(line))
    lex_filenum = parse(Int, popfirst!(line))
    synset_type = popfirst!(line)[1]

    n_words = parse_int_hex(popfirst!(line))  # hex!
    word_counts = Dict{String, Int}()
    for _ in 1:n_words
        k = popfirst!(line)
        word_counts[k] = parse_int_hex(popfirst!(line))
    end

    n_pointers = parse(Int, popfirst!(line))
    pointers = [
        Pointer(
            string(popfirst!(line)[1]),
            parse(UInt, popfirst!(line)),
            popfirst!(line)[1],
            popfirst!(line)
        )
        for _ in 1:n_pointers
    ]

    Synset(
        offset,
        lex_filenum,
        word_counts,
        synset_type,
        pos,
        pointers,
        gloss
    )
end

word_count(synset::Synset) = length(synset.word_counts)
words(synset::Synset) = keys(synset.word_counts)

function Base.show(io::IO, synset::Synset)
    ws = join([replace(word, "_" => " ") for word in words(synset)], ", ")
    print(io, "($(synset.synset_type)) $ws ($(synset.gloss))")
end


text:
using ThresholdStability
using Test

include("indicatort.jl")
include("automatont.jl")
include("cksvart.jl")
include("veroneset.jl")
include("stabilityt.jl")


text:
using DifferentialEquations
using Plots
using LinearAlgebra
using MatrixEquations

const M = 1
const m = 0.7
const l = 0.15
const D_x = 0.01
const D_θ = 0.01
const g = 9.80665

function fxu(x::Vector{T}, u::T) where T
    #dx = Vector{T}(undef, 4)
    dx1 = x[2]
    dx2 = (3*D_θ*cos(x[3])*x[4]/(4*M*l + 3*l*m*sin(x[3])^2 + l*m)
    - 4*D_x*l*x[2]/(4*M*l + 3*l*m*sin(x[3])^2 + l*m)
    - 3*g*l*m*sin(2*x[3])/(8*M*l + 6*l*m*sin(x[3])^2 + 2*l*m)
    + 4*l^2*m*sin(x[3])*x[4]^2/(4*M*l + 3*l*m*sin(x[3])^2 + l*m)
    + 4*l*u/(4*M*l + 3*l*m*sin(x[3])^2 + l*m))
    println("dx2 = ", dx2)
    dx3 = x[4]
    dx4 = (-3*D_θ*M*x[4]/(4*M*l^2*m + 3*l^2*m^2*sin(x[3])^2 + l^2*m^2)
    - 3*D_θ*m*x[4]/(4*M*l^2*m + 3*l^2*m^2*sin(x[3])^2 + l^2*m^2)
    + 3*D_x*l*m*cos(x[3])*x[2]/(4*M*l^2*m + 3*l^2*m^2*sin(x[3])^2 + l^2*m^2)
    + 3*M*g*l*m*sin(x[3])/(4*M*l^2*m + 3*l^2*m^2*sin(x[3])^2 + l^2*m^2)
    + 3*g*l*m^2*sin(x[3])/(4*M*l^2*m + 3*l^2*m^2*sin(x[3])^2 + l^2*m^2)
    - 3*l^2*m^2*sin(2*x[3])*x[4]^2/(8*M*l^2*m + 6*l^2*m^2*sin(x[3])^2 + 2*l^2*m^2)
    - 3*l*m*u*cos(x[3])/(4*M*l^2*m + 3*l^2*m^2*sin(x[3])^2 + l^2*m^2))
    dx = [dx1, dx2, dx3, dx4]
    println("x = ", x)
    println("dx = ", dx)
    return dx
end

function A(x::Vector{T}) where T
    out = Matrix{T}(undef, 4, 4)
    out[1,1] = 0
    out[1,2] = 1
    out[1,3] = 0
    out[1,4] = 0
    out[2,1] = 0
    out[2,2] = -4*D_x/(4*M + 3*m*sin(x[3])^2 + m)
    out[2,3] = (-6*m*(3*D_θ*cos(x[3])*x[4] - 4*D_x*l*x[2] - 3*g*l*m*sin(2*x[3])/2
    + 4*l^2*m*sin(x[3])*x[4]^2)*sin(x[3])*cos(x[3])/(l*(4*M
    + 3*m*sin(x[3])^2 + m)^2) + (-3*D_θ*sin(x[3])*x[4] - 3*g*l*m*cos(2*x[3])
    + 4*l^2*m*cos(x[3])*x[4]^2)/(l*(4*M + 3*m*sin(x[3])^2 + m)))
    out[2,4] = (3*D_θ*cos(x[3]) + 8*l^2*m*sin(x[3])*x[4])/(l*(4*M + 3*m*sin(x[3])^2 + m))
    out[3,1] = 0
    out[3,2] = 0
    out[3,3] = 0
    out[3,4] = 1
    out[4,1] = 0
    out[4,2] = 3*D_x*cos(x[3])/(l*(4*M + 3*m*sin(x[3])^2 + m))
    out[4,3] = (-18*(-D_θ*M*x[4] - D_θ*m*x[4] + D_x*l*m*cos(x[3])*x[2] + M*g*l*m*sin(x[3])
    + g*l*m^2*sin(x[3]) - l^2*m^2*sin(2*x[3])*x[4]^2/2)*sin(x[3])*cos(x[3])/(l^2*(4*M + 3*m*sin(x[3])^2 + m)^2)
    + 3*(-D_x*l*m*sin(x[3])*x[2] + M*g*l*m*cos(x[3]) + g*l*m^2*cos(x[3])
    - l^2*m^2*cos(2*x[3])*x[4]^2)/(l^2*m*(4*M + 3*m*sin(x[3])^2 + m)))
    out[4,4] = 3*(-D_θ*M - D_θ*m - l^2*m^2*sin(2*x[3])*x[4])/(l^2*m*(4*M + 3*m*sin(x[3])^2 + m))
    out
end

function B(x::Vector{T}) where T
    out = Vector{T}(undef, 4)
    out[1] = 0
    out[2] = 4*l/(4*M*l + 3*l*m*(sin(x[3])^2) + l*m)
    out[3] = 0
    out[4] = -3*l*m*cos(x[3])/(4*M*l^2*m + 3*l^2*m^2*(sin(x[3])^2) + l^2*m^2)
    out
end

struct PID{T}
    Kp::T
    Ki::T
    Kd::T
end

function input(p::PID{T}, x::Vector{T}) where T
    0
end

function dx!(dx::Vector{T}, x::Vector{T}, p::PID{T}, t::T) where T
    fxu!(dx, x, input(p, x))
end

function run_simuration(p::PID{T}, t_span::Tuple{T, T}, Δt::T, x₀::Vector{T}) where T
    prob = ODEProblem(dx!, x₀, t_span, p)
    sol = solve(prob)
    sol
end

mutable struct LQR{T}
    Q::Matrix{T}
    R::Matrix{T}
    F::Matrix{T}
end

function init!(p::LQR{T}, A::Matrix{T}, B::Vector{T}) where T
    S = zero(B)
    P, _, _ = arec(A, B, p.R, p.Q, S)
    p.F = inv(p.R) * B' * P
end

function input(p::LQR{T}, x::Vector{T}) where T
    (p.F * x)[1,1]
end

function dx(x::Vector{T}, p::LQR{T}, t::T) where T
    println("t = ", t)
    u = input(p, x)
    println("u = ", u)
    fxu(x, u)
end

function run_simuration(p::LQR{T}, t_span::Tuple{T, T}, Δt::T, x₀::Vector{T}) where T
    init!(p, A(zeros(T, 4)), B(zeros(T, 4)))
    prob = ODEProblem(dx, x₀, t_span, p)
    sol = solve(prob)
    sol
end

x₀ = [0, 0, π/6, 0]
p = LQR(diagm([10.0, 0, 100, 0]), Matrix{Float64}(I, 1, 1).*0.01, Matrix{Float64}(undef, 4, 4))
sol = run_simuration(p, (0.0, 5.0), 0.01, x₀)
plot(sol)


text:
#Adapted from IndirectArray

struct Vec{T, N} <: AbstractArray{T,1}
    cpp_object
    data::AbstractArray{T, 1}
    cpp_allocated::Bool
    @inline function Vec{T, N}(obj) where {T, N}

        new{T, N}(obj, Base.unsafe_wrap(Array{T, 1}, Ptr{T}(obj.cpp_object), N), true)
    end

    @inline function Vec{T, N}(data_raw::AbstractArray{T, 1}) where {T, N}
        if size(data_raw, 1) != N
            throw("Array is improper Size for Vec declared")
        end
        new{T, N}(nothing, data_raw, false)
    end
end

function Base.deepcopy_internal(x::Vec{T,N}, y::IdDict) where {T, N}
    if haskey(y, x)
        return y[x]
    end
    ret = Base.copy(x)
    y[x] = ret
    return ret
end

Base.size(A::Vec) = Base.size(A.data)
Base.axes(A::Vec) = Base.axes(A.data)
Base.IndexStyle(::Type{Vec{T,N}}) where {T, N} = IndexLinear()

Base.strides(A::Vec{T,N}) where {T, N} = (1)
function Base.copy(A::Vec{T,N}) where {T, N}
    return Vec{T, N}(copy(A.data))
end
Base.pointer(A::Vec) = Base.pointer(A.data)

Base.unsafe_convert(::Type{Ptr{T}}, A::Vec{S, N}) where {T, S, N} = Base.unsafe_convert(Ptr{T}, A.data)

@inline function Base.getindex(A::Vec{T,N}, I::Int) where {T, N}
    @boundscheck checkbounds(A.data, I)
    return A.data[I]
end

@inline function Base.setindex!(A::Vec, x, I::Int)
    @boundscheck checkbounds(A.data, I)
    A.data[I] = x
    return A
end


text:
using QuNet

# println("-----Running maxpaths.jl-----")
# include("maxpaths.jl")
# println("-----Running percolations.jl-----")
# include("percolations.jl")
# println("-----Running gridsize.jl-----")
# include("gridsize.jl")
# println("-----Running userpairs.jl-----")
# include("userpairs.jl")
# println("-----Running timedepth.jl-----")
# include("timedepth.jl")
println("Running bandwidth_userpairs.jl")
include("new_bandwidth_userpairs.jl")
println("-----Running bandwidth_memory.jl")
include("new_bandwidth_memory.jl")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
"""
Streamable Random Number Generator, see `AbstractStreamableRNG`
"""

module RDST

using Random

import Base: rand, show, copy

export AbstractRNGStream, AbstractStreamableRNG
export copy
export checkseed, MRG32k3a, rand, srand, reset_stream!, reset_substream!, next_substream!, MRG32k3aGen, show, get_state, next_stream
export Xoshiro256p, Xoshiro256plusGen
export Rigged

abstract type AbstractRNGStream end
abstract type AbstractStreamableRNG <: AbstractRNG end # object of subtype of AbstractStreamableRNG are StreamableRNG
                                                       #currently just mrg32k3a and xoshiro256plus
###mrg32k3a
include("mrg32k3a/mrg32k3a.jl")
include("mrg32k3a/mrg32k3a_types.jl")
include("mrg32k3a/stolen.jl")
#xoshiro256+

include("xoshiro/xoshiro256plus.jl")
include("xoshiro/xoshiro256plus_types.jl")

#Rigged

include("Rigged/Rigged.jl")
end # module


text:
using PlusMinus
using Test

@test ±(1,2) == (1-2,1+2)

@test (±).(collect(1:3),1) == [(i-1,i+1) for i in 1:3]

@test (±).(ones(3),collect(1:3)) == [(1-i,1+i) for i in 1:3]

@test (±).(ones(3,2),collect(1:3)) == [(1-i,1+i) for i in 1:3, y in 1:2]

@test plusminus(1.5,2) == (1.5-2,1.5+2)

@test plusminus.(collect(1:3),1) == [(i-1,i+1) for i in 1:3]


text:
module ReviseLoop

using Dates

export edit_pos_num, more_than, looper

@noinline function more_than(x)
    pos_num = 1
    x + pos_num
end

function looper(f)
    while true
        sleep(5)
        i = rand(1:9)
        println((Dates.now(), i => f(i)))
    end
end

function edit_pos_num(n)
    open(@__FILE__, "r+") do io
        s = read(io, String)
        ns = replace(s, r"(\s+pos_num =) \d+"m => SubstitutionString("\\1 $n"))
        seekstart(io)
        write(io, ns)
    end
end

end # module


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using DifferentialEquations
using LinearAlgebra
using Plots

f(u,p,t) = 0.98*u

u0 = 1.0
tspan = (0.0, 1.0)

prob = ODEProblem(f, u0, tspan)

sol = solve(prob)

plot(sol; label="")

sol(0.8)

function lorenz(u, p, t)
    σ, ρ, β = p
    x_t = σ*(u[2]-u[1])
    y_t = u[1]*(ρ-u[3]) - u[2]
    z_t = u[1]*u[2] - β*u[3]
    return [x_t; y_t; z_t]
end

u0 = [1.0; 0.0; 0.0]
p = [10; 28; 8/3]

tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan, p)

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
abstract type AbstractRecorder end

################################################################################
### ConfigRecorder
################################################################################

"""
    ConfigRecorder(mc, model[; rate = 10])

Creates a `ConfigRecorder` object which accumulates configurations during the
simulation.

ConfigRecorder pushed to a `ConfigRecorder` objects are forwarded to
`compress(mc, model, conf)` for potential compression. When replaying a
simulation `decompress(mc, model, conf)` is called to undo this action. You may
define both to save on memory/disk space.

See also [`BufferedConfigRecorder`](@ref), [`Discarder`](@ref)
"""
struct ConfigRecorder{CT} <: AbstractRecorder
    configs::Vector{CT}
    rate::Int64
end
ConfigRecorder{CT}(rate::Integer = 10) where CT = ConfigRecorder{CT}(Vector{CT}(), rate)
function ConfigRecorder(MC::Type, M::Type, rate::Integer = 10)
    ConfigRecorder{compressed_conf_type(MC, M)}(rate)
end
function Base.push!(c::ConfigRecorder, mc, model, sweep)
    (sweep % c.rate == 0) && push!(c.configs, compress(mc, model, conf(mc)))
    nothing
end

# Base extensions
Base.length(c::ConfigRecorder) = length(c.configs)
Base.isempty(c::ConfigRecorder) = isempty(c.configs)
Base.lastindex(c::ConfigRecorder) = lastindex(c.configs)
Base.iterate(c::ConfigRecorder, i=1) = iterate(c.configs, i)
Base.getindex(c::ConfigRecorder, i) = c.configs[i]

# Compression
compress(mc, model, conf) = copy(conf)
decompress(mc, model, conf) = conf

function _save(file::JLDFile, cs::ConfigRecorder, entryname::String="configs")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/tag", "ConfigRecorder")
    write(file, entryname * "/type", typeof(cs))
    write(file, entryname * "/data", cs.configs)
    write(file, entryname * "/rate", cs.rate)
end
function _load(data, ::Val{:ConfigRecorder})
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load $T version $(data["VERSION"])"))
    end
    data["type"](data["data"], data["rate"])
end
to_tag(::Type{<: ConfigRecorder}) = Val(:ConfigRecorder)

################################################################################
### BufferedConfigRecorder
################################################################################

struct FilePath
    is_relative::Bool
    relative_path::String
    absolute_path::String
end
RelativePath(filename) = FilePath(true, filename, joinpath(pwd(), filename))
AbsolutePath(filename) = FilePath(false, filename, filename)

"""
    BufferedConfigRecorder(mc, model, filename[; rate = 10, chunk_size = 1000])

Creates a `BufferedConfigRecorder` object which accumulates configurations
during the simulation in an in-memory buffer. The buffer has a set `chunk_size`.
Every time it becomes it is commited to a JLD2 file `filename` and reset.

This uses the same configuration compression used in ConfigRecorder. I.e. every
`push!`ed configuration is passed to `compress(mc, model, conf)`.

See also [`ConfigRecorder`](@ref), [`Discarder`](@ref)
"""
mutable struct BufferedConfigRecorder{CT} <: AbstractRecorder
    filename::FilePath
    link_id::String
    buffer::Vector{CT}
    rate::Int64
    idx::Int64
    chunk::Int64
    total_length::Int64
    save_idx::Int64
end
function BufferedConfigRecorder(fn::String, link_id, buffer, rate, idx, chunk, N, sidx)
    BufferedConfigRecorder(AbsolutePath(fn), link_id, buffer, rate, idx, chunk, N, sidx)
end
function BufferedConfigRecorder{CT}(filename, rate = 10, chunk_size = 1000) where CT
    link_id = string(rand(UInt128))
    BufferedConfigRecorder(filename, link_id, Vector{CT}(undef, chunk_size), rate, 1, 1, 0, -1)
end
function BufferedConfigRecorder(MC::Type, M::Type, filename; rate = 10, chunk_size = 1000)
    BufferedConfigRecorder{compressed_conf_type(MC, M)}(filename, rate, chunk_size)
end
function Base.push!(cr::BufferedConfigRecorder, mc, model, sweep)
    if (sweep % cr.rate == 0)
        _push!(cr, compress(mc, model, conf(mc)))
    end
    nothing
end
function _push!(cr::BufferedConfigRecorder, data)
    if cr.idx == -1
        # The recorder is not working on the right chunk
        if cr.save_idx != cr.total_length
            @error "Stored data length does not match expected data length. Reseting head to last checkpoint."
            cr.total_length = cr.save_idx
        end
        # get chunk, idx of next writeable location
        chunk, idx = idx2chunk_idx(cr, cr.total_length+1)
        # load chunk if not loaded and required
        chunk != cr.chunk && idx != 1 && load_chunk!(cr, chunk)
        cr.chunk = chunk
        cr.idx = idx
    elseif cr.idx > length(cr.buffer)
        # The buffer is full - commit to file
        save_final_chunk!(cr)
        cr.chunk += 1
        cr.idx = 1
    end
    cr.buffer[cr.idx] = data
    cr.idx += 1
    cr.total_length += 1
    nothing
end

# Base extensions
Base.length(cr::BufferedConfigRecorder) = cr.total_length
Base.isempty(cr::BufferedConfigRecorder) = cr.total_length == 0
Base.lastindex(cr::BufferedConfigRecorder) = cr.total_length
Base.iterate(cr::BufferedConfigRecorder, i=1) = iterate(cr.configs, i)

function Base.getindex(cr::BufferedConfigRecorder, i)
    @boundscheck 1 <= i <= cr.total_length
    chunk, idx = idx2chunk_idx(cr, i)
    if chunk != cr.chunk
        # Save unsaved data in buffer
        cr.save_idx < cr.total_length && save_final_chunk!(cr)
        load_chunk!(cr, chunk)
    end
    return cr.buffer[idx]
end

# maps index to (chunk_idx, buffer_idx)
function idx2chunk_idx(cr::BufferedConfigRecorder, i)
    chunk_size = length(cr.buffer)
    chunk = div(i-1, chunk_size) + 1
    idx = mod1(i, chunk_size)
    return chunk, idx
end

# chunk loading and saving
function save_final_chunk!(cr::BufferedConfigRecorder)
    if cr.total_length != (cr.chunk - 1) * length(cr.buffer) + cr.idx - 1
        error("Attempting to save non-final chunk $(cr.chunk)!") # TODO this error triggered?
    end
    save_chunk!(cr)
    cr.save_idx = cr.total_length
    nothing
end

function save_chunk!(cr::BufferedConfigRecorder, chunk = cr.chunk)
    @boundscheck chunk > 0 && (chunk-1) * length(cr.buffer) < cr.total_length
    k = string(chunk)
    JLD2.jldopen(cr.filename.absolute_path, "a+", compress = true) do file
        if haskey(file, k)
            @debug "Replacing chunk $k"
            delete!(file, k)
        end
        file[k] = cr.buffer
    end
    nothing
end

function load_chunk!(cr::BufferedConfigRecorder, chunk)
    @boundscheck (chunk-1) * length(cr.buffer) < cr.save_idx
    JLD2.jldopen(cr.filename.absolute_path, "r") do file
        copyto!(cr.buffer, file[string(chunk)])
    end
    cr.chunk = chunk
    cr.idx = -1
end

function update_filepath!(cr::BufferedConfigRecorder, parent_path)
    # this only messes with relative filepaths
    # - update absolute path if necessary (parent path changed)
    # - move and rename on collision if link_id doesn't match
    # - replace if link_id matches

    if cr.filename.is_relative
        # get dir of parent savefile
        path, _ = splitdir(parent_path)

        # get adjusted absolute path, removing potential preceding /
        rp = cr.filename.relative_path
        filepath = joinpath(path, startswith(rp, '/') ? rp[2:end] : rp)

        # if absolute_path does not point to a file the parent has been moved to
        # a different system. In this case we assume that the file will be moved
        # as well.
        # if there is a file at absolute_path and the path differs from filepath
        # we should move the file.
        if isfile(cr.filename.absolute_path) && filepath != cr.filename.absolute_path

            # if there is already a file at filepath we may need to replace it.
            # If the id's don't match it belongs to a different simulation and
            # we rename this file, otherwise we replace.
            if isfile(filepath)
                target_link_id = try
                    JLD2.jldopen(filepath, "r") do f
                        get(f, "link_id", "N/A") # maybe default to current?
                    end
                catch e # if we can't read the link_id we do not want to touch the file
                    "ERROR"
                end

                if target_link_id != cr.link_id
                    new_filepath = _generate_unique_filename(filepath)
                    @warn(
                        "There already exists an independent file at " *
                        "$filepath. Renaming to $new_filepath to avoid " *
                        "overwriting other data."
                    )
                    filepath = new_filepath
                else
                    # we already have a matching file at the new location, so
                    # we don't need to replace anything
                    cr.filename = FilePath(true, rp, filepath)
                    return nothing
                end
            end

            mv(cr.filename.absolute_path, filepath, force = true)
        end

        # update filepath
        cr.filename = FilePath(true, string(filepath[length(path)+2:end]), filepath)
    end

    nothing
end

function _save(file, cr::BufferedConfigRecorder, entryname::String="configs")
    # save link_id
    JLD2.jldopen(cr.filename.absolute_path, "a+") do file
        if !haskey(file, "link_id")
            file["link_id"] = cr.link_id
        end
    end

    # save current buffer
    if !(cr.chunk == -1 && cr.idx == -1) # not in unitialized state
        if cr.save_idx != cr.total_length || !isfile(cr.filename.absolute_path)
            save_final_chunk!(cr)
        end
    end

    # adjust relative FilePath
    update_filepath!(cr, file.path)

    # main save information
    write(file, entryname * "/VERSION", 3)
    write(file, entryname * "/tag", "BufferedConfigRecorder")
    write(file, entryname * "/filename", cr.filename)
    write(file, entryname * "/link_id", cr.link_id)
    write(file, entryname * "/buffer", cr.buffer)
    write(file, entryname * "/rate", cr.rate)
    write(file, entryname * "/total_length", cr.total_length)
    write(file, entryname * "/save_idx", cr.total_length)
    nothing
end

function _load(data, ::Val{:BufferedConfigRecorder})
    if !(data["VERSION"] <= 3)
        throw(ErrorException("Failed to load BufferedConfigRecorder version $(data["VERSION"])"))
    end

    link_id = get(data, "link_id", "N/A")

    cr = BufferedConfigRecorder(
        data["filename"], link_id, data["buffer"], data["rate"], -1, -1,
        data["total_length"], data["save_idx"]
    )

    # adjust relative FilePath
    update_filepath!(cr, data.path)

    # if link_id unknown get it from file or generate new
    if link_id == "N/A" && isfile(cr.filename.absolute_path)
        cr.link_id = JLD2.jldopen(cr.filename.absolute_path, "r") do file
            haskey(file, "link_id") ? file["link_id"] : string(rand(UInt128))
        end
    end

    return cr
end

################################################################################
### Discarder
################################################################################

"""
    Discarder(args...; kwargs...)

Creates an object that discards all configurations pushed to it.

See also: [`BufferedConfigRecorder`](@ref), [`ConfigRecorder`](@ref)
"""
struct Discarder <: AbstractRecorder end
Discarder(args...; kwargs...) = Discarder()
Base.push!(::Discarder, args...) = nothing
Base.length(::Discarder) = 0
Base.isempty(::Discarder) = true
Base.getindex(v::Discarder, i) = BoundsError(v, i)
Base.iterate(c::Discarder, i=1) = nothing

function _save(file::JLDFile, ::Discarder, entryname::String="configs")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/tag", "Discarder")
end
_load(data, ::Val{:Discarder}) = Discarder()
to_tag(::Type{<: Discarder}) = Val(:Discarder)

################################################################################
### Utility
################################################################################

function Base.merge!(target::BufferedConfigRecorder, source::ConfigRecorder; rate = 1)
    for i in eachindex(source.configs)
        if i % rate == 0
            _push!(target, source.configs[i])
        end
    end
end


text:
using MarketData

facts("constructors") do

  context("inner constructor enforces constraints") do
      @fact 1 => roughly(1.0000000001)
  end

  context("from Blotter and FinancialTimeSeries") do
      @fact 1 => roughly(1.0000000001)
  end
end

facts("array of Trades") do

  context("closed trades") do
    @fact_throws log(-1)
  end

  context("opening trades") do
    @fact_throws log(-1)
  end

  context("create an array of trades") do
    @fact_throws log(-1)
  end
end


text:
"""
    Calculates the corrected satellite time
    $SIGNATURES

    ´sat_state´: satellite state, combining decoded data, code- and carrierphase

    Provides the implementation of the clock computations
    from raw to corrected

    Following instructions in IS-GPS-200K:
    20.3.3.3.3.1 User Algorithm for SV Clock Correction.
    Equation (1)
"""
function calc_corrected_time(sat_state::SatelliteState)

    t_sv = calc_uncorrected_time(sat_state)
    Δt_sv = code_phase_offset(sat_state.decoder_state, t_sv)
    t = t_sv - Δt_sv
end

"""
    Calculates the code Phase offset for time correction
    $SIGNATURES

    ´decoder_state´: Decoder containing ephemeris data of satellite
    ´t´: Current time(i.e. raw time without correction)

    Provides the implementation of the code phase clock
    correction term

    Following instructions in IS-GPS-200K:
    20.3.3.3.3.1 User Algorithm for SV Clock Correction.
    Equation (2)
"""
function code_phase_offset(decoder_state::GNSSDecoderState, t)
    Δtr = relativistic_correction(decoder_state, t)
    a_f0 = decoder_state.data.a_f0
    a_f1 = decoder_state.data.a_f1
    a_f2 = decoder_state.data.a_f2
    t_oc = decoder_state.data.t_oc
    Δt_sv = a_f0 + a_f1 * (t - t_oc) + a_f2 * (t - t_oc)^2 + Δtr
    Δt_sv = L1_correction(decoder_state, Δt_sv)
end

"""
    Provides the implementation of the relativistic clock
    correction term

    $SIGNATURES

    ´decoder_state´: Decoder containing ephemeris data of satellite
    ´t´: Current time(i.e. raw time without correction)

    Provides the implementation of the relativistic clock
    correction term

    Following instructions in IS-GPS-200K:
    # 20.3.3.3.3.1 User Algorithm for SV Clock Correction
"""
function relativistic_correction(decoder_state::GNSSDecoderState, t, dt = 0)

    F = decoder_state.constants.F
    e = decoder_state.data.e
    sqrtA = decoder_state.data.sqrt_A
    Ek = calc_eccentric_anomaly(decoder_state, t, dt)
    Δtr = F * e * sqrtA * sin(Ek)
end

"""
    Calculates raw time data by observed times
    $SIGNATURES

    ´sat_state´: satellite state, combining decoded data, code- and carrierphase

    Calculates time by number of bits since TOW, TOW and code and carrier phases
"""
function calc_uncorrected_time(sat_state::SatelliteState)
    gpsl1 = GNSSSignals.GPSL1()
    t_tow = sat_state.decoder_state.data.TOW * 6
    t_bits = sat_state.decoder_state.num_bits_buffered / GNSSSignals.get_data_frequency(gpsl1) * Hz
    t_code_phase = sat_state.code_phase / GNSSSignals.get_code_frequency(gpsl1) * Hz
    t_carrier_phase = sat_state.carrier_phase / GNSSSignals.get_center_frequency(gpsl1) * Hz

    t = t_tow + t_bits + t_code_phase + t_carrier_phase
end

"""
    Provides the implementation of the Group Delay correction
    $SIGNATURES
    ´decoder_state´: Decoder
    ´Δt_sv´: code phase offset term

    Following instructions in IS-GPS-200K:
    20.3.3.3.3.2 L1 or L2 Correction.
"""
function L1_correction(decoder_state::GNSSDecoderState, Δt_sv)
    return Δt_sv - decoder_state.data.T_GD
end

"""
Calculates the eccentric anomaly
$SIGNATURES
´decoder_state´: Decoder storing ephemeris data
´t_sv´: time of satellite vehicle
´dt_sv´: time correction for SV

This function implements the calculation of the Eccentric
anomaly, which is required for both, the calculation of the
satellite position and the time correction.

Following instructions in IS-GPS-200K:
Table 20-IV
"""
function calc_eccentric_anomaly(decoder_state, t_sv, dt_sv = 0 )

    t = t_sv - dt_sv

    tk = check_crossover(t - decoder_state.data.t_oe)
    μ = decoder_state.constants.μ
    A = decoder_state.data.sqrt_A^2

    n_0 = sqrt(μ / (A^3))
    n = n_0 + decoder_state.data.Δn
    M = decoder_state.data.M_0 + n * tk

    e = decoder_state.data.e
    E = M
    for k = 1:30
        Et = E
        E = M + e * sin(E)
        if abs(E - Et) <= 1e-12
            break;
        end
    end
    return E
end

"""
Correction for end of week crossovers
$SIGNATURES

´t´: Time to be controlled

Following instructions in IS-GPS-200K:
20.3.3.3.3.1 User Algorithm for SV Clock Correction.
"""
function check_crossover(t)

    half_week = 302400.0  #Half of Week in Seconds
    if (t > half_week)
        t = t - 2.0 * half_week
    elseif (t < -half_week)
        t = t + 2.0 * half_week
    end
    return t
end