text:
@testset "Floats" begin

res = Parsers.xparse(Float64, "1x")

testcases = [
    (str="", x=0.0, code=(INVALID | EOF), len=0, tot=0),
    (str="-", x=0.0, code=(INVALID | EOF), len=1, tot=1),
    (str=".", x=0.0, code=(INVALID | EOF), len=1, tot=1),
    (str="1", x=1.0, code=(OK | EOF), len=1, tot=1),
    (str="+1", x=1.0, code=(OK | EOF), len=2, tot=2),
    (str="1.1", x=1.1, code=(OK | EOF), len=3, tot=3),
    (str="1.", x=1.0, code=(OK | EOF), len=2, tot=2),
    (str="-1.", x=-1.0, code=(OK | EOF), len=3, tot=3),
    (str="-1. ", x=-1.0, code=(OK | EOF), len=4, tot=4),
    (str="1.a", x=1.0, code=(OK | EOF | INVALID_DELIMITER), len=2, tot=2),
    (str="1.1.", x=1.1, code=(OK | EOF | INVALID_DELIMITER), len=3, tot=3),
    (str="1e23", x=1e23, code=(OK | EOF), len=4, tot=4),
    (str="1E23", x=1e23, code=(OK | EOF), len=4, tot=4),
    (str="1f23", x=1e23, code=(OK | EOF), len=4, tot=4),
    (str="1F23", x=1e23, code=(OK | EOF), len=4, tot=4),
    (str="428.E+03", x=428000.0, code=(OK | EOF), len=8, tot=8),
    (str="1.0e", x=0.0, code=(INVALID | EOF), len=4, tot=4),
    (str="1.0ea", x=0.0, code=(INVALID | EOF | INVALID_DELIMITER), len=5, tot=5),
    (str="100000000000000000000000", x=1e23, code=(OK | EOF), len=24, tot=24),
    (str="1e-100", x=1e-100, code=(OK | EOF), len=6, tot=6),
    (str="123456700", x=1.234567e+08, code=(OK | EOF), len=9, tot=9),
    (str="99999999999999974834176", x=9.999999999999997e+22, code=(OK | EOF), len=23, tot=23),
    (str="100000000000000000000001", x=1.0000000000000001e+23, code=(OK | EOF), len=24, tot=24),
    (str="100000000000000008388608", x=1.0000000000000001e+23, code=(OK | EOF), len=24, tot=24),
    (str="100000000000000016777215", x=1.0000000000000001e+23, code=(OK | EOF), len=24, tot=24),
    (str="100000000000000016777216", x=1.0000000000000003e+23, code=(OK | EOF), len=0, tot=0),

    (str="-1", x=-1.0, code=(OK | EOF), len=0, tot=0),
    (str="-0.1", x=-0.1, code=(OK | EOF), len=0, tot=0),
    (str="-0", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="1e-20", x=1e-20, code=(OK | EOF), len=0, tot=0),
    (str="625e-3", x=0.625, code=(OK | EOF), len=0, tot=0),

    # zeros
    (str="0", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="0e0", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="-0e0", x=-0.0, code=(OK | EOF), len=0, tot=0),
    (str="+0e0", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="0e-0", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="-0e-0", x=-0.0, code=(OK | EOF), len=0, tot=0),
    (str="+0e-0", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="0e+0", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="-0e+0", x=-0.0, code=(OK | EOF), len=0, tot=0),
    (str="+0e+0", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="0e+01234567890123456789", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="0.00e-01234567890123456789", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="-0e+01234567890123456789", x=-0.0, code=(OK | EOF), len=0, tot=0),
    (str="-0.00e-01234567890123456789", x=-0.0, code=(OK | EOF), len=0, tot=0),
    (str="0e291", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="0e292", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="0e347", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="0e348", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="-0e291", x=-0.0, code=(OK | EOF), len=0, tot=0),
    (str="-0e292", x=-0.0, code=(OK | EOF), len=0, tot=0),
    (str="-0e347", x=-0.0, code=(OK | EOF), len=0, tot=0),
    (str="-0e348", x=-0.0, code=(OK | EOF), len=0, tot=0),

    # Infs/invalid NaNs
    (str="n", x=0.0, code=(INVALID | EOF), len=1, tot=1),
    (str="na", x=0.0, code=(INVALID | EOF), len=2, tot=2),
    (str="n", x=0.0, code=(INVALID | EOF), len=1, tot=1),
    (str="inf", x=Inf, code=(OK | EOF), len=0, tot=0),
    (str="infinity", x=Inf, code=(OK | EOF), len=0, tot=0),
    (str="-inf", x=-Inf, code=(OK | EOF), len=0, tot=0),
    (str="-Inf", x=-Inf, code=(OK | EOF), len=0, tot=0),
    (str="-infinity", x=-Inf, code=(OK | EOF), len=0, tot=0),
    (str="-INFINITY", x=-Inf, code=(OK | EOF), len=0, tot=0),
    (str="+inf", x=Inf, code=(OK | EOF), len=0, tot=0),
    (str="i", x=0.0, code=(INVALID | EOF), len=1, tot=1),
    (str="in", x=0.0, code=(INVALID | EOF), len=2, tot=2),
    (str="infi", x=Inf, code=(OK | EOF), len=1, tot=1),
    (str="infin", x=Inf, code=(OK | EOF), len=1, tot=1),
    (str="infini", x=Inf, code=(OK | EOF), len=1, tot=1),
    (str="infinit", x=Inf, code=(OK | EOF), len=1, tot=1),
    (str="i,", x=0.0, code=(INVALID | DELIMITED | INVALID_DELIMITER), len=1, tot=1),
    (str="in,", x=0.0, code=(INVALID | DELIMITED | INVALID_DELIMITER), len=2, tot=2),
    (str="infi,", x=Inf, code=(OK | DELIMITED), len=1, tot=1),
    (str="infin,", x=Inf, code=(OK | DELIMITED), len=1, tot=1),
    (str="infini,", x=Inf, code=(OK | DELIMITED), len=1, tot=1),
    (str="infinit,", x=Inf, code=(OK | DELIMITED), len=1, tot=1),

    # largest float64
    (str="1.7976931348623157e308", x=1.7976931348623157e+308, code=(OK | EOF), len=0, tot=0),
    (str="-1.7976931348623157e308", x=-1.7976931348623157e+308, code=(OK | EOF), len=0, tot=0),
    # next float64 - too large
    (str="1.7976931348623159e308", x=+Inf, code=(OK | EOF), len=0, tot=0),
    (str="-1.7976931348623159e308", x=-Inf, code=(OK | EOF), len=0, tot=0),
    # the border is ...158079
    # borderline - okay
    (str="1.7976931348623158e308", x=1.7976931348623157e+308, code=(OK | EOF), len=0, tot=0),
    (str="-1.7976931348623158e308", x=-1.7976931348623157e+308, code=(OK | EOF), len=0, tot=0),
    # borderline - too large
    (str="1.797693134862315808e308", x=+Inf, code=(OK | EOF), len=0, tot=0),
    (str="-1.797693134862315808e308", x=-Inf, code=(OK | EOF), len=0, tot=0),

    # a little too large
    (str="1e308", x=1e+308, code=(OK | EOF), len=0, tot=0),
    (str="2e308", x=+Inf, code=(OK | EOF), len=0, tot=0),
    (str="1e309", x=+Inf, code=(OK | EOF), len=0, tot=0),

    # way too large
    (str="1e310", x=+Inf, code=(OK | EOF), len=0, tot=0),
    (str="-1e310", x=-Inf, code=(OK | EOF), len=0, tot=0),
    (str="1e400", x=+Inf, code=(OK | EOF), len=0, tot=0),
    (str="-1e400", x=-Inf, code=(OK | EOF), len=0, tot=0),
    (str="1e400000", x=+Inf, code=(OK | EOF), len=0, tot=0),
    (str="-1e400000", x=-Inf, code=(OK | EOF), len=0, tot=0),

    # denormalized
    (str="1e-305", x=1e-305, code=(OK | EOF), len=0, tot=0),
    (str="1e-306", x=1e-306, code=(OK | EOF), len=0, tot=0),
    (str="1e-307", x=1e-307, code=(OK | EOF), len=0, tot=0),
    (str="1e-308", x=1e-308, code=(OK | EOF), len=0, tot=0),
    (str="1e-309", x=1e-309, code=(OK | EOF), len=0, tot=0),
    (str="1e-310", x=1e-310, code=(OK | EOF), len=0, tot=0),
    (str="1e-322", x=1e-322, code=(OK | EOF), len=0, tot=0),
    # smallest denormal
    (str="5e-324", x=5e-324, code=(OK | EOF), len=0, tot=0),
    (str="4e-324", x=5e-324, code=(OK | EOF), len=0, tot=0),
    (str="3e-324", x=5e-324, code=(OK | EOF), len=0, tot=0),
    # too small
    (str="2e-324", x=0.0, code=(OK | EOF), len=0, tot=0),
    # way too small
    (str="1e-350", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="1e-400000", x=0.0, code=(OK | EOF), len=0, tot=0),

    # try to overflow exponent
    (str="1e-4294967296", x=0.0, code=(OK | EOF), len=0, tot=0),
    (str="1e+4294967296", x=+Inf, code=(OK | EOF), len=0, tot=0),
    (str="1e-18446744073709551616", x=0.0, code=(OK | EOF), len=0, tot=0),

    (str="1e+18446744073709551616", x=+Inf, code=(OK | EOF), len=0, tot=0),

    # Parse errors
    # (str="1e", x=missing, code=(INVALID | EOF), len=0, tot=0),
    # (str="1e-", x=missing, code=(INVALID | EOF), len=0, tot=0),
    # (str=".e-1", x=missing, code=(INVALID | EOF), len=0, tot=0),
    # (str="1\x00.2", x=1.0, code=(OK), len=0, tot=0),

    # http:#www.exploringbinary.com/java-hangs-when-converting-2-2250738585072012e-308/
    (str="2.2250738585072012e-308", x=2.2250738585072014e-308, code=(OK | EOF), len=0, tot=0),
    # http://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/
    (str="2.2250738585072011e-308", x=2.225073858507201e-308, code=(OK | EOF), len=0, tot=0),

    # A very large number (initially wrongly parsed by the fast algorithm).
    (str="4.630813248087435e+307", x=4.630813248087435e+307, code=(OK | EOF), len=0, tot=0),

    # A different kind of very large number.
    (str="22.222222222222222", x=22.22222222222222, code=(OK | EOF), len=0, tot=0),

    # Exactly halfway between 1 and math.Nextafter(1, 2).
    # Round to even (down).
    (str="1.00000000000000011102230246251565404236316680908203125", x=1.0, code=(OK | EOF), len=0, tot=0),
    # # Slightly lower; still round down.
    (str="1.00000000000000011102230246251565404236316680908203124", x=1.0, code=(OK | EOF), len=0, tot=0),
    # # Slightly higher; round up.
    (str="1.00000000000000011102230246251565404236316680908203126", x=1.0000000000000002, code=(OK | EOF), len=0, tot=0),

    (str="-5.871153289887625082e-01", x=-5.871153289887625082e-01, code=(OK | EOF), len=0, tot=0),
    (str="8.095032986136727615e-01", x=8.095032986136727615e-01, code=(OK | EOF), len=0, tot=0),
    (str="9.900000000000006573e-01", x=9.900000000000006573e-01, code=(OK | EOF), len=0, tot=0),
    (str="9.900000000000006573e-01", x=9.900000000000006573e-01, code=(OK | EOF), len=0, tot=0),
    (str="-9.866066418838319585e-01", x=-9.866066418838319585e-01, code=(OK | EOF), len=0, tot=0),
    (str="-3.138907529596844714e+00", x=-3.138907529596844714e+00, code=(OK | EOF), len=0, tot=0),
    (str="-5.129218106887371675e+00", x=-5.129218106887371675e+00, code=(OK | EOF), len=0, tot=0),
    (str="-4.803915800699462224e+00", x=-4.803915800699462224e+00, code=(OK | EOF), len=0, tot=0),

    # issue #18
    (str=".24409E+03", x=244.09, code=(OK | EOF), len=0, tot=0),
    (str=".24409E+03 ", x=244.09, code=(OK | EOF), len=0, tot=0),
    (str="-.2", x=-0.2, code=(OK | EOF), len=0, tot=0),
    (str=".2", x=0.2, code=(OK | EOF), len=0, tot=0),

    # from https://www.icir.org/vern/papers/testbase-report.pdf
    (str=string(5.0 * exp10(+125)), x=(5.0 * exp10(+125)), code=(OK | EOF), len=0, tot=0),
    (str=string(69.0 * exp10(+267)), x=(69.0 * exp10(+267)), code=(OK | EOF), len=0, tot=0),
    (str=string(999.0 * exp10(-26)), x=(999.0 * exp10(-26)), code=(OK | EOF), len=0, tot=0),
    (str=string(7861.0 * exp10(-34)), x=(7861.0 * exp10(-34)), code=(OK | EOF), len=0, tot=0),
    (str=string(75569.0 * exp10(-254)), x=(75569.0 * exp10(-254)), code=(OK | EOF), len=0, tot=0),
    (str=string(928609.0 * exp10(-261)), x=(928609.0 * exp10(-261)), code=(OK | EOF), len=0, tot=0),
    (str=string(9210917.0 * exp10(+80)), x=(9210917.0 * exp10(+80)), code=(OK | EOF), len=0, tot=0),
    (str=string(84863171.0 * exp10(+114)), x=(84863171.0 * exp10(+114)), code=(OK | EOF), len=0, tot=0),
    (str=string(653777767.0 * exp10(+273)), x=(653777767.0 * exp10(+273)), code=(OK | EOF), len=0, tot=0),
    (str=string(5232604057.0 * exp10(-298)), x=(5232604057.0 * exp10(-298)), code=(OK | EOF), len=0, tot=0),
    (str=string(27235667517.0 * exp10(-109)), x=(27235667517.0 * exp10(-109)), code=(OK | EOF), len=0, tot=0),
    (str=string(653532977297.0 * exp10(-123)), x=(653532977297.0 * exp10(-123)), code=(OK | EOF), len=0, tot=0),
    (str=string(3142213164987.0 * exp10(-294)), x=(3142213164987.0 * exp10(-294)), code=(OK | EOF), len=0, tot=0),
    (str=string(46202199371337.0 * exp10(-72)), x=(46202199371337.0 * exp10(-72)), code=(OK | EOF), len=0, tot=0),
    (str=string(231010996856685.0 * exp10(-73)), x=(231010996856685.0 * exp10(-73)), code=(OK | EOF), len=0, tot=0),
    (str=string(9324754620109615.0 * exp10(+212)), x=(9324754620109615.0 * exp10(+212)), code=(OK | EOF), len=0, tot=0),
    (str=string(78459735791271921.0 * exp10(+49)), x=(78459735791271921.0 * exp10(+49)), code=(OK | EOF), len=0, tot=0),
    (str=string(272104041512242479.0 * exp10(+200)), x=(272104041512242479.0 * exp10(+200)), code=(OK | EOF), len=0, tot=0),
    (str=string(6802601037806061975.0 * exp10(+198)), x=(6802601037806061975.0 * exp10(+198)), code=(OK | EOF), len=0, tot=0),
    (str=string(20505426358836677347.0 * exp10(-221)), x=(20505426358836677347.0 * exp10(-221)), code=(OK | EOF), len=0, tot=0),
    (str=string(836168422905420598437.0 * exp10(-234)), x=(836168422905420598437.0 * exp10(-234)), code=(OK | EOF), len=0, tot=0),
    (str=string(4891559871276714924261.0 * exp10(+222)), x=(4891559871276714924261.0 * exp10(+222)), code=(OK | EOF), len=0, tot=0),

    (str=string(9.0 * exp10(-265)), x=(9.0 * exp10(-265)), code=(OK | EOF), len=0, tot=0),
    (str=string(85.0 * exp10(-037)), x=(85.0 * exp10(-037)), code=(OK | EOF), len=0, tot=0),
    (str=string(623.0 * exp10(+100)), x=(623.0 * exp10(+100)), code=(OK | EOF), len=0, tot=0),
    (str=string(3571.0 * exp10(+263)), x=(3571.0 * exp10(+263)), code=(OK | EOF), len=0, tot=0),
    (str=string(81661.0 * exp10(+153)), x=(81661.0 * exp10(+153)), code=(OK | EOF), len=0, tot=0),
    (str=string(920657.0 * exp10(-023)), x=(920657.0 * exp10(-023)), code=(OK | EOF), len=0, tot=0),
    (str=string(4603285.0 * exp10(-024)), x=(4603285.0 * exp10(-024)), code=(OK | EOF), len=0, tot=0),
    (str=string(87575437.0 * exp10(-309)), x=(87575437.0 * exp10(-309)), code=(OK | EOF), len=0, tot=0),
    (str=string(245540327.0 * exp10(+122)), x=(245540327.0 * exp10(+122)), code=(OK | EOF), len=0, tot=0),
    (str=string(6138508175.0 * exp10(+120)), x=(6138508175.0 * exp10(+120)), code=(OK | EOF), len=0, tot=0),
    (str=string(83356057653.0 * exp10(+193)), x=(83356057653.0 * exp10(+193)), code=(OK | EOF), len=0, tot=0),
    (str=string(619534293513.0 * exp10(+124)), x=(619534293513.0 * exp10(+124)), code=(OK | EOF), len=0, tot=0),
    (str=string(2335141086879.0 * exp10(+218)), x=(2335141086879.0 * exp10(+218)), code=(OK | EOF), len=0, tot=0),
    (str=string(36167929443327.0 * exp10(-159)), x=(36167929443327.0 * exp10(-159)), code=(OK | EOF), len=0, tot=0),
    (str=string(609610927149051.0 * exp10(-255)), x=(609610927149051.0 * exp10(-255)), code=(OK | EOF), len=0, tot=0),
    (str=string(3743626360493413.0 * exp10(-165)), x=(3743626360493413.0 * exp10(-165)), code=(OK | EOF), len=0, tot=0),
    (str=string(94080055902682397.0 * exp10(-242)), x=(94080055902682397.0 * exp10(-242)), code=(OK | EOF), len=0, tot=0),
    (str=string(899810892172646163.0 * exp10(+283)), x=(899810892172646163.0 * exp10(+283)), code=(OK | EOF), len=0, tot=0),
    (str=string(7120190517612959703.0 * exp10(+120)), x=(7120190517612959703.0 * exp10(+120)), code=(OK | EOF), len=0, tot=0),
    (str=string(25188282901709339043.0 * exp10(-252)), x=(25188282901709339043.0 * exp10(-252)), code=(OK | EOF), len=0, tot=0),
    (str=string(308984926168550152811.0 * exp10(-052)), x=(308984926168550152811.0 * exp10(-052)), code=(OK | EOF), len=0, tot=0),
    (str=string(6372891218502368041059.0 * exp10(+064)), x=(6372891218502368041059.0 * exp10(+064)), code=(OK | EOF), len=0, tot=0),

    (str=string(8511030020275656.0 * exp2(-0342)), x=(8511030020275656.0 * exp2(-0342)), code=(OK | EOF), len=0, tot=0),
    (str=string(5201988407066741.0 * exp2(-0824)), x=(5201988407066741.0 * exp2(-0824)), code=(OK | EOF), len=0, tot=0),
    (str=string(6406892948269899.0 * exp2(+0237)), x=(6406892948269899.0 * exp2(+0237)), code=(OK | EOF), len=0, tot=0),
    (str=string(8431154198732492.0 * exp2(+0072)), x=(8431154198732492.0 * exp2(+0072)), code=(OK | EOF), len=0, tot=0),
    (str=string(6475049196144587.0 * exp2(+0099)), x=(6475049196144587.0 * exp2(+0099)), code=(OK | EOF), len=0, tot=0),
    (str=string(8274307542972842.0 * exp2(+0726)), x=(8274307542972842.0 * exp2(+0726)), code=(OK | EOF), len=0, tot=0),
    (str=string(5381065484265332.0 * exp2(-0456)), x=(5381065484265332.0 * exp2(-0456)), code=(OK | EOF), len=0, tot=0),
    (str=string(6761728585499734.0 * exp2(-1057)), x=(6761728585499734.0 * exp2(-1057)), code=(OK | EOF), len=0, tot=0),
    (str=string(7976538478610756.0 * exp2(+0376)), x=(7976538478610756.0 * exp2(+0376)), code=(OK | EOF), len=0, tot=0),
    (str=string(5982403858958067.0 * exp2(+0377)), x=(5982403858958067.0 * exp2(+0377)), code=(OK | EOF), len=0, tot=0),
    (str=string(5536995190630837.0 * exp2(+0093)), x=(5536995190630837.0 * exp2(+0093)), code=(OK | EOF), len=0, tot=0),
    (str=string(7225450889282194.0 * exp2(+0710)), x=(7225450889282194.0 * exp2(+0710)), code=(OK | EOF), len=0, tot=0),
    (str=string(7225450889282194.0 * exp2(+0709)), x=(7225450889282194.0 * exp2(+0709)), code=(OK | EOF), len=0, tot=0),
    (str=string(8703372741147379.0 * exp2(+0117)), x=(8703372741147379.0 * exp2(+0117)), code=(OK | EOF), len=0, tot=0),
    (str=string(8944262675275217.0 * exp2(-1001)), x=(8944262675275217.0 * exp2(-1001)), code=(OK | EOF), len=0, tot=0),
    (str=string(7459803696087692.0 * exp2(-0707)), x=(7459803696087692.0 * exp2(-0707)), code=(OK | EOF), len=0, tot=0),
    (str=string(6080469016670379.0 * exp2(-0381)), x=(6080469016670379.0 * exp2(-0381)), code=(OK | EOF), len=0, tot=0),
    (str=string(8385515147034757.0 * exp2(+0721)), x=(8385515147034757.0 * exp2(+0721)), code=(OK | EOF), len=0, tot=0),
    (str=string(7514216811389786.0 * exp2(-0828)), x=(7514216811389786.0 * exp2(-0828)), code=(OK | EOF), len=0, tot=0),
    (str=string(8397297803260511.0 * exp2(-0345)), x=(8397297803260511.0 * exp2(-0345)), code=(OK | EOF), len=0, tot=0),
    (str=string(6733459239310543.0 * exp2(+0202)), x=(6733459239310543.0 * exp2(+0202)), code=(OK | EOF), len=0, tot=0),
    (str=string(8091450587292794.0 * exp2(-0473)), x=(8091450587292794.0 * exp2(-0473)), code=(OK | EOF), len=0, tot=0),

    (str=string(6567258882077402.0 * exp2(+952)), x=(6567258882077402.0 * exp2(+952)), code=(OK | EOF), len=0, tot=0),
    (str=string(6712731423444934.0 * exp2(+535)), x=(6712731423444934.0 * exp2(+535)), code=(OK | EOF), len=0, tot=0),
    (str=string(6712731423444934.0 * exp2(+534)), x=(6712731423444934.0 * exp2(+534)), code=(OK | EOF), len=0, tot=0),
    (str=string(5298405411573037.0 * exp2(-957)), x=(5298405411573037.0 * exp2(-957)), code=(OK | EOF), len=0, tot=0),
    (str=string(5137311167659507.0 * exp2(-144)), x=(5137311167659507.0 * exp2(-144)), code=(OK | EOF), len=0, tot=0),
    (str=string(6722280709661868.0 * exp2(+363)), x=(6722280709661868.0 * exp2(+363)), code=(OK | EOF), len=0, tot=0),
    (str=string(5344436398034927.0 * exp2(-169)), x=(5344436398034927.0 * exp2(-169)), code=(OK | EOF), len=0, tot=0),
    (str=string(8369123604277281.0 * exp2(-853)), x=(8369123604277281.0 * exp2(-853)), code=(OK | EOF), len=0, tot=0),
    (str=string(8995822108487663.0 * exp2(-780)), x=(8995822108487663.0 * exp2(-780)), code=(OK | EOF), len=0, tot=0),
    (str=string(8942832835564782.0 * exp2(-383)), x=(8942832835564782.0 * exp2(-383)), code=(OK | EOF), len=0, tot=0),
    (str=string(8942832835564782.0 * exp2(-384)), x=(8942832835564782.0 * exp2(-384)), code=(OK | EOF), len=0, tot=0),
    (str=string(8942832835564782.0 * exp2(-385)), x=(8942832835564782.0 * exp2(-385)), code=(OK | EOF), len=0, tot=0),
    (str=string(6965949469487146.0 * exp2(-249)), x=(6965949469487146.0 * exp2(-249)), code=(OK | EOF), len=0, tot=0),
    (str=string(6965949469487146.0 * exp2(-250)), x=(6965949469487146.0 * exp2(-250)), code=(OK | EOF), len=0, tot=0),
    (str=string(6965949469487146.0 * exp2(-251)), x=(6965949469487146.0 * exp2(-251)), code=(OK | EOF), len=0, tot=0),
    (str=string(7487252720986826.0 * exp2(+548)), x=(7487252720986826.0 * exp2(+548)), code=(OK | EOF), len=0, tot=0),
    (str=string(5592117679628511.0 * exp2(+164)), x=(5592117679628511.0 * exp2(+164)), code=(OK | EOF), len=0, tot=0),
    (str=string(8887055249355788.0 * exp2(+665)), x=(8887055249355788.0 * exp2(+665)), code=(OK | EOF), len=0, tot=0),
    (str=string(6994187472632449.0 * exp2(+690)), x=(6994187472632449.0 * exp2(+690)), code=(OK | EOF), len=0, tot=0),
    (str=string(8797576579012143.0 * exp2(+588)), x=(8797576579012143.0 * exp2(+588)), code=(OK | EOF), len=0, tot=0),
    (str=string(7363326733505337.0 * exp2(+272)), x=(7363326733505337.0 * exp2(+272)), code=(OK | EOF), len=0, tot=0),
    (str=string(8549497411294502.0 * exp2(-448)), x=(8549497411294502.0 * exp2(-448)), code=(OK | EOF), len=0, tot=0),

    # from lemire's test suite
    (str="7.3177701707893310e+15", x=7.3177701707893310e+15, code=(OK | EOF), len=0, tot=0),
    (str="9007199254740995", x=9.007199254740996e15, code=(OK | EOF), len=0, tot=0),
    (str="7e23", x=7e23, code=(OK | EOF), len=0, tot=0),
    (str="-65.613616999999977", x=-65.613616999999977, code=(OK | EOF), len=0, tot=0),
    (str="7.2057594037927933e+16", x=7.2057594037927933e+16, code=(OK | EOF), len=0, tot=0),
    (str="0.1e-308", x=0.1e-308, code=(OK | EOF), len=0, tot=0),
    (str="0.01e-307", x=0.01e-307, code=(OK | EOF), len=0, tot=0),
    (str="1.79769e+308", x=1.79769e+308, code=(OK | EOF), len=0, tot=0),
    (str="2.22507e-308", x=2.22507e-308, code=(OK | EOF), len=0, tot=0),
    (str="-1.79769e+308", x=-1.79769e+308, code=(OK | EOF), len=0, tot=0),
    (str="-2.22507e-308", x=-2.22507e-308, code=(OK | EOF), len=0, tot=0),
];

