text:
module mcompshelldyn0

# Simply supported square isotropic homogeneous plate.

using Arpack
using LinearAlgebra: norm, Transpose, mul!, I, Symmetric
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test_homogeneous()
    formul = FEMMShellT3FFModule

    ax = ay = 1000*phun("mm")
    nx = ny = 8
    E = 1000.0*phun("MPa")
    nu = 0.396
    rho = 2.0*phun("kg/m^3")
    thickness = ax/1000;
    tolerance = ax/nx/100
    D = E*thickness^3/12/(1-nu^2)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # No in plane displacements
    # l1 = collect(1:count(fens))
    # for i in [1, 2,]
    #     setebc!(dchi, l1, true, i)
    # end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1, 2, 3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = formul.mass(femm, geom0, dchi);

    # Solve
    neigvs = 8
    d, v, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    # @show nconv
    fs = real(sqrt.(complex(d))) / (2 * pi)
    # From Blevins, Table 5.3, square plate SSSS
    # for m in 1:3
    #     for n in 1:3
    #         @show pi^2 * (m^2 + n^2) * sqrt(D/rho/thickness) / (2 * pi * ax^2)
    #     end
    # end

    # vectors = []
    # for i in 1:neigvs
    #     scattersysvec!(dchi, v[:, i])
    #     push!(vectors, ("mode_$i", deepcopy(dchi.values[:, 1:3])))
    # end
    # vtkwrite("plate-modes.vtu", fens, fes; vectors = vectors)
    # @test maximum(dchi.values[:, 3]) ./phun("mm") ≈ 0.22004349767718365
    @test norm(fs - [21.822774909379287, 54.52203720717488, 54.60475772077202, 86.5480437899711, 108.65349048567792, 108.78050954575491, 138.29506871044424, 140.89370172063866]) < 1.0e-3 * norm(fs)
    true
end

