text:
using Parsers, Test, Dates

import Parsers: INVALID, OK, SENTINEL, QUOTED, DELIMITED, NEWLINE, EOF, INVALID_QUOTED_FIELD, INVALID_DELIMITER, OVERFLOW, ESCAPED_STRING

@testset "Parsers" begin

@testset "Core Parsers.xparse" begin

sentinels = ["NANA", "NAN", "NA"]

testcases = [
    (str="", kwargs=(), x=0, code=(INVALID | EOF), vpos=1, vlen=0, tlen=0),
    (str="", kwargs=(sentinel=missing,), x=0, code=(SENTINEL | EOF), vpos=1, vlen=0, tlen=0),
    (str=" ", kwargs=(), x=0, code=(INVALID | EOF), vpos=1, vlen=0, tlen=1),
    (str=" ", kwargs=(sentinel=missing,), x=0, code=(INVALID | EOF), vpos=1, vlen=0, tlen=1),
    (str=" -", kwargs=(sentinel=missing,), x=0, code=(INVALID | EOF), vpos=1, vlen=2, tlen=2),
    (str=" +", kwargs=(sentinel=missing,), x=0, code=(INVALID | EOF), vpos=1, vlen=2, tlen=2),
    (str="-", kwargs=(sentinel=missing,), x=0, code=(INVALID | EOF), vpos=1, vlen=1, tlen=1),
    (str=" {-", kwargs=(sentinel=missing,), x=0, code=(QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=3, vlen=1, tlen=3),
    (str="{+", kwargs=(sentinel=missing,), x=0, code=(QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=2, vlen=1, tlen=2),
    (str=" {+,", kwargs=(sentinel=missing,), x=0, code=(QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=3, vlen=1, tlen=4),
    (str="{-,", kwargs=(sentinel=missing,), x=0, code=(QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=2, vlen=1, tlen=3),
    (str="+,", kwargs=(sentinel=missing,), x=0, code=(INVALID | DELIMITED), vpos=1, vlen=1, tlen=2),
    (str="-,", kwargs=(sentinel=missing,), x=0, code=(INVALID | DELIMITED), vpos=1, vlen=1, tlen=2),
    (str=" {-},", kwargs=(sentinel=missing,), x=0, code=(INVALID | QUOTED | DELIMITED), vpos=3, vlen=1, tlen=5),
    (str="{+} ,", kwargs=(sentinel=missing,), x=0, code=(INVALID | QUOTED | DELIMITED), vpos=2, vlen=1, tlen=5),
    (str="{", kwargs=(), x=0, code=(QUOTED | INVALID_QUOTED_FIELD), vpos=2, vlen=-1, tlen=1),
    (str="{}", kwargs=(), x=0, code=(INVALID | QUOTED | EOF), vpos=2, vlen=0, tlen=2),
    (str=" {", kwargs=(), x=0, code=(QUOTED | INVALID_QUOTED_FIELD), vpos=3, vlen=-2, tlen=2),
    (str=" {\\\\", kwargs=(), x=0, code=(QUOTED | INVALID_QUOTED_FIELD | ESCAPED_STRING | EOF), vpos=3, vlen=0, tlen=4),
    (str=" {\\}} ", kwargs=(), x=0, code=(QUOTED | INVALID | ESCAPED_STRING | EOF), vpos=3, vlen=2, tlen=6),
    (str=" {\\\\}", kwargs=(), x=0, code=(INVALID | QUOTED | ESCAPED_STRING | EOF), vpos=3, vlen=2, tlen=5),

    (str=" {}", kwargs=(), x=0, code=(INVALID | QUOTED | EOF), vpos=3, vlen=0, tlen=3),
    (str=" { }", kwargs=(), x=0, code=(INVALID | QUOTED | EOF), vpos=3, vlen=1, tlen=4),
    (str=" {,} ", kwargs=(), x=0, code=(INVALID | QUOTED | EOF), vpos=3, vlen=1, tlen=5),
    (str=" { } ", kwargs=(), x=0, code=(INVALID | QUOTED | EOF), vpos=3, vlen=1, tlen=5),
    (str=" {\t", kwargs=(), x=0, code=(INVALID_QUOTED_FIELD | QUOTED | EOF), vpos=3, vlen=-2, tlen=3),

    (str="NA", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | EOF), vpos=1, vlen=2, tlen=2),
    (str="£", kwargs=(sentinel=["£"],), x=0, code=(SENTINEL | EOF), vpos=1, vlen=2, tlen=2),
    (str="NA2", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | EOF | INVALID_DELIMITER), vpos=1, vlen=3, tlen=3),
    (str="NAN", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | EOF), vpos=1, vlen=3, tlen=3),
    (str="NANA", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | EOF), vpos=1, vlen=4, tlen=4),
    (str=" NA", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | EOF), vpos=1, vlen=3, tlen=3),
    (str="{NAN", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=2, vlen=3, tlen=4),
    (str="{NANA}", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | QUOTED | EOF), vpos=2, vlen=4, tlen=6),
    (str=" {NA", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=3, vlen=2, tlen=4),
    (str=" {NANA}", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=4, tlen=7),
    (str=" { NAN}", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=4, tlen=7),
    (str=" {NAN }", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=4, tlen=7),
    (str=" {NA} ", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=2, tlen=6),
    (str=" { NANA} ", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=5, tlen=9),
    (str=" {\tNA", kwargs=(sentinel=sentinels,), x=0, code=(SENTINEL | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=3, vlen=3, tlen=5),

    (str="-", kwargs=(sentinel=["-"],), x=0, code=(SENTINEL | EOF), vpos=1, vlen=1, tlen=1),
    (str=" +", kwargs=(sentinel=["+"],), x=0, code=(SENTINEL | EOF), vpos=1, vlen=2, tlen=2),
    (str="+1 ", kwargs=(sentinel=["+1"],), x=1, code=(SENTINEL | EOF), vpos=1, vlen=3, tlen=3),
    (str=" -1 ", kwargs=(sentinel=["-1"],), x=-1, code=(SENTINEL | EOF), vpos=1, vlen=4, tlen=4),
    (str="{1", kwargs=(sentinel=["1"],), x=1, code=(SENTINEL | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=2, vlen=1, tlen=2),
    (str="{1 ", kwargs=(sentinel=["1"],), x=1, code=(SENTINEL | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=2, vlen=2, tlen=3),
    (str="{-1}", kwargs=(sentinel=["-1"],), x=-1, code=(SENTINEL | QUOTED | EOF), vpos=2, vlen=2, tlen=4),
    (str=" {+1", kwargs=(sentinel=["+1"],), x=1, code=(SENTINEL | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=3, vlen=2, tlen=4),
    (str=" {-}", kwargs=(sentinel=["-"],), x=0, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=1, tlen=4),
    (str=" { +}", kwargs=(sentinel=["+"],), x=0, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=2, tlen=5),
    (str=" {-1 }", kwargs=(sentinel=["-1"],), x=-1, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=3, tlen=6),
    (str=" {+1} ", kwargs=(sentinel=["+1"],), x=1, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=2, tlen=6),
    (str=" { 1} ", kwargs=(sentinel=["1"],), x=1, code=(SENTINEL | QUOTED | EOF), vpos=3, vlen=2, tlen=6),
    (str=" {\t-", kwargs=(sentinel=["-"],), x=0, code=(SENTINEL | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=3, vlen=2, tlen=4),

    (str="+a ", kwargs=(sentinel=["+1"],), x=0, code=(INVALID | EOF | INVALID_DELIMITER), vpos=1, vlen=3, tlen=3),
    (str=" -a ", kwargs=(sentinel=["-1"],), x=0, code=(INVALID | EOF | INVALID_DELIMITER), vpos=1, vlen=4, tlen=4),
    (str="{a", kwargs=(sentinel=["1"],), x=0, code=(INVALID | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=2, vlen=0, tlen=2),
    (str="{-a}", kwargs=(sentinel=["-1"],), x=0, code=(INVALID | EOF | QUOTED), vpos=2, vlen=2, tlen=4),
    (str=" {+1a", kwargs=(sentinel=["+1"],), x=1, code=(SENTINEL | QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=3, vlen=2, tlen=5),
    (str=" {-1a }", kwargs=(sentinel=["-1"],), x=-1, code=(SENTINEL | QUOTED | INVALID | EOF), vpos=3, vlen=4, tlen=7),
    (str=" {+a} ", kwargs=(sentinel=["+1"],), x=0, code=(INVALID | QUOTED | EOF), vpos=3, vlen=2, tlen=6),
    (str=" { a} ", kwargs=(sentinel=["1"],), x=0, code=(INVALID | QUOTED | EOF), vpos=3, vlen=2, tlen=6),

    (str="-", kwargs=(), x=0, code=(INVALID | EOF), vpos=1, vlen=1, tlen=1),
    (str="0", kwargs=(), x=0, code=(OK | EOF), vpos=1, vlen=1, tlen=1),
    (str=" +", kwargs=(), x=0, code=(INVALID | EOF), vpos=1, vlen=2, tlen=2),
    (str=" +1", kwargs=(), x=1, code=(OK | EOF), vpos=1, vlen=3, tlen=3),
    (str="-1 ", kwargs=(), x=-1, code=(OK | EOF), vpos=1, vlen=3, tlen=3),
    (str=" +1 ", kwargs=(), x=1, code=(OK | EOF), vpos=1, vlen=4, tlen=4),
    (str="{1", kwargs=(), x=1, code=(OK | QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=2, vlen=1, tlen=2),
    (str="{1 ", kwargs=(), x=1, code=(OK | QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=2, vlen=2, tlen=3),
    (str="{-1}", kwargs=(), x=-1, code=(OK | QUOTED | EOF), vpos=2, vlen=2, tlen=4),
    (str=" {+1", kwargs=(), x=1, code=(OK | QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=3, vlen=2, tlen=4),
    (str=" {-}", kwargs=(), x=0, code=(INVALID | QUOTED | EOF), vpos=3, vlen=1, tlen=4),
    (str=" { +}", kwargs=(), x=0, code=(INVALID | QUOTED | EOF), vpos=3, vlen=2, tlen=5),
    (str=" {-1 }", kwargs=(), x=-1, code=(OK | QUOTED | EOF), vpos=3, vlen=3, tlen=6),
    (str=" {+1} ", kwargs=(), x=1, code=(OK | QUOTED | EOF), vpos=3, vlen=2, tlen=6),
    (str=" { 1} ", kwargs=(), x=1, code=(OK | QUOTED | EOF), vpos=3, vlen=2, tlen=6),
    (str=" {\t-", kwargs=(), x=0, code=(QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=3, vlen=2, tlen=4),
    (str="{1ab\\\\c}", kwargs=(), x=1, code=(OK | QUOTED | ESCAPED_STRING | EOF | INVALID), vpos=2, vlen=6, tlen=8),
    (str="{1\\\\abc,", kwargs=(), x=1, code=(OK | QUOTED | ESCAPED_STRING | EOF | INVALID_QUOTED_FIELD), vpos=2, vlen=6, tlen=8),
    (str="{1abc},", kwargs=(sentinel=["1abc"],), x=1, code=(SENTINEL | QUOTED | DELIMITED), vpos=2, vlen=4, tlen=7),
    (str="{1abc", kwargs=(sentinel=["1abc"],), x=1, code=(SENTINEL | QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=2, vlen=4, tlen=5),

    (str="922337203,", kwargs=(sentinel=["92233"],), x=922337203, code=(OK | DELIMITED), vpos=1, vlen=9, tlen=10),
    (str="92233,", kwargs=(sentinel=["92233"],), x=92233, code=(SENTINEL | DELIMITED), vpos=1, vlen=5, tlen=6),
    (str="92233  ", kwargs=(sentinel=["92233"],), x=92233, code=(SENTINEL | EOF), vpos=1, vlen=7, tlen=7),
    (str="{92233  ,", kwargs=(sentinel=["92233"],), x=92233, code=(SENTINEL | QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=2, vlen=7, tlen=9),
    (str="{92233} ,", kwargs=(sentinel=["92233"],), x=92233, code=(SENTINEL | QUOTED | DELIMITED), vpos=2, vlen=5, tlen=9),
    (str=" { 92233 },", kwargs=(sentinel=["92233"],), x=92233, code=(SENTINEL | QUOTED | DELIMITED), vpos=3, vlen=7, tlen=11),
    (str="922337203685477580", kwargs=(), x=922337203685477580, code=(OK | EOF), vpos=1, vlen=18, tlen=18),
    (str="9223372036854775808", kwargs=(), x=-9223372036854775808, code=(OVERFLOW | EOF), vpos=1, vlen=19, tlen=19),
    (str="9223372036854775808a", kwargs=(sentinel=["9223372036854775808a"],), x=-9223372036854775808, code=(SENTINEL | EOF), vpos=1, vlen=20, tlen=20),
    (str="{9223372036854775808a", kwargs=(sentinel=["9223372036854775808a"],), x=-9223372036854775808, code=(SENTINEL | QUOTED | EOF | INVALID_QUOTED_FIELD), vpos=2, vlen=20, tlen=21),

    (str="{9223372036854775808", kwargs=(), x=-9223372036854775808, code=(OVERFLOW | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=2, vlen=19, tlen=20),
    (str="{9223372036854775807", kwargs=(), x=9223372036854775807, code=(OK | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=2, vlen=19, tlen=20),
    (str="{9223372036854775807a", kwargs=(), x=9223372036854775807, code=(OK | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=2, vlen=19, tlen=21),
    (str="9223372036854775807a", kwargs=(), x=9223372036854775807, code=(OK | EOF | INVALID_DELIMITER), vpos=1, vlen=20, tlen=20),
    (str="9223372036854775807,", kwargs=(), x=9223372036854775807, code=(OK | DELIMITED), vpos=1, vlen=19, tlen=20),
    (str="9223372036854775807", kwargs=(sentinel=["9223372036854775807"],), x=9223372036854775807, code=(SENTINEL | EOF), vpos=1, vlen=19, tlen=19),
    (str="{9223372036854775807a", kwargs=(sentinel=["9223372036854775807a"],), x=9223372036854775807, code=(SENTINEL | EOF | QUOTED | INVALID_QUOTED_FIELD), vpos=2, vlen=20, tlen=21),

    (str="9223372036854775808", kwargs=(sentinel=["92233"],), x=-9223372036854775808, code=(OVERFLOW | EOF), vpos=1, vlen=19, tlen=19),
    (str="922337203685477580", kwargs=(sentinel=["92233"],), x=922337203685477580, code=(OK | EOF), vpos=1, vlen=18, tlen=18),
    (str="922337203685477580,", kwargs=(), x=922337203685477580, code=(OK | DELIMITED), vpos=1, vlen=18, tlen=19),
    (str="922337203685477580,", kwargs=(sentinel=["92233"],), x=922337203685477580, code=(OK | DELIMITED), vpos=1, vlen=18, tlen=19),
    (str="9223372036854775800000,", kwargs=(sentinel=["92233"],), x=-80, code=(DELIMITED | INVALID_DELIMITER | OVERFLOW), vpos=1, vlen=22, tlen=23),

    (str="1;", kwargs=(), x=1, code=(OK | EOF | INVALID_DELIMITER), vpos=1, vlen=2, tlen=2),
    (str="1;", kwargs=(delim=UInt8(';'),), x=1, code=(OK | DELIMITED), vpos=1, vlen=1, tlen=2),
    (str="1abc;", kwargs=(delim=UInt8(';'),), x=1, code=(OK | DELIMITED | INVALID_DELIMITER), vpos=1, vlen=4, tlen=5),
    (str="1\n", kwargs=(), x=1, code=(OK | NEWLINE | EOF), vpos=1, vlen=1, tlen=2),
    (str="1\r", kwargs=(), x=1, code=(OK | NEWLINE | EOF), vpos=1, vlen=1, tlen=2),
    (str="1\r\n", kwargs=(), x=1, code=(OK | NEWLINE | EOF), vpos=1, vlen=1, tlen=3),
    (str="1\n2", kwargs=(), x=1, code=(OK | NEWLINE), vpos=1, vlen=1, tlen=2),
    (str="1\r2", kwargs=(), x=1, code=(OK | NEWLINE), vpos=1, vlen=1, tlen=2),
    (str="1\r\n2", kwargs=(), x=1, code=(OK | NEWLINE), vpos=1, vlen=1, tlen=3),
    (str="1a\n", kwargs=(), x=1, code=(OK | NEWLINE | EOF | INVALID_DELIMITER), vpos=1, vlen=2, tlen=3),
    (str="1a\r", kwargs=(), x=1, code=(OK | NEWLINE | EOF | INVALID_DELIMITER), vpos=1, vlen=2, tlen=3),
    (str="1a\r\n", kwargs=(), x=1, code=(OK | NEWLINE | EOF | INVALID_DELIMITER), vpos=1, vlen=2, tlen=4),
    (str="1a\n2", kwargs=(), x=1, code=(OK | NEWLINE | INVALID_DELIMITER), vpos=1, vlen=2, tlen=3),
    (str="1a\r2", kwargs=(), x=1, code=(OK | NEWLINE | INVALID_DELIMITER), vpos=1, vlen=2, tlen=3),
    (str="1a\r\n2", kwargs=(), x=1, code=(OK | NEWLINE | INVALID_DELIMITER), vpos=1, vlen=2, tlen=4),

    (str="1,,,2,null,4", kwargs=(), x=1, code=(OK | DELIMITED), vpos=1, vlen=1, tlen=2),
    (str="1,,,2,null,4", kwargs=(ignorerepeated=true,), x=1, code=(OK | DELIMITED), vpos=1, vlen=1, tlen=4),
    (str="1,", kwargs=(ignorerepeated=true,), x=1, code=(OK | DELIMITED), vpos=1, vlen=1, tlen=2),
    (str="1,,", kwargs=(ignorerepeated=true,), x=1, code=(OK | DELIMITED), vpos=1, vlen=1, tlen=3),
    (str="1,,,", kwargs=(ignorerepeated=true,), x=1, code=(OK | DELIMITED), vpos=1, vlen=1, tlen=4),
    (str="1::2", kwargs=(delim="::",), x=1, code=(OK | DELIMITED), vpos=1, vlen=1, tlen=3),
    (str="1::::2", kwargs=(ignorerepeated=true, delim="::"), x=1, code=(OK | DELIMITED), vpos=1, vlen=1, tlen=5),
    (str="1a::::2", kwargs=(ignorerepeated=true, delim="::"), x=1, code=(OK | DELIMITED | INVALID_DELIMITER), vpos=1, vlen=2, tlen=6),
    (str="1[][]", kwargs=(delim="[]", ignorerepeated = true), x = 1, code=(OK | DELIMITED), vpos=1, vlen=1, tlen=5),
    (str="1a[][]", kwargs=(delim="[]", ignorerepeated = true), x = 1, code=(OK | DELIMITED | INVALID_DELIMITER), vpos=1, vlen=2, tlen=6),
    (str="1a[][]", kwargs=(delim="[]",), x = 1, code=(OK | DELIMITED | INVALID_DELIMITER), vpos=1, vlen=2, tlen=4),
    # ignorerepeated
    (str="1a,,", kwargs=(ignorerepeated=true,), x=1, code=(OK | DELIMITED | INVALID_DELIMITER), vpos=1, vlen=2, tlen=4),
    (str="1a,,2", kwargs=(ignorerepeated=true,), x=1, code=(OK | DELIMITED | INVALID_DELIMITER), vpos=1, vlen=2, tlen=4),
];

for useio in (false, true)
    for (oq, cq, e) in ((UInt8('"'), UInt8('"'), UInt8('"')), (UInt8('"'), UInt8('"'), UInt8('\\')), (UInt8('{'), UInt8('}'), UInt8('\\')))
        for (i, case) in enumerate(testcases)
            str = replace(replace(replace(case.str, '{'=>Char(oq)), '}'=>Char(cq)), '\\'=>Char(e))
            source = useio ? IOBuffer(str) : str
            x, code, vpos, vlen, tlen = Parsers.xparse(Int64, source; openquotechar=oq, closequotechar=cq, escapechar=e, case.kwargs...)
            if x != case.x || code != case.code || vpos != case.vpos || vlen != case.vlen || tlen != case.tlen
                println("ERROR on useio=$useio, case=$i, oq='$(Char(oq))', cq='$(Char(cq))', e='$(Char(e))', str=$(escape_string(str)), $case")
                source = useio ? IOBuffer(str) : str
                x, code, vpos, vlen, tlen = Parsers.xparse(Int64, source; debug=true, openquotechar=oq, closequotechar=cq, escapechar=e, case.kwargs...)
            end
            @test x == case.x
            @test code == case.code
            @test vpos == case.vpos
            @test vlen == case.vlen
            @test tlen == case.tlen
        end
    end
end

# strings
oq = UInt8('{')
cq = UInt8('}')
e = UInt8('\\')
for (i, case) in enumerate(testcases)
    str = replace(replace(replace(case.str, '{'=>Char(oq)), '}'=>Char(cq)), '\\'=>Char(e))
    x, code, vpos, vlen, tlen = Parsers.xparse(String, case.str; openquotechar=oq, closequotechar=cq, escapechar=e, case.kwargs...)
    if vpos != case.vpos || vlen != case.vlen || tlen != case.tlen
        println("ERROR on case=$i, oq='$(Char(oq))', cq='$(Char(cq))', e='$(Char(e))', str=$(escape_string(str)), $case")
        x, code, vpos, vlen, tlen = Parsers.xparse(String, case.str; debug=true, openquotechar=oq, closequotechar=cq, escapechar=e, case.kwargs...)
    end
    if !Parsers.invalidquotedfield(code)
        @test vpos == case.vpos
        @test vlen == case.vlen
        @test tlen == case.tlen
    end
end

end # @testset "Core Parsers.xparse"

@testset "ints" begin

# test lots of ints
@time for i in typemin(Int64):100_000_000_000_000:typemax(Int64)
    str = string(i)
    x, code, vpos, vlen, tlen = Parsers.xparse(Int64, str)
    @test string(x) == str
    @test code == OK | EOF
    @test vlen == length(str)
end

end # @testset "ints"

@testset "bools" begin

trues1 = ["T"]
falses1 = ["F"]
trues2 = ["truthy"]
falses2 = ["falsy"]

testcases = [
    (str="", kwargs=(), x=false, code=(INVALID | EOF), vpos=1, vlen=0),
    (str="t", kwargs=(), x=false, code=(INVALID | EOF), vpos=1, vlen=1),
    (str="tr", kwargs=(), x=false, code=(INVALID | EOF), vpos=1, vlen=2),
    (str="tru", kwargs=(), x=false, code=(INVALID | EOF), vpos=1, vlen=3),
    (str="true", kwargs=(), x=true, code=(OK | EOF), vpos=1, vlen=4),

    (str="f", kwargs=(), x=false, code=(INVALID | EOF), vpos=1, vlen=1),
    (str="fa", kwargs=(), x=false, code=(INVALID | EOF), vpos=1, vlen=2),
    (str="fal", kwargs=(), x=false, code=(INVALID | EOF), vpos=1, vlen=3),
    (str="fals", kwargs=(), x=false, code=(INVALID | EOF), vpos=1, vlen=4),
    (str="false", kwargs=(), x=false, code=(OK | EOF), vpos=1, vlen=5),

    (str="t,", kwargs=(), x=false, code=(INVALID | DELIMITED), vpos=1, vlen=1),
    (str="tr,", kwargs=(), x=false, code=(INVALID | DELIMITED), vpos=1, vlen=2),
    (str="tru,", kwargs=(), x=false, code=(INVALID | DELIMITED), vpos=1, vlen=3),
    (str="true,", kwargs=(), x=true, code=(OK | DELIMITED), vpos=1, vlen=4),

    (str="f,", kwargs=(), x=false, code=(INVALID | DELIMITED), vpos=1, vlen=1),
    (str="fa,", kwargs=(), x=false, code=(INVALID | DELIMITED), vpos=1, vlen=2),
    (str="fal,", kwargs=(), x=false, code=(INVALID | DELIMITED), vpos=1, vlen=3),
    (str="fals,", kwargs=(), x=false, code=(INVALID | DELIMITED), vpos=1, vlen=4),
    (str="false,", kwargs=(), x=false, code=(OK | DELIMITED), vpos=1, vlen=5),

    (str="0", kwargs=(), x=false, code=(OK | EOF), vpos=1, vlen=1),
    (str="1", kwargs=(), x=true, code=(OK | EOF), vpos=1, vlen=1),
    (str="001", kwargs=(), x=true, code=(OK | EOF), vpos=1, vlen=3),
    (str="2", kwargs=(), x=false, code=(INVALID | EOF | INVALID_DELIMITER), vpos=1, vlen=1),

    (str="t", kwargs=(trues=trues1, falses=falses1,), x=false, code=(INVALID | EOF | INVALID_DELIMITER), vpos=1, vlen=1),
    (str="T", kwargs=(trues=trues1, falses=falses1,), x=true, code=(OK | EOF), vpos=1, vlen=1),
    (str="T,", kwargs=(trues=trues1, falses=falses1,), x=true, code=(OK | DELIMITED), vpos=1, vlen=1),
    (str="Tr", kwargs=(trues=trues1, falses=falses1,), x=true, code=(OK | EOF | INVALID_DELIMITER), vpos=1, vlen=2),
    (str="F", kwargs=(trues=trues1, falses=falses1,), x=false, code=(OK | EOF), vpos=1, vlen=1),
    (str="truthy", kwargs=(trues=trues2, falses=falses2,), x=true, code=(OK | EOF), vpos=1, vlen=6),
    (str="truthyfalsy", kwargs=(trues=trues2, falses=falses2,), x=true, code=(OK | EOF | INVALID_DELIMITER), vpos=1, vlen=11),
    (str="falsytruthy", kwargs=(trues=trues2, falses=falses2,), x=false, code=(OK | EOF | INVALID_DELIMITER), vpos=1, vlen=11),
];

for useio in (false, true)
    for (i, case) in enumerate(testcases)
        x, code, vpos, vlen, tot = Parsers.xparse(Bool, useio ? IOBuffer(case.str) : case.str; case.kwargs...)
        if x != case.x || code != case.code || vlen != case.vlen
            println("error for case = $i: useio=$useio, $case")
            println(Parsers.codes(code))
        end
        @test x == case.x
        @test code == case.code
        @test vlen == case.vlen
    end
end

end # @testset "bools"

@testset "misc" begin

# additional tests for full xparse branch coverage
oq = UInt8('{')
cq = UInt8('}')
e = UInt8('\\')
str=" {\\"
x, code, vpos, vlen, tlen = Parsers.xparse(Int64, str; openquotechar=oq, closequotechar=cq, escapechar=e)
@test x == 0
@test code == QUOTED | EOF | INVALID_QUOTED_FIELD
@test tlen == 3

@test Parsers.parse(Int, "101") === 101
@test Parsers.parse(Float64, "101,101", Parsers.Options(decimal=',')) === 101.101
@test Parsers.parse(Bool, IOBuffer("true")) === true
@test_throws Parsers.Error Parsers.parse(Int, "abc")

@test Parsers.tryparse(Int, "abc") === nothing
@test Parsers.tryparse(Float32, IOBuffer("101,101"), Parsers.Options(decimal=',')) === Float32(101.101)
@test Parsers.parse(Date, "01/20/2018", Parsers.Options(dateformat="mm/dd/yyyy")) === Date(2018, 1, 20)

# https://github.com/JuliaData/CSV.jl/issues/345
x, code, vpos, vlen, tlen = Parsers.xparse(String, "\"DALLAS BLACK DANCE THEATRE\",")
@test vpos == 2
@test vlen == 26
@test tlen == 29
@test code == OK | QUOTED | DELIMITED

# https://github.com/JuliaData/CSV.jl/issues/344
str = "1,2,null,4"
pos = 1
x, code, vpos, vlen, tlen = Parsers.xparse(Int, str; pos=pos, sentinel=["null"])
pos += tlen
@test x === 1
@test code === OK | DELIMITED
@test vpos == 1
@test vlen == 1
@test pos == 3
x, code, vpos, vlen, tlen = Parsers.xparse(Int, str; pos=pos, sentinel=["null"])
pos += tlen
@test x === 2
@test code === OK | DELIMITED
@test vpos == 3
@test vlen == 1
@test pos == 5
x, code, vpos, vlen, tlen = Parsers.xparse(Int, str; pos=pos, sentinel=["null"])
pos += tlen
@test x == 0
@test code === SENTINEL | DELIMITED
@test vpos == 5
@test vlen == 4
@test pos == 10
x, code, vpos, vlen, tlen = Parsers.xparse(Int, str; pos=pos, sentinel=["null"])
pos += tlen
@test x === 4
@test code === OK | EOF
@test vpos == 10
@test vlen == 1
@test pos == 11

@test Parsers.parse(Int, SubString("101")) === 101
@test Parsers.parse(Float64, SubString("101,101"), Parsers.Options(decimal=',')) === 101.101

@test Parsers.asciival(' ')

@test_throws ArgumentError Parsers.Options(sentinel=[" "])
@test_throws ArgumentError Parsers.Options(sentinel=["\""])
@test_throws ArgumentError Parsers.Options(sentinel=[","], delim=',')
@test_throws ArgumentError Parsers.Options(sentinel=[","], delim=",")

@test Parsers.default(Int32) === Int32(0)
@test Parsers.default(Float32) === Float32(0)
@test Parsers.xparse(Int64, "10", 1, 2, Parsers.XOPTIONS) == (Int64(10), OK | EOF, 1, 2, 2)
@test Parsers.xparse(Int64, "a", 1, 1, Parsers.Options(sentinel=missing)) == (Int64(0), SENTINEL, 1, 0, 0)

@test Parsers.checkdelim!(UInt8[], 1, 0, Parsers.OPTIONS) == 1
@test Parsers.checkdelim!(codeunits(","), 1, 1, Parsers.XOPTIONS) == 2
@test Parsers.checkdelim!(codeunits("::"), 1, 2, Parsers.Options(delim="::")) == 3
@test Parsers.checkdelim!(codeunits(",,"), 1, 2, Parsers.Options(ignorerepeated=true, delim=',')) == 3
@test Parsers.checkdelim!(codeunits("::::"), 1, 4, Parsers.Options(delim="::", ignorerepeated=true)) == 5

e = Parsers.Error(Vector{UInt8}("hey"), Int64, INVALID | EOF, 1, 3)
io = IOBuffer()
showerror(io, e)
@test String(take!(io)) == "Parsers.Error (INVALID: EOF ):\ninitial value parsing failed, reached EOF\nattempted to parse Int64 from: \"hey\"\n"
e2 = Parsers.Error(IOBuffer("hey"), Int64, INVALID | EOF, 1, 3)
showerror(io, e2)
@test String(take!(io)) == "Parsers.Error (INVALID: EOF ):\ninitial value parsing failed, reached EOF\nattempted to parse Int64 from: \"hey\"\n"

@test Parsers.invalid(INVALID_DELIMITER)
@test Parsers.sentinel(SENTINEL)
@test !Parsers.sentinel(OK)
@test Parsers.quoted(QUOTED)
@test !Parsers.quoted(INVALID_DELIMITER)
@test Parsers.delimited(DELIMITED)
@test !Parsers.delimited(OK)
@test Parsers.newline(NEWLINE)
@test !Parsers.newline(DELIMITED)
@test Parsers.escapedstring(ESCAPED_STRING)
@test !Parsers.escapedstring(OK)
@test Parsers.invalidquotedfield(INVALID_QUOTED_FIELD)
@test !Parsers.invalidquotedfield(INVALID_DELIMITER)
@test Parsers.invaliddelimiter(INVALID_DELIMITER)
@test !Parsers.invaliddelimiter(INVALID_QUOTED_FIELD)
@test Parsers.overflow(OVERFLOW)
@test !Parsers.overflow(OK)
@test Parsers.quotednotescaped(QUOTED)
@test !Parsers.quotednotescaped(QUOTED | ESCAPED_STRING)

# https://github.com/JuliaData/CSV.jl/issues/454
x, code, vpos, vlen, tlen = Parsers.xparse(Float64, "\"\"", 1, 2)
@test Parsers.sentinel(code)

x, code, vpos, vlen, tlen = Parsers.xparse(String, "\"\"", 1, 2)
@test Parsers.sentinel(code)

end # @testset "misc"

include("floats.jl")
include("dates.jl")

end # @testset "Parsers"


text:
include("../get_home_dir.jl")

using DataFrames
using CSV
using PyPlot
using Random
using StatsBase
using Statistics

df_raw = CSV.read(string(PROJECT_HOME,"/data/moflux_land_data_skipyear_hourly2.csv"), DataFrame);
df_raw[!,"RAIN"] *= 2;

include(string(PROJECT_HOME,"/simulation_code/sim_vary_new_stomata_med_varyLayer.jl"));
include("time_averaging.jl")
include("mironov.jl")
include("tau_omega_funs.jl")

N = 24*365*1
istart = 24*365*3 + 1;
soil0 = 0.30;

deltaT = FT(60*60);
alpha = FT(50);
#alpha = FT(20);

alpha = FT(5);
nsoil = FT(1.2);
#nsoil = FT(1.4);

function run_sim_3layer(vcmax_par::FT, k_frac::FT, k_plant::FT,  z_soil::FT, weibB::FT, vol_factor::FT, g1::FT)
	return run_sim_varyB(vcmax_par, k_frac, weibB, FT(2), k_plant, FT(1e-4), z_soil, istart, N, soil0, vol_factor, 1e-5, 3, df_raw, g1, deltaT, alpha, nsoil,0,0);
end

c1 = 1

sim_res1 = run_sim_3layer(FT(60),FT(0.5), FT(12), FT(2000),FT(5),FT(0.4),FT(250));

factor_list = [1/3,1/2,1/1.5,1.5,2,3];

sim_list = [];

#push all of these to sim_list

parnames = ["g1","p63","Kplant","Capacitance","Zsoil"];

println("Set 0 of 4")
#change g1
partial_list = [run_sim_3layer(FT(60),FT(0.5), FT(12), FT(2000),FT(5),FT(0.4),FT(250*factor_list[i])) for i in 1:length(factor_list)];
append!(sim_list, partial_list);

println("Set 1 of 4")
#change stomatal and xylem p50
partial_list = [run_sim_3layer(FT(60),FT(0.5), FT(12), FT(2000),FT(5*factor_list[i]),FT(0.4),FT(250)) for i in 1:length(factor_list)];
append!(sim_list, partial_list);

println("Set 2 of 4")
#change conductance
partial_list = [run_sim_3layer(FT(60),FT(0.5), FT(12*factor_list[i]), FT(2000),FT(5),FT(0.4),FT(250)) for i in 1:length(factor_list)];
append!(sim_list, partial_list);

println("Set 3 of 4")
#change cap
partial_list = [run_sim_3layer(FT(60),FT(0.5), FT(12), FT(2000),FT(5),FT(0.4*factor_list[i]),FT(250)) for i in 1:length(factor_list)];
append!(sim_list, partial_list);

#change cap and cond
#sim_list = [run_sim_3layer(FT(60),FT(0.5), FT(12*factor_list[i]), FT(2000),FT(5),FT(0.4*factor_list[i]),FT(250)) for i in 1:length(factor_list)];

println("Set 4 of 4")
#change soil depth
partial_list = [run_sim_3layer(FT(60),FT(0.5), FT(12), FT(2000*factor_list[i]),FT(5),FT(0.4),FT(250)) for i in 1:length(factor_list)];
append!(sim_list, partial_list);

cm1 = mean(sim_res1[6][:,1:3],dims=2);
cmx = [mean(x[6][:,1:3],dims=2) for x in sim_list];

pdLWP = Float64.(replace(sim_res1[1].LWP_predawn, missing => NaN));

tsoil = sim_res1[1].T_SOIL .+ 273.15;
tair = sim_res1[1].T_AIR .+ 273.15;
tcan = (sim_res1[1].LW_OUT / 0.97 / K_STEFAN ) .^ 0.25;
laiM = sim_res1[1].LAI_modis;

using Distributions
using LabelledArrays
using LinearAlgebra

vodA = 0.067;
vodB = 0.82;
vodC = 0.051;

tb1 = get_TB_2(sim_res1[2][:,1], tsoil, tcan, cm1, laiM, vodA, vodB, vodC);

function optim_vod_mask(obsHi, obsVi, leafpot, smc_mod, Tsoil, Tcan, laiM, k, alpha, beta, maxiter, stop_crit,step_size, nskip, mask)
	return optim_vod_stages2(obsHi[mask], obsVi[mask], leafpot[mask], smc_mod[mask], Tsoil[mask], Tcan[mask], laiM[mask], k, alpha, beta, maxiter, stop_crit,step_size, nskip)
end

hodlist = (1:length(tb1[1])) .% 24;

mask6 = ((hodlist .== 7) + (hodlist .== (7+12))) .== 1;
mask1 = ((hodlist .== 2) + (hodlist .== (2+12))) .== 1;

tbx_all = []
tbx_6 = []
tbx_1 = []

for i in 1:length(sim_list)
	println(i)
	sim_res2 = sim_list[i];
	cm2 = cmx[i]
	ores2_all = optim_vod_stages2(tb1[1], tb1[2], cm2, sim_res2[2][:,1], tsoil, tcan, laiM, vodA, vodB, vodC, 10^5, 1e-6, 1e-5, 24);
	ores2_6 = optim_vod_mask(tb1[1], tb1[2], cm2, sim_res2[2][:,1], tsoil, tcan, laiM, vodA, vodB, vodC, 10^5, 1e-6, 1e-5, 2, mask6);
	ores2_1 = optim_vod_mask(tb1[1], tb1[2], cm2, sim_res2[2][:,1], tsoil, tcan, laiM, vodA, vodB, vodC, 10^5, 1e-6, 1e-5, 2, mask1);

	tb2_all = get_TB_2(sim_res2[2][:,1], tsoil, tcan, cm2, laiM, ores2_all[1:3]...);
	tb2_1 = get_TB_2(sim_res2[2][:,1], tsoil, tcan, cm2, laiM, ores2_1[1:3]...);
	tb2_6 = get_TB_2(sim_res2[2][:,1], tsoil, tcan, cm2, laiM, ores2_6[1:3]...);

	push!(tbx_all,tb2_all)
	push!(tbx_6,tb2_1)
	push!(tbx_1,tb2_6)
end

function cat_temp(x)
	return vcat(x[1],x[2])
end

tdiff_all = [sqrt.(mean((cat_temp(tb1)-cat_temp(x)) .^ 2)) for x in tbx_all];
tdiff_6 = [sqrt.(mean((cat_temp(tb1)-cat_temp(x))[vcat(mask6,mask6)] .^ 2)) for x in tbx_6];
tdiff_1 = [sqrt.(mean((cat_temp(tb1)-cat_temp(x))[vcat(mask1,mask1)] .^ 2)) for x in tbx_1];

factor_list4 = repeat(factor_list, outer=5);

significant_rmse_1yr = 1.3*sqrt(quantile.(Chisq.(365), 0.95)/(365) - 1);
significant_rmse_10yr = 1.3*sqrt(quantile.(Chisq.(365*10), 0.95)/(365*10) - 1);

figure()

subplot(1,3,1)
plot(reshape(factor_list4,(6,5)),reshape(tdiff_all,(6,5)),"o-",label=parnames)
#xlabel("Parameter multiplier")
#ylabel("Tb RMSE (K)")
xscale("log")
ylim(0,0.9)
title("Full diurnal")
legend()
xticks(factor_list,round.(factor_list,digits=2))
plot([1/3,3],significant_rmse_1yr*[1,1],"k--")
plot([1/3,3],significant_rmse_10yr*[1,1],"k--")

subplot(1,3,2)
plot(reshape(factor_list4,(6,5)),reshape(tdiff_1,(6,5)),"o-",label=parnames)
xscale("log")
ylim(0,0.9)
title("1 AM/PM")
xticks(factor_list,round.(factor_list,digits=2))
plot([1/3,3],significant_rmse_1yr*[1,1],"k--")
plot([1/3,3],significant_rmse_10yr*[1,1],"k--")

subplot(1,3,3)
plot(reshape(factor_list4,(6,5)),reshape(tdiff_6,(6,5)),"o-",label=parnames)
xscale("log")
ylim(0,0.9)
title("6 AM/PM")
xticks(factor_list,round.(factor_list,digits=2))
plot([1/3,3],significant_rmse_1yr*[1,1],"k--")
plot([1/3,3],significant_rmse_10yr*[1,1],"k--")

figure()
plot(factor_list4, tdiff_all,"o",label="All")
plot(factor_list4, tdiff_6,"o",label="6 AM/PM")
plot(factor_list4, tdiff_1,"o",label="1 AM/PM")
plot([1],[0],"ko")

significant_rmse_1yr = 1.3*sqrt(quantile.(Chisq.(365), 0.95)/(365) - 1);
plot([factor_list[1],factor_list[end]],[1,1]*significant_rmse_1yr,"k--")

significant_rmse_10yr = 1.3*sqrt(quantile.(Chisq.(365*10), 0.95)/(365*10) - 1);
plot([factor_list[1],factor_list[end]],[1,1]*significant_rmse_10yr,"k--")

xlabel("Parameter change")
ylabel("Tb RMSE")
xscale("log")

#npoint_list = 120:10:(360*10);
#chilist = [1.3*sqrt(quantile.(Chisq.(i), 0.95)/i - 1) for i in npoint_list];

ETdiff = [sqrt.(mean((sim_res1[1].ETmod - x[1].ETmod) .^ 2)) for x in sim_list]*18/1000*60*60*24;
GPPdiff = [sqrt.(mean((sim_res1[1].GPP - x[1].GPP) .^ 2)) for x in sim_list];
LWPdiff = [sqrt.(mean((cm1 - x) .^ 2)) for x in cmx];
LWPcor = [cor(cm1,x)[1,1] for x in cmx];

stdET = std(sim_res1[1].ETmod)*18/1000*60*60*24;
stdGPP = std(sim_res1[1].GPP);
stdLWP = std(cm1);

figure()
plot(factor_list4, 1 .- ETdiff ./ stdET, "o", label="ET")
plot(factor_list4, 1 .- GPPdiff ./ stdGPP, "o", label="GPP")
plot(factor_list4, 1 .- LWPdiff ./ stdLWP, "o", label="LWP")
plot([1],[1],"ko")
xlabel("Parameter change")
ylabel("R2")
legend()
title("Hourly")
xscale("log")

ETdiffD = [sqrt.(mean(get_daily(sim_res1[1].ETmod - x[1].ETmod, 24) .^ 2)) for x in sim_list]*18/1000*60*60*24;
GPPdiffD = [sqrt.(mean(get_daily(sim_res1[1].GPP - x[1].GPP, 24) .^ 2)) for x in sim_list];
LWPdiffD = [sqrt.(mean(get_daily(cm1 - x, 24) .^ 2)) for x in cmx];

DstdET = std(get_daily(sim_res1[1].ETmod,24))*18/1000*60*60*24;
DstdGPP = std(get_daily(sim_res1[1].GPP,24));
DstdLWP = std(get_daily(cm1,24));
DstdSMC = std(get_daily(sim_res1[1].ColumnSMC,24));

figure()
plot(factor_list4, 1 .- ETdiffD ./ DstdET, "o", label="ET")
plot(factor_list4, 1 .- GPPdiffD ./ DstdGPP, "o", label="GPP")
plot(factor_list4, 1 .- LWPdiffD ./ DstdLWP, "o", label="LWP")
plot([1],[1],"ko")
xlabel("Parameter change")
ylabel("R2")
legend()
title("Daily")
xscale("log")

figure()
plot(tdiff_all, LWPdiffD,"o",label="Full diurnal")
#plot(tdiff_1, LWPdiffD,"o",label="1 AM/PM")
#plot(tdiff_6, LWPdiffD,"o",label="6 AM/PM")
plot([0],[0],"ko")
xlabel("RMSE of Tb (K)")
ylabel("RMSE of LWP (MPa)")
plot([1,1]*significant_rmse_10yr,[0,maximum(LWPdiffD)],"k--")
plot([1,1]*significant_rmse_1yr,[0,maximum(LWPdiffD)],"k--")
text(significant_rmse_10yr+0.01,maximum(LWPdiffD)*0.75,"Significant\nover 10 years",rotation="vertical")
text(significant_rmse_1yr+0.01,maximum(LWPdiffD)*0.75,"Significant\nover 1 year",rotation="vertical")
#legend()

figure()
plot(tdiff_all, ETdiffD,"o",label="Full diurnal")
#plot(tdiff_1, ETdiffD,"o",label="1 AM/PM")
#plot(tdiff_6, ETdiffD,"o",label="6 AM/PM")
plot([0],[0],"ko")
xlabel("RMSE of Tb (K)")
ylabel("RMSE of ET (mm/day)")
plot([1,1]*significant_rmse_10yr,[0,maximum(ETdiffD)],"k--")
plot([1,1]*significant_rmse_1yr,[0,maximum(ETdiffD)],"k--")
text(significant_rmse_10yr+0.01,maximum(ETdiffD)*0.75,"Significant\nover 10 years",rotation="vertical")
text(significant_rmse_1yr+0.01,maximum(ETdiffD)*0.75,"Significant\nover 1 year",rotation="vertical")
#legend()

SMCdiffD = [sqrt.(mean(get_daily(sim_res1[1].ColumnSMC - x[1].ColumnSMC, 24) .^ 2)) for x in sim_list];

sSMCdiffD = [sqrt.(mean(get_daily(sim_res1[2][:,1] - x[2][:,1], 24) .^ 2)) for x in sim_list];

myY = SMCdiffD;

figure()
plot(tdiff_all, myY,"o",label="Full diurnal")
#plot(tdiff_1, myY,"o",label="1 AM/PM")
#plot(tdiff_6, myY,"o",label="6 AM/PM")
plot([0],[0],"ko")
xlabel("RMSE of Tb (K)")
ylabel("RMSE of Column SMC")
plot([1,1]*significant_rmse_10yr,[0,maximum(myY)],"k--")
plot([1,1]*significant_rmse_1yr,[0,maximum(myY)],"k--")
text(significant_rmse_10yr+0.01,maximum(myY)*0.75,"Significant\nover 10 years",rotation="vertical")
text(significant_rmse_1yr+0.01,maximum(myY)*0.75,"Significant\nover 1 year",rotation="vertical")
#legend()

figure()
plot(tdiff_all,tdiff_1,"o",label="1 AM/PM")
plot(tdiff_all,tdiff_6,"o",label="6 AM/PM")
plot([0,1],[0,1],"k--")
xlabel("Full diurnal RMSE (K)")
ylabel("Twice-daily RMSE (K)")
legend()

#ETdiff_obs = [sqrt.(mean((sim_res1[1].ETmod - sim_res1[1].LE/44200) .^ 2)) for x in sim_list]*18/1000*60*60;

figure()

subplot(1,3,1)
plot(reshape(factor_list4,(6,5)),reshape(LWPdiffD,(6,5)),"o-",label=parnames)
#xlabel("Parameter multiplier")
#ylabel("Tb RMSE (K)")
xscale("log")
#ylim(0,0.9)
title("LWP (MPa)")
legend()
xticks(factor_list,round.(factor_list,digits=2))

subplot(1,3,2)
plot(reshape(factor_list4,(6,5)),reshape(ETdiffD,(6,5)),"o-",label=parnames)
xscale("log")
title("ET (mm/day)")
xticks(factor_list,round.(factor_list,digits=2))

subplot(1,3,3)
plot(reshape(factor_list4,(6,5)),reshape(SMCdiffD,(6,5)),"o-",label=parnames)
xscale("log")
title("Column SMC")
xticks(factor_list,round.(factor_list,digits=2))

function diurnal_pattern(x)
	z1 = get_diurnal(x,24);
	return (z1 .- mean(z1)) / std(z1);
end
shapeLWPdiff = [sqrt.(mean((diurnal_pattern(cm1) - diurnal_pattern(x)) .^ 2)) for x in cmx];

all_diurnal = [diurnal_pattern(x) for x in cmx];

figure()
plot(hcat(all_diurnal...),alpha=0.25)
plot(diurnal_pattern(cm1),"k")
xlabel("Time of day")
ylabel("Normalized LWP")

figure()

subplot(1,2,1)
plot(reshape(factor_list4,(6,5)),reshape(LWPdiffD,(6,5)),"o-",label=parnames)
#xlabel("Parameter multiplier")
#ylabel("Tb RMSE (K)")
xscale("log")
#ylim(0,0.9)
title("Daily LWP (MPa)")
legend()
xticks(factor_list,round.(factor_list,digits=2))

subplot(1,2,2)
plot(reshape(factor_list4,(6,5)),reshape(shapeLWPdiff,(6,5)),"o-",label=parnames)
#xlabel("Parameter multiplier")
#ylabel("Tb RMSE (K)")
xscale("log")
#ylim(0,0.9)
title("LWP diurnal shape")
xticks(factor_list,round.(factor_list,digits=2))

tbx_H = []

for i in 1:length(sim_list)
	println(i)
	sim_res2 = sim_list[i];
	cm2 = cmx[i]

	tb2_all = get_TB_2(sim_res2[2][:,1], tsoil, tcan, cm2, laiM, vodA, vodB, vodC);

	push!(tbx_H,tb2_all)
end

tdiff_allH = [sqrt.(mean((cat_temp(tb1)-cat_temp(x)) .^ 2)) for x in tbx_H];
tdiff_6H = [sqrt.(mean((cat_temp(tb1)-cat_temp(x))[vcat(mask6,mask6)] .^ 2)) for x in tbx_H];
tdiff_1H = [sqrt.(mean((cat_temp(tb1)-cat_temp(x))[vcat(mask1,mask1)] .^ 2)) for x in tbx_H];

figure()
plot(tdiff_all, shapeLWPdiff,"o",label="Full diurnal")
#plot(tdiff_1, shapeLWPdiff, "o",label="1 AM/PM")
#plot(tdiff_6, shapeLWPdiff, "o",label="6 AM/PM")
plot([0],[0],"ko")
xlabel("RMSE of Tb (K)")
ylabel("RMSE of LWP diurnal shape")
plot([1,1]*significant_rmse_10yr,[0,maximum(myY)],"k--")
plot([1,1]*significant_rmse_1yr,[0,maximum(myY)],"k--")
text(significant_rmse_10yr+0.01,maximum(myY)*0.75,"Significant\nover 10 years",rotation="vertical")
text(significant_rmse_1yr+0.01,maximum(myY)*0.75,"Significant\nover 1 year",rotation="vertical")
#legend()

figure()
plot(tdiff_all, 1 .- LWPdiffD / DstdLWP,"o",label="LWP")
plot(tdiff_all, 1 .- ETdiffD / DstdET,"o",label="ET")
plot(tdiff_all, 1 .- SMCdiffD / DstdSMC,"o",label="SMC")
plot([0],[1],"ko")
xlabel("RMSE of Tb (K)")
ylabel("R2 of model state")
plot([1,1]*significant_rmse_10yr,[-2,1],"k--")
plot([1,1]*significant_rmse_1yr,[-2,1],"k--")
plot([0,maximum(tdiff_all)], [0,0],"k")
text(significant_rmse_10yr+0.01,-1,"Significant\nover 10 years",rotation="vertical")
text(significant_rmse_1yr+0.01,-1,"Significant\nover 1 year",rotation="vertical")
legend()

ext: jl


text:
macro mmatrace()
    esc(quote
        if tracing
            dt = Dict()
            if model.extended_trace
                dt["x"] = copy(x)
                dt["g(x)"] = copy(∇f_x)
                dt["λ"] = copy(results.minimizer)
            end
            update!(tr,
                    iter,
                    f_x,
                    gr_residual,
                    dt,
                    model.store_trace,
                    model.show_trace)
        end
    end)
end


text:
# StringFields
# ============
#
# UTF8-encoded mutable string type.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module StringFields

export StringField

import BufferedStreams

"""
A simplistic mutable, utf8 encoded string.
"""
type StringField <: AbstractString
    data::Vector{UInt8}

    # position of used data in `data` in bytes (not characters)
    part::UnitRange{Int}
end

function StringField()
    return StringField(UInt8[], 1:0)
end

function StringField(data::Vector{UInt8})
    return StringField(data, 1:length(data))
end

# From base unicode/utf8.jl
function Base.endof(s::StringField)
    d = s.data
    i = s.part.stop
    i == 0 && return i
    while Base.is_valid_continuation(d[i])
        i -= 1
    end
    @assert i >= s.part.start
    return i - s.part.start + 1
end

# From base unicode/utf8.jl
function Base.isvalid(s::StringField, i::Integer)
    return (1 <= i <= endof(s.data)) &&
        !Base.is_valid_continuation(s.data[s.part.start + i - 1])
end

# From base unicode/utf8.jl
function Base.next(s::StringField, i::Int)
    d = s.data
    b = d[s.part.start + i - 1]
    if Base.is_valid_continuation(b)
        throw(UnicodeError(Base.UTF_ERR_INVALID_INDEX, i, d[s.part.start + i - 1]))
    end

    trailing = Base.utf8_trailing[b+1]
    if length(s.part) < i + trailing
        return '\ufffd', i+1
    end
    c::UInt32 = 0
    for j = 1:trailing+1
        c <<= 6
        c += d[s.part.start + i - 1]
        i += 1
    end
    c -= Base.utf8_offset[trailing+1]
    return Char(c), i
end

function Base.copy!(field::StringField, data::Vector{UInt8},
               start::Integer, stop::Integer)
    if length(field.data) < length(data)
        resize!(field.data, length(data))
    end
    n = stop - start + 1
    copy!(field.data, 1, data, start, n)
    field.part = 1:n
    return n
end

function Base.empty!(field::StringField)
    field.part = 1:0
    return field
end

function Base.isempty(field::StringField)
    return field.part.start > field.part.stop
end

function Base.convert(::Type{StringField}, str::String)
    return StringField(copy(str.data), 1:length(str.data))
end

function Base.convert(::Type{String}, field::StringField)
    return String(field.data[field.part])
end

function Base.convert(::Type{AbstractString}, field::StringField)
    return convert(String, field)
end

function Base.write(io::IO, field::StringField)
    write(io, convert(String, field))
end

function Base.show(io::IO, field::StringField)
    show(io, convert(String, field))
end

function Base.copy(field::StringField)
    data = field.data[field.part]
    return StringField(data, 1:length(data))
end

# From Base.hash over strings in hashing2.jl
function Base.hash(field::StringField, h::UInt64)
    h += Base.memhash_seed
    return ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32),
                 pointer(field.data, field.part.start),
                 length(field.part), h % UInt32) + h
end

function Base.:(==)(a::StringField, b::StringField)
    return length(a) == length(b) &&
        ccall(:memcmp, Cint, (Ptr{Void}, Ptr{Void}, Csize_t),
              pointer(a.data, a.part.start), pointer(b.data, b.part.start), length(a)) == 0
end

function Base.:(==)(a::StringField, b::BufferedStreams.BufferedOutputStream)
    if a === b
        return true
    elseif length(a) == length(b)
        return ccall(:memcmp, Cint, (Ptr{Void}, Ptr{Void}, Csize_t),
                     pointer(a.data, a.part.start), b.buffer, length(a.data)) == 0
    else
        return false
    end
end

end # module StringFields


text:
##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
# M
##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
mutable struct highDimMat{T, S<:AbstractVecOrMat}
    M::S
end

# constructor for highDimMat M,
function highDimMat(M::AbstractVecOrMat{T}) where {T}
    if size(M,1) >= size(M,2)
        @warn "n >= p in high dimensional matrix"
    end
    highDimMat{T, typeof(M)}(M)
end

"""
# constructor for highDimMat M,
function highDimMat(M::AbstractMatrix)
    if size(M,1) >= size(M,2)
        @warn "n >= p in high dimensional matrix"
    end
    return highDimMat{eltype(M), typeof(M)}(M)
end
"""

Base.copyto!(A::highDimMat{T}, src::AbstractVecOrMat{T}) where {T} = copyto!(A.M, src)

Base.eltype(::highDimMat{T}) where {T} = T

Base.getindex(A::highDimMat{T}, i::Int, j::Int) where {T} = getindex(A.M, i, j)

Base.size(A::highDimMat{T}) where {T} = size(A.M)

Base.size(A::highDimMat{T}, i::Integer) where {T} = size(A.M, i)
##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
# X
##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
mutable struct XMat{T, S<:AbstractVecOrMat}
    X::S
end

# constructor for XMat X,
function XMat(X::AbstractVecOrMat{T}) where {T}
    if rank(X) < size(X,2)
        @warn "fixed effect matrix is not of full rank"
    end

    if size(X,1) < size(X,2)
        @warn "n < p in covariate matrix X"
    end
    XMat{T, typeof(X)}(X)
end

Base.copyto!(A::XMat{T}, src::AbstractVecOrMat{T}) where {T} = copyto!(A.X, src)

Base.eltype(::XMat{T}) where {T} = T

LinearAlgebra.rank(A::XMat{T}) where {T} = rank(A.X)

isfullrank(A::XMat{T}) where {T} = rank(A) == size(A.X,2)

Base.getindex(A::XMat{T}, i::Int, j::Int) where {T} = getindex(A.X, i, j)

#Base.adjoint(A::XMat{T}) = Adjoint(A)

Base.size(A::XMat{T}) where {T} = size(A.X)

Base.size(A::XMat{T}, i::Integer) where {T} = size(A.X, i)

## define new infix binary operator?

function Base.:*(A::highDimMat{T}, B::XMat{T}) where{T}
    A.M*B.X
end

#*(A::highDimMat{T, AbstractMatrix{T}}, B::XMat{T, AbstractMatrix{T}}) where{T} = A.M*B.X

*(A::highDimMat{Int64, Matrix{Int64}}, B::XMat{Int64, Matrix{Int64}}) where{T} = A.M*B.X

##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
# Z
##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
abstract type AbstractReMat{T} <: AbstractMatrix{T} end

mutable struct ReMat{T,S} <: AbstractMatrix{T}
    #trm # the grouping factor as a `StatsModels.CategoricalTerm`   ##????
    Z::Matrix{T}
end

# constructor for XMat X,
function ReMat(Z::AbstractMatrix{T}) where {T}
    ReMat{T, typeof(Z)}(Z)
end

LinearAlgebra.rank(A::ReMat) = rank(A.Z)

isfullrank(A::ReMat) = rank(A) == size(A.Z,2)

Base.getindex(A::ReMat, i::Int, j::Int) = getindex(A.Z, i, j)

Base.size(A::ReMat) = size(A.Z)

##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
# highDimMixedModel
##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
"""
    highDimMixedModel
High dim mixed-effects model representation
## Fields
* `formula`: the formula for the model
* `optsum`: an [`OptSummary`](@ref) object
## Properties
* `θ` or `theta`: the covariance parameter vector used to form λ
* `β` or `beta`: the fixed-effects coefficient vector
* `λ` or `lambda`: a vector of lower triangular matrices repeated on the diagonal blocks of `Λ`
* `σ` or `sigma`: current value of the standard deviation of the per-observation noise
* `b`: random effects on the original scale, as a vector of matrices
* `u`: random effects on the orthogonal scale, as a vector of matrices
* `lowerbd`: lower bounds on the elements of θ
* `X`: the fixed-effects model matrix
* `y`: the response vector
"""
mutable struct highDimMixedModel{T<:AbstractFloat}  <: MixedModel{T}
    formula::FormulaTerm
    M::highDimMat{T}
    X::XMat{T}
    Z::ReMat{T}
    y::Vector{T}
    optsum::Union{OptSummary{T}, Nothing}
end

""" (1)
    highDimMixedModel
dealing with the number of columns in M and X scenario, assume intercept alwalys exists (@formula: y ~ 1 + ... )
# @ numOfHDM: number of HDM variables
# @ numOfHDM: number of X variables

This function assumes the variables in dataframe is arranged in order HDM, X , Z(ReMat).
 The issue here is matrix of random effect is without intercept (set to zero constrain)

"""
function highDimMixedModel(
    f::FormulaTerm,
    df::DataFrame,
    contrasts::Dict{Symbol, UnionAll},
    numOfHDM::Int64,
    numOfXMat::Int64
) where{T}
    for i in 1:size(df,2)
        if(isa(df[1,i], Number))
            df[!,i] = convert(Vector{Float64},df[:,i])
        end
    end
    sch = schema(df,contrasts)
    form = apply_schema(f, sch)
    y, pred = modelcols(form, df);
    MXZ = pred[:,2:size(pred,2)]  ## get rid of intercept
    M = highDimMat(MXZ[:,1:numOfHDM])
    intercept = pred[:,1]
    X = XMat(hcat(reshape(intercept, size(intercept,1),1), MXZ[:, (numOfHDM + 1):(numOfHDM + numOfXMat)]))  ## concatenate intercept with X
    Z = ReMat(MXZ[:, (numOfHDM + numOfXMat + 1):size(MXZ,2)])

    #return highDimMixedModel{T<:AbstractFloat}(form, M, X, Z, y)
    return highDimMixedModel{Float64}(form, M, X, Z, y, nothing)
end

""" (2)
    highDimMixedModel(...)
   Private constructor for a highDimMixedModel.
To construct a model, you only need the formula (`f`), data(`df`), contrasts and indices of M,X,Z

# the ID given here must be coincide with the ID in formula
# The issue here is matrix of random effect is without intercept (set to zero constrain)

# @ idOfHDM: id of variables corresponding to matrix M, NOTE: it's the ID in hte formula! # [1,2]
# @ idOfXMat: same as before
# @ idOfReMat: same as before
"""

""" not useful right now
private help function for select indices for remaining random effect matrix
    preTerms is Vector(1:length(form.rhs.terms))

_filt(x, preTerms) = !(x in preTerms)

preTerms = Vector(1:length(form.rhs.terms))

"""

function highDimMixedModel(
    f::FormulaTerm,
    df::DataFrame,
    contrasts::Dict{Symbol, UnionAll},
    idOfHDM::Union{Int,AbstractArray{Int64,1}}, # [1,2]
    idOfXMat::Union{Int,AbstractArray{Int64,1}},
    idOfReMat::Union{Int,AbstractArray{Int64,1}}
) where {T}
    for i in 1:size(df,2)
        if(isa(df[1,i], Number))
            df[!,i] = convert(Vector{Float64},df[:,i])
        end
    end
    sch = schema(df, contrasts)
    form = apply_schema(f, sch)
    y, pred = modelcols(form, df);
    terms = form.rhs.terms
    M = highDimMat(modelmatrix(terms[idOfHDM],df))
    X = XMat(modelmatrix(terms[idOfXMat],df))

    #preTerms = Vector(1:length(form.rhs.terms))
    #idOfReMat = preTerms[preTerms .∉ Ref(vcat(idOfHDM, idOfXmat))]
    Z = ReMat(modelmatrix(terms[idOfReMat],df))

    return highDimMixedModel{Float64}(form, M, X, Z, y, nothing)

end

""" (3)
    highDimMixedModel(...)
   Private constructor for a highDimMixedModel.
To construct a model, you only need the formula (`f`), data(`df`), contrasts and names of columns of M,X,Z

!!! note
    This method requires the order of terms in formula corresponds to the order of names in df
"""
function highDimMixedModel(
    f::FormulaTerm,
    df::DataFrame,
    contrasts::Dict{Symbol, UnionAll},
    nameOfHDM::Union{AbstractString,AbstractArray{<:AbstractString,1}}, # ["a","b"]
    nameOfXMat::Union{AbstractString,AbstractArray{<:AbstractString,1}},
    nameOfReMat::Union{AbstractString,AbstractArray{<:AbstractString,1}}
) where {T}
    """
    for i in 1:size(df,2)
        if(eltype(df[:,i]) == Int64)
            df[!,i] = convert(Vector{Float64},df[:,i])
        end
    end
    """
    for i in 1:size(df,2)
        if(isa(df[1,i], Number))
            df[!,i] = convert(Vector{Float64},df[:,i])
        end
    end
    sch = schema(df, contrasts)
    form = apply_schema(f, sch)
    y, pred = modelcols(form, df);
    terms = form.rhs.terms

    if(typeof(form.rhs.terms[1]) <: InterceptTerm)
        ### This line is due to 0 intercept term in formula
        terms = terms[2:length(terms)]
    end

    ## get id by names
    namesOfVar = names(df)[2:length(names(df))] ## extract variable names
    if typeof(nameOfHDM) == String idOfHDM = findall(x -> x == nameOfHDM, namesOfVar)
    else idOfHDM = findall(x -> x in nameOfHDM, namesOfVar) ; end

    if typeof(nameOfXMat) == String idOfXMat = findall(x -> x == nameOfXMat, namesOfVar)
    else idOfXMat = findall(x -> x in nameOfXMat, namesOfVar) ; end

    if typeof(nameOfReMat) == String idOfReMat = findall(x -> x == nameOfReMat, namesOfVar)
    else idOfReMat = findall(x -> x in nameOfReMat, namesOfVar) ; end

    M = highDimMat(modelmatrix(terms[idOfHDM],df))
    X = XMat(modelmatrix(terms[idOfXMat],df))

    #preTerms = Vector(1:length(form.rhs.terms))
    #idOfReMat = preTerms[preTerms .∉ Ref(vcat(idOfHDM, idOfXmat))]
    Z = ReMat(modelmatrix(terms[idOfReMat],df))

    return highDimMixedModel{Float64}(form, M, X, Z, y, nothing)

end

""" (4)
    highDimMixedModel(...)
   Private constructor for a highDimMixedModel.
To construct a model, you only need the formula (`f`), data(`df`), names of columns of M,X,
this function utilize the apply_schema() from MixedModel, so it can parse the y ~ 0 + a + b + (1|c) as mixed model
The modelcols/modelmatrix in StatsModels.jl cannot parse (1|c), need more info to distinguish between matrix M and X

!!! note
    This method requires the order of terms in formula corresponds to the order of names in df
"""
function highDimMixedModel(
    f::FormulaTerm,
    df::DataFrame,
    numOfHDM::Int64,
    #numOfXMat::Int64
) where {T}
    """
    for i in 1:size(df,2)
        if(eltype(df[:,i]) == Int64)
            df[!,i] = convert(Vector{Float64},df[:,i])
        end
    end
    """
    for i in 1:size(df,2)
        if(isa(df[1,i], Number))
            df[!,i] = convert(Vector{Float64},df[:,i])
        end
    end

    ## below utilize apply_schema() from MixedModel, the reason to keep above is we still need variable names
    ## to identify M and X
    ## use apply_schema to get random effect matrix
    sch = schema(df)
    form = apply_schema(f, sch, highDimMixedModel)
    #preTerms = Vector(1:length(form.rhs.terms))
    #idOfReMat = preTerms[preTerms .∉ Ref(vcat(idOfHDM, idOfXmat))]
    y, pre = modelcols(form,df)
    FixedEffectMatrix,ZMatrix = pre  #modelmatrix(form,df)
    M = highDimMat(FixedEffectMatrix[:,1:numOfHDM])
    X = XMat(FixedEffectMatrix[:, (numOfHDM + 1):size(FixedEffectMatrix,2)])
    Z = ReMat(ZMatrix)

    return highDimMixedModel{Float64}(form, M, X, Z, y, nothing)

end

""" (5)
    highDimMixedModel(...)
   Private constructor for a highDimMixedModel.
   This function is a generalization of constructor when we have high dimensional data, in this setting, we won't need formula(since it's too long),
        we will need id of each random variable, then we will directly use id of variable of M and X make the matrix M and X,
        use tricks in (4) to make Z

# the ID given here must be coincide with the ID in formula
# The issue in (1) is matrix of random effect is without intercept (set to zero constrain), we try to tackle it here using tricks in (4)

# @ nameOfy: the variable name of response for constructing formula # ["a", "b"]
# @ idOfReMat: same as before

# @ idOfHDM: id of variables corresponding to matrix M, NOTE: it's the ID in the df! # [1,2]
# @ idOfXMat: same as before

"""

function highDimMixedModel(
    df::DataFrame,
    nameOfy::Union{AbstractString,AbstractArray{<:AbstractString,1}}, # ["a", "b"]
    nameOfReMat::Union{AbstractString,AbstractArray{<:AbstractString,1}},
    idOfHDM::Union{Int,AbstractArray{Int64,1}}, # [1,2]
    idOfXMat::Union{Int,AbstractArray{Int64,1}},
    #contrasts = Dict(term(nameOfReMat) => CategoricalTerm)
) where {T}
    for i in 1:size(df,2)
        if(isa(df[1,i], Number))
            df[!,i] = convert(Vector{Float64},df[:,i])
        end
    end

    M = highDimMat(Matrix(df[:,idOfHDM]))
    X = XMat(Matrix(df[:,idOfXMat]))

    # The function doesn't need formula to construct M and X, so the formula constucted here is only for extracting Z,
    # doesn't matter we have intercept or not
    # need to transfer nameOfy and nameOfReMat to term due to error: KeyError: key nameOfReMat not found
    # https://juliastats.org/StatsModels.jl/stable/api/#StatsModels.term
    # need to transform nameOfReMat from string to categorical ?

    #transform!(df, :trt => categorical, renamecols=false) # this is need to avoid error in term(nameOfReMat)
    f = @formula(term(nameOfy) ~ (1|term(nameOfReMat)))
    print(nameOfy)
    #sch = schema(df, contrasts)
    sch = schema(df)
    form = apply_schema(f, sch, highDimMixedModel)

    y, pre = modelcols(form,df)
    FixedEffectMatrix,ZMatrix = pre
    Z = ReMat(ZMatrix)

    return highDimMixedModel{Float64}(form, M, X, Z, y, nothing)

end

##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
# fit
##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
"""
function fit(
    ::Type{highDimMixedModel},
    f::FormulaTerm,
    df::DataFrame;
    contrasts,
    progress,
    REML,
    idOfHDM::Union{Int,AbstractArray{Int64,1}}, # [1,2]
    idOfXMat::Union{Int,AbstractArray{Int64,1}},
    idOfReMat::Union{Int,AbstractArray{Int64,1}}
) return fit!(highDimMixedModel(f, df, contrasts, idOfHDM, idOfXMat, idOfReMat), progress, REML)
end
"""
function fit(
    HMM::highDimMixedModel{T};
    verbose::Bool=true,
    REML::Bool=true,
    alg = :LN_COBYLA,
) where{T}
    return fit!(HMM, verbose = verbose, REML = REML, alg = alg)
end

function fit!(HMM::highDimMixedModel{T}; verbose::Bool=true, REML::Bool=true, alg = :LN_COBYLA) where {T}
    n = size(HMM.M, 1)
    A = hcat(HMM.M.M, HMM.X.X)
    P = I - A*inv(transpose(A)*A)*transpose(A)
    u,s,v = svd(P)
    r = size(HMM.M,2) + size(HMM.X,2)  # simplify: assume fixed effect full rank
    #C = randn((n-r),n)
    C = transpose(u[:,1:(n-r)])
    K = C*P
    Z = HMM.Z.Z
    y = HMM.y
    # C can be any full rank matrix with size n,r, e.g. randn(n,r)

    ## add optsum
    # init para
    sigma = [2200.0,14.0]
    lbd = [0.0; 0.0]
    optsum = OptSummary(Float64.(sigma), lbd, alg; ftol_rel=T(1.0e-15), ftol_abs=T(1.0e-8), xtol_rel = 1e-10)
    optsum.REML = REML

    ## init opt based on optsum
    opt = Opt(optsum)

    function negLogLik(sigma::Vector{Float64}, g::Vector{Float64})
        n = length(y)
        Sigma = sigma[1]*Z*transpose(Z) + sigma[2]*diagm(ones(n))
        negLog = 1/2*log(det(K*Sigma*transpose(K))) + 1/2*transpose(y)*transpose(K)*inv(K*Sigma*transpose(K))*K*y
        #println("OPT: parameter $(sigma) || objective eval $(negLog)")
        @show sigma
        @show negLog

        return negLog
    end

    println("The initial object value is $(negLogLik(Float64.(sigma), [1.0,1.0]))")

    opt.min_objective = negLogLik
    optsum.finitial = negLogLik(optsum.initial, [1.0,1.0])  # the second field not useful right now

    if verbose println("OPTBL: starting point $(sigma)") ; end    # to stdout

    (optf,optx,ret) = optimize(opt, sigma)

    Sigma = optx[1]*Z*transpose(Z) + optx[2]*diagm(ones(n))
    beta = inv(transpose(A)*inv(Sigma)*A)*transpose(A)*inv(Sigma)*y
    betaM = beta[1:size(HMM.M,2)]
    betaX = beta[(size(HMM.M,2) + 1): size(beta)[1]]

    optsum.feval = opt.numevals
    optsum.final = optx
    optsum.fmin = optf
    optsum.returnvalue = ret

    if verbose println("got $(round(optf, digits=5)) at $(round.(optx, digits=5)) after $(opt.numevals) iterations (returned $(ret))") ; end

    HMM.optsum = optsum

    return optx, betaM, betaX, optsum

end

##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
# Base.getproperty
##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============

function Base.getproperty(HMM::highDimMixedModel{T}, s::Symbol) where{T}
    if s == :σ || s == :sigma
        getσ(HMM)
    elseif s == :β || s ==:beta
        getβ(HMM)
    elseif s == :βM || s == :betaM
        getβM(HMM)
    elseif s == :βX || s == :betaX
        getβX(HMM)
    else
        getfield(HMM,s)
    end
end

function getσ(HMM::highDimMixedModel{T}) where{T}
    return HMM.optsum.final
end

function getβ(HMM::highDimMixedModel{T}) where{T}
    n = size(HMM.M, 1)
    A = hcat(HMM.M.M, HMM.X.X)
    Z = HMM.Z.Z
    y = HMM.y
    optx = HMM.sigma
    Sigma = optx[1]*Z*transpose(Z) + optx[2]*diagm(ones(n))
    beta = inv(transpose(A)*inv(Sigma)*A)*transpose(A)*inv(Sigma)*y
    return beta
end

function getβM(HMM::highDimMixedModel{T}) where{T}
    beta = getβ(HMM)[1:size(HMM.M,2)]
    return beta[1:size(HMM.M,2)]
end

function getβX(HMM::highDimMixedModel{T}) where{T}
    beta = getβ(HMM)[1:size(HMM.M,2)]
    return beta[(size(HMM.M,2) + 1): size(beta)[1]]
end

##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============
# Base.show
##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============##==============

function Base.show(io::IO, ::MIME"text/plain", HMM::highDimMixedModel)
    if isnothing(HMM.optsum)
        return
    end
    REML = HMM.optsum.REML
    println(io, "Linear mixed model fit by ", REML ? "REML" : "not implemented yet")
    println(io, " ", HMM.formula)
    obj = HMM.optsum.fmin
    if REML
        println(io, " REML criterion at convergence: ", obj)
    else
        println("Only implemented REML at this time")
    end

    println("σ_z : $(HMM.sigma[1])")
    println("σ : $(HMM.sigma[2])")
    println("β_M : $(HMM.betaM)")
    println("β_X : $(HMM.betaX)")

end

max_stars_repo_name: solislemuslab/HighDimMixedModels

max_issues_repo_name: solislemuslab/HighDimMixedModels

max_forks_repo_name: solislemuslab/HighDimMixedModels


text:
#! /usr/bin/julia

# Rosetta Code, Sorting algorithms/Cocktail sort

function coctailsort{T<:Real}(a::Array{T,1})
    b = copy(a)
    isordered = false
    lo = one(Int)
    hi = length(b)
    while !isordered && hi > lo
        isordered = true
        for i in lo+1:hi
            if b[i] < b[i-1]
                t = b[i]
                b[i] = b[i-1]
                b[i-1] = t
                isordered = false
            end
        end
        hi -= 1
        if isordered || hi <= lo
            break
        end
        for i in hi:-1:lo+1
            if b[i-1] > b[i]
                t = b[i]
                b[i] = b[i-1]
                b[i-1] = t
                isordered = false
            end
        end
        lo += 1
    end
    return b
end

a = [rand(-2^10:2^10) for i in 1:20]
println("Before Sort:")
println(a)
a = coctailsort(a)
println("\nAfter Sort:")
println(a)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

#     CurvedSurface{RealT<:Real}
#
# Contains the data needed to represent a curve with data points (x,y) as a Lagrange polynomial
# interpolant written in barycentric form at a given set of nodes.
struct CurvedSurface{RealT<:Real}
  nodes               ::Vector{RealT}
  barycentric_weights ::Vector{RealT}
  coordinates         ::Array{RealT, 2} #[nnodes, ndims]
end

# evalute the Gamma curve interpolant at a particular point s and return the (x,y) coordinate
function evaluate_at(s, boundary_curve::CurvedSurface)

   @unpack nodes, barycentric_weights, coordinates = boundary_curve

   x_coordinate_at_s_on_boundary_curve = lagrange_interpolation(s, nodes, view(coordinates, :, 1),
                                                                barycentric_weights)
   y_coordinate_at_s_on_boundary_curve = lagrange_interpolation(s, nodes, view(coordinates, :, 2),
                                                                barycentric_weights)

   return x_coordinate_at_s_on_boundary_curve, y_coordinate_at_s_on_boundary_curve
end

# evalute the derivative of a Gamma curve interpolant at a particular point s
# and return the (x,y) coordinate
function derivative_at(s, boundary_curve::CurvedSurface)

   @unpack nodes, barycentric_weights, coordinates = boundary_curve

   x_coordinate_at_s_on_boundary_curve_prime = lagrange_interpolation_derivative(s, nodes,
                                                                                view(coordinates, :, 1),
                                                                                barycentric_weights)
   y_coordinate_at_s_on_boundary_curve_prime = lagrange_interpolation_derivative(s, nodes,
                                                                                view(coordinates, :, 2),
                                                                                barycentric_weights)
   return x_coordinate_at_s_on_boundary_curve_prime, y_coordinate_at_s_on_boundary_curve_prime
end

# Chebysehv-Gauss-Lobatto nodes and weights for use with curved boundaries
function chebyshev_gauss_lobatto_nodes_weights(n_nodes::Integer)

  # Initialize output
  nodes   = zeros(n_nodes)
  weights = zeros(n_nodes)

  # Get polynomial degree for convenience
  N = n_nodes - 1

  for j in 1:n_nodes
    nodes[j]   = -cospi( (j-1) / N )
    weights[j] = pi / N
  end
  weights[1]   = 0.5 * weights[1]
  weights[end] = 0.5 * weights[end]

  return nodes, weights
end

# Calculate Lagrange interpolating polynomial of a function f(x) at a given point x for a given
# node distribution.
function lagrange_interpolation(x, nodes, fvals, wbary)
# Barycentric two formulation of Lagrange interpolant
  numerator   = zero(eltype(fvals))
  denominator = zero(eltype(fvals))

  for j in eachindex(nodes)
    if isapprox(x, nodes[j], rtol=eps(x))
      return fvals[j]
    end
    t            = wbary[j] / ( x - nodes[j] )
    numerator   += t * fvals[j]
    denominator += t
  end

  return numerator/denominator
end

# Calculate derivative of a Lagrange interpolating polynomial of a function f(x) at a given
# point x for a given node distribution.
function lagrange_interpolation_derivative(x, nodes, fvals, wbary)

  at_node   = false
  numerator = zero(eltype(fvals))
  i         = 0

  for j in eachindex(nodes)
    if isapprox(x, nodes[j])
      at_node     = true
      p           = fvals[j]
      denominator = -wbary[j]
      i           = j
    end
  end

  if at_node
    for j in eachindex(nodes)
      if j != i
        numerator += wbary[j] * ( p - fvals[j] ) / ( x - nodes[j] )
      end
    end
  else
    denominator = zero(eltype(fvals))
    p = lagrange_interpolation(x, nodes, fvals, wbary)
    for j in eachindex(nodes)
      t            = wbary[j] / (x - nodes[j])
      numerator   += t * ( p - fvals[j] ) / ( x - nodes[j] )
      denominator += t
    end
  end

  return numerator/denominator # p_prime
end

end # @muladd


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 7420c480-910c-11eb-2338-6bbd68005f0a
begin
    using Pkg
    Pkg.add("PlutoUI")
    Pkg.add("FreqTables")
    Pkg.add("RDatasets")
    Pkg.add("StatsBase")
    Pkg.add("StatsPlots")
    Pkg.add("Distributions")
    Pkg.add("Statistics")
    Pkg.add("LaTeXStrings")
	Pkg.add("VegaLite")
end

# ╔═╡ b0f0b9a4-938c-11eb-3929-331bdd555c4c
using DataFrames,JSON,HTTP,JSONTables,Dates,TableIO,CSV,VegaLite

# ╔═╡ 298ff14c-910d-11eb-1088-efa7ab843b4f
using PlutoUI,FreqTables,RDatasets,StatsBase,StatsPlots,Distributions,Statistics,LaTeXStrings,QuadGK

# ╔═╡ 4b759868-9494-11eb-1f4e-3d0ab484ec88

# ╔═╡ 131772d8-917f-11eb-1add-1f94cfb9bf63
begin
	Pkg.add("QuadGK")
	Pkg.add("JSON")
	Pkg.add("HTTP")
	Pkg.add("JSONTables")
	Pkg.add("CSV")
	Pkg.add("TableIO")
end

# ╔═╡ Cell order:
# ╠═7420c480-910c-11eb-2338-6bbd68005f0a
# ╠═4b759868-9494-11eb-1f4e-3d0ab484ec88
# ╠═131772d8-917f-11eb-1add-1f94cfb9bf63
# ╠═b0f0b9a4-938c-11eb-3929-331bdd555c4c
# ╠═298ff14c-910d-11eb-1088-efa7ab843b4f


text:
using Printf

simulation = nothing
selectionIndex = -1
simulationActive = true

include("controls.jl")

function setupSimulation(simulation)
  electron = Electron(Vector3(0.0, 3.0, 2.0), Vector3(0.008, 0.0, -0.003))
  proton = Proton(Vector3(0.0, 0.0, 2.0), Vector3(-0.0008, 0.0, 0.0003))
  neutron = Neutron(Vector3(0.0, 0.0, 5.0), Vector3(0.0, 0.0, 0.0))

  addParticle!(simulation, electron)
  addParticle!(simulation, proton)
  addParticle!(simulation, neutron)
end

bottomMenu = nothing
entityLabel = nothing
timeLabel = nothing
pauseButton = nothing

inspectorMenu = nothing
inspectorLabel = nothing
kineticEnergyLabel = nothing

function pause(button)
  if simulationActive == true
    global simulationActive = false
    setButtonText(pauseButton, "Resume")
  else
    global simulationActive = true
    setButtonText(pauseButton, "Pause")
  end
  return
end

function setupMenus()
  menus = []

  global bottomMenu = createMenu()

  setMenuOrientation(bottomMenu, Horizontal)
  setMenuPosition(bottomMenu, Vector2(0.0, 0.0))
  setMenuSize(bottomMenu, Vector2(1280.0, 20.0))
  setMenuColor(bottomMenu, Vector4(0.3, 0.3, 0.3, 1.0))

  global entityLabel = createLabel()
  setLabelZ(entityLabel, 1)
  setLabelSize(entityLabel, Vector2(100.0, 15.0))
  setLabelPosition(entityLabel, Vector2(0.0, 2.5))
  setLabelColor(entityLabel, Vector4(1.0, 1.0, 1.0, 1.0))
  setLabelTextHeight(entityLabel, 15)
  setLabelText(entityLabel, "Entities: N")

  addControlToMenu(bottomMenu, entityLabel)

  global timeLabel = createLabel()
  setLabelZ(timeLabel, 1)
  setLabelSize(timeLabel, Vector2(100.0, 15.0))
  setLabelPosition(timeLabel, Vector2(100.0, 2.5))
  setLabelColor(timeLabel, Vector4(1.0, 1.0, 1.0, 1.0))
  setLabelTextHeight(timeLabel, 15)
  setLabelText(timeLabel, "Time: 0")

  addControlToMenu(bottomMenu, timeLabel)

  push!(menus, bottomMenu)

  global pauseButton = createButton()
  setButtonZ(pauseButton, 1)
  setButtonSize(pauseButton, Vector2(100.0, 15.0))
  setButtonPosition(pauseButton, Vector2(200.0, 2.5))
  #setButtonColor(pauseButton, Vector4(0.0, 0.0, 0.0, 0.0))
  setButtonTextHeight(pauseButton, 15)
  setButtonText(pauseButton, "Pause")
  setButtonTextColor(pauseButton, Vector4(1.0, 1.0, 1.0, 1.0))
  setButtonAction(pauseButton, pause)

  addControlToMenu(bottomMenu, pauseButton)

  global inspectorMenu = createMenu()

  setMenuOrientation(inspectorMenu, Vertical)
  setMenuPosition(inspectorMenu, Vector2(20.0, 200.0))
  setMenuSize(inspectorMenu, Vector2(200.0, 400.0))
  setMenuColor(inspectorMenu, Vector4(0.3, 0.3, 0.3, 1.0))

  push!(menus, inspectorMenu)

  global inspectorLabel = createLabel()
  setLabelZ(inspectorLabel, 1)
  setLabelSize(inspectorLabel, Vector2(180.0, 15.0))
  setLabelPosition(inspectorLabel, Vector2(10.0, 380.0))
  setLabelColor(inspectorLabel, Vector4(1.0, 1.0, 1.0, 1.0))
  setLabelTextHeight(inspectorLabel, 15)
  setLabelText(inspectorLabel, "Inspector")

  addControlToMenu(inspectorMenu, inspectorLabel)

  global positionLabel = createLabel()
  setLabelZ(positionLabel, 1)
  setLabelSize(positionLabel, Vector2(180.0, 15.0))
  setLabelPosition(positionLabel, Vector2(10.0, 360.0))
  setLabelColor(positionLabel, Vector4(1.0, 1.0, 1.0, 1.0))
  setLabelTextHeight(positionLabel, 15)
  setLabelText(positionLabel, "Position: ( , , )")

  addControlToMenu(inspectorMenu, positionLabel)

  global kineticEnergyLabel = createLabel()
  setLabelZ(kineticEnergyLabel, 1)
  setLabelSize(kineticEnergyLabel, Vector2(180.0, 15.0))
  setLabelPosition(kineticEnergyLabel, Vector2(10.0, 340.0))
  setLabelColor(kineticEnergyLabel, Vector4(1.0, 1.0, 1.0, 1.0))
  setLabelTextHeight(kineticEnergyLabel, 15)
  setLabelText(kineticEnergyLabel, "Kinetic Energy: ")

  addControlToMenu(inspectorMenu, kineticEnergyLabel)

  global particleLabel = createLabel()
  setLabelZ(particleLabel, 1)
  setLabelSize(particleLabel, Vector2(180.0, 15.0))
  setLabelPosition(particleLabel, Vector2(10.0, 320.0))
  setLabelColor(particleLabel, Vector4(1.0, 1.0, 1.0, 1.0))
  setLabelTextHeight(particleLabel, 15)
  setLabelText(particleLabel, "Particle: ")

  addControlToMenu(inspectorMenu, particleLabel)

  return menus
end

function resize(renderer, width, height)
  setMenuSize(bottomMenu, Vector2(convert(Float64, width), 20.0))
end

# Main simulation loop
function simulate()
  # Create the simulation
  global simulation = newSimulation()

  # Create the rendering setup
  windowWidth = 1280
  windowHeight = 720
  windowTitle = "Simulation"

  Renderer.init()
  window = Renderer.createWindow(windowWidth, windowHeight, windowTitle)
  renderer = Renderer.createRenderer(window)
  Renderer.setResizeCallback(renderer, resize)
  camera = Renderer.getCamera(renderer)

  registerControls(renderer, camera)

  function esc_callback(renderer, key)
    Renderer.closeWindow(window)
  end

  Renderer.addKeyCallback(renderer, esc_callback, Renderer.Escape, Renderer.Press)

  menus = setupMenus()

  setupSimulation(simulation)

  step = 0
  while Renderer.windowCloseStatus(window) != true
    stepsperframe = 10
    for i = 1:stepsperframe
      if simulationActive
        simulationStep!(simulation)
      end
    end

    entityCount = length(simulation.objects)
    setLabelText(entityLabel, "Entities: $entityCount")
    time = simulation.time
    setLabelText(timeLabel, "Time: $time")

    if selectionIndex > 0
      inspectorReadout(simulation.objects[selectionIndex])
    end

    Renderer.render(renderer, simulation, menus)
  end

  Renderer.terminate(renderer)
end

function inspectorReadout(particle)
  setLabelText(positionLabel, Printf.@sprintf("Position: (%0.2f, %0.2f, %0.2f)", particle.position.x, particle.position.y, particle.position.z))
  setLabelText(kineticEnergyLabel, Printf.@sprintf("Kinetic Energy: %0.2f", kineticEnergy(particle) * 1000))
  setLabelText(particleLabel, Printf.@sprintf("Particle: %s", particleName(particle)))
end

# Contains a particle and the necessary rendering handle
mutable struct ParticleContainer
  particle::Particle
  renderObject::Ptr{Cvoid}
end

# Simulation object
"An object that encapsulates the state of a simulation"
mutable struct Simulation
  time::Float64
  objects::Array{Particle, 1}
end

"Creates a blank simulation"
function newSimulation()
  return Simulation(0, [])
end

"Runs a single step of a simulation"
function simulationStep!(simulation)
  simulation.time += 1
  forces = Array{Vector3}(undef, length(simulation.objects))
  for i in 1:length(simulation.objects)
    forces[i] = Vector3(0.0, 0.0, 0.0)
  end
  for i in 1:length(simulation.objects)
    for j in (i + 1):length(simulation.objects)
      (f1, f2) = coulombs(simulation.objects[i], simulation.objects[j])
      forces[i] = forces[i] + f1
      forces[j] = forces[j] + f2
    end
  end
  for i in 1:length(simulation.objects)
    particle = simulation.objects[i]
    acceleration = applyForce(particle, forces[i], 1)
    applyAcceleration(particle, acceleration, 1)
    applyVelocity(particle, 1)
  end
end

"Adds a particle to a simulation"
function addParticle!(simulation::Simulation, particle::Particle)
  push!(simulation.objects, particle)
end

"Removes a particle from a simulation"
function removeParticle!(simulation::Simulation, particle::Particle)
  filter!(a -> a != particle, simulation.objects)
end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
using Rheologies ; const R = Rheologies
using JuAFEM ; const J = JuAFEM
using Test
using SafeTestsets

@testset "variables interpolations" begin
    @test R.get_variables_interpolation((:u,), (2,), Cell{2,3,3}) ==
        Lagrange{2,RefTetrahedron,2}()
    @test R.get_variables_interpolation((:u,), (1,), Cell{2,3,3}) ==
        Lagrange{2,RefTetrahedron,1}()
    @test R.get_variables_interpolation((:u,:p), (2,1), Cell{2,3,3}) ==
        (Lagrange{2,RefTetrahedron,2}(), Lagrange{2,RefTetrahedron,1}())
    @test R.get_variables_interpolation((:u,:p), (2,1), Cell{2,4,4}) ==
        (Lagrange{2,RefCube,2}(), Lagrange{2,RefCube,1}())
end

@testset "get quad rules" begin
    @test typeof(R.get_quadrature_rules(2, :legendre, Cell{2,4,4})[1]) ==
        typeof(QuadratureRule{2,RefCube}(2))
    @test typeof(R.get_quadrature_rules(2, :legendre, Cell{2,4,4})[2]) ==
        typeof(QuadratureRule{1,RefCube}(2))
end

@testset "create values" begin
    el_geom = Cell{2,3,3}
    qr, qr_face = R.get_quadrature_rules(3, :legendre, el_geom)
    interp_u1 = R.get_variables_interpolation((:u,), (2,), el_geom)
    interp_u2 = R.get_variables_interpolation((:u,:p), (2,1), el_geom)
    interp_geom = J.default_interpolation(el_geom)
    @test length(R.create_values(qr, qr_face, interp_geom, interp_u1)) == 2
    @test length(R.create_values(qr, qr_face, interp_geom, interp_u2...)) == 3
end

@testset "create material properties" begin
    func_1D(x) = x[1]
    func_2D(x) = x[1] + x[2]
    elasticity1D = IncompressibleElasticity(E = func_1D, ν =0.4999)
    elasticity2D = IncompressibleElasticity(E = func_2D, ν =0.4999)

    r1D = Rheology(nothing,nothing,elasticity1D,nothing)
    r2D = Rheology(nothing,nothing,elasticity2D,nothing)

    grid1D = generate_grid(Line,(2,), Vec(0.0), Vec(10.0))
    grid2D = generate_grid(Quadrilateral,(2, 2), Vec(0.0,0.0), Vec(10.0,10.0))

    mp1D = create_material_properties(grid1D, r1D)
    mp2D = create_material_properties(grid1D, r1D)

    @test r1D(1.0).elasticity.E == 1.0
    @test r2D([1.0,2.0]).elasticity.E == 3.0

end


text:
@testset "$(format("Chained ADMM: Bernoulli"))"  begin
    let
        Random.seed!(65536)
        ROW = 500
        COL = 500
        for input_rank in collect(80:10:120)
            for input_sample in collect(50:5:90)
                try
                    @show(Pair(input_rank, input_sample))
                    truth_matrix        = rand([(FixedRankMatrix(Distributions.Bernoulli(0.5), rank = input_rank), ROW, COL)])
                    sample_model        = provide(Sampler{BernoulliModel}(), rate = input_sample / 100)
                    input_matrix        = sample_model.draw(truth_matrix)
                    manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
                    manual_type_matrix .= :Bernoulli
                    let
                        RESULTS_DIR    = GLOBAL_SIMULATION_RESULTS_DIR *
                            "bernoulli/small(" * string(ROW) * "x" * string(COL) * ")" * "(vary_missing_200_iter_max)/" *
                            "rank" * string(input_rank) * "/"  *
                            "sample" * string(input_sample) * "/oneshot/"
                        LOG_FILE_NAME  = "io.log"
                        DATA_FILE_NAME = "saved_variables.h5"
                        LOG_FILE_PATH  = RESULTS_DIR * LOG_FILE_NAME
                        DATA_FILE_PATH = RESULTS_DIR * DATA_FILE_NAME
                        Base.Filesystem.mkpath(RESULTS_DIR)
                        io = open(LOG_FILE_PATH, "w")
                        completed_matrix, type_tracker, tracker, imputed = complete(OneShotADMM(),
                                                                                A               = input_matrix,
                                                                                maxiter         = 1000,
                                                                                ρ               = 0.3,
                                                                                use_autodiff    = false,
                                                                                gd_iter         = 3,
                                                                                io              = io,
                                                                                debug_mode      = false,
                                                                                project_rank    = nothing,
                                                                                type_assignment = manual_type_matrix)
                        predicted_matrix = predict(MatrixCompletionModel(),
                                                   completed_matrix = completed_matrix,
                                                   type_tracker     = type_tracker)

                        summary_object  = summary(MatrixCompletionModel(),
                                                  predicted_matrix = predicted_matrix,
                                                  truth_matrix     = truth_matrix,
                                                  type_tracker     = type_tracker,
                                                  tracker          = tracker)
                        pickle(DATA_FILE_PATH,
                               "missing_idx"      => type_tracker[:Missing],
                               "completed_matrix" => completed_matrix,
                               "predicted_matrix" => predicted_matrix,
                               "truth_matrix"     => truth_matrix,
                               "summary"          => summary_object)
                        print(io, JSON.json(summary_object, 4))
                        close(io)
                    end

                    let
                        RESULTS_DIR    = GLOBAL_SIMULATION_RESULTS_DIR *
                            "bernoulli/small(" * string(ROW) * "x" * string(COL) * ")" * "(vary_missing_200_iter_max)/" *
                            "rank" * string(input_rank) * "/"  *
                            "sample" * string(input_sample) * "/chained/"
                        LOG_FILE_NAME  = "io.log"
                        DATA_FILE_NAME = "saved_variables.h5"
                        LOG_FILE_PATH  = RESULTS_DIR * LOG_FILE_NAME
                        DATA_FILE_PATH = RESULTS_DIR * DATA_FILE_NAME
                        Base.Filesystem.mkpath(RESULTS_DIR)
                        io = open(LOG_FILE_PATH, "w")

                        completed_matrix, type_tracker, tracker, imputed = complete(ChainedADMM(),
                                                                                A               = deepcopy(input_matrix),
                                                                                maxiter         = 200,
                                                                                ρ               = 0.3,
                                                                                use_autodiff    = false,
                                                                                gd_iter         = 3,
                                                                                imputation_round = 5,
                                                                                io              = io,
                                                                                debug_mode      = false,
                                                                                project_rank    = nothing,
                                                                                type_assignment = deepcopy(manual_type_matrix))
                        predicted_matrix = predict(MatrixCompletionModel(),
                                                   completed_matrix = completed_matrix,
                                                   type_tracker     = type_tracker)

                        summary_object  = summary(MatrixCompletionModel(),
                                                  predicted_matrix = predicted_matrix,
                                                  truth_matrix     = truth_matrix,
                                                  type_tracker     = type_tracker,
                                                  tracker          = tracker)
                        pickle(DATA_FILE_PATH,
                               "missing_idx"      => type_tracker[:Missing],
                               "completed_matrix" => completed_matrix,
                               "predicted_matrix" => predicted_matrix,
                               "truth_matrix"     => truth_matrix,
                               "summary"          => summary_object)
                        print(io, JSON.json(summary_object, 4))
                        close(io)
                    end
                catch
                    # nothing
                    @info("got exception not log")
                end
            end
        end
    end
end


text:
validate_field_grid(grid, field) = grid === field.grid

validate_field_grid(grid, field_tuple::NamedTuple) =
    all(validate_field_grid(grid, field) for field in field_tuple)

"""
    validate_field_tuple_grid(tuple_name, field_tuple, arch, grid, bcs)

Validates the grids associated with grids in the (possibly nested) `field_tuple`,
and returns `field_tuple` if validation succeeds.
"""
function validate_field_tuple_grid(tuple_name, field_tuple, grid)

    all(validate_field_grid(grid, field) for field in field_tuple) ||
        throw(ArgumentError("Model grid and $tuple_name grid are not identical! " *
                            "Check that the grid used to construct $tuple_name has the correct halo size."))

    return nothing
end


text:
#=
int arf_addmul(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)
Performs a fused multiply-add z=z+x⋅y, updating z in-place.

int arf_submul(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)
Performs a fused multiply-subtract z=z−x⋅y, updating z in-place.
=#

function muladd(x::ArbFloat{P}, y::ArbFloat{P}, z::ArbFloat{P}, roundingmode::RoundingMode) where {P}
    rounding = match_rounding_mode(roundingmode)
    zz = copy(z)
    res = ccall(@libarb(arf_addmul), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint),
                 zz, x, y, P, rounding)
    return zz
end

@inline muladd(x::ArbFloat{P}, y::ArbFloat{P}, z::ArbFloat{P}) where {P} =
    muladd(x, y, z, RoundNearest)

@inline fma(x::ArbFloat{P}, y::ArbFloat{P}, z::ArbFloat{P}) where {P} =
    muladd(x, y, z, RoundNearest)

#=
void arb_addmul(arb_t z, const arb_t x, const arb_t y, slong prec)
void arb_addmul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)
void arb_addmul_si(arb_t z, const arb_t x, slong y, slong prec)
void arb_addmul_ui(arb_t z, const arb_t x, ulong y, slong prec)
void arb_addmul_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec)¶
Sets z=z+x⋅y, rounded to prec bits.

void arb_submul(arb_t z, const arb_t x, const arb_t y, slong prec)¶
void arb_submul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)
void arb_submul_si(arb_t z, const arb_t x, slong y, slong prec)
void arb_submul_ui(arb_t z, const arb_t x, ulong y, slong prec)
void arb_submul_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec)
Sets z=z−x⋅y, rounded to prec bits. The precision can be ARF_PREC_EXACT provided that the result fits in memory.
=#

function muladd(x::ArbReal{P}, y::ArbReal{P}, z::ArbReal{P}) where {P}
    zz = copy(z)
    res = ccall(@libarb(arb_addmul), Cint, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong),
                 zz, x, y, P)
    return zz
end

function fma(x::ArbReal{P}, y::ArbReal{P}, z::ArbReal{P}) where {P}
    xx = ArbFloat{P}(x)
    yy = ArbFloat{P}(y)
    zz = ArbFloat{P}(z)
    zz = fma(xx, yy, zz)
    zzz = copy(z)
    res = ccall(@libarb(arb_addmul), Cint, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong),
                 zzz, x, y, P)
    return setball(zzz, radius(z))
end

#=
void acb_addmul(acb_t z, const acb_t x, const acb_t y, slong prec)
void acb_addmul_ui(acb_t z, const acb_t x, ulong y, slong prec)
void acb_addmul_si(acb_t z, const acb_t x, slong y, slong prec)
void acb_addmul_fmpz(acb_t z, const acb_t x, const fmpz_t y, slong prec)
void acb_addmul_arb(acb_t z, const acb_t x, const arb_t y, slong prec)
Sets z to z plus the product of x and y.

void acb_submul(acb_t z, const acb_t x, const acb_t y, slong prec)
void acb_submul_ui(acb_t z, const acb_t x, ulong y, slong prec)
void acb_submul_si(acb_t z, const acb_t x, slong y, slong prec)
void acb_submul_fmpz(acb_t z, const acb_t x, const fmpz_t y, slong prec)
void acb_submul_arb(acb_t z, const acb_t x, const arb_t y, slong prec)
Sets z to z minus the product of x and y.
=#

function muladd(x::ArbComplex{P}, y::ArbComplex{P}, z::ArbComplex{P}) where {P}
    zz = copy(z)
    res = ccall(@libarb(acb_addmul), Cint, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong),
                 zz, x, y, P)
    return zz
end


text:
using MarketData # to provide mark_price

txn = Input(res[1].value)

cash = foldl(100_000.00, txn) do init, t
    init - t.quantity * t.basis
end

A = foldl(Stock(:AAPL, 0, 0.0), txn) do asset, t
    transact(asset, t.quantity, t.basis)
end

mark_price = Input(cl.values[1])

Pval = foldl(0.0, A, mark_price) do init, asset, closing_price
    asset.quantity * closing_price
end

Aval = foldl(0.0, cash, Pval) do init, c, p
    c + p
end


text:
struct QuasiArray{T,N,AXES<:NTuple{N,AbstractVector}} <: AbstractQuasiArray{T,N}
    parent::Array{T,N}
    axes::AXES

    function QuasiArray{T,N,AXES}(par::AbstractArray{T,N}, axes::AXES) where {T,N,AXES<:NTuple{N,AbstractVector}}
        size(par) == length.(axes) || throw(ArgumentError("Axes must be compatible with parent dimensions"))
        new{T,N,AXES}(convert(Array{T,N},par),axes)
    end
end

QuasiArray{T,N,AXES}(par::AbstractArray{<:Any,N}, axes::AXES) where {T,N,AXES<:NTuple{N,AbstractVector}} =
    QuasiArray{T,N,AXES}(convert(AbstractArray{T,N}, par)::AbstractArray{T,N}, axes)

const QuasiMatrix{T,AXES<:NTuple{2,AbstractVector}} = QuasiArray{T,2,AXES}
const QuasiVector{T,AXES<:Tuple{AbstractVector}} = QuasiArray{T,1,AXES}

QuasiArray{T,N}(::UndefInitializer, axes::NTuple{N,AbstractVector}) where {T,N} =
    QuasiArray(Array{T}(undef, map(length,axes)), axes)
QuasiArray{T,N}(::UndefInitializer, axes::Vararg{AbstractVector,N}) where {T,N} =
    QuasiArray{T,N}(undef, axes)
QuasiVector(::UndefInitializer, axes::AbstractVector) where T =
    QuasiArray(Vector(undef,length(axes)), (axes,))
QuasiMatrix(::UndefInitializer, ax1::AbstractVector, ax2::AbstractVector) where T =
    QuasiArray(Matrix(undef,length(ax1),length(ax2)), (ax1,ax2))

QuasiArray(par::AbstractArray{T,N}, axes::NTuple{N,AbstractVector}) where {T,N} =
    QuasiArray{T,N,typeof(axes)}(par, axes)
QuasiMatrix(par::AbstractMatrix{T}, axes::NTuple{2,AbstractVector}) where T =
    QuasiArray{T,2,typeof(axes)}(par, axes)
QuasiVector(par::AbstractVector{T}, axes::Tuple{AbstractVector}) where T =
    QuasiArray{T,1,typeof(axes)}(par, axes)

QuasiArray{T}(par::AbstractArray{<:Any,N}, axes::NTuple{N,AbstractVector}) where {T,N} =
    QuasiArray{T,N,typeof(axes)}(par, axes)
QuasiArray{T}(par::AbstractArray{T,N}, axes::NTuple{N,AbstractVector}) where {T,N} =
    QuasiArray{T,N,typeof(axes)}(par, axes)
QuasiMatrix{T}(par::AbstractMatrix, axes::NTuple{2,AbstractVector}) where T =
    QuasiMatrix{T,typeof(axes)}(par, axes)
QuasiVector{T}(par::AbstractVector, axes::Tuple{AbstractVector}) where T =
    QuasiVector{T,typeof(axes)}(par, axes)

QuasiArray(par::AbstractArray{T,N}, axes::Vararg{AbstractVector,N}) where {T,N} =
    QuasiArray(par, axes)
QuasiMatrix(par::AbstractMatrix{T}, axes::Vararg{AbstractVector,2}) where T =
    QuasiMatrix(par, axes)
QuasiVector(par::AbstractVector{T}, axes::AbstractVector) where T =
    QuasiVector(par, (axes,))

QuasiArray{T}(par::AbstractArray{<:Any,N}, axes::Vararg{AbstractVector,N}) where {T,N} =
    QuasiArray{T}(par, axes)
QuasiMatrix{T}(par::AbstractMatrix, axes::Vararg{AbstractVector,2}) where T =
    QuasiMatrix{T}(par, axes)
QuasiVector{T}(par::AbstractVector, axes::AbstractVector) where T =
    QuasiVector{T}(par, (axes,))

QuasiMatrix(λ::UniformScaling, ax::NTuple{2,AbstractVector}) = QuasiMatrix(Matrix(λ,map(length,ax)...), ax)
QuasiMatrix{T}(λ::UniformScaling, ax::NTuple{2,AbstractVector}) where T = QuasiMatrix{T}(Matrix(λ,map(length,ax)...), ax)
QuasiMatrix(λ::UniformScaling, ax1::AbstractVector, ax2::AbstractVector) = QuasiMatrix(λ, (ax1, ax2))
QuasiMatrix{T}(λ::UniformScaling, ax1::AbstractVector, ax2::AbstractVector) where T = QuasiMatrix{T}(λ, (ax1, ax2))

QuasiArray(par::AbstractArray) = QuasiArray(par, axes(par))
QuasiMatrix(par::AbstractArray) = QuasiMatrix(par, axes(par))
QuasiVector(par::AbstractArray) = QuasiVector(par, axes(par))

QuasiArray{T,N,AXES}(par::AbstractArray{T,N}, axes::NTuple{N,AbstractQuasiOrVector}) where {T,N,AXES} =
    QuasiArray{T,N,AXES}(par, map(domain,axes))
QuasiArray{T,N}(par::AbstractArray{T,N}, axes::NTuple{N,AbstractQuasiOrVector}) where {T,N} =
    QuasiArray{T,N}(par, map(domain,axes))
QuasiArray{T}(par::AbstractArray{T,N}, axes::NTuple{N,AbstractQuasiOrVector}) where {T,N} =
    QuasiArray{T}(par, map(domain,axes))
QuasiArray(par::AbstractArray{T,N}, axes::NTuple{N,AbstractQuasiOrVector}) where {T,N} =
    QuasiArray(par, map(domain,axes))
QuasiMatrix(par::AbstractMatrix{T}, axes::NTuple{2,AbstractQuasiOrVector}) where T =
    QuasiMatrix(par, map(domain,axes))
QuasiVector(par::AbstractVector{T}, axes::Tuple{AbstractQuasiOrVector}) where T =
    QuasiVector(par, map(domain,axes))

QuasiVector(par::AbstractVector{T}, axes::AbstractArray) where {T} =
    QuasiVector(par, (axes,))

QuasiArray(a::AbstractQuasiArray{T}) where T = QuasiArray{T}(a[map(collect,axes(a))...], axes(a))
QuasiArray{T}(a::AbstractQuasiArray) where T = QuasiArray(convert(AbstractArray{T},a[map(collect,axes(a))...]), axes(a))
QuasiArray{T,N}(a::AbstractQuasiArray{<:Any,N}) where {T,N} = QuasiArray(convert(AbstractArray{T,N},a[map(collect,axes(a))...]), axes(a))
QuasiArray{T,N,AXES}(a::AbstractQuasiArray{<:Any,N}) where {T,N,AXES} = QuasiArray{T,N,AXES}(Array{T}(a), axes(a))
QuasiMatrix(a::AbstractQuasiMatrix) = QuasiArray(a)
QuasiVector(a::AbstractQuasiVector) = QuasiArray(a)

convert(::Type{T}, a::AbstractQuasiArray) where {T<:QuasiArray} = a isa T ? a : T(a)
convert(::Type{QuasiArray{T,N}}, a::AbstractQuasiArray) where {T,N} = a isa QuasiArray{T,N} ? a : QuasiArray{T,N}(a)
convert(::Type{QuasiArray{T}}, a::AbstractQuasiArray) where T = a isa QuasiArray{T} ? a : QuasiArray{T}(a)
convert(::Type{AbstractQuasiArray{T,N}}, a::QuasiArray) where {T,N} = convert(QuasiArray{T,N}, a)
convert(::Type{AbstractQuasiArray{T}}, a::QuasiArray) where T = convert(QuasiArray{T}, a)
convert(::Type{AbstractQuasiArray{T,N}}, a::QuasiArray{T,N}) where {T,N} = a
convert(::Type{AbstractQuasiArray{T}}, a::QuasiArray{T}) where T = a

AbstractQuasiArray{T,N}(a::QuasiArray) where {T,N} = QuasiArray{T,N}(a)
AbstractQuasiArray{T}(a::QuasiArray) where T = QuasiArray{T}(a)

_inclusion(d::Slice) = d
_inclusion(d::OneTo) = d
_inclusion(d::IdentityUnitRange{<:Integer}) = d
_inclusion(d::AbstractUnitRange{<:Integer}) = IdentityUnitRange(d)
_inclusion(d) = Inclusion(d)

axes(A::QuasiArray) = map(_inclusion,A.axes)
parent(A::QuasiArray) = A.parent

@propagate_inbounds @inline function _getindex(::Type{IND}, A::QuasiArray, I::IND) where IND
    @boundscheck checkbounds(A, I...)
    A.parent[findfirst.(isequal.(I), A.axes)...]
end

@propagate_inbounds @inline function _setindex!(::Type{IND}, A::QuasiArray, v, I::IND) where IND
    @boundscheck checkbounds(A, I...)
    @inbounds A.parent[findfirst.(isequal.(I), A.axes)...] = v
    A
end

function _quasimatrix_pow(A, p)
    axes(A,1) == axes(A,2) || throw(DimensionMismatch("axes must match"))
    QuasiArray(A.parent^p, A.axes)
end

^(A::QuasiMatrix, p::Number) = _quasimatrix_pow(A, p)
^(A::QuasiMatrix, p::Integer) = _quasimatrix_pow(A, p)

==(A::QuasiArray{T,N,NTuple{N,OneTo{Int}}}, B::AbstractArray{V,N}) where {T,V,N} =
    axes(A) == axes(B) && A.parent == B
==(B::AbstractArray{V,N}, A::QuasiArray{T,N,NTuple{N,OneTo{Int}}}) where {T,V,N} =
    A == B

function reshape(A::QuasiVector, ax::Tuple{Any,OneTo{Int}})
    @assert ax == (axes(A,1),Base.OneTo(1))
    QuasiMatrix(reshape(A.parent,size(A.parent,1),1), (A.axes[1], Base.OneTo(1)))
end


text:
using Foo;
Foo.command_main();


text:
using Enzyme
using Test

@testset "ABI & Calling convention" begin

    f(x) = x

    # GhostType -> Nothing
    res = autodiff(f, Const(nothing))
    @test res === nothing

    # ConstType -> Type{Int}
    res = autodiff(f, Const(Int))
    @test res === nothing

    cres,  = Enzyme.autodiff(f, Active(1.5 + 0.7im))
    @test cres ≈ 1.0 + 0.0im

    unused(_, y) = y
    res0, = autodiff(unused, Const(nothing), Active(2.0))
    @test res0 ≈ 1.0

    # Multi arg => sret
    mul(x, y) = x * y
    pair = autodiff(mul, Active(2.0), Active(3.0))
    @test pair[1] ≈ 3.0
    @test pair[2] ≈ 2.0

    # SequentialType
    struct Foo
        baz::Int
        qux::Float64
    end

    g(x) = x.qux
    res2,  = autodiff(g, Active(Foo(3, 1.2)))
    @test res2.qux ≈ 1.0

    unused2(_, y) = y.qux
    resF, = autodiff(unused2, Const(nothing), Active(Foo(3, 2.0)))
    @test resF.qux ≈ 1.0

    h(x, y) = x.qux * y.qux
    res3 = autodiff(h, Active(Foo(3, 1.2)), Active(Foo(5, 3.4)))
    @test res3[1].qux ≈ 3.4
    @test res3[2].qux ≈ 1.2

    caller(f, x) = f(x)
    res4, = autodiff(caller, (x)->x, Active(3.0))
    @test res4 ≈ 1.0

    struct LList
        next::Union{LList,Nothing}
        val::Float64
    end

    function sumlist(n::LList)
        sum = 0.0
        while n !== nothing
            sum += n.val
            n = n.next
        end
        sum
    end

    regular = LList(LList(nothing, 1.0), 2.0)
    shadow  = LList(LList(nothing, 0.0), 0.0)
    ad = autodiff(sumlist, Duplicated(regular, shadow))
    @test ad === nothing
    @test shadow.val ≈ 1.0 && shadow.next.val ≈ 1.0

    mulr(x, y) = x[] * y[]
    x = Ref(2.0)
    y = Ref(3.0)
    dx = Ref(0.0)
    dy = Ref(0.0)
    n = autodiff(mulr, Duplicated(x, dx), Duplicated(y, dy))
    @test n === nothing
    @test dx[] ≈ 3.0
    @test dy[] ≈ 2.0

    mid, = Enzyme.autodiff((fs, x) -> fs[1](x), (x->x*x,), Active(2.0))
    @test mid ≈ 4.0

    mid, = Enzyme.autodiff((fs, x) -> fs[1](x), [x->x*x], Active(2.0))
    @test mid ≈ 4.0

    # deserves_argbox yes and no
    struct Bar
        r::Ref{Int}
    end

    # ConstType

    # primitive type Int128, Float64, Float128

    # returns: sret, const/ghost, !deserve_retbox
end

@testset "Callable ABI" begin
    function method(f, x)
        return f(x)
    end

    struct AFoo
       x::Float64
    end

    function (f::AFoo)(x::Float64)
       return f.x * x
    end

    @test Enzyme.autodiff(method, AFoo(2.0), Active(3.0))[1]≈ 2.0
    @test Enzyme.autodiff(AFoo(2.0), Active(3.0))[1]≈ 2.0

    struct ABar
    end

    function (f::ABar)(x::Float64)
       return 2.0 * x
    end

    @test Enzyme.autodiff(method, ABar(), Active(3.0))[1]≈ 2.0
    @test Enzyme.autodiff(ABar(), Active(3.0))[1]≈ 2.0
end


text:
module MyModule
function complexfunc(dims::Int, iters::Int)::Float64
    x::Array{Float64, 2} = ones(Float64, (dims, dims))
    for i in 1:iters
        x .+= 1.0
    end

    z::Float64 = 0.0
    for j in 1:dims
        z += x[j, j]
    end

    z
end
end

max_stars_repo_name: gobuki/jlrs

max_issues_repo_name: gobuki/jlrs

max_forks_repo_name: gobuki/jlrs


text:
module AdversarialPrediction

# import for expression.jl
import Base.sqrt, Base.log, Base.exp

# packages used in metric.jl
using LinearAlgebra
using LBFGSB

# packages used in nn.jl
using Flux.Zygote
using Flux.Zygote: @adjoint
using Requires

include("expression.jl")
include("projection.jl")
include("metric.jl")
include("nn.jl")

# common metrics
include("common_metrics/CommonMetrics.jl")

# if Gurobi or CuArrays are loaded
function __init__()
    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("nncuda.jl")
end

export ConfusionMatrix, CM_Value, PerformanceMetric
export @metric, define, constraint
export special_case_positive!, special_case_negative!, cs_special_case_positive!, cs_special_case_negative!
export compute_metric, compute_constraints, objective
export ap_objective

end # module