for useio in (false, true)
    for (i, case) in enumerate(testcases)
        # println("testing case = $i, str = $(case.str)")
        res = Parsers.xparse(Float64, useio ? IOBuffer(case.str) : case.str)
        x, code, tlen = res.val, res.code, res.tlen
        if !Parsers.invalid(code)
            @test x == case.x
        end
        @test code == case.code
    end
end

# NaNs
case = (str="nan", x=NaN, code=(OK | EOF), len=3, tot=3)
res = Parsers.xparse(Float64, case.str)
x, code, tlen = res.val, res.code, res.tlen
@test isnan(x)
@test code == case.code
@test tlen == 3
case = (str="NaN", x=NaN, code=(OK | EOF), len=3, tot=3)
res = Parsers.xparse(Float64, case.str)
x, code, tlen = res.val, res.code, res.tlen
@test isnan(x)
@test code == case.code
@test tlen == 3
case = (str="NAN", x=NaN, code=(OK | EOF), len=3, tot=3)
res = Parsers.xparse(Float64, case.str)
x, code, tlen = res.val, res.code, res.tlen
@test isnan(x)
@test code == case.code
@test tlen == 3
case = (str="nAN", x=NaN, code=(OK | EOF), len=3, tot=3)
res = Parsers.xparse(Float64, case.str)
x, code, tlen = res.val, res.code, res.tlen
@test isnan(x)
@test code == case.code
@test tlen == 3
case = (str="nAN,", x=NaN, code=(OK | DELIMITED), len=4, tot=4)
res = Parsers.xparse(Float64, case.str)
x, code, tlen = res.val, res.code, res.tlen
@test isnan(x)
@test code == case.code
@test tlen == 4