function test_composite()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 1000*phun("mm")
    nx = ny = 8
    E = 1000.0*phun("MPa")
    nu = 0.396
    rho = 2.0*phun("kg/m^3")
    thickness = ax/1000;
    tolerance = ax/nx/100
    D = E*thickness^3/12/(1-nu^2)

    mater = CM.lamina_material(rho, E, nu)
    # Note that we model the isotropic homogeneous plate with multiple plies.
    plies = CM.Ply[]
    push!(plies, CM.Ply("ply_1", mater, thickness/4, 45))
    push!(plies, CM.Ply("ply_2", mater, thickness/4, -45))
    push!(plies, CM.Ply("ply_3", mater, thickness/4, 45))
    push!(plies, CM.Ply("ply_4", mater, thickness/4, -45))

    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("petyt_9.2", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Pin one of the corners
    l1 = selectnode(fens; box = Float64[ax ax 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,2, 3,]
        setebc!(dchi, l1, true, i)
    end
    # Roller at the other
    l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1, 2, 3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = formul.mass(femm, geom0, dchi);

    # Solve
    neigvs = 8
    d, v, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    # @show nconv
    fs = real(sqrt.(complex(d))) / (2 * pi)
    # oms = real(sqrt.(complex(d)))
    # @show (oms*ax^2*sqrt(rho/E/thickness^2))
    # vectors = []
    # for i in 1:neigvs
    #     scattersysvec!(dchi, v[:, i])
    #     push!(vectors, ("mode_$i", deepcopy(dchi.values[:, 1:3])))
    # end
    # vtkwrite("plate-modes.vtu", fens, fes; vectors = vectors)
    @test norm(fs - [21.822774909379287, 54.52203720717488, 54.60475772077202, 86.5480437899711, 108.65349048567792, 108.78050954575491, 138.29506871044424, 140.89370172063866]) < 1.0e-3 * norm(fs)
    true
end
end
using .mcompshelldyn0
mcompshelldyn0.test_homogeneous()
mcompshelldyn0.test_composite()

# module mcompshelldyn1

# # Introduction to Finite Element Vibration Analysis 2nd Edition
# # by Maurice Petyt
# #  Publisher: Cambridge University Press; 2nd edition (August 1, 2010)
# # ISBN-10: 0521191602
# # ISBN-13: 978-0521191609

# # Table 9.2. Non-dimensional natural frequencies of a simply supported
# # square laminated plate (a/h = 10) [45/−45/45/−45]
# # Natural frequencies, non dimensional
# # 18.46, 34.87, 54.27, 67.17, 75.58, 82.84

# using Arpack
# using LinearAlgebra: norm, Transpose, mul!, I, Symmetric
# using FinEtools
# using FinEtoolsDeforLinear
# using FinEtoolsFlexStructures.CompositeLayupModule
# using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
# using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
# using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
# using FinEtools.MeshExportModule.VTKWrite: vtkwrite
# using Test

# function test()
#     formul = FEMMShellT3FFCompModule
#     CM = CompositeLayupModule

#     ax = ay = 1000*phun("mm")
#     nx = ny = 50
#     E1 = 1000.0*phun("MPa")
#     E2 = E1/40
#     G12 = E2*0.6
#     G13 = E2*0.5
#     nu12 = 0.25;
#     nu23 = 0.396
#     G23 = E2*0.5
#     rho = 10.0*phun("kg/m^3")
#     thickness = ax/1000;
#     tolerance = ax/nx/100
#     CM = CompositeLayupModule

#     mater = CM.lamina_material(rho, E1, E2, nu12, G12, G13, G23)
#     plies = CM.Ply[]
#     push!(plies, CM.Ply("ply_1", mater, thickness/4, 45))
#     push!(plies, CM.Ply("ply_2", mater, thickness/4, -45))
#     push!(plies, CM.Ply("ply_3", mater, thickness/4, 45))
#     push!(plies, CM.Ply("ply_4", mater, thickness/4, -45))

#     mcsys = CM.cartesian_csys((1, 2, 3))
#     layup = CM.CompositeLayup("petyt_9.2", plies, mcsys)

#     fens, fes = T3block(ax,ay,nx,ny);
#     fens.xyz = xyz3(fens)
#     sfes = FESetShellT3()
#     accepttodelegate(fes, sfes)

#     femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

#     # Construct the requisite fields, geometry and displacement
#     # Initialize configuration variables
#     geom0 = NodalField(fens.xyz)
#     u0 = NodalField(zeros(size(fens.xyz,1), 3))
#     Rfield0 = initial_Rfield(fens)
#     dchi = NodalField(zeros(size(fens.xyz,1), 6))

#     # Apply EBC's
#     # Pin one of the corners
#     l1 = selectnode(fens; box = Float64[ax ax 0 0 -Inf Inf], inflate = tolerance)
#     for i in [1,2, 3,]
#         setebc!(dchi, l1, true, i)
#     end
#     # Roller at the other
#     l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
#     for i in [2,]
#         setebc!(dchi, l1, true, i)
#     end
#     # Simple support
#     l1 = connectednodes(meshboundary(fes))
#     for i in [1, 2, 3]
#         setebc!(dchi, l1, true, i)
#     end
#     applyebc!(dchi)
#     numberdofs!(dchi);

#     # Assemble the system matrix
#     formul.associategeometry!(femm, geom0)
#     K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);
#     M = formul.mass(femm, geom0, dchi);

#     # Solve
#     neigvs = 8
#     d, v, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
#     @show nconv
#     oms = real(sqrt.(complex(d)))
#     @show (oms*ax^2*sqrt(rho/E2/thickness^2))
#     vectors = []
#     for i in 1:neigvs
#         scattersysvec!(dchi, v[:, i])
#         push!(vectors, ("mode_$i", deepcopy(dchi.values[:, 1:3])))
#     end
#     vtkwrite("plate-modes.vtu", fens, fes; vectors = vectors)
#     # @test maximum(dchi.values[:, 3]) ./phun("mm") ≈ 0.22004349767718365
#     true
# end
# end
# using .mcompshelldyn1
# mcompshelldyn1.test()

module mcompshelldyn2

# Free vibration analysis of composite sandwich plates based on Reddy’s
# higher-order theory
# A.K. Nayak, S.S.J. Moy, R.A. Shenoi
# Section 4.3, Free vibration analysis of a simply supported square
# plate (orthotropic)

# From Table 2, nondimensional frequencies
# (1,1)    0.0474
# (1,2)    0.1033
# (2,1)    0.1188
# (2,2)    0.1694
# (1,3)    0.1888
# (3,1)    0.2180
# (2,3)    0.2475
# (3,2)    0.2624
# (1,4)    0.2969

using Arpack
using LinearAlgebra: norm, Transpose, mul!, I, Symmetric
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test(nplies = 10, axes = (1, 2, 3))
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 100*phun("mm")
    nx = ny = 9
    # aragonite crystals
    E1 = 143.52*phun("GPa")
    E2 = 75.38*phun("GPa")
    G12 = 42.03*phun("GPa")
    G13 = 25.56*phun("GPa")
    nu12 = 0.44;
    nu23 = 0.23
    G23 = 42.65*phun("GPa")
    rho = 1500.0*phun("kg/m^3") # only a guess
    C11 = 159.85*phun("GPa")
    thickness = ax/10;
    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(rho, E1, E2, nu12, G12, G13, G23)
    plies = [CM.Ply("ply_$i", mater, thickness/nplies, 0) for i in 1:nplies]

    mcsys = CM.cartesian_csys(axes)
    layup = CM.CompositeLayup("Nayak 4.3", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Pin one of the corners
    l1 = selectnode(fens; box = Float64[ax ax 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,2, 3,]
        setebc!(dchi, l1, true, i)
    end
    # Roller at the other
    l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1, 2, 3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = formul.mass(femm, geom0, dchi);

    # Solve
    neigvs = 9
    d, v, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    @test nconv == neigvs
    fs = real(sqrt.(complex(d))) / (2 * pi)

    @test norm([0.04571652264814249, 0.10096323319412286, 0.11467782865238098, 0.16264669289730896, 0.18683613787902217, 0.20655080934826295, 0.23954509827380233, 0.2476225295798577, 0.2718218525817157] - thickness * sqrt(rho/C11) .* (2 * pi * fs)) < 1.0e-13

    vectors = []
    for i in 1:neigvs
        scattersysvec!(dchi, v[:, i])
        push!(vectors, ("mode_$i", deepcopy(dchi.values[:, 1:3])))
    end
    vtkwrite("plate-modes.vtu", fens, fes; vectors = vectors)

    true
end
end
using .mcompshelldyn2
mcompshelldyn2.test()
mcompshelldyn2.test(3, (2, -1, 3))
mcompshelldyn2.test(5, (-1, -2, 3))
mcompshelldyn2.test(4, (-2, 1, 3))

module mcompshelldyn3

# Free vibration analysis of composite sandwich plates based on Reddy’s
# higher-order theory
# A.K. Nayak, S.S.J. Moy, R.A. Shenoi
# Section 4.4, Free vibration analysis of a simply supported cross ply
# (0/90/90/0) square plate (symmetric cross-ply)

# Better source of data:

# STABILITY AND VIBRATION OF ISOTROPIC, ORTHOTROPIC AND LAMINATED PLATES
# ACCORDING TO A HIGHER-ORDER SHEAR DEFORMATION THEORY J. N. REDDY AND N. D.
# PHAN
# TABLE 3  Non-dimensionalized fundamental frequencies of cross-ply square
# plates

using Arpack
using LinearAlgebra: norm, Transpose, mul!, I, Symmetric
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test(tL_ratio = 1/100, axes = (1, 2, 3))
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 100*phun("mm")
    nx = ny = 32
    # aragonite crystals
    E1 = 143.52*phun("GPa")
    E2 = E1/40
    G12 = E2*0.6
    G13 = E2*0.6
    nu12 = 0.25;
    G23 = E2*0.5
    rho = 1500.0*phun("kg/m^3") # only a guess
    thickness = ax*tL_ratio;
    nd_fundamental = Dict(1/5 => 10.989, 1/10 => 15.270, 1/100 => 18.755)
    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(rho, E1, E2, nu12, G12, G13, G23)
    nplies = 4
    plies = [
    CM.Ply("ply_1", mater, thickness/nplies, 0),
    CM.Ply("ply_2", mater, thickness/nplies, 90),
    CM.Ply("ply_3", mater, thickness/nplies, 90),
    CM.Ply("ply_4", mater, thickness/nplies, 0), ]

    mcsys = CM.cartesian_csys(axes)
    layup = CM.CompositeLayup("Nayak 4.4", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Pin one of the corners
    l1 = selectnode(fens; box = Float64[ax ax 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,2, 3,]
        setebc!(dchi, l1, true, i)
    end
    # Roller at the other
    l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1, 2, 3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = formul.mass(femm, geom0, dchi);

    # Solve
    neigvs = 9
    d, v, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    @test nconv == neigvs
    fs = real(sqrt.(complex(d))) / (2 * pi)

    # @show nd_fundamental[tL_ratio]
    # @show 2*pi*fs[1]*ax^2/thickness*sqrt(rho/E2)
    @test abs(nd_fundamental[tL_ratio] - 2*pi*fs[1]*ax^2/thickness*sqrt(rho/E2)) / nd_fundamental[tL_ratio] < 3.0e-2
    # vectors = []
    # for i in 1:neigvs
    #     scattersysvec!(dchi, v[:, i])
    #     push!(vectors, ("mode_$i", deepcopy(dchi.values[:, 1:3])))
    # end
    # vtkwrite("plate-modes.vtu", fens, fes; vectors = vectors)

    true
end
end
using .mcompshelldyn3
mcompshelldyn3.test(1/5, (1, 2, 3))
mcompshelldyn3.test(1/10, (1, 2, 3))
mcompshelldyn3.test(1/100, (1, 2, 3))
mcompshelldyn3.test(1/100, (2, -1, 3))
mcompshelldyn3.test(1/100, (-2, 1, 3))

module mcompshelldyn4

# Journal of Sound and Vibration (1985) 98(2), 157-170 STABILITY AND VIBRATION
# OF ISOTROPIC, ORTHOTROPIC AND LAMINATED PLATES ACCORDING TO A HIGHER-ORDER
# SHEAR DEFORMATION THEORY J. N. REDDY AND N. D. PHAN

# TABLE 3  Non-dimensionalized fundamental frequencies of cross-ply square
# plates

using Arpack
using LinearAlgebra: norm, Transpose, mul!, I, Symmetric
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test(tL_ratio = 1/100, axes = (1, 2, 3))
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    a = 100*phun("mm")
    nx = ny = 8
    # aragonite crystals
    E1 = 143.52*phun("GPa")
    E2 = E1/40
    G12 = E2*0.6
    G13 = E2*0.6
    nu12 = 0.25;
    G23 = E2*0.5
    rho = 1500.0*phun("kg/m^3") # only a guess
    thickness = a*tL_ratio;
    nd_fundamental = Dict(1/5 => 9.01, 1/10 => 10.449, 1/100 => 11.156)
    approx_nd_fundamental = Dict(1/5 => 8.754289808606417, 1/10 => 10.35796013315713, 1/100 => 11.160996567154124)
    tolerance = a/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(rho, E1, E2, nu12, G12, G13, G23)
    nplies = 2
    plies = [
    CM.Ply("ply_1", mater, thickness/nplies, 0),
    CM.Ply("ply_2", mater, thickness/nplies, 90),  ]

    mcsys = CM.cartesian_csys(axes)
    layup = CM.CompositeLayup("Nayak 4.4", plies, mcsys)

    fens, fes = T3block(a,a,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), CM.thickness(layup)), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Simple support
        # The boundary conditions are a bit peculiar: refer to the paper
        l1 = selectnode(fens; box = Float64[0 0 -Inf Inf 0 0], inflate = tolerance)
        for i in [2, 4]
            setebc!(dchi, l1, true, i)
        end
        l1 = selectnode(fens; box = Float64[a a -Inf Inf 0 0], inflate = tolerance)
        for i in [2, 4]
            setebc!(dchi, l1, true, i)
        end
        l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 0 0], inflate = tolerance)
        for i in [1, 5]
            setebc!(dchi, l1, true, i)
        end
        l1 = selectnode(fens; box = Float64[-Inf Inf a a 0 0], inflate = tolerance)
        for i in [1, 5]
            setebc!(dchi, l1, true, i)
        end
        # Simple support
        l1 = connectednodes(meshboundary(fes))
        for i in [3]
            setebc!(dchi, l1, true, i)
        end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = formul.mass(femm, geom0, dchi);

    # Solve
    neigvs = 9
    d, v, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    @test nconv == neigvs
    fs = real(sqrt.(complex(d))) / (2 * pi)

    # @show nd_fundamental[tL_ratio]
    # @show 2*pi*fs[1]*a^2/thickness*sqrt(rho/E2)
    @test 2*pi*fs[1]*a^2/thickness*sqrt(rho/E2) ≈ approx_nd_fundamental[tL_ratio]
    # @test abs(nd_fundamental[tL_ratio] - 2*pi*fs[1]*ax^2/thickness*sqrt(rho/E2)) / nd_fundamental[tL_ratio] < 3.0e-2
    # vectors = []
    # for i in 1:neigvs
    #     scattersysvec!(dchi, v[:, i])
    #     push!(vectors, ("mode_$i", deepcopy(dchi.values[:, 1:3])))
    # end
    # vtkwrite("plate-modes.vtu", fens, fes; vectors = vectors)

    true
end
end
using .mcompshelldyn4
mcompshelldyn4.test(1/5, (1, 2, 3))
mcompshelldyn4.test(1/10, (1, 2, 3))
mcompshelldyn4.test(1/100, (1, 2, 3))
mcompshelldyn4.test(1/100, (2, -1, 3))
mcompshelldyn4.test(1/100, (-2, 1, 3))

module mcompshelldyn6
# Similar definition of the plate as in:
# From Barbero's Finite Element Analysis using Abaqus ... book Example 3.1
# Free vibration problem, fundamental frequency with Abaqus 42.656 Hz
using Arpack
using LinearAlgebra: norm, Transpose, mul!, I, Symmetric
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 1000*phun("mm") # note the changed dimension
    nx = ny = 50
    # ASFD/9310
    rho = 1500*phun("kg/m^3")
    E1 = 133860*phun("MPa")
    E2 = 7706*phun("MPa")
    G12 = 4306*phun("MPa")
    G13 = G12
    nu12 = 0.301;
    nu23 = 0.396
    G23 = 2760*phun("MPa")

    npairs = 1
    thickness = 10*phun("mm");
    # With these inputs, the Abaqus-verified solution is 42.656 Hz

    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(rho, E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[]
    for p in 1:npairs
        push!(plies, CM.Ply("ply_0_$p", mater, thickness/npairs/2, 0))
        push!(plies, CM.Ply("ply_90_$p", mater, thickness/npairs/2, 90))
    end
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("example_3.1", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1, 2, 3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = formul.mass(femm, geom0, dchi);

        # Solve
    neigvs = 9
    d, v, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    @test nconv == neigvs
    fs = real(sqrt.(complex(d))) / (2 * pi)
    @test abs(fs[1] - 42.62) / 42.62 < 1.0e-2
    # vtkwrite("plate-$npairs-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    # @test minimum(dchi.values[:, 3]) ./phun("mm") ≈ -1502.0324040447633
    true
end
end
using .mcompshelldyn6
mcompshelldyn6.test()

module mcompshelldyn7
# Similar definition of the plate as in mcompshell6
# but the layup is [-45/45]
# Free vibration problem, fundamental frequency with Abaqus 42.656 Hz
using Arpack
using LinearAlgebra: norm, Transpose, mul!, I, Symmetric
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test(axes = (1, 2, 3))
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 1000*phun("mm") # note the changed dimension
    nx = ny = 50
    # ASFD/9310
    rho = 1500*phun("kg/m^3")
    E1 = 133860*phun("MPa")
    E2 = 7706*phun("MPa")
    G12 = 4306*phun("MPa")
    G13 = G12
    nu12 = 0.301;
    nu23 = 0.396
    G23 = 2760*phun("MPa")

    npairs = 1
    thickness = 10*phun("mm");
    # With these inputs, the Abaqus-verified solution is 47.03 Hz

    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(rho, E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[]
    for p in 1:npairs
        push!(plies, CM.Ply("ply_-45_$p", mater, thickness/npairs/2, -45))
        push!(plies, CM.Ply("ply_+45_$p", mater, thickness/npairs/2, +45))
    end
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("example_3.1", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1, 2, 3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = formul.mass(femm, geom0, dchi);

        # Solve
    neigvs = 9
    d, v, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, explicittransform=:none)
    @test nconv == neigvs
    fs = real(sqrt.(complex(d))) / (2 * pi)
    @test abs(fs[1] - 47.86476186783638) / 47.86476186783638 < 1.0e-2
    # vtkwrite("plate-$npairs-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    # @test minimum(dchi.values[:, 3]) ./phun("mm") ≈ -1502.0324040447633
    true
end
end
using .mcompshelldyn7
mcompshelldyn7.test()
mcompshelldyn7.test((-2, 1, 3))


text:
function minimg(p1 , p2 , side)

    delta_x = (p1[1] - p2[1]) % side
    if delta_x > side/2
        delta_x -= side
    elseif delta_x < -side/2
        delta_x += side
    end

    delta_y = (p1[2] - p2[2]) % side
    if delta_y > side/2
        delta_y -= side
    elseif delta_y < -side/2
        delta_y += side
    end
    minimg = sqrt(delta_x^2+delta_y^2)
    return (minimg, delta_x, delta_y)
end

ext: jl


text:
# Utils
# #+HTML: Source
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Utils][Utils:1]]
export ⊗, realize

abstract type Direction end

struct Left  <: Direction end # Often useful to dispatch on direction an algorithm is going
struct Right <: Direction end

const left  = Left()
const right = Right()

A ⊗ B = kron(A, B)

realize(x::Number) = error("Unrecognized numerical type")
realize(x::Real) = x
function realize(x::Complex; ϵ=1e-10)
    abs(imag(x)) < ϵ || error("Non-zero imaginary component, $(imag(x))")
    real(x)
end

dg(M::Array{T, 4}) where {T} = permutedims(conj.(M), (2, 1, 3, 4))
dg(M::Array{T, 3}) where {T} = permutedims(conj.(M), (2, 1, 3))

not(x) = ~x
# Utils:1 ends here


text:
### Fetch required packages and reaction networks ###
using Catalyst, Random, Statistics, StochasticDiffEq, Test
include("test_networks.jl")

### Compares to the manually calcualted function ###
identical_networks = Vector{Pair}()

function real_f_1(du,u,p,t)
    X1,X2,X3 = u
    p,k1,k2,k3,d = p
    du[1] = 2*p - k1*X1
    du[2] = k1*X1 - k2*X2 - k3*X2
    du[3] = k2*X2 + k3*X2 - d*X3
end
function real_g_1(du,u,p,t)
    X1,X2,X3 = u
    p,k1,k2,k3,d = p
    du[1,1] = 2*sqrt(p)
    du[1,2] = -sqrt(k1*X1)
    du[1,3] = 0
    du[1,4] = 0
    du[1,5] = 0
    du[2,1] = 0
    du[2,2] = sqrt(k1*X1)
    du[2,3] = -sqrt(k2*X2)
    du[2,4] = -sqrt(k3*X2)
    du[2,5] = 0
    du[3,1] = 0
    du[3,2] = 0
    du[3,3] = sqrt(k2*X2)
    du[3,4] = sqrt(k3*X2)
    du[3,5] = -sqrt(d*X3)
end
push!(identical_networks, reaction_networks_standard[8] => (real_f_1,real_g_1,zeros(3,5)))

function real_f_2(du,u,p,t)
    X1, = u
    v,K,n,d = p
    du[1] = v/10+v*X1^n/(X1^n+K^n) - d*X1
end
function real_g_2(du,u,p,t)
    X1, = u
    v,K,n,d = p
    du[1,1] = sqrt(v/10+v*X1^n/(X1^n+K^n))
    du[1,2] = -sqrt(d*X1)
end
push!(identical_networks, reaction_networks_hill[6] => (real_f_2,real_g_2,zeros(1,2)))

function real_f_3(du,u,p,t)
    X1,X2,X3,X4,X5,X6,X7 = u
    k1,k2,k3,k4,k5,k6 = p
    du[1] = -k1*X1*X2 + k2*X3
    du[2] = -k1*X1*X2 + k2*X3
    du[3] = k1*X1*X2 - k2*X3 - k3*X3*X4 + k4*X5
    du[4] = -k3*X3*X4 + k4*X5
    du[5] = k3*X3*X4 - k4*X5 - k5*X5*X6 + k6*X7
    du[6] = -k5*X5*X6 + k6*X7
    du[7] = k5*X5*X6 - k6*X7
end
function real_g_3(du,u,p,t)
    X1,X2,X3,X4,X5,X6,X7 = u
    k1,k2,k3,k4,k5,k6 = p
    du[1,1] = -sqrt(k1*X1*X2)
    du[1,2] = sqrt(k2*X3)
    du[1,3] = 0; du[1,4] = 0; du[1,5] = 0; du[1,6] = 0;
    du[2,1] = -sqrt(k1*X1*X2)
    du[2,2] = sqrt(k2*X3)
    du[2,3] = 0; du[2,4] = 0; du[2,5] = 0; du[2,6] = 0;
    du[3,1] = sqrt(k1*X1*X2)
    du[3,2] = -sqrt(k2*X3)
    du[3,3] = -sqrt(k3*X3*X4)
    du[3,4] = sqrt(k4*X5)
    du[3,5] = 0; du[3,6] = 0;
    du[4,3] = -sqrt(k3*X3*X4)
    du[4,4] = sqrt(k4*X5)
    du[4,1] = 0; du[4,2] = 0; du[4,5] = 0; du[4,6] = 0;
    du[5,3] = sqrt(k3*X3*X4)
    du[5,4] = -sqrt(k4*X5)
    du[5,5] = -sqrt(k5*X5*X6)
    du[5,6] = sqrt(k6*X7)
    du[5,1] = 0; du[5,2] = 0;
    du[6,5] = -sqrt(k5*X5*X6)
    du[6,6] = sqrt(k6*X7)
    du[6,1] = 0; du[6,2] = 0; du[6,3] = 0; du[6,4] = 0;
    du[7,5] = sqrt(k5*X5*X6)
    du[7,6] = -sqrt(k6*X7)
    du[7,1] = 0; du[7,2] = 0; du[7,3] = 0; du[7,4] = 0;
end
push!(identical_networks, reaction_networks_constraint[9] => (real_f_3,real_g_3,zeros(7,6)))

for (i,networks) in enumerate(identical_networks)
    for factor in [1e-1, 1e0, 1e1], repeat in 1:3
        u0 = 100. .+ factor*rand(length(networks[1].states))
        p = 0.01 .+ factor*rand(length(networks[1].ps))
        (i==2) && (u0[1] += 1000.)
        (i==3) ? (p[2:2:6] .*= 1000.; u0 .+= 1000) : (p[1] += 500.)
        prob1 = SDEProblem(networks[1],u0,(0.,100.),p)
        sol1 = solve(prob1,ImplicitEM(),saveat=0.01,maxiters=1e7)
        prob2 = SDEProblem(networks[2][1],networks[2][2],u0,(0.,100.),p,noise_rate_prototype=networks[2][3])
        sol2 = solve(prob2,ImplicitEM(),saveat=0.01,maxiters=1e7)
        for i = 1:length(u0)
            vals1 = getindex.(sol1.u[1000:end],i);
            vals2 = getindex.(sol1.u[1000:end],i);
            @test 0.8 < mean(vals1)/mean(vals2) < 1.25
            @test 0.8 < std(vals1)/std(vals2) < 1.25
        end
    end
end

### Compares level of noise with noise scalling. ###

# Tests with a single noise scaling parameter.
noise_scaling_network = @reaction_network begin
    (k1,k2), X1 ↔ X2
end k1 k2
for repeat = 1:5
    p = 1. .+ rand(2)
    u0 = 10000*(1. .+ rand(2))
    sol001 = solve(SDEProblem(noise_scaling_network,u0,(0.,1000.),vcat(p,0.01),noise_scaling=(@variables η1)[1]),ImplicitEM())
    sol01 = solve(SDEProblem(noise_scaling_network,u0,(0.,1000.),vcat(p,0.1),noise_scaling=(@variables η1)[1]),ImplicitEM())
    sol1 = solve(SDEProblem(noise_scaling_network,u0,(0.,1000.),vcat(p,1.),noise_scaling=(@variables η2)[1]),ImplicitEM())
    sol10 = solve(SDEProblem(noise_scaling_network,u0,(0.,1000.),vcat(p,10.),noise_scaling=(@variables η3)[1]),ImplicitEM())
    @test 2*std(first.(sol001.u)[100:end]) < std(first.(sol01.u)[100:end])
    @test 2*std(last.(sol001.u)[100:end]) < std(last.(sol01.u)[100:end])
    @test 2*std(first.(sol01.u)[100:end]) < std(first.(sol1.u)[100:end])
    @test 2*std(last.(sol01.u)[100:end]) < std(last.(sol1.u)[100:end])
    @test 2*std(first.(sol1.u)[100:end]) < std(first.(sol10.u)[100:end])
    @test 2*std(last.(sol1.u)[100:end]) < std(last.(sol10.u)[100:end])
end

### Tries to create a large number of problem, ensuring there are no errors (cannot solve as solution likely to go into negatives). ###
for reaction_network in reaction_networks_all
    for factor in [1e-2, 1e-1, 1e0, 1e1]
        u0 = factor*rand(length(reaction_network.states))
        p = factor*rand(length(reaction_network.ps))
        prob = SDEProblem(reaction_network,u0,(0.,1.),p)
    end
end

### No parameter test ###

no_param_network = @reaction_network begin
    (1.2,5), X1 ↔ X2
end
for factor in [1e3, 1e4]
    u0 = factor*(1. .+ rand(length(no_param_network.states)))
    prob = SDEProblem(no_param_network,u0,(0.,1000.))
    sol = solve(prob,ImplicitEM())
    vals1 = getindex.(sol.u[1:end],1)
    vals2 = getindex.(sol.u[1:end],2)
    @test mean(vals1) > mean(vals2)
end


text:
using SobolSearches, Test

@testset "simple consistency checks" begin
    y = [0.2, 0.7]
    f(x) = sum(abs2, x .- y)
    s = sobol_search(f, [-2, -1.0], [3, 1.0], 10, 10000)
    @test s[1].value == f(s[1].position)
    @test s[1].position ≈ y atol = 0.02
    @test IndexStyle(typeof(s)) ≡ IndexLinear()
    @test size(s) == (10, )
    @test length(s) == 10
end


text:
using DataValues

A = [DataValue(9), DataValue(8), DataValue(15)]

map(i -> isna(i) ? false : get(i) % 3 == 0, A)

f(i) = isna(i) ? false : get(i) % 3 == 0

f.(A)


text:
if basename(pwd()) == "aoc"
    cd("2019/5")
end
include("../IntCodeMachine/src/IntCodeMachine.jl")

part1(filename::AbstractString) = IntCodeMachine.execute!(IntCodeMachine.load(filename))
part1("input.txt")
# type 1 for part 1, 5 for part 2

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
#####
##### PMEM initialization
#####

"""
$(SIGNATURES)

Configure the PMM allocator in ngraph to allocated PMM files in `dir`.
Default: `/mnt/public`
"""
function setup_pmem(dir = "/mnt/public/")
    if isdir(dir)
        for file in readdir(dir)
            rm(joinpath(dir, file); recursive = true)
        end
    end

    manager = nGraph.Lib.getinstance()
    nGraph.Lib.set_pool_dir(manager, dir)
    return nothing
end

#####
##### Setup Affinities
#####

# For whatever reason, we have to set up the number of threads before the first
# compilation of an nGraph executable.
#
# After that, changes in the number of threads won't be seen by nGraph

"""
$(SIGNATURES)

Set OMP thread affinities be setting appropriate environmental variables.

* `omp_num_threads`: The number of threads to reserve. Default: `24`
* `reserved_cores`: The number of physical CPU cores to reserve. Default: `omp_num_threads`
* `threads_per_core`: If `1`, essentially disables hyperthreading. If `2`, enables
    hyperthreading. Default: `1`.

The effects of this function are only valid if the ngraph JIT compiler has not yet been
invoked. Afterwards, no changes will take place.

This method is called during Module initialization, but may be called as many times as
desired before the first ngraph compilation.

See also, [`teardown_affinities`](@ref)
"""
function setup_affinities(;
        omp_num_threads = 24,
        reserved_cores = omp_num_threads,
        threads_per_core = 1
    )

    if !in(threads_per_core, (1, 2))
        throw(ArgumentError("""
            Expected key word `threads_per_core` to be either 1 or 2.
            Instead, got $threads_per_core.
            """
        ))
    end

    if nGraph.have_compiled() == true
        @error """
        The nGraph compiler has already been invoked.

        The number of threads WILL NOT change.
        """
        return nothing
    end

    ENV["KMP_AFFINITY"] = "compact,granularity=fine"

    # Use the first 24 cores - 1 threads for each core
    # Send to numa-node 1 for a hopefully more quiet system
    #
    # See docs/src/runner/kmp.md for syntax documentation
    ENV["KMP_HW_SUBSET"] = "1s@1,$(reserved_cores)c,$(threads_per_core)t"

    # 1 Threads for each core
    ENV["OMP_NUM_THREADS"] = omp_num_threads
    ENV["OMP_MAX_ACTIVE_LEVELS"] = 2

    return nothing
end

_env_context() = (
    kmp_hw_subset = ENV["KMP_HW_SUBSET"]::String,
    omp_num_threads = ENV["OMP_NUM_THREADS"]::String,
)

"""
$(SIGNATURES)

Turn on ngraph's JIT feature.
"""
setup_codegen() = nGraph.enable_codegen()

"""
$(SIGNATURES)

Enable the ngraph compiler to reuse memory for intermediate tensors based on liveness
analysis.
"""
setup_passes() = nGraph.set_pass_attributes(nGraph.ReuseMemory())

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
"""
    causality(source::AbstractVector, target::AbstractVector, test::CausalityTest)

Test for a causal influence from `source` to `target` using the provided causality `test`.

## Examples

```julia
x, y = rand(300), rand(300)

# Define some causality tests and apply them to `x` and `y`.
test_ccm = ConvergentCrossMappingTest(timeseries_lengths = [45, 50], n_reps = 20)
test_cm = CrossMappingTest(n_reps = 10)
test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)
test_tog = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
test_nn = NearestNeighbourMITest(ηs = 1:5)
test_jdd = JointDistanceDistributionTest()
test_jddt = JointDistanceDistributionTTest()
test_pa = PredictiveAsymmetryTest(
    VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5))
test_pan = NormalisedPredictiveAsymmetryTest(f = 1.0,
    VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5))

causality(x, y, test_ccm)
causality(x, y, test_cm)
causality(x, y, test_vf)
causality(x, y, test_nn)
causality(x, y, test_tog)
causality(x, y, test_jdd)
causality(x, y, test_jddt)
causality(x, y, test_pa)
causality(x, y, test_pan)
```
"""
function causality(source::AbstractVector, target::AbstractVector, test::CausalityTest) end

"""
    causality(source::AbstractVector, target::AbstractVector, cond::AbstractVector,
        test::CausalityTest)

Test for a causal influence from `source` to `target` conditioned on `cond`,
using the provided causality `test`.

Conditional analysis is applicable to the following tests:

- [`VisitationFrequencyTest`](@ref)
- [`TransferOperatorGridTest`](@ref)
- [`PredictiveAsymmetryTest`](@ref)
- [`NormalisedPredictiveAsymmetryTest`](@ref)

## Examples

```julia
x, y, z = rand(300), rand(300), rand(300)

# Transfer entropy from `x` to `y`, conditioned on `z`. Return one
# transfer entropy value per prediction lag `η`.
test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)

causality(x, y, z, test_vf)
```
"""
function causality(source::AbstractVector, target::AbstractVector, cond::AbstractVector, test::CausalityTest) end

function summarise(test::CausalityTest)
    _type = typeof(test)

    strs = ["$fn = $(test.:($fn))" for fn in fieldnames(_type)]
    return "$_type" * "(" * join(strs, ", ") * ")"
end

Base.show(io::IO, test::CausalityTest) = print(io, summarise(test))


text:
function dcs_greedy{T}(G::Array{SparseMatrixCSC{T, Int64}, 1})

  for i=n:-1:2

  end
end

function charikar_greedy{T}(G::SparseMatrixCSC{T, Int64})
  n = size(G, 1)
  deg = sum()

end


text:
#This REPL requires:
# - a Slack Outgoing WebHooks integration, and
# - a Slack Incoming WebHooks integration

#Print debugging information to console?
DEBUG = false

#Replace this with your Slack token for your outgoing webhook
TOKEN = ""

#Replace this with your Webhook URL for your incoming webhook
INHOOK = "https://hooks.slack.com/services/"

#Add custom JSON entries to the returning payload to the incoming webhook
DEFAULTPAYLOAD = Dict()

#Load data from config.jl if present
#Search for config.jl in current working directory
configfile = joinpath(pwd(), "config.jl")
if isfile(configfile)
    include(configfile)
    DEBUG && info(string("Loading configuration from ", configfile))
end

isa(TOKEN, String) && (TOKEN = [TOKEN])

KOALA = """
　　　　　　　::.;;|
　　　　　　　::;;;|
　　　 r´""'''ヽー--､,,,,,.....__
　　　（:::::,,,..::;　　　 　　 .,,;;;:,..ヽ
　　　　ヾ,,.　　　　　　　 ';;;..:.:;; ）
　　　　　　| ｀ﾟ　l⌒l　ﾟ´ i.,,,,.ノ⌒"''''⌒ヽ
　　　　　　ト 　 l.,,,,.l　 ノY　　　　　 ,;;;;, ）
___　,ｨ,ｒ,ｧ　! ｀;ｰ‐--イ 　|　 ﾟ.ｒ‐t. ･ ｒ,ノ
ト ｀"　ｲ-┴''´""　;;;（"^''ヽ, ゝｲ　.ノ
ゝー-､,,　　　　　　　 ヾ､　　""'⌒ヽ､
　　　　 "''ゝ　　　　　　　ヽ　　 　 　 |
　　　　　 /'""　　　　　 　 |　　　　 　|
　　　　　|　　　　　　 　 r ､(　　　　　.j
　　　　<ヽi　　　　　 　 ゝ　　　　　 ノ
　　　　ゝ　　　　　 　　　 しイ⌒'""
　　　　　ヽ,,,ノ;;人,,,,,__,,......ノ
　　　　　　::;;;:|
　　　　　　:::;;;;|
"""
#####################

using Requests
using Morsel

app = Morsel.app()

#For debugging purposes, add an escape hatch to terminate current Julia
#instance if its internal state gets too hairy
DEBUG && route(app, GET | POST | PUT, "/reboot") do req, res
    exit(1)
end

route(app, GET | POST | PUT, "/") do req, res
    mycmd = ""
    channelname = "general"
    output, color, username = try
        DEBUG && info(string("Received HTTP packet with data\n", repr(req.state[:data]), "\n"))
        data = req.state[:data]

        (haskey(data, "token") && any(data["token"] .== TOKEN)) || error("Invalid Slack token")
        channelname = bytestring(data["channel_name"])
        username = bytestring(data["user_name"])

        if haskey(data, "text")
            cmdstart = if haskey(data, "trigger_word")
                length(data["trigger_word"])+2
            else#if haskey(data, "command")
                1
	        end
            mycmd = bytestring(strip(data["text"][cmdstart:end]))

            if length(mycmd)>=6 && mycmd[1:3] == mycmd[end-2:end] == "```"
                #Strip triple backquotes from start and end if present
                DEBUG && info("Removing triple backquotes")
                mycmd = mycmd[4:end-3]
            elseif length(mycmd)>=2 && mycmd[1] == mycmd[end] == '`'
                #Strip single backquotes from start and end if present
                DEBUG && info("Removing single backquotes")
                mycmd = mycmd[2:end-1]
            end

            #Unescape HTML entities produced by Slack
            mycmd = replace(mycmd, ">", ">")
            mycmd = replace(mycmd, "<", "<")
            mycmd = replace(mycmd, "&", "&")

            #Replace smart quotes
            mycmd = replace(mycmd, "â\u80\u009c", "\"")
            mycmd = replace(mycmd, "â\u80\u9d", "\"")
            mycmd = replace(mycmd, "â\u80\u0098", "'")
            mycmd = replace(mycmd, "â\u80\u99", "'")

            status = "good"
            cmdout = ""
            try
                cmdout = string(eval(parse(mycmd)))
            catch exc
                io = IOBuffer()
                Base.showerror(io, exc, catch_backtrace())
                status = "danger"
                cmdout = takebuf_string(io)
            end

            mycmd == "koala" && (cmdout = KOALA)
            (cmdout, status, username)
        else
            ("Could not recognize input", "danger", username)
        end
    catch exc
        io = IOBuffer()
        Base.showerror(io, exc, catch_backtrace())
        (takebuf_string(io), "danger")
    end

    payload = merge(DEFAULTPAYLOAD, Dict("channel"=>"#"*channelname,
           "attachments" => [
              Dict("title"=>"Julia input from "*username, "text"=>mycmd, "fallback"=>mycmd),
              Dict("title"=>"Julia output", "color"=>color, "text"=>output, "fallback"=>output)
    ]))

    DEBUG && info(string("Posting to URL ", INHOOK, " the JSON payload:\n", payload))
    Requests.post(INHOOK; json=payload)
end

start(app, 8000)


text:
include("ArrayFire/AF.jl")

module GenEvo

using Reexport
@reexport using AF

import AF: release!, copy, eval!

include("utils.jl")
include("genome/index.jl")
include("algorithms/index.jl")
include("networks/index.jl")
include("data/index.jl");

end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using Documenter, BilevelJuMP

makedocs(
    modules = [BilevelJuMP],
    doctest  = false,
    clean    = true,
    format   = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename = "BilevelJuMP.jl",
    authors  = "Joaquim Garcia",
    pages   = [
        "Home" => "index.md",
        "manual.md",
        "examples.md"
    ]
)

deploydocs(
    repo = "github.com/joaquimg/BilevelJuMP.jl.git",
)


text:
module Profiler

using Profile, JSON, FlameGraphs

import ..Atom:  expandpath, @msg, handle

function tojson(node::FlameGraphs.Node, root = false)
  name, path = expandpath(string(node.data.sf.file))
  classes = []

  node.data.status & 0x01 ≠ 0 && push!(classes, "dynamic-dispatch")
  node.data.status & 0x02 ≠ 0 && push!(classes, "garbage-collection")

  Dict(
    :path => path,
    :location => name,
    :func => node.data.sf.func,
    :line => node.data.sf.line,
    :count => root ? sum(length(c.data.span) for c in node) : length(node.data.span),
    :classes => classes,
    :children => [tojson(c) for c in node]
  )
end

function profiler(data = Profile.fetch() ;kwargs...)
  graph = FlameGraphs.flamegraph(data ;kwargs...)
  graph === nothing && return
  pruneinternal!(graph)
  prunetask!(graph)
  @msg profile(tojson(graph, true))
end

function prunetask!(node)
  for c in node
    if (
          (
            c.data.sf.func == :task_done_hook &&
            endswith(string(c.data.sf.file), "task.jl")
          ) ||
          (
            # REPL
            c.data.sf.func == :eval &&
            FlameGraphs.isleaf(c) &&
            endswith(string(c.data.sf.file), "boot.jl")
          )
       )
      FlameGraphs.prunebranch!(c)
    end
  end
end

function pruneinternal!(node)
  for c in node
    if (
        (
          # REPL evaluation
          c.data.sf.func == :eval &&
          c.parent.data.sf.func == :repleval &&
          endswith(string(c.parent.data.sf.file), "repl.jl")
        ) ||
        (
          # inline evaluation
          c.data.sf.func == :include_string &&
          c.parent.data.sf.func == :include_string &&
          endswith(string(c.parent.parent.data.sf.file), "eval.jl")
        )
       )

      # Add children directly to the root node
      root = node
      replchild = c
      while !FlameGraphs.isroot(root)
          root = root.parent
          replchild = replchild.parent
      end
      FlameGraphs.graftchildren!(root, c)
      # Eliminate all nodes in between. This might include some that don't
      # call REPL code, but as this is also internal it seems OK.
      FlameGraphs.prunebranch!(replchild)
      return true
    else
      pruneinternal!(c) && return true
    end
  end
  return false
end

handle("loadProfileTrace") do path
  if path === nothing
    return
  end

  json = try
    JSON.parse(String(open(read, first(path))))
  catch e
    @error "Error reading profile trace file at $path."
    nothing
  end
  json !== nothing && @msg profile(json)
end

handle("saveProfileTrace") do path, data
  if path === nothing
    return
  end

  try
    write(path, JSON.json(data))
  catch e
    @error "Error writing profile trace file at $path."
  end
end

end


text:
"""
    OceanGyre <: AbstractSimpleBoxProblem

Container structure for a simple box problem with wind-stress, coriolis force, and temperature forcing.
Lˣ = zonal (east-west) length
Lʸ = meridional (north-south) length
H  = height of the ocean
τₒ = maximum value of wind-stress (amplitude)
λʳ = temperature relaxation penetration constant (meters / second)
θᴱ = maximum surface temperature
"""
struct OceanGyre{T, BC} <: AbstractSimpleBoxProblem
    Lˣ::T
    Lʸ::T
    H::T
    τₒ::T
    λʳ::T
    θᴱ::T
    boundary_conditions::BC
    function OceanGyre{FT}(
        Lˣ,                  # m
        Lʸ,                  # m
        H;                   # m
        τₒ = FT(1e-1),       # N/m²
        λʳ = FT(4 // 86400), # m/s
        θᴱ = FT(10),         # K
        BC = (
            OceanBC(Impenetrable(NoSlip()), Insulating()),
            OceanBC(Impenetrable(NoSlip()), Insulating()),
            OceanBC(Penetrable(KinematicStress()), TemperatureFlux()),
        ),
    ) where {FT <: AbstractFloat}
        return new{FT, typeof(BC)}(Lˣ, Lʸ, H, τₒ, λʳ, θᴱ, BC)
    end
end

"""
    ocean_init_state!(::OceanGyre)

initialize u,v,η with 0 and θ linearly distributed between 9 at z=0 and 1 at z=H

# Arguments
- `p`: OceanGyre problem object, used to dispatch on and obtain ocean height H
- `Q`: state vector
- `A`: auxiliary state vector, not used
- `localgeo`: the local geometry information
- `t`: time to evaluate at, not used
"""
function ocean_init_state!(
    ::Union{HBModel, OceanModel},
    p::OceanGyre,
    Q,
    A,
    localgeo,
    t,
)
    coords = localgeo.coord
    @inbounds y = coords[2]
    @inbounds z = coords[3]
    @inbounds H = p.H

    Q.u = @SVector [-0, -0]
    Q.η = -0
    Q.θ = (5 + 4 * cos(y * π / p.Lʸ)) * (1 + z / H)

    return nothing
end

function ocean_init_state!(
    ::Union{SWModel, BarotropicModel},
    ::OceanGyre,
    Q,
    A,
    localgeo,
    t,
)
    Q.U = @SVector [-0, -0]
    Q.η = -0

    return nothing
end

"""
    kinematic_stress(::OceanGyre)

jet stream like windstress

# Arguments
- `p`: problem object to dispatch on and get additional parameters
- `y`: y-coordinate in the box
"""
@inline kinematic_stress(p::OceanGyre, y, ρ) =
    @SVector [(p.τₒ / ρ) * cos(y * π / p.Lʸ), -0]

@inline kinematic_stress(p::OceanGyre, y) =
    @SVector [-p.τₒ * cos(π * y / p.Lʸ), -0]

"""
    surface_flux(::OceanGyre)

cool-warm north-south linear temperature gradient

# Arguments
- `p`: problem object to dispatch on and get additional parameters
- `y`: y-coordinate in the box
- `θ`: temperature within element on boundary
"""
@inline function surface_flux(p::OceanGyre, y, θ)
    Lʸ = p.Lʸ
    θᴱ = p.θᴱ
    λʳ = p.λʳ

    θʳ = θᴱ * (1 - y / Lʸ)
    return λʳ * (θ - θʳ)
end


text:
"""
Time-lock Puzzle
================

This module reimplements Ron Rivest's "LCS35 Time Capsule" given at the 35th birthday of MIT's
Laboratory for Computer Science.

# What are missing in the code:

1. Estimate wall time needed to crack the puzzle according to Moore's law ("The free lunch is over",
   anyway).
2. Error checking during computation.
3. Attach information that allows one to factor `n`.

# References

- LCS35 Time Capsule Crypto-Puzzle (http://people.csail.mit.edu/rivest/lcs35-puzzle-description.txt)
- Time-lock encryption (http://www.gwern.net/Self-decrypting-files)
"""
module TimeLockPuzzle

using Primes: nextprime

PRIME_LENGTH = 1024
TWO_POWER = big(1) << PRIME_LENGTH

# The `Puzzle` struct can be (de)serialized but the code is omitted here.
"The output of the generated puzzle"
struct Puzzle
    "The modulo"
    n::BigInt

    "The time (steps of computation) expected to crack the puzzle"
    t::BigInt

    "Encrypted secret message, encoded as `BigInt`"
    enc_message::BigInt
end

# The string encoding is UTF-8. For secret messages written in Eastern-asian languages, UTF-16 is
# preferred, but the Julia community discriminates against the poor UTF-16 and puts it into
# "LegacyStrings.jl".

# Steps to create the puzzle:
#
# 1. Encrypt message `M` with key `K` using a symmetric encryption algorithm
#
#         C_M = Encrypt(K, M)
#
# 2. Encrypt the key
#
#         C_K = K + a ^ (2^t)  (mod n)
#
#     By Euler's theorem, to compute this efficiently, one can first compute
#
#         e = 2 ^ t  (mod phi(n))
#
#     and then
#
#         b = a ^ e  (mod n)
#
# 3. The output of the puzzle is
#
#         (n, a, t, C_M, C_K)
#
# --------------------------------------------------------------------------------------------------
#
# The implementation is slightly different from the description in the paper. Here `a` is always 2,
# `C_K` is always 0, and the encrytion algorithm is just XOR. Therefore the actual output tuple is
# `(n, t, C_M)`.
function create_puzzle(p_seed::BigInt,
                       q_seed::BigInt,
                       t::BigInt,
                       secret::String)::Union{Puzzle, Nothing}
    @assert p_seed != q_seed "Two identical seeds"

    # 5 has maximal order modulo 2^k (Let's assume it's right).
    p = nextprime(powermod(big(5), p_seed, TWO_POWER))
    q = nextprime(powermod(big(5), q_seed, TWO_POWER))

    n = p * q
    # Euler's totient function
    phi = (p - 1) * (q - 1)

    # Generates the puzzle
    u = powermod(big(2), t, phi)
    w = powermod(big(2), u, n)

    secret_bytes = Vector{UInt8}(secret)
    secret_big = convert(BigInt, secret_bytes)
    if secret_big > n
        return nothing
    end

    enc_secret_big = xor(secret_big, w)
    Puzzle(n, t, enc_secret_big)
end

# The puzzle can be solved by performing t successive squarings modulo n, beginning with 2.
#
#     W(0) = 2
#     W(i+1) = W(i)^2  (mod n)    for i>0
#
# There is no known way to perform this more quickly without knowing the factorization of n.
function solve_puzzle(puzzle::Puzzle)::String
    t = puzzle.t
    n = puzzle.n

    k = big(2)
    for _i in 1:t
        k = (k * k) % n
    end

    dec_message = xor(puzzle.enc_message, k)
    String(convert(Vector{UInt8}, dec_message))
end

import Base: convert

"`BigInt` -> `Vector{UInt8}` conversion"
function convert(::Type{Vector{UInt8}}, x::BigInt)
    sz = div(ndigits(x, base=2) + 8 - 1, 8)
    bytes = UInt8[]
    sizehint!(bytes, sz)
    while x != 0
        push!(bytes, x & 255)
        x >>= 8
    end
    reverse!(bytes)
end

"`Vector{UInt8}` -> `BigInt` conversion"
function convert(::Type{BigInt}, v::Vector{UInt8})
    foldl((acc, b) -> (acc << 8) + b, v; init=big(0))
end

function run()
    t = big(1048576)  # A distant future
    maybe_puzzle = create_puzzle(big(17070415), big(17830918), t,
        "I have proven Euler's sum of powers conjecture.")  # Some nonsense
    if maybe_puzzle == nothing
        error("Cannot create the puzzle")
    end

    dec_msg = solve_puzzle(maybe_puzzle)
    println("The secret message is: \"$dec_msg\"")
end

end  # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
################
# Content Type #
################

@ghdef mutable struct Content
    typ::Union{String, Nothing}
    filename::Union{String, Nothing}
    name::Union{String, Nothing}
    path::Union{String, Nothing}
    target::Union{String, Nothing}
    encoding::Union{String, Nothing}
    content::Union{String, Nothing}
    sha::Union{String, Nothing}
    url::Union{HTTP.URI, Nothing}
    git_url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    download_url::Union{HTTP.URI, Nothing}
    size::Union{Int, Nothing}
end

Content(path::AbstractString) = Content(Dict("path" => path))

namefield(content::Content) = content.path

###############
# API Methods #
###############

@api_default function file(api::GitHubAPI, repo, path; options...)
    result = gh_get_json(api, content_uri(repo, path); options...)
    return Content(result)
end

@api_default function directory(api::GitHubAPI, repo, path; options...)
    results, page_data = gh_get_paged_json(api, content_uri(repo, path); options...)
    return map(Content, results), page_data
end

@api_default function create_file(api::GitHubAPI, repo, path; options...)
    result = gh_put_json(api, content_uri(repo, path); options...)
    return build_content_response(result)
end

@api_default function update_file(api::GitHubAPI, repo, path; options...)
    result = gh_put_json(api, content_uri(repo, path); options...)
    return build_content_response(result)
end

@api_default function delete_file(api::GitHubAPI, repo, path; options...)
    result = gh_delete_json(api, content_uri(repo, path); options...)
    return build_content_response(result)
end

@api_default function readme(api::GitHubAPI, repo; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/readme"; options...)
    return Content(result)
end

function permalink(content::Content, commit)
    url = string(content.html_url)
    prefix = something(content.typ, "") == "file" ? "blob" : "tree"
    rgx = Regex(string("/", prefix, "/.*?/"))
    replacement = string("/", prefix, "/", name(commit), "/")
    return HTTP.URI(replace(url, rgx => replacement))
end

###########################
# Content Utility Methods #
###########################

content_uri(repo, path) = "/repos/$(name(repo))/contents/$(name(path))"

function build_content_response(json::Dict)
    results = Dict()
    haskey(json, "commit") && setindex!(results, Commit(json["commit"]), "commit")
    haskey(json, "content") && setindex!(results, isnothing(json["content"]) ? nothing : Content(json["content"]), "content")
    return results
end


text:
load("core.jl")

srand(1)
means = [0.1, 0.1, 0.1, 0.1, 0.9]
n = length(means)
means = means[randperm(n)]
current_arms = map(mu -> BernoulliArm(mu), means)
println("Best arm is $(ind_max(means))")

total_results = zeros(0, 6)

for eta in [0.1:0.1:0.5]
  algo = Hedge(n, eta)
  results = test_algorithm(algo, current_arms, 5000, 250)
  results = hcat(repmat([eta], size(results, 1), 1), results)
  total_results = vcat(total_results, results)
end

csvwrite("algorithms/hedge/hedge_results.csv", total_results)

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
const ARGS_TOKENS = Dict{Char, Vector{Pair{TokenFinder, Symbol}}}(
    '\\' => [
        forward_match("\\\"") => :SKIP,
    ],
    '"' => [
        forward_match("\"", ['"'], false)     => :SINGLE_QUOTE,
        forward_match("\"\"\"", ['"'], false) => :TRIPLE_QUOTE
    ]
)


text:
mutable struct DEStats
  nf::Int
  nf2::Int
  nw::Int
  nsolve::Int
  njacs::Int
  nnonliniter::Int
  nnonlinconvfail::Int
  ncondition::Int
  naccept::Int
  nreject::Int
  maxeig::Float64
end

DEStats(x=-1) = DEStats(ntuple(_->x, 10)..., 0.0)

function Base.show(io::IO, s::DEStats)
  println(io, summary(s))
  println(io, "Number of function 1 evaluations: $(s.nf)")
  println(io, "Number of function 2 evaluations: $(s.nf2)")
  println(io, "Number of W matrix evaluations: $(s.nw)")
  println(io, "Number of linear solves:        $(s.nsolve)")
  println(io, "Number of Jacobians created:        $(s.njacs)")
  println(io, "Number of nonlinear solver iterations:        $(s.nnonliniter)")
  println(io, "Number of nonlinear solver convergence failures:        $(s.nnonlinconvfail)")
  println(io, "Number of rootfind condition calls:        $(s.ncondition)")
  println(io, "Number of accepted steps:       $(s.naccept)")
  print(io,   "Number of rejected steps:       $(s.nreject)")
  iszero(s.maxeig) || print(io, "\nMaximum eigenvalue recorded:    $(s.maxeig)")
end


text:
module CodecXz

export
    XzCompressor,
    XzCompressorStream,
    XzDecompressor,
    XzDecompressorStream

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error,
    initialize,
    finalize,
    splitkwargs
using Libdl
using XZ_jll

include("liblzma.jl")
include("compression.jl")
include("decompression.jl")

end # module


text:
animals = [
    ("fly", "I don't know why she swallowed a fly, perhaps she'll die."),
    ("spider", "It wiggled and jiggled and tickled inside her."),
    ("bird", "How absurd, to swallow a bird."),
    ("cat", "Imagine that, she swallowed a cat."),
    ("dog", "What a hog, to swallow a dog."),
    ("goat", "She just opened her throat and swallowed a goat."),
    ("cow", "I don't know how she swallowed a cow."),
    ("horse", "She's dead, of course.")]

for (i, (animal, lyric)) in enumerate(animals)
    println("There was an old lady who swallowed a $animal.\n$lyric")

    if animal == "horse" break end

    for ((predator, _), (prey, _)) in zip(animals[i:-1:1], animals[i-1:-1:1])
        println("\tShe swallowed the $predator to catch the $prey")
    end

    if animal != "fly" println(animals[1][2]) end  # fly lyric
    println()  # new line
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestGBBaseLearner

using FactCheck
importall GradientBoost.GBBaseLearner
importall GradientBoost.LossFunctions

type DummyLearner; end
type StubLearner; end

function GBBaseLearner.learner_fit(lf::LossFunction, learner::StubLearner,
  instances, labels)

  return (instance) -> mean(instance)
end

function GBBaseLearner.learner_predict(lf::LossFunction, learner::StubLearner,
  model, instances)

  pred_func = model
  num_instances = size(instances, 1)

  predictions = [pred_func(instances[i,:]) for i = 1:num_instances]
end

facts("GB Learner") do
  context("not implemented functions throw an error") do
    dl = DummyLearner()
    emp_mat = Array(Any, 1, 1)
    emp_vec = Array[]
    dummy_model = emp_vec

    @fact_throws learner_fit(
      LeastSquares(),
      dgb,
      emp_mat,
      emp_vec
    )
    @fact_throws learner_predict(
      LeastSquares(),
      dgb,
      dummy_model,
      emp_mat
    )
  end

  context("build_base_func works") do
    sl = StubLearner()
    gb = GBBL(sl)
    instances = [
      1 1;
      2 4;
    ]
    labels = [
      2;
      6;
    ]

    prev_func_pred =
      fill(minimizing_scalar(gb.loss_function, labels), size(instances, 1))
    psuedo = negative_gradient(
      gb.loss_function,
      labels,
      gb.learning_rate .* prev_func_pred
    )

    base_func = build_base_func(
      gb,
      instances,
      labels,
      prev_func_pred,
      psuedo
    )

    predictions = base_func(instances)
    expected = { 1.0, 3.0 }
    @fact predictions => roughly(expected)
  end

  context("LeastSquares fit_best_constant works") do
    lf = LeastSquares()
    dummy_vec = [0.0,0.0,0.0,0.0]
    expected = 1.0

    actual = GBBaseLearner.fit_best_constant(
      lf, dummy_vec, dummy_vec, dummy_vec, dummy_vec
    )
    @fact actual => expected
  end
  context("LeastAbsoluteDeviation fit_best_constant works") do
    lf = LeastAbsoluteDeviation()
    dummy_vec = [0.0,0.0,0.0,0.0]
    labels = [0.0,1.0,2.0,3.0]
    psuedo_pred = [3.0,2.0,1.0,0.0]
    prev_func_pred = [1.0,0.0,1.0,0.0]
    expected = -0.333333

    actual = GBBaseLearner.fit_best_constant(
      lf, labels, dummy_vec, psuedo_pred, prev_func_pred
    )
    @fact actual => roughly(expected)
  end
  context("BinomialDeviance fit_best_constant throws error") do
    lf = BinomialDeviance()
    dummy_vec = [0.0,0.0,0.0,0.0]

    @fact_throws GBBaseLearner.fit_best_constant(
      lf, dummy_vec, dummy_vec, dummy_vec, dummy_vec
    )
  end
end

end # module


text:
@testset "Internal Methods" begin
    m = InfiniteModel()
    @infinite_parameter(m, x in Normal())
    @infinite_parameter(m, xi[1:2] in MvNormal(ones(2), [1. 0.; 0. 2.]))
    dx = dispatch_variable_ref(x)
    dxi = dispatch_variable_ref.(xi)
    set1 = UniDistributionSet(Normal())
    set2 = MultiDistributionSet(MvNormal(ones(2), [1. 0.; 0. 2.]))
    set3 = CollectionSet([set1, set1])
    # test _has_distribution_set
    @testset "_has_distribution_set" begin
        @test IOMT._has_distribution_set(dx, set1)
        @test IOMT._has_distribution_set(dxi[1], set2)
        @test IOMT._has_distribution_set(dxi[1], set3)
        @test !IOMT._has_distribution_set(dx, BadSet())
    end
    # test _expect_coeffs
    @testset "_expect_coeffs" begin
        @test IOMT._expect_coeffs([1, 2]) == [0.5, 0.5]
    end
end

# TODO test expect with JuMP container of infinite parameters

@testset "Expect" begin
    m = InfiniteModel()
    @infinite_parameter(m, x in Normal())
    @infinite_parameter(m, xi[1:2] in MvNormal(ones(2), [1. 0.; 0. 2.]))
    @infinite_parameter(m, y in [0, 1])
    @infinite_variable(m, inf(x, xi, y))

    @test InfiniteOpt._index_type(expect(inf, x)) == MeasureIndex
    warn = "Cannot specify a nonzero `min_num_supports` for individual " *
           "dependent parameters."
    @test_logs (:warn, warn) expect(inf, xi[1], min_num_supports = 3)
    @test !InfiniteOpt._is_expect(InfiniteOpt._core_variable_object(expect(inf, y)).data)
    @test 𝔼(inf, x) isa GeneralVariableRef
end

@testset "Macro" begin
    m = InfiniteModel()
    @infinite_parameter(m, x in Normal())
    @infinite_parameter(m, xi[1:2] in MvNormal(ones(2), [1. 0.; 0. 2.]))
    @infinite_variable(m, inf(x, xi))

    @test InfiniteOpt._index_type(@expect(inf, x)) == MeasureIndex
    @test InfiniteOpt._index_type(@expect(inf, xi, min_num_supports = 20)) == MeasureIndex
    @test_macro_throws ErrorException @expect(inf, x, 1)
    @test_macro_throws ErrorException @expect(inf, x, bob = 35)
    @test @𝔼(inf, x) isa GeneralVariableRef
end


text:
using BenchmarkTools

const SUITE = BenchmarkGroup()

include("bm_basis.jl")
include("bm_linear.jl")
SUITE["basis"] = basis_suite
SUITE["linear"] = linear_suite


text:
module PropertyDicts

"""
    PropertyDict(::AbstactDict)

Wraps an `AbstractDict` to add `getproperty` support for `Symbol` and
`AbstractString` keys.

e.g.
```
> julia d = PropertyDict(Dict("foo" => 1, :bar => 2))
PropertyDict with 2 entries:
  :bar  => 2
  "foo" => 1

julia> d.foo, d.bar, d."foo"
(1, 2, 1)

julia> d."bar"
ERROR: KeyError: key "bar" not found
```
"""
struct PropertyDict{K, V, T <: AbstractDict{K, V}} <: AbstractDict{K, V}
    d::T
    PropertyDict(d::T) where {T <: AbstractDict} =
        new{keytype(d), valtype(d), T}(d)
end

unwrap(d::PropertyDict) = getfield(d, :d)

Base.getproperty(d::PropertyDict{AbstractString}, n::Symbol) =
    getindex(d, String(n))

function Base.getproperty(d::PropertyDict, n::Symbol)
    v = get(d, n, Base.secret_table_token)
    if v != Base.secret_table_token
        return v
    end
    return getindex(d, String(n))
end

Base.getproperty(d::PropertyDict, n) = getindex(d, n)

Base.IteratorSize(::Type{PropertyDict{K,V,T}}) where {K,V,T} = Base.IteratorSize(T)
Base.IteratorEltype(::Type{PropertyDict{K,V,T}}) where {K,V,T} = Base.IteratorEltype(T)
Base.getindex(d::PropertyDict, i) = getindex(unwrap(d), i)
Base.get(d::PropertyDict, k, default) = get(unwrap(d), k, default)
Base.length(d::PropertyDict) = length(unwrap(d))
Base.iterate(d::PropertyDict) = iterate(unwrap(d))
Base.iterate(d::PropertyDict, i) = iterate(unwrap(d), i)

Base.string(d::PropertyDict) = string(unwrap(d))

Base.convert(::Type{Any}, d::PropertyDict) = d
Base.convert(::Type{PropertyDict{K,V,T}}, d::PropertyDict{K,V,T}) where {K,V,T<:AbstractDict{K,V}} = d
Base.convert(::Type{T}, d::PropertyDict) where T <: AbstractDict =
    T === AbstractDict ? d : convert(T, PropertyDicts.unwrap(d))
Base.convert(::Type{T}, d::PropertyDict) where T = convert(T, PropertyDicts.unwrap(d))

end # module PropertyDicts


text:
# Example from JuliaOpt/SumOfSquares.jl
# COSMO supports MOI/JuMP and therefore also works with JuMP extensions
# like SumOfSquares.jl

using SumOfSquares, COSMO
using DynamicPolynomials
@polyvar x y
motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2

model = SOSModel(with_optimizer(COSMO.Optimizer, rho = 1e-5);)
@constraint(model, (x^2 + y^2) * motzkin >= 0)
optimize!(model)

primal_status(model)

# The full example can be found at:
# https://github.com/JuliaOpt/SumOfSquares.jl/blob/master/examples/motzkin.ipynb


text:
using BrainFlow

# specify logging library to use
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params)

BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)
sleep(5)
BrainFlow.stop_stream(board_shim)
data = BrainFlow.get_board_data(board_shim)
println(data)
data = BrainFlow.get_board_data(10, board_shim)
println(data)
data = BrainFlow.get_current_board_data(10, board_shim)
println(data)
BrainFlow.release_session(board_shim)


text:
##
## function get_polyA_prefixes(fasta_record::FASTA.Record,
##    minimum_not_polyA::Int64,
##    minimum_polyA_length::Int64)::Array{String,1}
##
minimum_not_polyA = 20
minimum_polyA_length = 10
re = Regex("([ATGC]{$minimum_not_polyA})A{$minimum_polyA_length,}")

# Test set 1 just a simple nt real example
seq = dna"ATCCTTGGCTTTCTTCCGGGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGAAAAAAAAAAAAAAAAAAAAAAAAAAA"
output = ["ATCCTTGGCTTTCTTCCGGG","CAAAAAAAAAAAAAAAAAAA", "GGGGGGGGGGGGGGGGGGGG"]
ref_seq=FASTA.Record("seq",seq)
@test get_polyA_prefixes(ref_seq,20,10,re)  == output

# Test set 2 no AAAA
seq = dna"ATCCTTGGCTTTCTTCCGGGCGGGGGGGGGGGGGGGGGGGGAAA"
output = []
ref_seq=FASTA.Record("seq",seq)
@test get_polyA_prefixes(ref_seq,20,10,re)  == output

# Test set 3
seq = dna"ATCCTTGGCTTTCTTCCGGGCGGGGGGGGGGGGGGGGAGAGAGAAAAAAAAAGAAAAAAAAAA"
output = ["GGGGAGAGAGAAAAAAAAAG"]
ref_seq=FASTA.Record("seq",seq)
@test get_polyA_prefixes(ref_seq,20,10,re)  == output

# Test set 4 Only AA...
seq = dna"AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA"
output = ["AAAAAAAAAAAAAAAAAAAA"]
ref_seq=FASTA.Record("seq",seq)
@test get_polyA_prefixes(ref_seq,20,10,re)  == output

# Test set 5 Only AA... but too short
seq = dna"AAAAAAAAAAA"
output = []
ref_seq=FASTA.Record("seq",seq)
@test get_polyA_prefixes(ref_seq,20,10,re)  == output

# Test set 6
seq = dna"AAAAAAAAAAGGGGGGGGGGAAAAAAAAAA"
output = ["AAAAAAAAAAGGGGGGGGGG"]
ref_seq=FASTA.Record("seq",seq)
@test get_polyA_prefixes(ref_seq,20,10,re)  == output

# Test set 7 too short
seq = dna"GTCGTCGTC"
output = []
ref_seq=FASTA.Record("seq",seq)
@test get_polyA_prefixes(ref_seq,20,10,re)  == output

max_stars_repo_name: laurabiggins/polyAterminus

max_issues_repo_name: laurabiggins/polyAterminus

max_forks_repo_name: laurabiggins/polyAterminus


text:
@testset "constructors" begin
    node = Node()
    @test isa(node, Node) && node.value == "" && node.class == "" &&
        node.entr == 0 && node.ig == 1 && node.children == []

    node2 = Node("", "", 0, 1, [])
    @test isa(node2, Node) && node2.value == "" && node.class == "" &&
        node2.entr == 0 && node2.ig == 1 && node2.children == []

    node3 = Node("value")
    @test isa(node3, Node) && node3.value == "value" && node.class == "" &&
        node3.entr == 0 && node3.ig == 1 && node3.children == []

    node4 = Node("", [node3])
    @test isa(node4, Node) && node4.value == "" && node.class == "" &&
        node4.entr == 0 && node4.ig == 1 && node4.children == [node3]
end

@testset "accessors" begin
    node = Node("", "", 0, 1, [])
    @test value(node) == node.value
    @test class(node) == node.class
    @test entropy(node) == node.entr
    @test information_gain(node) == node.ig
    @test children(node) == node.children
    @test nchildren(node) == 0
end

@testset "properties" begin
    node = Node("", "", 0, 1, [])
    @test nchildren(node) == 0
end

@testset "operators" begin
    node = Node("", "", 0, 1, [])
    node2 = Node("", "", 0, 1, [])
    @test node == node2
end

@testset "mutations" begin
    root = Node("", "", 0, 1, [])
    node = Node("", "", 0, 1, [])
    prev_nchildren = length(root.children)
    push!(root, node)
    @test length(root.children) == prev_nchildren + 1

    root = Node("", "", 0, 1, [])
    node = Node("", "", 0, 1, [])
    prev_nchildren = length(root.children)
    append!(root, [node, node, node])
    @test length(root.children) == prev_nchildren + 3
end

@testset "show" begin
    node = Node("", "", 0, 1, [])
    @test sprint(show, node) == "()"

    node2 = Node("2", "", 0, 1, [])
    @test sprint(show, node2) == "[2] (entr=0.0000, ig=1.0000)"

    node3 = Node("3", "", 0, 1, [])
    node4 = Node("4", "", 0, 1, [])
    push!(node, node2)
    push!(node2, node4)
    push!(node2, node4)
    push!(node, node3)
    @test sprint(show, node) == """()
    |-- [2] (entr=0.0000, ig=1.0000)
    |   |-- [4] (entr=0.0000, ig=1.0000)
    |   `-- [4] (entr=0.0000, ig=1.0000)
    `-- [3] (entr=0.0000, ig=1.0000)"""
end


text:
using StatsModels
using StatsBase
using StatsModels: width
using DataStructures

@testset "Temporal Terms" begin
    @testset "Lag" begin
        @testset "Basic use" begin
            df = (y=1:10, x = 1:10)
            f = @formula(y ~ lag(x, 0) + lag(x, 1) + lag(x, 3) + lag(x, 11))
            f = apply_schema(f, schema(f, df))
            resp, pred = modelcols(f, df)

            @test isequal(pred[:, 1], 1.0:10)
            @test isequal(pred[:, 2], [missing; 1.0:9])
            @test isequal(pred[:, 3], [missing; missing; missing; 1.0:7])
            @test isequal(pred[:, 4], fill(missing, 10))

            @test coefnames(f)[2] == ["x_lag0", "x_lag1", "x_lag3", "x_lag11"]
        end

        @testset "1 arg form" begin
            df = (y=1:10, x = 1:10)
            f = @formula(y ~ lag(x))
            f = apply_schema(f, schema(f, df))
            resp, pred = modelcols(f, df)

            @test isequal(pred[:, 1], [missing; 1.0:9])
            @test coefnames(f)[2] == "x_lag1"
        end

        @testset "Row Table" begin
            rowdata = [(y=i, x=2i) for i in 1:10]
            f = @formula(y ~ lag(x))
            f = apply_schema(f, schema(f, rowdata))
            resp, pred = modelcols(f, rowdata)
            @test isequal(pred[:, 1], [missing; 2.0; 4.0; 6.0; 8.0; 10.0; 12.0; 14.0; 16.0; 18.0])
        end

        @testset "Nested Use" begin
            df = (y=1:10, x = 1:10)
            f = @formula(y ~ lag(lag(x, 1), 2))  # equiv to `lag(x, 3)`
            f = apply_schema(f, schema(f, df))
            resp, pred = modelcols(f, df);

            @test isequal(pred[:, 1], [missing; missing; missing; 1.0:7])
        end

        @testset "Negative lag" begin
            df = (y=1:10, x = 1:10)
            neg_f = @formula(y ~ lag(x, -2))
            neg_f = apply_schema(neg_f, schema(neg_f, df))
            resp, pred = modelcols(neg_f, df);

            @test isequal(pred[:, 1], [3.0:10; missing; missing])
            @test coefnames(neg_f)[2] == "x_lag-2"
        end

        @testset "Categorical Term use" begin
            df = (y=1:4, x = ["A", "B", "A", "C"])
            f = @formula(y ~ lag(x, 2))
            f = apply_schema(f, schema(f, df))
            resp, pred = modelcols(f, df)

            # Note the even though "C" is lagged out of the data, we still get 2 columns
            @test isequal(pred[:, 1], [missing; missing; 0; 1])
            @test isequal(pred[:, 2], [missing; missing; 0; 0])

            @test coefnames(f)[2] == ["x: B_lag2", "x: C_lag2"]
        end

        @testset "Diff Demo" begin
            df = (y=1:10, x = 1:10)
            f = @formula(y ~ (x - lag(x)))
            f = apply_schema(f, schema(f, df))
            # Broken because of: https://github.com/JuliaStats/StatsModels.jl/issues/114
            @test_broken resp, pred = modelcols(f, df);
            @test_broken isequal(pred[:, 1], [missing; fill(1, 9)])
        end

        @testset "Unhappy path" begin
            @testset "Variable lag" begin
                df = (y=1:5, x = 1:5, offset=[0, 1, 0, 2, 1])
                bad_f = @formula(y ~ lag(x, offset))
                @test_throws ArgumentError apply_schema(bad_f, schema(bad_f, df))
            end

            @testset "Fractional lag" begin
                df = (y=1:10, x = 1:10)
                bad_f = @formula(y ~ lag(x, 1.5))
                @test_throws InexactError apply_schema(bad_f, schema(bad_f, df))
            end
        end # Unhappy Path testset
    end # Lag testset

    # The code for lag and lead is basically the same, as we tested lag comprehensively above
    # the tests for lead are more sparse.
    @testset "Lead" begin
        @testset "Basic use" begin
            df = (y=1:10, x = 1:10)
            f = @formula(y ~ lead(x, 0) + lead(x, 1) + lead(x, 3) + lead(x, 11))
            f = apply_schema(f, schema(f, df))
            resp, pred = modelcols(f, df)

            @test isequal(pred[:, 1], 1.0:10)
            @test isequal(pred[:, 2], [2.0:10; missing])
            @test isequal(pred[:, 3], [4.0:10; missing; missing; missing])
            @test isequal(pred[:, 4], fill(missing, 10))

            @test coefnames(f)[2] == ["x_lead0", "x_lead1", "x_lead3", "x_lead11"]
        end
    end
end


text:
module LeapSeconds

using ERFA
using OptionalData
using RemoteFiles

export leapseconds, LSK, LSK_FILE, LSK_DATA

# Constants for calculating the offset between TAI and UTC for
# dates between 1960-01-01 and 1972-01-01

const EPOCHS = (
    2.4369345e6,
    2.4373005e6,
    2.4375125e6,
    2.4376655e6,
    2.4383345e6,
    2.4383955e6,
    2.4384865e6,
    2.4386395e6,
    2.4387615e6,
    2.4388205e6,
    2.4389425e6,
    2.4390045e6,
    2.4391265e6,
    2.4398875e6,
)

const OFFSETS = (
    1.417818,
    1.422818,
    1.372818,
    1.845858,
    1.945858,
    3.240130,
    3.340130,
    3.440130,
    3.540130,
    3.640130,
    3.740130,
    3.840130,
    4.313170,
    4.213170,
)

const DRIFT_EPOCHS = (
    2.4373005e6,
    2.4373005e6,
    2.4373005e6,
    2.4376655e6,
    2.4376655e6,
    2.4387615e6,
    2.4387615e6,
    2.4387615e6,
    2.4387615e6,
    2.4387615e6,
    2.4387615e6,
    2.4387615e6,
    2.4391265e6,
    2.4391265e6,
)

const DRIFT_RATES = (
    0.0012960,
    0.0012960,
    0.0012960,
    0.0011232,
    0.0011232,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0025920,
    0.0025920,
)

struct LSK
    t::Vector{Float64}
    leapseconds::Vector{Float64}
end

function LSK(file)
    t = Vector{Float64}()
    leapseconds = Vector{Float64}()
    re = r"(?<dat>[0-9]{2}),\s+@(?<date>[0-9]{4}-[A-Z]{3}-[0-9])"
    lines = open(readlines, file)
    for line in lines
        s = string(line)
        if ismatch(re, s)
            m = match(re, s)
            push!(leapseconds, float(m["dat"]))
            push!(t, Dates.datetime2julian(DateTime(m["date"], "y-u-d")))
        end
    end
    LSK(t, leapseconds)
end

@RemoteFile LSK_FILE "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls"
@OptionalData LSK_DATA LSK "Run 'AstroTime.update()' to load it."

function leapseconds(lsk::LSK, jd)
    # Before 1960-01-01
    if jd < 2.4369345e6
        return 0.0
    elseif jd < lsk.t[1]
        idx = findlast(jd .>= EPOCHS)
        return OFFSETS[idx] + (jd - DRIFT_EPOCHS[idx]) * DRIFT_RATES[idx]
    else
        return lsk.leapseconds[findlast(jd .>= lsk.t)]
    end
end
leapseconds(jd) = leapseconds(get(LSK_DATA), jd)
end


text:
@time @testset "Macrokernels" begin
    m = 20
    n = 30
    k = 40
    A1 = rand(Float64, m, k)
    B1 = rand(Float64, k, n)
    C1 = rand(Float64, m, n)
    A2 = deepcopy(A1)
    B2 = deepcopy(B1)
    C2 = deepcopy(C1)
    α = Float64(2.0)
    β = Float64(2.0)
    Octavian.loopmul!(VectorizationBase.zstridedpointer(C1), VectorizationBase.zstridedpointer(A1), VectorizationBase.zstridedpointer(B1), α, β, m, k, n)
    C2 = α*A2*B2 + β*C2
    @test C1 ≈ C2
end