text:
const testterm = treat(:g, TR, PR)

function valid_didargs(d::Type{TestDID}, ::AbstractTreatment, ::AbstractParallel,
        args::Dict{Symbol,Any})
    name = get(args, :name, "")
    ns = (setdiff([keys(args)...], [:name, :d])...,)
    args = NamedTuple{ns}(map(n->args[n], ns))
    return name, d, args
end

@testset "DiffinDiffsEstimator" begin
    d = DefaultDID()
    @test length(d) == 0
    @test eltype(d) == StatsStep
    @test_throws BoundsError d[1]
    @test collect(d) == StatsStep[]
    @test iterate(d) === nothing

    d = TestDID()
    @test length(d) == 2
    @test eltype(d) == StatsStep
    @test d[1] == TestStep()
    @test d[1:2] == [TestStep(), TestNextStep()]
    @test d[[2,1]] == [TestNextStep(), TestStep()]
    @test_throws BoundsError d[3]
    @test_throws MethodError d[:a]
    @test collect(d) == StatsStep[TestStep(), TestNextStep()]
    @test iterate(d) == (TestStep(), 2)
    @test iterate(d, 2) === (TestNextStep(), 3)
    @test iterate(d, 3) === nothing
end

@testset "parse_didargs!" begin
    args = Dict{Symbol,Any}(:xterms=>(term(:x),), :treatintterms=>:z)
    _totermset!(args, :xterms)
    @test args[:xterms] == TermSet(term(:x))
    _totermset!(args, :treatintterms)
    @test args[:treatintterms] == TermSet(term(:z))

    @test parse_didargs!(Any["test"], Dict{Symbol,Any}()) == Dict{Symbol,Any}(:name=>"test")

    args = parse_didargs!([TestDID, TR, PR], Dict{Symbol,Any}(:a=>1, :b=>2))
    @test args == Dict{Symbol,Any}(pairs((d=TestDID, tr=TR, pr=PR, a=1, b=2)))

    args = parse_didargs!(["test", testterm, TestDID], Dict{Symbol,Any}())
    @test args == Dict{Symbol,Any}(pairs((d=TestDID, tr=TR, pr=PR, name="test", treatname=:g)))

    args0 = parse_didargs!([TestDID, term(:y) ~ testterm, "test"], Dict{Symbol,Any}())
    @test args0 == Dict{Symbol,Any}(pairs((d=TestDID, tr=TR, pr=PR, name="test", yterm=term(:y), treatname=:g, treatintterms=TermSet(), xterms=TermSet())))

    args1 = parse_didargs!([TestDID, "test", @formula(y ~ treat(g, ttreat(t, 0), tpara(0)))],
        Dict{Symbol,Any}())
    @test args1 == args0

    args0 = parse_didargs!([TestDID, term(:y) ~ testterm & term(:z) + term(:x)],
        Dict{Symbol,Any}())
    @test args0 == Dict{Symbol,Any}(pairs((d=TestDID, tr=TR, pr=PR, yterm=term(:y),
        treatname=:g, treatintterms=TermSet(term(:z)), xterms=TermSet(term(:x)))))

    args1 = parse_didargs!([TestDID,
        @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & z + x)], Dict{Symbol,Any}())
    @test args1 == args0

    @test_throws ArgumentError parse_didargs!(['a', :a, 1], Dict{Symbol,Any}())
    @test parse_didargs!(Any[TestDID, DefaultDID], Dict{Symbol,Any}()) == Dict{Symbol,Any}(:d=>DefaultDID)
end

@testset "valid_didargs" begin
    arg = Dict{Symbol,Any}(pairs((d=TestDID, tr=TR, pr=PR)))
    @test valid_didargs(arg) == ("", TestDID, (tr=TR, pr=PR))
    arg = Dict{Symbol,Any}(pairs((name="name", d=TestDID, tr=TR, pr=PR)))
    @test valid_didargs(arg) == ("name", TestDID, (tr=TR, pr=PR))
    arg = Dict{Symbol,Any}(pairs((d=TestDID, tr=TR, a=1)))
    @test_throws ArgumentError valid_didargs(arg)
    arg = Dict{Symbol,Any}(pairs((tr=TR, pr=PR)))
    @test_throws ErrorException valid_didargs(arg)
    arg = Dict{Symbol,Any}(pairs((d=NotImplemented, tr=TR, pr=PR)))
    @test_throws ErrorException valid_didargs(arg)
end

@testset "StatsSpec" begin
    @testset "== ≊" begin
        sp1 = StatsSpec("", DefaultDID, NamedTuple())
        sp2 = StatsSpec("", DefaultDID, NamedTuple())
        @test sp1 == sp2
        @test sp1 ≊ sp2

        sp1 = StatsSpec("", DefaultDID, (tr=TR, pr=PR, a=1, b=2))
        sp2 = StatsSpec("", DefaultDID, (pr=PR, tr=TR, b=2, a=1))
        @test sp1 != sp2
        @test sp1 ≊ sp2

        sp2 = StatsSpec("", DefaultDID, (tr=TR, pr=PR, a=1.0, b=2.0))
        @test sp1 == sp2
        @test sp1 ≊ sp2

        sp2 = StatsSpec("", DefaultDID, (tr=TR, pr=PR, a=1, b=2, c=3))
        @test !(sp1 ≊ sp2)

        sp2 = StatsSpec("", DefaultDID, (tr=TR, pr=PR, a=1, b=1))
        @test !(sp1 ≊ sp2)
    end

    @testset "show" begin
        sp = StatsSpec("", DefaultDID, NamedTuple())
        @test sprint(show, sp) == "unnamed"
        @test sprint(show, MIME("text/plain"), sp) == "unnamed (StatsSpec for DefaultDID)"

        sp = StatsSpec("name", DefaultDID, NamedTuple())
        @test sprint(show, sp) == "name"
        @test sprint(show, MIME("text/plain"), sp) == "name (StatsSpec for DefaultDID)"

        sp = StatsSpec("", TestDID, (tr=dynamic(:time,-1), pr=nevertreated(-1)))
        @test sprint(show, sp) == "unnamed"
        @test sprint(show, MIME("text/plain"), sp) == """
            unnamed (StatsSpec for TestDID):
              Dynamic{S}(-1)
              NeverTreated{U,P}(-1)"""

        sp = StatsSpec("", TestDID, (tr=dynamic(:time,-1),))
        @test sprint(show, sp) == "unnamed"
        @test sprint(show, MIME("text/plain"), sp) == """
            unnamed (StatsSpec for TestDID):
              Dynamic{S}(-1)"""

        sp = StatsSpec("name", TestDID, (pr=notyettreated(-1),))
        @test sprint(show, sp) == "name"
        @test sprint(show, MIME("text/plain"), sp) == """
            name (StatsSpec for TestDID):
              NotYetTreated{U,P}(-1)"""
    end
end

@testset "didspec @did" begin
    @test_throws ArgumentError @did [noproceed]
    @test_throws ArgumentError didspec()

    sp0 = StatsSpec("", TestDID, (tr=TR, pr=PR, a=1, b=2))
    sp1 = didspec(TestDID, TR, PR, a=1, b=2)
    @test sp1 ≊ sp0
    @test sp0 ≊ @did [noproceed] TR a=1 b=2 PR TestDID

    sp2 = StatsSpec("name", TestDID, (tr=TR, pr=PR, a=1, b=2))
    sp3 = didspec("name", TR, PR, TestDID, b=2, a=1)
    @test sp2 ≊ sp1
    @test sp3 ≊ sp2
    @test sp3 ≊ @did [noproceed] TR PR TestDID "name" b=2 a=1

    sp4 = StatsSpec("name", TestDID, (tr=TR, pr=PR, treatname=:g))
    sp5 = didspec(TestDID, testterm)
    @test sp5 ≊ sp4
    @test sp4 ≊ @did [noproceed] TestDID testterm "name"

    sp6 = StatsSpec("", TestDID, (tr=TR, pr=PR, yterm=term(:y), treatname=:g,
        treatintterms=TermSet(term(:z)), xterms=TermSet(term(:x))))
    sp7 = didspec(TestDID, term(:y) ~ testterm & term(:z) + term(:x))
    sp8 = didspec(TestDID, @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & z + x))
    @test sp7 ≊ sp6
    @test sp8 ≊ sp7
    @test sp6 ≊ @did [noproceed] TestDID term(:y) ~ testterm & term(:z) + term(:x)
    @test sp6 ≊ @did [noproceed] TestDID @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & z + x)
end

@testset "did @did" begin
    r = result(TestDID, NamedTuple()).result
    d0 = @did TestDID TR PR
    d1 = @did PR TR TestDID
    @test d0 == r
    @test d1 == d0
    @test did(TestDID, TR, PR) == r
    @test did(PR, TR, TestDID) == r

    d = @did [keepall] TestDID testterm
    @test d ≊ (tr=TR, pr=PR, treatname=:g, str=sprint(show, TR), spr=sprint(show, PR),
        next="next"*sprint(show, TR), result=r)
    @test did(TestDID, testterm; keepall=true) == d

    d0 = @did [keepall] TestDID term(:y) ~ testterm
    @test d0 ≊ merge(d, (yterm=term(:y), treatintterms=TermSet(), xterms=TermSet()))
    @test did(TestDID, term(:y) ~ testterm; keepall=true) == d0
    d1 = @did [keepall] TestDID @formula(y ~ treat(g, ttreat(t, 0), tpara(0))) "test"
    @test d1 ≊ d0
    @test did(TestDID, @formula(y ~ treat(g, ttreat(t, 0), tpara(0))); keepall=true) == d1

    d0 = @did [keepall] TestDID term(:y) ~ testterm & term(:z) + term(:x)
    @test d0 ≊ merge(d, (yterm=term(:y), treatintterms=TermSet(term(:z)),
        xterms=TermSet(term(:x))))
    @test did(TestDID, term(:y) ~ testterm & term(:z) + term(:x); keepall=true) == d0
    d1 = @did [keepall] "test" TestDID @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & z + x)
    @test d1 ≊ d0
    @test did(TestDID, @formula(y ~ treat(g, ttreat(t, 0), tpara(0)) & z + x); keepall=true) == d1

    d = @did [keep=:treatname] TestDID testterm
    @test d ≊ (treatname=:g, result=r)
    @test did(TestDID, testterm; keep=:treatname) == d
    d = @did [keep=[:treatname,:tr]] TestDID testterm
    @test d ≊ (treatname=:g, tr=TR, result=r)
    @test did(TestDID, testterm; keep=[:treatname,:tr]) == d

    @test_throws ArgumentError @did TestDID
    @test_throws ArgumentError did(TestDID)
    @test_throws ArgumentError @did TestDID TR
    @test_throws ArgumentError did(TestDID, TR)
    @test_throws ArgumentError @did TestDID TR PR 1
    @test_throws ArgumentError did(TestDID, TR, PR, 1)
    @test_throws ErrorException @did NotImplemented TR PR
    @test_throws ErrorException did(NotImplemented, TR, PR)

    @test_throws ArgumentError @did [keep] TestDID testterm
    @test_throws ArgumentError did(TestDID, testterm, keep=true)
    @test_throws ArgumentError @did [keep=1] TestDID testterm
    @test_throws ArgumentError did(TestDID, testterm, keep=1)
    @test_throws ArgumentError @did [keep="treatname"] TestDID testterm
    @test_throws ArgumentError did(TestDID, testterm, keep="treatname")
end

@testset "@specset @did" begin
    s = @specset [noproceed] begin
        @did TestDID TR PR
    end
    @test s == [didspec(TestDID, TR, PR)]
    @test proceed(s) == [@did TestDID TR PR]

    s = @specset [noproceed] begin
        @did TestDID TR PR
        @did PR TR TestDID
    end
    @test s == StatsSpec[didspec(TestDID, TR, PR), didspec(PR, TR, TestDID)]

    s = @specset [noproceed] DefaultDID PR a=1 begin
        @did TR PR TestDID a=2
        @did TR TestDID
    end
    @test s[1] ≊ didspec(TestDID, TR, PR; a=2)
    @test s[2] ≊ didspec(TestDID, TR, PR; a=1)

    s = @specset [noproceed] DefaultDID PR begin
        @did [keep=:treatname verbose] TestDID testterm
        @did [noproceed] TestDID @formula(y ~ treat(g, ttreat(t, 0), tpara(0)))
    end
    @test s[1] ≊ didspec(TestDID, TR, PR; treatname=:g)
    @test s[2] ≊ didspec(TestDID, TR, PR; yterm=term(:y), treatname=:g,
        treatintterms=TermSet(), xterms=TermSet())

    s = @specset [noproceed] for i in 1:2
        @did "name"*"$i" TestDID TR PR a=i
    end
    @test s[1] ≊ didspec("name1", TestDID, TR, PR; a=1)
    @test s[2] ≊ didspec("name2", TestDID, TR, PR; a=2)

    s = @specset [noproceed] DefaultDID PR a=1 begin
        @did "name" TR PR TestDID a=0
        for i in 1:2
            @did "name"*"$i" TR TestDID a=i
        end
    end
    @test s[1] ≊ didspec("name", TestDID, TR, PR; a=0)
    @test s[2] ≊ didspec("name1", TestDID, TR, PR; a=1)
    @test s[3] ≊ didspec("name2", TestDID, TR, PR; a=2)
end