# #25
case = (str="74810199.033988851037472901827191090834", x=7.481019903398885e7, code=(OK | EOF))
res = Parsers.xparse(Float64, case.str)
x, code, tlen = res.val, res.code, res.tlen
@test x == case.x
@test code == case.code
@test tlen == 39

# https://github.com/quinnj/JSON3.jl/issues/57
@test_throws Parsers.Error Parsers.parse(Float64,join(rand(1:9, 2000), ""))
@test Parsers.parse(BigFloat, "3.14") == BigFloat("3.14")
@test Parsers.parse(BigFloat, "3.14 ") == BigFloat("3.14")

# https://github.com/JuliaData/Parsers.jl/issues/53
bytes = codeunits("a,b,c\n.,1,3")
res = Parsers.xparse(Float64, bytes, 7, 11)
x, code, tlen = res.val, res.code, res.tlen
@test code == (INVALID | DELIMITED)

# https://github.com/JuliaData/CSV.jl/issues/710
# problematic really small floats on various threshold boundaries
@test Parsers.parse(Float64, "9.88e-324") === 1.0e-323
@test Parsers.parse(Float64, "4.94e-324") === 5.0e-324
@test Parsers.parse(Float64, "8.40e-323") === 8.4e-323
@test Parsers.parse(Float64, "3091.") === 3091.0

# https://github.com/JuliaData/CSV.jl/issues/769
@test Parsers.parse(Float64, "9223372036854775808") === 9.223372036854776e18

# random sampling of Float16 testing
for _ = 1:1000
    x = rand(Float16)
    str = string(x)
    @test Base.parse(Float16, str) === Parsers.parse(Float16, str)
end

# discovered from JSON tests
@test Parsers.tryparse(Float64, "0e+") === nothing

# https://github.com/JuliaData/CSV.jl/issues/916
@test  Parsers.parse(Float64, "0.44311945001372019574271437679879349172") === 0.4431194500137202
@test Parsers.parse(BigFloat, "0.44311945001372019574271437679879349172") == BigFloat("0.44311945001372019574271437679879349172")

# https://github.com/JuliaData/Parsers.jl/issues/83
@test Parsers.parse(Float64,"9.3494547075363499E-311") === 9.3494547075363e-311
@test Parsers.parse(Float64,"9.349454707536349999E-311") === 9.3494547075363e-311

# https://github.com/JuliaData/CSV.jl/issues/938
v1 = Parsers.parse(BigFloat, "0.994322841995507271075540969506")
v2 = Parsers.parse(BigFloat, "0.794322841995507271075540969506")
@test v1 !== v2

# parsing larger than Int64 with dangling non-digits
@test_throws Parsers.Error Parsers.parse(Float64, "3130307457683247493156627A")

end # @testset


text:
using PyCall
anim = pyimport("matplotlib.animation")
using DataFrames
using DataFramesMeta
import Interpolations: interpolate, extrapolate, BSpline, Constant, Flat
import ImageFiltering: imfilter, Kernel

import MapImages: nearest_pixel, nearest, grid, km_to_deg
import Geodesy: LLA, nad27, wgs84, UTMZfromLLA, UTMfromLLA
import KernelDensity.kde

kde(df::DataFrame, xcol::Symbol, ycol::Symbol; kwargs...) =
    kde((df[!, xcol], df[!, ycol]); kwargs...)
kde(df::DataFrame, xcol::Symbol, ycol::Symbol, weightcol::Symbol; kwargs...) =
    kde((df[!, xcol], df[!, ycol]); weights = df[!, weightcol], kwargs...)

ccrs = pyimport("cartopy.crs")
shpreader = pyimport("cartopy.io.shapereader")
sgeom = pyimport("shapely.geometry")
cfeature = pyimport("cartopy.feature")

function project(lons, lats; datum = nad27, zone = 13, north = true)
    txutm = UTMZfromLLA(datum)
    # txutm = UTMfromLLA(zone, north, datum)
    pts = txutm.(LLA.(lats, lons, 0))
    # Note: these will have units of meters
    return [p.x for p in pts], [p.y for p in pts]
end

project(df::DataFrame, xcol::Symbol, ycol::Symbol; datum = nad27) =
    project(df[!, xcol], df[!, ycol], datum = datum)

function kde_image(
    df,
    xcol,
    ycol,
    weightcol;
    unnormalize = true,
    extent = nothing,
    boundary = nothing,
    do_project = true,
)
    xx, yy = do_project ? project(df, xcol, ycol) : (df[!, xcol], df[!, ycol])

    if isnothing(boundary)
        boundary = isnothing(extent) ? (extrema(xx), extrema(yy)) :
            ((extent[1], extent[2]), (extent[3], extent[4]))
    end

    # k = kde(df, xcol, ycol, weightcol; boundary=boundary)
    #
    weights = df[!, weightcol]
    k = kde((xx, yy), weights = weights; boundary = boundary)
    # k = kde((yy, xx), weights=weights)

    xs, ys, img = k.x, k.y, permutedims(k.density)

    # Note on units:
    # https://github.com/JuliaStats/KernelDensity.jl/blob/master/src/bivariate.jl#L58
    # library does `1 / (sum(weights)*(sx*sy)^2)` to make units of the pdf `probabilty/area`,
    # so we mult by the weight sum to make it `weightunit / area`
    # E.G. barrels / meter^2 if we have projected
    img = unnormalize ? img .* sum(df[!, weightcol]) : img
    # If we have projecte4d and are using meters, multiply by 1e6 to make it / km^2
    # If we're using degrees... divide by 10656 to make it / km^2
    # julia> MapImages.latlon_to_dist((ys[1], xs[1]), [1, 0] + [ys[1], xs[1]])  # 111.31709969219834
    # julia> MapImages.latlon_to_dist((ys[1], xs[1]), [0, 1] + [ys[1], xs[1]])  # 96.42150149618836
    # julia> 111*96 10656
    scale = do_project ? 1e6 : 1.0 / 10656
    img .*= scale

    return img, xs, ys
end

# function kde_stack(df, years, xcol=:LatNAD27, ycol=:LatNAD27; kwargs...)
# end

function plot_kde(
    df,
    xcol,
    ycol,
    weightcol;
    unnormalize = true,
    extent = nothing,
    do_project = true,
    cmap = "Reds",
    vmax = nothing,
)
    img, xs, ys = kde_image(
        df,
        xcol,
        ycol,
        weightcol,
        unnormalize = unnormalize,
        extent = extent,
        do_project = do_project,
    )

    fig, ax = plt.subplots()
    axim = ax.imshow(
        img,
        origin = "lower",
        extent = Iterators.flatten(extent),
        cmap = cmap,
        vmax = vmax,
    )
    cbar = fig.colorbar(axim)
    cbar.set_label("Units/sq. km", rotation = 270)
end

function plotstuff()
    imgs = [kde_image(oil_prod, xc, yc, Symbol(y), do_project = false)[1] for y = 2014:2017]
end

### Gridding Functions for summing CSVs into bins ###
#
oob(row, col, out) = row < 1 || row > size(out, 1) || col < 1 || col > size(out, 2)

function bin_vals(
    demrsc,
    df;
    valcol = :sum15_17,
    loncol = :LongNAD27,
    latcol = :LatNAD27,
    smoothing::Union{Nothing,<:Real} = nothing,
)
    out = zeros(size(demrsc))
    for (lon, lat, val) in eachrow(df[:, [loncol, latcol, valcol]])
        row, col = nearest_pixel(demrsc, lat, lon)
        oob(row, col, out) && continue
        out[row, col] += coalesce(val, 0)
    end

    return (!isnothing(smoothing) && smoothing > 0) ?
           imfilter(out, Kernel.gaussian(smoothing)) : out
end

function coarse_bin_vals(
    demrsc,
    df;
    step_km = nothing,
    digits = nothing,
    valcol = :sum15_17,
    loncol = :LongNAD27,
    latcol = :LatNAD27,
    samesize = true,
    sum_vals = true,
)
    if !isnothing(step_km)
        lons, lats = coarse_grid_km(demrsc, step_km)
    elseif !isnothing(digits)
        lons, lats = coarse_grid_deg(demrsc, digits)
    else
        error("Need step_km or digits")
    end

    out = zeros(size(lats, 1), size(lons, 1))
    for (lon, lat, val) in eachrow(df[:, [loncol, latcol, valcol]])
        row = nearest(lats, lat)
        col = nearest(lons, lon)
        oob(row, col, out) && continue
        out[row, col] += (sum_vals ? val : 1)
    end
    return samesize ? resize_nearest(out, size(demrsc)...) : out
end

function resize_nearest(img, nrows, ncols)
    intp = interpolate(img, BSpline(Constant()))
    etpf = extrapolate(intp, Flat())
    return etpf(
        range(0.5, stop = size(img, 1) + 0.5, length = nrows),
        range(0.5, stop = size(img, 2) + 0.5, length = ncols),
    )
end

