text:
using NiLangCore
using Test
using Base.Threads

@testset "to_standard_format" begin
    for (OP, FUNC) in [(:+=, PlusEq), (:-=, MinusEq), (:*=, MulEq), (:/=, DivEq), (:⊻=, XorEq)]
        @test NiLangCore.to_standard_format(Expr(OP, :x, :y)) == :($FUNC(identity)(x, y))
        @test NiLangCore.to_standard_format(Expr(OP, :x, :(sin(y; z=3)))) == :($FUNC(sin)(x, y; z=3))

        OPD = Symbol(:., OP)
        @test NiLangCore.to_standard_format(Expr(OPD, :x, :y)) == :($FUNC(identity).(x, y))
        @test NiLangCore.to_standard_format(Expr(OPD, :x, :(sin.(y)))) == :($FUNC(sin).(x, y))
        @test NiLangCore.to_standard_format(Expr(OPD, :x, :(y .* z))) == :($FUNC(*).(x, y, z))
    end
    @test NiLangCore.to_standard_format(Expr(:⊻=, :x, :(y && z))) == :($XorEq($(NiLangCore.logical_and))(x, y, z))
    @test NiLangCore.to_standard_format(Expr(:⊻=, :x, :(y || z))) == :($XorEq($(NiLangCore.logical_or))(x, y, z))
end

@testset "i" begin
    @i function test1(a::T, b, out) where T<:Number
        add(a, b)
        out += a * b
    end

    @i function tt(a, b)
        out ← 0.0
        test1(a, b, out)
        (~test1)(a, b, out)
        a += b
        out → 0.0
    end

    # compute (a+b)*b -> out
    x = 3.0
    y = 4.0
    out = 0.0
    @test isreversible(test1, Tuple{Number, Any, Any})
    @test check_inv(test1, (x, y, out))
    @test check_inv(tt, (x, y))
    @test check_inv(tt, (x, y))
end

@testset "if statement 1" begin
    # compute (a+b)*b -> out
    @i function test1(a, b, out)
        add(a, b)
        if (a > 8, a > 8)
            out += a*b
        else
        end
    end

    x = 3
    y = 4
    out = 0
    @instr test1(x, y, out)
    @test out==0
    @test x==7
    @instr (~test1)(x, y, out)
    @test out==0
    @test x==3
end

@testset "if statement error" begin
    x = 3
    y = 4
    out = 0

    # compute (a+b)*b -> out
    @i function test1(a, b, out)
        add(a, b)
        if (out < 4,  out < 4)
            out += a*b
        else
        end
    end

    @test_throws InvertibilityError test1(x, y, out)
end

@testset "if statement 3" begin
    x = 3
    y = 4
    out = 0
    @i @inline function test1(a, b, out)
        add(a, b)
        if (a > 2, a > 2)
            out += a*b
        else
        end
    end

    x = 3
    y = 4
    out = 0
    @instr test1(x, y, out)
    @test out==28
    @instr (~test1)(x, y, out)
    @test out==0
end

@testset "if statement 4" begin
    @i function test1(a, b, out)
        add(a, b)
        if a > 8.0
            out += a*b
        end
    end
    @test test1(1.0, 8.0, 0.0)[3] == 72.0

    @i function test2(a, b)
        add(a, b)
        if a > 8.0
            a -= b^2
        end
    end
    @test_throws InvertibilityError test2(1.0, 8.0)

    @test_throws LoadError macroexpand(Main, :(@i function test3(a, b)
        add(a, b)
        if a > 8.0
            a -= b*b
        end
    end))
end

@testset "for" begin
    @i function looper(x, y, k)
        for i=1:1:k
            x += y
        end
    end
    x = 0.0
    y = 1.0
    k = 3
    @instr looper(x, y, k)
    @test x == 3
    @instr (~looper)(x, y, k)
    @test x == 0.0

    shiba = 18
    @i function looper2(x, y, k)
        for i=1:1:k
            k += shiba
            x += y
        end
    end
    @test_throws InvertibilityError looper2(x, y, k)
end

@testset "while" begin
    @i function looper(x, y)
        while (x<100, x>0)
            x += y
        end
    end
    x = 0.0
    y = 9
    @instr looper(x, y)
    @test x == 108
    @instr (~looper)(x, y)
    @test x == 0.0

    @i function looper2(x, y)
        while (x<100, x>-10)
            x += y
        end
    end
    @test_throws InvertibilityError looper2(x, y)

   @test_throws LoadError macroexpand(@__MODULE__, :(@i function looper3(x, y)
        while (x<100, x>0)
            z ← 0
            x += y
            z += 1
        end
    end))
end

@testset "ancilla" begin
    one, ten = 1, 10
    @i function looper(x, y)
        z ← 0
        x += y
        z += one
        z -= one
        z → 0
    end
    x = 0.0
    y = 9
    @instr looper(x, y)
    @test x[] == 9
    @instr (~looper)(x, y)
    @test x[] == 0.0

    @i function looper2(x, y)
        z ← 0
        x += y
        z += one
        z -= ten
        z → 0
    end
    x = 0.0
    y = 9
    @test_throws InvertibilityError looper2(x, y)
end

@testset "broadcast" begin
    # compute (a+b)*b -> out
    @i function test1(a, b)
        a .+= b
    end
    x = [3, 1.0]
    y = [4, 2.0]
    @instr test1(x, y)
    @test x == [7, 3.0]
    @instr (~test1)(x, y)
    @test x == [3, 1.0]

    @i function test2(a, b, out)
        a .+= identity.(b)
        out .+= (a .* b)
    end

    x = Array([3, 1.0])
    y = [4, 2.0]
    out = Array([0.0, 1.0])
    @instr test2(x, y, out)
    @test out==[28, 7]
    @test check_inv(test2, (x, y, out))
end

@testset "broadcast arr" begin
    @i function f5(x, y, z, a, b)
        x += y + z
        b += a + x
    end
    @i function f4(x, y, z, a)
        x += y + z
        a += y + x
    end
    @i function f3(x, y, z)
        y += x + z
    end
    @i function f2(x, y)
        y += x
    end
    @i function f1(x)
        l ← zero(x)
        l += x
        x -= 2 * l
        l += x
        l → zero(x)
    end
    a = randn(10)
    b = randn(10)
    c = randn(10)
    d = randn(10)
    e = randn(10)
    aa = copy(a)
    @instr f1.(aa)
    @test aa ≈ -a
    aa = copy(a)
    bb = copy(b)
    @instr f2.(aa, bb)
    @test aa ≈ a
    @test bb ≈ b + a

    aa = copy(a)
    bb = copy(b)
    cc = copy(c)
    @instr f3.(aa, bb, cc)
    @test aa ≈ a
    @test bb ≈ b + a + c
    @test cc ≈ c

    aa = copy(a)
    bb = copy(b)
    cc = copy(c)
    dd = copy(d)
    @instr f4.(aa, bb, cc, dd)
    @test aa ≈ a + b + c
    @test bb ≈ b
    @test cc ≈ c
    @test dd ≈ a + 2b + c + d

    aa = copy(a)
    bb = copy(b)
    cc = copy(c)
    dd = copy(d)
    ee = copy(e)
    @instr f5.(aa, bb, cc, dd, ee)
    @test aa ≈ a + b + c
    @test bb ≈ b
    @test cc ≈ c
    @test dd ≈ d
    @test ee ≈ a + b + c + d + e

    x = randn(5)
    @test_throws AssertionError @instr x .+= c
end

@testset "broadcast tuple" begin
    @i function f5(x, y, z, a, b)
        x += y + z
        b += a + x
    end
    @i function f4(x, y, z, a)
        x += y + z
        a += y + x
    end
    @i function f3(x, y, z)
        y += x + z
    end
    @i function f2(x, y)
        y += x
    end
    @i function f1(x)
        l ← zero(x)
        l += x
        x -= 2 * l
        l += x
        l → zero(x)
    end
    a = (1,2)
    b = (3,1)
    c = (6,7)
    d = (1,11)
    e = (4,1)
    aa = a
    @instr f1.(aa)
    @test aa == -1 .* a
    aa = a
    bb = b
    @instr f2.(aa, bb)
    @test aa == a
    @test bb == b .+ a

    aa = a
    bb = b
    cc = c
    @instr f3.(aa, bb, cc)
    @test aa == a
    @test bb == b .+ a .+ c
    @test cc == c

    aa = a
    bb = b
    cc = c
    dd = d
    @instr f4.(aa, bb, cc, dd)
    @test aa == a .+ b .+ c
    @test bb == b
    @test cc == c
    @test dd == a .+ 2 .* b .+ c .+ d

    aa = a
    bb = b
    cc = c
    dd = d
    ee = e
    @instr f5.(aa, bb, cc, dd, ee)
    @test aa == a .+ b .+ c
    @test bb == b
    @test cc == c
    @test dd == d
    @test ee == a .+ b .+ c .+ d .+ e

    x = (2,1,5)
    @test_throws AssertionError @instr x .+= c
end

@testset "broadcast 2" begin
    # compute (a+b)*b -> out
    @i function test1(a, b)
        a += b
    end
    x = [3, 1.0]
    y = [4, 2.0]
    @instr test1.(x, y)
    @test x == [7, 3.0]
    @instr (~test1).(x, y)
    @test x == [3, 1.0]

    @i function test2(a, b, out)
        add(a, b)
        out += (a * b)
    end

    x = [3, 1.0]
    y = [4, 2.0]
    out = [0.0, 1.0]
    @instr test2.(x, y, out)
    @test out==[28, 7]
    @instr (~test2).(x, y, out)
    @test out==[0, 1.0]
    args = (x, y, out)
    @instr test2.(args...)
    @test args[3]==[28, 7]
end

@testset "neg sign" begin
    @i function test(out, x, y)
        out += x * (-y)
    end
    @test check_inv(test, (0.1, 2.0, -2.5); verbose=true)
end

@testset "@ibounds" begin
    @i function test(x, y)
        for i=1:length(x)
            @inbounds x[i] += y[i]
        end
    end
    @test test([1,2], [2,3]) == ([3,5], [2,3])
end

@testset "kwargs" begin
    @i function test(out, x; y)
        out += x * (-y)
    end
    @test check_inv(test, (0.1, 2.0); y=0.5, verbose=true)
end

@testset "routines" begin
    @i function test(out, x)
        @routine begin
            out += x
        end
        ~@routine
    end
    out, x = 0.0, 1.0
    @instr test(out, x)
    @test out == 0.0
end

@testset "inverse a prog" begin
    @i function test(out, x)
        ~(out += x;
        out += x)
        ~for i=1:3
            out += x
        end
    end
    out, x = 0.0, 1.0
    @test check_inv(test, (out, x))
    @instr test(out, x)
    @test out == -5.0
end

@testset "invcheck" begin
    @i function test(out, x)
        anc ← 0
        @invcheckoff for i=1:x[]
            x[] -= 1
        end
        @invcheckoff while (anc<3, anc<3)
            anc += 1
        end
        out += anc
        @invcheckoff anc → 0
    end
    res = test(0, Ref(7))
    @test res[1] == 3
    @test res[2][] == 0
end

@testset "nilang ir" begin
    ex = :(
        @inline function f(x!::T, y) where T
            anc ← zero(T)
            @routine anc += x!
            x! += y * anc
            ~@routine
            anc → zero(T)
        end
    )
    ex2 = :(
    @inline function f(x!::T, y) where T
          anc ← zero(T)
          anc += identity(x!)
          x! += y * anc
          anc -= identity(x!)
          anc → zero(T)
    end)

    ex3 = :(
    @inline function (~f)(x!::T, y) where T
          anc ← zero(T)
          anc += identity(x!)
          x! -= y * anc
          anc -= identity(x!)
          anc → zero(T)
    end)
    @test nilang_ir(@__MODULE__, ex) |> NiLangCore.rmlines == ex2 |> NiLangCore.rmlines
    @test nilang_ir(@__MODULE__, ex; reversed=true) |> NiLangCore.rmlines == ex3 |> NiLangCore.rmlines
end

@testset "protectf" begin
    struct C<:Function end
    # protected
    @i function (a::C)(x)
        @safe @show a
        if (protectf(a) isa Inv, ~)
            add(x, 1.0)
        else
            sub(x, 1.0)
        end
    end
    a = C()
    @test (~a)(a(1.0)) == 1.0
    # not protected
    @i function (a::C)(x)
        @safe @show a
        if (a isa Inv, ~)
            add(x, 1.0)
        else
            sub(x, 1.0)
        end
    end
    @test (~a)(a(1.0)) == -1.0