@testset "DIDResult" begin
    r = TestResult(2, 2)

    @test coef(r) == r.coef
    @test coef(r, 1) == 1
    @test coef(r, "rel: 1 & c: 1") == 1
    @test coef(r, :c5) == 5
    @test coef(r, 1:2) == [1.0, 2.0]
    @test coef(r, 5:-1:3) == [5.0, 4.0, 3.0]
    @test coef(r, (1, :c5, "rel: 1 & c: 1")) == [1.0, 5.0, 1.0]
    @test coef(r, [1 "rel: 1 & c: 1" :c5]) == [1.0 1.0 5.0]
    @test coef(r, :rel=>x->x==1) == [1.0, 2.0]
    @test coef(r, :rel=>x->x==1, :c=>x->x==1) == [1.0]
    @test coef(r, (:rel, :c)=>(x, y)->x==1) == [1.0, 2.0]

    @test vcov(r) == r.vcov
    @test vcov(r, 1) == 1
    @test vcov(r, 1, 2) == 2
    @test vcov(r, "rel: 1 & c: 1") == 1
    @test vcov(r, :c5) == 25
    @test vcov(r, "rel: 1 & c: 1", :c5) == 5
    @test vcov(r, 1:2) == r.vcov[1:2, 1:2]
    @test vcov(r, 5:-1:3) == r.vcov[5:-1:3, 5:-1:3]
    @test vcov(r, (1, :c5, "rel: 1 & c: 1")) == r.vcov[[1,5,1], [1,5,1]]
    @test vcov(r, [1 :c5 "rel: 1 & c: 1"]) == r.vcov[[1,5,1], [1,5,1]]
    @test vcov(r, :rel=>x->x==1) == r.vcov[[1,2], [1,2]]
    @test vcov(r, :rel=>x->x==1, :c=>x->x==1) == reshape([1.0], 1, 1)
    @test vcov(r, (:rel, :c)=>(x, y)->x==1) == r.vcov[[1,2], [1,2]]

    @test vce(r) === nothing
    cfint = confint(r)
    @test length(cfint) == 2
    @test all(cfint[1] + cfint[2] .≈ 2 * coef(r))

    @test treatment(r) === TR
    @test nobs(r) == 6
    @test outcomename(r) == "y"
    @test coefnames(r) == r.coefnames
    @test treatcells(r) == r.treatcells
    @test weights(r) == :w
    @test ntreatcoef(r) == 4
    @test treatcoef(r) == r.coef[1:4]
    @test treatvcov(r) == r.vcov[1:4, 1:4]
    @test treatnames(r) == r.coefnames[1:4]
    @test parent(r) === r
    @test dof_residual(r) == 5
    @test responsename(r) == "y"
    @test coefinds(r) == r.coefinds
    @test ncovariate(r) == 2

    @test_throws ErrorException agg(r)

    w = VERSION < v"1.6.0" ? " " : ""
    @test sprint(show, r) == """
        ─────────────────────────────────────────────────────────────────────────
                       Estimate  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
        ─────────────────────────────────────────────────────────────────────────
        rel: 1 & c: 1       1.0         1.0  1.00    0.3632   -1.57058    3.57058
        rel: 1 & c: 2       2.0         2.0  1.00    0.3632   -3.14116    7.14116
        rel: 2 & c: 1       3.0         3.0  1.00    0.3632   -4.71175   10.7117$w
        rel: 2 & c: 2       4.0         4.0  1.00    0.3632   -6.28233   14.2823$w
        c5                  5.0         5.0  1.00    0.3632   -7.85291   17.8529$w
        c6                  6.0         6.0  1.00    0.3632   -9.42349   21.4235$w
        ─────────────────────────────────────────────────────────────────────────"""

    r = TestResultBARE(2, 2)
    @test dof_residual(r) === nothing
    @test coefinds(r) === nothing

    w = VERSION < v"1.6.0" ? " " : ""
    @test sprint(show, r) == """
        ─────────────────────────────────────────────────────────────────────────
                       Estimate  Std. Error     z  Pr(>|z|)  Lower 95%  Upper 95%
        ─────────────────────────────────────────────────────────────────────────
        rel: 1 & c: 1       1.0         1.0  1.00    0.3173  -0.959964    2.95996
        rel: 1 & c: 2       2.0         2.0  1.00    0.3173  -1.91993     5.91993
        rel: 2 & c: 1       3.0         3.0  1.00    0.3173  -2.87989     8.87989
        rel: 2 & c: 2       4.0         4.0  1.00    0.3173  -3.83986    11.8399$w
        ─────────────────────────────────────────────────────────────────────────"""
end

@testset "_treatnames" begin
    t = VecColumnTable((rel=[1, 2],))
    @test _treatnames(t) == ["rel: 1", "rel: 2"]
    r = TestResult(2, 2)
    @test _treatnames(r.treatcells) == ["rel: $a & c: $b" for a in 1:2 for b in 1:2]
end

@testset "_parse_bycells!" begin
    cells = VecColumnTable((rel=repeat(1:2, inner=2), c=repeat(1:2, outer=2)))
    bycols = copy(getfield(cells, :columns))
    _parse_bycells!(bycols, cells, :rel=>isodd)
    @test bycols[1] == isodd.(cells.rel)
    @test bycols[2] === cells.c

    bycols = copy(getfield(cells, :columns))
    _parse_bycells!(bycols, cells, :rel=>:c=>isodd)
    @test bycols[1] == isodd.(cells.c)
    @test bycols[2] === cells.c

    bycols = copy(getfield(cells, :columns))
    _parse_bycells!(bycols, cells, (:rel=>1=>isodd, 2=>1=>isodd))
    @test bycols[1] == isodd.(cells.rel)
    @test bycols[2] == isodd.(cells.rel)

    bycols = copy(getfield(cells, :columns))
    _parse_bycells!(bycols, cells, nothing)
    @test bycols[1] == cells.rel
    @test bycols[2] == cells.c

    @test_throws ArgumentError _parse_bycells!(bycols, cells, (:rel,))
end

@testset "_parse_subset _nselected" begin
    r = TestResult(2, 2)
    @test _parse_subset(r, 1:4, true) == 1:4
    @test _parse_subset(r, collect(1:4), true) == 1:4
    @test _parse_subset(r, trues(4), true) == trues(4)
    @test _parse_subset(r, :rel=>isodd, true) == ((1:6) .< 3)
    @test _parse_subset(r, :rel=>isodd, false) == ((1:4) .< 3)
    @test _parse_subset(r, (:rel=>isodd, :c=>isodd), true) == ((1:6) .< 2)
    @test _parse_subset(r, (:rel=>isodd, :c=>isodd), false) == ((1:4) .<2)
    @test _parse_subset(r, :, false) == 1:4
    @test _parse_subset(r, :, true) == 1:6

    @test _nselected(1:10) == 10
    @test _nselected(collect(1:10)) == 10
    @test _nselected(isodd.(1:10)) == 5
    @test _nselected([true, false]) == 1
    @test_throws ArgumentError _nselected(:)
end

@testset "treatindex checktreatindex" begin
    @test treatindex(10, :) == 1:10
    @test treatindex(10, 1) == 1
    @test length(treatindex(10, 11)) == 0
    @test treatindex(10, 1:5) == 1:5
    @test treatindex(10, collect(1:20)) == 1:10
    @test treatindex(10, isodd.(1:20)) == isodd.(1:10)
    @test_throws ArgumentError treatindex(10, :a)

    @test checktreatindex([3,2,1,5,4], [1,2,3])
    @test checktreatindex(1:10, 1:5)
    @test_throws ArgumentError checktreatindex([3,2,1,5,4], [3,5])
    @test_throws ArgumentError checktreatindex(10:-1:1, 1:5)
    @test checktreatindex(trues(10), trues(5))
    @test checktreatindex(:, 1:5)
    @test checktreatindex(5, 1:0)
end

@testset "SubDIDResult" begin
    r = TestResult(2, 2)
    sr = view(r, :)
    @test coef(sr) == coef(r)
    @test vcov(sr) == vcov(r)
    @test vce(sr) === nothing
    @test treatment(sr) === TR
    @test nobs(sr) == nobs(r)
    @test outcomename(sr) == outcomename(r)
    @test coefnames(sr) == coefnames(r)
    @test treatcells(sr) == treatcells(r)
    @test weights(sr) == weights(r)
    @test ntreatcoef(sr) == ntreatcoef(r)
    @test treatcoef(sr) == treatcoef(r)
    @test treatvcov(sr) == treatvcov(r)
    @test treatnames(sr) == treatnames(r)
    @test parent(sr) === r
    @test dof_residual(sr) == dof_residual(r)
    @test responsename(sr) == responsename(r)
    @test coefinds(sr) == coefinds(r)
    @test ncovariate(sr) == 2

    sr = view(r, isodd.(1:6))
    @test coef(sr) == coef(r)[[1,3,5]]
    @test vcov(sr) == vcov(r)[[1,3,5],[1,3,5]]
    @test vce(sr) === nothing
    @test treatment(sr) === TR
    @test nobs(sr) == nobs(r)
    @test outcomename(sr) == outcomename(r)
    @test coefnames(sr) == coefnames(r)[[1,3,5]]
    @test treatcells(sr) == view(treatcells(r), [1,3])
    @test weights(sr) == weights(r)
    @test ntreatcoef(sr) == 2
    @test treatcoef(sr) == treatcoef(r)[[1,3]]
    @test treatvcov(sr) == treatvcov(r)[[1,3],[1,3]]
    @test treatnames(sr) == treatnames(r)[[1,3]]
    @test parent(sr) === r
    @test dof_residual(sr) == dof_residual(r)
    @test responsename(sr) == responsename(r)
    @test coefinds(sr) == Dict("rel: 1 & c: 1"=>1, "rel: 2 & c: 1"=>2, "c5"=>3)
    @test ncovariate(sr) == 1

    @test_throws BoundsError view(r, 1:7)
    @test_throws ArgumentError view(r, [6,1])
end

@testset "TransformedDIDResult" begin
    r = TestResult(2, 2)
    m = reshape(1:36, 6, 6)
    tr = TransformedDIDResult(r, m, copy(r.coef), copy(r.vcov))

    @test coef(tr) === tr.coef
    @test vcov(tr) === tr.vcov
    @test vce(tr) === nothing
    @test treatment(tr) === TR
    @test nobs(tr) == nobs(r)
    @test outcomename(tr) == outcomename(r)
    @test coefnames(tr) === coefnames(r)
    @test treatcells(tr) === treatcells(r)
    @test weights(tr) == weights(r)
    @test ntreatcoef(tr) == 4
    @test treatcoef(tr) == tr.coef[1:4]
    @test treatvcov(tr) == tr.vcov[1:4,1:4]
    @test treatnames(tr) == treatnames(r)
    @test parent(tr) === r
    @test dof_residual(tr) == dof_residual(r)
    @test responsename(tr) == responsename(r)
    @test coefinds(tr) === coefinds(r)
    @test ncovariate(tr) == 2
end

@testset "TransSubDIDResult" begin
    r = TestResult(2, 2)
    m = reshape(1:18, 3, 6)
    inds = [2,1,6]
    tr = TransSubDIDResult(r, m, r.coef[inds], r.vcov[inds, inds], inds)

    @test coef(tr) === tr.coef
    @test vcov(tr) === tr.vcov
    @test vce(tr) === nothing
    @test treatment(tr) === TR
    @test nobs(tr) == nobs(r)
    @test outcomename(tr) == outcomename(r)
    @test coefnames(tr) == coefnames(r)[inds]
    @test treatcells(tr) == view(treatcells(r), 2:-1:1)
    @test weights(tr) == weights(r)
    @test ntreatcoef(tr) == 2
    @test treatcoef(tr) == tr.coef[1:2]
    @test treatvcov(tr) == tr.vcov[1:2,1:2]
    @test treatnames(tr) == treatnames(r)[2:-1:1]
    @test parent(tr) === r
    @test dof_residual(tr) == dof_residual(r)
    @test responsename(tr) == responsename(r)
    @test coefinds(tr) === tr.coefinds
    @test ncovariate(tr) == 1

    @test_throws BoundsError TransSubDIDResult(r, m, r.coef[1:2], r.vcov[1:2,1:2], 1:7)
    @test_throws ArgumentError TransSubDIDResult(r, m, r.coef[1:2], r.vcov[1:2,1:2], [6,1])
end

@testset "lincom" begin
    r = TestResult(2, 2)
    m = reshape(1:36, 6, 6)
    tr = lincom(r, m)
    @test coef(tr) == m * r.coef
    @test vcov(tr) == m * r.vcov * m'

    tr = lincom(r, m, nothing)
    @test coef(tr) == m * r.coef
    @test vcov(tr) == m * r.vcov * m'
    @test tr isa TransformedDIDResult

    m = reshape(1:18, 3, 6)
    inds = [1,2,6]
    tr = lincom(r, m, inds)
    @test coef(tr) == m * r.coef
    @test vcov(tr) == m * r.vcov * m'

    m1 = reshape(1:15, 3, 5)
    @test_throws DimensionMismatch lincom(r, m1)
    @test_throws DimensionMismatch lincom(r, m1, 1:3)

    @test_throws ArgumentError lincom(r, m)
    @test_throws ArgumentError lincom(r, m, 1:6)
end

@testset "rescale" begin
    r = TestResult(2, 2)
    scale = fill(2, 6)
    m = Diagonal(scale)
    tr = rescale(r, scale)
    @test coef(tr) == m * r.coef
    @test vcov(tr) == m * r.vcov * m'

    tr = rescale(r, scale, nothing)
    @test coef(tr) == m * r.coef
    @test vcov(tr) == m * r.vcov * m'
    @test tr isa TransformedDIDResult

    scale = fill(2, 3)
    m = Diagonal(scale)
    inds = [1,2,6]
    tr = rescale(r, scale, inds)
    @test coef(tr) == m * r.coef[inds]
    @test vcov(tr) == m * r.vcov[inds,inds] * m'

    tr = rescale(r, :rel=>identity)
    @test coef(tr) == r.treatcells.rel.*r.coef[1:4]
    m = Diagonal(r.treatcells.rel)
    @test vcov(tr) == m * r.vcov[1:4,1:4] * m'

    tr = rescale(r, :rel=>identity, 1:3)
    @test coef(tr) == r.treatcells.rel[1:3].*r.coef[1:3]
    m = Diagonal(r.treatcells.rel[1:3])
    @test vcov(tr) == m * r.vcov[1:3,1:3] * m'
end

@testset "post!" begin
    @test getexportformat() == DefaultExportFormat[1]
    setexportformat!(StataPostHDF())
    @test DefaultExportFormat[1] == StataPostHDF()

    f = Dict{String,Any}()
    r = TestResult(2, 2)
    post!(f, r)
    @test f["model"] == "TestResult{TestTreatment}"
    @test f["b"] == coef(r)
    @test f["V"] == vcov(r)
    @test f["vce"] == "nothing"
    @test f["N"] == nobs(r)
    @test f["depvar"] == "y"
    @test f["coefnames"][1] == "rel: 1 & c: 1"
    @test f["weights"] == "w"
    @test f["ntreatcoef"] == 4

    f = Dict{String,Any}()
    fds = Union{Symbol, Pair{String,Symbol}}[Symbol("extra$i") for i in 1:8]
    fds[1] = "e1"=>:extra1
    post!(f, StataPostHDF(), r, model="test", fields=fds, at=1:6)
    @test f["e1"] == 1
    @test f["extra2"] == "a"
    @test f["extra3"] == "a"
    @test f["extra4"] == [1.0]
    @test f["extra5"] == ["a"]
    @test f["extra6"] == [1.0 2.0]
    @test f["extra7"] == ["a"]
    @test f["extra8"] == ""

    f = Dict{String,Any}()
    @test_throws ArgumentError post!(f, StataPostHDF(), r, fields=[:extra9])
    f = Dict{String,Any}()
    @test_throws ArgumentError post!(f, StataPostHDF(), r, at=false)
    f = Dict{String,Any}()
    @test_throws ArgumentError post!(f, StataPostHDF(), r, at=1:2)

    r1 = TestResultBARE(2, 2)
    f = Dict{String,Any}()
    post!(f, StataPostHDF(), r1, at=true)
    @test f["at"] == treatcells(r1).rel
    f = Dict{String,Any}()
    post!(f, StataPostHDF(), r1, at=false)
    @test !haskey(f, "at")
end


text:
using LinearAlgebra

function readmat(n::Int)
    m = zeros(Int, n, n)
    for i = 1:n
        m[i, :] = [parse(Int, x) for x in split(readline())]
    end
    m
end

function solution()
    n_cases = parse(Int, readline())
    for x=1:n_cases
        n = parse(Int, readline())
        m = readmat(n)
        k = tr(m)
        r = n - count(allunique(m[i, :]) for i in 1:n)
        c = n - count(allunique(m[:, i]) for i in 1:n)
        println("Case #$x: $k $r $c")
    end
end

solution()

ext: jl


text:
export DefaultExtremeRef, DefaultProgressiveRef, AddExtremeConstraint, AddExtremeCollection,
       AddProgressiveCollection, AddProgressiveConstraint

"""
    @enum IterationOutcome

Has values Dominating, Improving, or Unsuccessful. Corresponding to the three
iteration outcomes in MADS-PB.
"""
@enum IterationOutcome begin
    Dominating
    Improving
    Unsuccessful
end

"""
    @enum ConstraintOutcome

Has the values `Feasible`, `WeakInfeasible`, or `StrongInfeasible` to classify
the outcome of the constraint evaluations of a single point.

A `Feasible` point meets the requirement of all constraints with no relaxation.

A `WeakInfeasible` outcome has at least one relaxable constraint violated
but no unrelaxable constraints violated.

A `StrongInfeasible` outcome indicates at least one unrelaxable constraint has
been violated or the relaxable constraint violation is reater than hmax.
"""
@enum ConstraintOutcome begin
    Feasible
    WeakInfeasible
    StrongInfeasible
end

"""
    CollectionIndex

An `Int` wrapper that is used for indexing constraint collections within a
`Constraints` object.
"""
struct CollectionIndex
    value::Int
end

"""
    CollectionIndex

An `Int` wrapper that is used for indexing constraints within a
`ConstraintCollection` object.
"""
struct ConstraintIndex
    value::Int
end

"""
    DefaultExtremeRef

The collection index that refers to the default location of extreme barrier
constraints.
"""
const DefaultExtremeRef = CollectionIndex(1)

"""
    DefaultProgressiveRef

The collection index that refers to the default location of progressive barrier
constraints.
"""
const DefaultProgressiveRef = CollectionIndex(2)

"""
    ConstraintCollection{FT,C}(h_max::FT,
                               h_max_update::Function,
                               aggregator::Function
                              ) where {FT<:AbstractFloat,C<:AbstractConstraint}

Contains multiple constraints of the same type that have the same settings
applied to them.

`h_max` is the initial hmax value.
`h_max_update` is a function that should update `h_max` given an
`IterationOutcome` value.
`aggregator` is a function that will bring all constraint violations of a
collection into a single `h` result.

Defaults for each of these values are set in the `AddProgressiveCollection` and
`AddExtremeCollection` functions.
"""
mutable struct ConstraintCollection{FT<:AbstractFloat, C<:AbstractConstraint}
    constraints::Vector{C}
    count::Int
    ignore::Bool
    h_max::FT
    result_aggregate::Function
    violation::FT

    function ConstraintCollection{FT, C}(h_max::FT, aggregator::Function
                                        ) where {FT<:AbstractFloat, C<:AbstractConstraint}
        c = new{FT, C}()
        c.constraints = []
        c.count = 0
        c.ignore = false
        c.h_max = h_max
        c.result_aggregate = aggregator
        c.violation = zero(FT)
        return c
    end
end

"""
    AbstractProgressiveConstraint <: AbstractConstraint

Parent type for progressive contraints.
"""
abstract type AbstractProgressiveConstraint <: AbstractConstraint end

"""
    ExtremeConstraint(f::Function)

Create an extreme barrier constraint. Function `f` should take a vector argument
and return `true` or `false` to indicate if the vector meets the constraint.
"""
mutable struct ExtremeConstraint <: AbstractConstraint
    f::Function
    ignore::Bool
    ExtremeConstraint(f::Function) = new(f, false)
end

"""
    ProgressiveConstraint(f::Function)

Create a progressive barrier constraint.

Argument `f` is a function that should take a single vector argument and
return a value that gives the amount the constraint function has been
violated.

A value greater than 0 indicates the function has been violated, 0 shows
that the input lies on the constraint, and negative numbers show a feasible
value.

Negative numbers may be truncated to 0 without affecting the algorithm.
"""
mutable struct ProgressiveConstraint <: AbstractProgressiveConstraint
    f::Function
    ignore::Bool
    ProgressiveConstraint(f::Function) = new(f, false)
end

"""
    ConstraintCache{FT<:AbstractFloat}

Store constraint information on an iteration-by-iteration basis.

Stores the final infeasible point hmax value of the previous iteration.
Also stores the collection hmax values computed for each point in an iteration.
"""
mutable struct ConstraintCache{FT<:AbstractFloat}
    OldHmax::FT
    hmax_map::Dict{Vector{FT}, Vector{FT}}

    function ConstraintCache{FT}() where {FT<:AbstractFloat}
        c = new{FT}()
        c.OldHmax = FT(Inf)
        c.hmax_map = Dict{Vector{FT}, Vector{FT}}()
        return c
    end
end

"""
    Constraints{FT<:AbstractFloat}()

Create an object that constains multiple `ConstraintCollection` objects and their
corresponding `ConstraintCache`.

Upon creation `Constraints` is automaticvally populated with two constraint
collections, an `ExtremeCollection` and a `ProgressiveCollection`.
"""
mutable struct Constraints{FT<:AbstractFloat}
    collections::Vector{ConstraintCollection}
    count::Int
    cache::ConstraintCache{FT}

    function Constraints{FT}() where {FT<:AbstractFloat}
        c = new{FT}()
        c.count = 0
        c.collections = []
        c.cache = ConstraintCache{FT}()

        AddExtremeCollection(c)
        AddProgressiveCollection(c)
        return c
    end
end

"""
    CollectionTypeCount(c::Constraints{T}, C::AbstractConstraint)::Int where T

Return the total number of constraints of type `C` that are stored in all collections.
"""
(CollectionTypeCount(c::Constraints{T}, C::AbstractConstraint)::Int) where T =
                sum([col.count for col in c.collections if typeof(col) ==
                     ConstraintCollection{T, C}])

"""
    ConstraintUpdate!(c::Constraints, result::IterationOutcome)

Perform all necessary actions required to update a problem's constraints
between iterations and clear the constraint cache.
"""
function UpdateConstraints(c::Constraints, h_max, result::IterationOutcome, feasible, infeasible)

    # Internal hmax should only be updated when the result is improving
    if result == Improving
        # hmax values per collection for the infeasible point's overall hmax value
        collections_hmax = c.cache.hmax_map[infeasible]

        # collection hmax values for all considered points
        all_points_hmax = values(c.cache.hmax_map)

        #Update internal hmax for each collection
        for (i, collection) in enumerate(c.collections)
            #The hmax value from the collection for each considered point
            hmax_set = [c[i] for c in all_points_hmax]

            if collections_hmax[i] == minimum(hmax_set)
                # if the infeasible point's hmax is minimum considered, then use that
                SetCollectionHmax(collection, collections_hmax[i])
            else
                # otherwise use the largest value less than the the infeasible point's
                SetCollectionHmax(collection, maximum(filter(x->x<collections_hmax[i], hmax_set)))
            end
        end
        SetOldHmax(c, h_max)
    end

    # Cache is only maintained over a single iteration
    ClearCache(c)
end

"""
    ConstraintCachePush(c::Constraints{T}, x::Vector{T}, i::Int, h::T) where T

For point `x` and constraint collection `i` push the violation function result
`h` to the cache.
"""
function ConstraintCachePush(p::AbstractProblem, x::Vector{T}, i::Int, h::T) where T
    if p.config.max_simultaneous_evaluations > 1
        lock(() -> constraint_cache_push(p.constraints, x, i, h), p.config.parallel_lock)
    else
        constraint_cache_push(p.constraints, x, i, h)
    end
end

function constraint_cache_push(c::Constraints{T}, x::Vector{T}, i::Int, h::T) where T
    if !haskey(c.cache.hmax_map, x)
        c.cache.hmax_map[x] = zeros(T, c.count)
    end
    c.cache.hmax_map[x][i] = h
end

function SetCollectionHmax(c::ConstraintCollection, new_hmax)
    c.h_max = new_hmax
end

function ClearCache(c::Constraints{T}) where T
    c.cache.hmax_map = Dict{Vector{T},Vector{T}}()
end

function SetOldHmax(c::Constraints{T}, hmax::T) where T
    c.cache.OldHmax = hmax
end

"""
    ConstraintEvaluation(constraints::Constraints{T}, p::Vector{T})::Tuple{ConstraintOutcome,T} where T

Evaluate point `p` over all constraint collections in `constraints`. Returns a
ConstraintOutcome indicating the result of the evaluation:

`Feasible`: `p` evaluated as feasible for all constraints (extreme and progressive barrier)

`WeakInfeasible`: `p` evaluated as feasible for all extreme barrier constraints, and had no
progressive barrier constraint violations greater than h_max

`StrongInfeasible`: At least one extreme barrier constraint evaluated as infeasible, or at least
one progressive barrier constraint had a violation greater than h_max

The second returned value is the sum of h_max values evaluated during the constraint checks.
"""
function ConstraintEvaluation(p::AbstractProblem, point::Vector{T})::ConstraintOutcome where T
    # Initially define result as feasible
    eval_result = Feasible
    for (i,collection) in enumerate(p.constraints.collections)
        collection.ignore && continue

        result = ConstraintCollectionEvaluation(p, collection, point)
        ConstraintCachePush(p, point, i, collection.violation)

        if result == WeakInfeasible
            eval_result = WeakInfeasible
        end

        # A single completely infeasible result (h(x) > 0 for extreme, or h(x) > h_max for prog) means invalid
        result == StrongInfeasible && return StrongInfeasible
    end

    return eval_result
end

"""
    GetViolationSum(c::Constraints, p::Vector{T})::T where T

Return the sum of all cached `h` values for point `p`.

If `p` hasn't been evaluated (which generally shouldn't happen) then
return ``\\infty``.
"""
(GetViolationSum(c::Constraints{T}, p::Vector{T})::T) where T =
                                            sum(get(c.cache.hmax_map, p, Inf))

(GetViolationSumParallel(p::AbstractProblem, point::Vector{T})::T) where T =
                                            lock(() -> GetViolationSum(p.constraints, point), p.config.parallel_lock)

(GetOldHmaxSum(c::Constraints{T})::T) where T = c.cache.OldHmax