"""Create a grid in DemRsc area with spacing `step_km`"""
function coarse_grid_km(demrsc, step_km::Union{Number, Nothing} = (10 / 0.62))
    lons, lats = grid(demrsc, sparse = true)
    step = km_to_deg(step_km)
    return lons[1]:step:lons[end], lats[1]:(-step):lats[end]
end

"""Create a grid from a DemRsc on the 10^(-`digits`) grid lines"""
function coarse_grid_deg(demrsc, digits::Union{Int,Nothing} = 2)
    step = 10.0^(-digits)
    lons, lats = grid(demrsc, sparse = true)
    lat1, lat2 = round.(extrema(lats), digits = digits)
    lon1, lon2 = round.(extrema(lons), digits = digits)
    return lon1:step:lon2, lat2:(-step):lat1
end

# Make a new column in `df` with the sum of columns `years`
function sum_years(df, years...; outcol = :yearsum)
    df_out = copy(df)
    yearsyms = [Symbol(y) for y in years]
    df_out[!, outcol] = zeros(size(df_out, 1))

    for i = 1:length(yearsyms)
        df_out[!, outcol] += coalesce.(df[:, yearsyms[i]], 0)
    end
    return df_out
end

funccols(df, f = sum) = [f(col) for col in eachcol(curdf)]
funccols(df, cols, f = sum) = [f(df[!, col]) for col in cols]

function _levels_colors(maxlevel = 100)
    levels = [0, 1, 5, 10, 15, 25, 50, maxlevel]
    colors = broadcast(
        t -> t ./ 256,
        [
            (69, 117, 199, 256),
            (145, 191, 219, 256),
            (224, 243, 248, 256),
            (255, 255, 191, 255),
            (254, 224, 144, 256),
            (252, 141, 89, 256),
            (215, 48, 39, 256),
        ],
    )
    return levels, colors
end

dfmt = dateformat"y-m-d H:M:S.s"
readdf(fname) = coalesce.(CSV.read(fname, dateformat = dfmt), 0)
function plot_wells_per_mi(demrsc; years = [2015, 2016, 2017], outname = nothing)
    oil_prod = CSV.read("oil_production.csv")

    outcol = :yearsum
    df = sum_years(oil_prod, years...; outcol = outcol)

    oil_per_mi =
        coarse_bin_vals(demrsc, df, sum_vals = false, step_km = 1 / 0.62, valcol = outcol)

    levels, colors = _levels_colors(maximum(oil_per_mi))
    fig, ax = plt.subplots()
    axim = ax.contourf(
        oil_per_mi,
        colors = colors,
        levels = levels,
        origin = "image",
        vmax = 45,
        extend = "max",
    )
    fig.colorbar(axim)

    cmap_, norm_ =
        plt.matplotlib.colors.from_levels_and_colors(levels = levels, colors = colors)
    # ax.imshow(oil_per_mi, norm_(oil_per_mi))#, cmap=cmap_)

    if !isnothing(outname)
        save_img_geotiff(outname, oil_per_mi, demrsc, levels, colors)
    end
    return oil_per_mi, fig, ax
end

function save_img_geotiff(outfile::AbstractString, img, demrsc, levels, colors)
    save_img_figure("tmp.png", img, levels, colors)
    kml.create_geotiff(rsc_data = Dict(demrsc), img_filename = "tmp.png", outfile = outfile)
end

function plot_oil_water_bar()
    oil_prod = readdf("oil_production.csv")
    water_inj = readdf("water_injection.csv")

    xs = 2009:2017
    oilyears = funccols(oil_prod, [Symbol(y) for y in xs])
    wateryears = funccols(water_inj, [Symbol(y) for y in xs])

    fig, axes = plt.subplots(1, 2)
    bg_color = "b"
    active_color = "r"
    colors = fill(bg_color, length(xs))
    colors[end] = active_color

    xticks = [2009, 2013, 2017]
    # yticks=[0, 0.5, 1, 1.5, 2, 2.5]

    titles = ["oil", "water"]
    for (title, ax, values) in zip(titles, axes, [oilyears, wateryears])
        ax.bar(xs, values, color = colors)
        ax.set_xticks(xticks)
        # ax.set_yticks(yticks)
        ax.set_ylabel("Million bbl/day")

        ax.set_title(title)
    end
end

function plot_oil_volume_grid(demrsc; outname = nothing)
    oil_prod = readdf("oil_production.csv")
    valcol = Symbol(2017)
    oil_prod17 = @where(oil_prod, cols(valcol) .> 0)

    grid_size_mi = 10
    oil_grid = coarse_bin_vals(
        demrsc,
        oil_prod17,
        sum_vals = true,
        step_km = grid_size_mi / 0.62,
        valcol = valcol,
    )

    oil_grid[oil_grid.==0] .= NaN

    levels = [0, 0.01, 1, 2, 4, 7, 10, 15, 20] .* 1e6
    colors = broadcast(
        t -> t ./ 256,
        [
            (0, 0, 0, 0),
            (69, 117, 199, 200),
            (145, 191, 219, 256),
            (224, 243, 248, 256),
            (255, 255, 191, 255),
            (254, 224, 144, 256),
            (252, 141, 89, 256),
            (215, 48, 39, 256),
        ],
    )
    fig, ax = plt.subplots()
    axim = ax.contourf(
        oil_grid,
        colors = colors,
        levels = levels,
        origin = "image",
        vmax = maximum(levels),
        extend = "max",
    )
    fig.colorbar(axim)

    cmap_, norm_ =
        plt.matplotlib.colors.from_levels_and_colors(levels = levels, colors = colors)

    if !isnothing(outname)
        save_img_geotiff(outname, oil_grid, demrsc, levels, colors)
    end
    return oil_grid, fig, ax
end

function plot_water_volume_grid(demrsc; outname = nothing, grid_size_mi = 5)
    water_inj = readdf("water_injection.csv")
    valcol = Symbol(2017)
    water_inj17 = @where(water_inj, cols(valcol) .> 0)

    water_grid = coarse_bin_vals(
        demrsc,
        water_inj17,
        sum_vals = true,
        step_km = grid_size_mi / 0.62,
        valcol = valcol,
    )

    water_grid[water_grid.==0] .= NaN

    # levels = [0, .01, 1, 2, 4, 7, 10, 15, 20] .* 1e6  # For 10 mi grid
    # levels_adjust = 10 / grid_size_mi / 2
    # levels .*= levels_adjust
    levels = [0, 0.01, 1, 5, 10, 20, 40, 80, 180] .* 1e6
    colors = broadcast(
        t -> t ./ 256,
        [
            (0, 0, 0, 0),
            (69, 117, 199, 200),
            (145, 191, 219, 256),
            (224, 243, 248, 256),
            (255, 255, 191, 255),
            (254, 224, 144, 256),
            (252, 141, 89, 256),
            (215, 48, 39, 256),
        ],
    )
    fig, ax = plt.subplots()
    axim = ax.contourf(
        water_grid,
        colors = colors,
        levels = levels,
        origin = "image",
        vmax = maximum(levels),
        extend = "max",
    )
    fig.colorbar(axim)

    cmap_, norm_ =
        plt.matplotlib.colors.from_levels_and_colors(levels = levels, colors = colors)

    if !isnothing(outname)
        save_img_geotiff(outname, water_grid, demrsc, levels, colors)
    end
    return water_grid, fig, ax
end

macro name(arg)
    x = string(arg)
    quote
        $x
    end
end

depth1(df, depthcol = :BotPerfDepth, lowdepth = 3281, middepth = 6562) =
    @where(df, cols(depthcol) .< lowdepth)
depth2(df, depthcol = :BotPerfDepth, lowdepth = 3281, middepth = 6562) =
    @where(df, cols(depthcol) .> lowdepth, cols(depthcol) .< middepth)
depth3(df, depthcol = :BotPerfDepth, lowdepth = 3281, middepth = 6562) =
    @where(df, cols(depthcol) .> middepth)

function plot_water_bubbles_grouped(demrsc; outname = nothing)
    loncol, latcol, depthcol = :LongNAD27, :LatNAD27, :BotPerfDepth

    delaware =
        sgeom.shape(JSON.parsefile("delaware_shape.geojson")["features"][1]["geometry"])
    midland =
        sgeom.shape(JSON.parsefile("midland_shape.geojson")["features"][1]["geometry"])
    central_basin =
        sgeom.shape(JSON.parsefile("central_basin_shape.geojson")["features"][1]["geometry"])

    water_inj = readdf("water_injection.csv")
    valcol = Symbol(2017)
    water_inj17 = @where(water_inj, cols(valcol) .> 0)
    pts = sgeom.Point.(collect(zip(water_inj17[!, loncol], water_inj17[!, latcol])))

    delaware_pts = water_inj17[delaware.contains.(pts), :]
    central_basin_pts = water_inj17[central_basin.contains.(pts), :]
    midland_pts = water_inj17[midland.contains.(pts), :]

    rows = []
    for p in [delaware_pts, central_basin_pts, midland_pts]
        for d in [depth1, depth2, depth3]
            curdf = d(p)[!, [loncol, latcol, depthcol, valcol]]
            curout = [f(col) for (f, col) in zip([mean, mean, mean, sum], eachcol(curdf))]
            push!(rows, curout)
            println("$d: $curout")
        end
    end
    outdf = DataFrame(permutedims(hcat(rows...)), [:Lon, :Lat, :DepthBot, :Volume2017])
    CSV.write("water_inj_depth_avgs.csv", outdf)
    return outdf
end
# E.g.
# m = rand(5, 5);
# stack = cat([m .+ .2i for i in 1:10]..., dims=3)
function animate_stack(
    stack,
    date_list;
    plot_map = true,
    bigfont = true,
    loncol = :LongNAD27,
    latcol = :LatNAD27,
    xc = loncol,
    yc = latcol,
    sizecol = :mag,
    vm = maximum(stack),
    vmin = nothing,
    vmax = nothing,
    twoway = true,
    outname = "test.gif",
    delay = 500,
    cmap = "seismic_wide_y",
    transparent = false,
    extent = nothing,
    titles = nothing,
)

    if bigfont
        rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
        rcParams["font.size"] = 16
        rcParams["font.weight"] = "bold"
    end

    # (-104.0, -103.001666673056, 30.901666673336, 31.90000000028)
    extent = isnothing(extent) ? MapImages.grid_extent(stack) : extent
    xmin, xmax, ymin, ymax = extent

    fig, ax = plot_map ? plot_states(extent) : plt.subplots()

    # titles = !isnothing(geolist) ? string(geolist) : [string(i) for i in 1:size(stack,3)]
    titles = isnothing(titles) ? string.(date_list) : titles

    for i = 1:size(stack, 3)
        ax.clear()

        fig, ax = plot_map ? plot_states(extent, fig, ax) : (fig, ax)

        img = stack[:, :, i]
        vmin, vmax = _get_vminmax(img, vm, vmin, vmax, twoway = twoway)
        img = transparent ? _cmap_transp(img, plt.cm.get_cmap(cmap), vmin, vmax) : img

        if plot_map
            ax.imshow(
                img,
                vmax = vmax,
                vmin = vmin,
                cmap = cmap,
                extent = extent,
                zorder = 2.5,
                transform = ccrs.PlateCarree(),
            )
        else
            ax.imshow(img, vmax = vmax, vmin = vmin, cmap = cmap, extent = extent)
        end

        ax.set_title(titles[i])
        fname = @sprintf("testgif_%04d", i)
        println("Saving $(titles[i]) as $fname")
        fig.savefig(fname, bbox_inches = "tight")
    end
    run(`magick -delay $(delay/10) -loop 0 testgif_\*.png $outname`)
    rm.(glob("testgif_*.png"))
    return fig, ax

end

# Example:
# xc, yc = :LongNAD27, :LatNAD27
# date_list = collect(2009:2017);
# imgs_oil = cat([kde_image(oil_prod, xc, yc, Symbol(y), do_project=false, extent=extent)[1] for y in date_list]..., dims=3);
# animate_stack(imgs, date_list, delay=900, vmax=5e4, vmin=0, cmap="Reds", bigfont=false, extent=extent, transparent=true, titles=["Oil Production density in $d" for d in date_list])

# TODO: do this instead to CMAP, and not to the image, so that we can still use a colorbar (that isn't 0 to 1)
# https://stackoverflow.com/questions/25408393/getting-individual-colors-from-a-color-map-in-matplotlib
Normalize = pyimport("matplotlib.colors").Normalize
function _cmap_transp(img, cmap, vmin = 0, vmax = maximum(img))
    # Any value whose absolute value is > .0001 will have zero transparency
    # alphas = Normalize(0, .3, clip=true)(abs.(img))
    alphas = Normalize(0, vmax, clip = true)(abs.(img))
    alphas = clamp.(alphas, 0.2, 1)  # alpha value clipped at the bottom at .2

    colors = Normalize(vmin, vmax, clip = true)(img)
    colors = cmap(colors)
    colors[:, :, end] = alphas
    return colors
