text:
using Test, BenchmarkTools, StaticArrays, LinearAlgebra, KiteUtils
using KiteModels, KitePodModels

if ! @isdefined kcu
    const kcu = KCU(se())
end
if ! @isdefined kps4
    const kps4 = KPS4(kcu)
end

pos, vel = nothing, nothing
poss, vels = nothing, nothing

@testset verbose = true "KPS4 tests...." begin

function set_defaults()
    KiteModels.clear!(kps4)
    kps4.set.l_tether = 150.0
    kps4.set.elevation = 60.0
    kps4.set.area = 20.0
    kps4.set.rel_side_area = 50.0
    kps4.set.v_wind = 8.0
    kps4.set.mass = 11.4
    kps4.set.damping =  2 * 473.0
    kps4.set.alpha = 1.0/7
    kps4.set.c_s = 0.6

end

function init_392()
    KiteModels.clear!(kps4)
    kps4.set.l_tether = 392.0
    kps4.set.elevation = 70.0
    kps4.set.area = 10.0
    kps4.set.rel_side_area = 50.0
    kps4.set.v_wind = 9.1
    kps4.set.mass = 6.2
    kps4.set.c_s = 0.6
end

function init_150()
    KiteModels.clear!(kps4)
    kps4.set.l_tether = 150.0
    kps4.set.elevation = 70.0
    kps4.set.area = 10.18
    kps4.set.rel_side_area = 30.6
    kps4.set.v_wind = 9.1
    kps4.set.mass = 6.21
    kps4.set.c_s = 0.6
    kps4.set.damping = 473.0     # unit damping coefficient
    kps4.set.c_spring = 614600.0 # unit spring coefficent
    kps4.set.width = 4.9622
end

function init3()
    kps4.set.alpha =  0.08163
    KiteModels.clear!(kps4)
    kps4.set.l_tether = 150.0 # - kps4.set.height_k - kps4.set.h_bridle
    kps4.set.area = 10.18
    kps4.set.rel_side_area = 30.6
    kps4.set.mass = 6.21
    kps4.set.c_s = 0.6
    kps4.set.damping = 473.0     # unit damping coefficient
    kps4.set.c_spring = 614600.0 # unit spring coefficent
    kps4.set.width = 4.9622
    kps4.set.elevation = 70.7
    kps4.set.profile_law = Int(EXPLOG)
    pos, vel = KiteModels.init_pos_vel(kps4)
    posd = copy(vel)
    veld = zero(vel)
    height = 134.14733504839947
    kps4.v_wind .= kps4.v_wind_gnd * calc_wind_factor(kps4.am, height)
    kps4.stiffness_factor = 1.0
    KiteModels.init_springs!(kps4)
    return pos, vel, posd, veld
end

set_defaults()

@testset "calc_rho              " begin
    @test isapprox(calc_rho(kps4.am, 0.0), 1.225, atol=1e-5)
    @test isapprox(calc_rho(kps4.am, 100.0), 1.210756, atol=1e-5)
end

@testset "initial_kite_ref_frame" begin
    vec_c    = [-15., 0., -25.98076211]
    v_app    = [10.4855, 0, -3.08324]
    x, y, z = KiteModels.initial_kite_ref_frame(vec_c, v_app)
    @test x == [-0.8660254037549957, 0.0, 0.5000000000509968]
    @test y == [0.0, 1.0, 0.0]
    @test z == [-0.5000000000509968, 0.0, -0.8660254037549957]
end

@testset "get_particles         " begin
    init_150()
    particles = KiteUtils.get_particles(kps4.set.height_k, kps4.set.h_bridle, kps4.set.width, kps4.set.m_k)
    @test particles[1] == zeros(3)
    @test particles[2] == [  75.,              0.    ,       129.90381057]   # P_KCU
    @test particles[3] == [76.590521748547275, 0.    , 134.64355504845008]   # pos_A
    @test particles[4] == [78.565000000363611, 0.    , 136.07857169877312]   # pos_B (top)
    @test particles[5] == [77.450000000249887, 2.4811, 134.14733504839947]   # pos_C (left  as seen from GS)
    @test particles[6] == [77.450000000249887,-2.4811, 134.14733504839947]   # pos_D (right as seen from GS)
end

@testset "init_springs!          " begin
    init_150()
    sp = KiteModels.init_springs!(kps4)
    # test springs
    @test length(sp) == 6 + KiteModels.KITE_SPRINGS
    @test sp[1].p1 == 1
    @test sp[1].p2 == 2
    for i in 1:6
        @test sp[i].length   ≈ 150.0/6
        @test sp[i].c_spring ≈ 2.76460154e+04
        @test sp[i].damping  ≈ 1.89200000e+01
    end
    @test sp[6].p1 == 6
    @test sp[6].p2 == 7
    @test sp[7].p1 == 7 # bridle segment S1
    @test sp[7].p2 == 8
    @test sp[7].length ≈ 4.998493790987047
    @test sp[7].c_spring ≈ 54012.39452466341
    @test sp[7].damping ≈ 94.62850606174248
    @test sp[8].length ≈ 2.671692035300166       # S2
    @test sp[9].length ≈ 4.96120756              # S3
    @test sp[10].length ≈ 3.3353120022370204
    @test sp[11].length ≈ 5.4912468596622306
    @test sp[11].c_spring ≈ 49165.631334315847
    @test sp[15].length ≈  2.440379941126399     # S9 p11 p8
    # TODO also test spring 13 .. 15
end

@testset "init_masses!           " begin
    init_150()
    m = KiteModels.init_masses!(kps4)
    @test m[1] ≈ 0.1137256540599505
    for i in 2:6
        @test m[i] ≈ 0.227451308119901
    end
    @test m[7] ≈ 8.5137256540599502
    @test m[8] ≈ 2.9187
    @test m[9] ≈ 1.31652
    @test m[11] ≈ 0.98739
    @test m[11] ≈ 0.98739
end

@testset "calc_particle_forces!  " begin
    init_150()
    pos1 = KVec3(1.0, 2.0, 3.0)
    pos2 = KVec3(2.0, 3.0, 4.0)
    vel1 = KVec3(3.0, 4.0, 5.0)
    vel2 = KVec3(4.0, 5.0, 6.0)
    rho = kps4.set.rho_0
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        kps4.forces[i] .= zeros(3)
    end
    bytes = 0
    for i in 1:length(kps4.springs)
        spring = kps4.springs[i]
        kps4.stiffness_factor = 0.5
        kps4.v_wind_tether .= KVec3(8.0, 0.1, 0.0)
        bytes = @allocated KiteModels.calc_particle_forces!(kps4, pos1, pos2, vel1, vel2, spring, se().segments, se().d_tether/1000.0, rho, i)
    end
    # @test bytes == 0
    # Python output
    res=[[ 18550.4729309395152086  18550.6132232745367219  18550.6305627766196267]
         [    -0.1986161147506209      0.0819685552924057      0.1166475594582153]
         [    -0.1986161147506209      0.0819685552924057      0.1166475594582153]
         [    -0.1986161147506209      0.0819685552924057      0.1166475594582153]
         [    -0.1986161147506209      0.0819685552924057      0.1166475594582153]
         [    -0.1986161147506209      0.0819685552924057      0.1166475594582153]
         [-32417.6381463687685027 -32417.0769770286824496 -32417.0076190203544684]
         [-20528.0512582440096594 -20527.4900889039272442 -20527.420730895599263 ]
         [ 12986.35257788861054    12986.7734548936750798  12986.8254733999201562]
         [ 23289.9810739697131794  23290.5422433097955945  23290.6116013181235758]
         [ -1883.1033393325606085  -1882.5421699924754648  -1882.4728119841502121]]
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        @test all(res[i,:] .≈ kps4.forces[i])
    end
end

@testset "init                  " begin
    init_150()
    kps4.set.elevation = 60.0
    pos, vel, acc = KiteModels.init_pos_vel_acc(kps4, old=true, delta = 1e-6)
    pos1 = [[  -0.                    0.000001             -0.                ]
            [  12.5000000000000036    0.000001             21.6506350946109656]
            [  25.0000000000000071    0.000001             43.3012701892219312]
            [  37.5000000000000071    0.000001             64.9519052838329003]
            [  50.0000000000000142    0.000001             86.6025403784438623]
            [  62.5000000000000142    0.000001            108.2531754730548244]
            [  75.0000000000000142    0.000001            129.9038105676658006]
            [  76.5905217485472747    0.                  134.6435550484500823]
            [  78.5650000003636109    0.                  136.0785716987731178]
            [  77.4500000002498865    2.4810999999999996  134.1473350483994693]
            [  77.4500000002498865   -2.4810999999999996  134.1473350483994693]
           ]
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        if ! all(pos1[i,:] .≈ pos[i])
            println(i, " ", pos1[i,:], " ", pos[i])
        else
            @test all(pos1[i,:] .≈ pos[i])
        end
    end
end

function split_res(res)
    particles = se().segments+KiteModels.KITE_PARTICLES
    pos=res[1:3*particles]
    vel=res[3*particles+1:6*particles]
    pos, vel
end

@testset "initial_residual      " begin
    init3()
    res = zeros(MVector{6*(se().segments+KiteModels.KITE_PARTICLES), SimFloat})
    X00 = zeros(SimFloat, 2*(se().segments+KiteModels.KITE_PARTICLES-1)+2)
    y0, yd0 = KiteModels.init(kps4, X00)
    residual!(res, yd0, y0, kps4, 0.0)
    res_pos, res_vel = split_res(res)
    @test res_pos == zeros(length(res_pos))
    # in the second test we check if the norm of the accelerations of the tether particles is low
    # println((norm(res_vel[1:15])))
    @test_broken (norm(res_vel[1:15])) < 15.0
end

@testset "inner_loop!            " begin
    kps4.set.alpha = 1.0/7.0
    init_150()
    kps4.set.elevation = 60.0
    kps4.set.profile_law = Int(EXP)
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        kps4.forces[i] .= zeros(3)
    end
    pos, vel, acc = KiteModels.init_pos_vel_acc(kps4, old=true, delta = 1e-6)
    v_wind_gnd = KVec3(7.0, 0.1, 0.0)
    kps4.stiffness_factor = 0.5
    segments = kps4.set.segments
    d_tether = kps4.set.d_tether/1000.0
    KiteModels.inner_loop!(kps4, pos, vel, v_wind_gnd, segments, d_tether)
    forces =  [[ -1.1039795506035208  -0.0210281466470539   0.6374018106640786]
               [ -2.6112444243501161  -0.0497374193345597   1.5075837513184493]
               [ -3.2469329482256093  -0.0618458809786162   1.8746176116987205]
               [ -3.6500726231166869  -0.0695247045888339   2.1073704115748102]
               [ -3.9578038144544192  -0.0753862315039757   2.2850390976806585]
               [ -4.2101657193446407  -0.0801931096212583   2.4307403112207719]
               [-36.5076651709692399  -0.0567130032327103 -65.4202414843858833]
               [-31.8997735210856632  -0.0144731948993382  20.0026639693058534]
               [ 39.0565328469252862  -0.0092205388684263  47.4026845127173857]
               [  9.7854297156670409  84.4721933009509485   0.5027899265154456]
               [  9.7877962918810688 -84.5072746928884726   0.5210774922270867]]
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        @test isapprox(forces[i,:], kps4.forces[i], rtol=1e-4)
    end
end

@testset "calc_aero_forces!      " begin
    kps4.set.alpha = 1.0/7.0
    init_150()
    kps4.set.elevation = 60.0
    kps4.set.profile_law = Int(EXP)
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        kps4.forces[i] .= zeros(3)
    end
    pos, vel, acc = KiteModels.init_pos_vel_acc(kps4, old=true, delta = 1e-6)
    rho = 1.25
    kps4.v_wind .= KVec3(8.0, 0.2, 0.0)
    alpha_depower = 0.1
    rel_steering = -0.1
    kps4.set.alpha_zero = 5.0
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        kps4.forces[i] .= zeros(3)
    end
    KiteModels.calc_aero_forces!(kps4, pos, vel, rho, alpha_depower, rel_steering)
    forces=[[   0.                    0.                    0.                ]
            [   0.                    0.                    0.                ]
            [   0.                    0.                    0.                ]
            [   0.                    0.                    0.                ]
            [   0.                    0.                    0.                ]
            [   0.                    0.                    0.                ]
            [   0.                    0.                    0.                ]
            [   0.                    0.                    0.                ]
            [ -81.3257203383301146   -2.0331316776625457 -454.15328374043645  ]
            [  -9.9385080719600989  -68.915067335201357    -0.9904916722429121]
            [ -11.4093091631036021   53.2874848115847612    0.7727700100708267]]
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        @test all(forces[i,:] .≈ kps4.forces[i])
    end
end

function init2()
    kps4.set.alpha = 1.0/7.0
    init_150()
    kps4.set.elevation = 60.0
    kps4.set.profile_law = Int(EXP)
    pos, vel, acc = KiteModels.init_pos_vel_acc(kps4,old=true, delta = 1e-6)
    posd = copy(vel)
    veld = zero(vel)
    kps4.v_wind_gnd .= [7.0, 0.1, 0.0]
    height = 134.14733504839947
    kps4.v_wind .= kps4.v_wind_gnd * calc_wind_factor(kps4.am, height)
    kps4.stiffness_factor = 0.5
    kps4.set.alpha = 1.0/7.0
    length = 150.0
    kps4.segment_length = length/se().segments
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        kps4.forces[i] .= zeros(3)
    end
    KiteModels.init_springs!(kps4)
    return pos, vel, posd, veld
end