end

@testset "ifelse statement" begin
    @i function f(x, y)
        if (x > 0, ~)
            y += 1
        elseif (x < 0, ~)
            y += 2
        else
            y += 3
        end
    end
    @test f(1, 0) == (1, 1)
    @test f(-2, 0) == (-2, 2)
    @test f(0, 0) == (0, 3)

    @i function f2(x, y)
        if (x > 0, x < 0)
            y += 1
        elseif (x < 0, x < 0)
            y += 2
        else
            y += 3
        end
    end
    @test_throws InvertibilityError f2(-1, 0)
end

@testset "skip!" begin
    x = 0.4
    @instr (@skip! 3) += x
    @test x == 0.4
    y = 0.3
    @instr x += @const y
    @test x == 0.7
    @test y == 0.3
end

@testset "for x in range" begin
    @i function f(x, y)
        for item in y
            x += item
        end
    end
    @test check_inv(f, (0.0, [1,2,5]))
end

@testset "@simd and @threads" begin
    @i function f(x)
        @threads for i=1:length(x)
            x[i] += 1
        end
    end
    x = [1,2,3]
    @test f(x) == [2,3,4]
    @i function f2(x)
        @simd for i=1:length(x)
            x[i] += 1
        end
    end
    x = [1,2,3]
    @test f2(x) == [2,3,4]
end

@testset "xor over ||" begin
    x = false
    @instr x ⊻= true || false
    @test x
    @instr x ⊻= true && false
    @test x
end

macro zeros(T, x, y)
    esc(:($x ← zero($T); $y ← zero($T)))
end

@testset "macro" begin
    @i function f(x)
        @zeros Float64 a b
        x += a * b
        ~@zeros Float64 a b
    end
    @test f(3.0) == 3.0
end

@testset "allow nothing pass" begin
    @i function f(x)
        nothing
    end
    @test f(2) == 2
end

@testset "ancilla check" begin
    ex1 = :(@i function f(x)
        x ← 0
    end)
    @test_throws LoadError macroexpand(Main, ex1)

    ex2 = :(@i function f(x)
        y ← 0
        y ← 0
    end)
    @test_throws LoadError macroexpand(Main, ex2)

    ex3 = :(@i function f(x)
        y ← 0
        y → 0
    end)
    @test macroexpand(Main, ex3) isa Expr

    ex4 = :(@i function f(x; y=5)
        y ← 0
    end)
    @test_throws LoadError macroexpand(Main, ex4)

    ex5 = :(@i function f(x)
        y → 0
    end)
    @test_throws LoadError macroexpand(Main, ex5)

    ex6 = :(@i function f(x::Int)
        y ← 0
        y → 0
    end)
    @test macroexpand(Main, ex6) isa Expr

    ex7 = :(@i function f(x::Int)
        if x>3
            y ← 0
            y → 0
        elseif x<-3
            y ← 0
            y → 0
        else
            y ← 0
            y → 0
        end
    end)
    @test macroexpand(Main, ex7) isa Expr

    ex8 = :(@i function f(x; y=5)
        z ← 0
        z → 0
    end)
    @test macroexpand(Main, ex8) isa Expr

    ex9 = :(@i function f(x; y)
        z ← 0
        z → 0
    end)
    @test macroexpand(Main, ex9) isa Expr

    ex10 = :(@i function f(x; y)
        begin
            z ← 0
        end
        ~begin
            z ← 0
        end
    end)
    @test macroexpand(Main, ex10) isa Expr
end

@testset "dict access" begin
    d = Dict(3=>4)
    @instr d[3] → 4
    @instr d[4] ← 3
    @test d == Dict(4=>3)
    @test_throws InvertibilityError @instr d[4] → 5
    @test (@instr @invcheckoff d[8] → 5; true)
    @test_throws InvertibilityError @instr d[4] ← 5
    @test (@instr @invcheckoff d[4] ← 5; true)
end

@testset "@routine,~@routine" begin
    @test_throws LoadError macroexpand(Main, :(@i function f(x)
        @routine begin
        end
    end))
    @test_throws LoadError macroexpand(Main, :(@i function f(x)
        ~@routine
    end))
    @test macroexpand(Main, :(@i function f(x)
        @routine begin end
        ~@routine
    end)) !== nothing
end

@testset "@from post while pre" begin
    @i function f()
       x ← 5
       z ← 0
       @from z==0 while x > 0
           x -= 1
           z += 1
       end
       z → 5
       x → 0
    end

    @test f() == ()
    @test (~f)() == ()
end

@testset "argument with function call" begin
    @test_throws LoadError @macroexpand @i function f(x, y)
        x += sin(exp(y))
    end
    @i function f(x, y)
        x += sin(exp(0.4)) + y
    end
end

@testset "allocation multiple vars" begin
    info = NiLangCore.PreInfo()
    @test NiLangCore.precom_ex(NiLangCore, :(x,y ← var), info) == :((x, y) ← var)
    @test NiLangCore.precom_ex(NiLangCore, :(x,y → var), info) == :((x, y) → var)
    @test NiLangCore.precom_ex(NiLangCore, :((x,y) ↔ (a, b)), info) == :((x,y) ↔ (a,b))
    @test (@code_reverse (x,y) ← var) == :((x, y) → var)
    @test (@code_reverse (x,y) → var) == :((x, y) ← var)
    @test (@code_julia (x,y) ← var) == :((x, y) = var)
    @test (@code_julia (x,y) → var) == :(try
        $(NiLangCore.deanc)((x, y), var)
    catch e
        $(:(println("deallocate fail `$($(QuoteNode(:((x, y))))) → $(:var)`")) |> NiLangCore.rmlines)
        throw(e)
    end)

    x = randn(2,4)
    @i function f(y, x)
        m, n ← size(x)
        (l, k) ← size(x)
        y += m*n
        y += l*k
        (l, k) → size(x)
        m, n → size(x)
    end
    twosize = f(0, x)[1]
    @test  twosize == 16
    @test (~f)(twosize, x)[1] == 0

    @i function g(x)
        (m, n) ← size(x)
        (m, n) → (7, 5)
    end

    @test_throws InvertibilityError g(x)
    @test_throws InvertibilityError (~g)(x)
end

@testset "argument without argname" begin
    @i function f(::Complex)
    end
    @test f(1+2im) == 1+2im
end

@testset "tuple input" begin
    @i function f(x::Tuple{<:Tuple, <:Real})
        f(x.:1)
        (x.:1).:1 += x.:2
    end
    @i function f(x::Tuple{<:Real, <:Real})
        x.:1 += x.:2
    end
    @i function g(data)
        f(((data.:1, data.:2), data.:3))
    end
    @test g((1,2,3)) == (6,2,3)
end

@testset "single argument" begin
    @i function f(x)
        neg(x)
    end
    @i function g(x::Vector)
        neg.(x)
    end
    @test f(3) == -3
    @test g([3, 2]) == [-3, -2]
    x = (3,)
    @instr f(x...)
    @test x == (-3,)
    x = ([3, 4],)
    @instr f.(x...)
    @test x == ([-3, -4],)
end

@testset "type constructor" begin
    @i function f(x, y, a, b)
        add(Complex{}(x, y), Complex{}(a, b))
    end
    @test f(1,2, 3, 4) == (4, 6, 3, 4)
    @test_throws LoadError macroexpand(NiLangCore, :(@i function f(x, y, a, b)
        add(Complex(x, y), Complex{}(a, b))
    end))
    @i function g(x::Inv, y::Inv)
        add(x.f, y.f)
    end
    @i function g(x, y)
        g(Inv{}(x), Inv{}(y))
    end
    @test g(2, 3) == (5, 3)
end

@testset "variable_analysis" begin
    # kwargs should not be assigned
    @test_throws LoadError macroexpand(@__MODULE__, :(@i function f1(x; y=4)
        y ← 5
        y → 5
    end))
    # deallocated variables should not be used
    @test_throws LoadError macroexpand(@__MODULE__, :(@i function f1(x; y=4)
        z ← 5
        z → 5
        x += 2 * z
    end))
    # deallocated variables should not be used in local scope
    @test_throws LoadError macroexpand(@__MODULE__, :(@i function f1(x; y=4)
        z ← 5
        z → 5
        for i=1:10
            x += 2 * z
        end
    end))
end

@testset "boolean" begin
    @i function f1(x, y, z)
        x ⊻= true
        y .⊻= z
    end
    @test f1(false, [true, false], [true, false]) == (true, [false, false], [true, false])

    @i function f2(x, y, z)
        z[2] ⊻= true && y[1]
        z[1] ⊻= z[2] || x
    end
    @test f2(false, [true, false], [true, false]) == (false, [true, false], [false, true])
end

@testset "swap ↔" begin
    @i function f1(x, y)
        j::∅ ↔ k::∅   # dummy swap
        a::∅ ↔ x
        a ↔ y
        a ↔ x::∅   # ↔ is symmetric
    end
    @test f1(2, 3) == (3, 2)
    @test check_inv(f1, (2, 3))

    # stack
    @i function f2(x, y)
        x[end+1] ↔ y
        y ← 2
    end
    @test f2([1,2,3], 4) == ([1,2,3,4], 2)
    @test check_inv(f2, ([1,2,3], 3))

    @i function f4(x, y)
        y ↔ x[end+1]
        y ← 2
    end
    @test f4([1,2,3], 4) == ([1,2,3,4], 2)
    @test check_inv(f4, ([1,2,3], 3))

    @i function f3(x, y::TY, s) where TY
        y → _zero(TY)
        x[end] ↔ (y::TY)::∅
        @safe @show x[2], s
        x[2] ↔ s
    end
    @test f3(Float32[1,2,3], 0.0, 4f0) == (Float32[1,4], 3.0, 2f0)
    @test check_inv(f3, (Float32[1,2,3], 0.0, 4f0))
end

@testset "feed tuple and types" begin
    @i function f3(a, d::Complex)
        a.:1 += d.re
        d.re ↔ d.im
    end
    @i function f4(a, b, c, d, e)
        f3((a, b, c), Complex{}(d, e))
    end
    @test f4(1,2,3,4,5) == (5,2,3,5,4)
    @test check_inv(f4, (1,2,3,4,5))
end

@testset "exchange tuple and fields" begin
    @i function f1(x, y, z)
        (x, y) ↔ @fields z
    end
    @test f1(1,2, 3+4im) == (3,4,1+2im)

    @i function f2(re, x)
        r, i ← @fields x
        re += r
        r, i → @fields x
    end
    @test f2(0.0, 3.0+2im) == (3.0, 3.0 + 2.0im)

    @i function f3(x, y, z)
        (@fields z) ↔ (x, y)
    end
    @test f3(1,2, 3+4im) == (3,4,1+2im)

    @test_throws LoadError macroexpand(@__MODULE__, :(@i function f3(x, y, z)
        (x, y) ↔ (z, j)
    end))
    @i function f4(x, y, z, j)
        (x, y) ↔ (z, j)
    end
    @test f4(1,2, 3, 4) == (3,4,1,2)

    @i function swap_fields(obj::Complex)
        (x, y)::∅ ↔ @fields obj
        x += y
        (x, y) ↔ (@fields obj)::∅
    end

    @test swap_fields(1+2im) == (3+2im)
end


text:
function integrant_volume_to_surface_area(f)
    f_x(u) = ForwardDiff.derivative(x -> f(x, u[2]), u[1])
    f_z(u) = ForwardDiff.derivative(z -> f(u[1], z), u[2])
    u -> sqrt(f_x(u)^2 + f_z(u)^2 + 1)
end

function integrate_xonly(f, xmin, xmax)
    res = hquadrature(x -> f(SA[x, bd.zmax-bd.hb]), xmin, xmax)
    res[1]
end
function integrate_limit(f, f_zlimit; xmin, xmax, zmin = bd.zmax - bd.hb, zmax = bd.zmax)
    function hcfun(u)
        _x, _z = u
        _z > f_zlimit(_x) ? 0.0 : f(u)
    end
    a = SA[xmin, zmin]
    b = SA[xmax, zmax]
    res = @test_nowarn hcubature(hcfun, a, b, maxevals = 1_000_000)
    # println(res[2])
    res[1]
end
# Start with the middle part
integrate_middle(f = bd.fi2m; kwargs...) =
    integrate_limit(f, (x) -> bd.zmax; xmin = bd.xmin0, xmax = bd.xmax0, kwargs...)
# dziob - psi2, fi2f
integrate_front(f = bd.fi2f; kwargs...) =
    integrate_limit(f, bd.psi2; xmin = bd.xmax0, xmax = bd.xmax, kwargs...)