end

function animate_imgs_pts(
    df,
    date_list,
    stack::Union{MapImage,Nothing} = nothing;
    loncol = :lon,
    latcol = :lat,
    sizecol = :mag,
    datecol = :datetime,
    size_scale = 4,
    outname = "test.gif",
    alpha = 0.4,
    vm = 6,
    c = "r",
    delay = 200,
    cmap = "seismic_wide_y",
    extent = nothing,
    titles = nothing,
    bigfont = true,
)
    # (-104.0, -103.001666673056, 30.901666673336, 31.90000000028)
    extent = isnothing(stack) ? extent : MapImages.grid_extent(stack)
    xmin, xmax, ymin, ymax = extent

    # fig, ax = plt.subplots()
    fig, ax = plot_states(extent)
    # titles = !isnothing(geolist) ? string(geolist) : [string(i) for i in 1:size(stack,3)]
    titles = isnothing(titles) ? string(date_list) : titles

    if bigfont
        rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
        rcParams["font.size"] = 16
        rcParams["font.weight"] = "bold"
    end

    prevdate = Date(1, 1, 1)
    # for i=1:size(stack, 3)
    for i = 1:length(date_list)
        # for i=1:4
        curdate = date_list[i]
        df_sub = df[(df[:, datecol].<curdate).&(df[:, datecol].>prevdate), :]
        prevdate = curdate

        (lons, lats, sizes) = (df_sub[:, loncol], df_sub[:, latcol], df_sub[:, sizecol])
        sizes .*= size_scale

        ax.clear()
        fig, ax = plot_states(extent, fig, ax)

        ax.scatter(
            lons,
            lats,
            sizes,
            c = c,
            alpha = alpha,
            edgecolor = "none",
            zorder = 2.5,
            transform = ccrs.PlateCarree(),
        )
        !isnothing(stack) &&
        ax.imshow(stack[:, :, i], vmax = vm, vmin = -vm, cmap = cmap, extent = extent)
        ax.set_title(titles[i])
        fname = @sprintf("testgif_%04d", i)
        println("Saving $(titles[i]) as $fname")
        fig.savefig(fname, bbox_inches = "tight")
    end
    run(`magick -delay $(delay/10) -loop 0 testgif_\*.png $outname`)
    rm.(glob("testgif_*.png"))
    return fig, ax
end
# E.G.
# eqs_txar = readdf("TXAR_events.csv"); eqs_txar_big = eqs_txar[(eqs_txar[!, :mag ] .> 0.7) .& (eqs_txar[!, :QP] .< 3), :];
# animate_imgs_pts(eqs_txar_big, Date(2001,1,1):Dates.Year(1):Date(2018,1,1), extent=(-105., -100., 29., 33.), delay=650, size_scale=8, titles=["TXAR earthquakes in $(string(year(d)))" for d in date_list])
#
# animate_imgs_pts(eqs_txar_big, Date(2001,1,1):Dates.Year(1):Date(2018,1,1), extent=MapImages.grid_extent(demrsc), delay=550, size_scale=8, titles=[string(year(d)) for d in date_list])
# animate_imgs_pts(eqs_txar, Date(2001,1,1):Dates.Year(1):Date(2017,1,1), extent=MapImages.grid_extent(demrsc), delay=350)
# animate_imgs_pts(eqs_txar, geolist, stack78)

# plt.figure(); plt.contourf(oil_per_mi, colors=colors, levels=[0, 1, 5, 10, 15, 25, 50, maximum(oil_per_mi)], origin="image", vmax=45, extend="max"); plt.colorbar()

function plot_states(
    extent,
    fig = nothing,
    ax = nothing;
    add_counties = true,
    add_basins = true,
)
    fig = isnothing(fig) ? plt.figure() : fig
    ax = isnothing(ax) ? fig.add_axes([0, 0, 1, 1], projection = ccrs.LambertConformal()) :
        ax

    # extent = isnothing(stack) ? extent : MapImages.grid_extent(stack)
    ax.set_extent(extent, ccrs.Geodetic())

    # to get the effect of having just the states without a map "background"
    # turn off the outline and background patches
    ax.background_patch.set_visible(false)
    ax.outline_patch.set_visible(false)

    # shapename = "admin_1_states_provinces_lakes_shp"
    # states_shp = shpreader.natural_earth(resolution="110m", category="cultural", name=shapename)
    # states_shp = shpreader.natural_earth(resolution="50m", category="cultural", name=shapename)
    # states_shp = shpreader.natural_earth(resolution="10m", category="cultural", name=shapename)

    # records = shpreader.Reader(states_shp).records()
    # states = shpreader.Reader(states_shp).geometries()
    # ff = shpreader.Reader("/Users/scott/.local/share/cartopy/shapefiles/natural_earth/cultural/ne_10m_admin_1_states_provinces_lakes_shp.shp")
    ff =
        shpreader.Reader("/Users/scott/Downloads/cb_2018_us_state_5m/cb_2018_us_state_5m.shp")
    records = ff.records()
    states = ff.geometries()
    # records = ff.records()
    for (record, state) in zip(records, states)
        # for state in states
        # pick a default color for the land with a black outline,
        # this will change if the storm intersects with our track
        facecolor =
            record.attributes["NAME"] == "Texas" ? [0.9375, 0.9375, 0.859375] : [1, 1, 1, 0]
        # facecolor = [0.9375, 0.9375, 0.859375]
        edgecolor = "black"

        # if state.intersects(track)
        #     facecolor = "red"
        # elseif state.intersects(track_buffer):
        #     facecolor = "#FF7E00"
        # end

        ax.add_geometries(
            [state],
            ccrs.PlateCarree(),
            facecolor = facecolor,
            edgecolor = edgecolor,
        )
    end
    pecos = (-103.49283, 31.4243)
    ax.plot(pecos..., "bo", ms = 10, transform = ccrs.PlateCarree(), zorder = 3)
    plt.text(
        ((-.09, 0.09) .+ pecos)...,
        "Pecos",
        fontdict = Dict("size" => 14, "weight" => "bold"),
        horizontalalignment = "right",
        transform = ccrs.Geodetic(),
    )

    if add_counties
        reader = shpreader.Reader("./countyl010g_shp_nt00964/countyl010g.shp")
        counties = collect(reader.geometries())
        COUNTIES = cfeature.ShapelyFeature(counties, ccrs.PlateCarree())
        ax.add_feature(COUNTIES, facecolor = "none", edgecolor = "gray")
    end

    if add_basins
        reader =
            shpreader.Reader("./repermianbasinshapefile/Permian_Boundary_2018_TORA_Dutton.shp")
        basins = collect(reader.geometries())[[3, 4, 6]]
        # ax.add_geometries(basins, ccrs.epsg(3665), facecolor="gray", edgecolor="gray")
        BASINS = cfeature.ShapelyFeature(basins, ccrs.epsg(3665))
        ax.add_feature(BASINS, facecolor = "none", edgecolor = "black", linewidth = 3)
    end

    return fig, ax
end

# oil_totals = CSV.read("PermianTotalByYear.csv"); oil_totals = oil_totals[1:end-1, :]
function animate_bar(
    xs,
    ys;
    outname = "test.gif",
    delay = 200,
    bg_color = "b",
    active_color = "r",
    titles = string.(xs),
    bigfont = true,
    xticks = [2000, 2009, 2018],
    yticks = [0, 0.5, 1, 1.5, 2, 2.5],
)
    # years, values = oil_totals[:, :year], oil_totals[:, :oil_bbl]
    fig, ax = plt.subplots()
    if bigfont
        rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
        rcParams["font.size"] = 16
        rcParams["font.weight"] = "bold"
    end

    for i = 1:length(xs)
        ax.clear()

        colors = fill(bg_color, length(xs))
        colors[i] = active_color

        ax.bar(xs, values, color = colors)

        ax.set_xticks(xticks)
        ax.set_yticks(yticks)
        ax.set_ylabel("Million bbl/day")

        ax.set_title(titles[i])
        fname = @sprintf("testgif_%04d", i)
        println("Saving $(titles[i]) as $fname")
        fig.savefig(fname, bbox_inches = "tight")
    end
    run(`magick -delay $(delay/10) -loop 0 testgif_\*.png $outname`)
    rm.(glob("testgif_*.png"))
    return fig, ax
end

# function make_frame(i)
#     plt.imshow(stack[:,:,i+1], vmax=vm, vmin=-vm)
# end
# myanim = anim.FuncAnimation(fig, make_frame, frames=size(stack,3), interval=20)
# myanim[:save]("test2.mp4", bitrate=-1, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])

# TODO Data cleanup
# water_inj[!, :BotPerfDepth] .== 0
# water_inj17 = @transform(water_inj17, depthdiff = :BotPerfDepth .- :TopPerfDepth);
# water_inj17[water_inj17.depthdiff .< 0, :]

ext: jl


text: @test problem1(1000) == 233168

max_stars_repo_name: halleysfifthinc/ProjectEuler

max_issues_repo_name: halleysfifthinc/ProjectEuler

max_forks_repo_name: halleysfifthinc/ProjectEuler


text:
wall = Lar.svg2lar("./test/svg/assembly/wall.svg", flag=false);
openings = Lar.svg2lar("./test/svg/assembly/svg/openings.svg", flag=false);
rectangle = Lar.svg2lar("./test/svg/assembly/rectangle.svg", flag=false);
box = Lar.svg2lar("./test/svg/assembly/box.svg", flag=false);

assembly = Lar.Struct([wall, openings, rectangle, box]);
V,EV = Lar.struct2lar(assembly);
if todisplay  GL.VIEW([ GL.GLGrid(V,EV, GL.COLORS[1],1), GL.GLFrame2 ]);  end

model = CAGD.Model(V);
CAGD.addModelCells!(model, 1, EV);

A_model = CAGD.planar_arrangement(model);

A_model.G
Matrix(A_model.T[1])
Matrix(A_model.T[2])

# Display Models

if todisplay
    VV = [[k] for k=1:size(A_model.G, 2)];
    triangulated_faces = CAGD.triangulate2D(A_model);
    FVs = convert(Array{Lar.Cells}, triangulated_faces);

    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.0, 1.0, 1.0, 99, 1));
    GL.VIEW(GL.GLExplode(A_model.G, FVs, 4.0, 4.0, 4.0, 99, 1));
end

