text:
@testset "4.4.7 (d trig)^m (a+b (c cot)^n)^p" begin
    (A, C, a, b, c, d, x, ) = @variables A C a b c d x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(b*cot(e+f*x)^n)^p=#

    #= ::Title:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+b*cot(e+f*x)^n)^p=#

    #= ::Title:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cot(e+f*x)^n)^p=#

    #= ::Title:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cot(e+f*x)^n)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cot(e+f*x)^2)^p=#

    @test_int [(A + C*cot(c + d*x)^2)/sqrt(b*tan(c + d*x)), x, 15, -(((A - C)*atan(1 - (sqrt(2)*sqrt(b*tan(c + d*x)))/sqrt(b)))/(sqrt(2)*sqrt(b)*d)) + ((A - C)*atan(1 + (sqrt(2)*sqrt(b*tan(c + d*x)))/sqrt(b)))/(sqrt(2)*sqrt(b)*d) - ((A - C)*log(sqrt(b) + sqrt(b)*tan(c + d*x) - sqrt(2)*sqrt(b*tan(c + d*x))))/(2*sqrt(2)*sqrt(b)*d) + ((A - C)*log(sqrt(b) + sqrt(b)*tan(c + d*x) + sqrt(2)*sqrt(b*tan(c + d*x))))/(2*sqrt(2)*sqrt(b)*d) - (2*b*C)/(3*d*(b*tan(c + d*x))^(3/2))]

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*cot(e+f*x)^m*(a+b*cot(e+f*x)^n)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cot(e+f*x)^m*(a+b*cot(e+f*x)^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cot(e+f*x)^m*(a+b*cot(e+f*x)^2)^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [(a + b*cot(c + d*x)^2)^1, x, 3, a*x - b*x - (b*cot(c + d*x))/d]

    @test_int [(a + b*cot(c + d*x)^2)^2, x, 4, (a - b)^2*x - ((2*a - b)*b*cot(c + d*x))/d - (b^2*cot(c + d*x)^3)/(3*d)]

    @test_int [(a + b*cot(c + d*x)^2)^3, x, 4, (a - b)^3*x - (b*(3*a^2 - 3*a*b + b^2)*cot(c + d*x))/d - ((3*a - b)*b^2*cot(c + d*x)^3)/(3*d) - (b^3*cot(c + d*x)^5)/(5*d)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [1/(a + b*cot(c + d*x)^2)^1, x, 3, x/(a - b) + (sqrt(b)*atan((sqrt(b)*cot(c + d*x))/sqrt(a)))/(sqrt(a)*(a - b)*d)]

    @test_int [1/(a + b*cot(c + d*x)^2)^2, x, 5, x/(a - b)^2 + ((3*a - b)*sqrt(b)*atan((sqrt(b)*cot(c + d*x))/sqrt(a)))/(2*a^(3/2)*(a - b)^2*d) + (b*cot(c + d*x))/(2*a*(a - b)*d*(a + b*cot(c + d*x)^2))]

    @test_int [1/(a + b*cot(c + d*x)^2)^3, x, 6, x/(a - b)^3 + (sqrt(b)*(15*a^2 - 10*a*b + 3*b^2)*atan((sqrt(b)*cot(c + d*x))/sqrt(a)))/(8*a^(5/2)*(a - b)^3*d) + (b*cot(c + d*x))/(4*a*(a - b)*d*(a + b*cot(c + d*x)^2)^2) + ((7*a - 3*b)*b*cot(c + d*x))/(8*a^2*(a - b)^2*d*(a + b*cot(c + d*x)^2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cot(e+f*x)^m*(a+b*cot(e+f*x)^2)^(p/2)*with*a-b=0=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [(1 + cot(x)^2)^(3/2), x, 4, (-(1/2))*asinh(cot(x)) - (1/2)*cot(x)*sqrt(csc(x)^2)]
    @test_int [sqrt(1 + cot(x)^2), x, 3, -asinh(cot(x))]
    @test_int [1/sqrt(1 + cot(x)^2), x, 3, -(cot(x)/sqrt(csc(x)^2))]

    @test_int [(-1 - cot(x)^2)^(3/2), x, 5, (-(1/2))*atan(cot(x)/sqrt(-csc(x)^2)) + (1/2)*cot(x)*sqrt(-csc(x)^2)]
    @test_int [sqrt(-1 - cot(x)^2), x, 4, atan(cot(x)/sqrt(-csc(x)^2))]
    @test_int [1/sqrt(-1 - cot(x)^2), x, 3, -(cot(x)/sqrt(-csc(x)^2))]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [cot(x)^3/sqrt(a + a*cot(x)^2), x, 4, -(1/sqrt(a*csc(x)^2)) - sqrt(a*csc(x)^2)/a]
    @test_int [cot(x)^2/sqrt(a + a*cot(x)^2), x, 5, cot(x)/sqrt(a*csc(x)^2) - (atanh(cos(x))*csc(x))/sqrt(a*csc(x)^2)]
    @test_int [cot(x)^1/sqrt(a + a*cot(x)^2), x, 3, 1/sqrt(a*csc(x)^2)]
    @test_int [tan(x)^1/sqrt(a + a*cot(x)^2), x, 5, atanh(sqrt(a*csc(x)^2)/sqrt(a))/sqrt(a) - 1/sqrt(a*csc(x)^2)]
    @test_int [tan(x)^2/sqrt(a + a*cot(x)^2), x, 5, cot(x)/sqrt(a*csc(x)^2) + (csc(x)*sec(x))/sqrt(a*csc(x)^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cot(e+f*x)^m*(a+b*cot(e+f*x)^2)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [cot(x)^3*sqrt(a + b*cot(x)^2), x, 6, (-sqrt(a - b))*atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b)) + sqrt(a + b*cot(x)^2) - (a + b*cot(x)^2)^(3/2)/(3*b)]
    @test_int [cot(x)^1*sqrt(a + b*cot(x)^2), x, 5, sqrt(a - b)*atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b)) - sqrt(a + b*cot(x)^2)]
    @test_int [tan(x)^1*sqrt(a + b*cot(x)^2), x, 7, sqrt(a)*atanh(sqrt(a + b*cot(x)^2)/sqrt(a)) - sqrt(a - b)*atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))]

    @test_int [cot(x)^2*sqrt(a + b*cot(x)^2), x, 7, sqrt(a - b)*atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2)) - ((a - 2*b)*atanh((sqrt(b)*cot(x))/sqrt(a + b*cot(x)^2)))/(2*sqrt(b)) - (1/2)*cot(x)*sqrt(a + b*cot(x)^2)]
    @test_int [cot(x)^0*sqrt(a + b*cot(x)^2), x, 6, (-sqrt(a - b))*atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2)) - sqrt(b)*atanh((sqrt(b)*cot(x))/sqrt(a + b*cot(x)^2))]
    @test_int [tan(x)^2*sqrt(a + b*cot(x)^2), x, 5, sqrt(a - b)*atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2)) + sqrt(a + b*cot(x)^2)*tan(x)]
    @test_int [tan(x)^4*sqrt(a + b*cot(x)^2), x, 6, (-sqrt(a - b))*atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2)) - ((3*a - b)*sqrt(a + b*cot(x)^2)*tan(x))/(3*a) + (1/3)*sqrt(a + b*cot(x)^2)*tan(x)^3]

    @test_int [cot(x)^3*(a + b*cot(x)^2)^(3/2), x, 7, (-(a - b)^(3/2))*atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b)) + (a - b)*sqrt(a + b*cot(x)^2) + (1/3)*(a + b*cot(x)^2)^(3/2) - (a + b*cot(x)^2)^(5/2)/(5*b)]
    @test_int [cot(x)^2*(a + b*cot(x)^2)^(3/2), x, 8, (a - b)^(3/2)*atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2)) - ((3*a^2 - 12*a*b + 8*b^2)*atanh((sqrt(b)*cot(x))/sqrt(a + b*cot(x)^2)))/(8*sqrt(b)) - (1/8)*(5*a - 4*b)*cot(x)*sqrt(a + b*cot(x)^2) - (1/4)*b*cot(x)^3*sqrt(a + b*cot(x)^2)]
    @test_int [cot(x)^1*(a + b*cot(x)^2)^(3/2), x, 6, (a - b)^(3/2)*atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b)) - (a - b)*sqrt(a + b*cot(x)^2) - (1/3)*(a + b*cot(x)^2)^(3/2)]
    @test_int [tan(x)^1*(a + b*cot(x)^2)^(3/2), x, 8, a^(3/2)*atanh(sqrt(a + b*cot(x)^2)/sqrt(a)) - (a - b)^(3/2)*atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b)) - b*sqrt(a + b*cot(x)^2)]
    @test_int [tan(x)^2*(a + b*cot(x)^2)^(3/2), x, 7, (a - b)^(3/2)*atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2)) - b^(3/2)*atanh((sqrt(b)*cot(x))/sqrt(a + b*cot(x)^2)) + a*sqrt(a + b*cot(x)^2)*tan(x)]

    @test_int [(a + b*cot(c + d*x)^2)^(5/2), x, 8, -(((a - b)^(5/2)*atan((sqrt(a - b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2)))/d) - (sqrt(b)*(15*a^2 - 20*a*b + 8*b^2)*atanh((sqrt(b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2)))/(8*d) - ((7*a - 4*b)*b*cot(c + d*x)*sqrt(a + b*cot(c + d*x)^2))/(8*d) - (b*cot(c + d*x)*(a + b*cot(c + d*x)^2)^(3/2))/(4*d)]
    @test_int [(a + b*cot(c + d*x)^2)^(3/2), x, 7, -(((a - b)^(3/2)*atan((sqrt(a - b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2)))/d) - ((3*a - 2*b)*sqrt(b)*atanh((sqrt(b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2)))/(2*d) - (b*cot(c + d*x)*sqrt(a + b*cot(c + d*x)^2))/(2*d)]
    @test_int [(a + b*cot(c + d*x)^2)^(1/2), x, 6, -((sqrt(a - b)*atan((sqrt(a - b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2)))/d) - (sqrt(b)*atanh((sqrt(b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2)))/d]
    @test_int [1/(a + b*cot(c + d*x)^2)^(1/2), x, 3, -(atan((sqrt(a - b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2))/(sqrt(a - b)*d))]
    @test_int [1/(a + b*cot(c + d*x)^2)^(3/2), x, 4, -(atan((sqrt(a - b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2))/((a - b)^(3/2)*d)) + (b*cot(c + d*x))/(a*(a - b)*d*sqrt(a + b*cot(c + d*x)^2))]
    @test_int [1/(a + b*cot(c + d*x)^2)^(5/2), x, 6, -(atan((sqrt(a - b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2))/((a - b)^(5/2)*d)) + (b*cot(c + d*x))/(3*a*(a - b)*d*(a + b*cot(c + d*x)^2)^(3/2)) + ((5*a - 2*b)*b*cot(c + d*x))/(3*a^2*(a - b)^2*d*sqrt(a + b*cot(c + d*x)^2))]
    @test_int [1/(a + b*cot(c + d*x)^2)^(7/2), x, 7, -(atan((sqrt(a - b)*cot(c + d*x))/sqrt(a + b*cot(c + d*x)^2))/((a - b)^(7/2)*d)) + (b*cot(c + d*x))/(5*a*(a - b)*d*(a + b*cot(c + d*x)^2)^(5/2)) + ((9*a - 4*b)*b*cot(c + d*x))/(15*a^2*(a - b)^2*d*(a + b*cot(c + d*x)^2)^(3/2)) + (b*(33*a^2 - 26*a*b + 8*b^2)*cot(c + d*x))/(15*a^3*(a - b)^3*d*sqrt(a + b*cot(c + d*x)^2))]

    @test_int [(1 - cot(x)^2)^(3/2), x, 6, (5/2)*asin(cot(x)) - 2*sqrt(2)*atan((sqrt(2)*cot(x))/sqrt(1 - cot(x)^2)) + (1/2)*cot(x)*sqrt(1 - cot(x)^2)]
    @test_int [sqrt(1 - cot(x)^2), x, 5, asin(cot(x)) - sqrt(2)*atan((sqrt(2)*cot(x))/sqrt(1 - cot(x)^2))]
    @test_int [1/sqrt(1 - cot(x)^2), x, 3, -(atan((sqrt(2)*cot(x))/sqrt(1 - cot(x)^2))/sqrt(2))]

    @test_int [(-1 + cot(x)^2)^(3/2), x, 7, (5/2)*atanh(cot(x)/sqrt(-1 + cot(x)^2)) - 2*sqrt(2)*atanh((sqrt(2)*cot(x))/sqrt(-1 + cot(x)^2)) - (1/2)*cot(x)*sqrt(-1 + cot(x)^2)]
    @test_int [sqrt(-1 + cot(x)^2), x, 6, -atanh(cot(x)/sqrt(-1 + cot(x)^2)) + sqrt(2)*atanh((sqrt(2)*cot(x))/sqrt(-1 + cot(x)^2))]
    @test_int [1/sqrt(-1 + cot(x)^2), x, 3, -(atanh((sqrt(2)*cot(x))/sqrt(-1 + cot(x)^2))/sqrt(2))]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [cot(x)^3/sqrt(a + b*cot(x)^2), x, 5, -(atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))/sqrt(a - b)) - sqrt(a + b*cot(x)^2)/b]
    @test_int [cot(x)^2/sqrt(a + b*cot(x)^2), x, 6, atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2))/sqrt(a - b) - atanh((sqrt(b)*cot(x))/sqrt(a + b*cot(x)^2))/sqrt(b)]
    @test_int [cot(x)^1/sqrt(a + b*cot(x)^2), x, 4, atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))/sqrt(a - b)]
    @test_int [tan(x)^1/sqrt(a + b*cot(x)^2), x, 7, atanh(sqrt(a + b*cot(x)^2)/sqrt(a))/sqrt(a) - atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))/sqrt(a - b)]
    @test_int [tan(x)^2/sqrt(a + b*cot(x)^2), x, 5, atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2))/sqrt(a - b) + (sqrt(a + b*cot(x)^2)*tan(x))/a]

    @test_int [cot(x)^3/(a + b*cot(x)^2)^(3/2), x, 5, -(atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))/(a - b)^(3/2)) + a/((a - b)*b*sqrt(a + b*cot(x)^2))]
    @test_int [cot(x)^2/(a + b*cot(x)^2)^(3/2), x, 4, atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2))/(a - b)^(3/2) - cot(x)/((a - b)*sqrt(a + b*cot(x)^2))]
    @test_int [cot(x)^1/(a + b*cot(x)^2)^(3/2), x, 5, atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))/(a - b)^(3/2) - 1/((a - b)*sqrt(a + b*cot(x)^2))]
    @test_int [tan(x)^1/(a + b*cot(x)^2)^(3/2), x, 8, atanh(sqrt(a + b*cot(x)^2)/sqrt(a))/a^(3/2) - atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))/(a - b)^(3/2) + b/(a*(a - b)*sqrt(a + b*cot(x)^2))]
    @test_int [tan(x)^2/(a + b*cot(x)^2)^(3/2), x, 6, atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2))/(a - b)^(3/2) + (b*tan(x))/(a*(a - b)*sqrt(a + b*cot(x)^2)) + ((a - 2*b)*sqrt(a + b*cot(x)^2)*tan(x))/(a^2*(a - b))]

    @test_int [cot(x)^3/(a + b*cot(x)^2)^(5/2), x, 6, -(atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))/(a - b)^(5/2)) + a/(3*(a - b)*b*(a + b*cot(x)^2)^(3/2)) + 1/((a - b)^2*sqrt(a + b*cot(x)^2))]
    @test_int [cot(x)^2/(a + b*cot(x)^2)^(5/2), x, 6, atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2))/(a - b)^(5/2) - cot(x)/(3*(a - b)*(a + b*cot(x)^2)^(3/2)) - ((2*a + b)*cot(x))/(3*a*(a - b)^2*sqrt(a + b*cot(x)^2))]
    @test_int [cot(x)^1/(a + b*cot(x)^2)^(5/2), x, 6, atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))/(a - b)^(5/2) - 1/(3*(a - b)*(a + b*cot(x)^2)^(3/2)) - 1/((a - b)^2*sqrt(a + b*cot(x)^2))]
    @test_int [tan(x)^1/(a + b*cot(x)^2)^(5/2), x, 9, atanh(sqrt(a + b*cot(x)^2)/sqrt(a))/a^(5/2) - atanh(sqrt(a + b*cot(x)^2)/sqrt(a - b))/(a - b)^(5/2) + b/(3*a*(a - b)*(a + b*cot(x)^2)^(3/2)) + ((2*a - b)*b)/(a^2*(a - b)^2*sqrt(a + b*cot(x)^2))]
    @test_int [tan(x)^2/(a + b*cot(x)^2)^(5/2), x, 7, atan((sqrt(a - b)*cot(x))/sqrt(a + b*cot(x)^2))/(a - b)^(5/2) + (b*tan(x))/(3*a*(a - b)*(a + b*cot(x)^2)^(3/2)) + ((7*a - 4*b)*b*tan(x))/(3*a^2*(a - b)^2*sqrt(a + b*cot(x)^2)) + ((a - 4*b)*(3*a - 2*b)*sqrt(a + b*cot(x)^2)*tan(x))/(3*a^3*(a - b)^2)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cot(e+f*x)^m*(a+b*cot(e+f*x)^3)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*cot(e+f*x)^3)^n=#

    @test_int [1/(1 + cot(x)^3), x, 7, x/2 - (1/6)*log(1 + cot(x)) + (1/3)*log(1 - cot(x) + cot(x)^2) + (1/2)*log(sin(x))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cot(e+f*x)^m*(a+b*cot(e+f*x)^4)^p=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*cot(e+f*x)^m*(a+b*cot(e+f*x)^4)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cot(e+f*x)^m*(a+b*cot(e+f*x)^4)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [cot(x)*sqrt(a + b*cot(x)^4), x, 8, (1/2)*sqrt(b)*atanh((sqrt(b)*cot(x)^2)/sqrt(a + b*cot(x)^4)) + (1/2)*sqrt(a + b)*atanh((a - b*cot(x)^2)/(sqrt(a + b)*sqrt(a + b*cot(x)^4))) - (1/2)*sqrt(a + b*cot(x)^4)]

    @test_int [cot(x)*(a + b*cot(x)^4)^(3/2), x, 9, (1/4)*sqrt(b)*(3*a + 2*b)*atanh((sqrt(b)*cot(x)^2)/sqrt(a + b*cot(x)^4)) + (1/2)*(a + b)^(3/2)*atanh((a - b*cot(x)^2)/(sqrt(a + b)*sqrt(a + b*cot(x)^4))) - (1/4)*(2*(a + b) - b*cot(x)^2)*sqrt(a + b*cot(x)^4) - (1/6)*(a + b*cot(x)^4)^(3/2)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [cot(x)/sqrt(a + b*cot(x)^4), x, 4, atanh((a - b*cot(x)^2)/(sqrt(a + b)*sqrt(a + b*cot(x)^4)))/(2*sqrt(a + b))]

    @test_int [cot(x)/(a + b*cot(x)^4)^(3/2), x, 6, atanh((a - b*cot(x)^2)/(sqrt(a + b)*sqrt(a + b*cot(x)^4)))/(2*(a + b)^(3/2)) - (a + b*cot(x)^2)/(2*a*(a + b)*sqrt(a + b*cot(x)^4))]

    @test_int [cot(x)/(a + b*cot(x)^4)^(5/2), x, 7, atanh((a - b*cot(x)^2)/(sqrt(a + b)*sqrt(a + b*cot(x)^4)))/(2*(a + b)^(5/2)) - (a + b*cot(x)^2)/(6*a*(a + b)*(a + b*cot(x)^4)^(3/2)) - (3*a^2 + b*(5*a + 2*b)*cot(x)^2)/(6*a^2*(a + b)^2*sqrt(a + b*cot(x)^4))]
    end

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.7 (d trig)^m (a+b (c cot)^n)^p.jl


text:
function simplify!(ir)
  for (v, st) in ir
    ex = st.expr
    if iscall(ex, randn) && length(ex.args) == 1
      ir[v] = xcall(rand, Normal(0, 1))
    elseif iscall(ex, rand) && length(ex.args) == 1
      ir[v] = xcall(rand, Uniform(0, 1))
    elseif iscall(ex, rand) && length(ex.args) == 2 && ex.args[2] == Bool
      ir[v] = xcall(rand, Bernoulli(1//2))
    end
  end
  return ir
end

trace(Ts...) = simplify!(Abstract.trace(Ts...))

ext: jl


text:
using SequentialLOB
using Plots
pyplot()
using PyPlot
pygui(true)

num_paths = 1; M = 800 ; T = 2299 ; p₀ = 238.75 ; L = 200 ; Δx = L / M ; λ = 1.0
D = 0.5 ; Δx = L / M ; Δt = (Δx^2) / (2.0 * D) ; σ = 1.5 ; μ = 0.03 ; nu = 0.5
α = 0.0 ; α_slob = 500.0

slob_model = SLOB(num_paths, T, p₀, M, L, D, σ, nu, α, SourceTerm(λ, μ))

lob_densities, raw_price_paths, sample_price_paths, P⁺s, P⁻s, Ps = slob_model(true, 5487)

num_xs = Int(M / 8 + 1)
time_steps_range = 0:3000:30000

x_s = repeat(collect(range(p₀ - L / 2, stop=p₀ + L / 2, length=num_xs)), inner=11)
y_s = repeat(collect(time_steps_range), outer=num_xs)
z_s = vec(transpose(lob_densities[Int.(range(1, M + 1, length=num_xs)), (time_steps_range) .+ 1, 1]))

plt = Plots.plot(
    x_s, y_s, z_s, st=:surface,
    size=(1000, 800),
        xlabel="Price" , ylabel="Time",
        zlabel="Density (φ)", camera=(45, 15),
        xguidefontsize=14,yguidefontsize=14,zguidefontsize=14,
        xtickfontsize=8,ytickfontsize=8,ztickfontsize=8,
        dpi=200, zguidefontrotation=90, formatter=:plain,
        c=:redsblues, legend=:none
)

Plots.savefig(plt, "$(dirname(@__FILE__))/lob_density_values_sample.pdf")

slob_model = SLOB(num_paths, T, p₀, M, L, D, σ, nu, α_slob, SourceTerm(λ, μ))

lob_densities, raw_price_paths, sample_price_paths, P⁺s, P⁻s, Ps = slob_model(true, 173781)

num_xs = Int(M / 8 + 1)
time_steps_range = 0:3000:30000

x_s = repeat(collect(range(p₀ - L / 2, stop=p₀ + L / 2, length=num_xs)), inner=11)
y_s = repeat(collect(time_steps_range), outer=num_xs)
z_s = vec(transpose(lob_densities[Int.(range(1, M + 1, length=num_xs)), (time_steps_range) .+ 1, 1]))

plt = Plots.plot(
    x_s, y_s, z_s, st=:surface,
    size=(1000, 800),
        xlabel="Price" , ylabel="Time",
        zlabel="Density (φ)", camera=(35, 15),
        xguidefontsize=14,yguidefontsize=14,zguidefontsize=14,
        xtickfontsize=8,ytickfontsize=8,ztickfontsize=8,
        dpi=200, zguidefontrotation=90, formatter=:plain,
        c=:redsblues, legend=:none
)

Plots.savefig(plt, "$(dirname(@__FILE__))/slob_density_values_sample.pdf")

anim = @animate for i = 1:40
    plt1 = Plots.plot(
        x_s, y_s, z_s, st=:surface,
        size=(1000, 800), title="Camera Angle=($(i * 5), 15)",
            xlabel="Price" , ylabel="Time",
            zlabel="Density (φ)", camera=(i * 5, 15),
            xguidefontsize=14,yguidefontsize=14,zguidefontsize=14,
            xtickfontsize=8,ytickfontsize=8,ztickfontsize=8,
            dpi=100, zguidefontrotation=90, formatter=:plain,
            c=:redsblues, legend=:none
    )
end
gif(anim, "$(dirname(@__FILE__))/slob_density_values.gif", fps=8)

max_stars_repo_name: GantZA/CALOBMTMD

max_issues_repo_name: GantZA/CALOBMTMD

max_forks_repo_name: GantZA/CALOBMTMD


text:
using LinearAlgebra
using HCubature
using SpecialFunctions
using PyPlot
using SparseArrays
using ProgressBars
using StatsBase
using AlphaStableDistributions
using AdvancedHMC
using Dates
using MAT
using Optim
using MCMCDiagnostics
using Dates

function difference1(X)
    N = length(X);
    M = zeros(N,N);

    for i  = 2:N
       M[i,i] = -1;
       M[i,i-1] = 1;
    end
    M[1,1] = 1;
    return M
end

function  difference2(X)
    N = length(X);
    M = zeros(N,N);

    for i  = 3:N
       M[i,i] = -1;
       M[i,i-1] = 2;
       M[i,i-2] = -1;
    end
    M[1,1] = 1; M[2,1] = 1; M[2,2] = -1;
    return M

end

function smatrix(X0,a,b)
    N = length(X0);
    dx = X0[2]-X0[1];

    M = zeros(N,N);
    for i=2:N-1
        M[i,i-1] = -a/dx^2;
        M[i,i+1] = -a/dx^2;
        M[i,i] = 2*a/dx^2 +b;

    end
    #M[1,:] = [b+2*a/dx^2 -a/dx^2 zeros(1,N-3) -a/dx^2];
    #M[N,:] = [-a/dx^2 zeros(1,N-3) -a/dx^2 b+2*a/dx^2];

    M[1,1] = 1;
    #M(N,1) = 1/dx;
    #M(N,2) = -1/dx;

    M[N,N] = 1
    return M

end

function measurementmatrix(X,MX,kernel)
    s = X[2]-X[1];
    F = zeros(size(MX)[1],size(X)[1])
    for i = 1:size(MX)[1]
        F[i,:] = s.*kernel.(X,MX[i]);
    end
    return F
end

function divide(N,NT)
    even = mod(N,NT) == 0
    k = Int64(floor(N/NT))
    t = 0
    if(even)
        q = Vector{Vector{Int64}}(undef,k)
        t = k

    else
        q = Vector{Vector{Int64}}(undef,k+1)
        t = k +1
    end

    for i = 1:t
        q[i] = Vector(NT*(i-1)+1:min(NT*i,N))
    end

    return q

end

function partialreplace(orig,replace,ix)
    a = copy(orig)
    a[ix] = replace
    return a
end

function  logpcauchy(W,arguments)
    Mat = arguments.Mat;
    noisesigma = arguments.noisesigma
    scale = arguments.scale; y  = arguments.y; #F = arguments.F
    G = 1;

    r = Mat*(W);
    #r = F*U;

    logp = -0.5/noisesigma^2*sum((r-y).^2);
    logp = logp + sum(log.(scale./(pi*(scale^2.0.+W.^2))));

    return logp;

end

function logpcauchygradi(W,arguments;both=false)
    logp = 0.0
    Mat = arguments.Mat;
    noisesigma = arguments.noisesigma
    scale = arguments.scale; y  = arguments.y; #F = arguments.F

    r = Mat*(W);
    #r = F*U;

    if(both)
        logp = -0.5/noisesigma^2*sum((r-y).^2);
        logp = logp + sum(log.(scale./(pi*(scale^2.0.+W.^2))));
    end

    Gr = Mat

    nt = length(y); N =size(W)[1];
    G = zeros(N,);
    G =  -1/noisesigma^2*Mat'*(r-y)
    # for i = 1:nt
    #     q =(-((r[i]-y[i])./noisesigma.^2))*Gr[i,:]';
    #     G = G + dropdims(q, dims = tuple(findall(size(q) .== 1)...))
    # end

   G = G  -2*W./(scale^2 .+ W.^2);
   return logp,G

end

function  logpdiff(U,args,cache)
    M = args.Mat;
    noisesigma = args.noisesigma
    scale = args.scale; y  = args.y; F = args.F

    Fxprop = cache.Fxprop
    Dxprop = cache.Dxprop
    res = cache.residual
    Fxprop .= F*U;
    Dxprop .= M*U;

    res .= Fxprop-y

    logp = -0.5/noisesigma^2*dotitself(res);
    logp = logp + sum(log.(scale./((scale^2 .+ Dxprop.^2))));

    return logp

end

function  logpdiffgradi(U,arguments,cache;both=true)
    M = arguments.Mat;
    noisesigma = arguments.noisesigma
    scale = arguments.scale; y  = arguments.y; F = arguments.F

    Fxprop = cache.Fxprop
    Dxprop = cache.Dxprop
    res = cache.residual
    Fxprop .= F*U;
    Dxprop .= M*U;

    res .= Fxprop-y
    logp = -0.5/noisesigma^2*dotitself(res);
    logp = logp + sum(log.(scale./((scale^2 .+ Dxprop.^2))));

    #nt = length(y); N = length(U);
    #for i = 1:nt
            #println(Gr[i,:]')
         #q =  (-((r[i]-y[i])./noisesigma.^2))*Gr[i,:]';
         #G = G + dropdims(q, dims = tuple(findall(size(q) .== 1)...))
     #end
     G= F'*(-((res)./noisesigma.^2))

    Gd =   M'*(-2.0*Dxprop./(scale^2 .+ Dxprop.^2));
    #Gd = Gd'*M;
    G = G + Gd;

    return logp,G

end

function  logpdiffpartial(ix,U,args,cache,current)
    D = @view args.Mat[:,ix];
    noisesigma = args.noisesigma
    scale = args.scale; y  = args.y;
    F  = @view args.F[:,ix]

    logp = 0.0

    currx = @view current.xcurr[ix]
    Fxcurr = current.Fxcurr
    Dxcurr =  current.Dxcurr

    xprop = cache.xprop
    xprop[ix] = U

    compx = U-currx
    Fxprop = cache.Fxprop
    Fxprop .= Fxcurr
    Fxprop .+= F*compx

    Dxprop = cache.Dxprop
    Dxprop .= Dxcurr
    Dxprop .+= D*compx

    res = cache.residual
    res .= Fxprop
    res .-=  y

    logp = -0.5/noisesigma^2*dotitself(res);
    logp = logp + sum(log.(scale./((scale^2 .+ Dxprop.^2))));

    return logp

end

function  logpdiffpartialgradi(ix,U,arguments,cache,current;both=true)

    D = @view arguments.Mat[:,ix];
    noisesigma = arguments.noisesigma
    scale = arguments.scale; y  = arguments.y; F = @view arguments.F[:,ix]

    #(Fxcurr=(F*X0),Dxcurr=(Md1*X0), xcurr = copy(X0))

    currx = @view current.xcurr[ix]
    Fxcurr = current.Fxcurr
    Dxcurr =  current.Dxcurr

    xprop = cache.xprop
    xprop[ix] = U

    compx = U-currx
    Fxprop = cache.Fxprop
    Fxprop .= Fxcurr
    Fxprop .+= F*compx

    Dxprop = cache.Dxprop
    Dxprop .= Dxcurr
    Dxprop .+= D*compx

    res = cache.residual
    res .= Fxprop
    res .-=  y

    logp = 0.0

    if (both)
        logp = -0.5/noisesigma^2*dotitself(res);
        logp = logp + sum(log.(scale./((scale^2 .+ Dxprop.^2))));
    end

    G = cache.gradiprop

    G[ix] = -F'*res/noisesigma^2

    G[ix] = G[ix] + D'*(-2.0*Dxprop./(scale^2 .+ Dxprop.^2));

    return logp,G[ix]

end

function refreshcurrent!(cache,current;mode=1)
    if (mode ==1)
        current.Mxcurr .= cache.Mxprop
        current.noisecurr .= cache.noiseprop
    elseif (mode==2)
        current.xcurr .= cache.xprop
        current.Fxcurr .= cache.Fxprop
        current.Dxcurr .= cache.Dxprop
    end
end

function repatpartial(lpdf,placeholder,ix,x0,args,cache,current,tuning;mode=1,Np=1,Nruns=1,T=-Inf,extra=copy(x0))::Vector{Float64}
    x = copy(x0[ix])
    dim = length(x)

    logpdf(xp) = lpdf(ix,xp,args,cache,current)
    #logpdf(xp) = lpdf(partialreplace(x0,xp,ix),args,cache)
    density = logpdf(x)

    z = copy(extra[ix])
    u = randn(dim)
    densityz = logpdf(z)

    for i = 1:Nruns
        R = tuning.Cho.L/sqrt(2)
		u .= randn(dim)
		xtilde = x + R * u
		densitytilde = logpdf(xtilde)
		while (log(rand()) >=  density - densitytilde)
			u .= randn(dim)
			xtilde .= x + R * u
			densitytilde = logpdf(xtilde)
		end

		u .= randn(dim)
		xstar = xtilde + R * u
		densitystar = logpdf(xstar)

		while (log(rand()) >=  densitystar - densitytilde)
			u .= randn(dim)
			xstar .= x + R * u
			densitystar = logpdf(xstar)
		end

		u .= randn(dim)
		zstar = xstar + R * u
		densityzstar = logpdf(zstar)

		while (log(rand()) >=  densitystar - densityzstar)
			u .= randn(dim)
			zstar .= xstar + R * u
			densityzstar = logpdf(zstar)
		end

        ratio = densitystar + min(0,density-densityz) - density - min(0,densitystar-densityzstar)
        if (log(rand()) <= ratio)
            x .= xstar
            z .= zstar
			extra[ix] = zstar
			density = densitystar
			densityz = densityzstar

        end
        tuning.n = tuning.n + 1

        if (tuning.n <= tuning.Nadapt)
            xmp = copy(tuning.xm)
            tuning.xm .= tuning.xm + (x-tuning.xm)/tuning.n
            xm = tuning.xm
            tuning.C .= ((tuning.n-1)*tuning.C + (x-xm)*(x-xmp)')/tuning.n
            tuning.Cho = cholesky(Hermitian(2.38^2/dim*tuning.C + 1e-11*I(dim)))
        end

    end
    logpdf(x)
    refreshcurrent!(cache,current;mode=mode)
    return x
end

function mwgpartial(lpdf,gradilpdf,ix,x00,args,cache,current,tuning;Nruns=5,Np=3,T=5.0,mode=1,extra=nothing)
    x0 = repeat(x00[ix],Np)
    x = copy(x0)
    Nall = length(x)

    xnew = similar(x)
    pnew = similar(xnew)
    p = similar(xnew)

    logpdf_can(xp) = lpdf(ix,xp,args,cache,current)
    #logpdf_can(xp) = lpdf(partialreplace(x00,xp,ix),args,cache)

    logpdf(x) = totallogpdf(x,logpdf_can,Np,T)

    glpdf_can(xp) = gradilpdf(ix,xp,args,cache,current;both=false)[2]
    #glpdf_can(xp) = gradilpdf(partialreplace(x00,xp,ix),args,cache)[2][ix]

    glpdf(x) = totalloggrad(x,logpdf_can,glpdf_can,Np,T)[2]
    glpdfboth(x) = totalloggrad(x,logpdf_can,glpdf_can,Np,T)

    ll = logpdf(x)

    for _ = 1:Nruns

        tuning.n = tuning.n + 1

        V=randn(Nall)
        xnew .= x+ tuning.Cho.L*V

        lnew = logpdf(xnew)

        U = rand()
        apr = min(1.0,exp(-ll + lnew))

        if (U < apr)
            ll = lnew
            x .= xnew
            tuning.acc = tuning.acc +1
        end

        if (tuning.n <= tuning.Nadapt)
            xmp = copy(tuning.xm)
            tuning.xm .= tuning.xm + (x-tuning.xm)/tuning.n
            xm = tuning.xm
            tuning.C .= ((tuning.n-1)*tuning.C + (x-xm)*(x-xmp)')/tuning.n
            tuning.Cho = cholesky(Hermitian(2.38^2/Nall*tuning.C + 1e-11*I(Nall)))

            # atarg = 0.3
            # V = V/norm(V)
            # z = sqrt(tuning.n^(-0.7) *abs(tuning.acc/tuning.n - atarg)) *  tuning.Cho.L*V
            # #println(z)
            # if (tuning.acc/tuning.n >= atarg)
            #     lowrankupdate!(tuning.Cho, z)

            # else
            #     lowrankdowndate!(tuning.Cho, z)
            # end

        end

    end

    xret = reweight(x,Np,logpdf_can,T)
    logpdf_can(xret)
    refreshcurrent!(cache,current;mode=mode)

    return xret

end

mutable struct rram
    n::Int64
    acc::Int64
    Nadapt::Int64
    C::Array{Float64,2}
    Cho::Cholesky{Float64,Array{Float64,2}}
    xm::Vector{Float64}
end

function totalsample(lpdf,glpdf,x0,meas,arg,ix;N=1,Nadapt=1000,Nruns=5,thinning=10,algo=repatpartial,mode=1,Np=2,T=2.0)
    x = copy(x0)
    Npar = length(x0)
    chain = zeros(Npar, Int(floor(N / thinning)))
    Nix = length(ix)

    if (mode==2)
        cache=(xprop=copy(x0),Fxprop=(arg.F*x0),Dxprop=(arg.Mat*x0), residual=similar(meas),gradiprop=similar(x0))
        current = (Fxcurr=(arg.F*x0),Dxcurr=(arg.Mat*x0), xcurr = copy(x0))
    else
        cache = []
        current = []
    end

    if algo == repatpartial
        extra = copy(x0)
        tuning = Vector{rram}(undef,Nix)
        for k = 1:Nix
            n = length(ix[k])
            C = 0.1*Array(I(n))
            Cho = cholesky(C)
            xm = copy(x0[ix[k]])
            tuning[k] = rram(1,0,Nadapt,C,Cho,xm)
        end

    elseif (algo==mwgpartial)
        extra = []
        tuning = Vector{rram}(undef,Nix)
        #tuning = Vector{dahmc2}(undef,Nix)
        #tuning = Vector{ahmc}(undef,Nix)
        for k = 1:Nix

            n = length(ix[k])*Np
            C = 5*Array(I(n))
            Cho = cholesky(C)
            xm = repeat(copy(x0[ix[k]]),Np)
            tuning[k] = rram(1,0,Nadapt,C,Cho,xm)

        end

    pb = ProgressBar(1:N)

    for i in pb
        for k = 1:Nix
            x[ix[k]] = algo(lpdf,glpdf,ix[k],x,arg,cache,current,tuning[k];mode=mode,Nruns=Nruns,extra=extra,Np=Np,T=T)#lpdf,ix,x0,args,cache,current,tuning
            if (i % thinning == 0)
                chain[ix[k], Int(i / thinning)] = x[ix[k]]
            end

        end
    end

    return chain,tuning
end

function savechain(c,prefix)
    p = pwd()
    f = p*"/chains/"*prefix*"_"*Dates.format(Dates.now(), "yyyy_mm_dd_HH_MM_SS")*".mat"
    matwrite(f,Dict("chain"=>c))

end

N = 200;
x = Vector(range(-0.0,stop=1.0,length=N));    dx = x[2]-x[1];
xi =  copy(x)
measxi = copy(xi[1:3:end])

triangle(x) = 1.0*(x+1)*(x<=0)*(x>=-1) + 1.0*(-x+1)*(x>0)*(x<=1);
heavi(x) = 1.0*(x>=0);
tf(x) = 0+ 0.2*heavi(x-0.1)*heavi(-x+0.15) - 0.5*heavi(x-0.65)*heavi(-x+0.7)+ 0*exp(-60*abs(x-0.2)) - 0*exp(-180*abs(x-0.8)) + 0*triangle(10*(x-0.8)) + 1*heavi(x-0.3)*heavi(-x+0.6);
tf(x) = heavi(x-0.75)*heavi(-x+0.9)+ triangle(10*(x-0.15)) + triangle(10*(x-0.55))*heavi(x-0.55) + 0*triangle(30*(x-0.65))*heavi(-x+0.65) + exp(-70*abs(x-0.4)) # heavi(x-0.45)*heavi(-x+0.55)
sigma = 0.1;  # Standard deviation of the measurement noise.
cw = 500;
kernel(x,y) = sqrt(cw)/sqrt(π)*exp(-cw*(x-y)^2)
tfc(x) = quadgk(y ->  tf(y)*kernel(x,y), 0, 1, rtol=1e-9)[1]
tfc(x) = HCubature.hquadrature(y ->  tf(y)*kernel(x,y), 0.0, 1.0; rtol=1e-6,atol=1e-6,initdiv=7)[1]

F = measurementmatrix(xi,measxi,kernel); # Theory matrix.
gt = tf.(xi);
gtc =  tfc.(measxi); #F*gt
gtc2 = F*gt

Random.seed!(100)

meas = gtc .+ sigma.*randn(size(measxi)); # Simulate the measurements.

lcauchy = 0.015;
acauchy = lcauchy^2; bcauchy = 1;
cauchypriorscale = 0.5;

cd1scale = 0.05;
cd2scale = 0.01;

Md1=difference1(xi);

Md2=difference2(xi);

Mcauchy = smatrix(x,acauchy,bcauchy); MN = size(Mcauchy)[1];

measnoisescaling = 1.0
diff1arg = (F=F, y=meas,Mat=sparse(Md1),scale=cd1scale,noisesigma=measnoisescaling*sigma)
diff2arg = (F=F, y=meas,Mat=sparse(Md2),scale=cd2scale,noisesigma=measnoisescaling*sigma)
varg = (Mat=Mcauchy,F=F,y=meas,scale=cauchypriorscale,noisesigma=measnoisescaling*sigma)

cached=(xprop=similar(xi),Fxprop=(F*xi),Dxprop=(Md1*xi), residual=similar(meas),gradiprop=similar(xi))
cache2d=(xprop=similar(xi),Fxprop=(F*xi),Dxprop=(Md2*xi), residual=similar(meas),gradiprop=similar(xi))
cachev =(xprop = similar(xi), Fxprop = F*xi, Dxprop = Mcauchy*xi,residual=similar(meas),gradiprop=similar(xi))

targetdiff1(w) =-logpdiff(w,diff1arg,cached)
targetdiffgrad1(w) = -logpdiffgradi(w,diff1arg,cached;both=false)[2]#-logpcauchygradi(w,cauchyarg)[2]#
res1 = Optim.optimize(targetdiff1, targetdiffgrad1,0*randn(N,),Optim.LBFGS(), Optim.Options(allow_f_increases=true,show_trace=false,iterations=1000); inplace = false)
MAP_diff1 = res1.minimizer

targetdiff2(w) =-logpdiff(w,diff2arg,cache2d)
targetdiffgrad2(w) = -logpdiffgradi(w,diff2arg,cached;both=false)[2]#-logpcauchygradi(w,cauchyarg)[2]#
res2 = Optim.optimize(targetdiff2, targetdiffgrad2,0*randn(N,),Optim.LBFGS(), Optim.Options(allow_f_increases=true,show_trace=false,iterations=30000); inplace = false)
MAP_diff2 = res2.minimizer

targetv(w) = -logpdiff(w,varg,cachev)
targetvgrad(w) = -logpdiffgradi(w,varg,cachev;both=false)[2]#-logpcauchygradi(w,cauchyarg)[2]#
res2 = Optim.optimize(targetv, targetvgrad,0.00*randn(N,),Optim.LBFGS(), Optim.Options(allow_f_increases=true,show_trace=false,iterations=30000); inplace = false)
MAP_v = res2.minimizer

Ntot = length(xi)

function NUTS_sample(N,arg,cac)
    n_samples, n_adapts = 100, 30
    ss = Vector{Any}(undef,N)
    st = Vector{Any}(undef,N)
    Np = 1
    T = 3.0
    for i = 1:N

        Ntot = length(xi)*Np
        Random.seed!(i)
        initial_θ = repeat(MAP_v,Np)
        global c

        targetdiff(w) = logpdiff(w,arg,cac)
        targetdiffgrad(w) = logpdiffgradi(w,arg,cac;both=true)

        metric = DiagEuclideanMetric(Ntot)# DenseEuclideanMetric(Ntot)# DiagEuclideanMetric(Ntot)
        hamiltonian = Hamiltonian(metric, targetdiff, targetdiffgrad)

        initial_ϵ = find_good_stepsize(hamiltonian, initial_θ)
        integrator =  Leapfrog(initial_ϵ)#  TemperedLeapfrog(initial_ϵ, 1.05)# Leapfrog(initial_ϵ)
        proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator,10,1000.0)
        adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))
        c, stats = sample(hamiltonian, proposal, initial_θ, n_samples, adaptor, n_adapts; progress=true)
        ss[i] = hcat(c...)
        # ss[i] = c
        st[i] = stats
        pr="NUTS_"*string("diff1_NUTS_")
        savechain(c,pr)

    end
    return ss, st
end

#c,b = NUTS_sample(10,diff1arg,cached)

ix = divide(N,1)
algo =  repatpartial#mwgpartial
grfunk = logpdiffpartialgradi#  logpcauchyreppartialgradi#
funk =  logpdiffpartial#logpcauchyreppartial#logpdiffpartial#
par = diff1arg
mode = 2# 1
x000  = copy(MAP_diff1) #randn(N,)

cached=(xprop=copy(x000),Fxprop=(F*x000),Dxprop=(Md1*x000), residual=copy(meas),gradiprop=copy(gt))
currentd = (Fxcurr=(F*MAP_diff1),Dxcurr=(Md1*MAP_diff1), xcurr = copy(MAP_diff1))

for t=1:0
    Random.seed!(t)
    xz = copy(MAP_v)#0.1*randn(N,)

    Nsamples = 1500; Nadapt = 1250; thinning = 20
    Np = 1; Nruns = 5; T = 1.0

    c,tun=totalsample(funk,grfunk,xz,meas,par,ix;N=Nsamples,Nadapt=Nadapt,thinning=thinning,algo=algo,mode=mode,Np=Np,T=T,Nruns=Nruns)
    savechain(c,"diff1_"*string(nameof(algo))*"_"*string(nameof(funk))*"_N"*string(Nsamples)*"_Nadapt"*string(Nadapt)*"_thin"*string(thinning)*"_Np"*string(Np)*"_Nruns"*string(Nruns)*"_T"*string(T))
    println(effective_sample_size(c[1,:]))
    global c
    global tun
end

max_stars_repo_name: suurj/markovcauchy

max_issues_repo_name: suurj/markovcauchy

max_forks_repo_name: suurj/markovcauchy


text:
module VarStructs

if VERSION <= v"1.1"
    isnothing(x) = x == nothing
end
if VERSION <= v"1.2.0"
    @warn("Please upgrade to Julia 1.4 for a better experience")
    Base.print(io, x::Nothing) = Base.show(io, x)
end

export @var, @shared_var, Unset, isunset

include("sugar.jl")

"""
    Unset

The value of a VarStruct field when it is not set yet.

# Examples
```julia
mystruct = @var struct MyStruct
    x::Int64
end

typeof(mystruct.x) # Unset

```
"""
struct Unset end

"""
  isunset(x)

Returns true if the value of a field is not set yet.

# Examples
```julia
mystruct = @var struct MyStruct
    x::Int64
end

isunset(mystruct.x) # true
```

```julia
isunset( Unset() ) # true
```
"""
isunset(x) = isa(x, Unset)

"""
type of the field names

# Example
```julia
f = :Amin
typeof(f) <: VarStructs.Name # true
```
"""
Name = Union{Symbol, String}

"""
A struct to hold the information about a VarStruct field

# Example
```julia
fp = Props(20.0, Float64)
```
"""
mutable struct Props{T}
    value::Union{Unset, T} # Unset when the value is not provided yet, but the type is known
    type::Type{T}
    # check::Union{Function, Nothing}
end

function Props(value::Union{Unset, T}) where T
    return Props{T}(value, T)
end

"""
A Dict that maps the field names to field properties

# Examples
```julia
ft = FieldTable( :Amin => Props(20.0, Float64) )
```
"""
FieldTable = Dict{Name, Props}

################################################################
include("parse.jl")
include("create.jl")

################################################################
"""
Create shared var struct instances (calling constructor will mutate the old instances).
"""
macro shared_var(expr)
    # __module__ = @__MODULE__ # for functions debuging.
    expr = macroexpand(__module__, expr) # to expand literal macros and @static
    return var(expr, __module__, __source__,  true)
end

macro var(expr)
    # __module__ = @__MODULE__ # for functions debuging.
    expr = macroexpand(__module__, expr) # to expand literal macros and @static
    return var(expr, __module__, __source__,  false)
end

function var(expr, __module__, __source__,  shred_var::Bool = false)
    #  check if @var is used before struct
    if expr isa Expr && expr.head == :struct

        # expr.args[3] # arguments
        expr.args[3], args_field, args_defaultvalue, args_type, args_param, args_check, is_struct_mutable, T = varstruct_struct_parse(__module__, expr)

    #  check if @var  is used before call
    elseif expr isa Expr && expr.head == :call

        expr, args_field, args_defaultvalue, args_type, args_param, T = varstruct_call_parse(__module__, expr)

    else
        error("Invalid usage of @var")
    end

    # check if the field name is a defined type
    for arg_field in args_field
        if isdefined(__module__, arg_field) && isa(getfield(__module__, arg_field), Type)
            @error "Change the field name `$arg_field` in struct `$T` to something else. The name conflicts with an already defined type name." _module = __module__  _file =__source__.file  _line = __source__.line
        end
    end

    return varstruct_construct(__module__, T, args_field, args_defaultvalue, args_type, shred_var)
end

end


text:
export visualize_displacement, visualize_von_mises_stress, visualize_mesh, visualize_boundary,
visualize_scalar_on_scoped_body, visualize_total_deformation_on_scoped_body, visualize_von_mises_stress_on_scoped_body,
visualize_x_deformation_on_scoped_body, visualize_y_deformation_on_scoped_body, visualize_scalar_on_undeformed_body

"""
    visualize_von_mises_stress(domain::Domain; frames::Int64 = 20, kwargs...)

Animation of von Mises stress tensors.

![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/NNFEM/visualize_von_mises_stress.gif?raw=true)
"""
function visualize_von_mises_stress(domain::Domain; frames::Int64 = 20, kwargs...)
    NT = length(domain.history["stress"])
    if NT == 0
        error(ArgumentError("history[\"stress\"] is empty."))
    end
    visualize_von_mises_stress_on_scoped_body(zeros(NT+1, domain.nnodes*2), domain; frames = frames, kwargs...)
end

"""
    visualize_von_mises_stress(domain::Domain, t_step::Int64; kwargs...)

Plot of von Mises stress tensors at time step `t_step`.

![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/NNFEM/visualize_von_mises_stress_50.png?raw=true)
"""
function visualize_von_mises_stress(domain::Domain, t_step::Int64; kwargs...)
    domain2 = copy(domain)
    domain2.history["stress"] = [domain.history["stress"][t_step]]
    visualize_von_mises_stress_on_scoped_body(zeros(1, domain2.nnodes*2), domain2; frames=1, kwargs...)
end

"""
    visualize_displacement(domain::Domain)
    visualize_displacement(u::Array{Float64, 2}, domain::Domain)

Animation of displacements using points.

![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/NNFEM/visualize_displacement.gif?raw=true)
"""
function visualize_displacement(domain::Domain; scale_factor::Float64 = 1.0)
    u = hcat(domain.history["state"]...) * scale_factor
    visualize_displacement(u, domain)
end

function visualize_displacement(u::Array{Float64, 2}, domain::Domain)
    if size(u,2)==2domain.nnodes
        u = Array(u')
    end
    X0, Y0 = domain.nodes[1:domain.nnodes], domain.nodes[domain.nnodes+1:end]
    NT = size(u, 2)
    U0 = zeros(domain.nnodes, NT)
    V0 = zeros(domain.nnodes, NT)
    for i = 1:NT
        U0[:,i] = X0 + u[1:domain.nnodes, i]
        V0[:,i] = Y0 + u[domain.nnodes+1:end, i]
    end
    xmin, xmax = minimum(U0), maximum(U0)
    ymin, ymax = minimum(V0), maximum(V0)
    l1 = (xmax-xmin)/10
    l2 = (ymax-ymin)/10
    xmin -= l1
    xmax += l1
    ymin -= l2
    ymax += l2
    close("all")
    p, = plot([], [], ".", markersize=5)
    t = title("Snapshot = 0")
    xlim(xmin, xmax)
    ylim(ymin, ymax)
    axis("scaled")
    function update(frame)
        p.set_data(U0[:,frame], V0[:,frame])
        t.set_text("Snapshot = $frame")
    end
    out = animate(update, Int64.(round.(LinRange(1, NT, 20))))
    xlabel("x")
    ylabel("y")
    gca().invert_yaxis()
    out
end

"""
    visualize_displacement(domain::Domain)
    visualize_displacement(nodes::Array{Float64,2}, elems::Array{Int64, 2})

Visualizes the mesh.
![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/NNFEM/visualize_mesh.png?raw=true)
"""
function visualize_mesh(nodes::Array{Float64,2}, elems::Array{Int64, 2})
    patches = PyObject[]
    for i = 1:size(elems,1)
        e = elems[i,:]
        p = plt.Polygon(nodes[e,:],edgecolor="k",lw=1,fc=nothing,fill=false)
        push!(patches, p)
    end
    p = matplotlib.collections.PatchCollection(patches, match_original=true)
    gca().add_collection(p)
    axis("scaled")
    xlabel("x")
    ylabel("y")
    gca().invert_yaxis()
end

function visualize_mesh(domain::Domain)
    elements = zeros(Int64, length(domain.elements), 4)
    for (k,e) in enumerate(domain.elements)
        elements[k,:] = e.elnodes
    end
    visualize_mesh(domain.nodes, elements)
end

"""
    visualize_boundary(domain::Domain, direction::String="x")

Visualizes the boundary conditions. The boundary configuration is shown in the direction `direction`.

![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/NNFEM/visualize_boundary.png?raw=true)
"""
function visualize_boundary(domain::Domain, direction::String="x")
    visualize_mesh(domain)
    direction = direction == "x" ? 1 : 2;

    function _helper(idx, marker, Label)
        if length(idx)==0
            return
        end
        x, y = domain.nodes[idx,1], domain.nodes[idx,2]
        plot(x, y, marker, markersize=10, label="$Label")
    end
    _helper(findall(domain.EBC[:,direction] .== -1), ".", "Fixed Dirichlet")
    _helper(findall(domain.EBC[:,direction] .== -2), ".", "Time-dependent Dirichlet")
    _helper(findall(domain.FBC[:,direction] .== -1), "+", "Fixed Neumann")
    _helper(findall(domain.FBC[:,direction] .== -2), "+", "Time-dependent Neumann")

    if size(domain.edge_traction_data,1)>0
        ids = unique(domain.edge_traction_data[:,3])
        for (k, ids_) in enumerate(ids)
            data = domain.edge_traction_data[domain.edge_traction_data[:,3].==ids_,1:2]

            for i = 1:size(data, 1)
                elem = domain.elements[data[i,1]]
                if data[i,2]==1
                    idx = [1;2]
                elseif data[i,2]==2
                    idx = [2;3]
                elseif data[i,2]==3
                    idx = [3;4]
                elseif data[i,2]==4
                    idx = [4;1]
                end
                x = elem.coords[idx,1]
                y = elem.coords[idx,2]
                if i==1
                    plot(x, y, "--",linewidth=2, color="C$(k+3)", label="Edge Traction (ID=$ids_)")
                else
                    plot(x, y, "--",linewidth=2, color="C$(k+3)")
                end
            end
        end

    end
    legend()
    # gca().invert_yaxis()

end

"""
    visualize_scalar_on_scoped_body(d::Array{Float64}, domain::Domain, kwargs...)
    visualize_scalar_on_scoped_body(s::Array{Float64, 1}, d::Array{Float64,1}, domain::Domain;
        scale_factor::Float64 = -1.0, kwargs...)

Plot the scalar on scoped body. For example, `s` can be the von Mises stress tensor.
"""
function visualize_scalar_on_scoped_body(s::Array{Float64, 1}, d::Array{Float64,1}, domain::Domain;
        scale_factor::Real = -1.0, kwargs...)
    local sv
    if scale_factor<0.0

        scale_factor = min(
            0.1 * (maximum(domain.nodes[:,1]) - minimum(domain.nodes[:,1]))/maximum(abs.(d[1:domain.nnodes])),
            0.1 * (maximum(domain.nodes[:,2]) - minimum(domain.nodes[:,2]))/maximum(abs.(d[domain.nnodes+1:end]))
        )
        if isinf(scale_factor)
            scale_factor = 1.0
        end
        @info "scale_factor automatically set to $scale_factor"
    end
    vmin = haskey(kwargs, :vmin) ? kwargs[:vmin] : minimum(s)
    vmax = haskey(kwargs, :vmax) ? kwargs[:vmax] : maximum(s)
    n = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
    cm = matplotlib.cm
    cmap = cm.jet
    m = cm.ScalarMappable(norm=n, cmap=cmap)
    elements = domain.elements
    patches = []
    sv_id = 0
    for i = 1:size(elements,1)
        e = elements[i].elnodes
        if length(s)==domain.nnodes
            sv = mean(s[e])
        elseif length(s)==length(elements)
            sv = s[i]
        elseif length(s)==getNGauss(domain)
            ngpt = length(elements[i].weights)
            sv = mean(s[sv_id+1:sv_id+ngpt])
            sv_id = sv_id + ngpt
        else
            error(ArgumentError("Dimension of the scalar function must be #nodes, #elements, or #guass_points."))
        end
        dv = [d[e] d[e .+ domain.nnodes]]
        p = plt.Polygon(elements[i].coords + scale_factor * dv,edgecolor="k",lw=1,fill=true,
            fc=m.to_rgba(sv))
        push!(patches, p)
    end
    p = matplotlib.collections.PatchCollection(patches, match_original=true)
    gca().add_collection(p)
    colorbar(m)
    gca().invert_yaxis()
    axis("scaled")
    xlabel("x")
    ylabel("y")
end

function visualize_scalar_on_scoped_body(s_all::Array{Float64, 2}, d_all::Array{Float64,2}, domain::Domain;
    scale_factor::Real = -1.0, frames::Int64 = 20, kwargs...)
    local sv
    if scale_factor<0.0
        scale_factor = min(
            0.1 * (maximum(domain.nodes[:,1]) - minimum(domain.nodes[:,1]))/maximum(abs.(d_all[:, 1:domain.nnodes])),
            0.1 * (maximum(domain.nodes[:,2]) - minimum(domain.nodes[:,2]))/maximum(abs.(d_all[:, domain.nnodes+1:end]))
        )
        if isinf(scale_factor)
            scale_factor = 1.0
        end
        @info "scale_factor automatically set to $scale_factor"
    end

    if frames==1
         visualize_scalar_on_scoped_body(s_all[1,:], d_all[1,:], domain; scale_factor = scale_factor, kwargs...)
         return nothing
    end
    vmin = haskey(kwargs, :vmin) ? kwargs[:vmin] : minimum(s_all)
    vmax = haskey(kwargs, :vmax) ? kwargs[:vmax] : maximum(s_all)
    n = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
    n = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
    cm = matplotlib.cm
    cmap = cm.jet
    m = cm.ScalarMappable(norm=n, cmap=cmap)
    elements = domain.elements
    close("all")
    # Find the bounding box: the domain takes up 80%
    D = zeros(size(d_all,1), domain.nnodes)
    for i = 1:size(d_all,1)
        D[i,:] = domain.nodes[:,1] + scale_factor*d_all[i,1:domain.nnodes]
    end
    a1, b1 = minimum(D), maximum(D)
    h1 = (b1-a1)*0.1

    D = zeros(size(d_all,1), domain.nnodes)
    for i = 1:size(d_all,1)
        D[i,:] = domain.nodes[:,2] + scale_factor*d_all[i,domain.nnodes+1:end]
    end
    a2, b2 = minimum(D), maximum(D)
    h2 = (b2-a2)*0.1

    function update(frame)
        s = s_all[frame,:]
        d = d_all[frame,:]
        gca().clear()
        patches = []
        sv_id = 0
        for i = 1:size(elements,1)
            e = elements[i].elnodes
            if length(s)==domain.nnodes
                sv = mean(s[e])
            elseif length(s)==length(elements)
                sv = s[i]
            elseif length(s)==getNGauss(domain)
                ngpt = length(elements[i].weight)
                sv = mean(s[sv_id+1:sv_id+ngpt])
                sv_id = sv_id + ngpt
            else
                error(ArgumentError("Dimension of the scalar function must be #nodes, #elements, or #guass_points."))
            end
            dv = [d[e] d[e .+ domain.nnodes]]
            p = plt.Polygon(elements[i].coords + scale_factor * dv,edgecolor="k",lw=1,fill=true,
                fc=m.to_rgba(sv))
            push!(patches, p)
        end
        p = matplotlib.collections.PatchCollection(patches, match_original=true)
        gca().add_collection(p)
        xlabel("x")
        ylabel("y")
        title("Snapshot = $frame")
        axis("scaled")
        xlim(a1-h1, b1+h1)
        ylim(a2-h2, b2+h2)
        gca().invert_yaxis()
    end
    update(1)
    t = title("Snapshot = 1")
    cb = colorbar(m)
    animate(update, Int64.(round.(LinRange(div(size(s_all,1), frames), size(s_all,1),frames))))
end

"""
    visualize_scalar_on_undeformed_body(s::Array{Float64, 1}, domain::Domain; kwargs...)
    visualize_scalar_on_undeformed_body(s::Array{Float64, 2}, domain::Domain; frames::Int64 = 20, kwargs...)

Plots or animates scalar values `s` on the domain `domain`

"""
function visualize_scalar_on_undeformed_body(s::Array{Float64, 1}, domain::Domain; kwargs...)
    visualize_scalar_on_scoped_body(s, zeros(2*domain.nnodes), domain)
end

function visualize_scalar_on_undeformed_body(s::Array{Float64, 2}, domain::Domain; frames::Int64 = 20, kwargs...)
    visualize_scalar_on_scoped_body(s, zeros(2*domain.nnodes), domain; frames = frames, kwargs...)
end

@doc raw"""
    visualize_total_deformation_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
    scale_factor::Float64 = 1.0, frames::Int64 = 20)

Visualizes the total deformation

$$\sqrt{u_x^2 + u_y^2}$$

![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/NNFEM/visualize_total_deformation_on_scoped_body.gif?raw=true)
"""
function visualize_total_deformation_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
    scale_factor::Real = 1.0, frames::Int64 = 20, kwargs...)
    s_all = d_all
    S = zeros(size(s_all,1), domain.nnodes)
    for i = 1:size(s_all,1)
        S[i,:] = @. sqrt( s_all[i,1:domain.nnodes]^2 + s_all[i,domain.nnodes+1:end]^2 )
    end
    visualize_scalar_on_scoped_body(S, d_all, domain; scale_factor = scale_factor, frames=frames, kwargs...)
end

@doc raw"""
    visualize_x_deformation_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
    scale_factor::Real = 1.0, frames::Int64 = 20, kwargs...)

Visualizes the $x$ directional displacement.
"""
function visualize_x_deformation_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
    scale_factor::Real = 1.0, frames::Int64 = 20, kwargs...)
    visualize_scalar_on_scoped_body(d_all[:, 1:domain.nnodes], d_all, domain; scale_factor = scale_factor, frames=frames, kwargs...)
end

function visualize_x_deformation_on_scoped_body(d_all::Array{Float64,1}, domain::Domain;
    scale_factor::Real = 1.0, frames::Int64 = 20, kwargs...)
    visualize_scalar_on_scoped_body(d_all[1:domain.nnodes], d_all, domain; scale_factor = scale_factor, kwargs...)
end

@doc raw"""
    visualize_y_deformation_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
    scale_factor::Real = 1.0, frames::Int64 = 20, kwargs...)

Visualizes the $y$ directional displacement.
"""
function visualize_y_deformation_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
    scale_factor::Real = 1.0, frames::Int64 = 20, kwargs...)
    visualize_scalar_on_scoped_body(d_all[:, domain.nnodes+1:end], d_all, domain; scale_factor = scale_factor, frames=frames, kwargs...)
end

"""
    visualize_von_mises_stress_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
    scale_factor::Real = -1.0, frames::Int64 = 20)

Similar to [`visualize_von_mises_stress`](@ref), but the domain can be a deformed body.

![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/NNFEM/visualize_von_mises_stress_on_scoped_body.gif?raw=true)
"""
function visualize_von_mises_stress_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
    scale_factor::Real = -1.0, frames::Int64 = 20, kwargs...)
    stress = domain.history["stress"]
    S = zeros(length(stress), length(domain.elements))
    x = zeros(length(domain.elements))
    y = zeros(length(domain.elements))
    for t = 1:length(stress)
        cnt = 1
        for (k,e) in enumerate(domain.elements)
            ct = mean(domain.elements[k].coords, dims=1)
            x[k], y[k] = ct[1,1], ct[1,2]

                ss = Float64[]
                nstress = length(e.mat)
                for p = 1:nstress
                    push!(ss, postprocess_stress(stress[t][cnt, :] ,"vonMises"))
                    cnt += 1
                end
                S[t, k] = mean(ss)
        end
    end

    visualize_scalar_on_scoped_body(S, d_all, domain; scale_factor = scale_factor, frames=frames, kwargs...)