# rufa - psi1, fi2r
integrate_rear(f = bd.fi2r; kwargs...) =
    integrate_limit(f, bd.psi1; xmin = bd.xmin, xmax = bd.xmin0, kwargs...)

@testset "Use adaptive cubature to approximate boat mass" begin
    # kadlub
    m_m_b = integrate_middle(integrant_volume_to_surface_area(bd.fi2m))
    m_f_b = integrate_front(integrant_volume_to_surface_area(bd.fi2f))
    m_r_b = integrate_rear(integrant_volume_to_surface_area(bd.fi2r))
    snb = 2 * (m_m_b + m_f_b + m_r_b)
    m_hull = snb * bd.rhob * bd.deltab
    # @test snb ≈ bd.snb
    @test 2m_m_b ≈ bd.snbmv
    @test m_f_b ≈ bd.iinbfvv rtol = 1e-4
    @test m_r_b ≈ bd.iinbrvv rtol = 1e-4

    # poklad
    m_m_d = integrate_xonly(u -> bd.fi2m(u...), bd.xmin0, bd.xmax0)
    m_f_d = integrate_xonly(u -> bd.fi2f(u...), bd.xmax0, bd.xmax)
    m_r_d = integrate_xonly(u -> bd.fi2r(u...), bd.xmin, bd.xmin0)
    sbd = 2 * (m_m_d + m_f_d + m_r_d)
    sbde = bd.bbe * (bd.xmaxe - bd.xmine)
    m_deck_full = sbd * bd.rhod * bd.deltad
    m_deck_hole = bd.rhod * bd.deltad * sbde
    @test sbd ≈ bd.sbd
    @test sbde ≈ bd.sbde
    # podłoga
    sbdg = bd.bbg * (bd.xmaxg - bd.xming)
    m_ground = bd.rhog * bd.deltag * sbdg
    @test sbdg ≈ bd.sbdg

    m = m_hull + m_deck_full - m_deck_hole + m_ground
    @test m ≈ bd.m_B rtol = 1e-5
end

@testset "Use adaptive cubature to approximate boat static inertia sxmb" begin
    # kadlub
    sxfun(f) = u -> u[1] * f(u)
    sxfun_area(f) = sxfun(integrant_volume_to_surface_area(f))
    sx_m_b = integrate_middle(sxfun_area(bd.fi2m))
    sx_f_b = integrate_front(sxfun_area(bd.fi2f))
    sx_r_b = integrate_rear(sxfun_area(bd.fi2r))
    x1snb = 2 * (sx_m_b + sx_f_b + sx_r_b)
    sx_hull = x1snb * bd.rhob * bd.deltab
    # @test x1snb ≈ bd.x1snb
    @test 2sx_m_b ≈ bd.x1snbmv
    @test sx_f_b ≈ bd.ix1inbfvv rtol = 1e-3
    @test_broken sx_r_b ≈ bd.ix1inbrvv rtol = 1e-3

    # poklad
    sxfun2(f) = sxfun(u -> f(u[1], u[2]))
    sx_m_d = integrate_xonly(sxfun2(bd.fi2m), bd.xmin0, bd.xmax0)
    sx_f_d = integrate_xonly(sxfun2(bd.fi2f), bd.xmax0, bd.xmax)
    sx_r_d = integrate_xonly(sxfun2(bd.fi2r), bd.xmin, bd.xmin0)
    x1sbd = 2 * (sx_m_d + sx_f_d + sx_r_d)
    x1sbde = 0.5bd.bbe * (bd.xmaxe - bd.xmine)^2
    sx_deck_full = x1sbd * bd.rhod * bd.deltad
    sx_deck_hole = bd.rhod * bd.deltad * x1sbde
    @test x1sbd ≈ bd.x1sbd
    @test_broken x1sbde ≈ bd.x1sbde
    # podłoga
    x1sbdg = 0.5bd.bbg * (bd.xmaxg - bd.xming)^2
    sx_ground = bd.rhog * bd.deltag * x1sbdg
    @test_broken x1sbdg ≈ bd.x1sbdg

    sx = sx_hull + sx_deck_full - sx_deck_hole + sx_ground
    @test_broken sx ≈ bd.sxmb rtol = 1e-5
    # println((sx / bd.m_B, bd.xo))
end

@testset "Use adaptive cubature to approximate boat static inertia szmb" begin
    # kadlub
    szfun(f) = u -> u[2] * f(u)
    szfun_area(f) = szfun(integrant_volume_to_surface_area(f))
    sz_m_b = integrate_middle(szfun_area(bd.fi2m))
    sz_f_b = integrate_front(szfun_area(bd.fi2f))
    sz_r_b = integrate_rear(szfun_area(bd.fi2r))
    z1snb = 2 * (sz_m_b + sz_f_b + sz_r_b)
    sz_hull = z1snb * bd.rhob * bd.deltab
    # @test x1snb ≈ bd.x1snb
    @test 2sz_m_b ≈ bd.z1snbmnv
    @test sz_f_b ≈ bd.iz1inbfvv rtol = 8e-3
    @test sz_r_b ≈ bd.iz1inbrvv rtol = 1e-2

    # poklad
    m_m_d = integrate_xonly(u -> bd.fi2m(u...), bd.xmin0, bd.xmax0)
    m_f_d = integrate_xonly(u -> bd.fi2f(u...), bd.xmax0, bd.xmax)
    m_r_d = integrate_xonly(u -> bd.fi2r(u...), bd.xmin, bd.xmin0)
    sbd = 2 * (m_m_d + m_f_d + m_r_d)
    z1sbd = sbd * (bd.zmax - bd.hb)
    @test z1sbd ≈ bd.z1sbd
    z1sbde = bd.bbe * (bd.xmaxe - bd.xmine) * (bd.zmax - bd.hb)
    @test z1sbde ≈ bd.z1sbde
    sz_deck_full = z1sbd * bd.rhod * bd.deltad
    sz_deck_hole = bd.rhod * bd.deltad * z1sbde
    # podłoga
    z1sbdg = bd.bbg * (bd.xmaxg - bd.xming) * (bd.zmax - bd.hb + bd.zg)
    @test z1sbdg ≈ bd.z1sbdg
    sz_ground = bd.rhog * bd.deltag * z1sbdg

    sz = sz_hull + sz_deck_full - sz_deck_hole + sz_ground
    @test sz ≈ bd.szmb rtol = 5e-3
    # println((sz / bd.m_B, bd.zo))
end

@testset "Use adaptive cubature to approximate boat static inertia j_By" begin
    # kadlub
    x2fun(f) = u -> u[1]^2 * f(u)
    z2fun(f) = u -> u[2]^2 * f(u)
    x2fun_area(f) = x2fun(integrant_volume_to_surface_area(f))
    z2fun_area(f) = z2fun(integrant_volume_to_surface_area(f))
    x2_m_b = integrate_middle(x2fun_area(bd.fi2m))
    x2_f_b = integrate_front(x2fun_area(bd.fi2f))
    x2_r_b = integrate_rear(x2fun_area(bd.fi2r))
    z2_m_b = integrate_middle(z2fun_area(bd.fi2m))
    z2_f_b = integrate_front(z2fun_area(bd.fi2f))
    z2_r_b = integrate_rear(z2fun_area(bd.fi2r))
    x2snb = 2 * (x2_m_b + x2_f_b + x2_r_b)
    z2snb = 2 * (z2_m_b + z2_f_b + z2_r_b)
    J_hull = (x2snb + z2snb) * bd.rhob * bd.deltab
    @test 2x2_m_b ≈ bd.x2snbmv
    @test x2_f_b ≈ bd.ix2inbfvv rtol = 4e-4
    @test_broken x2_r_b ≈ bd.ix2inbrvv rtol = 1e-2

    @test 2z2_m_b ≈ bd.z2snbmnv
    @test z2_f_b ≈ bd.iz2inbfvv rtol = 1e-2
    @test z2_r_b ≈ bd.iz2inbrvv rtol = 1e-2

    # poklad
    m_m_d = integrate_xonly(u -> bd.fi2m(u...), bd.xmin0, bd.xmax0)
    m_f_d = integrate_xonly(u -> bd.fi2f(u...), bd.xmax0, bd.xmax)
    m_r_d = integrate_xonly(u -> bd.fi2r(u...), bd.xmin, bd.xmin0)
    sbd = 2 * (m_m_d + m_f_d + m_r_d)
    z2sbd = sbd * (bd.zmax - bd.hb)^2
    @test z2sbd ≈ bd.z2sbd
    x2_m_d = integrate_xonly(u -> u[1]^2 * bd.fi2m(u...), bd.xmin0, bd.xmax0)
    x2_f_d = integrate_xonly(u -> u[1]^2 * bd.fi2f(u...), bd.xmax0, bd.xmax)
    x2_r_d = integrate_xonly(u -> u[1]^2 * bd.fi2r(u...), bd.xmin, bd.xmin0)
    x2sbd = 2 * (x2_m_d + x2_f_d + x2_r_d)
    @test x2sbd ≈ bd.x2sbd
    z2sbde = bd.bbe * (bd.xmaxe - bd.xmine)^2 * (bd.zmax - bd.hb)
    @test z2sbde ≈ bd.z2sbde
    x2sbde = bd.bbe * (bd.xmaxe^3 - bd.xmine^3) / 3
    @test x2sbde ≈ bd.x2sbde
    J_deck_full = (x2sbd + z2sbd) * bd.rhod * bd.deltad
    J_deck_hole = bd.rhod * bd.deltad * (x2sbde + z2sbde)
    # podłoga
    x2sbdg = bd.bbg * (bd.xmaxg^3 - bd.xming^3) / 3
    z2sbdg = bd.bbg * (bd.xmaxg - bd.xming) * (bd.zmax - bd.hb + bd.zg)^2
    @test x2sbdg ≈ bd.x2sbdg
    @test z2sbdg ≈ bd.z2sbdg
    J_ground = bd.rhog * bd.deltag * (x2sbdg + z2sbdg)

    J_st = bd.m_B * (bd.xo^2 + bd.zo^2)
    J = J_hull + J_deck_full - J_deck_hole + J_ground - J_st
    @test_broken J ≈ bd.J_By rtol = 5e-3
    println((J, bd.J_By))
end

@testset "Check if approximation of vh2o is correct" begin
    tol = eps()
    @test abs(bd.vh2o(0.0, bd.zmax - tol)) < tol
    # vh2o(0.0, 0.0) - should be the full boat volume!
    V_m_b = integrate_middle(u -> bd.fi2m(u...))
    V_f_b = integrate_front(u -> bd.fi2f(u...))
    V_r_b = integrate_rear(u -> bd.fi2r(u...))
    @test bd.vh2o(0.0, 0.0) ≈ 2 * (V_m_b + V_f_b + V_r_b) rtol = 1e-6
    # compute for some bh2o value
    bh2o = 0.111245
    V_m_b = integrate_middle(u -> bd.fi2m(u...), zmin = bh2o)
    V_f_b = integrate_front(u -> bd.fi2f(u...), zmin = bh2o)
    V_r_b = integrate_rear(u -> bd.fi2r(u...), zmin = bh2o)
    @test bd.vh2o(0.0, bh2o) ≈ 2 * (V_m_b + V_f_b + V_r_b) rtol = 1e-6
    # the = 0.003
    # R = SA[cos(the) -sin(the); sin(the) cos(the)]
    # V_m_b = integrate_middle(u -> bd.fi2m((R'*u)...), zmin = bh2o)
    # @test bd.vh2om(the, bh2o) ≈ 2 * (V_m_b) rtol = 1e-6
    # V_f_b = integrate_front(u -> bd.fi2f(u...), zmin = bh2o)
    # V_r_b = integrate_rear(u -> bd.fi2r(u...), zmin = bh2o)
    # @test bd.vh2o(the, bh2o) ≈ 2 * (V_m_b + V_f_b + V_r_b) rtol = 1e-6
end

ext: jl


text:
#   This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2016 Nathan Smith
export MExpr,
       @m_str,
       parse,
       mcall,
       convert,
       error,
       MaximaError,
       MaximaSyntaxError,
       ==,
       getindex

import Base: parse,
             convert,
             error,
             ==,
             getindex,
             *,
             split

struct MaximaError <: Exception
    errstr::String
end

Base.showerror(io::IO, err::MaximaError) = print(io, err.errstr)

struct MaximaSyntaxError <: Exception
    errstr::String
end

Base.showerror(io::IO, err::MaximaSyntaxError) = print(io, err.errstr)

const infix_ops = [:+, :-, :*, :/, :^]