innerpoints = Lar.internalpoints2d(A_model.G', A_model.T[1], A_model.T[2][1:end, :])
# associate internal points to 2-cells
listOfModels = Lar.evalStruct(assembly)
inputfacenumbers = [length(listOfModels[k][2]) for k=1:length(listOfModels)];
# test input data for containment of reference points
boolmatrix = sparse(BitArray(undef, length(innerpoints), length(listOfModels)));
containmenttest = Lar.testinternalpoint2d(listOfModels);
for (k,point) in enumerate(innerpoints) # k runs on columns
    cells = containmenttest(point); # contents of columns
    #println(k," ",faces)
    for l in cells
        boolmatrix[k,l] = 1;
    end
end

boolmatrix = sparse(Bool[
    true  false false true
    true  false true  true
    false false false true
    true  true  false true
    false true  false true
    false false false true
    true  true  false true
    false false false true
    false true  false false
    false false true  true
    true  false false true
    false true  false true
    true  false false true
    true  false false true
    false true  false true
    false false  true  true
    false false true  false
]);

Matrix(boolmatrix)

wall = boolmatrix[:, 1]
Lshape = boolmatrix[:, 2]
rect = boolmatrix[:, 3]
box = boolmatrix[:, 4]

GL.VIEW(GL.GLExplode(A_model.G, FVs[wall.nzind],   1.0, 1.0, 1.0, 99, 1));
GL.VIEW(GL.GLExplode(A_model.G, FVs[Lshape.nzind], 1.0, 1.0, 1.0, 99, 1));
GL.VIEW(GL.GLExplode(A_model.G, FVs[rect.nzind],   1.0, 1.0, 1.0, 99, 1));
GL.VIEW(GL.GLExplode(A_model.G, FVs[box.nzind],    1.0, 1.0, 1.0, 99, 1));

shells = Array{SparseVector{Bool,Int64},1}()
push!(shells, .|(wall, Lshape, rect));
push!(shells, .&(wall, .!Lshape, .!rect));
push!(shells, .&(box, .!shells[2]))

if todisplay
    for i = 1 : length(shells)
    GL.VIEW(GL.GLExplode(A_model.G, FVs[shells[i].nzind], 1.0, 1.0, 1.0, 1, 1));
    GL.VIEW(GL.GLExplode(A_model.G, FVs[shells[i].nzind], 1.0, 1.0, 1.0, 99, 1));
    end
end


text:
__precompile__(false)

module WindowManager

using GLFW

using CodeGeneration
using SystemManager

###############################################################################

CodeGeneration.module_listeners(@__MODULE__)

###############################################################################

DRIVER_CHOICE = Symbol(GLFW)
#macro __DRIVER__() :($DRIVER_CHOICE) end

macro DRIVER() Meta.quot(SystemManager.@DRIVER DRIVER_CHOICE) end
macro DRIVER_INIT() Meta.quot(SystemManager.@DRIVER_INIT DRIVER_CHOICE) end
macro DRIVER_CHECK() Meta.quot(SystemManager.@DRIVER_CHECK DRIVER_CHOICE) end
#macro DRIVER_ENABLED() Meta.quot(SystemManager.@DRIVER_ENABLED DRIVER_CHOICE) end

###############################################################################

@DRIVER

function driver_init()
  GLFW.Terminate()  # remove previous Window
  GLFW.Init() # OS X-specific GLFW hints to initialize the correct version of OpenGL
  @DRIVER_INIT
end

const WindowRef = GLFW.Window
export WindowRef

mutable struct Window
  id::Symbol
  ref::Union{Nothing, WindowRef}
  pos::Tuple{Int32,Int32}
  size::Tuple{UInt32,UInt32}
  title::String
  focus::Bool
  fullScreen::Bool

  listenList::Dict{Symbol, Tuple{Symbol, Function}}
  events::Array{Array{Any,1},1}

  lastPos::AbstractVector
  lastCursorPos::AbstractVector
  lastSize::Tuple{Number,Number}

  function Window(id::Union{Nothing,Symbol}=nothing)
    this=new(generateID(id), #random id
    nothing,(0,0),(1,1),"",true,false,
    Dict(),[],[],[],(1,1))
    WindowList[this.id]=this
    this
  end
end
export Window

WindowList = Dict{Symbol,Window}()
WindowRefList = Dict{Symbol,Window}()

lastKey=Number(0)
KEYS=Dict{Number,Number}()

function cleanUp()
  global WindowList, WindowRefList
  callOnListeners(:cleanUp)
  WindowList = typeof(WindowList)()
  WindowRefList = typeof(WindowRefList)()
  GLFW.Terminate()
end

getWindow(ref::WindowRef)::Window = WindowRefList[Symbol(ref)]
export getWindow

function setWindow(this::Window, ref::WindowRef)
  this.ref = ref
  WindowRefList[Symbol(ref)] = this
end

function remove(this::Window)
  global WindowRefList, WindowList
  filter!((x)->x[1]!=this.id,WindowList)
  filter!((x)->x[1]!=Symbol(this.ref),WindowRefList)
end

# Create a window and its OpenGL context
"""
TODO
"""
function open(this::Window, options::Array{Pair{Symbol,T},1}=[]) where T
    if !SystemManager.driver_check(DRIVER_CHOICE;soft=true) driver_init() end

    for (id,val) in options
      GLFW.WindowHint(eval(:(GLFW.$id)),val)
    end

    (width,height) = this.size

    setWindow(this, GLFW.CreateWindow(width, height, this.title))

		pos = GLFW.GetWindowPos(this.ref)
    this.pos = (pos.x,pos.y)

  	GLFW.MakeContextCurrent(this.ref)
    GLFW.SwapInterval(0)
    #showError(true)
    GLFW.ShowWindow(this.ref)

    # OpenGL Version
    #GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR,4)
    #GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR,6)

    #GLFW.SetWindowSize(window, w, h) # Seems to be necessary to guarantee that window > 0

    # Window Callbacks
    GLFW.SetWindowIconifyCallback(this.ref, OnWindowIconify)
    GLFW.SetWindowSizeCallback(this.ref, OnWindowResize)
    GLFW.SetWindowCloseCallback(this.ref, OnWindowClose)
    GLFW.SetWindowFocusCallback(this.ref, OnWindowFocus)
    GLFW.SetWindowPosCallback(this.ref, OnWindowPos)
    GLFW.SetWindowRefreshCallback(this.ref, OnWindowRefresh)
    GLFW.SetFramebufferSizeCallback(this.ref, OnFramebufferResize)
    GLFW.SetCursorPosCallback(this.ref, OnCursorPos)
    GLFW.SetCursorEnterCallback(this.ref, OnCursorEnter)
    GLFW.SetDropCallback(this.ref, OnDroppedFiles)

	   # Input Callbacks
    GLFW.SetKeyCallback(this.ref, OnKey)
    GLFW.SetMouseButtonCallback(this.ref, OnMouseKey)
    GLFW.SetScrollCallback(this.ref, OnScroll)
    GLFW.SetCharCallback(this.ref, OnUnicodeChar)
    GLFW.SetCharModsCallback(this.ref, OnCharMods)
end

function resize(this::Window, size::Tuple{Integer,Integer})
  (width,height) = this.size = size
  GLFW.SetWindowSize(this.ref, width, height) # Seems to be necessary to guarantee that window > 0
end

################################################################################
# Window Events

"""
TODO
"""
getEventNames() = [
  :OnWindowResize=>[Window,Tuple{Number,Number}],
  :OnWindowIconify=>[Window,Number],
  :OnWindowClose=>[Window],
  :OnWindowFocus=>[Window,Number],
  :OnWindowRefresh=>[Window],
  :OnFramebufferResize=>[Window,Tuple{Number,Number}],
  :OnCursorEnter=>[Window,Number],
  :OnDroppedFiles=>[Window,AbstractArray],
  :OnWindowPos=>[Window,Tuple{Number,Number}],
  :OnWindowMove=>[Window,Tuple{Number,Number}],
  :OnMousePos=>[Window,Tuple{Number,Number}],
  :OnMouseMove=>[Window,Tuple{Number,Number}],
  :OnScroll=>[Window,:x=>Number,:y=>Number],
  :OnCharMods=>[Window,:code=>Char,:mods=>Number],
  :OnKey=>[Window,:key=>Number,:scancode=>Number,:action=>Number,:mods=>Number],
  :OnMouseKey=>[Window,:key=>Number,:action=>Number,:mods=>Number],
]

"""
TODO
"""
function OnWindowResize(this::WindowRef, width::Number, height::Number)
	window = getWindow(this)
	window.size = (width,height)
	OnEvent(window, :OnWindowResize, [width,height])
end

"""
TODO
"""
OnWindowIconify(this::WindowRef, iconified::Number) = OnEvent(getWindow(this), :OnWindowIconify, iconified)

"""
TODO
"""
OnWindowClose(this::WindowRef) = OnEvent(getWindow(this), :OnWindowClose)

"""
TODO
"""
function OnWindowFocus(this::WindowRef, focused::Number)
  window = getWindow(this)
  window.focus = focused
  OnEvent(window, :OnWindowFocus, focused)
end

"""
TODO
"""
OnWindowRefresh(this::WindowRef) = OnEvent(getWindow(this), :OnWindowRefresh)

"""
TODO
"""
OnFramebufferResize(this::WindowRef, width::Number, height::Number) =  OnEvent(getWindow(this), :OnFramebufferResize, [width,height])

"""
TODO
"""
OnCursorEnter(this::WindowRef, entered::Number) = OnEvent(getWindow(this), :OnCursorEnter, entered)

"""
TODO
"""
OnDroppedFiles(this::WindowRef, files::AbstractArray) = OnEvent(getWindow(this), :OnDroppedFiles, files)

"""
TODO
"""
function OnWindowPos(this::WindowRef, x::Number, y::Number)
  window = getWindow(this)
  t = [-x,y]
  OnEvent(window, :OnWindowPos, (t[1],t[2])) #
  l = window.lastPos; if l == [] l = t end
  n = l - t
  window.lastPos = t
  OnEvent(window, :OnWindowMove, (n[1],n[2])) #
end

"""
TODO
"""
function OnCursorPos(this::WindowRef, x::Number, y::Number)
  window = getWindow(this)
  t = [x,y]
  OnEvent(window, :OnMousePos, (t[1],t[2]))
  l = window.lastCursorPos; if l == [] l = t end
  n = l - t
  window.lastCursorPos = t
  OnEvent(window, :OnMouseMove, (n[1],n[2]))
end

# Input Events
"""
TODO
"""
OnScroll(this::WindowRef, x::Number, y::Number) = OnEvent(getWindow(this), :OnScroll, [x,y])

"""
TODO
"""
OnCharMods(this::WindowRef, code::Char, mods::Number) = OnEvent(getWindow(this), :OnCharMods, code, mods)

"""
TODO
"""
function OnKey(this::WindowRef, key::GLFW.Key, scancode::Number, action::GLFW.Action, mods::Number)
	window = getWindow(this)
  key, scancode, action, mods, unicode = UInt32(key), UInt32(scancode), UInt32(action), UInt32(mods), UInt32(0)
	global lastKey = key
	push!(window.events, [ :OnKey, key, scancode, action, mods, unicode ])
end

"""
TODO
"""
OnMouseKey(this::WindowRef, key::GLFW.MouseButton, action::GLFW.Action, mods::Number) = OnEvent(getWindow(this), :OnMouseKey, UInt32(key), UInt32(action), UInt32(mods))

"""
TODO
"""
function OnUnicodeChar(this::WindowRef, unicode::Char)
	window = getWindow(this)
	if lastKey > 0 KEYS[lastKey]=Number(unicode) end
	global lastKey=0
end

################################################################################

"""
TODO
"""
function OnUpdateEvents()
  GLFW.PollEvents()  # Poll for and process events
	for (_,window) in WindowRefList
		for e in window.events
			if length(e) > 5 && e[1] == :OnKey
				key=e[2]
				if haskey(KEYS, key) e[6]=KEYS[key] end
			end
			OnEvent(window, e...)
		end
		if length(window.events) > 0	window.events	= [] end
	end
end

"""
TODO
"""
function OnEvent(this::Window, eventName::Symbol, args...)
 #a = values(d) # convert Dict Values to Array
  for (_,listener) in this.listenList
		if listener[1] == eventName listener[2](this, args...) end
    #s = listener.storage
		#if haskey(s.events, eventName) s.events[eventName](args...) end
  end
end

################################################################################

"""
TODO
"""
addListener(this::Window, eventName::Symbol, callback::Function; id::Union{Nothing,Symbol}=nothing) =
	this.listenList[generateID(id)] = (eventName, callback)

"""
TODO
"""
function close(this::Window)
  remove(this)
	if length(WindowRefList) == 0 GLFW.Terminate() end
end

"""
TODO
"""
swap(this::Window) = GLFW.SwapBuffers(this.ref)   # Swap front and back buffers

#function GLFW_SetWindowMonitor(window::GLFW.Window, monitor::GLFW.Monitor, xpos, ypos, width, height, refreshRate)
#  ccall((:glfwSetWindowMonitor, GLFW.lib), Nothing, (GLFW.Window, GLFW.Monitor, Cint, Cint, Cint, Cint, Cint), window, monitor, xpos, ypos, width, height, refreshRate)
  # GLFW.SetWindowMonitor
#end

"""
TODO
"""
function fullscreen(this::Window, full::Bool)
	monitor = GLFW.GetPrimaryMonitor()
	mode = GLFW.GetVideoMode(monitor)
	pos = GLFW.GetMonitorPos(monitor)

	if !full
		pos = this.pos
		(width,height) = this.size = this.lastSize
		mode = GLFW.VidMode(width, height, mode.redbits, mode.greenbits, mode.bluebits, mode.refreshrate)
		monitor = GLFW.Monitor(C_NULL)
	else
		this.lastSize = this.size
		this.pos = GLFW.GetWindowPos(this.ref)
		this.size = GLFW.GetWindowSize(this.ref)
	end

	this.fullScreen = full
	GLFW.SetWindowMonitor(this.ref, monitor, pos[1], pos[2], mode.width, mode.height, mode.refreshrate)
end

"""
TODO
"""
cursor(this::Window, mode::Symbol) = GLFW.SetInputMode(this.ref, GLFW.CURSOR, eval(:(GLFW.$mode)))
#GLFW.STICKY_KEYS or GLFW.STICKY_MOUSE_BUTTONS

"""
TODO
"""
function title(this::Window, title::String)
  this.title = title
  GLFW.SetWindowTitle(this.ref, title)
end

"""
TODO
"""
function showError(debugging::Bool)
  @static if Sys.isapple() return end # not working on apple
  GLFW.WindowHint(GLFW.OPENGL_DEBUG_CONTEXT, convert(Cint, debugging))
end

"""
TODO
"""
isOpen(this::Window) = !GLFW.WindowShouldClose(this.ref)

"""
TODO
"""
hasFocus(this::Window) = this.focus

"""
TODO
"""
function update(this::Window)
  if hasFocus(this) swap(this) end
	OnUpdateEvents()
end

"""
Loop until the user closes the window
"""
function loop(this::Window, repeat::Function)
	while isOpen(this)
    if !hasFocus(this) sleep(0.1)
    else repeat(this)
    end
    update(this) # update window
	end
	close(this)
end

end #WindowManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
@testset "SingleVector" begin
    a = SingleVector(3)
    @test length(a) == 1
    @test a[1] == 3
    @test a[] == 3
    @test_throws BoundsError a[2]
end


text:
## Initialization Filters

using LinearAlgebra, FFTW

generatefilters(type,dims;form=:list) = _generatefilters(Val(type),dims,Val(form))

# DCT
function _generatefilters(::Val{:DCT},dims,::Val{:list})
    H = _generatefilters(Val(:DCT),dims,Val(:matrix))
    return Array.(_filterlist(H,dims))
end
function _generatefilters(::Val{:DCT},dims,::Val{:matrix})
    @assert length(dims) == 2           "Only 2D DCT is implemented"      # todo
    @assert all(i->i==first(dims),dims) "Only square DCT is implemented"  # todo

    temp = dct(Matrix(I,dims),1)
    return kron(temp,temp)' / sqrt(prod(dims))
end


text:
using StatisticalRethinking, CmdStan, LinearAlgebra
#gr(size=(600,600));

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

df2 = filter(row -> row[:age] >= 18, df);
first(df2, 5)

!isfile(joinpath(ProjDir, "m4.1s.jls")) && include(joinpath(ProjDir, "m4.1s.jl"))

chn = deserialize(joinpath(ProjDir, "m4.1s.jls"))

MCMCChains.describe(chn)

density(chn, lab="All heights", xlab="height [cm]", ylab="density")

mu_sigma = hcat(chn.value[:, 2, 1], chn.value[:,1, 1])
LinearAlgebra.diag(cov(mu_sigma))

cor(mu_sigma)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
#------------------------------
# functions to simulate data
#------------------------------

function generate_xs(n, p, true_u0, sol_group1, sol_group2; t_start=1.5, t_end=10, maxntps = 10, dt=0.1, σ_var=0.1, σ_ind=0.5)

    # (1) generate artifical group labels
    groups = zeros(n)
    groups[randperm(n)[1:Int(floor(n/2))]] .= 1 # for groups of the same size
    group1 = findall(x -> x==1, groups)
    group2 = findall(x -> x==0, groups)

    # (2) generate artificial time stamps
    ntps = rand(1:maxntps, n)
    tvals = [sort(rand(t_start:dt:t_end,ntps[i])) for i in 1:n]

    # (3) obtain true values as solutions of the ODEs at the initial time point and the drawn second time point
    # check for equal number of variables:
    if p%2 != 0
        error("Please select an even number of variables")
    end
    # true starting point
    z_t0_p1 = true_u0[1] # for variables 1-(p/2)
    z_t0_p2 = true_u0[2] # for variables (p/2+1)-p
    z_t0 = repeat([z_t0_p1, z_t0_p2], inner=Int(p/2))

    # now use ODE solutions to obtain true temporal development value
    # for all individuals in both variables u1 and u2
    z_later_ts = collect((i ∈ group1) ? (Array(sol_group1)[:,Int.(round.(tvals[i].*(1 ./dt)).+1)]) : (Array(sol_group2)[:,Int.(round.(tvals[i].*(1 ./dt)).+1)]) for i in 1:n)

    # (4) sample variable- specific and individual-specific errors at both time points
    # variable specific random effect (general difficulty measuring that specific variable)
    us = rand(Normal(0,σ_var),p)

    xs = []
    for i in 1:n
        # make time series structure, should have shape (p x ntps[i])
        cur_timeseries = zeros(p, ntps[i]+1)
        for j in 1:p
            cur_timeseries[j,1] = z_t0[j] + us[j] + randn() .* σ_ind
            for tp in 1:ntps[i]
                if j <= Int(p/2)
                    cur_timeseries[j,tp+1] = z_later_ts[i][1,tp] + us[j] + randn() .* σ_ind
                else
                    cur_timeseries[j,tp+1] = z_later_ts[i][2,tp] + us[j] + randn() .* σ_ind
                end
            end
        end
        push!(xs, cur_timeseries)
    end

    return xs, tvals, group1, group2
end

function generate_baseline(n, q, q_info, group1; σ_info=1, σ_noise=1)
    zs = fill(1,(n,1))
    zs[group1].= -1
    means = fill(0,n,q)
    means[:,1:q_info] .= zs

    vars=fill(σ_noise,q)
    vars[1:q_info] .= σ_info
    x_params = [cat([rand(Normal(means[i,j],vars[j])) for j in 1:q]..., dims=1) for i in 1:n]
    return x_params
end

function generate_baseline(n, q, q_info, group1, true_odeparams_group1, true_odeparams_group2; σ_info=0.1, σ_noise=0.1)
    signs = fill(1,(n,1))
    signs[group1] .= -1
    z1s = zeros(n,1)
    z1s[signs .== -1] .= true_odeparams_group1[1]
    z1s[signs .== 1] .= true_odeparams_group2[1]
    z2s = zeros(n,1)
    z2s[signs .== -1] .= true_odeparams_group1[4]
    z2s[signs .== 1] .= true_odeparams_group2[4]

    means = zeros(n,q)
    means[:,1:Int(floor(q_info/2))] .= z1s
    means[:,Int(floor(q_info/2))+1:q_info] .=z2s

    vars=fill(σ_noise,q)
    vars[1:q_info] .= σ_info
    x_params = [cat([rand(Normal(means[i,j],vars[j])) for j in 1:q]..., dims=1) for i in 1:n]
    return x_params
end

function generate_baseline_4params(n, q, q_info, group1, true_odeparams_group1, true_odeparams_group2; σ_info=0.1, σ_noise=0.1)
    signs = fill(1,n)
    signs[group1] .= -1
    zs = cat(repeat([true_odeparams_group2], n)..., dims=2)
    zs[:,findall(signs .==-1)] .= true_odeparams_group1
    zs = zs'

    means = zeros(n,q)
    means[:,1:Int(floor(q_info/4))] .= zs[:,1]
    means[:,Int(floor(q_info/4))+1:Int(floor(q_info/2))] .= zs[:,2]
    means[:,Int(floor(q_info/2))+1:Int(floor(q_info*3/4))] .= zs[:,3]
    means[:,Int(floor(q_info*3/4))+1:q_info] .= zs[:,4]

    vars=fill(σ_noise, q)
    vars[1:q_info] .= σ_info
    x_params = [[rand(Normal(means[i,j],vars[j])) for j in 1:q] for i in 1:n]
    return x_params
end


text:
using Catalyst, ReactionNetworkImporters, SparseArrays

# version giving parameters and rates
rs = @reaction_network begin
    k1, 2S1 --> S2
    k2, S2 --> 2S1
    k3, S1 + S2 --> S3
    k4, S3 --> S1 + S2
    k5, 3S3 --> 3S1
end k1 k2 k3 k4 k5

@parameters k1 k2 k3 k4 k5
pars = [k1,k2,k3,k4,k5]
substoich =[2 0 0;
            0 1 0;
            1 1 0;
            0 0 1;
            0 0 3]'
prodstoich = [0 1 0;
              2 0 0;
              0 0 1;
              1 1 0;
              3 0 0]'
prn = loadrxnetwork(MatrixNetwork(), pars, substoich, prodstoich; params=pars)
@test rs == prn.rn

# sparse version
prn = loadrxnetwork(MatrixNetwork(), pars, sparse(substoich), sparse(prodstoich); params=pars)
@test rs == prn.rn

# version with hard coded rates (no parameter symbols)
rs = @reaction_network begin
    1., 2S1 --> S2
    2., S2 --> 2S1
    3., S1 + S2 --> S3
    4., S3 --> S1 + S2
    5., 3S3 --> 3S1
end
prn = loadrxnetwork(MatrixNetwork(), convert.(Float64,1:5), substoich, prodstoich)
@test rs == prn.rn

# sparse version
prn = loadrxnetwork(MatrixNetwork(), convert.(Float64,1:5), sparse(substoich), sparse(prodstoich))
@test rs == prn.rn

# version with species names, rate functions and symbolic parameters
rs = @reaction_network begin
    k1*A, 2A --> B
    k2, B --> 2A
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
end k1 k2 k3 k4 k5
@parameters t
@variables A(t) B(t) C(t)
species = [A,B,C]
rates = [k1*A,k2,k3,k4,k5]
prn = loadrxnetwork(MatrixNetwork(), rates, substoich, prodstoich; species=species, params=pars)
@test rs == prn.rn

# sparse version
prn = loadrxnetwork(MatrixNetwork(), rates, sparse(substoich), sparse(prodstoich); species=species, params=pars)
@test rs == prn.rn


text:
#A partir de PDFs con tablas de la SSA, devuelve un CSV.

#Paquetería para leer pdfs y manipular fechas
using PDFIO
using Dates
using DelimitedFiles
using ArgParse

#Devuelve el texto del pdf como un string
function texto_pdf(archivo)

  #Abre el documento.
  documento = pdDocOpen(archivo)
  #Genera objetos que representan las páginas.
  pages = pdDocGetPageRange(documento, 1:pdDocGetPageCount(documento))
  #Creamos un buffer para compilar los datos.
  datos = IOBuffer()

  #Escribimos el texto al buffer.
  #Entre cada página forzamos una línea nueva.
  for page in pages
    pdPageExtractText(datos, page)
    write(datos, " \n")
  end

  #Cerramos el documento.
  pdDocClose(documento)

  #Devolvemos el texto en un string.
  return String(take!(datos))
end

#Separa el string por líneas y elimina aquellas que no corresponden a casos.
function eliminar_nocasos(string)

  #Separa por línea
  rows = split(string, "\n")
  #Revisa que el inicio de la línea sea de la forma {espacios}Número y descarta la línea si no.
  filter!(row -> occursin(r"^ +\d", row), rows)
  #Transforma los substrings en strings:
  rows = String.(rows)

  return rows
end

function procesa_fecha(string)

  longitud_campo = length(string)

  if longitud_campo == 10

    try
      return Dates.format(Date(string, "dd/mm/yyyy"), "yyyy-mm-dd")
    catch e
      return ""
    end

  elseif longitud_campo == 5

    #Son fechas representadas en un formato de Excel: con el conteo de días iniciando en 1900-01-01. (Día 1)
    #https://www.covid19in.mx/docs/datos/tablas-casos/normalizacion/fecha/
    #Sin embargo, por retrocompatibilidad, se considera erróneamente que 1900 es un año bisiesto,
    #por lo que para convertir las fechas a partir del 1900-02-28 se requiere restar dos al número mostrado,
    #una unidad para tener un intervalo de tiempo respecto de 1900-01-01 y otra unidad para corregir por el día
    #extra considerado.
    #Gracias a Juan Claudio Toledo Roy por estas correciones.

    Δt = parse(Int, string) - 2  #Nuestras fechas son posteriores a 1900-02-28.

    fecha = Date("1900-01-01") + Day(Δt) |> Base.string

    return fecha

  else

    return ""
  end
end

function procesa_fila(string, index_fechas)
  out = replace(string, r"^\s+" => "")  # espacios al inicio
  out = replace(out, r"\s+$" => "")     # espacios al final

  # Correciones para normalizar una tabla como la de Abril 6
  correcciones = Dict("DISTRITO FEDERAL" => "CIUDAD DE MÉXICO",
                      "MEXICO" => "MÉXICO",
                      "MICHOACAN" => "MICHOACÁN",
                      "NUEVO LEON" => "NUEVO LEÓN",
                      "QUERETARO" => "QUERÉTARO",
                      "SAN LUIS POTOSI" => "SAN LUIS POTOSÍ",
                      "YUCATAN" => "YUCATÁN",
                      "MASCULINO" => "M",
                      "FEMENINO" => "F")

  for (k, v) in correcciones
    out = replace(out, k => v)
  end

  out = split(out, r"\s{2,}")           # mas de dos espacios define las entradas

  for i in index_fechas
    out[i] = procesa_fecha(out[i])  # fechas en formato ISO (de ser posible)
  end

  return out
end

#Función principal que toma el pdf y escribe csv correspondiente
function scraping(archivo_pdf, archivo_csv;
                  procedencia=false,
                  fecha_llegada=false,
                  index_fechas=[5])

  #Obtenemos los casos en un array:
  casos = procesa_fila.(eliminar_nocasos(texto_pdf(archivo_pdf)), Ref(index_fechas))
  # Ref evita el broadcasting con el segundo argumento (la lista de las columnas
  # que se tienen que considerar fechas)

  header = ["Número_caso", "Estado", "Sexo", "Edad", "Fecha_síntomas",
            "Situación"]
  # TODO: evitar utilizar acentos para nombres de las columnas?

  # Esto se puede cambiar para modificar el header
  procedencia ? push!(header, "País_fuente") : nothing
  fecha_llegada ? push!(header, "Fecha_llegada") : nothing

  #Escribe el archivo
  open(archivo_csv, "w") do io
    writedlm(io, [header], ',')
    writedlm(io, casos, ',')  # esto escribe todas las filas con separadores y \n
  end

  return "Done"
end

# Si el nombre del csv no se proporciona, se utiliza la misma base
function scraping(archivo_pdf; kwargs...)
  archivo_csv = splitext(basename(archivo_pdf))[1] * ".csv"
  scraping(archivo_pdf, archivo_csv; kwargs...)
end

function parse_commandline()
  settings = ArgParseSettings("convierte PDF con casos de covid-19 a archivo csv",
                              suppress_warnings=true)

  @add_arg_table! settings begin
    "pdf"
        help = "el archivo PDF con la tabla"
        required = true
    "--outfile", "-o"
        help = "el archivo csv donde se va a escribir el resultado (default: el nombre del PDF)"
        arg_type = String
        default = ""
    "--procedencia"
        help = "utilizar la columna `Procedencia` (default: falso)"
        action = :store_true
        default = false
    "--llegada"
        help = "utilizar la columna `Fecha de llegada a México` (default: falso)"
        action = :store_true
        default = false
  end

  return parse_args(ARGS, settings)
end

function main()
  parsed_args = parse_commandline()

  archivo_pdf = parsed_args["pdf"]
  @assert endswith(archivo_pdf, ".pdf")

  archivo_csv = parsed_args["outfile"]
  procedencia = parsed_args["procedencia"]
  fecha_llegada = parsed_args["llegada"]

  kwargs = Dict{Symbol,Any}()
  kwargs[:procedencia] = procedencia
  kwargs[:fecha_llegada] = fecha_llegada
  # trata de adivinar que fecha_llegada es una columna de fecha
  fecha_llegada ? kwargs[:index_fechas] = [5,8] : nothing

  if length(archivo_csv) == 0
    status = scraping(archivo_pdf; kwargs...)
  else
    @assert endswith(archivo_csv, ".csv")
    status = scraping(archivo_pdf, archivo_csv; kwargs...)
  end

  println(status)
end

if abspath(PROGRAM_FILE) == @__FILE__
  main()
end

#Función que hace el scraping de las dos tablas de datos diarias:
#fecha en formato "yyyymmdd"
function daily_scraping(fecha)

  año_mes = fecha[1:6]

  #Casos positivos:
  scraping("Documents/Tablas casos/$(año_mes)/positivos_$(fecha).pdf")
  #Casos sospechosos:
  scraping("Documents/Tablas casos/$(año_mes)/sospechosos_$(fecha).pdf")

  #Mueve los archivos a la carpeta respectiva:
  mv("Documents/Tablas casos/$(año_mes)/positivos_$(fecha).csv", "Daily data/$(año_mes)/positivos_$(fecha).csv")
  mv("Documents/Tablas casos/$(año_mes)/sospechosos_$(fecha).csv", "Daily data/$(año_mes)/sospechosos_$(fecha).csv")

  return "Done"
end


text:
##################### compute loss per pattern ############################

# type
struct SemFIML{
        INV <: AbstractArray,
        C <: AbstractArray,
        L <: AbstractArray,
        M <: AbstractArray,
        IM <: AbstractArray,
        I <: AbstractArray,
        T <: AbstractArray,
        U,
        V,
        W} <: LossFunction
    inverses::INV #preallocated inverses of imp_cov
    choleskys::C #preallocated choleskys
    logdets::L #logdets of implied covmats
    imp_mean::IM
    meandiff::M
    imp_inv::I
    mult::T # is this type known?
    objective::U
    grad::V
    interaction::W
end

function SemFIML(observed::O where {O <: SemObs},
    imply::I where {I <: Imply}, objective, grad)

    if isnothing(imply.imp_mean)
        error("A model implied meanstructure is needed for FIML")
    end

    inverses = broadcast(x -> zeros(x, x), Int64.(observed.pattern_nvar_obs))
    choleskys = Array{Cholesky{Float64,Array{Float64,2}},1}(undef, length(inverses))

    n_patterns = size(observed.rows, 1)
    logdets = zeros(n_patterns)

    imp_mean = zeros.(Int64.(observed.pattern_nvar_obs))
    meandiff = zeros.(Int64.(observed.pattern_nvar_obs))

    imp_inv = zeros(size(observed.data, 2), size(observed.data, 2))
    mult = similar.(inverses)

    return SemFIML(
    inverses,
    choleskys,
    logdets,
    imp_mean,
    meandiff,
    imp_inv,
    mult,
    copy(objective),
    copy(grad),
    nothing
    )
end

function F_one_pattern(imp_mean, meandiff, inverse, obs_cov, obs_mean, logdet, N)
    F = logdet
    @. meandiff = obs_mean - imp_mean
    F += dot(meandiff, inverse, meandiff)
    if N > one(N)
        F += dot(obs_cov, inverse)
    end
    F = N*F
    return F
end

function F_FIML(F, rows, semfiml, model)
    for i = 1:size(rows, 1)
        let (imp_mean, meandiff, inverse, obs_cov, obs_mean, logdet, N) =
            (semfiml.imp_mean[i],
            semfiml.meandiff[i],
            semfiml.inverses[i],
            model.observed.obs_cov[i],
            model.observed.obs_mean[i],
            semfiml.logdets[i],
            model.observed.pattern_n_obs[i])

            F += F_one_pattern(imp_mean, meandiff, inverse, obs_cov, obs_mean, logdet, N)
        end
    end
    return F
end

function (semfiml::SemFIML)(
    par,
    model::Sem{O, I, L, D}
    ) where
    {O <: SemObs,
    L <: Loss,
    I <: Imply,
    D <: SemDiff}

    if !check_fiml_2(semfiml, model) return Inf end

    copy_per_pattern!(semfiml, model)
    #batch_cholesky!(semfiml, model)
    batch_sym_inv_update!(semfiml, model)
    semfiml.logdets .= -logdet.(semfiml.inverses)

    F = zero(eltype(par))
    F = F_FIML(F, model.observed.rows, semfiml, model)
    return F
end

#= function (semfiml::SemFIML)(
    par,
    model::Sem{O, I, L, D}
    ) where
    {O <: SemObs,
    L <: Loss,
    I <: Imply,
    D <: SemDiff}

    if !check_fiml(semfiml, model) return Inf end

    copy_per_pattern!(semfiml, model)
    batch_cholesky!(semfiml, model)
    semfiml.logdets .= logdet.(semfiml.choleskys)
    batch_inv!(semfiml, model)

    F = zero(eltype(par))
    F = F_FIML(F, model.observed.rows, semfiml, model)
    return F