end


text:
module CSV

# stdlib
using Mmap, Dates, Unicode
using Parsers, Tables
using PooledArrays, CategoricalArrays, WeakRefStrings, DataFrames

function validate(fullpath::Union{AbstractString,IO}; kwargs...)
    Base.depwarn("`CSV.validate` is deprecated. `CSV.read` now prints warnings on misshapen files.", :validate)
    Tables.columns(File(fullpath; kwargs...))
    return
end

include("utils.jl")
include("filedetection.jl")

struct Error <: Exception
    msg::String
end

Base.showerror(io::IO, e::Error) = println(io, e.msg)

struct File
    name::String
    names::Vector{Symbol}
    types::Vector{Type}
    rows::Int64
    cols::Int64
    e::UInt8
    categorical::Bool
    refs::Vector{Vector{String}}
    buf::Vector{UInt8}
    tapes::Vector{Vector{UInt64}}
end

getname(f::File) = getfield(f, :name)
getnames(f::File) = getfield(f, :names)
gettypes(f::File) = getfield(f, :types)
getrows(f::File) = getfield(f, :rows)
getcols(f::File) = getfield(f, :cols)
gete(f::File) = getfield(f, :e)
getcategorical(f::File) = getfield(f, :categorical)
function getrefs(f::File, col)
    @inbounds r = getfield(f, :refs)[col]
    return r