isinfix(args) = args[1] in infix_ops && length(args) > 2
show_expr(io::IO, ex) = print(io, ex)

function show_expr(io::IO, expr::Expr)
    if expr.head == :call
        seperator = isinfix(expr.args) ? " $(expr.args[1]) " : ", "
        !isinfix(expr.args) && show_expr(io, expr.args[1])
        print(io, "(")
        args = expr.args[2:end]
        for (i, arg) in enumerate(args)
            show_expr(io, arg)
            i != lastindex(args) ? print(io, seperator) : print(io, ")")
        end
    elseif expr.head == :(=)
        show_expr(io,expr.args[1])
        print(io,": ")
        show_expr(io,expr.args[2])
    elseif expr.head == :function
        show_expr(io,expr.args[1])
        print(io," := block([], ")
        args = unparse(expr.args[2])
        for (i, arg) in enumerate(args)
            print(io,arg)
            i != endof(args) ? print(io, ", ") : print(io, ")")
        end
    elseif expr.head == :line
        nothing
    else
      error("Nested :$(expr.head) block structure not supported by Maxima.jl")
    end
end

function unparse(expr::Expr)
    str = String[]
    io = IOBuffer();
    if expr.head == :block
        for line ∈ expr.args
            show_expr(io,line)
            push!(str,takebuf_string(io))
        end
        return mtrim(str)
    else
        show_expr(io, expr)
        return push!(str,String(take!(copy(io))))
    end
end

"""

A Maxima expression

## Summary:

type MExpr <: Any

## Fields:

str :: String
"""
struct MExpr
    str::Vector{String}
    MExpr(m::Vector{String}) = new(m)
end

MExpr(m::Array{SubString{String},1}) = MExpr(convert(Array{String, 1}, m))
MExpr(str::String) = MExpr(push!(String[], str))
MExpr(m::Any) = MExpr("$m")

macro m_str(str)
    MExpr(str)
end

*(x::MExpr, y::String) = MExpr(push!(deepcopy(x.str), y))
*(x::String, y::MExpr) = MExpr(unshift!(deepcopy(y.str), x))
*(x::MExpr, y::MExpr) = MExpr(vcat(x.str..., y.str...))

function mtrim(m::Vector{String})
    n = String[]
    for h ∈ 1:length(m)
        !isempty(m[h]) && push!(n, m[h])
    end
    return n
end

function split(m::MExpr)
    n = String[]
    for i in 1:length(m.str)
        p = split(replace(m.str[i], r"\$" => ";"), ';')
        for j in 1:length(p)
            push!(n, p[i])
        end
    end
    return MExpr(n)
end

const m_to_jl = Dict(
    "%e"     => "ℯ",
    "%pi"    =>  "π",
    "%i"     =>  "im",
    "%gamma" => "eulergamma",
    "%phi"   =>  "φ",
    "inf"    =>  "Inf",
    "minf"   =>  "-Inf"
)

const jl_to_m = Dict(
    "ℯ"          => "%e",
    "pi"         => "%pi",
    "π"          => "%pi",
    "γ"          => "%gamma",
    "eulergamma" => "%gamma",
    "golden"     => "%phi",
    "φ"          => "%phi",
    "im"         => "%i",
    "Inf"        => "inf"
)

_subst(a, b, expr) = "subst($a, $b, '($expr))" |> MExpr |> mcall

"""
    MExpr(expr::Expr)

Convert Julia expression to Maxima expression

## Examples
```julia
julia> MExpr(:(sin(x*im) + cos(y*φ)))

                           cos(%phi y) + %i sinh(x)

```
"""
function MExpr(expr::Expr)
    # str = "$expr"
    str = unparse(expr)
    for h in 1:length(str)
        for key in keys(jl_to_m)
            str[h] = _subst(jl_to_m[key], key, str[h])
        end
    end
    return MExpr(str)
end

"""
    parse(mexpr::MExpr)

Parse a Maxima expression into a Julia expression

## Examples
```julia
julia> parse(m\"sin(%i*x)\")
:(im * sinh(x))

```
"""
function parse(m::MExpr)
    pexpr = Any[]
    for subexpr in split(m).str
        for key in keys(m_to_jl)
            subexpr = replace(subexpr, key => m_to_jl[key])
        end
        if occursin( ":=", subexpr)
            sp = split(subexpr, ":=")
            push!(pexpr, Expr(:function, parse(sp[1]), sp[2] |> String |> MExpr |> parse))
        elseif occursin( "block([], ", subexpr)
            rp = replace(subexpr, "block([],", "") |> chop
            sp = split(rp, ",")
            ep = Vector{Any}(length(sp))
            for j in 1:length(sp)
                ep[j] = sp[j] |> String |> MExpr |> parse
            end
            push!(pexpr,Expr(:block,ep...))
        elseif occursin( ":", subexpr)
            sp = split(subexpr, ":")
            push!(pexpr, Expr(:(=), parse(sp[1]), sp[2] |> String |> MExpr |> parse))
        else
            push!(pexpr, Meta.parse(subexpr))
        end
    end
    return length(pexpr) == 1 ? pexpr[1] : Expr(:block, pexpr...)
end

convert(::Type{MExpr}, m::MExpr) = m
convert(::Type{Array{String,1}}, m::MExpr) = m.str
convert(::Type{String}, m::MExpr) = join(m.str, "; ")
convert(::Type{T}, m::MExpr) where T = T <: Number ? eval(parse(m)) : parse(m)

"""
    mcall(m::MExpr)

Evaluate a Maxima expression.

## Examples
```julia
julia> m\"integrate(sin(x), x)\"

                             integrate(sin(x), x)

julia> mcall(ans)

                                   - cos(x)

```
"""
function mcall(m::MExpr)
    write(ms, replace(convert(String, m), r";" => "; print(ascii(3))\$ "))
    output = read(ms)
    if occursin(maxerr, output)
            write(ms.input, "errormsg()\$")
            write(ms.input, "print(ascii(4))\$")
            message = read(ms)
            throw(MaximaError(message))
    elseif occursin(synerr, output)
        throw(MaximaSyntaxError(output))
    else
        sp = split(output, '\x03')
        for k in 1:length(sp)
          sp[k] = replace(sp[k], '\n' => "")
          sp[k] = replace(sp[k], ' ' => "")
        end
        return MExpr(sp)
    end
end

"""
	mcall(expr::T) where T

Evaluate a Julia expression or string using the Maxima interpretor and convert
output back into the input type

## Examples
```julia
julia> mcall(\"integrate(sin(y)^2, y)\")
\"(y-sin(2*y)/2)/2\"

julia> mcall(:(integrate(1/(1+x^2), x)))
:(atan(x))

```
"""
function mcall(expr::T) where T
    mexpr = MExpr(expr)
    return convert(T, mcall(mexpr))
end

function ==(m::MExpr, n::MExpr)
    r = split(m).str
    s = split(n).str
    l = length(r)
    l != length(s) && (return false)
    b = true
    for j ∈ 1:l
        out = mcall("is($(r[j]) = $(s[j]))")
        b &= !occursin( "false", out)
    end
    return b
end

function getindex(m::MExpr, i)
    return MExpr("$m[$i]") |> mcall
end


text:
import TreeLas

data_dir(paths...) =
    joinpath(dirname(pathof(TreeLas)), "..", "..", "data", paths...) |> realpath

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
# This file contains all the wrapper functions that expose the C API of the cmark
# or cmark-gfm shared libraries.
#
# Each takes a libcmark :: CMarkLibrary function as an argument, defaulting to
# libcmarkgfm.

"""
    struct CMarkLibrary{LIBCM,LIBCMX}

Can be used to pass different `cmark` or `cmark-gfm` shared library products to the C calls.

It stores the full string paths to the shared libraries as type parameters. This way it is
possible to use the `ccall((:foo, LIBCM), ...)` syntax semi-dynamically, while still making
sure that the `LIBCM` is a compile-time constant and therefore that `ccall` is valid.

`LIBCM` is the path to the `libcmark` library, whereas the `LIBCMX` is the path to the
`libcmark-extensions` library, which can optionally be `nothing`, in which case extensions
are not allowed.
"""
struct CMarkLibrary{LIBCM,LIBCMX} end
function CMarkLibrary(libcmark::AbstractString, libcmark_extensions::Union{AbstractString,Nothing} = nothing)
    isfile(libcmark) || error("Invalid libcmark path: $(libcmark)")
    libcm = Symbol(libcmark)
    libcmx = if isnothing(libcmark_extensions)
        nothing
    else
        isfile(libcmark_extensions) || error("Invalid libcmark_extensions path: $(libcmark_extensions)")
        Symbol(libcmark_extensions)
    end
    CMarkLibrary{libcm,libcmx}()
end

function init!(libcmark::CMarkLibrary{LIBCM,LIBCMX}) where {LIBCM, LIBCMX}
    # If the extension shared library is not provided, we assume that there are no extensions.
    isnothing(LIBCMX) && return

    cmark_gfm_core_extensions_ensure_registered(libcmark)
    exts = findsyntaxextension.(extensions)
    any(exts .== C_NULL) && error("Failed to load GFM extensions, $exts")

    # Not entirely sure about the RTLD_* options here, but I am borrowing them from the JLL
    # library.
    h = Libdl.dlopen(LIBCMX, RTLD_LAZY | RTLD_DEEPBIND)
    CMARK_NODE_STRIKETHROUGH[] = unsafe_load(Ptr{Cint}(dlsym(h, :CMARK_NODE_STRIKETHROUGH)))
    CMARK_NODE_TABLE[] = unsafe_load(Ptr{Cint}(dlsym(h, :CMARK_NODE_TABLE)))
    CMARK_NODE_TABLE_ROW[] = unsafe_load(Ptr{Cint}(dlsym(h, :CMARK_NODE_TABLE_ROW)))
    CMARK_NODE_TABLE_CELL[] = unsafe_load(Ptr{Cint}(dlsym(h, :CMARK_NODE_TABLE_CELL)))
end

"""
    const libcmarkgfm :: CMarkLibrary

Default [`CMarkLibrary`](@ref) instance, pointing to the cmark-gfm library JLL that is a
binary dependency of this package.
"""
const libcmarkgfm = CMarkLibrary(
    cmark_gfm_jll.libcmark_gfm,
    cmark_gfm_jll.libcmark_gfm_extensions,
)

"""
    @cmarkapi <cmark ccall wrapper definition>

A decorator that defines an extra method for the `libcmark` library wrappers that uses the
default shared library ([`libcmarkgfm`](@ref)).

The functions are expected to be defined with the signature

```julia
function cmark_fn(::CMarkLibrary{LIBCM, LIBCMX}, args...; kwargs...) where {LIBCM, LIBCMX}
    ...
end
```

for some value of `args` and `kwargs`. The macro then defines an additional method with
the signature

```julia
cmark_fn(args...; kwargs...) = cmark_fn(libcmarkgfm, args...; kwargs...)
```
"""
macro cmarkapi(expr)
    # The expression should be a function definition that looks something like
    #
    #   function cmark_fn(::CMarkLibrary{LIBCM, LIBCMX}, args...) where {LIBCM, LIBCMX} ... end
    #
    # This leads to the AST that looks something like:
    #
    # (:function, [
    #     (:where, [
    #         (:call, [
    #             :cmark_fn,
    #             ::CMarkLibrary{LIBCM, LIBCMX},
    #             args...
    #         ]),
    #         :PATH
    #     ]),
    #     <function body>
    # ])
    #
    # The only thing we really need is to extract the function name:
    @assert expr.head === :function
    @assert expr.args[1].head === :where
    @assert expr.args[1].args[1].head === :call
    @assert expr.args[1].args[1].args[1] isa Symbol
    fn_name = expr.args[1].args[1].args[1]
    quote
        $(esc(expr))
        function $(esc(fn_name))(args...; kwargs...)
            $(esc(fn_name))($(esc(:libcmarkgfm)), args...; kwargs...)
        end
    end
end

# libcmark C API global constants:

const CMARK_OPT_DEFAULT = 0