function ConstraintCollectionEvaluation(p::AbstractProblem, collection::ConstraintCollection,
                                        x::Vector{FT})::ConstraintOutcome where {FT<:AbstractFloat}
    if p.config.max_simultaneous_evaluations > 1
        lock(() -> constraint_collection_evaluation(collection, x), p.config.parallel_lock)
    else
        constraint_collection_evaluation(collection, x)
    end
end

"""
    constraint_collection_evaluation(collection::ConstraintCollection{T,ProgressiveConstraint},
                                   x::Vector{T})::ConstraintOutcome where T

Evalute every constraint within progressive constraint collection `collection`
for point `x`.

If the aggregate value of the constraint evaluations exceeds the collection's
h_max then a `StrongInfeasible` is returned. If the value is less than or equal
to 0.0 then `Feasible` is returned. Otherwise `WeakInfeasible` is returned.
"""
function constraint_collection_evaluation(collection::ConstraintCollection{T,ProgressiveConstraint},
                                          x::Vector{T})::ConstraintOutcome where T
    sum = 0.0
    for c in collection.constraints
        c.ignore && continue
        sum += collection.result_aggregate(c.f(x))
    end
    collection.violation = sum
    if sum <= collection.h_max
        return sum == 0.0 ? Feasible : WeakInfeasible
    else
        return StrongInfeasible
    end
end

"""
    constraint_collection_evaluation(collection::ConstraintCollection{FT,ExtremeConstraint},
                                   x::Vector{FT})::ConstraintOutcome where {FT<:AbstractFloat}

Evalute every constraint within extreme constraint collection `collection` for
point `x`.

If any constraint returns false  or a value greater than 0 then a
`StrongInfeasible` result is returned. Otherwise a `Feasible` result is returned.
"""
function constraint_collection_evaluation(collection::ConstraintCollection{FT,ExtremeConstraint},
                                          x::Vector{FT})::ConstraintOutcome where {FT<:AbstractFloat}
    for c in collection.constraints
        c.ignore && continue
        v = c.f(x)
        v < zero(FT) || v == true || return StrongInfeasible
    end

    return Feasible
end

"""
    AddExtremeConstraint(p::AbstractProblem, c::Function
                         index::CollectionIndex=CollectionIndex(1)
                        )::ConstraintIndex where T

Register a single function that defines an extreme barrier constraint. Return
a constraint index.

The provided function should take a vector input and return a boolean or numeric
value indicating if the constraint has been met or not. `true` or less than or
equal to 0 indicates the constraint has been met. `false` or greater than 0
shows the constraint has been violated.

The `index` argument can be specified to give a collection to add the constraint
to. The specified collection must exist, and must be able to accept extreme
barrier constraints. If `index` is not specified then it is added to collection
1, the default extreme constraint collection.
"""
AddExtremeConstraint(p::AbstractProblem, f::Function; index::CollectionIndex=CollectionIndex(1)
                    ) = AddExtremeConstraint(p.constraints, f, index=index)

function AddExtremeConstraint(p::Constraints{T}, f::Function;
                              index::CollectionIndex=CollectionIndex(1))::ConstraintIndex where T
    i = index.value
    i < 0 && error("Collection indices must be positive")
    i > p.count && error("Invalid CollectionIndex")

    c = ExtremeConstraint(f)
    push!(p.collections[i].constraints, c)
    p.collections[i].count += 1
    return ConstraintIndex(p.collections[i].count)
end

"""
    AddExtremeConstraint(p::AbstractProblem, c::Vector{Function}; index::CollectionIndex=CollectionIndex(1))

Register a group of functions that define extreme barrier constraints. Calls
[`AddExtremeConstraint`](@ref) on each function individually.
"""
AddExtremeConstraint(p::AbstractProblem, f::Vector{Function}; index::CollectionIndex=CollectionIndex(1)
                    ) = AddExtremeConstraint(p.constraints, f, index=index)

function AddExtremeConstraint(p::Constraints, f::Vector{Function};
                              index::CollectionIndex=CollectionIndex(1))::Vector{ConstraintIndex}
    return [AddExtremeConstraint(p, fnc, index=index) for fnc in f]
end

"""
    AddProgressiveConstraint(p::AbstractProblem, c::Function; index::CollectionIndex=CollectionIndex(2))

Register a single function that defines a progressive barrier constraint. Return
an index that refers to the constraint.

The provided function should take a vector input and return a numeric value
indicating if the constraint has been met or not. Less than or
equal to 0 indicates the constraint has been met. 0 shows the constraint has
been violated.

The `index` argument can be specified to give a collection to add the constraint
to. The specified collection must exist, and must be able to accept progressive
barrier constraints. If `index` is not specified then it is added to collection
2, the default progressive barrier constraint collection.
"""
AddProgressiveConstraint(p::AbstractProblem, f::Function; index::CollectionIndex=CollectionIndex(2)
                        ) = AddProgressiveConstraint(p.constraints, f, index=index)

function AddProgressiveConstraint(p::Constraints, f::Function;
                                  index::CollectionIndex=CollectionIndex(2))::ConstraintIndex
    i = index.value
    i < 0 && error("Collection indices must be positive")
    i > p.count && error("Invalid CollectionIndex")

    c = ProgressiveConstraint(f)
    push!(p.collections[i].constraints, c)
    p.collections[i].count += 1
    return ConstraintIndex(p.collections[i].count)
end

"""
    AddProgressiveConstraint(p::AbstractProblem, c::Vector{Function})::Vector{Int}

Register a group of functions that define progressive barrier constraints. Calls
[`AddProgressiveConstraint`](@ref) on each function individually.
"""
AddProgressiveConstraint(p::AbstractProblem, f::Vector{Function}; index::CollectionIndex=CollectionIndex(2)
                        ) = AddProgressiveConstraint(p.constraints, f, index=index)

function AddProgressiveConstraint(p::Constraints, c::Vector{Function};
                                  index::CollectionIndex=CollectionIndex(2))::Vector{ConstraintIndex}
    return map(f -> AddProgressiveConstraint(p, f, index=index), c)
end

"""
    AddProgressiveCollection(p::Constraints{FT}; h_max=Inf, h_max_update::Function=h_max_update,
                             aggregator::Function=x->max(0,x)^2)::CollectionIndex where {FT<:AbstractFloat}

Instantiate a new constraint collection within the problem. Returns an index that refers to this
new collection.

The default constraint settings match those from Audet & Dennis 2009:

`h_max`: Begins as infinity

`h_max_update`: Sets h_max to the largest valid h evaluation if an iteration is improving

`aggregator`: Creates h as ``\\sum k(x)`` where ``k=\\max(0,x)^2``
"""
AddProgressiveCollection(p::AbstractProblem; kwargs...)::CollectionIndex = AddProgressiveCollection(p.constraints; kwargs...)

function AddProgressiveCollection(p::Constraints{FT}; h_max=FT(Inf),
                                  aggregator::Function=x->max(0,x)^2)::CollectionIndex where {FT<:AbstractFloat}
    push!(p.collections,
          ConstraintCollection{FT,ProgressiveConstraint}(h_max, aggregator))
    p.count += 1

    return CollectionIndex(p.count)
end

"""
    AddExtremeCollection(p::Constraints{FT})::CollectionIndex where {FT<:AbstractFloat}

Instantiate a new constraint collection for extreme constraints. Returns an index that
refers to the new collection.
"""
(AddExtremeCollection(p::AbstractProblem)::CollectionIndex) = AddExtremeCollection(p.constraints)

function AddExtremeCollection(p::Constraints{FT})::CollectionIndex where {FT<:AbstractFloat}
    f() = error("Should not be calling aggregator update for an extreme constraint collection")
    push!(p.collections, ConstraintCollection{FT,ExtremeConstraint}(zero(FT), f))
    p.count += 1
    return CollectionIndex(p.count)
end

#DisableConstraint(p::AbstractProblem, constraint_index::ConstraintIndex;
#                  collection_index::CollectionIndex = CollectionIndex(
#                 ) = DisableConstraint(p.constraints, i)
#function DisableConstraint(c::Constraints, i::Int)
#    i > c.count && error("Invalid constraint index, largest known index is $(c.count)")
#    c.constraints[i].ignore == true
#end
#
#EnableConstraint(p::AbstractProblem, i::Int) = EnableConstraint(p.constraints, i)
#function EnableConstraint(c::ConstraintCollection, i::Int)
#    i > c.count && error("Invalid constraint index, largest known index is $(c.count)")
#    c.constraints[i].ignore == false
#end


text:
# If an integer is passed down to the function, it means that the user is requesting an exact number of iterations.
function bisection_it(f,ansatz::Tuple{Float64,Float64}; N::Int64=100)
    (lower,upper) = ansatz
    mid = 0
    for _ in 1:N

    mid = (lower + upper)/2

        if f(lower)*f(mid) < 0
            upper = mid
        elseif f(lower)*f(mid) > 0
            lower = mid
        else
            return mid
        end

    end
    return mid
end

# f(x) = sin(x)
# bisection(f,(2.,4.))

#If then a Float is passed to the function, then it means it must mean an objective relative error.
function bisection_err(f,ansatz::Tuple{Float64,Float64}; error::Float64=0.01)
    mid = (lower + upper)/2;
    ea = 1; #Initialize
    #We do the first iteration, otherwise, there is no way to have a relative error to measure
    if f(lower)*f(mid) < 0
        upper = mid
    elseif f(lower)*f(mid) > 0
        lower = mid
    else
        return mid
    end
    #Now we start doing the iterations until it is done
    while ea > error
        last = mid #We store the last approximation
        mid = (lower + upper)/2
        ea = abs((mid-last)/mid);

        if f(lower)*f(mid) < 0
            upper = mid
        elseif f(lower)*f(mid) > 0
            lower = mid
        else
            return mid
        end

    end
    return mid
end


text:
using JlrsReflect
using Test
include("SingleFieldBits.jl")
include("MultiFieldBits.jl")
include("BitsWithCustom.jl")
include("BitsWithTuples.jl")
include("WithBuiltinFields.jl")
include("WithGenericFields.jl")
include("WithBitsUnion.jl")
include("WithNonBitsUnion.jl")
include("ZeroSized.jl")


text:
@testset "Random trees" begin
	tree = randtree(2, true)
	@test n_tip(tree) == 2
	@test n_node(tree) == 3
	@test isrooted(tree)
	tree = randtree(3)
	@test n_tip(tree) == 3
	@test n_node(tree) == 4
	@test ! isrooted(tree)
	tree = randtree(3, true)
	@test n_tip(tree) == 3
	@test n_node(tree) == 5
	@test isrooted(tree)
	for _ in 1:100
		tree = randtree(50)
		@test n_tip(tree) == 50
		@test n_node(tree) == 98
		@test ! isrooted(tree)
	end
	for _ in 1:100
		tree = randtree(50, true)
		@test n_tip(tree) == 50
		@test n_node(tree) == 99
		@test isrooted(tree)
	end
	tree = randtree(string.(1:34))
	@test n_tip(tree) == 34
	@test n_node(tree) == 66
	tree = randtree(SpeciesDirectory(string.(1:19)))
	@test n_tip(tree) == 19
	@test n_node(tree) == 36
end


text:
import Base.string
import Base.print
import Base.show

export GenieType, GenieController, Controller

abstract type GenieType end
string(io::IO, t::T) where {T<:GenieType} = "$(typeof(t)) <: $(super(typeof(t)))"
print(io::IO, t::T) where {T<:GenieType} = print(io, "$(typeof(t)) <: $(super(typeof(t)))")
show(io::IO, t::T) where {T<:GenieType} = print(io, genietype_to_print(t))

mutable struct GenieController <: GenieType
end

const Controller = GenieController

mutable struct GenieChannel <: GenieType
end

# const Channel = GenieChannel

import Millboard

"""
    genietype_to_print{T<:GenieType}(m::T) :: String

Pretty printing of Genie types.
"""
function genietype_to_print(m::T) :: String where {T<:GenieType}
  output = "\n" * "$(typeof(m))" * "\n"
  output *= string(config.log_formatted ? Millboard.table(to_string_dict(m)) : to_string_dict(m) ) * "\n"

  output
end

"""
    to_dict(m::Any) :: Dict{String,Any}
    to_string_dict(m::Any; all_output::Bool = false) :: Dict{String,String}
    to_string_dict(m::Any, fields::Array{Symbol,1}; all_output::Bool = false) :: Dict{String,String}

Creates a `Dict` using the fields and the values of `m`.
"""
function to_dict(m::Any) :: Dict{String,Any}
  Dict(string(f) => getfield(m, Symbol(f)) for f in fieldnames(typeof(m)))
end
function to_string_dict(m::Any; all_output::Bool = false) :: Dict{String,String}
  to_string_dict(m, fieldnames(m), all_output = all_output)
end
function to_string_dict(m::Any, fields::Array{Symbol,1}; all_output::Bool = false) :: Dict{String,String}
  response = Dict{String,String}()
  for f in fields
    response[string(f)] = string(getfield(m, Symbol(f)))
  end

  response
end


text:
using Statistics

export blackScholesFormula

"Price of a vanilla option under the Black-Scholes model"
function blackScholesFormula(isCall::Bool, strike::T, spot::T, variance::T, driftDf::T, discountDf::T) where {T}

    sign = 1.0
    if !isCall
        sign = -1.0
    end
    forward = spot / driftDf
    if variance < eps()
        price = discountDf * max(sign * (forward - strike), 0.0)
        return price
    elseif spot < eps()
        if isCall
            return 0.0
        else
            return discountDf * strike
        end
    elseif strike < eps()
        if isCall
            return discountDf * forward
        else
            return 0.0
        end
    else
        sqrtVar = sqrt(variance)
        d1 = 1.0 / sqrtVar * log(forward / strike) + 0.5 * sqrtVar
        d2 = d1 - sqrtVar
        normal = Normal()
        nd1 = cdf(normal, sign * d1)
        nd2 = cdf(normal, sign * d2)
        price = sign * discountDf * (forward * nd1 - strike * nd2)
        return price
    end
end


text:
module Datasets
using Random
using DelimitedFiles: readdlm

include("datasets.jl")
export load_iris

include("generators.jl")
export make_spiral,
       make_poly,
       make_sin

end


text:
function add_additional_state_provider(obj::AbstractPropagator, arg0::AdditionalStateProvider)
    return jcall(obj, "addAdditionalStateProvider", void, (AdditionalStateProvider,), arg0)
end

function add_event_detector(obj::AbstractAnalyticalPropagator, arg0::EventDetector)
    return jcall(obj, "addEventDetector", void, (EventDetector,), arg0)
end

function clear_events_detectors(obj::AbstractAnalyticalPropagator)
    return jcall(obj, "clearEventsDetectors", void, ())
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_additional_state_providers(obj::AbstractPropagator)
    return jcall(obj, "getAdditionalStateProviders", List, ())
end

function get_attitude_provider(obj::AbstractPropagator)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_events_detectors(obj::AbstractAnalyticalPropagator)
    return jcall(obj, "getEventsDetectors", Collection, ())
end

function get_frame(obj::AbstractPropagator)
    return jcall(obj, "getFrame", Frame, ())
end

function get_generated_ephemeris(obj::AbstractAnalyticalPropagator)
    return jcall(obj, "getGeneratedEphemeris", BoundedPropagator, ())
end

function get_initial_state(obj::EphemerisSegmentPropagator)
    return jcall(obj, "getInitialState", SpacecraftState, ())
end

function get_managed_additional_states(obj::AbstractPropagator)
    return jcall(obj, "getManagedAdditionalStates", Vector{JString}, ())
end

function get_max_date(obj::EphemerisSegmentPropagator)
    return jcall(obj, "getMaxDate", AbsoluteDate, ())
end

function get_min_date(obj::EphemerisSegmentPropagator)
    return jcall(obj, "getMinDate", AbsoluteDate, ())
end

function get_mode(obj::AbstractPropagator)
    return jcall(obj, "getMode", jint, ())
end

function get_pv_coordinates(obj::EphemerisSegmentPropagator, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function get_pv_provider(obj::AbstractAnalyticalPropagator)
    return jcall(obj, "getPvProvider", PVCoordinatesProvider, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function is_additional_state_managed(obj::AbstractPropagator, arg0::JString)
    return jcall(obj, "isAdditionalStateManaged", jboolean, (JString,), arg0)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function propagate(obj::AbstractAnalyticalPropagator, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "propagate", SpacecraftState, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end

function propagate(obj::AbstractPropagator, arg0::AbsoluteDate)
    return jcall(obj, "propagate", SpacecraftState, (AbsoluteDate,), arg0)
end

function reset_initial_state(obj::EphemerisSegmentPropagator, arg0::SpacecraftState)
    return jcall(obj, "resetInitialState", void, (SpacecraftState,), arg0)
end

function set_attitude_provider(obj::AbstractPropagator, arg0::AttitudeProvider)
    return jcall(obj, "setAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function set_ephemeris_mode(obj::AbstractPropagator)
    return jcall(obj, "setEphemerisMode", void, ())
end

function set_ephemeris_mode(obj::AbstractPropagator, arg0::OrekitStepHandler)
    return jcall(obj, "setEphemerisMode", void, (OrekitStepHandler,), arg0)
end

function set_master_mode(obj::AbstractPropagator, arg0::OrekitStepHandler)
    return jcall(obj, "setMasterMode", void, (OrekitStepHandler,), arg0)
end

function set_master_mode(obj::AbstractPropagator, arg0::jdouble, arg1::OrekitFixedStepHandler)
    return jcall(obj, "setMasterMode", void, (jdouble, OrekitFixedStepHandler), arg0, arg1)
end

function set_slave_mode(obj::AbstractPropagator)
    return jcall(obj, "setSlaveMode", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using ControlSystems

"""
    ntf = synthesizeNTF(order=3, osr=64, opt=0, H_inf=1.5, f0=0)

Synthesize a noise transfer function for a delta-sigma modulator.
- `order`: order of the modulator
- `osr`: oversampling ratio
- `opt`: flag for optimized zeros
    - 0 -> not optimized
    - 1 -> optimized
    - 2 -> optimized with at least one zero at band-center
    - [z] -> zero locations in complex form
- `H_inf`: maximum NTF gain
- `f0`: center frequency (1->fs)

`ntf` is a TransferFunction object containing the zeros and poles of the NTF.
"""
function synthesizeNTF(order=3, osr=64.0, opt=0, H_inf=1.5, f0=0.0)
    if f0 > 0.5
        error("f0 must be less than 0.5")
    end
    if f0 != 0 && f0 < (0.25/osr)
        @warn "Creating a lowpass ntf"
        f0 = 0
    end
    if f0 != 0 && isodd(order)
        error("order must be even for a bandpass modulator")
    end
    if length(opt) > 1 && length(opt) != order
        error("The opt vector must be of length $order(=order)")
    end

    # Determine the zeros
    if f0 != 0  # Bandpass design-- halve the order temporarily
        order = order ÷ 2
        dw = pi / (2*osr)
    else
        dw = pi / osr
    end

    if length(opt) == 1
        if opt == 0
            z = zeros(order)
        else
            z = dw * ds_optzeros(order, opt)
            if isempty(z)
                return
            end
            if f0 != 0  # Bandpass design -- shift and replicate the zeros.
                order = order * 2
                z = z .+ 2*pi*f0
                z = [z'; -z'][:]
            end
        end
        z = exp.(im * z)
    else
        z = opt[:]
    end

    p = zeros(order)
    ntf = zpk(z, p, 1, 1)
    itn_limit = 100

    delta_x = 0
    fprev = 0

    # Iteratively determine the poles by finding the value of the x-parameter
    # which results in the desired H_inf.
    if f0 == 0  # Lowpass design
        HinfLimit = 2^order
        if H_inf >= HinfLimit
            @warn "Unable to achieve specified Hinf. Setting all NTF poles to zero."
            p = zeros(order)
        else
            x = 0.3 ^ (order-1)   # starting guess
            for itn = 1:itn_limit
                me2 = -0.5 * x^(2/order)
                w = (2*(1:order) .- 1) * pi / order
                mb2 = 1 .+ me2 * exp.(im * w)
                p = mb2 - sqrt.(mb2.^2 .- 1)
                out = findall(x -> abs(x) > 1, p)
                p[out] = 1 ./ p[out] # reflect poles to be inside the unit circle
                p = pairComplex(p)

                ntf = zpk(z, p, 1, 1)
                f = real(evalTF(ntf, -1)) - H_inf
                if itn == 1
                    delta_x = -f / 100
                else
                    delta_x = -f * delta_x / (f - fprev)
                end

                xplus = x + delta_x
                if xplus > 0
                    x = xplus
                else
                    x = x * 0.1
                end

                fprev = f
                if abs(f) < 1e-10 || abs(delta_x) < 1e-10
                    break
                end

                if x > 1e6
                    @warn "Unable to achieve specified Hinf. Setting all NTF poles to zero."
                    p = zeros(order)
                    break
                end

                if itn == itn_limit
                    @warn "Danger! Iteration limit exceeded."
                end
            end
        end
    else    # Bandpass design
        x = 0.3 ^ (order/2 - 1)   # starting guess (not very good for f0~0)
        if f0 > 0.25
            z_inf = 1
        else
            z_inf = -1
        end
        c2pif0 = cos(2 * pi * f0)

        for itn = 1:itn_limit
            e2 = 0.5 * x^(2/order)
            w = (2*(1:order) .- 1) * pi / order
            mb2 = c2pif0 .+ e2 * exp.(im * w)
            p = mb2 - sqrt.(mb2.^2 .- 1)
            out = findall(x -> abs(x) > 1, p)
            p[out] = 1 ./ p[out]    # reflect poles to be inside the unit circle.
            p = pairComplex(p)

            ntf = zpk(z, p, 1, 1)
            f = real(evalTF(ntf, z_inf)) - H_inf
            if itn == 1
                delta_x = -f / 100
            else
                delta_x = -f * delta_x / (f - fprev)
            end

            xplus = x + delta_x
            if xplus > 0
                x = xplus
            else
                x = x * 0.1
            end

            fprev = f
            if abs(f) < 1e-10 || abs(delta_x) < 1e-10
                break
            end

            if x > 1e6
                @warn "Unable to achieve specified Hinf. Setting all NTF poles to zero."
                p = zeros(order)
                break
            end

            if itn == itn_limit
                @warn "Danger! Iteration limit exceeded."
            end
        end
    end

    ntf = zpk(z, p, 1, 1)
    return ntf
end

"""
    pairComplex(a)

A helper function to pair up and order complex conjugate pairs in `a`
so it can be passed to `zpk`(). Based off MATLAB cplxpair().
"""
function pairComplex(a)
    a = sort(a, by=x->(real(x), imag(x)), rev=true)
    rtol = 100 * eps(Float64) # tolerance
    i = 1
    while i <= length(a)
        if isapprox(a[i], real(a[i]), rtol=rtol)
            # imaginary component is zero or very small, treat as purely real
            a[i] = real(a[i])
        else
            if isapprox(a[i], a[i+1]', rtol=rtol)
                # complex conjugate pairs with some error, turn to exact pairs
                # by setting real and imaginary to midpoint of the two
                a_real = (real(a[i]) + real(a[i+1])) / 2
                a_imag = (abs(imag(a[i])) + abs(imag(a[i+1]))) / 2
                a[i] = a_real + a_imag*im   # put positive pair first
                a[i+1] = a_real - a_imag*im
            else
                error("Complex numbers can't be paired")
            end
            i += 1
        end
        i += 1
    end
    return a
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
function load_celestial_body(obj::CelestialBodyLoader, arg0::JString)
    return jcall(obj, "loadCelestialBody", CelestialBody, (JString,), arg0)
end


text:
using SumOfSquares
using Documenter, Literate

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")

include(joinpath(EXAMPLES_DIR, "all_examples.jl"))
deleteat!(EXAMPLES, findfirst(isequal("chordal_sparsity.jl"), EXAMPLES))
deleteat!(EXAMPLES, findfirst(isequal("chordal_sparsity_with_domain.jl"), EXAMPLES))

for example in EXAMPLES
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR)
    Literate.notebook(example_filepath, OUTPUT_DIR)
    Literate.script(example_filepath, OUTPUT_DIR)
end

makedocs(
    sitename = "SumOfSquares",
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    # See https://github.com/jump-dev/JuMP.jl/issues/1576
    strict = true,
    pages = [
        "Index" => "index.md",
        "Sum-of-Squares Programming" => "sumofsquares.md",
        "Variables" => "variables.md",
        "Constraints" => "constraints.md",
        "Examples" => map(EXAMPLES) do jl_file
            # Need `string` as Documenter fails if `name` is a `SubString{String}`.
            name = string(split(jl_file, ".")[1])
            return name => "generated/$name.md"
        end
    ],
    # The following ensures that we only include the docstrings from
    # this module for functions define in Base that we overwrite.
    modules = [SumOfSquares, PolyJuMP]
)

deploydocs(
    repo   = "github.com/jump-dev/SumOfSquares.jl.git",
    push_preview = true,
)


text:
#included within module Time

export  DayCount, Actual360, Actual365, Thirty360, BondThirty360, EuroBondThirty360, ItalianThirty360, ISMAActualActual, ISDAActualActual, AFBActualActual,
		daycount, yearfraction

abstract DayCount

type Actual360 <:DayCount ; end
type Actual365 <: DayCount ; end
abstract Thirty360 <:DayCount

type BondThirty360 <: Thirty360; end
type EuroBondThirty360 <: Thirty360; end
type ItalianThirty360 <: Thirty360; end

typealias USAThirty360 BondThirty360
typealias EuroThirty360 EuroBondThirty360

abstract ActualActual <:DayCount

type ISMAActualActual <: ActualActual; end
type ISDAActualActual <: ActualActual; end
type AFBActualActual <: ActualActual; end

type BusinessDay252 <: DayCount
	bc::BusinessCalendar
end

string(::Actual360) = "Actual/360"
string(::Actual365) = "Actual/365 Fixed"
string(::BondThirty360) = "Thirty/360 Bond/USA basis"
string(::EuroBondThirty360) = "Thirty/360 EuroBond basis"
string(::ItalianThirty360) = "Thirty/360 Italian basis"

string(::ISMAActualActual) = "Actual/Actual ISMA basis"
string(::ISDAActualActual) = "Actual/Actual ISDA basis"
string(::AFBActualActual) = "Actual/Actual AFB basis"

#Assumes both CalendarTime objects represent midnight in the same tz.
#The rounding takes care of summer time and leap second
#Will fail if they represent arbitrary times of day
daycount(c::DayCount,  d_start::CalendarTime,  d_end::CalendarTime) = iround((d_end-d_start).millis/86400e3)

function daycount(c::BondThirty360,  d_start::CalendarTime,  d_end::CalendarTime)
	dd1 = day(d_start)
     dd2 = day(d_end)
     mm1 = month(d_start)
     mm2 = month(d_end)
     yy1 = year(d_start)
     yy2 = year(d_end)

    if dd2 == 31 && dd1 < 30
    	dd2 = 1
    	mm2 = mm2+1
    end

    return 360*(yy2-yy1) + 30*(mm2-mm1-1) + max(0, 30-dd1) + min(30, dd2)
end

function daycount(c::EuroBondThirty360,  d_start::CalendarTime,  d_end::CalendarTime)
	dd1 = day(d_start)
    dd2 = day(d_end)
    mm1 = month(d_start)
    mm2 = month(d_end)
    yy1 = year(d_start)
    yy2 = year(d_end)

   	return 360*(yy2-yy1) + 30*(mm2-mm1-1) + max(0, 30-dd1) + min(30, dd2)
end

function daycount(c::ItalianThirty360,  d_start::CalendarTime,  d_end::CalendarTime)
	dd1 = day(d_start)
    dd2 = day(d_end)
    mm1 = month(d_start)
    mm2 = month(d_end)
    yy1 = year(d_start)
    yy2 = year(d_end)

    if (mm1 == 2 && dd1 > 27)
        dd1 = 30
    end
    if (mm2 == 2 && dd2 > 27)
        dd2 = 30
    end

    return 360*(yy2-yy1) + 30*(mm2-mm1-1) + max(0, 30-dd1) + min(30, dd2)
end

daycount(c::BusinessDay252, d_start::CalendarTime, d_end::CalendarTime) = businessDaysBetween(c.bc, d_start, d_end)

yearfraction(c::Union(Actual360,Thirty360),  d_start::CalendarTime,  d_end::CalendarTime) = daycount(c, d_start, d_end) / 360

yearfraction(c::Actual365,  d_start::CalendarTime,  d_end::CalendarTime) = daycount(c, d_start, d_end) / 365

function yearfraction(c::ISDAActualActual,  d_start::CalendarTime,  d_end::CalendarTime )

    if (d_start == d_end ); return 0.0; end

    if (d_start > d_end)
        return -yearFraction(d_end, d_start)
    end

    y1 = year(d_start)
    y2 = year(d_end)
    dib1 = isleapyear(d_start) ? 366.0 : 365.0
    dib2 = isleapyear(d_end)   ? 366.0 : 365.0

    sum = y2 - y1 - 1

    #Days from start to starting of following year
    sum += (dib1 - dayofyear(d_start) + 1) / dib1
    #Days from beginning of year to the endDate
    sum += (dayofyear(d_end) - 1) / dib2
    return sum
end

function yearfraction(c::AFBActualActual,  d_start::CalendarTime,  d_end::CalendarTime )
	if (d_start == d_end ); return 0.0; end

    if (d_start > d_end)
        return -yearFraction(d_end, d_start)
    end

    newD2 = d_end
	temp = d_end
	sum = 0.0
	while (temp>d_start)
	    temp = newD2+years(-1)
	    if (day(temp) == 28 && month(temp) == 2 && isleapyear(temp))
	        temp + days(1)
	    end
	    if (temp>=(d_start))
	        sum += 1.0
	        newD2 = temp
	    end
	end

	den = 365.0

	if isleapyear(newD2)
	    if (newD2 > ymd(year(newD2), February, 29) && d_start < ymd(year(newD2), February, 29) )
	        den += 1.0
	    end
	elseif (isleapyear(d_start))
	    if (newD2 > (ymd(year(d_start), February, 29)) && d_start < (ymd(year(d_start), February, 29)))
	        den += 1.0
	    end
	end
	return sum + daycount(c, d_start, newD2) / den
end

yearfraction(c::ISMAActualActual,  d_start::CalendarTime,  d_end::CalendarTime ) = yearfraction(c, d_start, d_end, d_start, d_end)

function yearfraction(c::ISMAActualActual,  d_start::CalendarTime,  d_end::CalendarTime,  ref_start::CalendarTime,  ref_end::CalendarTime)

 	if (d_start == d_end ); return 0.0; end

    if (d_start > d_end)
        return -yearFraction(d_end, d_start, ref_start, ref_end)
    end

    ref_period_start =  ref_start
    ref_period_end   = ref_end

    @assert ref_period_end > ref_period_start && ref_period_end > d_start

    # estimate roughly the length in months of a period
    mnths = ifloor (0.5 + 12 * daycount(c,  ref_period_start, ref_period_end) / 365)

    # for short periods, take the reference period as 1 year from d_start
    if (mnths == 0)
        ref_period_start = d_start
        ref_period_end = d_start + year(1)
        mnths = 12
    end

    period = mnths / 12.0

    if (d_end<=ref_period_end)
        # here ref_period_end is a future (notional?) payment date
        if (d_start>=ref_period_start)
            # here ref_period_start is the last (maybe notional)
            # payment date.
            # ref_period_start <= d_start <= d_end <= ref_period_end
            # [maybe the equality should be enforced, since
            # ref_period_start < d_start <= d_end < ref_period_end
            # could give wrong results] ???
            return period * daycount(c,d_start, d_end) / daycount(c,ref_period_start, ref_period_end)
        else
            # here ref_period_start is the next (maybe notional)
            # payment date and ref_period_end is the second next
            # (maybe notional) payment date.
            # d_start < ref_period_start < ref_period_end
            # AND d_end <= ref_period_end
            # this case is long first coupon

            # the last notional payment date
            previousRef = ref_period_start + months(-mnths)
            if (d_end>ref_period_start)
                return yearfraction(c, d_start, ref_period_start, previousRef,ref_period_start) + yearfraction(c,ref_period_start, d_end, ref_period_start, ref_period_end)
            else
                return yearfraction(d_start, d_end, previousRef,ref_period_start)
            end
        end
    else
        # here ref_period_end is the last notional payment date
        # d_start < ref_period_end < d_end AND ref_period_start < ref_period_end
        @assert ref_period_start<=d_start
        # now it is: ref_period_start <= d_start < ref_period_end < d_end

        # the part from d_start to ref_period_end
         sum = yearfraction(c,d_start, ref_period_end, ref_period_start, ref_period_end)

        # the part from ref_period_end to d_end
        # count how many regular periods are in [ref_period_end, d_end],
        # then add the remaining time
        i = 0
        local newRefStart, newRefEnd
        while true
            newRefStart = ref_period_end + months(months * i)
            newRefEnd = ref_period_end+ months(months * (i + 1))
            if (d_end<newRefEnd)
                break
            else
                sum += period
                i += 1
            end
        end
        sum += yearFraction(c, newRefStart, d_end, newRefStart, newRefEnd)
        return sum
    end
end

yearfraction(c::BusinessDay252, d_start::CalendarTime, d_end::CalendarTime) = daycount(c, d_start, d_end) / 252


text:
@test IntInf == IntInf
@test IntInf != IntSup
@test IntInf == Float64(-Inf)
@test IntInf != 1
@test IntSup == IntSup
@test IntSup != IntInf
@test IntSup == BigFloat(Inf)
@test IntSup != 1

@test isless(IntInf, IntInf) == false
@test isless(IntInf, 0)
@test isless(IntInf, IntSup)
@test isless(IntInf, -Inf) == false
@test isless(IntInf, 0.0)
@test isless(IntInf, Inf)
@test isless(0, IntInf) == false
@test isless(IntSup, IntInf) == false
@test isless(-Inf, IntInf) == false
@test isless(0.0, IntInf) == false
@test isless(Inf, IntInf) == false

@test isless(IntSup, IntSup) == false
@test isless(IntSup, 0) == false
@test isless(IntSup, IntInf) == false
@test isless(IntSup, -Inf) == false
@test isless(IntSup, 0.0) == false
@test isless(IntSup, Inf) == false
@test isless(0, IntSup)
@test isless(-Inf, IntSup)
@test isless(0.0, IntSup)
@test isless(Inf, IntSup) == false

@test IntInf == typemin(Rational{Int})
@test IntSup == typemax(Rational{BigInt})
@test IntInf == typeinf(Rational)
@test IntSup == typesup(Rational)

@test IntInf + IntInf == IntInf
@test IntInf + 0 == IntInf
@test_throws ArgumentError IntInf + IntSup
@test 0 + IntInf == IntInf
@test_throws ArgumentError IntInf - IntInf
@test IntInf - 0 == IntInf
@test IntInf - IntSup == IntInf

@test IntSup + IntSup == IntSup
@test IntSup + 0 == IntSup
@test_throws ArgumentError IntSup + IntInf
@test 0 + IntSup == IntSup
@test_throws ArgumentError IntSup - IntSup
@test IntSup - 0 == IntSup
@test IntSup - IntInf == IntSup

@test -IntInf == IntSup
@test -IntSup == IntInf


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
const default_buildings = [
    "FastFoodRest",
    "FullServiceRest",
    "Hospital",
    "LargeHotel",
    "LargeOffice",
    "MediumOffice",
    "MidriseApartment",
    "Outpatient",
    "PrimarySchool",
    "RetailStore",
    "SecondarySchool",
    "SmallHotel",
    "SmallOffice",
    "StripMall",
    "Supermarket",
    "Warehouse",
    "FlatLoad",
]

const MMBTU_TO_KWH = 293.07107

function find_ashrae_zone_city(lat, lon)::String
    file_path = joinpath(dirname(@__FILE__), "..", "..", "data", "climate_cities.shp")
    table = Shapefile.Table(file_path)
    geoms = Shapefile.shapes(table)
    # TODO following for loop is relatively slow
    for (row, geo) in enumerate(geoms)
        g = length(geo.points)
        nodes = zeros(g, 2)
        edges = zeros(g, 2)
        for (i,p) in enumerate(geo.points)
            nodes[i,:] = [p.x, p.y]
            edges[i,:] = [i, i+1]
        end
        edges[g, :] = [g, 1]
        edges = convert(Array{Int64,2}, edges)
        # shapefiles have longitude as x, latitude as y
        if inpoly2([lon, lat], nodes, edges)[1]
            return table.city[row]
        end
        GC.gc()
    end
    @info "Could not find latitude/longitude in U.S. Using geometrically nearest city."
    cities = [
        (city="Miami", lat=25.761680, lon=-80.191790),
        (city="Houston", lat=29.760427, lon=-95.369803),
        (city="Phoenix", lat=33.448377, lon=-112.074037),
        (city="Atlanta", lat=33.748995, lon=-84.387982),
        (city="LasVegas", lat=36.1699, lon=-115.1398),
        (city="LosAngeles", lat=34.052234, lon=-118.243685),
        (city="SanFrancisco", lat=37.3382, lon=-121.8863),
        (city="Baltimore", lat=39.290385, lon=-76.612189),
        (city="Albuquerque", lat=35.085334, lon=-106.605553),
        (city="Seattle", lat=47.606209, lon=-122.332071),
        (city="Chicago", lat=41.878114, lon=-87.629798),
        (city="Boulder", lat=40.014986, lon=-105.270546),
        (city="Minneapolis", lat=44.977753, lon=-93.265011),
        (city="Helena", lat=46.588371, lon=-112.024505,),
        (city="Duluth", lat=46.786672, lon=-92.100485),
        (city="Fairbanks", lat=59.0397, lon=-158.4575),
    ]
    min_distance = 0.0
    nearest_city = ""
    for (i, c) in enumerate(cities)
        distance = sqrt((lat - c.lat)^2 + (lon - c.lon)^2)
        if i == 1
            min_distance = distance
            nearest_city = c.city
        elseif distance < min_distance
            min_distance = distance
            nearest_city = c.city
        end
    end
    return nearest_city
end

"""
    built_in_load(type::String, city::String, buildingtype::String,
        year::Int, annual_energy::Real, monthly_energies::AbstractArray{Real,1}
    )
Scale a normalized Commercial Reference Building according to inputs provided and return the 8760.
"""
function built_in_load(type::String, city::String, buildingtype::String,
    year::Int, annual_energy::Real, monthly_energies::AbstractArray{Real,1}
    )

    @assert type in ["electric", "domestic_hot_water", "space_heating"]
    monthly_scalers = ones(12)
    lib_path = joinpath(dirname(@__FILE__), "..", "..", "data", "load_profiles", type)

    profile_path = joinpath(lib_path, string("crb8760_norm_" * city * "_" * buildingtype * ".dat"))
    normalized_profile = vec(readdlm(profile_path, '\n', Float64, '\n'))

    if length(monthly_energies) == 12
        annual_energy = 1.0  # do not scale based on annual_energy
        t0 = 1
        for month in 1:12
            plus_hours = daysinmonth(Date(string(year) * "-" * string(month))) * 24
            if month == 2 && isleapyear(year)
                plus_hours -= 24
            end
            month_total = sum(normalized_profile[t0:t0+plus_hours-1])
            if month_total == 0.0  # avoid division by zero
                monthly_scalers[month] = 0.0
            else
                monthly_scalers[month] = monthly_energies[month] / month_total
            end
            t0 += plus_hours
        end
    end

    scaled_load = Float64[]
    boiler_efficiency = 1.0
    mmbtu_to_kwh = 1.0  # do not convert electric loads
    if type in ["domestic_hot_water", "space_heating"]
        # CRB thermal "loads" are in terms of energy input required (boiler fuel), not the actual energy demand.
        # So we multiply the fuel energy by the boiler_efficiency to get the actual energy demand.
        boiler_efficiency = EXISTING_BOILER_EFFICIENCY
        mmbtu_to_kwh = MMBTU_TO_KWH  # do convert thermal loads
    end
    datetime = DateTime(year, 1, 1, 1)
    for ld in normalized_profile
        month = Month(datetime).value
        push!(scaled_load, ld * annual_energy * monthly_scalers[month] * boiler_efficiency * mmbtu_to_kwh)
        datetime += Dates.Hour(1)
    end

    return scaled_load
end

"""
    blend_and_scale_doe_profiles(
        constructor,
        latitude::Float64,
        longitude::Float64,
        year::Int,
        blended_doe_reference_names::Array{String, 1},
        blended_doe_reference_percents::Array{<:Real,1},
        city::String = "",
        annual_energy::Union{Real, Nothing} = nothing,
        monthly_energies::Array{<:Real,1} = Real[],
    )

Given `blended_doe_reference_names` and `blended_doe_reference_percents` use the `constructor` function to load in DoE
    CRB profiles and create a single profile, where `constructor` is one of:
    - BuiltInElectricLoad
    - BuiltInDomesticHotWaterLoad
    - BuiltInSpaceHeatingLoad
"""
function blend_and_scale_doe_profiles(
    constructor,
    latitude::Float64,
    longitude::Float64,
    year::Int,
    blended_doe_reference_names::Array{String, 1},
    blended_doe_reference_percents::Array{<:Real,1},
    city::String = "",
    annual_energy::Union{Real, Nothing} = nothing,
    monthly_energies::Array{<:Real,1} = Real[],
    )

    @assert sum(blended_doe_reference_percents) ≈ 1 "The sum of the blended_doe_reference_percents must equal 1"
    if year != 2017
        @warn "Changing ElectricLoad.year to 2017 because DOE reference profiles start on a Sunday."
    end
    year = 2017
    if isempty(city)
        city = find_ashrae_zone_city(latitude, longitude)  # avoid redundant look-ups
    end
    profiles = Array[]  # collect the built in profiles
    for name in blended_doe_reference_names
        push!(profiles, constructor(city, name, latitude, longitude, year, annual_energy, monthly_energies))
    end
    if isnothing(annual_energy) # then annual_energy should be the sum of all the profiles' annual kwhs
        # we have to rescale the built in profiles to the total_kwh by normalizing them with their
        # own annual kwh and multiplying by the total kwh
        annual_kwhs = [sum(profile) for profile in profiles]
        total_kwh = sum(annual_kwhs)
        monthly_scaler = 1
        if length(monthly_energies) == 12
            monthly_scaler = length(blended_doe_reference_names)
        end
        for idx in 1:length(profiles)
            profiles[idx] .*= total_kwh / annual_kwhs[idx] / monthly_scaler
        end
    end
    for idx in 1:length(profiles)  # scale the profiles
        profiles[idx] .*= blended_doe_reference_percents[idx]
    end
    sum(profiles)
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
function barzanRegressTreeInvoke(model, features)
	#% Invokes the regression tree for the given set of features
	features=features''
	predictions = zeros(size(features,1),1)
	for i=1:size(features,1)
		predictions[i]= apply_tree(model,vec(features[i,:]))
	end
	return predictions
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
module Day3

export get_inputs, get_solution1, get_solution2

## Input getting
is_tree(char) = char=='#'
is_tree(str::String) = [is_tree(char) for char in str]

parse_input(f_name) = readlines(joinpath(@__DIR__, f_name)) .|> is_tree

function get_inputs()
    test_input1 = test_input2 = parse_input("test_input1.txt")

    test_output1 = 7
    test_output2 = 336

    data = parse_input("input.txt")

    return (; test_input1, test_input2, test_output1, test_output2, data)
end

## Solution functions
# Part 1
function get_solution1(data; slope=1//3)
    pos_right = 1
    pos_down = 1
    tree_count = 0
    period = length(data[1])

    while pos_down ≤ length(data)
        tree_count += data[pos_down][pos_right]
        pos_right = mod1(pos_right+slope.den, period)
        pos_down += slope.num
    end

    return tree_count
end

# Part 2
get_solution2(data) = prod(get_solution1(data; slope) for slope in [1//1, 1//3, 1//5, 1//7, 2//1])

end


text:
include("preamble.jl")

T = [0.0 0.0; 1.0 -1.0]
c = [1.0; 0.0]
model = BoundedFluidQueue(T,c,10.0)

hx = 1.0
nodes = 0.0:hx:model.b
t = 4.0
h = 0.005 # I think this should be less than 0.0087 to obey CFL-like stuff
x_lims = (-0.1,10.1)
n_evals = 10_001
norm_pdf(x) = Distributions.pdf(Distributions.Normal(2.5,0.5),x)
norm_cdf(x) = Distributions.cdf(Distributions.Normal(2.5,0.5),x)
orders = 1:2:21
CDF0s = [
    (z,i)->Float64(z>=0.0)*(i==1)*(z>=0.0);
    (z,i)->Float64(z>0.5)*(i==1)*(z>=0.0);
    (z,i)->((z<0.5)*2*z+(z>=0.5))*(i==1)*(z>=0.0);
    # (z,i)->(z>0.5)*(z^2/2-z/2+1/8)*8;
    (z,i)->(z*(z<1.0) + Float64(z>=1.0))*(i==1)*(z>=0.0);
    (z,i)->(z^3/10^3)*(i==1)*(z>=0.0)*(z<6.0) + (z>=6.0)*(i==1);
    # (z,i)->-2*z^3+3*z^2;
    (z,i)->(norm_cdf(z)-norm_cdf(0))/(norm_cdf(10)-norm_cdf(0))*(i==1)*(z>=0.0)#((norm_cdf(z)-norm_cdf(z-2.5))/(1-norm_cdf(z-2.5)))*(i==1)*(z>=0.0)*(z<6.0) + (z>=6.0)*(i==1);
    (z,i)->(2*(z/2 - sin(2*pi*z)/(4*pi))*(z<1.0) + Float64(z>=1.0))*(i==1)*(z>=0.0);
    (z,i)->((2*(z/2 - sin(2*pi*z)/(4*pi)))/10.0)*(i==1)*(z>=0.0)*(z<6.0) + (i==1)*(z>=0.0)*(z>=6.0);
    (z,i)->((cdf(ConcentratedMatrixExponential(5,mean=2.5))(z))./(cdf(ConcentratedMatrixExponential(5,mean=2.5))(10.0)))*(i==1)*(z>=0.0)*(z<6.0) + (z>=6.0)*(i==1)
    ]
PDF0s = [
    (z,i)->0.0;
    (z,i)->NaN;
    (z,i)->2*Float64(z<0.5)*(i==1)*(z>=0.0)*(z<6.0);
    # (z,i)->(z>0.5)*8*(z-0.5);
    (z,i)->1.0*(z<1.0)*(i==1)*(z>=0.0)*(z<6.0);
    (z,i)->z^2*3/10^3*(i==1)*(z>=0.0)*(z<6.0);
    # (z,i)->-6*z^2+6*z;
    (z,i)->(norm_pdf(z)/(norm_cdf(10)-norm_cdf(0)))*(i==1)*(z>=0.0)*(z<6.0);
    (z,i)->((cos(2*π*(z+0.5))+1)*(z<1.0))*(i==1)*(z>=0.0)*(z<6.0);
    (z,i)->(cos(2*π*(z+0.5))+1)/10*(i==1)*(z>=0.0);
    (z,i)->pdf(ConcentratedMatrixExponential(5,mean=2.5))(z)*(i==1)*(z>=0.0)./cdf(ConcentratedMatrixExponential(5,mean=2.5))(10.0)*(z<6.0)
    ]
CDFs = [ (z,j)->exp(T[2,2]*t)*(exp(-T[2,2]*z/c[1])-1)*(z<=c[1]*t) + (z>c[1]*t)*(1-exp(T[2,2]*t));
    [(z,j)->CDF0s[i](z-t,j) for i in 2:length(CDF0s)]]
CDFs[6] = (z,i)->(norm_cdf(z-t)-norm_cdf(0))/(norm_cdf(10)-norm_cdf(0))*(i==1)*(z>=t)*(z<10) + 1(z>=10)

PDFs = [(z,j)->exp(T[2,2]*t)*exp(-T[2,2]*z/c[1])*-T[2,2]/c[1]*(z<=c[1]*t);
    [(z,j)->PDF0s[i](z-t,j) for i in 2:length(PDF0s)]]
PDFs[6] = (z,i)->(norm_pdf(z-t)/(norm_cdf(10)-norm_cdf(0)))*(i==1)*(z>=t);

init_distns = [x->left_point_mass(2,x);
         x->interior_point_mass(0.5,1,x);
    [ x->SFMDistribution(PDF0s[i],x,TrapezoidRule;fun_evals=n_evals) for i in 3:length(PDFs)]]
    # x->SFMDistribution(PDFs[2],x,TrapezoidRule;fun_evals=2001);
    # x->SFMDistribution(PDFs[3],x,TrapezoidRule;fun_evals=2001);
    # [ x->SFMDistribution(PDFs[i],x) for i in 4:length(PDFs)]]
# x->SFMDistribution((z,i)->1.0,x),
# x->SFMDistribution((z,i)->2.0*z,x),
# x->SFMDistribution((z,i)->-6*z^2+6*z,x),
# x->SFMDistribution((z,i)->3*exp(-3*z)/(1-exp(-3)),x),
# x->SFMDistribution((z,i)->cos(2*π*(z+0.5))+1,x),
# x->SFMDistribution((z,i)->cos(4*π*(z+0.5))+1,x),
# x->SFMDistribution((z,i)->pdf(ConcentratedMatrixExponential(5,mean=0.5))(z),x),
# ]
pdf_ylims = [(-0.2,1.2); (-0.2,3.2); (-0.2,2.2); (-0.2,1.2); (-0.1,0.25); (-0.2,1.2); (-0.2,2.2); (-0.05,0.25); (-0.1,0.9)]

m_loop_vec = (DGMesh,DGMesh,DGMesh,FRAPMesh)

# coeff_matrix = Array{Any,3}(undef,length(CDFs),length(orders),length(m_loop_vec))
jldopen("./wave/coeffs_matrix.jld2","r") do f
    global coeff_matrix=f["coeff_matrix"]
end
for (func_count,init_dist_fun) in enumerate(init_distns)
    if func_count==6
        pth = mkpath((@__FILE__)[1:end-3]*"/func_count_"*string(func_count)*"/mesh_comp")
        mkpath(pth*"/data2")
        mkpath(pth*"/figs2")

        L1_cdf_errors = DataFrame(DG = Float64[], Order_1 = Float64[], DG_limiter = Float64[], QBDRAP = Float64[])
        ks_errors = DataFrame(DG = Float64[], Order_1 = Float64[], DG_limiter = Float64[], QBDRAP = Float64[])
        L1_pdf_errors = DataFrame(DG = Float64[], Order_1 = Float64[], DG_limiter = Float64[], QBDRAP = Float64[])
        L2_pdf_errors = DataFrame(DG = Float64[], Order_1 = Float64[], DG_limiter = Float64[], QBDRAP = Float64[])
        c_plt = 0
        for plt_type in (cdf,pdf)
            c_plt += 1
            p = plot(layout=(length(m_loop_vec),length(orders)),legend=false)
            linetypes = [:solid,:solid,:solid,:solid]
            c_order = 0
            for order in orders
                c_order += 1
                c_mesh = 0

                L1_cdf_error_row = zeros(1,length(m_loop_vec))
                ks_error_row = zeros(1,length(m_loop_vec))
                L1_pdf_error_row = zeros(1,length(m_loop_vec))
                L2_pdf_error_row = zeros(1,length(m_loop_vec))
                for mtype in m_loop_vec
                    c_mesh += 1
                    tmp_nodes = c_mesh==2 ? (nodes[1]:(hx/order):nodes[end]) : nodes
                    tmp_order = c_mesh==2 ? 1 : order
                    tmp_nodes = c_mesh==3 ? (nodes[1]:(hx/((order+1)/2)):nodes[end]) : tmp_nodes ##
                    tmp_order = c_mesh==3 ? 2 : tmp_order ##
                    mesh = mtype(tmp_nodes,tmp_order)

                    dq = DiscretisedFluidQueue(model,mesh)
                    B = build_full_generator(dq)

                    d0 = init_dist_fun(dq)

                    # rerun only those that need it
                    if (!isassigned(coeff_matrix,func_count,c_order,c_mesh))||(func_count==6)#||(c_mesh==3)#||(c_mesh==2) ##
                        # throw(DomainError("huh?"))
                        dt = ((mtype==DGMesh) && (c_mesh<3)) ? (limit_str=" NoLimit ";integrate_time(d0,B,t,StableRK4(h); limiter=NoLimiter)) : (limit_str=" ";integrate_time(d0,B,t,StableRK4(h)))
                        coeff_matrix[func_count,c_order,c_mesh] = dt.coeffs
                    else
                        dt = SFMDistribution(coeff_matrix[func_count,c_order,c_mesh],dq)
                    end
                    # global dt ##
                    # @show order ##
                    # @show tmp_order ##
                    # @show c_mesh ##
                    rec = (mtype==FRAPMesh) ? plt_type(dt, eval(Symbol(:normalised_closing_operator_,plt_type))) : plt_type(dt)

                    ## errors
                    (plt_type==cdf) && (L1_cdf_error_row[c_mesh] = log10(DiscretisedFluidQueues.Lp(x->rec(x,1),x->CDFs[func_count](x,1),range(-eps(),10+eps(),length=n_evals))))
                    (plt_type==cdf) && (ks_error_row[c_mesh] = log10(DiscretisedFluidQueues.kolmogorov_smirnov(x->rec(x,1),x->CDFs[func_count](x,1),range(-eps(),10+eps(),length=n_evals))))
                    (plt_type==pdf) && (L1_pdf_error_row[c_mesh] = log10(DiscretisedFluidQueues.Lp(x->rec(x,1),x->PDFs[func_count](x,1),range(-eps(),10+eps(),length=n_evals))))
                    (plt_type==pdf) && (L2_pdf_error_row[c_mesh] = log10(DiscretisedFluidQueues.Lp(x->rec(x,1),x->PDFs[func_count](x,1),range(-eps(),10+eps(),length=n_evals),2)))

                    ## plotting...
                    y_lim_vals = (plt_type==pdf) ? pdf_ylims[func_count] : (-0.2,1.2)
                    # y_tick_vals = ((plt_type==pdf) ? (y_lim_vals[1]:3:y_lim_vals[2]-1) : (0:0.5:1))
                    std_plot(args...; kwargs...) = plot!(
                        p.layout.grid[c_mesh,c_order],
                        x->rec(x,1),
                        nodes[1]+eps(), nodes[end]+eps();
                        ylim=y_lim_vals,
                        linestyle=linetypes[c_mesh],
                        xticks=false,
                        # yticks=false,
                        grid=false,
                        legend=false,
                        kwargs...,
                    )
                    if c_order == 1
                        std_plot(;label=false,ylabel=string(mtype)[1:findfirst("M",string(mtype))[1]-1],legend=false)
                        # yticks!(p.layout.grid[c_mesh,c_order])
                    elseif c_order == length(orders)
                        std_plot(;legend=false, yticks=false)#:outertopright)
                    else
                        std_plot(;label=false,legend=false, yticks=false)
                    end
                    if c_mesh == 1
                        plot!(p.layout.grid[c_mesh,c_order]; title=string(order),legend=false)
                    elseif c_mesh == length(m_loop_vec)
                        xticks!(p.layout.grid[c_mesh,c_order], 0:3:11)
                        plot!(p.layout.grid[c_mesh,c_order]; xlabel="x",legend=false)
                    end
                    std_plot2(args... ; kwargs...) = if plt_type==cdf
                        plot!(
                            p.layout.grid[c_mesh,c_order],
                            x->CDFs[func_count](x,1),
                            nodes[1], nodes[end],
                            args...;
                            linestyle=:dot,
                            legend=false,
                            kwargs...,
                        )
                    elseif plt_type==pdf
                        plot!(
                            p.layout.grid[c_mesh,c_order],
                            x->PDFs[func_count](x,1),
                            nodes[1], nodes[end],
                            args...;
                            subplot=c_order,
                            linestyle=:dot,
                            legend=false,
                            kwargs...
                        )
                    end
                    if c_order==length(orders)
                        std_plot2(;legend=false)#:outertopright)
                    else
                        std_plot2(;label=false)
                    end
                end
                (plt_type==cdf) && push!(L1_cdf_errors,L1_cdf_error_row)
                (plt_type==cdf) && push!(ks_errors,ks_error_row)
                (plt_type==pdf) && push!(L1_pdf_errors,L1_pdf_error_row)
                (plt_type==pdf) && push!(L2_pdf_errors,L2_pdf_error_row)
            end
            # display(p)
            # display(pth)
            savefig(p,pth*"/figs2/mesh_"*string(plt_type)*"_func_count_"*string(func_count)*".svg")
        end
        # #display(L1_cdf_errors)
        file = pth*"/data2/meshs_l1_cdf_"*"func_count_"*string(func_count)
        CSV.write(file*".csv",L1_cdf_errors)
        q = plot()
        for reconstruction in names(L1_cdf_errors)
            plot!(orders,L1_cdf_errors[:,reconstruction],label=reconstruction)
        end
        plot!(xlabel="Order", ylabel="log₁₀ Error",
            title="L¹ error between the true CDF and approximations",
            legend=:outertopright)
        #display(q)
        # #display(pth)
        file = pth*"/figs2/meshs_l1_cdf_"*"func_count_"*string(func_count)
        savefig(q,file*".svg")

        # #display(ks_errors)
        file = pth*"/data2/meshs_ks_"*"func_count_"*string(func_count)
        CSV.write(file*".csv",ks_errors)
        q = plot()
        for reconstruction in names(ks_errors)
            plot!(orders,ks_errors[:,reconstruction],label=reconstruction)
        end
        plot!(xlabel="Order", ylabel="log₁₀ Error",
            title="KS error between the true CDF and approximations",
            legend=:outertopright)
        #display(q)
        file = pth*"/figs2/meshs_ks_"*"func_count_"*string(func_count)
        savefig(q,file*".svg")

        # #display(L1_pdf_errors)
        file = pth*"/data2/meshs_l1_pdf_"*"func_count_"*string(func_count)
        CSV.write(file*".csv",L1_pdf_errors)
        q = plot()
        for reconstruction in names(L1_pdf_errors)
            plot!(orders,L1_pdf_errors[:,reconstruction],label=reconstruction)
        end
        plot!(xlabel="Order", ylabel="log₁₀ Error",
            title="L¹ error between the true PDF and approximations",
            legend=:outertopright)
        #display(q)
        file = pth*"/figs2/meshs_l1_pdf_"*"func_count_"*string(func_count)
        savefig(q,file*".svg")

        # #display(L2_pdf_errors)
        file = pth*"/data2/meshs_l2_pdf_"*"func_count_"*string(func_count)
        CSV.write(file*".csv",L2_pdf_errors)
        q = plot()
        for reconstruction in names(L2_pdf_errors)
            plot!(orders,L2_pdf_errors[:,reconstruction],label=reconstruction)
        end
        plot!(xlabel="Order", ylabel="log₁₀ Error",
            title="L² error between the true PDF and approximations",
            legend=:outertopright)
        #display(q)
        file = pth*"/figs2/meshs_l2_pdf_"*"func_count_"*string(func_count)
        savefig(q,file*".svg")
    end
end
file = (@__FILE__)[1:end-3]*"/coeffs_matrix.jld2"
# jldsave(file;coeff_matrix=coeff_matrix)
## TO READ
# tmp = copy(coeff_matrix)
# jldopen("./wave/coeffs_matrix.jld2","a") do f
#     coeff_matrix=f["coeff_matrix"]
#     coeff_matrix[6,:,:] = tmp[6,:,:]
# end


text:
module LLVM

const platform_name = "llvm_version"
const augment = """
    function augment_llvm!(platform::Platform)
        haskey(platform, "llvm_version") && return p

        llvm_version = Base.libllvm_version
        # does our LLVM build use assertions?
        llvm_assertions = try
            cglobal((:_ZN4llvm24DisableABIBreakingChecksE, Base.libllvm_path()), Cvoid)
            false
        catch
            true
        end
        platform["llvm_version"] = if llvm_assertions
            "\$(llvm_version.major).asserts"
        else
            "\$(llvm_version.major)"
        end
        return platform
    end
"""

function platform(llvm_version, llvm_assertions)
    if llvm_assertions
        return "$(llvm_version.major).asserts"
    else
        return "$(llvm_version.major)"
    end
end

end


text:
module Scenes

export HitRecord, Sphere, Scene, TriangleMesh, ray_intersect, create_triangles, Portal
#export has_uvs, has_normals, get_vertex, get_uv, get_normal

using LinearAlgebra

#push!(LOAD_PATH, pwd())
using ..GfxBase
using ..OBJMeshes
using ..Materials

#####################################
###### Generic Scene Data Type ######
#####################################
struct Scene
    id::Int64
    background::RGB{Float32}
    objects::Array{Any,1}
    lights::Array{Any,1}
end
""" Structure to store data about an intersection
of a ray with an object (a "hit")."""
mutable struct HitRecord
    t::Float64
    intersection::Vec3
    normal::Vec3
    uv::Union{Vec2,Nothing}
    object
end

# Abstract ray-object intersection function:
# Needs to be implemented for each type of object to be rendered
""" Intersect a ray with an object.
Returns a HitRecord with info about the intersetion, or nothing if
the ray doesn't intersect with the object. """
function ray_intersect(ray::Ray, object) end

##################
##### Sphere #####
##################

# Data type:
struct Sphere
    center::Vec3
    radius::Float64
    material::Material
end

mutable struct Portal
    mesh::Array{Any, 1}
    center::Vec3
    buddy::Union{Portal, Nothing}
end

""" Ray-sphere intersection. """
function ray_intersect(ray::Ray, object::Sphere)
    # define constants to use later
    cent = object.center
    d = ray.direction
    p = ray.origin
    r = object.radius

    # find the value of t
    a = dot(d,d)
    b = dot(p - cent,d)
    c = dot(p-cent, p-cent)-r^2
    discrim = b*b-a*c
    dx = p[1] - cent[1]
    dy = p[2] - cent[2]
    dz = p[3] - cent[3]

    # not solns if discriminant is < 0
    if discrim < 0
        return nothing
    else
        num = -b-sqrt(discrim)
        if num > 0
            t = num/a
        else
            num = -b+sqrt(discrim)
            if num > 0
                t = num/a
            else
                return nothing
            end
        end
    end
    # intersection points
    x = p[1] + t*d[1]
    y = p[2] + t*d[2]
    z = p[3] + t*d[3]

    # to use for the uv coords, intersection - center
    dist = normalize(Vec3(x,y,z)-cent)

    normal = normalize((Vec3(x,y,z)- cent)/r)
    if isnothing(object.material.texture)
        return HitRecord(t,Vec3(x,y,z), normal, nothing, object)
    else
        theta = atan(dist[1],dist[3])
        phi = asin(-dist[2]) # negative as image was upside down
        # shift by .5 for middle of the texel
        u = .5 + (theta/ (2*π))
        v = .5 - (phi/π)
        uv = Vec2(u,v)
        return HitRecord(t,Vec3(x,y,z), normal, uv, object)
    end
end

###########################
###### Triangle Mesh ######
###########################

""" Data type: stores the OBJTriangle, a reference to its Mesh
object, and the material it should be rendered with. """
struct Triangle
    geometry::OBJTriangle
    mesh::OBJMesh
    material
end

""" Return an Array of all Triangles belonging to the given mesh, assigning
each one the given material. """
function create_triangles(mesh::OBJMesh, material)
    [Triangle(f, mesh, material) for f in mesh.triangles]
end

""" Some helper functions that make for easier access to triangle data: """
function get_vertex(tri::Triangle, i)
    tri.mesh.positions[tri.geometry.positions[i]]
end

function has_uvs(tri::Triangle)
    length(tri.geometry.uvs) == 3
end

function get_uv(tri::Triangle, i)
    tri.mesh.uvs[tri.geometry.uvs[i]]
end

function has_normals(tri::Triangle)
    length(tri.geometry.normals) == 3
end

function get_normal(tri::Triangle, i)
    tri.mesh.normals[tri.geometry.normals[i]]
end

function ray_intersect(ray::Ray, object::Triangle)

    # vertices
    x = get_vertex(object, 1)
    y = get_vertex(object, 2)
    z = get_vertex(object, 3)

    dir = ray.direction
    org = ray.origin

    # letters from the text
    a = x[1]-y[1]
    b = x[2]-y[2]
    c = x[3]-y[3]
    g = dir[1]
    h = dir[2]
    i = dir[3]
    j = x[1] - org[1]
    k = x[2] - org[2]
    l = x[3] - org[3]
    d = x[1]-z[1]
    e = x[2]-z[2]
    f = x[3]-z[3]

    # the main matrix
    m = (a * ((e * i) - (h * f))) + (b * ((g * f) - (d * i))) + (c * ((d * h) - (e * g)))
    if m == 0 # no solns
        return nothing
    end
    # cramers rule
    t = -(((f * (a * k - j * b)) + (e * (j * c - a * l)) + (d * (b * l - k * c))) / m)
    γ = ((i * (a * k - j * b)) + (h * (j * c - a * l)) + (g * (b * l - k * c))) / m
    if (γ > 1) || (γ < 0)
        return nothing
    end
    β = (j * (e * i - h * f) + k * (g * f - d * i) + l * (d * h - e * g)) / m
    if (β < 0) || (β > (1 -γ))
        return nothing
    end
    # intersection point
    dir *= t
    inter = org + dir
    if !has_normals(object)
        norm = normalize(cross(y-x, z-x))
    else
        # barycentric interpolation for normals
        n1 = get_normal(object, 1)
        n2 = get_normal(object, 2)
        n3 = get_normal(object, 3)
        x = n1[1]+(β*(n2[1]-n1[1]))+(γ*(n3[1]-n1[1]))
        y = n1[2]+(β*(n2[2]-n1[2]))+(γ*(n3[2]-n1[2]))
        z = n1[3]+(β*(n2[3]-n1[3]))+(γ*(n3[3]-n1[3]))
        norm = normalize(Vec3(x,y,z))
    end
    if isnothing(object.material.texture)
        return HitRecord(t, inter, norm, nothing, object)
    else
        if has_uvs(object)
            uv1 = get_uv(object, 1)
            uv2 = get_uv(object, 2)
            uv3 = get_uv(object, 3)

            # barycentric interpolation
            u = uv1[1]+(β*(uv2[1]-uv1[1]))+(γ*(uv3[1]-uv1[1]))
            v = uv1[2]+(β*(uv2[2]-uv1[2]))+(γ*(uv3[2]-uv1[2]))
            uvs = Vec2(u,v)
            return HitRecord(t, inter, norm, uvs, object)
        end
    end
end
end # module Scenes