@testset "test_loop          " begin
    init2()
    pos, vel, posd, veld = init2()
    KiteModels.loop!(kps4, pos, vel, posd, veld)
    res1=[[-0.        0.000001 -0.      ]
          [ 0.        0.        0.      ]
          [ 0.        0.        0.      ]
          [ 0.        0.        0.      ]
          [ 0.        0.        0.      ]
          [ 0.        0.        0.      ]
          [ 0.        0.        0.      ]
          [ 0.        0.        0.      ]
          [ 0.        0.        0.      ]
          [ 0.        0.        0.      ]
          [ 0.        0.        0.      ]]
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        @test all(kps4.res1[i] .≈ res1[i,:])
    end
    forces=[[ -1.1039795506226362  -0.0210281466470539   0.6374018106309699]
            [ -2.6112444243501161  -0.0497374193345597   1.5075837513184491]
            [ -3.2469329482256093  -0.0618458809786162   1.8746176116987205]
            [ -3.6500726230975715  -0.0695247045888339   2.1073704116079188]
            [ -3.9578038144525074  -0.0753862315039757   2.2850390976839696]
            [ -4.2101657194038982  -0.0801931096212583   2.4307403111181349]
            [-36.507665170911892   -0.0567130032327103 -65.4202414842865778]
            [-31.8997735210856632  -0.0144731948993382  20.0026639693058534]
            [ 39.0565328469252862  -0.0092205388684263  47.4026845127173857]
            [  9.7854297156670409  84.4721933009509485   0.5027899265154456]
            [  9.7877962918810688 -84.5072746928884726   0.5210774922270867]]
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        @test isapprox(forces[i,:], kps4.forces[i], rtol=1e-4)
    end
    res2 = [[  0.000001             0.000001             0.000001          ]
            [ -9.4954342703640595  -0.1808633430347626  15.2921227320670887]
            [-11.8070289026385815  -0.2248941126495136  16.6267913152680755]
            [-13.2729913567184443  -0.2528171075957596  17.4731651331197071]
            [-14.3920138707363918  -0.2741317509235479  18.1192330824871632]
            [-15.3096935251050859  -0.2916113077136664  18.6490556767932176]
            [ -4.2761540463732821  -0.0066428115060276   2.1473070003763892]
            [-10.9294458221419344  -0.0049587812722576  16.663278503890723 ]
            [ 29.6664941261243911  -0.0070037210740637  45.816049670887935 ]
            [  9.9103998578748431  85.5509913012598417  10.3192110782116959]
            [  9.9127966577351092 -85.5865207191570505  10.3377321952086678]]
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        @test isapprox(res2[i,:], kps4.res2[i], rtol=1e-4)
    end
end

function unpack_add_origin(y::MVector{S, SimFloat},) where S
    T = S-2
    y_  = @view y[1:end-2]
    part  = reshape(SVector{T}(y_),  Size(3, div(T,6), 2))
    pos1, vel1 = part[:,:,1], part[:,:,2]
    pos = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(pos1[:,i-1]) end for i in 1:div(T,6)+1)
    vel = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(vel1[:,i-1]) end for i in 1:div(T,6)+1)
    pos, vel
end

function unpack(y::MVector{T, SimFloat},) where T
    part  = reshape(SVector{T}(y),  Size(3, div(T,6), 2))
    pos1, vel1 = part[:,:,1], part[:,:,2]
    pos = SVector{div(T,6)}(SVector(pos1[:,i]) for i in 1:div(T,6))
    vel = SVector{div(T,6)}(SVector(vel1[:,i]) for i in 1:div(T,6))
    pos, vel
end

@testset "test_residual!       " begin
    global pos, poss, vel, vels
    # init2()
    kps4.alpha_depower = -0.820659579962
    kps4.stiffness_factor = 0.04
    kps4.set.alpha_zero = 0.0
    res =  zeros(MVector{6*(kps4.set.segments+4), SimFloat})
    y0, yd0 = KiteModels.init(kps4; old=true, delta=1e-6)
    pos, vel = unpack_add_origin(y0)
    y0s =[  -0.                ,    0.000          ,   -0.                ,
         12.5000000000000036,    0.000001          ,   21.6506350946109656,
         25.0000000000000071,    0.000001          ,   43.3012701892219312,
         37.5000000000000071,    0.000001          ,   64.9519052838329003,
         50.0000000000000142,    0.000001          ,   86.6025403784438623,
         62.5000000000000142,    0.000001          ,  108.2531754730548244,
         75.0000000000000142,    0.000001          ,  129.9038105676658006,
         76.5905217485472747,    0.                ,  134.6435550484500823,
         78.5650000003636109,    0.                ,  136.0785716987731178,
         77.4500000002498865,    2.4810999999999996,  134.1473350483994693,
         77.4500000002498865,   -2.4810999999999996,  134.1473350483994693,
          0.00000           ,    0.00000           ,    0.00000           ,
          0.000001          ,    0.000001          ,   -0.                ,
          0.000001          ,    0.000001          ,   -0.                ,
          0.000001          ,    0.000001          ,   -0.                ,
          0.000001          ,    0.000001          ,   -0.                ,
          0.000001          ,    0.000001          ,   -0.                ,
          0.000001          ,    0.000001          ,   -0.                ,
          0.000001          ,    0.000001          ,    0.000001          ,
          0.000001          ,    0.000001          ,    0.000001          ,
          0.000001          ,    0.000001          ,    0.000001          ,
          0.000001          ,    0.000001          ,    0.000001          ,
        150.                ,    0.                ][1:end-2]
    poss, vels=unpack(MVector{66, Float64}(y0s))
    @test length(pos) == length(poss)
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        @test all(pos[i] .≈ poss[i])
        @test all(vel[i] .≈ vels[i])
    end
    time = 0.0
    @test length(res) == length(y0)-2
    residual!(res, yd0, y0, kps4, time)
    res1, res2 = kps4.res1, kps4.res2
    height = calc_height(kps4)
    @test height ≈ 136.07857169877312
    res1_= zeros(MVector{(kps4.set.segments+4+1), KVec3})
    @test res1 == res1_
    # reference values output of KPS4P_v2.py of FreeKiteSim
    res2_ =[[   0.                    0.                    0.                ]
            [  -9.4953865037803435   -0.1808625895902219    5.4822042566742599]
            [ -11.8070279026385805   -0.2248931126495136    6.816791315268075 ]
            [ -13.272990356870892    -0.2528161075957597    7.6631651328556583]
            [ -14.392012870751639    -0.2741307509235479    8.309233082460759 ]
            [ -15.3096925246324975   -0.2916103077136664    8.8390556776117659]
            [  -0.7543438292745689   -0.0077244009806506   -0.4074644786734876]
            [  -1.2766764661730894   -0.0049560787501444    0.6267331714897555]
            [-170.1485390054543245   -2.4632176272693642  327.4727436084622809]
            [ -19.5086166662826592 -108.1053016325064959   -0.6311974416990687]
            [ -22.7946534141140447  107.5032248949864737    1.2645875252505476]]

    @test res2_[1, :] == res2[1]
    @test all(res2_[2, :] .≈ res2[2])
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        if ! all(res2_[i, :] .≈ res2[i])
            # println(i, res2_[i, :], " ", res2[i])
            @test_broken all(res2_[i, :] .≈ res2[i])
        else
            @test all(res2_[i, :] .≈ res2[i])
        end
    end

end

@testset "test_getters" begin
    x, y, z = kite_ref_frame(kps4)
    @test all(x .≈ [-0.8660254037549963, 0.0, 0.5000000000509957])
    @test all(y .≈ [0.0, 1.0, 0.0])
    @test all(z .≈ [-0.5000000000509957, -0.0, -0.8660254037549963])
    @test all(orient_euler(kps4) .≈ [1.5707963267948966, -0.5235987756571836, 1.5707963267948966])
    @test all(pos_kite(kps4) .≈ [78.56500000036361, 0.0, 136.07857169877312])
    @test calc_elevation(kps4) ≈ 1.0471975512013534 # 60 degrees
    @test calc_azimuth(kps4) ≈ 0
    @test calc_heading(kps4) ≈ 0
    calc_course(kps4) # the course for vel_kite = zero is undefined
end

@testset "test_find_steady_state" begin
    init_392()
    clear!(kps4)
    KiteModels.set_depower_steering!(kps4, kps4.set.depower_offset/100.0, 0.0)
    height = sin(deg2rad(kps4.set.elevation)) * kps4.set.l_tether
    kps4.v_wind .= kps4.v_wind_gnd * calc_wind_factor(kps4.am, height)
    res1, res2 = find_steady_state!(kps4; stiffness_factor=0.035, prn=true)
    # TODO check why -9.81 appears in the residual
    @test sum(res2) ≈ -9.81*(se().segments+ KiteModels.KITE_PARTICLES) # velocity and acceleration must be near zero
    pre_tension = KiteModels.calc_pre_tension(kps4)
    @test pre_tension > 1.0001
    @test pre_tension < 1.01
    @test unstretched_length(kps4) ≈ 392.0              # initial, unstreched tether lenght
    println("length: ", tether_length(kps4))
    @test isapprox(tether_length(kps4), 406.4, rtol=1e-2) # real, streched tether length
#    @test winch_force(kps) ≈ 276.25776695110034        # initial force at the winch [N]
#    lift, drag = lift_drag(kps)
#    @test lift ≈ 443.63303000106197                    # initial lift force of the kite [N]
#    @test drag ≈ 94.25223134952152                     # initial drag force of the kite [N]
#    @test lift_over_drag(kps) ≈ 4.706870316479931      # initlal lift-over-drag
#    @test norm(v_wind_kite(kps)) ≈ 9.107670173739065   # inital wind speed at the height of the kite [m/s]
end

@testset "test_spring_forces    " begin
    init2()
    kps4.alpha_depower = -0.820659579962
    kps4.stiffness_factor = 0.04
    kps4.set.alpha_zero = 0.0
    res =  zeros(MVector{6*(kps4.set.segments+4), SimFloat})
    y0, yd0 = KiteModels.init(kps4)
    forces = spring_forces(kps4)
    ref_forces = [3.928735076156923e-12, 3.928735076156923e-12, 3.928735076156923e-12, 0.0, -3.928735076156923e-12, 1.1786205228470769e-11, 2.160277004750972, 2.1602770047433926, 2.1602770047441373, 2.1602770047074573, 2.1602770047483513, 2.1602770047483513, 2.1602770047074573, 2.1602770047433926, 2.160277004685822]
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        @test isapprox(forces[i], ref_forces[i], atol=1e-6, rtol=1e-4)
    end
end

function simulate(integrator, steps)
    start = integrator.p.iter
    for i in 1:steps
        KiteModels.next_step!(kps4, integrator)
    end
    (integrator.p.iter - start) / steps
end

@testset "test_simulate        " begin
    STEPS = 500
    kps4.set.depower = 23.6
    integrator = KiteModels.init_sim!(kps4; stiffness_factor=0.035, prn=false)
    println("\nStarting simulation...")
    simulate(integrator, 100)
    av_steps = simulate(integrator, STEPS-100)
    println(av_steps) #1102
    expected_steps = 300
    if Sys.isapple()
        expected_steps = 1000
    end
    @test isapprox(av_steps, expected_steps, rtol=0.6)
    lift, drag = KiteModels.lift_drag(kps4)
    println(lift, " ", drag) # 703.7699568972286 161.44746368100536
    @test isapprox(lift, 703.8, rtol=0.05)
    sys_state = SysState(kps4)
    # TODO Add testcase with varying reelout speed
end

# TODO Add test for winch_force
@testset "test_copy_examples" begin
    cd(tempdir())
    copy_examples()
    @test isdir("examples")
    cd("examples")
    @test isfile("plot2d.jl")
end

@testset "test_copy_bin" begin
    cd(tempdir())
    copy_bin()
    @test isdir("bin")
    cd("bin")
    @test isfile("create_sys_image")
end

end
nothing


text:
#=

  Asia medical problem in Turing.jl

  Port of *Church model https://github.com/stuhlmueller/jschurch/blob/master/tests/asia.church

  WebPPL model (exact probabilities)
  [ 'smoker', 0.6700973677390694 ],
  [ 'tb', 0.04072414942060647 ],
  [ 'cancer', 0.07774610343933955 ],
  [ 'tb_or_cancer', 0.11473103550405392 ],
  [ 'bronchitis', 1 ],
  [ 'xray', 0.15669986301877015 ],
  [ 'bronchitis', 1 ] ]

  This Turing model (using SMC())
  """
  Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

        smoker    0.6759    0.4681     0.0047    0.0073   4091.3819    1.0002     1207.2534
            tb    0.0400    0.1960     0.0020    0.0031   4124.4085    1.0004     1216.9987
        cancer    0.0719    0.2583     0.0026    0.0041   4258.4295    1.0001     1256.5446
  tb_or_cancer    0.1068    0.3089     0.0031    0.0052   4182.5913    1.0004     1234.1668
    bronchitis    1.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
          xray    0.1504    0.3575     0.0036    0.0057   4220.3910    1.0004     1245.3205
      dyspnoea    0.8080    0.3939     0.0039    0.0060   4220.4111    0.9999     1245.3264
  """

=#

using Turing
include("jl_utils.jl")

@model function asia()
    smoker ~ flip(0.5)
    tb ~ flip(smoker ? 0.05 : 0.01)
    cancer ~ flip(smoker ? 0.1 : 0.01)
    tb_or_cancer ~ Dirac(tb || cancer)
    bronchitis ~ flip(smoker ? 0.6 : 0.3)
    xray ~ flip(tb_or_cancer ? 0.98 : 0.05)
    dyspnoea ~ flip(tb_or_cancer ? bronchitis ? 0.9 : 0.7
                    : bronchitis ? 0.8 : 0.1)
    true ~ Dirac(bronchitis == true)
    true ~ Dirac(dyspnoea == true)