# Node type enum for libcmark
# @enum cmark_node_type begin
#     # Error status
#     CMARK_NODE_NONE
#
#     # Block
#     CMARK_NODE_DOCUMENT
#     CMARK_NODE_BLOCK_QUOTE
#     CMARK_NODE_LIST
#     CMARK_NODE_ITEM
#     CMARK_NODE_CODE_BLOCK
#     CMARK_NODE_HTML_BLOCK
#     CMARK_NODE_CUSTOM_BLOCK
#     CMARK_NODE_PARAGRAPH
#     CMARK_NODE_HEADING
#     CMARK_NODE_THEMATIC_BREAK
#
#     # CMARK_NODE_FIRST_BLOCK = CMARK_NODE_DOCUMENT
#     # CMARK_NODE_LAST_BLOCK = CMARK_NODE_THEMATIC_BREAK
#
#     # Inline
#     CMARK_NODE_TEXT
#     CMARK_NODE_SOFTBREAK
#     CMARK_NODE_LINEBREAK
#     CMARK_NODE_CODE
#     CMARK_NODE_HTML_INLINE
#     CMARK_NODE_CUSTOM_INLINE
#     CMARK_NODE_EMPH
#     CMARK_NODE_STRONG
#     CMARK_NODE_LINK
#     CMARK_NODE_IMAGE
#
#     # CMARK_NODE_FIRST_INLINE = CMARK_NODE_TEXT
#     # CMARK_NODE_LAST_INLINE = CMARK_NODE_IMAGE
# end

# src/cmark-gfm.h:34:#define CMARK_NODE_TYPE_PRESENT (0x8000)
const CMARK_NODE_TYPE_PRESENT = 0x8000
# src/cmark-gfm.h:35:#define CMARK_NODE_TYPE_BLOCK (CMARK_NODE_TYPE_PRESENT | 0x0000)
const CMARK_NODE_TYPE_BLOCK = (CMARK_NODE_TYPE_PRESENT | 0x0000)
# src/cmark-gfm.h:36:#define CMARK_NODE_TYPE_INLINE (CMARK_NODE_TYPE_PRESENT | 0x4000)
const CMARK_NODE_TYPE_INLINE  = (CMARK_NODE_TYPE_PRESENT | 0x4000)

@enum cmark_node_type begin
    # Error status
    CMARK_NODE_NONE = 0x0000

    # Block
    CMARK_NODE_DOCUMENT       = CMARK_NODE_TYPE_BLOCK | 0x0001
    CMARK_NODE_BLOCK_QUOTE    = CMARK_NODE_TYPE_BLOCK | 0x0002
    CMARK_NODE_LIST           = CMARK_NODE_TYPE_BLOCK | 0x0003
    CMARK_NODE_ITEM           = CMARK_NODE_TYPE_BLOCK | 0x0004
    CMARK_NODE_CODE_BLOCK     = CMARK_NODE_TYPE_BLOCK | 0x0005
    CMARK_NODE_HTML_BLOCK     = CMARK_NODE_TYPE_BLOCK | 0x0006
    CMARK_NODE_CUSTOM_BLOCK   = CMARK_NODE_TYPE_BLOCK | 0x0007
    CMARK_NODE_PARAGRAPH      = CMARK_NODE_TYPE_BLOCK | 0x0008
    CMARK_NODE_HEADING        = CMARK_NODE_TYPE_BLOCK | 0x0009
    CMARK_NODE_THEMATIC_BREAK = CMARK_NODE_TYPE_BLOCK | 0x000a
    CMARK_NODE_FOOTNOTE_DEFINITION = CMARK_NODE_TYPE_BLOCK | 0x000b

    # Inline
    CMARK_NODE_TEXT          = CMARK_NODE_TYPE_INLINE | 0x0001
    CMARK_NODE_SOFTBREAK     = CMARK_NODE_TYPE_INLINE | 0x0002
    CMARK_NODE_LINEBREAK     = CMARK_NODE_TYPE_INLINE | 0x0003
    CMARK_NODE_CODE          = CMARK_NODE_TYPE_INLINE | 0x0004
    CMARK_NODE_HTML_INLINE   = CMARK_NODE_TYPE_INLINE | 0x0005
    CMARK_NODE_CUSTOM_INLINE = CMARK_NODE_TYPE_INLINE | 0x0006
    CMARK_NODE_EMPH          = CMARK_NODE_TYPE_INLINE | 0x0007
    CMARK_NODE_STRONG        = CMARK_NODE_TYPE_INLINE | 0x0008
    CMARK_NODE_LINK          = CMARK_NODE_TYPE_INLINE | 0x0009
    CMARK_NODE_IMAGE         = CMARK_NODE_TYPE_INLINE | 0x000a
    CMARK_NODE_FOOTNOTE_REFERENCE = CMARK_NODE_TYPE_INLINE | 0x000b
end
const CMARK_NODE_STRIKETHROUGH = Ref{Int32}(-1)
const CMARK_NODE_TABLE         = Ref{Int32}(-1)
const CMARK_NODE_TABLE_ROW     = Ref{Int32}(-1)
const CMARK_NODE_TABLE_CELL    = Ref{Int32}(-1)

@enum cmark_list_type begin
    CMARK_NO_LIST
    CMARK_BULLET_LIST
    CMARK_ORDERED_LIST
end

# typedef enum {
# CMARK_NO_DELIM,
# CMARK_PERIOD_DELIM,
# CMARK_PAREN_DELIM
# } cmark_delim_type;

# libcmark C function ccall wrappers:

"""
    cmark_version([::CMarkLibrary]) -> Int

Returns the version of the libcmark library as an integer (wrapping the `cmark_version`
function from `libcmark`).
"""
function cmark_version end
@cmarkapi function cmark_version(::CMarkLibrary{LIBCM,LIBCMX}) where {LIBCM, LIBCMX}
    version = ccall((:cmark_version, LIBCM), Cint, ())
    return Int(version)
end

"""
    cmark_version_string([::CMarkLibrary]) -> String

Returns the version of the libcmark library as a string (wrapping the `cmark_version_string`
function from `libcmark`).
"""
function cmark_version_string end
@cmarkapi function cmark_version_string(::CMarkLibrary{LIBCM,LIBCMX}) where {LIBCM, LIBCMX}
    s = ccall((:cmark_version_string, LIBCM), Cstring, ())
    return unsafe_string(s)
end

"""
    cmark_markdown_to_html([::CMarkLibrary{LIBCM,LIBCMX},] markdown::AbstractString) -> String

Wraps the `cmark_markdown_to_html` function from `libcmark`.
"""
function cmark_markdown_to_html end
@cmarkapi function cmark_markdown_to_html(::CMarkLibrary{LIBCM,LIBCMX}, markdown::AbstractString) where {LIBCM, LIBCMX}
    s = ccall(
        (:cmark_markdown_to_html, LIBCM),
        Cstring, (Cstring, Csize_t, Cint),
        markdown, length(markdown), CMARK_OPT_DEFAULT
    )
    return unsafe_string(s)
end

"""
    cmark_parse_document([::CMarkLibrary{LIBCM,LIBCMX},] markdown::AbstractString) -> Ptr{Cvoid}

Wraps the `cmark_parse_document` function from `libcmark`:

> `cmark_node *cmark_parse_document(const char *buffer, size_t len, int options)`
>
> Parse a CommonMark document in 'buffer' of length 'len'.
> Returns a pointer to a tree of nodes.  The memory allocated for
> the node tree should be released using 'cmark_node_free'
> when it is no longer needed.
"""
function cmark_parse_document end
@cmarkapi function cmark_parse_document(::CMarkLibrary{LIBCM,LIBCMX}, markdown::AbstractString) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_parse_document, LIBCM),
        Ptr{Cvoid}, (Cstring, Csize_t, Cint),
        markdown, length(markdown), CMARK_OPT_DEFAULT
    )
end

"""
    cmark_node_get_type([::CMarkLibrary,] node::Ptr{Cvoid}) -> node_type::Cint

Wraps `cmark_node_get_type` from `libcmark`:

> `cmark_node_type cmark_node_get_type(cmark_node *node)`
>
> Returns the type of node, or `CMARK_NODE_NONE` on error.
"""
function cmark_node_get_type end
@cmarkapi function cmark_node_get_type(::CMarkLibrary{LIBCM,LIBCMX}, node::Ptr{Cvoid}) where {LIBCM, LIBCMX}
    t = ccall(
        (:cmark_node_get_type, LIBCM),
        Cint, (Ptr{Cvoid},),
        node
    )
end

"""
    cmark_node_next([::CMarkLibrary,] node::Ptr{Cvoid}) -> node::Ptr{Cvoid}

Wraps `cmark_node_next` from `libcmark`:

> `cmark_node * cmark_node_next(cmark_node *node)`
>
> Returns the next node in the sequence after `node`, or `NULL` if there is none.
"""
function cmark_node_next end
@cmarkapi function cmark_node_next(::CMarkLibrary{LIBCM,LIBCMX}, node::Ptr{Cvoid}) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_node_next, LIBCM),
        Ptr{Cvoid}, (Ptr{Cvoid},),
        node
    )
end

"""
    cmark_node_first_child([::CMarkLibrary,] node::Ptr{Cvoid}) -> node::Ptr{Cvoid}

Wraps `cmark_node_first_child` from `libcmark`:

> `cmark_node * cmark_node_first_child(cmark_node *node)`
>
> Returns the first child of `node`, or `NULL` if `node` has no children.
"""
function cmark_node_first_child end
@cmarkapi function cmark_node_first_child(::CMarkLibrary{LIBCM,LIBCMX}, node::Ptr{Cvoid}) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_node_first_child, LIBCM),
        Ptr{Cvoid}, (Ptr{Cvoid},),
        node
    )
end

"""
    cmark_node_get_literal([::CMarkLibrary,] node::Ptr{Cvoid}) -> s::Cstring

Wraps `cmark_node_get_literal` from `libcmark`:

> `const char * cmark_node_get_literal(cmark_node *node)`
>
> Returns the string contents of `node`, or an empty string if none is set. Returns `NULL` if
> called on a node that does not have string content.
"""
function cmark_node_get_literal end
@cmarkapi function cmark_node_get_literal(::CMarkLibrary{LIBCM,LIBCMX}, node::Ptr{Cvoid}) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_node_get_literal, LIBCM),
        Cstring, (Ptr{Cvoid},),
        node
    )
end

"""
   cmark_node_get_heading_level([::CMarkLibrary,] node::Ptr{Cvoid}) -> level::Cint

Wraps `cmark_node_get_heading_level` from `libcmark`:

> `int cmark_node_get_heading_level(cmark_node *node)`
>
> Returns the heading level of `node`, or `0` if `node` is not a heading.
"""
function cmark_node_get_heading_level end
@cmarkapi function cmark_node_get_heading_level(::CMarkLibrary{LIBCM,LIBCMX}, node::Ptr{Cvoid}) where {LIBCM, LIBCMX}
   ccall(
       (:cmark_node_get_heading_level, LIBCM),
       Cint, (Ptr{Cvoid},),
       node
   )
end

"""
   cmark_node_get_fence_info([::CMarkLibrary,] node::Ptr{Cvoid}) -> Cstring

Wraps `cmark_node_get_fence_info` from `libcmark`:

> `const char * cmark_node_get_fence_info(cmark_node *node)`
>
> Returns the info string from a fenced code block.
"""
function cmark_node_get_fence_info end
@cmarkapi function cmark_node_get_fence_info(::CMarkLibrary{LIBCM,LIBCMX}, node::Ptr{Cvoid}) where {LIBCM, LIBCMX}
   ccall(
       (:cmark_node_get_fence_info, LIBCM),
       Cstring, (Ptr{Cvoid},),
       node
   )
end

"""
   cmark_node_get_list_type([::CMarkLibrary,] node::Ptr{Cvoid}) -> cmark_list_type

Wraps `cmark_node_get_list_type` from `libcmark`:

> `cmark_list_type cmark_node_get_list_type(cmark_node *node)`
>
> Returns the list type of `node`, or `CMARK_NO_LIST` if `node` is not a list.
"""
function cmark_node_get_list_type end
@cmarkapi function cmark_node_get_list_type(::CMarkLibrary{LIBCM,LIBCMX}, node::Ptr{Cvoid}) where {LIBCM, LIBCMX}
   ccall(
       (:cmark_node_get_list_type, LIBCM),
       Cuint, (Ptr{Cvoid},),
       node
   )
end

"""
   cmark_node_get_url([::CMarkLibrary,] node::Ptr{Cvoid}) -> cmark_list_type

Wraps `cmark_node_get_url` from `libcmark`:

> `const char * cmark_node_get_url(cmark_node *node)`
>
> Returns the URL of a link or image node, or an empty string if no URL is set. Returns
> `NULL` if called on a node that is not a link or image.
"""
function cmark_node_get_url end
@cmarkapi function cmark_node_get_url(::CMarkLibrary{LIBCM,LIBCMX}, node::Ptr{Cvoid}) where {LIBCM, LIBCMX}
   ccall(
       (:cmark_node_get_url, LIBCM),
       Cstring, (Ptr{Cvoid},),
       node
   )