end
getbuf(f::File) = getfield(f, :buf)
function gettape(f::File, col)
    @inbounds t = getfield(f, :tapes)[col]
    return t
end

function Base.show(io::IO, f::File)
    println(io, "CSV.File(\"$(getname(f))\"):")
    println(io, "Size: $(getrows(f)) x $(getcols(f))")
    show(io, Tables.schema(f))
end

const EMPTY_POSITIONS = Int64[]
const EMPTY_TYPEMAP = Dict{TypeCode, TypeCode}()
const EMPTY_REFS = Vector{String}[]
const EMPTY_REFVALUES = String[]

"""
    CSV.File(source; kwargs...) => CSV.File

Read a csv input (a filename given as a String or FilePaths.jl type, or any other IO source), returning a `CSV.File` object.

Opens the file and uses passed arguments to detect the number of columns and column types.
The returned `CSV.File` object supports the [Tables.jl](https://github.com/JuliaData/Tables.jl) interface
and can iterate `CSV.Row`s. `CSV.Row` supports `propertynames` and `getproperty` to access individual row values.
Note that duplicate column names will be detected and adjusted to ensure uniqueness (duplicate column name `a` will become `a_1`).
For example, one could iterate over a csv file with column names `a`, `b`, and `c` by doing:

```julia
for row in CSV.File(file)
    println("a=\$(row.a), b=\$(row.b), c=\$(row.c)")
end
```

By supporting the Tables.jl interface, a `CSV.File` can also be a table input to any other table sink function. Like:

```julia
# materialize a csv file as a DataFrame
df = CSV.File(file) |> DataFrame!

# load a csv file directly into an sqlite database table
db = SQLite.DB()
tbl = CSV.File(file) |> SQLite.load!(db, "sqlite_table")
```

Supported keyword arguments include:
* File layout options:
  * `header=1`: the `header` argument can be an `Int`, indicating the row to parse for column names; or a `Range`, indicating a span of rows to be concatenated together as column names; or an entire `Vector{Symbol}` or `Vector{String}` to use as column names; if a file doesn't have column names, either provide them as a `Vector`, or set `header=0` or `header=false` and column names will be auto-generated (`Column1`, `Column2`, etc.)
  * `normalizenames=false`: whether column names should be "normalized" into valid Julia identifier symbols; useful when iterating rows and accessing column values of a row via `getproperty` (e.g. `row.col1`)
  * `datarow`: an `Int` argument to specify the row where the data starts in the csv file; by default, the next row after the `header` row is used. If `header=0`, then the 1st row is assumed to be the start of data
  * `skipto::Int`: similar to `datarow`, specifies the number of rows to skip before starting to read data
  * `footerskip::Int`: number of rows at the end of a file to skip parsing
  * `limit`: an `Int` to indicate a limited number of rows to parse in a csv file; use in combination with `skipto` to read a specific, contiguous chunk within a file
  * `transpose::Bool`: read a csv file "transposed", i.e. each column is parsed as a row
  * `comment`: rows that begin with this `String` will be skipped while parsing
  * `use_mmap::Bool=!Sys.iswindows()`: whether the file should be mmapped for reading, which in some cases can be faster
  * `ignoreemptylines::Bool=false`: whether empty rows/lines in a file should be ignored (if `false`, each column will be assigned `missing` for that empty row)
* Parsing options:
  * `missingstrings`, `missingstring`: either a `String`, or `Vector{String}` to use as sentinel values that will be parsed as `missing`; by default, only an empty field (two consecutive delimiters) is considered `missing`
  * `delim=','`: a `Char` or `String` that indicates how columns are delimited in a file; if no argument is provided, parsing will try to detect the most consistent delimiter on the first 10 rows of the file
  * `ignorerepeated::Bool=false`: whether repeated (consecutive) delimiters should be ignored while parsing; useful for fixed-width files with delimiter padding between cells
  * `quotechar='"'`, `openquotechar`, `closequotechar`: a `Char` (or different start and end characters) that indicate a quoted field which may contain textual delimiters or newline characters
  * `escapechar='"'`: the `Char` used to escape quote characters in a quoted field
  * `dateformat::Union{String, Dates.DateFormat, Nothing}`: a date format string to indicate how Date/DateTime columns are formatted for the entire file
  * `decimal='.'`: a `Char` indicating how decimals are separated in floats, i.e. `3.14` used '.', or `3,14` uses a comma ','
  * `truestrings`, `falsestrings`: `Vectors of Strings` that indicate how `true` or `false` values are represented; by default only `true` and `false` are treated as `Bool`
* Column Type Options:
  * `type`: a single type to use for parsing an entire file; i.e. all columns will be treated as the same type; useful for matrix-like data files
  * `types`: a Vector or Dict of types to be used for column types; a Dict can map column index `Int`, or name `Symbol` or `String` to type for a column, i.e. Dict(1=>Float64) will set the first column as a Float64, Dict(:column1=>Float64) will set the column named column1 to Float64 and, Dict("column1"=>Float64) will set the column1 to Float64; if a `Vector` if provided, it must match the # of columns provided or detected in `header`
  * `typemap::Dict{Type, Type}`: a mapping of a type that should be replaced in every instance with another type, i.e. `Dict(Float64=>String)` would change every detected `Float64` column to be parsed as `String`
  * `pool::Union{Bool, Float64}=0.1`: if `true`, *all* columns detected as `String` will be internally pooled; alternatively, the proportion of unique values below which `String` columns should be pooled (by default 0.1, meaning that if the # of unique strings in a column is under 10%, it will be pooled)
  * `categorical::Bool=false`: whether pooled columns should be copied as CategoricalArray instead of PooledArray
  * `strict::Bool=false`: whether invalid values should throw a parsing error or be replaced with `missing`
  * `silencewarnings::Bool=false`: if `strict=false`, whether invalid value warnings should be silenced
"""
function File(source;
    # file options
    # header can be a row number, range of rows, or actual string vector
    header::Union{Integer, Vector{Symbol}, Vector{String}, AbstractVector{<:Integer}}=1,
    normalizenames::Bool=false,
    # by default, data starts immediately after header or start of file
    datarow::Integer=-1,
    skipto::Union{Nothing, Integer}=nothing,
    footerskip::Integer=0,
    limit::Integer=typemax(Int64),
    transpose::Bool=false,
    comment::Union{String, Nothing}=nothing,
    use_mmap::Bool=!Sys.iswindows(),
    ignoreemptylines::Bool=false,
    # parsing options
    missingstrings=String[],
    missingstring="",
    delim::Union{Nothing, Char, String}=nothing,
    ignorerepeated::Bool=false,
    quotechar::Union{UInt8, Char}='"',
    openquotechar::Union{UInt8, Char, Nothing}=nothing,
    closequotechar::Union{UInt8, Char, Nothing}=nothing,
    escapechar::Union{UInt8, Char}='"',
    dateformat::Union{String, Dates.DateFormat, Nothing}=nothing,
    decimal::Union{UInt8, Char}=UInt8('.'),
    truestrings::Union{Vector{String}, Nothing}=nothing,
    falsestrings::Union{Vector{String}, Nothing}=nothing,
    # type options
    type=nothing,
    types=nothing,
    typemap::Dict=EMPTY_TYPEMAP,
    categorical::Union{Bool, Real}=false,
    pool::Union{Bool, Real}=0.1,
    strict::Bool=false,
    silencewarnings::Bool=false,
    debug::Bool=false,
    parsingdebug::Bool=false,
    allowmissing::Union{Nothing, Symbol}=nothing)
    file(source, header, normalizenames, datarow, skipto, footerskip,
        limit, transpose, comment, use_mmap, ignoreemptylines, missingstrings, missingstring,
        delim, ignorerepeated, quotechar, openquotechar, closequotechar,
        escapechar, dateformat, decimal, truestrings, falsestrings, type,
        types, typemap, categorical, pool, strict, silencewarnings, debug,
        parsingdebug, allowmissing)