end

model = asia()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, IS(), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, 4)

display(chns)

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using BenchmarkTools
using Distributions

function d1(S,K,τ,r,σ)
	(log(S/K) + (r + σ^2/2) * τ) / (σ * √(τ))
end

function d2(S,K,τ,r,σ)
	d1(S,K,τ,r,σ) - σ * √(τ)
end

function Call(S,K,τ,r,σ)
	N(x) = cdf(Normal(),x)
	d₁ = d1(S,K,τ,r,σ)
	d₂ = d2(S,K,τ,r,σ)
	return N(d₁)*S - N(d₂) * K * exp(-r*τ)

end

S = 300
K = 250
τ = 1
σ = 0.15
r = 0.03

# check that result is as expected
@assert Call(S,K,τ,r,σ) ≈ 58.81976813699322

b = @benchmark Call($S,$K,$τ,$r,$σ)

@show b

max_stars_repo_name: JeannotJeannot/Learn

max_issues_repo_name: JeannotJeannot/Learn

max_forks_repo_name: JeannotJeannot/Learn


text:
using Conda
import Pkg

@show Conda.ROOTENV

env   = "py38"
pyver = "3.8"

# Assign external Conda
if false
    ENV["CONDA_JL_HOME"] = joinpath(ENV["USERPROFILE"], "local", "mambaforge")
    Pkg.build("Conda")
end

# Back to Julia Conda
if false
    # ENV["JULIA_DEPOT_PATH"]???
    ENV["CONDA_JL_HOME"] = joinpath(ENV["USERPROFILE"], ".julia", "conda", "3")
    Pkg.build("Conda")
end

# Create new Conda env
if false
    Conda.runconda(`create --name $(env) python=$(pyver)`)
end

# Common Conda commands
if false
    Conda.runconda(`info`)
    Conda.runconda(`search python`)
    Conda.runconda(`env list`)
    Conda.runconda(`list`)
end

# Add package
if false
    Conda.add("vtk",   env)
    Conda.add("numpy", env)
    # or
    Conda.runconda(`install paraview -n $(env)`)
    Conda.runconda(`install vtk      -n $(env)`)
end

function learn_basic()
    #@show Conda.ROOTENV

    # From Conda document,
    # NOTE: If you are installing Python packages for use with PyCall,
    # you must use the root environment.
    #
    # This may allow using non-root-installed packages.
    #
    # https://github.com/JuliaPy/PyCall.jl/issues/730
    # ENV["PATH"] = Conda.bin_dir(Conda.ROOTENV) * ";" * ENV["PATH"]

    #Conda.runconda(`create --name py39 python=3.9`)
    #Conda.runconda(`create --name py38 python=3.8`)

    ## env names
    #env = "py38"
    #env = "py382"

    ## Common Conda commands
    #Conda.runconda(`search python`)
    #Conda.runconda(`env list`)
    #Conda.runconda(`info`)

    # Conda.add("vtk", :py38)
    # Conda.add("numpy", :py38)
    # Conda.runconda(`install paraview`)
    # Conda.runconda(`install vtk -n $env`)

    ## Revert back to default Conda
    # ENV["CONDA_JL_HOME"] = raw"C:\Users\Yi-Kun.Yang\.julia\conda\3"
    # Pkg.build("Conda")

    ## Assign external Conda
    # ENV["CONDA_JL_HOME"] = joinpath(ENV["USERPROFILE"], "local", "mambaforge")
    # Pkg.build("Conda")

    ## Assign external Conda
    # ENV["CONDA_JL_HOME"] = raw"C:\Users\Yi-Kun.Yang\local\conda3"
    # Pkg.build("Conda")

    ## Revert back to root python
    # ENV["PYTHON"] = ""
    # Pkg.build("PyCall")

    ## Assign non-root python
    # ENV["PYTHON"] = joinpath(Conda.ROOTENV, "envs", env, "python.exe")
    # Pkg.build("PyCall")

    #ENV["PYTHON"] = raw"C:\Users\Yi-Kun.Yang\.julia\conda\3\python.exe"
    #Pkg.build("PyCall")

    # ENV["PYTHON"] = joinpath(Conda.ROOTENV, "envs", "py382", "python.exe")
    # Pkg.build("PyCall")

    # ENV["PYTHON"] = joinpath(Conda.ROOTENV, "envs", env, "python.exe")
    # Pkg.build("PyCall")

end

learn_basic()

nothing


text:
using Base.Test
using MLDatasets

include("MNIST/runtests.jl")

# CIFAR10
x, y = CIFAR10.traindata()
x, y = CIFAR10.testdata()

# CIFAR100
x, y = CIFAR100.traindata()
x, y = CIFAR100.testdata()

# PTBLM
x, y = PTBLM.traindata()
x, y = PTBLM.testdata()

# UD_English
x = UD_English.traindata()
x = UD_English.devdata()
x = UD_English.testdata()
nothing


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgApisMetaV1APIResourceList <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    groupVersion::Any # spec type: Union{ Nothing, String } # spec name: groupVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    resources::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1APIResource} } # spec name: resources

    function IoK8sApimachineryPkgApisMetaV1APIResourceList(;apiVersion=nothing, groupVersion=nothing, kind=nothing, resources=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResourceList, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResourceList, Symbol("groupVersion"), groupVersion)
        setfield!(o, Symbol("groupVersion"), groupVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResourceList, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResourceList, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1APIResourceList

const _property_map_IoK8sApimachineryPkgApisMetaV1APIResourceList = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("groupVersion")=>Symbol("groupVersion"), Symbol("kind")=>Symbol("kind"), Symbol("resources")=>Symbol("resources"))
const _property_types_IoK8sApimachineryPkgApisMetaV1APIResourceList = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("groupVersion")=>"String", Symbol("kind")=>"String", Symbol("resources")=>"Vector{IoK8sApimachineryPkgApisMetaV1APIResource}")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1APIResourceList }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1APIResourceList))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1APIResourceList }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1APIResourceList[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1APIResourceList }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1APIResourceList[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1APIResourceList)
    (getproperty(o, Symbol("groupVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("resources")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1APIResourceList }, name::Symbol, val)
end


text:
# Exceptions
# ==========

module Exceptions

"""
    MissingFieldException <: Exception

An exception type thrown when a missing field of a record is accessed.
"""
immutable MissingFieldException <: Exception
    field::Symbol
end

# This shows a warning message?
function Base.showerror(io::IO, ex::MissingFieldException)
    print(io, "$(ex.field) is missing")
end

# Throws a MissingFieldException exception.
function missingerror(field::Symbol)
    throw(MissingFieldException(field))
end

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
module O

struct Program f::Function end
const TestSet = Dict{Symbol, Program}()

filenames = ("dir/A.jl", "dir/B.jl", "dir/C.jl")
rmext(x) = replace(x, r"\.[^.]*$"=>"")

for fn in filenames
    s = rmext(basename(fn))
    m = Symbol(s)
    @eval module $m
        using ..O
        f() = "f defined in the module O." * $s # or include(fn)
        O.TestSet[Symbol($s)] = O.Program(f)
    end
end

end

# %%
O.TestSet

# %%
O.TestSet[:A].f()

# %%
O.A.f()

# %%
O.TestSet[:B].f()

# %%
O.TestSet[:A].f |> dump

# %%
O.TestSet[:B].f |> dump

# %%
O.TestSet[:C].f |> dump

# %%
module P

struct Program f::Function end
const TestSet = Dict{Symbol, Program}()

filenames = ("dir/A.jl", "dir/B.jl", "dir/C.jl")
rmext(x) = replace(x, r"\.[^.]*$"=>"")

for fn in filenames
    let m = Symbol(rmext(basename(fn)))
        f() = "f defined in " * string(m) * ".jl" # or include(fn)
        TestSet[m] = Program(f)
    end
end

end

# %%
P.TestSet

# %%
P.TestSet[:A].f()

# %%
P.TestSet[:B].f()

# %%
P.TestSet[:A].f |> dump

# %%
P.TestSet[:B].f |> dump

# %%
P.TestSet[:C].f |> dump

# %%

max_stars_repo_path: 0016/How to create many modules in a module by for loop.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
px = convert(Ptr{$(elty)},x) + (ix-1)*sizeof($(elty))
ccall(($(string(fscal)),libblas), Void,
    (Ptr{BlasInt}, Ptr{$(elty)}, Ptr{$(elty)}, Ptr{BlasInt}),
    &(n), &(a), px, &(incx))


text:
import Ecto.MirrorTypes: MirrorModel, MirrorField, MirrorStructError
import Ecto.MirrorTypes: schema_to_mirrorstruct, plain_to_mirrorstruct, mirror_to_mirrorstruct
import Ecto.MirrorTypes: mirrorstruct_to_schema, mirrorstruct_to_plain, mirrorstruct_to_mirror

# plain
type User
   name::String
   age::Int
end

# schema
module UserSchema
using Ecto.Schema
schema("users") do
    field(:name, :string)
    field(:age, :integer)
end
end

# mirror
type _User <: MirrorModel
    name::typeof(MirrorField(:name, String, Dict()))
    age::typeof(MirrorField(:name, Int, Dict()))
end

using Base.Test

name = MirrorField(:name, String, Dict())
age = MirrorField(:age, Int, Dict())
_user = _User(name, age)
for (func, typ) in [(schema_to_mirrorstruct, UserSchema),
                    (plain_to_mirrorstruct, User),
                    (mirror_to_mirrorstruct, _user)
                   ]
    struct = func(typ)
    @test :User == struct.naming.plain
    @test :_User == struct.naming.mirror
    @test :UserSchema == struct.naming.schema
    @test [:id, :name, :age] == map(field->field.name, struct.fields)
    @test [Int, String, Int] == map(field->field.typ, struct.fields)
end

struct = plain_to_mirrorstruct(User)

schema = mirrorstruct_to_schema(struct)
@test isa(schema, Module)
@test isdefined(schema, :schema)

mirror = mirrorstruct_to_mirror(struct)
@test isa(mirror, getfield(Ecto.MirrorTypes, :_User))

plain = mirrorstruct_to_plain(struct)
@test plain == getfield(Ecto.MirrorTypes, :User)

@test_throws MirrorStructError schema_to_mirrorstruct(Base)


text:
@testset "Reduce dim" begin
    let rng = MersenneTwister(123456)
        # mapreducedim on a single-dimensional array should be consistent with mapreduce.
        x = Leaf(Tape(), [1.0, 2.0, 3.0, 4.0, 5.0])
        s = 5.0 * mapreduce(abs2, +, x, dims=1)
        @test ∇(s, oneslike(unbox(s)))[x] ≈ 5.0 * [2.0, 4.0, 6.0, 8.0, 10.0]

        # mapreduce on a two-dimensional array when reduced over a single dimension
        # should give different results to mapreduce over the same array.
        x2_ = reshape([1.0, 2.0, 3.0, 4.0,], (2, 2))
        x2 = Leaf(Tape(), x2_)
        s = mapreduce(abs2, +, x2, dims=1)
        @test ∇(s, oneslike(unbox(s)))[x2] ≈ 2.0 * x2_

        # mapreducedim under `exp` should trigger the first conditional in the ∇ impl.
        x3_ = randn(rng, 5, 4)
        x3 = Leaf(Tape(), x3_)
        s = mapreduce(exp, +, x3, dims=1)
        @test ∇(s, oneslike(unbox(s)))[x3] == exp.(x3_)

        # mapreducedim under an anonymous-function should trigger fmad.
        x4_ = randn(rng, 5, 4)
        x4 = Leaf(Tape(), x4_)
        s = mapreduce(x->x*x, +, x4, dims=2)
        @test ∇(s, oneslike(unbox(s)))[x4] == 2x4_

        # Check that `sum` and `mean` work correctly with `Node`s.
        x_sum = Leaf(Tape(), randn(rng, 5, 4, 3))
        @test unbox(sum(x_sum, dims=[2, 3])) == unbox(mapreduce(identity, +, x_sum, dims=[2, 3]))
        @test unbox(mean(x_sum, dims=[2, 3])) == mean(unbox(x_sum), dims=[2, 3])

        # Ensure that the underlying value is correct in the presence of keyword arguments
        x5_ = ones(5, 4, 3)
        x5 = Leaf(Tape(), x5_)
        @test unbox(sum(x5, dims=1)) == sum(x5_, dims=1) == fill(5.0, (1, 4, 3))
        @test unbox(sum(x5, dims=2)) == sum(x5_, dims=2) == fill(4.0, (5, 1, 3))
        @test unbox(sum(x5, dims=3)) == sum(x5_, dims=3) == fill(3.0, (5, 4, 1))
        @test unbox(sum(x5)) == 60.0
        @test getfield(sum(x5), :f) === sum

        x5_ = ones(5, 4, 3)
        x5 = Leaf(Tape(), x5_)
        @test unbox(mean(x5, dims=1)) == mean(x5_, dims=1) == fill(1.0, (1, 4, 3))
        @test unbox(mean(x5, dims=2)) == mean(x5_, dims=2) == fill(1.0, (5, 1, 3))
        @test unbox(mean(x5, dims=3)) == mean(x5_, dims=3) == fill(1.0, (5, 4, 1))
        @test unbox(mean(x5)) == 1.0
        @test getfield(mean(x5), :f) === mean
        @test getfield(mean(x5, dims=1), :f) === mean
        @test check_errs(mean, randn(rng), randn(rng, 10), randn(rng, 10))
        @test check_errs(x->mean(abs2, x), randn(rng), randn(rng, 10), randn(rng, 10))
        @test check_errs(x->mean(x, dims=2), randn(rng, 10, 1),
                         randn(rng, 10, 10), randn(rng, 10, 10))
        @test check_errs(x->mean(x, dims=(2,3)), randn(rng, 10, 1, 1),
                         randn(rng, 10, 10, 10), randn(rng, 10, 10, 10))

        # Issue #123
        x6_ = collect(1:10)
        tens = (fill(10.0, (10,)), fill(10.0, (10, 1)))
        @test ∇(x->sum(sum(x, dims=2)))(x6_) == (oneslike(x6_),)
        @test ∇((x, y)->sum(sum(x, dims=2) .+ sum(y, dims=2)'))(x6_, x6_) == tens
        @test ∇((x, y)->sum(x .+ y'))(x6_, x6_) == tens

        @test check_errs(sum, randn(rng), randn(rng, 10), randn(rng, 10))
        # Specialization on the function argument works
        @test check_errs(x->sum(abs2, x), randn(rng), randn(rng, 10), randn(rng, 10))
        # Function argument with no specialization also works
        @test check_errs(x->log(sum(exp, x)), randn(rng), randn(rng, 10), randn(rng, 10))
    end
end


text:
"""
    coordinates!(mtg; angle = 45; force = false)

Compute dummy 3d coordinates for the mtg nodes using an alterning phyllotaxy. Used when
coordinates are missing.
Coordinates are just node attributes with reserved names: :XX, :YY and :ZZ.

# Returns

Nothing, mutates the mtg in-place (adds :XX, :YY and :ZZ to nodes).

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
coordinates!(mtg)
DataFrame(mtg, [:XX, :YY, :ZZ])
```
"""
function coordinates!(mtg; angle = 45, force = false)

    coord_in_attributes = [:XX, :YY, :ZZ] .∈ Ref(get_attributes(mtg))

    if !force && any(coord_in_attributes)
        error(
            "Coordinates for $([:XX, :YY, :ZZ][coord_in_attributes]) are already present",
            " in the mtg attributes (at least for some nodes). Use `force = true` if you ",
            "want to overwrite their values."
        )
    end

    phyllotaxy = [-1]
    # This function adds a XX, YY and ZZ coordinates of each node
    append!(mtg, (XX = 0.0, YY = 0.0, ZZ = 0.0))
    append!(mtg[1], (XX = 0.0, YY = 0.1, ZZ = 0.0))
    traverse!(mtg[1][1], new_pos, angle, phyllotaxy)
end

function coordinates_parent!(mtg)
    coord_in_attributes = [:XX, :YY, :ZZ] .∈ Ref(get_attributes(mtg))
    if !all(coord_in_attributes)
        error("No coordinates found in MTG. Use `coordinates!(mtg)` first.")
    end

    transform!(
        mtg,
        :XX => (x -> 1.0) => :ZZ,
        (x -> get(ancestors(x, :XX, recursivity_level = 1), 1, nothing)) => :XX_from,
        (x -> get(ancestors(x, :YY, recursivity_level = 1), 1, nothing)) => :YY_from,
        (x -> get(ancestors(x, :ZZ, recursivity_level = 1), 1, nothing)) => :ZZ_from
    )
end

function new_pos(node, angle, phyllotaxy)
    parent_node = parent(node)
    great_parent_node = parent(parent_node)

    if great_parent_node === nothing
        great_parent_node_XX = 0.0
        great_parent_node_YY = -1.0
    else
        great_parent_node_XX = great_parent_node[:XX]
        great_parent_node_YY = great_parent_node[:YY]
    end

    if node.MTG.link == "/"
        extend_length = 0.1
    else
        extend_length = 1
    end

    point = extend_pos(great_parent_node_XX, great_parent_node_YY, parent_node[:XX], parent_node[:YY], extend_length)

    if node.MTG.link == "+"
        point =
            rotate_point(
                parent_node[:XX],
                parent_node[:YY],
                point[1],
                point[2],
                phyllotaxy[1] * angle
            )

        # Change phyllotaxy for next node:
        if phyllotaxy[1] == 1
            phyllotaxy[1] = -1
        else
            phyllotaxy[1] = 1
        end
    end
    node[:XX] = point[1]
    node[:YY] = point[2]
    node[:ZZ] = 0.0

    return
end

"""
Add a new point after (x1,y1) using same direction and length relative to it
"""
function extend_pos(x0, y0, x1, y1, extend_length)
    lengthAB = sqrt((x0 - x1)^2 + (y0 - y1)^2)
    x = x1 + (x1 - x0) / lengthAB * extend_length
    y = y1 + (y1 - y0) / lengthAB * extend_length
    return x, y
end

"""
Rotate a point (x1,y1) around (x0, y0) with `angle`.
"""
function rotate_point(x0, y0, x1, y1, angle)
    angle = -angle * pi / 180
    x1 = x1 - x0
    y1 = y1 - y0
    cos_a = cos(angle)
    sin_a = sin(angle)

    x = x1 * cos_a - y1 * sin_a + x0
    y = x1 * sin_a + y1 * cos_a + y0

    return x, y
end

"""
    mtg_coordinates_df(mtg; force = false)
    mtg_coordinates_df!(mtg; force = false)

Extract the coordinates of the nodes of the mtg and the
coordinates of their parents (:XX_from, :YY_from, :ZZ_from) and output
a DataFrame.

The coordinates are computed using [`coordinates!`](@ref) if missing, or if
`force = true`.
"""
function mtg_coordinates_df(mtg; force = false)
    mtg_coordinates_df!(deepcopy(mtg); force = force)
end

function mtg_coordinates_df!(mtg; force = false)
    coord_in_attributes = [:XX, :YY, :ZZ] .∈ Ref(get_attributes(mtg))
    if !all(coord_in_attributes) || force
        if coord_in_attributes == [true, true, false] && !force
            # Only :ZZ is missing, using a dummy value:
            transform!(
                mtg,
                :XX => (x -> 0.0) => :ZZ,
            )
        else
            # :XX and/or :YY missing, recompute them
            coordinates!(mtg; force = force)
        end
    end

    # Get the coordinates of the parent of the node to draw edges:
    coordinates_parent!(mtg)

    DataFrame(mtg, [:XX, :YY, :ZZ, :XX_from, :YY_from, :ZZ_from])
end


text:
using Test, UsingMerge

foo(x::Int)=2

module Bar
  export foo
  foo(x::Float64)=3
end

@usingmerge Bar

@test length(methods(foo))==2


text:
module AnimatedTwoCategories

    using Colors, Compose
    using Cairo, Fontconfig
    using Plots
    import FileIO
    import Cairo, Fontconfig

    using Compose: circle, rectangle

    export main, create_animation, two_morphism, one_morphism,
        universal, polar_linear, cartesian_linear, AnimationType,
        Disk

    include("Assisting Files/Disk.jl")
    include("keyFrames.jl")

    set_default_graphic_size(10cm, 10cm)

    main() = main_func()

    @enum AnimationType linear polarLinear

    function create_animation(save_file_path, image_file_path, start_index, end_index, extension)
        anim = @animate for i in start_index:end_index
            file_path = string(image_file_path, string(i), extension)
            image = FileIO.load(file_path)
            plot(image, axis = nothing, background_color=:transparent)
        end

        gif(anim, save_file_path, fps = 30)
    end

    function two_morphism(file_path, expression::Expr, keyframes, frames_per_keyframe, frames_per_animation, type1::AnimationType, type2::AnimationType, polar_shift1 = (0, 0), polar_shift2 = (0,0))
        for animation_index in 1:frames_per_animation
            frame = compose(context())
            for key_frame_index in 2:length(keyframes)
                for i in 2:frames_per_keyframe
                    animation1 = universal(type1, keyframes[key_frame_index-1], keyframes[key_frame_index], (i/frames_per_keyframe), polar_shift1)
                    animation2 = universal(type2, keyframes[key_frame_index-1], keyframes[key_frame_index], (i/frames_per_keyframe), polar_shift2)
                    final = cartesian_linear(animation1, animation2, (animation_index/frames_per_animation))
                    frame = compose(context(), disk_compose_single_base(expression, final, 0), frame)
                end
            end
            draw(PNG(string(file_path, string(animation_index), ".png"), 10cm, 10cm, dpi=250), frame)
        end
    end

    function one_morphism(expression::Expr, keyframes, frames_per_keyframe, type::AnimationType, file_path, polar_shift = (0,0))
        for key_frame_index in 2:length(keyframes)
            for i in 1:frames_per_keyframe
                frame = compose(context(), disk_compose_single_base(expression, universal(type, keyframes[key_frame_index-1], keyframes[key_frame_index], (i/frames_per_keyframe), polar_shift), 0))
                draw(PNG(string(file_path, string(((key_frame_index - 1)*frames_per_keyframe) + i), ".png"), 10cm, 10cm, dpi=250), frame)
            end
        end
    end

    function universal(type::AnimationType, frame1, frame2, weight, center = (0,0))
        if type == linear
            return cartesian_linear(frame1, frame2, weight)
        elseif type == polarLinear
            return polar_linear(frame1, frame2, weight, center)
        end
    end

    function polar_linear(frame1, frame2, weight, center)
        transition = deepcopy(frame1)

        for (key, value) in transition

            radius_v = Vector{Float64}()
            for r_index in 1:length(value.radius)
                push!(radius_v, (weight * frame1[key].radius[r_index]) + ((1-weight) * frame2[key].radius[r_index]))
            end

            parameter_v = Vector{Complex}()
            for p_index in 1:length(value.parameters)
                shifted1 = frame1[key].parameters[p_index] - complex(center[1], center[2])
                shifted2 = frame2[key].parameters[p_index] - complex(center[1], center[2])

                new_angle = (weight * angle(shifted1)) + ((1-weight) * angle(shifted2))
                new_radius = (weight * abs(shifted1)) + ((1-weight) * abs(shifted2))

                push!(parameter_v, (complex(new_radius*cos(new_angle),new_radius*sin(new_angle)) + complex(center[1],center[2])))
            end

            transition[key] = Disk(weighted_color_mean(weight, frame1[key].color, frame2[key].color), radius_v, parameter_v)
        end

        return transition
    end

    function cartesian_linear(frame1, frame2, weight)
        transition = deepcopy(frame1)

        for (key, value) in transition

            radius_v = Vector{Float64}()
            for r_index in 1:length(value.radius)
                push!(radius_v, (weight * frame1[key].radius[r_index]) + ((1-weight) * frame2[key].radius[r_index]))
            end

            parameter_v = Vector{Complex}()
            for p_index in 1:length(value.parameters)
                push!(parameter_v, (weight * frame1[key].parameters[p_index]) + ((1-weight) * frame2[key].parameters[p_index]))
            end

            transition[key] = Disk(weighted_color_mean(weight, frame1[key].color, frame2[key].color), radius_v, parameter_v)
        end

        return transition
    end

end # module


text:
@testset "Cross Check Properties" begin

	Broadcast.broadcastable(x::XAM.BAM.Record) = Ref(x) #TODO: consider moving to XAM.jl.
	Broadcast.broadcastable(x::XAM.SAM.Record) = Ref(x) #TODO: consider moving to XAM.jl

	function crosscheck(bam::BAM.Record, sam::SAM.Record, property::Symbol)

		bam_property = getproperty(XAM.BAM, property)
		sam_property = getproperty(XAM.SAM, property)

		if bam_property(bam) != sam_property(sam)
			@warn "$property result is not the same" bam_property(bam) sam_property(sam)
			return false
		end

		return true
	end

	samdir = path_of_format("SAM")
	bamdir = path_of_format("BAM")

	filenames = [
		"ce#1",
		"ce#2",
		"ce#5",
		"ce#5b",
		"ce#unmap",
		"ce#unmap1",
		"ce#unmap2",
	]

	properties = [
		:position,# POS
		:tempname,# QNAME
		:mappingquality,# MAPQ
		:cigar, # CIGAR
		:flag, # FLAG
		:sequence, # SEQ
		:nextposition, # PNEXT
		:templength, # TLEN
	]

	for filename in filenames

		records_bam = collect(open(BAM.Reader, joinpath(bamdir, filename * ".bam")))
		records_sam = collect(open(SAM.Reader, joinpath(samdir, filename * ".sam")))

		for (bam, sam) in zip(records_bam, records_sam)
			@test all(crosscheck.(bam, sam, properties)) == true
		end

	end

end # testset Crosscheck


text:
# Client and server for communicating with the RAFT cluster
using Sockets

const RAFT_API_VERSION = 1
const RAFT_PROTO_VERSION = 1
const RAFT_JULIA_CLIENT  = 0x1
const RAFT_PROTO_CLIENT  = 0x2

const COMMS_PORT = 2000

function process_preamble(h::UInt64)
    size    = UInt32(h & 0xffffffff)
    version = UInt16((h >> 32) & 0xffff)
    flags   = UInt16((h >> 48) & 0xffff)
    @debug ("size=$size version=$version flags=$flags", " ", h)
    (size, version, flags)
end

function check_version(version, flags, lookup_table)
    if (flags & RAFT_JULIA_CLIENT) == RAFT_JULIA_CLIENT
        (version != RAFT_API_VERSION) && throw(RavanaException("Unsupported API version $(version)"))
    elseif (flags & RAFT_PROTO_CLIENT) == RAFT_PROTO_CLIENT
        (version != RAFT_PROTO_VERSION) && throw(RavanaException("Unsupported proto version $(version)"))
    else
        throw(RavanaException("Bad protocol header. Unknown flag $(flags)"))
    end
end

# Disassemble protocol header
function get_opt(sock, lookup_table)
    (size::UInt32, version::UInt16, flags::UInt16) = process_preamble(read(sock, UInt64))
    check_version(version, flags, lookup_table)

    (op, argv) = array_to_type(read(sock, size))
    !haskey(lookup_table, op) && throw(RavanaException("Invalid op $(op)"))
    func = lookup_table[op]
    @debug ("op: ", op, "   argv: ", argv)
    return (op, func, argv)
end

"""

"""
function ravana_server(;address=IPv4(0), port=COMMS_PORT)
    @async begin
        server = 0 # Init server socket
        sockErr = true
        @debug ("In ravana_server")
        while (sockErr)
            try
                server = listen(address, port)
                sockErr = false
                @info ("Starting cluster communication server at $(address):$(port)")
            catch e
                @info ("Could not listen on port $(port). Trying $(port + 1)")
                port += 1  # Try next port
            end
        end

        while true
            sock = accept(server)
            if isopen(sock) != true
                throw(RaftException("Error! Socket not open"))
            end
            # Disassemble op and arguments
            @async begin
                try
                    (op, func, argv) = get_opt(sock, op_table)
                    # Execute on cluster
                    ret = raft_cluster_execute(op, func, argv)
                    # Return result to client
                    b = byte_array(ret)
                    write(sock, length(b), b)
                catch e
                    @error ("ravana_server(): Exception! ", e)
                    b = byte_array(e)
                    write(sock, length(b), b)
                end
                close(sock)
            end
        end
    end
end

"""
    ravana_client(address, port, op::Int32, argv...)
Low level function that can be called from a Julia prompt/program.
```jldoxctest
julia> Ravana.ravana_client(IPv4(0), 2000, Ravana.OP_GET_NODE_PARAMS)

```
"""
function ravana_client(address, port, op::Int32, argv...)
    bytes = byte_array((op, argv))
    size = UInt32(length(bytes))
    version = UInt16(RAFT_API_VERSION)
    flags = UInt16(RAFT_JULIA_CLIENT)
    client = connect(address, port)

    ret = write(client, size, version, flags, bytes)
    ret_size = read(client, Int)
    ret = array_to_type(read(client, ret_size))
    close(client)
    ret
end

# If leader execute op, otherwise redirect to leader
function raft_cluster_execute(op, func, argv)
    @debug ("In raft_cluster_execute")
    if op == OP_INIT_CLUSTER
        ret = raft_execute(op, func, argv) # Bootstrap bypasses append entries
    elseif  current_state == LEADER
        @debug("Leader before executing raft_run_command")
        ret = raft_run_command(op, func, argv) # Commits as per RAFT protocol
    else
        ret = ravana_client(leaderAddress, leaderPort, op, argv) # Redirect to leader
    end
    ret
end

function init_cluster(;address=IPv4(0), port=2000)
    ravana_client(address, port, OP_INIT_CLUSTER, nothing)
end

function get_cluster_config(node::raftNode)
    ravana_client(node.name, node.port, OP_GET_CLUSTER_CONF, nothing)
end

get_cluster_config() = get_cluster_config(raftNode(leaderAddress, COMMS_PORT, Int128(0), Int128(0)))


text:
######################################################################
# eid.jl
# Generates the expected information density (EID).
######################################################################

"""
`fisher(m::SearchDomain)`

Creates the Fisher information matrix. Vehicle positions are discretized to same level as possible jammer locations.
"""
function fisher2(m::SearchDomain)
	# TODO: multiply by measurement noise
	# TODO: what do we do about rad2deg business? Multiply by constant?
	# TODO: what do we do about xr = yr = 0??
	n = m.num_cells
	F = zeros(n, n, n, n, 2, 2)

	for vx = 1:n
		for vy = 1:n
			for theta_x = 1:n
				for theta_y = 1:n
					# Create the 2x2 matrix
					xr = theta_x - vx
					yr = theta_y - vy
					if xr == 0 && yr == 0
						xr = 1e-6
						yr = 1e-6
					end
					den = (xr^2 + yr^2)^2
					F[vx,vy,theta_x,theta_y,1,1] = -2*yr*xr / den
					F[vx,vy,theta_x,theta_y,2,2] = 2*yr*xr / den
					F[vx,vy,theta_x,theta_y,1,2] = (xr^2 - yr^2) / den
					F[vx,vy,theta_x,theta_y,2,1] = (xr^2 - yr^2) / den
				end
			end
		end
	end
	return F
end

function fisher(m::SearchDomain)
	# TODO: multiply by measurement noise
	# TODO: what do we do about rad2deg business? Multiply by constant?
	# TODO: what do we do about xr = yr = 0??
	n = m.num_cells
	F = zeros(n, n, n, n, 2, 2)

	for vx = 1:n
		for vy = 1:n
			for theta_x = 1:n
				for theta_y = 1:n
					# compute phi
					# Create the 2x2 matrix
					xr = theta_x - vx
					yr = theta_y - vy
					if xr == 0 && yr == 0
						xr = 1e-6
						yr = 1e-6
					end

					alpha = atan2(xr, yr)
					den = (xr^2 + yr^2)

					F[vx,vy,theta_x,theta_y,1,1] = cos(alpha)^2 / den
					F[vx,vy,theta_x,theta_y,2,2] = sin(alpha)^2 / den
					F[vx,vy,theta_x,theta_y,1,2] = -sin(2*alpha) / (2*den)
					F[vx,vy,theta_x,theta_y,2,1] = -sin(2*alpha) / (2*den)
				end
			end
		end
	end
	return F
end

"""
`EID(m::SearchDomain, b::Belief)`

Arguments:

* `m` is the search domain
* `b` is the belief

Returns the expected information density.
"""
function EID(m::SearchDomain, b::Belief)
	n = m.num_cells
	F = m.F
	ei = zeros(n, n, 2, 2)
	eid = zeros(n, n)
	for vx = 1:n
		for vy = 1:n
			for theta_x = 1:n
				for theta_y = 1:n
					p = b[theta_x, theta_y]
					ei[vx, vy, 1, 1] += F[vx, vy, theta_x, theta_y,1,1] * p
					ei[vx, vy, 2, 1] += F[vx, vy, theta_x, theta_y,2,1] * p
					ei[vx, vy, 1, 2] += F[vx, vy, theta_x, theta_y,1,2] * p
					ei[vx, vy, 2, 2] += F[vx, vy, theta_x, theta_y,2,2] * p
				end
			end
			ma = ei[vx, vy, 1, 1]
			mb = ei[vx, vy, 2, 1]
			mc = ei[vx, vy, 1, 2]
			md = ei[vx, vy, 2, 2]
			eid[vx, vy] = ma*md - mb*mc
		end
	end
	return eid
end

function EID(m::SearchDomain)
	return EID(m, m.b)
end

"""
`EID(m::SearchDomain, theta_x::Int64, theta_y::Int64)`

Arguments:

* `m` is the search domain
* `(theta_x, theta_y)` is the jammer location (0-indexed)

Returns the expected information density.
"""
function EID(m::SearchDomain, theta_x::Int, theta_y::Int)
	n = m.num_cells
	F = m.F
	eid = zeros(n, n)
	for vx = 1:n
		for vy = 1:n
			ma = F[vx, vy, theta_x+1, theta_y+1, 1, 1]
			mb = F[vx, vy, theta_x+1, theta_y+1, 1, 2]
			mc = F[vx, vy, theta_x+1, theta_y+1, 2, 1]
			md = F[vx, vy, theta_x+1, theta_y+1, 2, 2]

			eid[vx, vy] = ma*md - mb*mc
		end
	end
	return eid
end


text:
function tracked_files()
    setup()
    return lock(Main.MicroCoverage_lock) do
        all_tracked_files = Vector{String}(undef, 0)
        for lnn in Main.MicroCoverage_tracker_linenumbernodes
            push!(all_tracked_files, lnn.file)
        end
        unique!(all_tracked_files)
        return all_tracked_files
    end
end


text:
addprocs(2)

using Base.Test
using DistributedStatistics, DistributedArrays, StatsBase

@testset "Histograms" begin
    n = 1000
    x = randn(n)
    w = rand(n)
    xd = distribute(x)
    wd = distribute(w)
    @test fit(Histogram, x, closed=:right) == fit(Histogram, xd, closed=:right)
    @test fit(Histogram, x, weights(w), closed=:right).weights ≈ fit(Histogram, xd, weights(wd), closed=:right).weights
end


text:
println("Testing CFADistributionEM: loglikelihood...")
srand(10)
P = 40
K = 20
Sigma_X = spdiagm(rand(P) .+ 0.2)
A = sparse(1:P, Int64[ceil(i/2) for i in 1:P], ones(P))
Sigma_L = FactorAnalysis.randcor(K, 0.2)
d = CFADistributionEM(Sigma_X, A, Sigma_L)
N = 100000
X = rand(d, N)
S = X*X' / N
Base.cov2cor!(S, sqrt(diag(S)))
@test loglikelihood(d, S, N) < 0.0
@test loglikelihood(d, S, N) > loglikelihood(CFADistributionEM(Sigma_X, A, Sigma_L .+ eye(K)*10), S, N)
function ll_test(d, S, N)
    Theta = inv(d.Sigma_X + d.A*d.Sigma_L*d.A')
    logdet(Theta) - trace(S*Theta)
end
@test abs(loglikelihood(d, S, N) - ll_test(d, S, N)) < 1e-8

println("Testing CFADistributionEM: rand...")
d = CFADistributionEM(Sigma_X, A, Sigma_L)
@test size(rand(d, N)) == (P,N)

println("Testing CFADistributionEM: fit_mle...")
truthLL = loglikelihood(d, S, N)
dopt = fit_mle(CFADistributionEM, ones(K,K), spones(A), S, N, show_trace=false, iterations=15000)
@test loglikelihood(dopt, S, N) > truthLL
@test FactorAnalysis.area_under_pr(inv(Sigma_L), inv(dopt.Sigma_L)) > 0.9

# v = FactorAnalysis.area_under_pr(inv(Sigma_L), inv(dopt.Sigma_L))
# println(truthLL)
# println(v)
# a = collect(zip(abs(FactorAnalysis.upper(inv(Sigma_L))) .> 0.01, abs(FactorAnalysis.upper(inv(dopt.Sigma_L))), 1:length(FactorAnalysis.upper(inv(Sigma_L)))))
# for v in sort(a, by=x->x[2])
#     println(v)
# end

dopt2 = fit_mle(CFADistributionSigma, ones(K,K), spones(A), S, N, show_trace=false, iterations=1000)
#FactorAnalysis.area_under_pr(abs(FactorAnalysis.upper(inv(Sigma_L))) .> 0.01, abs(FactorAnalysis.upper(inv(dopt2.Sigma_L))))
@assert minimum(eig(dopt2.Sigma_L)[1]) > 0 "CFADistributionSigma found invalid Sigma_L so can't compare to CFADistributionEM"
@assert minimum(dopt2.Sigma_X.nzval) > 0 "CFADistributionSigma found invalid Sigma_X so can't compare to CFADistributionEM"
@test abs(loglikelihood(dopt, S, N) - loglikelihood(dopt2, S, N)) < 1e-6

println("Testing CFADistributionEM: fit_map...")
dopt3 = fit_map(CFADistributionEMRidge(1.1), CFADistributionEM, spones(A), S, N, show_trace=false, iterations=15000)
@test loglikelihood(dopt3, S, N) < loglikelihood(dopt, S, N)
@test loglikelihood(dopt3, S, N) > truthLL
@test FactorAnalysis.area_under_pr(inv(Sigma_L), inv(dopt3.Sigma_L)) > 0.6 # we over regularized here so we don't expect good performance

println("Testing CFADistributionEM: fit_mle with a fully connected latent variable...")
srand(10)
P = 40
K = 20
Sigma_X = spdiagm(rand(P) .+ 0.2)
A = sparse(1:P, Int64[ceil(i/2) for i in 1:P], ones(P))
Sigma_L = FactorAnalysis.randcor(K, 0.2)
d = CFADistributionEM(Sigma_X, A, Sigma_L)
N = 100000
X = rand(d, N)

# add a confounding latent factor
confounder = randn(P)
for i in 1:N
    X[:,i] .+= randn()*confounder
end

S = X*X' / N
Base.cov2cor!(S, sqrt(diag(S)))
truthLL = loglikelihood(d, S, N)

# create an independent component to handle the confounder
maskSigma_L = ones(K+1,K+1)
for i in 1:K
    maskSigma_L[K+1,i] = maskSigma_L[i,K+1] = 0
end
A2 = sparse([1:P; 1:P], [Int64[ceil(i/2) for i in 1:P]; ones(Int64, P)*(K+1)], ones(2*P))

dopt = fit_mle(CFADistributionEM, maskSigma_L, A2, S, N, show_trace=false, iterations=15000)
Th = inv(dopt.Sigma_X + dopt.A*dopt.Sigma_L*dopt.A')
# display(dopt.A)
# tmp = deepcopy(dopt.Sigma_L)
# dopt.Sigma_L = dopt.Sigma_L * (2*eye(K+1))
# #Base.cov2cor!(dopt.Sigma_L, 2*ones(K+1))
# FactorAnalysis.normalize_Sigma_L!(dopt::CFADistributionEM)
# display(dopt.A)
# display(tmp)
# println()
# display(dopt.Sigma_L)
# @test maximum(abs(tmp .- dopt.Sigma_L)) < 1e-8
#
# display(dopt.Sigma_X + dopt.A*dopt.Sigma_L*dopt.A')
# println()
# println("GGG2")
# display(S)
# println()
#
# println(logdet(Th) - trace(S*Th))
# println(logdet(inv(S)) - trace(S*inv(S)))
# # display(Sigma_X)
# display(dopt.Sigma_X)
# display(Sigma_L)
# println()
# display(dopt.Sigma_L)
# display(dopt.A)
@test loglikelihood(dopt, S, N) > truthLL
@test FactorAnalysis.area_under_pr(inv(Sigma_L), inv(dopt.Sigma_L[1:K,1:K])) > 0.9

# println(truthLL)
#println(v)
# println()
# v = FactorAnalysis.area_under_pr(abs(FactorAnalysis.upper(inv(Sigma_L))) .> 0.01, abs(FactorAnalysis.upper(inv(dopt.Sigma_L + 0.01*eye(K)))))
# println(v)

#display(Sigma_L)


text:
using Revise
using DifferentialEquations
using AbstractPlotting
using Makie
using Random
using Printf
using NPZ
using Infiltrator
using Bem2d

function vplot(v)
    vispower = 1.0 / 10.0
    return sign.(v) .* (abs.(v)).^(vispower)
end

function plotelementsmakie(els, subscene)
    names = unique(els.name)
    for j in 1:length(names) - 1
        idx = findall(x->x == names[j], els.name)
        for i in 1:length(idx)
            plot!(subscene, [els.x1[idx[i]], els.x2[idx[i]]], [els.y1[idx[i]], els.y2[idx[i]]], color = :black)
        end
    end
    return nothing
end

function derivsconst(u, p, t)
    intidx, partials, els, eta, thetalaw, dc, blockvxglobal, blockvyglobal = p
    nintidx = length(intidx)
    vxglobal = @. abs(u[1:3:end])
    vyglobal = u[2:3:end]
    theta = @. abs(u[3:3:end])
    dtracglobaldt =  partials["trac"]["fault"]["fault"] * [blockvxglobal .- vxglobal blockvyglobal .- vyglobal]'[:]
    vx, vy = multmatvec(els.rotmat[intidx, :, :], vxglobal, vyglobal)
    dtracxglobaldt, dtracyglobaldt = multmatvec(els.rotmat[intidx, :, :], dtracglobaldt[1:2:end], dtracglobaldt[2:2:end])

    dthetadt = zeros(nintidx)
    dvxdt = zeros(nintidx)
    dvydt = zeros(nintidx)
    for i in 1:length(intidx)
        dthetadt[i] = thetalaw(vx[i], theta[i], dc)
        dvxdt[i] = 1 / (eta / els.normalstress[intidx[i]] + els.a[intidx[i]] / vx[i]) * (dtracxglobaldt[i] / els.normalstress[intidx[i]] - els.b[intidx[i]] * dthetadt[i] / theta[i])
        dvydt[i] = 0
    end

    dvxglobaldt, dvyglobaldt = multmatvec(els.rotmatinv[intidx, :, :], dvxdt, dvydt)
    dudt = zeros(length(u))
    dudt[1:3:end] = dvxglobaldt
    dudt[2:3:end] = dvyglobaldt
    dudt[3:3:end] = dthetadt
    return dudt
end

function fig_qdmakie()
    # Fun things to play with
    nsteps = 500
    amplitude = 1.0
    nfault = 200

    # Constants
    siay = 365.25 * 24 * 60 * 60
    tspan = (0, siay * 100)
    abstol = 1e-4
    reltol = 1e-4
    mu = 3e10
    nu = 0.25
    eta = 0.25
    rho = 2700.0
    eta = mu / (2.0 * sqrt(mu / rho))
    dc = 0.05
    blockvelx = 1e-9
    blockvely = 0.0

    # Create fault elements
    els = Elements(Int(1e5))
    nnodes = 1 * nfault
    faultwidth = 10000
    x1, y1, x2, y2 = discretizedline(-faultwidth, 0, faultwidth, 0, nfault)
    # x1, y1, x2, y2 = discretizedline(-faultwidth, -faultwidth, faultwidth, faultwidth, nfault)

    # Spatially variable properties
    avec = collect(LinRange(0.010, 0.020, nfault))
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.a[els.endidx + i] = avec[i]
        els.b[els.endidx + i] = 0.020
        els.normalstress[els.endidx + i] = 50e6
        els.name[els.endidx + i] = "fault"
    end
    standardize_elements!(els)

    # Convenience data structures
    idx = getidxdict(els)
    partialsconst = initpartials(els)

    # Calculate slip to traction partials on the fault
    println("Calculating velocity to traction matrix")
    @time _, _, partialsconst["trac"]["fault"]["fault"] = partialsconstdispstress(slip2dispstress, els, idx["fault"], idx["fault"], mu, nu)

    # Set initial conditions
    ics = zeros(3 * nnodes)
    ics[1:3:end] = 1e-3 * blockvelx * ones(nnodes)
    ics[2:3:end] = 0.0 * blockvely * ones(nnodes)
    ics[3:3:end] = 1e8 * ones(nnodes)
    intidx = collect(1:1:els.endidx) # indices of elements to integrate
    p = (intidx, partialsconst, els, eta, thetaaginglaw, dc, blockvelx, blockvely)
    prob = ODEProblem(derivsconst, ics, tspan, p)
    integrator = init(prob, Vern7(), abstol = abstol, reltol = reltol)

    println("Step integrating")
    xplot = collect(1:1:nfault)
    textsize = 1.75
    vgloballimits = FRect(0, -2.0, nfault, 4.0)
    thetalimits = FRect(0, -1.0, nfault, 11.0)

    vxglobalupdate = Node(0.0 .* ones(nfault))
    vxglobalupdatemax = Node(0.0 .* ones(nfault))
    vxglobalupdatemaxvals = 0.0 .* ones(nfault)
    vyglobalupdate = Node(0.0 .* ones(nfault))
    vyglobalupdatemax = Node(0.0 .* ones(nfault))
    vyglobalupdatemaxvals = 0.0 .* ones(nfault)
    vxupdate = Node(0.0 .* ones(nfault))
    vxupdatemax = Node(0.0 .* ones(nfault))
    vxupdatemaxvals = 0.0 .* ones(nfault)
    vyupdate = Node(0.0 .* ones(nfault))
    vyupdatemax = Node(0.0 .* ones(nfault))
    vyupdatemaxvals = 0.0 .* ones(nfault)
    thetaupdate = Node(0.0 .* ones(nfault))
    thetaupdatemax = Node(0.0 .* ones(nfault))
    thetaupdatemaxvals = log10.(1e-10 .* ones(nfault))
    vglobalcurrent = Node(string(0))
    vcurrent = Node(string(0))
    thetacurrent = Node(string(0))

    scene = Scene(resolution=(3000, 2000))
    subscene1 = Scene(scene, IRect(50, 1050, 1350, 800))
    subscene2 = Scene(scene, IRect(1600, 1050, 1350, 800))
    subscene3 = Scene(scene, IRect(50, 50, 1350, 800))
    subscene4 = Scene(scene, IRect(1600, 50, 1350, 800))
    plotelementsmakie(els, subscene1)
    subscene1[Axis][:names][:axisnames] = ("x (m)", "y (m)")
    plot!(subscene2, xplot, vxglobalupdate, limits=vgloballimits, color = :red)
    plot!(subscene2, xplot, vxglobalupdatemax, limits=vgloballimits, color = :red, linestyle=:dot)
    plot!(subscene2, xplot, vyglobalupdate, limits=vgloballimits, color = :blue)
    plot!(subscene2, xplot, vyglobalupdatemax, limits=vgloballimits, color = :blue, linestyle=:dot)
    text!(subscene2, vglobalcurrent, position = (0.0, 2.0), align = (:left,  :center), textsize = textsize, limits=vgloballimits)
    subscene2[Axis][:names][:axisnames] = ("element index", "pow(v global)")
    plot!(subscene4, xplot, vxupdate, limits=vgloballimits, color = :red)
    plot!(subscene4, xplot, vxupdatemax, limits=vgloballimits, color = :red, linestyle=:dot)
    plot!(subscene4, xplot, vyupdate, limits=vgloballimits, color = :blue)
    plot!(subscene4, xplot, vyupdatemax, limits=vgloballimits, color = :blue, linestyle=:dot)
    text!(subscene4, vcurrent, position = (0.0, 2.0), align = (:left,  :center), textsize = textsize, limits=vgloballimits)
    subscene4[Axis][:names][:axisnames] = ("element index", "pow(v global)")
    plot!(subscene3, xplot, thetaupdate, limits=thetalimits, color = :green)
    plot!(subscene3, xplot, thetaupdatemax, limits=thetalimits, color = :green, linestyle=:dot)
    text!(subscene3, thetacurrent, position = (0.0, 10.0), align = (:left, :center), textsize = textsize, limits=thetalimits)
    subscene3[Axis][:names][:axisnames] = ("element index", "log10(θ)")
    display(scene)

    @time for i in 1:nsteps
        DifferentialEquations.step!(integrator)

        # Global velocities
        vxglobalupdate[] = vplot(integrator.u[1:3:end])
        vyglobalupdate[] = vplot(integrator.u[2:3:end])
        vxglobalupdatemaxvals = dropdims(findmax([integrator.u[1:3:end] vxglobalupdatemaxvals], dims=2)[1], dims=2)
        vxglobalupdatemax[] = vplot(vxglobalupdatemaxvals)
        vyglobalupdatemaxvals = dropdims(findmax([integrator.u[2:3:end] vyglobalupdatemaxvals], dims=2)[1], dims=2)
        vyglobalupdatemax[] = vplot(vyglobalupdatemaxvals)
        vglobalcurrent[] = @sprintf("t = %012.6f, n = %07d, max(vx) = %01.5f, min(vx) = %01.5f, max(vy) = %01.5f, min(vy) = %01.5f", integrator.t / siay, i, maximum(integrator.u[1:3:end]), minimum(integrator.u[1:3:end]), maximum(integrator.u[2:3:end]), minimum(integrator.u[2:3:end]))

        # Local velocities
        vx, vy = multmatvec(els.rotmat[1:els.endidx, :, :], integrator.u[1:3:end], integrator.u[2:3:end])
        vxupdate[] = vplot(vx)
        vyupdate[] = vplot(vy)
        vxupdatemaxvals = dropdims(findmax([integrator.u[1:3:end] vxupdatemaxvals], dims=2)[1], dims=2)
        vxupdatemax[] = vplot(vxupdatemaxvals)
        vyupdatemaxvals = dropdims(findmax([integrator.u[2:3:end] vyupdatemaxvals], dims=2)[1], dims=2)
        vyupdatemax[] = vplot(vyupdatemaxvals)
        vcurrent[] = @sprintf("t = %012.6f, n = %07d, max(vx) = %01.5f, min(vx) = %01.5f, max(vy) = %01.5f, min(vy) = %01.5f", integrator.t / siay, i, maximum(integrator.u[1:3:end]), minimum(integrator.u[1:3:end]), maximum(integrator.u[2:3:end]), minimum(integrator.u[2:3:end]))

        # State
        thetaupdate[] = log10.(integrator.u[3:3:end])
        thetaupdatemaxvals = dropdims(findmax([integrator.u[3:3:end] thetaupdatemaxvals], dims=2)[1], dims=2)
        thetaupdatemax[] = log10.(thetaupdatemaxvals)
        thetacurrent[] = @sprintf("t = %012.6f, n = %07d, max(θ) = %012.2f, min(θ) = %012.2f", integrator.t / siay, i, maximum(integrator.u[3:3:end]), minimum(integrator.u[3:3:end]))

        sleep(1e-10)
    end

end
fig_qdmakie()


text:
"""
x = admm(A, B; ...)

Solves non-negative least-squares problem by the Alternating Direction Method
of Multipliers (ADMM).

Optional arguments:
    ρ: penalty parameter (set heuristically by default)
    ε: tolerance for stopping (small number times sqrt[m*n] by default)

References:
    S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein (2011). Distributed
    Optimization and Statistical Learning via the Alternating Direction Method
    of Multipliers. Foundations and Trends in Machine Learning.

    http://stanford.edu/~eryu/nnlsqr.html
"""
function admm(A::Matrix{Float64},
	          B::Matrix{Float64};
	          ρ=max(0.1,vecnorm(A)^2/size(A,2)),
	          ε=sqrt(size(A,2)*size(B,2))*1e-15,
	          kwargs...)

	# Dimensions
	m,k = size(A)
	n = size(B,2)

	# Cache matrices
	AtB = A'*B
	AtAρ = A'*A + eye(k)*ρ

	# Cache cholesky factorization
	L = cholfact(AtAρ)

	# Matrix storing the solutions
	X = zeros(k,n)

	# Initialize variables
	Z,U = zeros(k,n),zeros(k,n)
	X = L \ (AtB+ρ*(Z-U))

	# Solve
	while vecnorm(X-Z) > ε
		Z = max(0,X+U)
		U = U+X-Z
		X = L \ (AtB+ρ*(Z-U))
	end

	# Z ≈ X, return Z because nonnegativity is strictly enforced
	return Z
end


text:
using AdventOfCode2021
using BenchmarkTools
using Printf

on_CI = haskey(ENV, "GITHUB_ACTIONS")
if on_CI
    tune = true
else
    # Select if benchmarks should be re-tuned when run locally
    tune = false
end

# Set up benchmarking suite
suite = BenchmarkGroup()
for day in 1:CURRENT_DAY
    _Day = Day(day)
    suite["$day"] = BenchmarkGroup(["load", "solve1", "solve2"])
    suite["$day"]["load"] = @benchmarkable data($_Day)
    _data = data(_Day)
    suite["$day"]["solve1"] = @benchmarkable solve1($_Day, $_data)
    suite["$day"]["solve2"] = @benchmarkable solve2($_Day, $_data)
end

loadparams!(suite, BenchmarkTools.load("benchmark/params.json")[1], :evals, :samples);

# Re-tune if selected
if tune
    print("Tuning benchmark suite... ")
    tune!(suite)
    BenchmarkTools.save("benchmark/params.json", params(suite))
    println("params saved.")
end

# Run bechmarks
println("Running benchmarks... ")
results = run(suite; verbose=true)

# Parse results to table-rows for readme
formattime(t) = @sprintf "%.2f μs" t / 1000 # where t = time in ns
formatmemory(m) = @sprintf "%.2f KiB" m / 1024 # where m = memory in bytes

function day2string(day::Integer)
    two_digit_day = @sprintf "%02d" day

    # Some hand-crafted strings
    url_task = "https://adventofcode.com/2021/day/$day"
    url_src = "https://github.com/adrhill/AdventOfCode2021.jl/blob/main/src/day$two_digit_day.jl"
    day_string = """[Day $day]($url_task)
[`src`]($url_src)"""
    functions_string = "`load`
`solve1`
`solve2`"

    # Strings from benchmarks
    res = results["$day"]
    ks = ["load", "solve1", "solve2"]
    times = join([formattime(mean(res[k]).time) for k in ks], "
")
    mems = join([formatmemory(mean(res[k]).memory) for k in ks], "
")
    allocs = join([mean(res[k]).allocs for k in ks], "
")
    return join([day_string, functions_string, times, mems, allocs], " | ")
end

# Handwritten README stuff
ci_badge = "https://github.com/adrhill/AdventOfCode2021.jl/actions/workflows/CI.yml/badge.svg?branch=main"
ci_url = "https://github.com/adrhill/AdventOfCode2021.jl/actions/workflows/CI.yml?query=branch%3Amain"
cov_badge = "https://codecov.io/gh/adrhill/AdventOfCode2021.jl/branch/main/graph/badge.svg"
cov_url = "https://codecov.io/gh/adrhill/AdventOfCode2021.jl"
readme_header = """
# AdventOfCode2021.jl

[![Build Status]($ci_badge)]($ci_url)
[![Coverage]($cov_badge)]($cov_url)

Solutions for the [Advent of Code 2021](https://adventofcode.com/2021) in Julia Base without dependencies.

## Benchmarks
Timings may vary across commits as they are automatically generated by GitHub Actions on different machines:

| Day | Part | Mean time | Memory | Allocs |
|:---:|:----:|----------:|-------:|-------:|"""

# Write table to file
open(joinpath(@__DIR__, "..", "README.md"), "w") do io
    write(io, join([readme_header, [day2string(d) for d in 1:CURRENT_DAY]...], "\n"))
end;


text:
module ExpmV

  # package code goes here
  include("normAm.jl")

  include("select_taylor_degree.jl")

  include("expmv_fun.jl")

end # module


text:
function isWorldmapData(dat::DataFrames.DataFrame, latlong = true)
  DataFrames.ncol(dat) == 5 || return false

  if eltypet(dat[:, 1]) <: AbstractString
    if eltypet(dat[:, 4]) <: Number
      if eltypet(dat[:, 5]) <: Number
        latlong || return true # risky
        if minimum(dropna(dat[:, 4])) > -181 && maximum(dropna(dat[:, 4])) < 181
          if minimum(dropna(dat[:, 5])) > -91 && maximum(dropna(dat[:, 5])) < 91
            return true
          end
        end
      end
    end
  end
  false
end

function parsesingleDataFrame(occ::DataFrames.DataFrame)
    if isWorldmapData(occ)
      println("Data format identified as Worldmap export file")
      coords = occ[!,4:5]
      coords[!,:sites] = createsitenames(coords)
      occ = DataFrame(site = coords[!,:sites], abu = ones(Int, DataFrames.nrow(occ)), species = occ[!,1])
      coords = unique(coords, :sites)
    else
        error("If coords is not specified, data must be in the worldmap format - of the column format [speciesnames, (ignored), (ignored), longitude, latitude]")
    end
  occ, coords
end

function guess_xycols(dat::DataFrames.DataFrame)
  numbers = map(x -> x<:Number, eltypest(dat))
  sum(.!numbers) == 1 || error("Site names cannot be numeric in the input matrix")
  ((findall(numbers)[1:2])...,)
end

function testbool(x::Int)
    x == 1 && return true
    x > 1 && error("Integer values > 1 used to create Boolean matrix")
    return false
end

testbool(x) = error("Value can not be interpreted as Boolean")
testbool(x::Missing) = false
testbool(x::Bool) = x
function testbool(x::Number)
  x == 0 && return false
  x == 1 && return true
  error("Value can not be interpreted as Boolean")
end

function dataFrametoSparseMatrix(dat::DataFrames.DataFrame, ::Type{T}) where T<:Bool
    is, js = Vector{Int}(), Vector{Int}()

    @inbounds for j in 1:DataFrames.ncol(dat)
        col = dat[:,j]
        for i in 1:DataFrames.nrow(dat)
            if testbool(col[i])
                push!(is, i)
                push!(js, j)
            end
        end
    end

    sparse(is, js, true, DataFrames.nrow(dat), DataFrames.ncol(dat))
end

function dataFrametoSparseMatrix(dat::DataFrames.DataFrame, ::Type{T}) where T<:Real
    is, js, vals = Vector{Int}(), Vector{Int}(), Vector{T}()

    @inbounds for j in 1:DataFrames.ncol(dat)
        col = dat[:,j]
        for i in 1:DataFrames.nrow(dat)
            if !ismissing(col[i]) && col[i] != 0
                push!(is, i)
                push!(js, j)
                push!(vals, col[i])
            end
        end
    end

    sparse(is, js, vals, DataFrames.nrow(dat), DataFrames.ncol(dat))
end

function match_commat_coords(occ::ComMatrix, coords::AbstractMatrix, sitestats::DataFrames.DataFrame)
  occ, coords, sitestats
 ## so far this does nothing TODO
end

function dropspecies!(occ::SpeciesData)
  occur = occurring(occ)
  occ.commatrix = occ.commatrix[occur, :]
  occ.traits = occ.traits[occur,:]
end

function dropbyindex!(site::Locations{PointData}, indicestokeep)
  site.coords = site.coords[indicestokeep,:]
  site.sitestats = site.sitestats[indicestokeep,:]
end

# these functions will be removed eventually
maxrange(x) = diff([extrema(x)...])[1]

# remember here - something wrong with the indices, make sure they are based from 1!

function dropbyindex!(site::Locations{GridData}, indicestokeep)
  site.indices = site.indices[indicestokeep,:]
  site.sitestats = site.sitestats[indicestokeep,:]
  site.grid.xmin = xrange(site.grid)[minimum(site.indices[:,1])]
  site.grid.ymin = yrange(site.grid)[minimum(site.indices[:,2])]
  site.grid.xcells = maxrange(site.indices[:,1]) + 1
  site.grid.ycells = maxrange(site.indices[:,2]) + 1
  site.indices = site.indices - minimum(site.indices) + 1
end

function dropsites!(occ::ComMatrix, site::SELocations)
  hasspecies = occupied(occ)
  occ.commatrix = occ.commatrix[:, hasspecies]
  dropbyindex!(site, hasspecies)
end

function createsitenames(coords::AbstractMatrix)
  size(coords, 2) == 2 || error("Only defined for matrices with two columns")
  mapslices(x->"$(x[1])_$(x[2])", coords, dims = 2)
end

function createsitenames(coords::DataFrames.DataFrame)
  size(coords, 2) == 2 || error("Only defined for matrices with two columns")
  ["$(coords[i,1])_$(coords[i,2])" for i in 1:DataFrames.nrow(coords)]
end

creategrid(coords::AbstractMatrix{<:Union{Number, Missing}}, tolerance = sqrt(eps())) =
    GridTopology(gridvar(coords[:,1], tolerance), gridvar(coords[:,2], tolerance))

# could allow for n-dimensional binning, using code from StatsBase.Histogram
function gridvar(x, tolerance = sqrt(eps()))
  sux = sort(unique(x))
  difx = diff(sux)
  length(difx) == 0 && error("Cannot make a grid with width 1 in the current implementation") #TODO
  rudifxmin, rudifxmax = extrema(unique(difx))
  if rudifxmax/rudifxmin < tolerance
    filter!(x->x>rudifxmax * tolerance, difx)
    rudifxmin, rudifxmax = extrema(unique(difx))
  end

  err1 = rudifxmax - rudifxmin
  if err1 > tolerance
    xx = [1, rudifxmax / rudifxmin]
    err2 = maximum(abs.(floor.(xx) - xx))
    err2  > tolerance && error("Cannot be converted to grid, as coordinate intervals are not constant. Try adjusting the tolerance (currently $tolerance)")
    difx = difx[difx .< rudifxmin + tolerance]
  end

  cellsize = mean(difx)
  min = minimum(sux)
  cellnumber = round(Int, maxrange(sux) / cellsize) + 1

  range(min, step = cellsize, length = cellnumber)
end

function getindices(coords::AbstractMatrix{<:Union{AbstractFloat, Missing}}, grid::GridTopology, tolerance = 2*sqrt(eps()))
  index1 = 1 .+ floor.(Int,(coords[:,1] .- xmin(grid) .+ 0.5xcellsize(grid)) ./ xcellsize(grid) .+ tolerance)
  index2 = 1 .+ floor.(Int,(coords[:,2] .- ymin(grid) .+ 0.5xcellsize(grid)) ./ ycellsize(grid) .+ tolerance)
  hcat(index1, index2)
end


text:
using FastTransforms, Compat
using Compat.Test
if VERSION ≥ v"0.7-"
    using FFTW
end

@testset "BigFloat FFT and DCT" begin

    c = collect(range(-big(1.0),stop=1,length=16))
    @test norm(fft(c) - fft(Float64.(c))) < 3Float64(norm(c))*eps()
    @test norm(ifft(c) - ifft(Float64.(c))) < 3Float64(norm(c))*eps()

    c = collect(range(-big(1.0),stop=1.0,length=201))
    @test norm(ifft(fft(c))-c) < 200norm(c)eps(BigFloat)

    p = plan_dct(c)
    @test norm(dct(c) - p*c) == 0

    pi = plan_idct!(c)
    @test norm(pi*dct(c) - c) < 1000norm(c)*eps(BigFloat)

    @test norm(dct(c)-dct(map(Float64,c)),Inf) < 10eps()

    cc = cis.(c)
    @test norm(dct(cc)-dct(map(Complex{Float64},cc)),Inf) < 10eps()

    c = big.(rand(100)) + im*big.(rand(100))
    @test norm(dct(c)-dct(map(ComplexF64,c)),Inf) < 10eps()
    @test norm(idct(c)-idct(map(ComplexF64,c)),Inf) < 10eps()
    @test norm(idct(dct(c))-c,Inf) < 1000eps(BigFloat)
    @test norm(dct(idct(c))-c,Inf) < 1000eps(BigFloat)
end


text:
@coeff function solve_pulsed!(reactor::Reactor, R_0_guess=nothing)
  LHS_this = log((1+epsilon_B)/(1-epsilon_B))/(1-epsilon_B)
  LHS_that = 4*epsilon_B/(1-epsilon_B^2)

  LHS = 1 - epsilon_B
  LHS -= B_0^2 / (4*mu_0*sigma_max) * (LHS_this+LHS_that)

  RHS_this = (
    1 + R_0^2 / ( ( kappa * a + b + c ) ^ 2 )
  ) ^ (1/4)

  RHS_that = sqrt( 1 - f_B + Tau )

  RHS = K_OH * RHS_this * RHS_that / ( B_0 * R_0 ^ 2 ) ^ (5/4)

  equation = real(LHS - RHS)

  root_function = function (tmp_R_0)
    tmp_error = subs(equation, symbols(:R_0) => tmp_R_0)

    if isa(nu_J, SymEngine.Basic)
      tmp_nu_J = calc_nu_J(reactor, tmp_R_0, B_0)
      tmp_error = subs(tmp_error, symbols(:nu_J) => tmp_nu_J)
    end

    float(real(tmp_error))
  end

  if isnothing(R_0_guess) || isnan(R_0_guess)
    valid_roots = []
  else
    work_R_0 = NaN

    try
      work_R_0 = find_zero(root_function, R_0_guess)
    catch cur_error
      isa(cur_error, DomainError) || rethrow(cur_error)
    end

    if is_valid(reactor, work_R_0, B_0)
      valid_roots = [work_R_0]
    else
      valid_roots = []
    end
  end

  if isempty(valid_roots)
    function_roots = find_zeros(root_function, R_0_min, R_0_max)

    for work_R_0 in function_roots
      is_valid(reactor, work_R_0, B_0) || continue
      push!(valid_roots, work_R_0)
    end
  end

  if isempty(valid_roots)
    reactor.R_0 = NaN
  else
    @assert length(valid_roots) == 1
    chosen_R_0 = first(valid_roots)

    reactor.R_0 = chosen_R_0

    if isa(nu_J, SymEngine.Basic)
      reactor.nu_J = calc_nu_J(reactor, chosen_R_0, B_0)
    end
  end
end


text:
using Documenter, NeutralLandscapes

makedocs(
    sitename="Neutral Landscapes",
    authors="Timothée Poisot",
    modules=[NeutralLandscapes],
    pages=[
        "Index" => "index.md",
        "Gallery" => "gallery.md"
        ]
)

deploydocs(
    deps=Deps.pip("pygments", "python-markdown-math"),
    repo="github.com/EcoJulia/NeutralLandscapes.jl.git",
    devbranch="main",
    push_preview=true
)


text:
using Tables, TableTraits, IteratorInterfaceExtensions

Tables.istable(::Type{<:AbstractDataFrame}) = true
Tables.columnaccess(::Type{<:AbstractDataFrame}) = true
Tables.columns(df::AbstractDataFrame) = df
Tables.rowaccess(::Type{DataFrame}) = true
Tables.rows(df::DataFrame) = Tables.rows(columntable(df))

Tables.schema(df::AbstractDataFrame) = Tables.Schema(names(df), eltypes(df))

getvector(x::AbstractVector) = x
getvector(x) = collect(x)
fromcolumns(x) = DataFrame(Any[getvector(c) for c in Tables.eachcolumn(x)], Index(collect(Symbol, propertynames(x))))

function DataFrame(x)
    if Tables.istable(typeof(x))
        return fromcolumns(Tables.columns(x))
    end
    if TableTraits.supports_get_columns_copy_using_missing(x)
        y = TableTraits.get_columns_copy_using_missing(x)
        return DataFrame(collect(y), collect(keys(y)))
    end
    it = TableTraits.isiterabletable(x)
    if it === true
        # Base.depwarn("constructing a DataFrame from an iterator is deprecated; $T should support the Tables.jl interface", nothing)
        y = IteratorInterfaceExtensions.getiterator(x)
        return fromcolumns(Tables.columns(Tables.DataValueUnwrapper(y)))
    elseif it === missing
        if x isa AbstractVector && all(col -> isa(col, AbstractVector), x)
            return DataFrame(Vector{AbstractVector}(x))
        end
        y = IteratorInterfaceExtensions.getiterator(x)
        # non-NamedTuple or EltypeUnknown
        return fromcolumns(Tables.buildcolumns(nothing, Tables.DataValueUnwrapper(y)))
    end
    if x isa AbstractVector && all(col -> isa(col, AbstractVector), x)
        return DataFrame(Vector{AbstractVector}(x))
    end
    throw(ArgumentError("unable to construct DataFrame from $(typeof(x))"))
end

Base.append!(df::DataFrame, x) = append!(df, DataFrame(x))

# This supports the Tables.RowTable type; needed to avoid ambiguities w/ another constructor
DataFrame(x::Vector{T}) where {T <: NamedTuple} = fromcolumns(Tables.columns(Tables.DataValueUnwrapper(x)))

IteratorInterfaceExtensions.getiterator(df::DataFrame) = Tables.datavaluerows(df)
IteratorInterfaceExtensions.isiterable(x::DataFrame) = true
TableTraits.isiterabletable(x::DataFrame) = true


text:
# Autogenerated wrapper script for Libassuan_jll for aarch64-linux-musl
export libassuan

using Libgpg_error_jll
JLLWrappers.@generate_wrapper_header("Libassuan")
JLLWrappers.@declare_library_product(libassuan, "libassuan.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Libgpg_error_jll)
    JLLWrappers.@init_library_product(
        libassuan,
        "lib/libassuan.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using PtFEM

data = Dict(
  # Beam(ndim, nst, nxe, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Beam(2, 1, 4, 1, :x, Line(2, 1), false),
  :properties => [1.0;],
  :x_coords => [0.0, 0.25, 0.5, 0.75, 1.0],
  :support => [
    (1, [0 1]),
    (5, [0 0])
    ],
  :limit => 100,
  :tol => 0.00001
)

data |> display
println()

@time m = p46(data)
println()

println("\nThe buckling load = $(m[1])")
m[4][1] = 0.0
println("\nThe buckling mode (iterations=$(m[2])):\n")
for i in 1:m[3]
  println("$i    $(m[4][m[5][:, i]+1])")
end


text:
# logging

using Logging

# FIXME: replace with an additional log level
macro trace(ex...)
    esc(:(@debug $(ex...)))
end

# fatal versions to `@error`, including a safe version for compile-time errors
# (to be used instead of `error`, which _emits_ an error)
macro fatal(ex...)
    esc(quote
        @error $(ex...)
        error("Fatal error occurred")
    end)
end

# define safe loggers for use in generated functions (where task switches are not allowed)
for level in [:trace, :debug, :info, :warn, :error, :fatal]
    @eval begin
        macro $(Symbol("safe_$level"))(ex...)
            macrocall = :(@placeholder $(ex...))
            # NOTE: `@placeholder` in order to avoid hard-coding @__LINE__ etc
            macrocall.args[1] = Symbol($"@$level")
            quote
                old_logger = global_logger()
                global_logger(Logging.ConsoleLogger(Core.stderr, old_logger.min_level))
                ret = $(esc(macrocall))
                global_logger(old_logger)
                ret
            end
        end
    end
end


text:
"""
    Map <: AbstractDict{MOI.ObjectiveFunction, AbstractBridge}

A mapping from a bridged objective function type to the bridge responsible for
bridging that type of objective function.
"""
mutable struct Map <: AbstractDict{MOI.ObjectiveFunction, AbstractBridge}
    bridges::Dict{MOI.ObjectiveFunction, AbstractBridge}
    function_type::Union{Nothing, Type{<:MOI.AbstractScalarFunction}}
end
function Map()
    return Map(Dict{MOI.ObjectiveFunction, AbstractBridge}(), nothing)
end

# Implementation of `AbstractDict` interface.

Base.isempty(map::Map) = isempty(map.bridges)
function Base.empty!(map::Map)
    empty!(map.bridges)
    map.function_type = nothing
end
function Base.haskey(map::Map, attr::MOI.ObjectiveFunction)
    return haskey(map.bridges, attr)
end
function Base.getindex(map::Map, attr::MOI.ObjectiveFunction)
    return map.bridges[attr]
end
Base.length(map::Map) = length(map.bridges)
Base.values(map::Map) = values(map.bridges)
Base.iterate(map::Map, args...) = iterate(map.bridges, args...)

# Custom interface for information needed by `AbstractBridgeOptimizer`s that is
# not part of the `AbstractDict` interface.

"""
    function_type(map::Map)

Return the function type of the [`root_bridge`](@ref) or `nothing` if `map` is
empty.
"""
function_type(map::Map) = map.function_type

"""
    root_bridge(map::Map)

Return the last bridge added.
"""
function root_bridge(map::Map)
    attr = MOI.ObjectiveFunction{function_type(map)}()
    return map[attr]
end

"""
    add_key_for_bridge(map::Map, bridge::AbstractBridge,
                       func::MOI.AbstractScalarFunction)

Stores the mapping `atttr => bridge` where `attr` is
`MOI.ObjectiveFunction{typeof(func)}()` and set [`function_type`](@ref) to
`typeof(func)`.
"""
function add_key_for_bridge(map::Map, bridge::AbstractBridge,
                            func::MOI.AbstractScalarFunction)
    attr = MOI.ObjectiveFunction{typeof(func)}()
    map.function_type = typeof(func)
    map.bridges[attr] = bridge
end

"""
    EmptyMap <: AbstractDict{MOI.ObjectiveFunction, AbstractBridge}

Empty version of [`Map`](@ref). It is used by
[`MathOptInterface.Bridges.Variable.SingleBridgeOptimizer`](@ref) and
[`MathOptInterface.Bridges.Constraint.SingleBridgeOptimizer`](@ref) as they do
not bridge any objective function.
"""
struct EmptyMap <: AbstractDict{MOI.ObjectiveFunction, AbstractBridge} end
Base.isempty(::EmptyMap) = true
function Base.empty!(::EmptyMap) end
Base.length(::EmptyMap) = 0
Base.haskey(::EmptyMap, ::MOI.ObjectiveFunction) = false
Base.values(::EmptyMap) = MOIU.EmptyVector{AbstractBridge}()
Base.iterate(::EmptyMap) = nothing
function_type(::EmptyMap) = nothing


text:
"""
	$(SIGNATURES)

Save an optimization history.
Since serialization is flaky, it is safest to save each "series" in a separate file. If one file cannot be read, the rest of the history may still be readable.
Since binary saving does not work reliably, use JSON.
"""
function save_mvhistory(histS :: MVHistory, fPath :: AbstractString)
    # Save the series that are available
    sNameV = save_meta(histS, fPath);

    # Save each series
    for sName in sNameV
        # d = ValueHistoriesLH.history_to_dict(histS, sName);
        sPath = series_path(fPath, sName);
        save_history(get_series(histS, sName), sPath);
    end
    return nothing
end

"""
	$(SIGNATURES)

Load MVHistory.
"""
function load_mvhistory(fPath :: AbstractString; notFoundError :: Bool = true)
    sNameV = load_meta(fPath; notFoundError = notFoundError);
    if isnothing(sNameV)
        return nothing
    end

    histS = MVHistory();
    for sName in sNameV
        h = load_history(fPath, sName; notFoundError = notFoundError);
        if isnothing(h)
            if notFoundError
                error("Series $sName not found");
            end
        else
            add_series!(histS, sName, h);
        end
    end
    return histS :: MVHistory
end

function save_meta(h :: MVHistory, fPath :: AbstractString)
    sNameV = collect(keys(h));
    d = vector_to_dict(sNameV);
    metaPath = meta_path(fPath);
    save_dict(d, metaPath);
    return sNameV
end

function load_meta(fPath :: AbstractString; notFoundError :: Bool = true)
    metaPath = meta_path(fPath);
    if !isfile(metaPath)
        if notFoundError
            showPath = path_to_show(fPath);
            error("Not found: $showPath");
        else
            return nothing
        end
    end
    d = load_dict(metaPath);
    _, v = dict_to_vector(d);
    v = Symbol.(v);
    return v :: Vector{Symbol}
end

meta_path(fPath :: AbstractString) = series_path(fPath, :meta);
file_ext() = ".json3";

function series_path(fPath :: AbstractString, sName :: Symbol)
    fDir, fName = splitdir(fPath);
    baseName, fExt = splitext(fName);
    return joinpath(fDir,  baseName * "_$sName" * file_ext())
end

function path_to_show(fPath :: AbstractString)
    fDir, fName = splitdir(fPath);
    v = splitpath(fDir);
    n = length(v);
    if n < 3
        return fPath
    else
        idxV = (n-1) : n;
        return joinpath(v[idxV]..., fName)
    end
end

function vector_to_dict(vVec :: Vector)
	d = Dict{String, eltype(vVec)}();
	for j = 1 : length(vVec)
		d[string(j)] = vVec[j];
	end
	return d
end

## ------------  Save uni values history

"""
	$(SIGNATURES)

Save a uni-value history.
"""
function save_history(h :: History, fPath :: AbstractString;  showMsg :: Bool = true)
    d = ValueHistoriesLH.history_to_dict(h);
    save_dict(d, fPath);
    if showMsg
        showPath = path_to_show(fPath);
        println("Saved history $showPath");
    end
    return nothing
end

# Save a `Dict` to a JSON file
function save_dict(d, fPath :: AbstractString)
    open(fPath, "w") do io
        JSON3.write(io, d);
    end
    return nothing
end

function load_dict(fPath :: AbstractString)
    @assert isfile(fPath)
    d = try
        open(fPath, "r") do io
            # This loads the `Dict` that was saved
            JSON3.read(io);
        end
    catch
        error("Could not open $fPath")
    end
    return d
end

"""
	$(SIGNATURES)

Load a uni-value history.
Due to a JSON limitation, matrices are flattened into Vectors.
"""
function load_history(fPath)
    d = load_dict(fPath);
    idxV, valueV = dict_to_vector(d);
    h = make_history(idxV, valueV);
    return h
end

function load_history(fPath :: AbstractString, sName :: Symbol; notFoundError :: Bool = true)
    sp = series_path(fPath, sName);
    if !isfile(sp)
        if notFoundError
            showPath = path_to_show(sp);
            error("Not found: $showPath");
        else
            return nothing
        end
    else
        return load_history(sp);
    end
end

# Convert a loaded `Dict` into a `Vector`.
# Assumes keys are sequential integers. But stored as `Symbol{String}` (b/c of JSON3).
# JSON arrays are made into ordinary arrays.
function dict_to_vector(d :: T1) where T1 <: AbstractDict
    n = length(d);
    keyV = sort(key_to_int.(collect(keys(d))));
    @assert all(x -> x > 0, keyV)
    # Because each value can be stored as a different type, start with `Vector{Any}`.
    idxV = Vector{Int}(undef, n);
    v = Vector{Any}(undef, n);
    for (j, k) in enumerate(keyV)
        v[j] = d[int_to_key(k)];
    end

    # Promote to common type
    vOut = collect(promote(v...));
    # Make JSON3 arrays into ordinary arrays
    vOut = convert_json_array(vOut);
	return keyV, vOut
end

# Convert a JSON key to an Int
key_to_int(k) = parse(Int, string(k));

int_to_key(j :: Integer) = string(j);

# Find the largest key in a `Dict`
function largest_key(d :: T1) where T1 <: AbstractDict
    kMax = -100;
    for k in keys(d)
        kInt = key_to_int(k);
        @assert kInt > 0
        @assert isa(kInt, Integer)
        (kInt > kMax)  &&  (kMax = kInt);
    end
    @assert kMax > 1
    return kMax
end

# Make JSON array into ordinary array
function convert_json_array(x :: JSON3.Array; vType = Float64)
    y = similar(x, vType);
    for j in 1 : length(x)
        y[j] = x[j];
    end
    return y
end

# Make Array{JSON3.Array} into ordinary Arrays
function convert_json_array(x :: Array{T1}; vType = Float64) where T1 <: JSON3.Array
    y1 = convert_json_array(x[1]);
    y = similar(x, typeof(y1));
    for j = 1 : length(x)
        y[j] = convert_json_array(x[j]; vType = vType);
    end
    return y
end

convert_json_array(x) = x;

# -------------