end

"""
    cmark_node_get_title([::CMarkLibrary,] node::Ptr{Cvoid}) -> Cstring

Wraps `cmark_node_get_title` from `libcmark`:

> `const char * cmark_node_get_title(cmark_node *node)`
>
> Returns the title of a link or image node, or an empty string if no title is set. Returns
> `NULL` if called on a node that is not a link or image.
"""
function cmark_node_get_title end
@cmarkapi function cmark_node_get_title(::CMarkLibrary{LIBCM,LIBCMX}, node::Ptr{Cvoid}) where {LIBCM, LIBCMX}
   ccall(
       (:cmark_node_get_title, LIBCM),
       Cstring, (Ptr{Cvoid},),
       node
   )
end

"""''
> cmark_syntax_extension *cmark_find_syntax_extension(const char *name);
"""
function cmark_find_syntax_extension end
@cmarkapi function cmark_find_syntax_extension(::CMarkLibrary{LIBCM,LIBCMX}, name::AbstractString) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_find_syntax_extension, LIBCM),
        Ptr{Cvoid}, (Cstring,),
        name
    )
end

"""
cmark_parser * cmark_parser_new(int options)

Creates a new parser object.
"""
function cmark_parser_new end
@cmarkapi function cmark_parser_new(::CMarkLibrary{LIBCM,LIBCMX}) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_parser_new, LIBCM),
        Ptr{Cvoid}, (Cint,),
        CMARK_OPT_DEFAULT
    )
end

"""
void cmark_parser_free(cmark_parser *parser)

Frees memory allocated for a parser object.
"""
function cmark_parser_free end
@cmarkapi function cmark_parser_free(::CMarkLibrary{LIBCM,LIBCMX}, p::Ptr{Cvoid}) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_parser_free, LIBCM),
        Cvoid, (Ptr{Cvoid},),
        p
    )
end

"""
void cmark_parser_feed(cmark_parser *parser, const char *buffer, size_t len)

Feeds a string of length len to parser.
"""
function cmark_parser_feed end
@cmarkapi function cmark_parser_feed(::CMarkLibrary{LIBCM,LIBCMX}, p::Ptr{Cvoid}, s::AbstractString) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_parser_feed, LIBCM),
        Cvoid, (Ptr{Cvoid}, Cstring, Csize_t),
        p, s, length(s)
    )
end

"""
cmark_node * cmark_parser_finish(cmark_parser *parser)

Finish parsing and return a pointer to a tree of nodes.
"""
function cmark_parser_finish end
@cmarkapi function cmark_parser_finish(::CMarkLibrary{LIBCM,LIBCMX}, p::Ptr{Cvoid}) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_parser_finish, LIBCM),
        Ptr{Cvoid}, (Ptr{Cvoid},),
        p
    )
end

# int cmark_parser_attach_syntax_extension(cmark_parser *parser, cmark_syntax_extension *extension);
@cmarkapi function cmark_parser_attach_syntax_extension(::CMarkLibrary{LIBCM,LIBCMX}, p::Ptr{Cvoid}, ext::Ptr{Cvoid}) where {LIBCM, LIBCMX}
    ccall(
        (:cmark_parser_attach_syntax_extension, LIBCM),
        Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
        p, ext
    )
end

# Extensions
#void cmark_gfm_core_extensions_ensure_registered(void);
@cmarkapi function cmark_gfm_core_extensions_ensure_registered(::CMarkLibrary{LIBCM,LIBCMX}) where {LIBCM, LIBCMX}
    ccall((:cmark_gfm_core_extensions_ensure_registered, LIBCMX), Cvoid, ())
end


text:
using CPLEX, Base.Test, MathOptInterface, MathOptInterface.Test, MathOptInterfaceCPLEX

const MOIT = MathOptInterface.Test
const MOICPX = MathOptInterfaceCPLEX

@testset "MathOptInterfaceCPLEX" begin
    @testset "Linear tests" begin
        linconfig = MOIT.TestConfig()
        solver = CPLEXOptimizer()
        MOIT.contlineartest(solver, linconfig, ["linear10","linear12","linear8a","linear8b","linear8c"])

        solver_nopresolve = CPLEXOptimizer(CPX_PARAM_SCRIND=0, CPX_PARAM_PREIND=0)
        MOIT.contlineartest(solver_nopresolve , linconfig, ["linear10","linear12","linear8b","linear8c"])

        linconfig_nocertificate = MOIT.TestConfig(infeas_certificates=false)
        MOIT.linear12test(solver, linconfig_nocertificate)
        MOIT.linear8ctest(solver, linconfig_nocertificate)
        MOIT.linear8btest(solver, linconfig_nocertificate)

        linconfig_noquery = MOIT.TestConfig(query=false)
        MOIT.linear10test(solver, linconfig_noquery)
    end

    @testset "Quadratic tests" begin
        quadconfig = MOIT.TestConfig(atol=1e-3, rtol=1e-4, duals=false, query=false)
        solver = CPLEXOptimizer()
        MOIT.contquadratictest(solver, quadconfig)
    end

    @testset "Linear Conic tests" begin
        linconfig = MOIT.TestConfig()
        solver = CPLEXOptimizer()
        MOIT.lintest(solver, linconfig, ["lin3","lin4"])

        solver_nopresolve = CPLEXOptimizer(CPX_PARAM_SCRIND=0, CPX_PARAM_PREIND=0)
        MOIT.lintest(solver_nopresolve, linconfig)
    end

    @testset "Integer Linear tests" begin
        intconfig = MOIT.TestConfig()
        solver = CPLEXOptimizer()
        MOIT.intlineartest(solver, intconfig)
    end

    @testset "ModelLike tests" begin
        intconfig = MOIT.TestConfig()
        solver = CPLEXOptimizer()
        MOIT.validtest(solver)
        MOIT.emptytest(solver)
        solver2 = CPLEXOptimizer()
        MOIT.copytest(solver,solver2)
    end
end
;


text:
export sampleTangent
export samplePoint
"""
    $SIGNATURES

Return a random sample as a tangent vector from a belief represented by coordinates on a manifold at point p.

Notes

"""
function sampleTangent end

# Sampling MKD
function sampleTangent(M::AbstractGroupManifold, x::ManifoldKernelDensity, p=mean(x))
    # get legacy matrix of coordinates and selected labels
    coords, lbls = sample(x.belief,1)
    X = hat(x.manifold, p, coords)
    return X
end

function sampleTangent(x::ManifoldKernelDensity, p=mean(x))
  return sampleTangent(x.manifold, x, p)
end

# Sampling Distributions
function sampleTangent(M::AbstractManifold, z::Distribution, p, basis::AbstractBasis)
  return get_vector(M, p, rand(z), basis)
end

function sampleTangent(M::AbstractGroupManifold, z::Distribution, p=identity_element(M))
  return hat(M, p, rand(z,1)[:]) #TODO find something better than (z,1)[:]
end

"""
    $SIGNATURES

Return a random sample point on a manifold from a belief represented by coordinates at point p.

Notes

"""
function samplePoint(M::AbstractManifold, sbelief, p, basis::AbstractBasis, retraction_method::AbstractRetractionMethod=ExponentialRetraction())
  X = sampleTangent(M, sbelief, p, basis)
  return retract(M, p, X, retraction_method)
end
function samplePoint(M::AbstractGroupManifold, sbelief, p=identity_element(M), retraction_method::AbstractRetractionMethod=ExponentialRetraction())
  X = sampleTangent(M, sbelief, p)
  return retract(M, p, X, retraction_method)
end

function samplePoint(M::AbstractGroupManifold, sbelief::ManifoldKernelDensity, p=identity_element(M, mean(sbelief)), retraction_method::AbstractRetractionMethod=ExponentialRetraction())
  X = sampleTangent(M, sbelief, p)
  return retract(M, p, X, retraction_method)
end

function samplePoint(x::ManifoldKernelDensity, p=mean(x))
  return samplePoint(x.manifold, x, p)
end

# FIXME: rather use manifolds
function samplePoint(distr::SamplableBelief)
  Base.depwarn("samplePoint(distr::SamplableBelief) should be replaced by samplePoint(M<:AbstractManifold, distr::SamplableBelief, ...)", :samplePoint)
  rand(distr,1)
end

## default getSample
# getSample(cf::CalcFactor{<:AbstractPrior}, N::Int=1) = ([samplePoint(getManifold(cf.factor), cf.factor.Z, ) for _=1:N], )
# getSample(cf::CalcFactor{<:AbstractRelative}, N::Int=1) = ([sampleTangent(getManifold(cf.factor), cf.factor.Z) for _=1:N], )


text:
function tcm_estimateMMNL(N, K, assortments, transaction_counts, numInitializations, time_limit, perTol)

	results_class_probs = Array{Array{Float64}, 1}(numInitializations)
	results_u = Array{Array{Float64}, 1}(numInitializations)
	results_loglik = zeros(numInitializations)

	#perTol = 1e-7;

	grand_start_time = time()

	for i in 1:numInitializations
		start_time = time()

		u = 10*(rand(K,N) - 0.5)
		u[:,N] = 0
		class_probs = 1/K * ones(K) #rand(K)
		#class_probs = class_probs / sum(class_probs)

		v, v_by_class = tcm_predictMMNL(N, K, class_probs, u, assortments)
		logv = log.(v)
		logv[ isinf.(logv)] = 0.0;
		loglik = dot( convert(Array{Float64,}, transaction_counts), logv)
		iter = 0;
		println("Init ", i, " -- iter ", iter, " MMNL iteration -- time elapsed = ", time() - start_time, " log likelihood = ",loglik)

		while ( (time() - start_time) < time_limit )
			iter += 1;
			h, v, v_by_class = tcm_MMNL_E(N, K, class_probs, u, assortments)
			u = tcm_MMNL_M_update_u(N, K, h, assortments, transaction_counts)
			class_probs = tcm_MMNL_M_update_class_probs(K, h, transaction_counts)

			v, v_by_class = tcm_predictMMNL(N, K, class_probs, u, assortments)
			logv = log.(v)
			logv[ isinf.(logv)] = 0.0;
			new_loglik = dot( convert(Array{Float64,}, transaction_counts), logv)

	        println("Init ", i, " -- iter ", iter, " MMNL iteration -- time elapsed = ", time() - start_time, " log likelihood = ",loglik)

	        if ( (new_loglik - loglik) / abs(loglik) < perTol )
	        	temp = (new_loglik - loglik) / abs(loglik)
	        	@show temp
	        	@show new_loglik
	        	@show loglik

	        	loglik = new_loglik
				break
			end

			loglik = new_loglik
		end

		results_loglik[i] = loglik
		results_u[i] = u
		results_class_probs[i] = class_probs
	end

	grand_end_time = time()

	best_i = indmax(results_loglik)
	u = results_u[best_i]
	class_probs = results_class_probs[best_i]
	loglik = results_loglik[best_i]

	grand_elapsed_time = grand_end_time - grand_start_time

	return u, class_probs, loglik, grand_elapsed_time

end

function tcm_MMNL_M_update_u(N, K, h, assortments, transaction_counts)
	u = zeros(K,N)

	for k in 1:K
		# @show h[k]
		# @show size(h[k])
		#@show size(h[k])
		#@show size(transaction_counts)
		weighted_transaction_counts = transaction_counts .* h[k]
		single_u, obj, elapsed_time = tcm_estimateMNL(N, assortments, weighted_transaction_counts)
		#@show single_u
		u[k,1:N-1] = single_u
	end

	return u
end

function tcm_MMNL_M_update_class_probs(K, h, transaction_counts)
	sum_of_h = zeros(K)

	for k in 1:K
		sum_of_h[k] = dot( h[k], convert(Array{Float64}, transaction_counts))
	end

	class_probs = sum_of_h ./ sum(sum_of_h)

	return class_probs
end

function tcm_MMNL_E(N, K, class_probs, u, assortments)

	v, v_by_class = tcm_predictMMNL(N, K, class_probs, u, assortments)

	#@show size(v)
	#@show size(v_by_class[1])

	h = Array{ Array{Float64,2}, 1}(K)
	for k in 1:K
		single_h = class_probs[k] * v_by_class[k] ./ v
		single_h[ isnan.(single_h)] = 0.0;

		h[k] = single_h
	end

	return h, v, v_by_class
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
function numerical_gradient(feval, x, dx=0.00001)

  g = zeros(size(x))

  for ii = 1:length(x)

    x[ii] = x[ii] + dx

    D1 = feval(x)

    x[ii] = x[ii] - 2*dx

    D2 = feval(x)

    x[ii] = x[ii] + dx

    g[ii] = (D1-D2)/(2*dx)

  end

  return g