end

function file(source,
    # file options
    # header can be a row number, range of rows, or actual string vector
    header=1,
    normalizenames=false,
    # by default, data starts immediately after header or start of file
    datarow=-1,
    skipto=nothing,
    footerskip=0,
    limit=typemax(Int64),
    transpose=false,
    comment=nothing,
    use_mmap=!Sys.iswindows(),
    ignoreemptylines=false,
    # parsing options
    missingstrings=String[],
    missingstring="",
    delim=nothing,
    ignorerepeated=false,
    quotechar='"',
    openquotechar=nothing,
    closequotechar=nothing,
    escapechar='"',
    dateformat=nothing,
    decimal=UInt8('.'),
    truestrings=nothing,
    falsestrings=nothing,
    # type options
    type=nothing,
    types=nothing,
    typemap=EMPTY_TYPEMAP,
    categorical=false,
    pool=0.1,
    strict=false,
    silencewarnings=false,
    debug=false,
    parsingdebug=false,
    allowmissing=nothing)

    # initial argument validation and adjustment
    !isa(source, IO) && !isa(source, Vector{UInt8}) && !isfile(source) && !isa(source, Cmd) && throw(ArgumentError("\"$source\" is not a valid file"))
    (types !== nothing && any(x->!isconcretetype(x) && !(x isa Union), types isa AbstractDict ? values(types) : types)) && throw(ArgumentError("Non-concrete types passed in `types` keyword argument, please provide concrete types for columns: $types"))
    if type !== nothing && typecode(type) == EMPTY
        throw(ArgumentError("$type isn't supported in the `type` keyword argument; must be one of: `Int64`, `Float64`, `Date`, `DateTime`, `Bool`, `Missing`, `PooledString`, `CategoricalString{UInt32}`, or `String`"))
    elseif types !== nothing && any(x->typecode(x) == EMPTY, types isa AbstractDict ? values(types) : types)
        T = nothing
        for x in (types isa AbstractDict ? values(types) : types)
            if typecode(x) == EMPTY
                T = x
                break
            end
        end
        throw(ArgumentError("unsupported type $T in the `types` keyword argument; must be one of: `Int64`, `Float64`, `Date`, `DateTime`, `Bool`, `Missing`, `PooledString`, `CategoricalString{UInt32}`, or `String`"))
    end
    delim !== nothing && ((delim isa Char && iscntrl(delim) && delim != '\t') || (delim isa String && any(iscntrl, delim) && !all(==('\t'), delim))) && throw(ArgumentError("invalid delim argument = '$(escape_string(string(delim)))', must be a non-control character or string without control characters"))
    ignorerepeated && delim === nothing && throw(ArgumentError("auto-delimiter detection not supported when `ignorerepeated=true`; please provide delimiter via `delim=','`"))
    allowmissing !== nothing && @warn "`allowmissing` is a deprecated keyword argument"
    if !(categorical isa Bool)
        @warn "categorical=$categorical is deprecated in favor of `pool=$categorical`; categorical is only used to determine CategoricalArray vs. PooledArrays"
        pool = categorical
        categorical = categorical > 0.0
    elseif categorical === true
        pool = categorical
    end
    header = (isa(header, Integer) && header == 1 && (datarow == 1 || skipto == 1)) ? -1 : header
    isa(header, Integer) && datarow != -1 && (datarow > header || throw(ArgumentError("data row ($datarow) must come after header row ($header)")))
    datarow = skipto !== nothing ? skipto : (datarow == -1 ? (isa(header, Vector{Symbol}) || isa(header, Vector{String}) ? 0 : last(header)) + 1 : datarow) # by default, data starts on line after header
    debug && println("header is: $header, datarow computed as: $datarow")
    # getsource will turn any input into a `Vector{UInt8}`
    buf = getsource(source, use_mmap)
    len = length(buf)
    # skip over initial BOM character, if present
    pos = consumeBOM!(buf)

    # we call guessnrows upfront to simultaneously guess how many rows are in a file (based on average # of bytes in first 10 rows), and to figure out our delimiter: if provided, we use that, otherwise, we auto-detect based on filename or detected common delimiters in first 10 rows
    oq = something(openquotechar, quotechar) % UInt8
    cq = something(closequotechar, quotechar) % UInt8
    eq = escapechar % UInt8
    cmt = comment === nothing ? nothing : (pointer(comment), sizeof(comment))
    IG = Val(ignoreemptylines)
    rowsguess, del = guessnrows(buf, oq, cq, eq, source, delim, cmt, IG, debug)
    debug && println("estimated rows: $rowsguess")
    debug && println("detected delimiter: \"$(escape_string(del isa UInt8 ? string(Char(del)) : del))\"")

    # build Parsers.Options w/ parsing arguments
    wh1 = del == UInt(' ') || delim == " " ? 0x00 : UInt8(' ')
    wh2 = del == UInt8('\t') || delim == "\t" ? 0x00 : UInt8('\t')
    trues = truestrings === nothing ? nothing : truestrings
    falses = falsestrings === nothing ? nothing : falsestrings
    sentinel = ((isempty(missingstrings) && missingstring == "") || (length(missingstrings) == 1 && missingstrings[1] == "")) ? missing : isempty(missingstrings) ? [missingstring] : missingstrings
    options = Parsers.Options(sentinel, wh1, wh2, oq, cq, eq, del, decimal, trues, falses, dateformat, ignorerepeated, true, parsingdebug, strict, silencewarnings)

    # determine column names and where the data starts in the file; for transpose, we note the starting byte position of each column
    if transpose
        rowsguess, names, positions = datalayout_transpose(header, buf, pos, len, options, datarow, normalizenames)
        datapos = isempty(positions) ? 0 : positions[1]
    else
        positions = EMPTY_POSITIONS
        names, datapos = datalayout(header, buf, pos, len, options, datarow, normalizenames, cmt, IG)
    end
    debug && println("column names detected: $names")
    debug && println("byte position of data computed at: $datapos")

    # deduce initial column types for parsing based on whether any user-provided types were provided or not
    T = type === nothing ? EMPTY : (typecode(type) | USER)
    if types isa Vector
        typecodes = TypeCode[typecode(T) | USER for T in types]
        categorical = categorical | any(x->x == CategoricalString{UInt32}, types)
    elseif types isa AbstractDict
        typecodes = initialtypes(T, types, names)
        categorical = categorical | any(x->x == CategoricalString{UInt32}, values(types))
    else
        typecodes = TypeCode[T for _ = 1:length(names)]
    end
    debug && println("computed typecodes are: $typecodes")

    # we now do our parsing pass over the file, starting at datapos
    # we fill in our "tape", which has two UInt64 slots for each cell in row-major order (linearly indexed)
    # the 1st UInt64 is used for noting the byte position, len, and other metadata of the field within the file:
        # leftmost bit indicates a sentinel value was detected while parsing, resulting cell value will be `missing`
        # 2nd leftmost bit indicates a cell initially parsed as Int (used if column later gets promoted to Float64)
        # 3rd leftmost bit indicates if a field was quoted and included escape chararacters (will have to be unescaped later)
        # 45 bits for position (allows for maximum file size of 35TB)
        # 16 bits for field length (allows for maximum field size of 65K)
    # the 2nd UInt64 is used for storing the raw bits of a parsed, typed value: Int64, Float64, Date, DateTime, Bool, or categorical/pooled UInt32 ref
    ncols = length(names)
    tapelen = rowsguess * 2
    tapes = Vector{UInt64}[Mmap.mmap(Vector{UInt64}, tapelen) for i = 1:ncols]
    pool = pool === true ? 1.0 : pool isa Float64 ? pool : 0.0
    refs = Vector{Dict{String, UInt64}}(undef, ncols)
    lastrefs = zeros(UInt64, ncols)
    t = Base.time()
    rows, tapes = parsetape(Val(transpose), IG, ncols, gettypecodes(typemap), tapes, tapelen, buf, datapos, len, limit, cmt, positions, pool, refs, lastrefs, rowsguess, typecodes, debug, options)
    debug && println("time for initial parsing to tape: $(Base.time() - t)")
    for i = 1:ncols
        typecodes[i] &= ~USER
    end
    finaltypes = Type[TYPECODES[T] for T in typecodes]
    debug && println("types after parsing: $finaltypes, pool = $pool")
    finalrefs = Vector{Vector{String}}(undef, ncols)
    if pool > 0.0
        for i = 1:ncols
            if isassigned(refs, i)
                finalrefs[i] = map(x->x[1], sort!(collect(refs[i]), by=x->x[2]))
            elseif typecodes[i] == POOL || typecodes[i] == (POOL | MISSING)
                # case where user manually specified types, but no rows were parsed
                # so the refs never got initialized; initialize them here to empty
                finalrefs[i] = Vector{String}[]
            end
        end
    end
    return File(getname(source), names, finaltypes, rows - footerskip, ncols, eq, categorical, finalrefs, buf, tapes)