end =#

# Helper functions
function copy_per_pattern!(inverses, source_inverses, means, source_means, patterns)
    @views for i = 1:size(patterns, 1)
        inverses[i] .=
            source_inverses[
                patterns[i],
                patterns[i]]
    end

    @views for i = 1:size(patterns, 1)
        means[i] .=
            source_means[patterns[i]]
    end
end

function copy_per_pattern!(inverses, source_inverses, means, source_means, patterns, which_source_pattern)
    @views for i = 1:size(patterns, 1)
        inverses[i] .=
            source_inverses[which_source_pattern[i]][
                patterns[i],
                patterns[i]]
    end

    @views for i = 1:size(patterns, 1)
        means[i] .=
            source_means[which_source_pattern[i]][patterns[i]]
    end
end

copy_per_pattern!(
    semfiml::Union{SemFIML, ∇SemFIML},
    model::M where {M <: AbstractSem}) =
    copy_per_pattern!(
        semfiml.inverses,
        model.imply.imp_cov,
        semfiml.imp_mean,
        model.imply.imp_mean,
        model.observed.patterns)

copy_per_pattern!(semfiml::Union{SemFIML, ∇SemFIML}, model::Sem{O, I, L, D}) where
    {O <: SemObsMissing, L , I <: ImplyDefinition, D} =
    copy_per_pattern!(
        semfiml.inverses,
        model.imply.imp_cov,
        semfiml.imp_mean,
        model.imply.imp_mean,
        semfiml.interaction.missing_patterns,
        semfiml.interaction.group_imp_per_comb)

function batch_cholesky!(semfiml::Union{SemFIML, ∇SemFIML}, model)
    for i = 1:size(semfiml.inverses, 1)
        semfiml.choleskys[i] = cholesky!(Hermitian(semfiml.inverses[i]))
    end
    return true
end

function batch_cholesky!(semfiml::Union{SemFIML, ∇SemFIML}, model::Sem{O, I, L, D}) where
    {O <: SemObsMissing, L, I <: ImplyDefinition, D}
    for i = 1:size(semfiml.inverses, 1)
        semfiml.choleskys[i] = cholesky!(Hermitian(semfiml.inverses[i]); check = false)
        if !isposdef(semfiml.choleskys[i]) return false end
    end
    return true
end

function check_fiml(semfiml, model)
    copyto!(semfiml.imp_inv, model.imply.imp_cov)
    a = cholesky!(Hermitian(semfiml.imp_inv); check = false)
    return isposdef(a)
end

function check_fiml_2(semfiml, model)
    copyto!(semfiml.imp_inv, model.imply.imp_cov)
    semfiml.choleskys[1] = cholesky!(Hermitian(semfiml.imp_inv); check = false)
    return isposdef(semfiml.choleskys[1])
end

##################### definition variables ############################
# interaction
struct SemInteractionFIMLDefition{A, B, C, D} <: SemInteraction
    rows_per_comb::A
    group_obs_per_comb::B
    group_imp_per_comb::C
    missing_patterns::D
end

# constructor
function SemFIML(
        observed::O where {O <: SemObsMissing},
        imply::I where {I <: ImplySymbolicDefinition},
        objective,
        grad)

        group_obs = findrow.(1:Int64(observed.n_obs), [observed.rows])
        group_imp = findrow.(1:Int64(observed.n_obs), [imply.rows])

        rows_per_comb = Vector{Vector{Int64}}()
        group_obs_per_comb = Vector{Int64}()
        group_imp_per_comb = Vector{Int64}()

        for i in 1:Int64(observed.n_obs)
            new = true
            for j in 1:length(rows_per_comb)
                if (group_obs[i] == group_obs_per_comb[j]) && (group_imp[i] == group_imp_per_comb[j])
                    push!(rows_per_comb[j], i)
                    new = false
                end
            end
            if new
                push!(rows_per_comb, [i])
                push!(group_obs_per_comb, group_obs[i])
                push!(group_imp_per_comb, group_imp[i])
            end
        end

        sorted_ind = sortperm(collect(zip(group_imp_per_comb, group_obs_per_comb)))
        rows_per_comb = rows_per_comb[sorted_ind]
        group_obs_per_comb = group_obs_per_comb[sorted_ind]
        group_imp_per_comb = group_imp_per_comb[sorted_ind]

        # inverses
        nvar_obs_per_comb = Int64.(observed.pattern_nvar_obs[group_obs_per_comb])
        inverses = zeros.(nvar_obs_per_comb, nvar_obs_per_comb)

        # choleskys
        choleskys = Vector{Cholesky{Float64, Matrix{Float64}}}(undef, length(inverses))

        n_patterns = size(group_obs_per_comb, 1)

        logdets = zeros(size(inverses, 1))

        imp_mean = zeros.(size.(inverses, 1))
        meandiff = zeros.(size.(inverses, 1))

        imp_inv = similar(imply.imp_cov)
        mult = similar.(inverses)

        # missing patterns
        missing_patterns = observed.patterns[group_obs_per_comb]

    return SemFIML(
    inverses,
    choleskys,
    logdets,
    imp_mean,
    meandiff,
    imp_inv,
    mult,
    copy(objective),
    copy(grad),
    SemInteractionFIMLDefition(rows_per_comb, group_obs_per_comb, group_imp_per_comb, missing_patterns)
    )
end

# loss function for defvars + fiml
function (semfiml::SemFIML)(par, model::Sem{O, I, L, D}) where
    {O <: SemObsMissing, L <: Loss, I <: ImplyDefinition, D <: SemFiniteDiff}

    copy_per_pattern!(semfiml, model)
    if !batch_cholesky!(semfiml, model) return Inf end
    semfiml.logdets .= logdet.(semfiml.choleskys)
    batch_inv!(semfiml, model)

    F = zero(eltype(par))
    F = F_FIML(F, semfiml.interaction.rows_per_comb, semfiml, model)
    return F
end