end


text:
using CharFuncPricing, TaylorSeries

@testset "CumulantsJoshi" begin
    κ = 4.0
    θ = 0.25
    σ = 1.0
    ρ = -0.5
    v0 = 0.01
    τ = 0.01
    params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc{HestonParams{Float64},Taylor1{Complex}}(params)
    t = Taylor1(Float64, 4)
    cft = CharFuncPricing.evaluateLogCharFuncCui(cf, t, τ)
    c1, c2, c4 = computeCumulants(params, τ)
    println(c1, " ", c2 / 2, " ", c4 / (2 * 3 * 4), " ", cft)
    @test isapprox(imag(cft.coeffs[2]), c1, atol = 1e-12)
    @test isapprox(-real(cft.coeffs[3]) * 2, c2, atol = 1e-12)
    @test isapprox(cft.coeffs[5] * 2 * 3 * 4, c4, atol = 1e-16)
end
@testset "CumulantsAlan" begin
    κ = 4.0
    θ = 0.25
    σ = 1.0
    ρ = -0.5
    v0 = 0.04
    τ = 1.0
    params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc{HestonParams{Float64},Taylor1{Complex}}(params)
    t = Taylor1(Float64, 4)
    cft = CharFuncPricing.evaluateLogCharFuncCui(cf, t, τ)
    c1, c2, c4 = computeCumulants(params, τ)
    println(c1, " ", c2 / 2, " ", c4 / (2 * 3 * 4), " ", cft)
    @test isapprox(imag(cft.coeffs[2]), c1, atol = 1e-12)
    @test isapprox(-real(cft.coeffs[3]) * 2, c2, atol = 1e-12)
    @test isapprox(cft.coeffs[5] * 2 * 3 * 4, c4, atol = 1e-15)
end
@testset "CumulantsNearZeroKappa" begin
    #Low accuracy with small kappa due to powers of kappa in the denominator.
    κ = 0.02
    θ = 0.25
    σ = 1.0
    ρ = -0.5
    v0 = 0.04
    τ = 1.0
    params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc{HestonParams{Float64},Taylor1{Complex}}(params)
    t = Taylor1(Float64, 4)
    cft = CharFuncPricing.evaluateLogCharFuncCui(cf, t, τ)
    c1, c2, c4 = computeCumulants(params, τ)
    println(c1, " ", c2 / 2, " ", c4 / (2 * 3 * 4), " ", cft)
    @test isapprox(imag(cft.coeffs[2]), c1, atol = 1e-12)
    @test isapprox(-real(cft.coeffs[3]) * 2, c2, atol = 1e-3)
    @test isapprox(cft.coeffs[5] * 2 * 3 * 4, c4, atol = 1e-3)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ fc872510-e7b7-11eb-0f83-ed2e18af538f
md"""
# Chapter 10
#### Nicetak
 $(import Dates; Dates.format(Dates.today(), Dates.DateFormat("U d, Y")))
"""

# ╔═╡ a6075125-eb1f-4486-9e51-185ff559dbcd
md"""
## 10.1
頂点数が $N$ の二分木において, 高さが $N-1$ となるような例をあげてください.
"""

# ╔═╡ 12155e72-f81e-4506-bb4a-087144245cee
md"""
**答.** 北斗神拳の継承者(一子相伝のため)
"""

# ╔═╡ a1da9aaa-f90b-4f67-95df-f79ff214dad5
md"""
## 10.2
空のヒープに3個の整数5, 6, 1をこの順に挿入して得られるヒープについて, 配列で表した時の様子を示してください.

"""

# ╔═╡ f05593a8-cdd7-4ae0-bfe9-6d70a17c356d
md"""
**答.** [6, 5, 1]
"""

# ╔═╡ 0a789e1a-f359-445a-bacf-f5bebdbf91f3
md"""
## 10.3
空のヒープに7個の整数5, 6, 1, 2, 7, 3, 4をこの順に挿入して得られるヒープについて, 配列で表したときの様子を示してください.
"""

# ╔═╡ 7a15992c-8c54-41cf-bd3f-4bb657b922ed
md"""
**答.** [7, 6, 4, 2, 5, 1, 3]
"""

# ╔═╡ 2c6c91c2-9c7d-4ae2-9d0d-f66e891845bd
md"""
## 10.4
強平衡二分木の高さが$h = O(\log N)$となることを示してください.
"""

# ╔═╡ e99a2952-e0be-4afb-8ee1-6438828ceb01
md"""
**答.** 高さ$h$段目に入る要素数は$2^h$であり, 強平衡二分木は各階層が全て埋まったら次の段を埋め始める. したがって,

$2^h \le N < 2^{h+1}$

よって, $h = O(\log N)$.
"""

# ╔═╡ 89478c0a-9a09-4bf6-96a4-d166d184e126
md"""
## 10.5
頂点数が $N$ の木の辺数が $N-1$であることを示してください.
"""

# ╔═╡ c015fddb-a32c-4ade-ba0c-336b0197fc00
md"""
**答.**
$N = 1$の時は自明. 頂点数 $N$ の木の辺数が $N- 1$ と仮定する.
任意の頂点数 $N+1$ の木は, 頂点数 $N$の木に一つの頂点を加えたものである.
この操作の際に増える辺の数は, 木の定義から, 1である.
したがって, 頂点数 $N+1$ の木の辺数は $N$ であり, 帰納法から命題は示された.
"""

# ╔═╡ Cell order:
# ╟─fc872510-e7b7-11eb-0f83-ed2e18af538f
# ╟─a6075125-eb1f-4486-9e51-185ff559dbcd
# ╟─12155e72-f81e-4506-bb4a-087144245cee
# ╟─a1da9aaa-f90b-4f67-95df-f79ff214dad5
# ╟─f05593a8-cdd7-4ae0-bfe9-6d70a17c356d
# ╟─0a789e1a-f359-445a-bacf-f5bebdbf91f3
# ╟─7a15992c-8c54-41cf-bd3f-4bb657b922ed
# ╟─2c6c91c2-9c7d-4ae2-9d0d-f66e891845bd
# ╟─e99a2952-e0be-4afb-8ee1-6438828ceb01
# ╟─89478c0a-9a09-4bf6-96a4-d166d184e126
# ╟─c015fddb-a32c-4ade-ba0c-336b0197fc00


text:
module Operators

using Compat
using Compat.Test
using EAGOSmoothMcCormickGrad
using IntervalArithmetic
using StaticArrays

function about(calc,val,tol)
    return (val - tol <= calc <= val + tol)
end

################################################################################
######################## Tests Standard McCormick Relaxations ##################
################################################################################
EAGOSmoothMcCormickGrad.set_diff_relax(0)