end

function parsetape(::Val{transpose}, ignoreemptylines, ncols, typemap, tapes, tapelen, buf, pos, len, limit, cmt, positions, pool, refs, lastrefs, rowsguess, typecodes, debug, options::Parsers.Options{ignorerepeated}) where {transpose, ignorerepeated}
    row = 0
    tapeidx = 1
    if pos <= len && len > 0
        while row < limit
            pos = consumecommentedline!(buf, pos, len, cmt, ignoreemptylines)
            if ignorerepeated
                pos = Parsers.checkdelim!(buf, pos, len, options)
            end
            pos > len && break
            row += 1
            for col = 1:ncols
                if transpose
                    @inbounds pos = positions[col]
                end
                @inbounds T = typecodes[col]
                @inbounds tape = tapes[col]
                type = typebits(T)
                if type === EMPTY
                    pos, code = detect(tape, tapeidx, buf, pos, len, options, row, col, typemap, pool, refs, lastrefs, debug, typecodes)
                elseif type === MISSINGTYPE
                    pos, code = detect(tape, tapeidx, buf, pos, len, options, row, col, typemap, pool, refs, lastrefs, debug, typecodes)
                elseif type === INT
                    pos, code = parseint!(T, tape, tapeidx, buf, pos, len, options, row, col, typecodes)
                elseif type === FLOAT
                    pos, code = parsevalue!(Float64, T, tape, tapeidx, buf, pos, len, options, row, col, typecodes)
                elseif type === DATE
                    pos, code = parsevalue!(Date, T, tape, tapeidx, buf, pos, len, options, row, col, typecodes)
                elseif type === DATETIME
                    pos, code = parsevalue!(DateTime, T, tape, tapeidx, buf, pos, len, options, row, col, typecodes)
                elseif type === TIME
                    pos, code = parsevalue!(Time, T, tape, tapeidx, buf, pos, len, options, row, col, typecodes)
                elseif type === BOOL
                    pos, code = parsevalue!(Bool, T, tape, tapeidx, buf, pos, len, options, row, col, typecodes)
                elseif type === POOL
                    pos, code = parsepooled!(T, tape, tapeidx, buf, pos, len, options, row, col, rowsguess, pool, refs, lastrefs, typecodes)
                else # STRING
                    pos, code = parsestring!(T, tape, tapeidx, buf, pos, len, options, row, col, typecodes)
                end
                if transpose
                    @inbounds positions[col] = pos
                else
                    if col < ncols
                        if Parsers.newline(code) || pos > len
                            options.silencewarnings || notenoughcolumns(col, ncols, row)
                            for j = (col + 1):ncols
                                # put in dummy missing values on the tape for missing columns
                                @inbounds tape = tapes[j]
                                tape[tapeidx] = MISSING_BIT
                                T = typecodes[j]
                                if T > MISSINGTYPE
                                    typecodes[j] |= MISSING
                                end
                            end
                            break # from for col = 1:ncols
                        end
                    else
                        if pos <= len && !Parsers.newline(code)
                            options.silencewarnings || toomanycolumns(ncols, row)
                            # ignore the rest of the line
                            pos = readline!(buf, pos, len, options)
                        end
                    end
                end
            end
            tapeidx += 2
            pos > len && break
            if tapeidx + 1 > tapelen
                debug && reallocatetape()
                oldtapes = tapes
                newtapelen = ceil(Int64, tapelen * 1.4)
                newtapes = Vector{UInt64}[Mmap.mmap(Vector{UInt64}, newtapelen) for i = 1:ncols]
                for i = 1:ncols
                    copyto!(newtapes[i], 1, oldtapes[i], 1, tapelen)
                end
                tapes = newtapes
                tapelen = newtapelen
                for i = 1:ncols
                    finalize(oldtapes[i])
                end
            end
        end
    end
    return row, tapes
end

@noinline reallocatetape() = println("warning: didn't pre-allocate enough tape while parsing, re-allocating...")
@noinline notenoughcolumns(cols, ncols, row) = println("warning: only found $cols / $ncols columns on data row: $row. Filling remaining columns with `missing`")
@noinline toomanycolumns(cols, row) = println("warning: parsed expected $cols columns, but didn't reach end of line on data row: $row. Ignoring any extra columns on this row")
@noinline stricterror(T, buf, pos, len, code, row, col) = throw(Error("error parsing $T on row = $row, col = $col: \"$(String(buf[pos:pos+len-1]))\", error=$(Parsers.codes(code))"))
@noinline warning(T, buf, pos, len, code, row, col) = println("warnings: error parsing $T on row = $row, col = $col: \"$(String(buf[pos:pos+len-1]))\", error=$(Parsers.codes(code))")
@noinline fatalerror(buf, pos, len, code, row, col) = throw(Error("fatal error, encountered an invalidly quoted field while parsing on row = $row, col = $col: \"$(String(buf[pos:pos+len-1]))\", error=$(Parsers.codes(code)), check your `quotechar` arguments or manually fix the field in the file itself"))

const INTVALUE = Val(true)
const NONINTVALUE = Val(false)

@inline function setposlen!(tape, tapeidx, code, pos, len, ::Val{IntValue}=NONINTVALUE) where {IntValue}
    pos = Core.bitcast(UInt64, pos) << 16
    pos |= ifelse(Parsers.sentinel(code), MISSING_BIT, UInt64(0))
    pos |= ifelse(Parsers.escapedstring(code), ESCAPE_BIT, UInt64(0))
    if IntValue
        pos |= INT_BIT
    end
    @inbounds tape[tapeidx] = pos | Core.bitcast(UInt64, len)
    return
end

function detect(tape, tapeidx, buf, pos, len, options, row, col, typemap, pool, refs, lastrefs, debug, typecodes)
    int, code, vpos, vlen, tlen = Parsers.xparse(Int64, buf, pos, len, options)
    if Parsers.invalidquotedfield(code)
        fatalerror(buf, pos, tlen, code, row, col)
    end
    if Parsers.sentinel(code) && code > 0
        setposlen!(tape, tapeidx, code, vpos, vlen)
        @inbounds typecodes[col] = MISSINGTYPE
        @goto done
    end
    @inbounds T = typecodes[col]
    if Parsers.ok(code) && !haskey(typemap, INT)
        setposlen!(tape, tapeidx, code, vpos, vlen, INTVALUE)
        @inbounds tape[tapeidx + 1] = uint64(int)
        @inbounds typecodes[col] = T == MISSINGTYPE ? (INT | MISSING) : INT
        @goto done
    end
    float, code, vpos, vlen, tlen = Parsers.xparse(Float64, buf, pos, len, options)
    if Parsers.ok(code) && !haskey(typemap, FLOAT)
        setposlen!(tape, tapeidx, code, vpos, vlen)
        @inbounds tape[tapeidx + 1] = uint64(float)
        @inbounds typecodes[col] = T == MISSINGTYPE ? (FLOAT | MISSING) : FLOAT
        @goto done
    end
    if options.dateformat === nothing
        try
            date, code, vpos, vlen, tlen = Parsers.xparse(Date, buf, pos, len, options)
            if Parsers.ok(code) && !haskey(typemap, DATE)
                setposlen!(tape, tapeidx, code, vpos, vlen)
                @inbounds tape[tapeidx + 1] = uint64(date)
                @inbounds typecodes[col] = T == MISSINGTYPE ? (DATE | MISSING) : DATE
                @goto done
            end
        catch e
        end
        try
            datetime, code, vpos, vlen, tlen = Parsers.xparse(DateTime, buf, pos, len, options)
            if Parsers.ok(code) && !haskey(typemap, DATETIME)
                setposlen!(tape, tapeidx, code, vpos, vlen)
                @inbounds tape[tapeidx + 1] = uint64(datetime)
                @inbounds typecodes[col] = T == MISSINGTYPE ? (DATETIME | MISSING) : DATETIME
                @goto done
            end
        catch e
        end
    else
        try
            # use user-provided dateformat
            DT = timetype(options.dateformat)
            dt, code, vpos, vlen, tlen = Parsers.xparse(DT, buf, pos, len, options)
            if Parsers.ok(code)
                setposlen!(tape, tapeidx, code, vpos, vlen)
                @inbounds tape[tapeidx + 1] = uint64(dt)
                @inbounds typecodes[col] = DT == Date ? (T == MISSINGTYPE ? (DATE | MISSING) : DATE) : DT == DateTime ? (T == MISSINGTYPE ? (DATETIME | MISSING) : DATETIME) : (T == MISSINGTYPE ? (TIME | MISSING) : TIME)
                @goto done
            end
        catch e
        end
    end
    bool, code, vpos, vlen, tlen = Parsers.xparse(Bool, buf, pos, len, options)
    if Parsers.ok(code) && !haskey(typemap, BOOL)
        setposlen!(tape, tapeidx, code, vpos, vlen)
        @inbounds tape[tapeidx + 1] = uint64(bool)
        @inbounds typecodes[col] = T == MISSINGTYPE ? (BOOL | MISSING) : BOOL
        @goto done
    end
    _, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
    setposlen!(tape, tapeidx, code, vpos, vlen)
    if pool > 0.0
        r = Dict{String, UInt64}()
        @inbounds refs[col] = r
        ref = getref!(r, PointerString(pointer(buf, vpos), vlen), lastrefs, col, code, options)
        @inbounds tape[tapeidx + 1] = ref
        @inbounds typecodes[col] = T == MISSINGTYPE ? (POOL | MISSING) : POOL
    else
        @inbounds typecodes[col] = T == MISSINGTYPE ? (STRING | MISSING) : STRING
    end
@label done
    return pos + tlen, code
end

@inline function parseint!(T, tape, tapeidx, buf, pos, len, options, row, col, typecodes)
    x, code, vpos, vlen, tlen = Parsers.xparse(Int64, buf, pos, len, options)
    if code > 0
        if !Parsers.sentinel(code)
            @inbounds tape[tapeidx + 1] = uint64(x)
            if !user(T)
                setposlen!(tape, tapeidx, code, vpos, vlen, INTVALUE)
            end
        else
            @inbounds typecodes[col] = INT | MISSING
            @inbounds tape[tapeidx] = MISSING_BIT
        end
    else
        if Parsers.invalidquotedfield(code)
            # this usually means parsing is borked because of an invalidly quoted field, hard error
            fatalerror(buf, pos, tlen, code, row, col)
        end
        if user(T)
            if !options.strict
                options.silencewarnings || warning(Int64, buf, pos, tlen, code, row, col)
                @inbounds typecodes[col] = INT | MISSING
                @inbounds tape[tapeidx] = MISSING_BIT
            else
                stricterror(Int64, buf, pos, tlen, code, row, col)
            end
        else
            y, code, vpos, vlen, tlen = Parsers.xparse(Float64, buf, pos, len, options)
            if code > 0
                @inbounds tape[tapeidx + 1] = uint64(y)
                @inbounds typecodes[col] = FLOAT | (missingtype(T) ? MISSING : EMPTY)
            else
                _, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
                @inbounds typecodes[col] = STRING | (missingtype(T) ? MISSING : EMPTY)
            end
            setposlen!(tape, tapeidx, code, vpos, vlen)
        end
    end
    return pos + tlen, code
end

function parsevalue!(::Type{type}, T, tape, tapeidx, buf, pos, len, options, row, col, typecodes) where {type}
    x, code, vpos, vlen, tlen = Parsers.xparse(type, buf, pos, len, options)
    if code > 0
        if !Parsers.sentinel(code)
            @inbounds tape[tapeidx + 1] = uint64(x)
        else
            @inbounds typecodes[col] = T | MISSING
            @inbounds tape[tapeidx] = MISSING_BIT
        end
    else
        if Parsers.invalidquotedfield(code)
            # this usually means parsing is borked because of an invalidly quoted field, hard error
            fatalerror(buf, pos, tlen, code, row, col)
        end
        if user(T)
            if !options.strict
                code |= Parsers.SENTINEL
                options.silencewarnings || warning(type, buf, pos, tlen, code, row, col)
                @inbounds typecodes[col] = T | MISSING
                @inbounds tape[tapeidx] = MISSING_BIT
            else
                stricterror(type, buf, pos, tlen, code, row, col)
            end
        else
            _, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
            @inbounds typecodes[col] = STRING | (missingtype(T) ? MISSING : EMPTY)
        end
    end
    if !user(T)
        setposlen!(tape, tapeidx, code, vpos, vlen)
    end
    return pos + tlen, code
end

@inline function parsestring!(T, tape, tapeidx, buf, pos, len, options, row, col, typecodes)
    x, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
    setposlen!(tape, tapeidx, code, vpos, vlen)
    if Parsers.invalidquotedfield(code)
        # this usually means parsing is borked because of an invalidly quoted field, hard error
        fatalerror(buf, pos, tlen, code, row, col)
    end
    if Parsers.sentinel(code)
        @inbounds typecodes[col] = STRING | MISSING
    end
    return pos + tlen, code
end

@inline function getref!(x::Dict, key::PointerString, lastrefs, col, code, options)
    if Parsers.escapedstring(code)
        key2 = unescape(key, options.e)
        index = Base.ht_keyindex2!(x, key2)
    else
        index = Base.ht_keyindex2!(x, key)
    end
    if index > 0
        @inbounds found_key = x.vals[index]
        return found_key::UInt64
    else
        @inbounds new = (lastrefs[col] += UInt64(1))
        @inbounds Base._setindex!(x, new, Parsers.escapedstring(code) ? key2 : String(key), -index)
        return new
    end
end

function parsepooled!(T, tape, tapeidx, buf, pos, len, options, row, col, rowsguess, pool, refs, lastrefs, typecodes)
    x, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
    setposlen!(tape, tapeidx, code, vpos, vlen)
    if Parsers.invalidquotedfield(code)
        # this usually means parsing is borked because of an invalidly quoted field, hard error
        fatalerror(buf, pos, tlen, code, row, col)
    end
    if Parsers.sentinel(code)
        @inbounds typecodes[col] = T | MISSING
        ref = UInt64(0)
    else
        if !isassigned(refs, col)
            r = Dict{String, UInt64}()
            @inbounds refs[col] = r
        else
            @inbounds r = refs[col]
        end
        ref = getref!(r, PointerString(pointer(buf, vpos), vlen), lastrefs, col, code, options)
    end
    if !user(T) && isassigned(refs, col) && (length(refs[col]) / rowsguess) > pool
        @inbounds typecodes[col] = STRING | (missingtype(typecodes[col]) ? MISSING : EMPTY)
    else
        @inbounds tape[tapeidx + 1] = ref
    end
    return pos + tlen, code
end

include("tables.jl")
include("iteration.jl")
include("rows.jl")
include("write.jl")

"""
`CSV.read(source; copycols::Bool=false, kwargs...)` => `DataFrame`

Parses a delimited file into a `DataFrame`. `copycols` determines whether a copy of columns should be made when creating the DataFrame; by default, no copy is made, and the DataFrame is built with immutable, read-only `CSV.Column` vectors. If mutable operations are needed on the DataFrame columns, set `copycols=true`.

`CSV.read` supports the same keyword arguments as [`CSV.File`](@ref).
"""
read(source; copycols::Bool=false, kwargs...) = DataFrame(CSV.File(source; kwargs...), copycols=copycols)

function __init__()
    Threads.resize_nthreads!(VALUE_BUFFERS)
    return
end

end # module


text:
"""
    Molecule

Multiple [`Localization`](@ref) objects combined into a single entity.
"""
mutable struct Molecule <: DataEntity
    index::Int
    group::LocalizationGroup
    accuracy::Float64
    Molecule(localization::Localization) =
        new(localization.index, LocalizationGroup(localization), localization.accuracy)
end

struct ResultMolecule <: DataEntity
    index::Int
    x::Float64
    y::Float64
    z::Float64
    accuracy::Float64

    ResultMolecule(mol::Molecule) = new(mol.index, mol.group.x, mol.group.y, mol.group.z, mol.accuracy)
end

import LocalizationMicroscopy.extractcoordinates
extractcoordinates(molecule::Molecule) = [molecule.group.x; molecule.group.y; molecule.group.z]

import Base.push!
function push!(molecule::Molecule, newlocalization::Localization...)
    push!(molecule.group, newlocalization...)
    calcaccuracy!(molecule)
end

function push!(molecule::Molecule, othermolecule::Molecule)
    append!(molecule.group, othermolecule.group)
end

import Base.append!
function append!(molecule::Molecule, newlocalizations::Vector{Localization})
    append!(molecule.group, newlocalizations)
    calcaccuracy!(molecule)
end

function calcaccuracy!(molecule::Molecule)
    count = length(molecule.group.localizations)
    molecule.accuracy = mean(l -> l.accuracy, molecule.group.localizations) / sqrt(count)
end

export Molecule, ResultMolecule, push!, append!


text:
for i in 1:10
    println("i is $i")
end

for name in ["Richard", "Sam", "Dawn"]
    println(name)
end

# nested loops
for i in 1:5, j in 1:3
    println(i, " / ", j)
end

# enumerating over a list
names = ("Richard", "Dawn", "Sam", "Eleanor", "Lizzy")
for (idx, val) in enumerate(names)
    println(idx, " / ", val)
end

# break or continue
for a in 1:10
    if isodd(a)
        continue
    elseif a > 8
        break
    else
        println("a is $a")
    end
end

# global values
b = 1
while b < 10
    print(b)
    global b += 1
end

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
using TrueSkill
using Test

@testset "TrueSkill.jl" begin
    # Write your tests here.
end


text:
module sbl

include("fileformats/pubchem_json.jl")

end


text:
# Autogenerated wrapper script for ghr_jll for x86_64-w64-mingw32
export ghr

JLLWrappers.@generate_wrapper_header("ghr")
JLLWrappers.@declare_executable_product(ghr)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        ghr,
        "bin\\ghr.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ContinuousWavelets, Plots, Wavelets, FFTW,Logging
global_logger(SimpleLogger(min_level=Logging.Error))
global_logger(Logging.SimpleLogger(stderr,Logging.Error))
n=2047;
t = range(0,n/1000,length=n); # 1kHz sampling rate
f = testfunction(n, "Doppler");
p1=plot(t, f,legend=false,title="Doppler",xticks=false)
c = wavelet(Morlet(π), β=2);
res = cwt(f, c)
# plotting
freqs = getMeanFreq(computeWavelets(n, c)[1])
freqs[1] = 0
p2=heatmap(t,freqs, abs.(res)', xlabel= "time (s)", ylabel="frequency (Hz)",colorbar=false)
l=@layout [a{.3h};b{.7h}]
plot(p1,p2,layout=l)
savefig("../docs/doppler.svg");#hide
DEFAULT_BOUNDARY; supertype(typeof(morl))
a,b=computeWavelets(n,c)
heatmap(a)
using LaTeXStrings
dRate = 4
waveType = Morlet()
Ψ1 = wavelet(waveType, s=8, β =dRate, averagingLength=2)
# sketch of how the frequencies are chosen
pyplot()
locs = ContinuousWavelets.polySpacing(8,Ψ1);
#Figure 3.1
scatter(1:length(locs), locs, legend=:bottomright, label="mean log frequency", xlabel="Wavelet Index (x)", ylabel= "log-Frequency (y)", color=:black)
scatter!(length(locs):length(locs), locs[end:end], markersize=10, markershape=:x, color=:black, label=:none)
firstW, lastW, stepS = ContinuousWavelets.genSamplePoints(8,Ψ1)
b = (dRate/Ψ1.Q)^(1 ./dRate)*(8+Ψ1.averagingLength)^((dRate-1)/dRate)
t= range(1,stop=length(locs),step=.1)
curve = b .*(range(firstW,stop=(locs[end]/b)^dRate,length=length(t))).^(1 / dRate)
plot!(t, curve, color=:blue, line=:dash, label=L"y=a(mx+x_0)^{^1/_\beta}", legend=:bottomright, legendfontsize=12, xrange=(0,length(locs)+3), xticks= [1; 5:5:1+length(locs)...], yrange=(minimum(locs)-1, maximum(locs)+1), yticks=(range(floor(Int,minimum(locs)), ceil(Int,maximum(locs)),step=2),[L"\alpha", (4:2:8)..., "N.Octaves"]))
x = range(15, stop=28, step=.5)
ycord(x)= locs[end] .+ b/dRate*24 .^(1/dRate-1).*(x .-length(locs))
plot!(x, ycord(x), c=:black,line=2,label=:none)
annotate!(length(locs)-1/8, locs[end]+7.5/16, Plots.text(L"\frac{dy}{dx}=^{1}/_{Q}", 11, :black, :center))
savefig("plotOfLogCentralFrequencies.svg")

n=2047
Ψ1 = wavelet(morl, s=8, β=1)
d1, ξ = computeWavelets(n,Ψ1)
Ψ2 = wavelet(morl, s=8, β =2)
d2, ξ = Wavelets.computeWavelets(n,Ψ2)
Ψ4 = wavelet(morl, s=8, β =4)
d4, ξ = Wavelets.computeWavelets(n,Ψ4)
matchingLimits = (minimum([d1 d2 d4]), maximum([d1 d2 d4]))# for
plot(heatmap(1:size(d1,2), ξ, d1, color=:Greys,
             yaxis = (L"\omega", ),
             xaxis = ("wavelet index", ),
             title=L"\beta=1"*" ("*L"\Psi1"*")", colorbar=false,
             clims=matchingLimits),
     heatmap(1:size(d2,2), ξ, d2, color=:Greys,
             yticks=[],
             xaxis = ("wavelet index", ),
             title=L"\beta=2"*" ("*L"\Psi2"*")", colorbar=false,
             clims=matchingLimits),
     heatmap(1:size(d4,2), ξ, d4,color=:Greys, yticks=[],
             xaxis = ("wavelet index", ),
             title=L"\beta=4"*" ("*L"\Psi4"*")"),
     layout=(1,3),  clims=matchingLimits,
     colorbar_title=L"\widehat{\psi_i}")#hide

savefig("changeBeta.png")
aasdf = randn(10,10)
plot(aasdf, aspect_ratio=1.0)
plot(aasdf)
l=@layout [grid(4,3) [a{.5w};b{.5w};c{.5w};d{.5w}]]
ps = [heatmap(i*j*rand(10,10),colorbar=false) for i=1:4,j=1:4]
plot(ps...,layout=l,size=(1100,800))
Q=8
eval(:Q)

using ContinuousWavelets, Plots, Wavelets, FFTW, Logging
using Plots; gr()
Plots.reset_defaults()
global_logger(Logging.SimpleLogger(stderr,Logging.Error))
n= 2047;
function mapTo(waveType, isReal=true,window=1:2047; d=1, kwargs...)
	if isReal
		c = wavelet(waveType; β=d, kwargs...)
		waves,ω = computeWavelets(n,c)
		return circshift(irfft(waves,2*n,1),(1024,0))[window,:]
	else
		c = wavelet(waveType; β=d, kwargs...)
		waves,ω = computeWavelets(n,c)
        waves = cat(waves, zeros(2047,size(waves,2)),dims=1)
		return circshift(ifft(waves,1),(1024,0))[window,:]
	end
end
tmp = mapTo(Morlet(π),false;averagingLength=-0.2)[:,2]
p1=plot([real.(tmp) imag.(tmp)], title="Morlet", labels=["real" "imaginary"], ticks=nothing,linewidth=5)
tmp = mapTo(paul2,false,averagingLength=-.5)[:,2]
p2=plot([real.(tmp) imag.(tmp)],title="Paul 2", labels=["real" "imaginary"],ticks=nothing,linewidth=5)
p3=plot(mapTo(dog2;averagingLength=-1.5)[:,2],title="derivative of gaussians (dog2)",legend=false,ticks=nothing,linewidth=5)
p4=plot(mapTo(cHaar,true; averagingLength=1)[:,2],title="Haar",legend=false,ticks=nothing,linewidth=5)
p5=plot(mapTo(cBeyl, true; d=1, averagingLength=-0)[:,2], title="Beylkyin", legend=false, ticks=nothing, linewidth=5)
p6=plot(mapTo(cVaid, true; d=1, averagingLength=-0)[:,2], title="Vaidyanthan", legend=false, ticks=nothing, linewidth=5)
p7=plot(mapTo(cDb2;d=1,averagingLength=-0)[:,2],title="Daubhechies 2",legend=false,ticks=nothing,linewidth=5)
p8=plot(mapTo(cCoif2, true;d=1, averagingLength=-0)[:,2],title="Coiflet 2",legend=false,ticks=nothing,linewidth=5)
p9=plot(mapTo(cSym4, true; d=1,averagingLength=-0)[:,2],title="Symlet 4",legend=false,ticks=nothing,linewidth=5)
k = 0600; p10=plot(mapTo(cBatt4, true, 1024-k:1024+k;d=1,averagingLength=-1)[:,2],title="Battle-Lemarie, 4",legend=false,ticks=nothing,linewidth=5)
plot(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,layout=(2,5),size=300 .*(5, 2.2))
savefig("mothers.svg")#hide
2047/2

c = wavelet(Morlet(π), averagingType=NoAve(), β=2);
exs = cat(testfunction(n, "Doppler"),testfunction(n,"Blocks"),testfunction(n,"Bumps"),testfunction(n,"HeaviSine"),dims=2)
res = cwt(exs, c)
p1=plot(plot(exs[:,1],legend=false,title="Doppler",yticks=[],xticks=[]), plot(exs[:,2],legend=false,title="Blocks",yticks=[],xticks=[]), plot(exs[:,3],legend=false,title="Bumps",yticks=[],xticks=[]), plot(exs[:,4],legend=false,title="HeaviSine",yticks=[],xticks=[]),layout=(1,4))
c = wavelet(cDb2, β=2,extraOctaves=-0);
res = cwt(exs, c)
p2=plot(heatmap(identity.(res[:,:,1])', xlabel= "time", xticks=false,yticks=[], c=:viridis, colorbar=false), heatmap(identity.(res[:,:,2])', xlabel= "time", xticks=false,yticks=[], c=:viridis, colorbar=false), heatmap(identity.(res[:,:,3])', xlabel= "time", xticks=false,yticks=[], c=:viridis, colorbar=false), heatmap(identity.(res[:,:,4])', xlabel= "time", xticks=false,yticks=[], c=:viridis, colorbar=false), layout=(1,4))
l=@layout [a{.3h};b{.7h}]
plot(p1,p2,layout=l)
savefig("multiEx.svg")
ContinuousWavelets.getNOctaves(n,wavelet(cCoif2))
log2(n)
n/2
c = wavelet(cCoif2,β=1,Q=1)
wave,ω = ContinuousWavelets.computeWavelets(n, wavelet(cCoif4,β=1,Q=1);space=true)
i=1; count(abs.(wave[:,i])/norm(wave[:,i]) .> 1e-7)
plot(wave[1000:1090,end])
ContinuousWavelets.calcDepth(c,n)
log2(n)
log2(length(qmf(c.waveType)))
log2(n)

gr()
n=2047
Ψ1 = wavelet(morl, s=8, β=1)
d1, ξ = computeWavelets(n,Ψ1)
Ψ2 = wavelet(morl, s=8, β =2)
d2, ξ = Wavelets.computeWavelets(n,Ψ2)
Ψ4 = wavelet(morl, s=8, β =4)
d4, ξ = Wavelets.computeWavelets(n,Ψ4)
matchingLimits = (minimum([d1 d2 d4]), maximum([d1 d2 d4]))# for
plot(heatmap(1:size(d1,2), ξ, d1, color=:Greys,
             yaxis = (L"\omega", ),
             xaxis = ("wavelet index", ),
             title=L"\beta=1"*" ("*L"\Psi1"*")", colorbar=false,
             clims=matchingLimits),
     heatmap(1:size(d2,2), ξ, d2, color=:Greys,
             yticks=[],
             xaxis = ("wavelet index", ),
             title=L"\beta=2"*" ("*L"\Psi2"*")", colorbar=false,
             clims=matchingLimits),
     heatmap(1:size(d4,2), ξ, d4,color=:Greys, yticks=[],
             xaxis = ("wavelet index", ),
             title=L"\beta=4"*" ("*L"\Psi4"*")"),
     layout=(1,3),  clims=matchingLimits,
     colorbar_title=L"\widehat{\psi_i}")

# readme examples
using ContinuousWavelets, Plots, Wavelets
n=2047;
f = testfunction(n, "Doppler");
p1=plot(f,legend=false,title="Doppler",xlims=(0,2000))
c = wavelet(Morlet(π), β=2);
res = cwt(f, c)
p2=heatmap(abs.(res)', xlabel= "time index",
	ylabel="frequency index",colorbar=false)
l=@layout [a{.3h};b{.7h}]
plot(p1,p2,layout=l)

f = testfunction(n, "Bumps");
p1=plot(f,legend=false,title="Bumps",xlims=(0,2000))
c = wavelet(dog2, β=2);
res = cwt(f, c)
# dropping the middle peaks
res[620:1100,:] .=0
# and smoothing the remaining peaks
res[:,10:29] .= 0
p2=heatmap(abs.(res)', xlabel= "time index", ylabel="frequency index",colorbar=false)
dropped = ContinuousWavelets.icwt(res,c,DualFrames())
p1 = plot([dropped f],legend=false, title="Smoothing and dropping bumps")
l=@layout [a{.3h};b{.7h}]
plot(p1,p2,layout=l)
savefig("../docs/bumps.svg")

n = 2047
f = testfunction(n, "Blocks") + .10randn(n)
f = testfunction(n, "HeaviSine")
wave = wavelet(morl)
fCWT = cwt(f,wave)
fRecon = ContinuousWavelets.icwt(fCWT,wave,DualFrames())
plot([real.(fRecon) f], labels=["reconstructed" "original"])
plot(fRecon)
deNoiseCWT = copy(fCWT)
deNoiseCWT[abs.(fCWT)/norm(fCWT,Inf) .< .07] .= 0
heatmap((abs.(fCWT)/norm(fCWT,Inf))' .< .01)
heatmap(fCWT', c=:viridis)
heatmap(deNoiseCWT', c=:viridis)
heatmap(fCWT' - deNoiseCWT')
fRecon = ContinuousWavelets.icwt(deNoiseCWT,wave,DualFrames())
plot([fRecon f], labels=["reconstructed" "original"])
plot(plot(f, title="with noise"), plot(fRecon, title="threshold and inverted"), layout=(2,1), legend=false)
fRecon = ContinuousWavelets.icwt(fCWT,wave,DualFrames())
plot([real.(fRecon) f])

c = wavelet(morl, β=2);
dualframe, dfNorm = getDualCoverage(n,c,DualFrames())
dualPenrose, dfNorm = getDualCoverage(n,c,PenroseDelta())
dualNaive, dfNorm = getDualCoverage(n,c,NaiveDelta())
plot([dualframe dualPenrose dualNaive], labels = ["Canonical" "Delta with least squares weights" "delta with naive weights"], legend=:topleft)

pyplot()
n=2047
Ψ1 = wavelet(morl, s=8, β=1)
d1, ξ = computeWavelets(n,Ψ1)
Ψ2 = wavelet(morl, s=8, β =2)
d2, ξ = computeWavelets(n,Ψ2)
Ψ4 = wavelet(morl, s=8, β =4)
d4, ξ = computeWavelets(n,Ψ4)
matchingLimits = (minimum([d1 d2 d4]), maximum([d1 d2 d4]))#hide
plot(heatmap(1:size(d1,2), ξ, d1, color=:Greys, yaxis = (L"\omega", ), xaxis = ("wavelet index", ), title=L"\beta=1"*" ("*L"\Psi1"*")", colorbar=false, clims=matchingLimits),  heatmap(1:size(d2,2), ξ, d2, color=:Greys, yticks=[], xaxis = ("wavelet index", ), title=L"\beta=2"*" ("*L"\Psi2"*")", colorbar=false, clims=matchingLimits),  heatmap(1:size(d4,2), ξ, d4,color=:Greys, yticks=[], xticks=[1, 5, 10, 14, 18], xaxis = ("wavelet index", ), title=L"\beta=4"*" ("*L"\Psi4"*")"), layout=(1,3), clims=matchingLimits, colorbar_title=L"\widehat{\psi_i}")
using LaTeXStrings
d4


text:
using Test, AxisKeys, BenchmarkTools

@testset "indexing & lookup" begin

    A = wrapdims(rand(2,3), 11.0:12.0, [:a, :b, :c])

    if VERSION >= v"1.2"
        # getindex
        @test 0 == @ballocated $A[1, 1]
        @test 272 >= @ballocated $A[1, :]
        @test (@inferred A[1, :]; true)
        @test (@inferred view(A, 1, :); true)

        # getkey
        @test 32 >= @ballocated $A(11.0, :a)
    end

    al_A = @ballocated view($A,1,:) # 96

    @test al_A == @ballocated $A(11.0,:)
    @test al_A == @ballocated $A(11.0)
    @test al_A == @ballocated $A(11)
    @test 0 == @ballocated AxisKeys.inferdim(11, $(axiskeys(A)))
    if VERSION >= v"1.2"
        @test (@inferred A(11); true)
    end
    @test al_A/2 >= @ballocated $A[1,:] .= 0 # dotview skips view of key vector

    # with names
    N = wrapdims(rand(2,3), row=11.0:12.0, col=[:a, :b, :c])

    if VERSION >= v"1.2"
        @test 0 == @ballocated $N[1, 1]
        @test 0 == @ballocated $N[col=1, row=1]
        @test 288 >= @ballocated $N[row=1]
        @test (@inferred N[row=1]; true)
    end

    # extraction
    @test 0 == @ballocated axiskeys($N)
    @test 0 == @ballocated axiskeys($N, 1)
    @test 0 == @ballocated axiskeys($N, :row)

    @test 0 == @ballocated dimnames($N)
    @test 0 == @ballocated dimnames($N, 1)

    @test 0 == @ballocated AxisKeys.hasnames($N)
    @test 0 == @ballocated AxisKeys.haskeys($N)

end
@testset "construction" begin

    M = rand(2,3);

    @test 64 >= @ballocated KeyedArray($M, ('a':'b', 10:10:30))
    @test 16 >= @ballocated NamedDimsArray($M, (:row, :col))
    @test (@inferred KeyedArray(M, ('a':'b', 10:10:30)); true)

    V = rand(3);
    @test 64 >= @ballocated KeyedArray($V, 'a':'c')

    # nested pair via keywords
    if VERSION >= v"1.3" # 144 alloc on 1.2, have not tried 1.1
        @test 80 >= @ballocated KeyedArray($M, row='a':'b', col=10:10:30) # 464 >=
        @test 80 >= @ballocated NamedDimsArray($M, row='a':'b', col=10:10:30) # 400 >=
    end
    @test 560 >= @ballocated wrapdims($M, row='a':'b', col=10:10:30) # 560 >=

    @test (@inferred KeyedArray(M, row='a':'b', col=10:10:30); true)
    @test (@inferred NamedDimsArray(M, row='a':'b', col=10:10:30); true)

end