######## tests division of same object ######
a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(1.0,5.0);Interval(-1.0,2.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(3.0,3.0,a,a,xIBox[1],false,xIBox,mBox)
out = X/X
@test out == 1.0

######## tests nonsmooth division ######
a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-3.0,4.0);Interval(-5.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-2.0,-2.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(-4.0,-4.0,b,b,xIBox[2],false,xIBox,mBox)
out = X/Y
@test about(out.cc,0.6,1E-6)
@test about(out.cv,0.41666666,1E-6)
@test about(out.cc_grad[1],-0.2,1E-6)
@test about(out.cc_grad[2],0.2,1E-6)
@test about(out.cv_grad[1],-0.333333,1E-6)
@test about(out.cv_grad[2],0.1875,1E-6)
@test about(out.Intv.lo,-1.33333333,1E-6)
@test about(out.Intv.hi,1.0,1E-6)

######## tests exponent on product ######
a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-3.0,4.0);Interval(-5.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-2.0,-2.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(-4.0,-4.0,b,b,xIBox[2],false,xIBox,mBox)
out = exp(X*Y)
@test about(out.cc,2.708614394334035e6,1E-1)
@test about(out.cv,1096.6331584284585,1E-1)
@test about(out.cc_grad[1],-2.80201e5,1E1)
@test about(out.cc_grad[2],-2.80201e5,1E1)
@test about(out.cv_grad[1],-5483.17,1E-1)
@test about(out.cv_grad[2],-3289.9,1E-1)
@test about(out.Intv.lo,2.06115e-09,1E-1)
@test about(out.Intv.hi,3.26902e+06,1E1)

######## tests log on product ######
a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(3.0,7.0);Interval(3.0,9.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.0,4.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(7.0,7.0,b,b,xIBox[2],false,xIBox,mBox)
out = log(X*Y)
@test about(out.cc,3.4011973816621555,1E-8)
@test about(out.cv,2.7377551742960287,1E-8)
@test about(out.cc_grad[1],0.3,1E-1)
@test about(out.cc_grad[2],0.1,1E-1)
@test about(out.cv_grad[1],0.108106,1E-5)
@test about(out.cv_grad[2],0.108106,1E-5)
@test about(out.Intv.lo,2.19722,1E-4)
@test about(out.Intv.hi,4.14314,1E-4)

# tests powers (square)
a = seed_g(Float64,1,2)
xIBox = [Interval(3.0,7.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.0,4.0,a,a,xIBox[1],false,xIBox,mBox)
out = X^2
@test about(out.cc,19,1E-8)
@test about(out.cv,16,1E-8)
@test about(out.cc_grad[1],10.0,1E-1)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],8.0,1E-5)
@test about(out.cv_grad[2],0.0,1E-5)
@test about(out.Intv.lo,9,1E-4)
@test about(out.Intv.hi,49,1E-4)

# tests powers (^-2 on positive domain)
a = seed_g(Float64,1,2)
xIBox = [Interval(3.0,7.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.0,4.0,a,a,xIBox[1],false,xIBox,mBox)
out = X^(-2)
@test about(out.cc,0.08843537414965986,1E-8)
@test about(out.cv,0.0625,1E-8)
@test about(out.cc_grad[1],-0.0226757,1E-1)
@test about(out.cc_grad[2],0.0,1E-4)
@test about(out.cv_grad[1],-0.03125,1E-5)
@test about(out.cv_grad[2],0.0,1E-4)
@test about(out.Intv.lo,0.0204081,1E-4)
@test about(out.Intv.hi,0.111112,1E-4)

# tests powers (^-2 on negative domain)
a = seed_g(Float64,1,2)
xIBox = [Interval(-8.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.5,-4.5,a,a,xIBox[1],false,xIBox,mBox)
out = X^(-2)
@test about(out.cc,0.08246527777777776,1E-8)
@test about(out.cv,0.04938271604938271,1E-8)
@test about(out.cc_grad[1],0.0190972,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],0.0219479,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,0.015625,1E-4)
@test about(out.Intv.hi,0.111112,1E-4)

# tests powers (^0)
a = seed_g(Float64,1,2)
xIBox = [Interval(-8.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.5,-4.5,a,a,xIBox[1],false,xIBox,mBox)
out = X^(0)
@test out == 1.0

# tests powers (^1)
a = seed_g(Float64,1,2)
xIBox = [Interval(-8.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.5,-4.5,a,a,xIBox[1],false,xIBox,mBox)
out = X^(1)
@test about(out.cc,-4.5,1E-8)
@test about(out.cv,-4.5,1E-8)
@test about(out.cc_grad[1],1.0,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],1.0,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,-8.0,1E-4)
@test about(out.Intv.hi,-3.0,1E-4)

# tests powers (^2)
a = seed_g(Float64,1,2)
xIBox = [Interval(-8.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.5,-4.5,a,a,xIBox[1],false,xIBox,mBox)
out = X^(2)
@test about(out.cc,25.5,1E-8)
@test about(out.cv,20.25,1E-8)
@test about(out.cc_grad[1],-11.0,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],-9.0,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,9.0,1E-4)
@test about(out.Intv.hi,64.0,1E-4)

# tests powers (^3)
a = seed_g(Float64,1,2)
xIBox = [Interval(-8.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.5,-4.5,a,a,xIBox[1],false,xIBox,mBox)
out = X^(3)
@test about(out.cc,-91.125,1E-8)
@test about(out.cv,-172.5,1E-8)
@test about(out.cc_grad[1],60.75,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],97.0,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,-512,1E-4)
@test about(out.Intv.hi,-27,1E-4)

# tests powers (^4)
a = seed_g(Float64,1,2)
xIBox = [Interval(-8.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.5,-4.5,a,a,xIBox[1],false,xIBox,mBox)
out = X^(4)
@test about(out.cc,1285.5,1E-8)
@test about(out.cv,410.0625,1E-8)
@test about(out.cc_grad[1],-803.0,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],-364.5,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,81,1E-4)
@test about(out.Intv.hi,4096,1E-4)

# tests sqrt
a = seed_g(Float64,1,2)
xIBox = [Interval(3.0,9.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.5,4.5,a,a,xIBox[1],false,xIBox,mBox)
out = sqrt(X)
@test about(out.cc,2.1213203435596424,1E-8)
@test about(out.cv,2.049038105676658,1E-8)
@test about(out.cc_grad[1],0.235702,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],0.211325,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,1.73205,1E-4)
@test about(out.Intv.hi,3,1E-4)

# tests powers (^3 greater than zero ISSUE WITH CC)
a = seed_g(Float64,1,2)
xIBox = [Interval(3.0,8.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.5,4.5,a,a,xIBox[1],false,xIBox,mBox)
out = X^(3)
@test about(out.cc,172.5,1E-8)
@test about(out.cv,91.125,1E-8)
@test about(out.cc_grad[1],97.0,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],60.75,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,27,1E-4)
@test about(out.Intv.hi,512,1E-4)

# tests powers (^4 greater than zero)

a = seed_g(Float64,1,2)
xIBox = [Interval(3.0,8.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.5,4.5,a,a,xIBox[1],false,xIBox,mBox)
out = X^(4)
@test about(out.cc,1285.5,1E-1)
@test about(out.cv,410.0625,1E-1)
@test about(out.cc_grad[1],803.0,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],364.5,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,81,1E-4)
@test about(out.Intv.hi,4096,1E-4)

# tests powers (^4 zero in range)

a = seed_g(Float64,1,2)
xIBox = [Interval(-8.0,3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.5,-4.5,a,a,xIBox[1],false,xIBox,mBox)
out = X^(4)
@test about(out.cc,2818.5,1)
@test about(out.cv,410.0625,1)
@test about(out.cc_grad[1],-365.0,1)
@test about(out.cc_grad[2],0.0,1)
@test about(out.cv_grad[1],-364.5,1)
@test about(out.cv_grad[2],0.0,1)
@test about(out.Intv.lo,0,1)
@test about(out.Intv.hi,4096,1)

a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-3.0,4.0);Interval(-5.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-2.0,-2.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(-4.0,-4.0,b,b,xIBox[2],false,xIBox,mBox)
out = exp2(X*Y)
@test about(out.cc,27150.62857159206,1E-1)
@test about(out.cv,128.0,1E-1)
@test about(out.cc_grad[1],-2808.69,1E1)
@test about(out.cc_grad[2],-2808.69,1E1)
@test about(out.cv_grad[1],-443.614,1E-1)
@test about(out.cv_grad[2],-266.169,1E-1)
@test about(out.Intv.lo,9.53674e-07,1E-1)
@test about(out.Intv.hi,32768,1E1)

a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-3.0,4.0);Interval(-5.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-2.0,-2.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(-4.0,-4.0,b,b,xIBox[2],false,xIBox,mBox)
out = exp10(X*Y)
@test about(out.cc,8.285714285714286e14,1E-1)
@test about(out.cv,1.0e7,1E-1)
@test about(out.cc_grad[1],-8.57143e13,1E8)
@test about(out.cc_grad[2],-8.57143e13,1E8)
@test about(out.cv_grad[1],-1.15129e8,1E3)
@test about(out.cv_grad[2],-6.90776e7,1E3)
@test about(out.Intv.lo,9.99999e-21,1E-1)
@test about(out.Intv.hi,1e+15,1E1)

a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(3.0,7.0);Interval(3.0,9.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.0,4.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(7.0,7.0,b,b,xIBox[2],false,xIBox,mBox)
out = log2(X*Y)
@test about(out.cc,4.906890595608519,1E-7)
@test about(out.cv,3.94974581312498,1E-7)
@test about(out.cc_grad[1],0.432809,1E-5)
@test about(out.cc_grad[2],0.14427,1E-5)
@test about(out.cv_grad[1],0.155964,1E-5)
@test about(out.cv_grad[2],0.155964,1E-5)
@test about(out.Intv.lo,3.16992,1E-5)
@test about(out.Intv.hi,5.97728,1E-5)

a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(3.0,7.0);Interval(3.0,9.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.0,4.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(7.0,7.0,b,b,xIBox[2],false,xIBox,mBox)
out = log10(X*Y)
@test about(out.cc,1.4771212547196624,1E-7)
@test about(out.cv,1.1889919649988407,1E-7)
@test about(out.cc_grad[1],0.130288,1E-5)
@test about(out.cc_grad[2],0.0434294,1E-5)
@test about(out.cv_grad[1],0.0469499,1E-5)
@test about(out.cv_grad[2],0.0469499,1E-5)
@test about(out.Intv.lo,0.954242,1E-5)
@test about(out.Intv.hi,1.79935,1E-5)

a = seed_g(Float64,1,2)
xIBox = [Interval(-3.0,8.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(4.5,4.5,a,a,xIBox[1],false,xIBox,mBox)
out = abs(X)
@test about(out.cc,6.409090909090908,1E-1)
@test about(out.cv,4.5,1E-1)
@test about(out.cc_grad[1],0.454545,1E-4)
@test about(out.cc_grad[2],0.0,1E-1)
@test about(out.cv_grad[1],1.0,1E-4)
@test about(out.cv_grad[2],0.0,1E-1)
@test about(out.Intv.lo,0,1E-4)
@test about(out.Intv.hi,8,1E-4)

end

# tests powers (^3 zero in range, ISSUE) TO DO
#=
a = seed_g(Float64,1,2)
xIBox = [Interval(-8.0,3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.5,-4.5,a,a,xIBox[1],false,xIBox,mBox)
out2 = X^(3)
ans1 = about(out.cc,-91.125,1E-8)
ans2 = about(out.cv,-172.5,1E-8)
ans3 = about(out.cc_grad[1],60.75,1E-4)
ans4 = about(out.cc_grad[2],0.0,1E-1)
ans5 = about(out.cv_grad[1],97.0,1E-4)
ans6 = about(out.cv_grad[2],0.0,1E-1)
ans7 = about(out.Intv.lo,-512,1E-4)
ans8 = about(out.Intv.hi,-27,1E-4)
=#


text:
import GDAL
using ArchGDAL, GDAL_jll, PROJ_jll

export rasterize, readraster, saveTIFF

function rasterize_AG(infile::String, outfile::String, options::Vector{<:AbstractString})
    ArchGDAL.read(infile) do dataset
        GDAL.close(GDAL.rasterize(
            outfile,
            Ptr{GDAL.GDALDatasetH}(C_NULL),
            dataset.ptr,
            GDAL.rasterizeoptionsnew(options, C_NULL), C_NULL))
    end
end

# works (creates the TIFF and saves it) but then crashes
function rasterize_AG2(infile::String, outfile::String, options::Vector{<:AbstractString})
    ArchGDAL.read(infile) do dataset
        ArchGDAL.unsafe_gdalrasterize(dataset, options, dest=outfile)
    end
end

# uses the command line version instead (gdal_rasterize)
# significantly faster for some reason, also gives a simple progress indication
function rasterize(infile::String, outfile::String, options::Vector{<:AbstractString}; sql::String="")
    gdal_rasterize_path() do gdal_rasterize
        if isempty(sql)
            run(`$gdal_rasterize $options $infile $outfile`)
        else
            run(`$gdal_rasterize $options -sql $sql $infile $outfile`)
        end
    end
end

function getextent(geotransform::Vector{Float64}, rastersize::Tuple{Int,Int})
    @assert length(geotransform) == 6 "A GeoTransform vector must have 6 elements."
    left, xres, _, top, _, yres = geotransform
    width, height = rastersize
    bottom, right = top+yres*height, left+xres*width
    return [left, bottom, right, top]
end

function read3draster(infile::String)
    ArchGDAL.read(infile) do dataset
        ArchGDAL.read(dataset)
    end
end

function readraster(infile::String, extentflag::Symbol, dim::Int=1)
    local raster, geotransform
    raster = ArchGDAL.read(infile) do dataset
        # display(ArchGDAL.getproj(dataset))
        geotransform = ArchGDAL.getgeotransform(dataset)
        ArchGDAL.read(dataset)[:,:,dim]
    end
    coordextent = getextent(geotransform, size(raster))
    if extentflag == :extend_to_full_globe
        left, bottom, right, top = coordextent
        xres, yres = geotransform[2], geotransform[6]
        newwidth, newheight = round.(Int, (360/xres, -180/yres))
        xindexes = 1+round(Int, (left-(-180))/xres):newwidth-round(Int, (right-180)/xres)
        yindexes = 1+round(Int, (top-90)/yres):newheight+round(Int, (bottom-(-90))/yres)
        adjusted = zeros(eltype(raster), (newwidth, newheight))
        adjusted[xindexes, yindexes] = raster
        return adjusted, coordextent
    else    # extentflag == :getextent
        return raster, coordextent
    end
end

readraster(infile::String, dim::Int=1) = readraster(infile, :none, dim)[1]

function saveTIFF(x::AbstractMatrix, filename::String, extent::Vector{Float64})
    # EPSG:4326 (switch to importEPSG later)
    # http://yeesian.com/ArchGDAL.jl/latest/projections/#Creating-Spatial-References-1
    wkt_string = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]"
    width, height = size(x)
    xres = (extent[3]-extent[1])/width
    yres = (extent[4]-extent[2])/height
    raster = ArchGDAL.create(
        filename,
        driver = ArchGDAL.getdriver("GTiff"),
        width = width,
        height = height,
        nbands = 1,
        dtype = eltype(x),
        options = ["COMPRESS=LZW"]
    )
    ## assign the projection and transformation parameters
    ArchGDAL.setgeotransform!(raster, [extent[1], xres, 0, extent[4], 0, -yres])
    ArchGDAL.setproj!(raster, wkt_string)

    ## write the raster
    ArchGDAL.write!(
        raster,
        x,      # image to "burn" into the raster
        1,      # update band 1
    )
    ArchGDAL.destroy(raster)
    nothing
end

saveTIFF(x::AbstractMatrix, filename::String) = saveTIFF(x, filename, [-180.0, -90.0, 180.0, 90.0])

# ArchGDAL tutorial: http://www.acgeospatial.co.uk/julia-prt3/

# ogrinfo -al -so C:/Stuff/Datasets/gadm36/gadm36.shp

# rasterize_AG("C:/Stuff/Datasets/gadm36/gadm36.shp", "testtest.tif", "-a ID_0 -ts 4000 2000 -ot Byte")
# shapefile2tif("C:/Stuff/Datasets/gadm36/gadm36.shp", "Europe", "ID_0", 4300, [-11, 34, 32, 72], ')

# ogr2ogr -f CSV C:/Stuff/Julia/gadmfields012.csv -sql "select uid,id_0,name_0,id_1,name_1,id_2,name_2 from gadm36" C:/Stuff/Datasets/gadm36/gadm36.shp
# gdal_rasterize -a UID -ot Int32 -ts 5000 2500 C:\Stuff\Datasets\gadm36\gadm36.shp C:/Stuff/Julia/globtest.tif
# gdal_rasterize -a UID -ot Int32 -ts 36000 18000 -co COMPRESS=LZW C:/Stuff/Datasets/gadm36/gadm36.shp C:/Users/niclas/Downloads/globtest.tif
# gdal_rasterize -a UID -ot Int32 -tr 0.01 0.01 -te -180 -90 180 90 -co COMPRESS=LZW C:/Stuff/Datasets/gadm36/gadm36.shp C:/Users/niclas/Downloads/globtest.tif

# run(`gdal_rasterize -a UID -ot Int32 -tr 0.01 0.01 -te -180 -90 180 90 -co COMPRESS=LZW C:/Stuff/Datasets/gadm36/gadm36.shp C:/Users/niclas/Downloads/globtest.tif`)
# rasterize_AG("C:/Stuff/Datasets/gadm36/gadm36.shp", "C:/Users/niclas/Downloads/globtest3.tif", split("-a ID_0 -ot Byte -tr 0.01 0.01 -te -180 -90 180 90 -co COMPRESS=LZW"))
# rasterize("C:/Stuff/Datasets/gadm36/gadm36.shp", "C:/Users/niclas/Downloads/globtest3.tif", split("-a ID_0 -ot Byte -tr 0.01 0.01 -te -180 -90 180 90 -co COMPRESS=LZW"))

# timemem-1.0 gdal_translate -r mode -tr 0.1 0.1 -co COMPRESS=LZW gadm.tif gadmsmall.tif

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS