text:
module All_EBC_2dir_examples
using FinEtools
using FinEtools.AlgoDeforLinearModule
using FinEtools.MeshUtilModule
using FinEtools.AlgoBaseModule
using LinearAlgebra

# Orthotropic material
E1 = 2.5e6*phun("PSI"); E2 = 1e6*phun("PSI"); E3 = E2;
G12 = 0.5e6*phun("PSI");  G13 = G12; G23 = 0.2e6*phun("PSI")
nu12 =  0.25; nu13 =  0.25; nu23 =  0.25;
# Coefficients of thermal expansion
CTE1 = CTE2 = CTE3 = 0.0

angles = vec([-15.0]);
nLayers = length(angles)
# dimensions of the plate
a = 100.0*phun("mm");
b = 100.0*phun("mm");
t = 100.0*phun("mm");

tolerance = 0.0001*t

Refinements = [2, 5, 10, 20, 40] # For the paper
Refinements = [2, 4, 8, 16] # for testing

# Here we define the layout and the thicknesses of the layers.
ts = t/nLayers * ones(nLayers); # layer thicknesses

# The material coordinate system function is defined as:
function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    rotmat3!(csmatout, angles[fe_label]/180.0*pi* [0.0; 1.0; 0.0]);
end

displacement_fun_x(x) = dot([0.0, 0.001, 0.0, 0.001, 0.0, -0.001, 0.0, 0.002, 0.0, 0.001], [1.0, x[1], x[2], x[3], x[1] * x[2], x[1] * x[3], x[2] * x[3], x[1]^2, x[2]^2, x[3]^2])
displacement_fun_y(x) = dot([0.0, 0.0, 0.0003, 0.0, 0.003, -0.002, 0.007, 0.0, 0.008, -0.007], [1.0, x[1], x[2], x[3], x[1] * x[2], x[1] * x[3], x[2] * x[3], x[1]^2, x[2]^2, x[3]^2])
displacement_fun_z(x) = dot([0.0, 0.0, 0.0, 0.00007, 0.008, 0.003, -0.015, 0.002, -0.003, 0.0001], [1.0, x[1], x[2], x[3], x[1] * x[2], x[1] * x[3], x[2] * x[3], x[1]^2, x[2]^2, x[3]^2])

function All_EBC_2dir_MST10_conv()
    elementtag = "MST10"
    println("""
    Fiber-reinforced block: $(elementtag)
    """)

    modeldatasequence = FDataDict[]
    for Refinement = Refinements

        MR = DeforModelRed3D
        material = MatDeforElastOrtho(MR, 0.0, E1, E2, E3, nu12, nu13, nu23, G12, G13, G23, CTE1, CTE2, CTE3)

        # Select how fine the mesh should be
        nts= Refinement * ones(Int, nLayers);# number of elements per layer
        tnts = sum(nts)
        na, nb = tnts, tnts;

        xs = collect(linearspace(0.0, a, na+1))
        ys = collect(linearspace(0.0, b, nb+1))
        fens,fes = T10layeredplatex(xs, ys, ts, nts)
        println("Mesh: na, nb, nts = $na, $nb, $nts")
        println("count(fens) = $(count(fens))")

        # The volume integrals are evaluated using this rule
        gr = SimplexRule(3, 4)

        # We will create one region per layer
        regions = FDataDict[]
        for layer = 1:nLayers
            rls = selectelem(fens, fes, label =  layer)
            push!(regions, FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegData(subset(fes, rls), gr), CSys(3, 3, updatecs!), material)))
        end

        # The essential boundary conditions: the entire surface
        bfes = meshboundary(fes)
        lx0 = connectednodes(bfes)
        eclamped1 = FDataDict("displacement"=> displacement_fun_x, "component"=> 1, "node_list"=>lx0)
        eclamped2 = FDataDict("displacement"=> displacement_fun_y, "component"=> 2, "node_list"=>lx0)
        eclamped3 = FDataDict("displacement"=> displacement_fun_z, "component"=> 3, "node_list"=>lx0)

        modeldata = FDataDict("fens"=>fens, "regions"=>regions, "essential_bcs"=>[eclamped1, eclamped2, eclamped3])
        modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

        modeldata["elementsize"] = t/Refinement
        modeldata["geometricaltolerance"] = tolerance
        push!(modeldatasequence, modeldata)
    end # for refinement

    for (extrap, nodevalmeth) = zip([:extrapmean, :extraptrend, :default], [:averaging, :averaging, :invdistance])
        filebase = "All_EBC_2dir_$(elementtag)_$(extrap)"
        for modeldata = modeldatasequence
            u = modeldata["u"]
            geom = modeldata["geom"]
            modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s", "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3), "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
            modeldata = AlgoDeforLinearModule.exportstress(modeldata)
        end

        println("")
        println("Stress RMS error")
        for md = modeldatasequence
            md["targetfields"] = [e["field"] for e in md["postprocessing"]["exported"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Stress" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        println("")
        println("Displacement RMS error")
        for md = modeldatasequence
            md["targetfields"] = [md["u"] for r in md["regions"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Displ" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")
    end
    println("Done")

end # All_EBC_2dir_MST10_conv

function All_EBC_2dir_MSH8_conv()
    elementtag = "MSH8"
    println("""
    Fiber-reinforced block: $(elementtag)
    """)

    modeldatasequence = FDataDict[]
    for Refinement = Refinements

        # This is the material  model
        MR = DeforModelRed3D
        material = MatDeforElastOrtho(MR, 0.0, E1, E2, E3, nu12, nu13, nu23, G12, G13, G23, CTE1, CTE2, CTE3)

        # Select how fine the mesh should be
        nts= Refinement * ones(Int, nLayers);# number of elements per layer
        tnts = sum(nts)
        na, nb = tnts, tnts;

        xs = collect(linearspace(0.0, a, na+1))
        ys = collect(linearspace(0.0, b, nb+1))
        fens,fes = H8layeredplatex(xs, ys, ts, nts)
        println("Mesh: na, nb, nts = $na, $nb, $nts")
        println("count(fens) = $(count(fens))")

        # The volume integrals are evaluated using this rule
        gr = GaussRule(3, 2)

        # We will create one region per layer
        regions = FDataDict[]
        for layer = 1:nLayers
            rls = selectelem(fens, fes, label =  layer)
            push!(regions, FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegData(subset(fes, rls), gr), CSys(3, 3, updatecs!), material)))
        end

        # The essential boundary conditions: the entire surface
        bfes = meshboundary(fes)
        lx0 = connectednodes(bfes)
        eclamped1 = FDataDict("displacement"=> displacement_fun_x, "component"=> 1, "node_list"=>lx0)
        eclamped2 = FDataDict("displacement"=> displacement_fun_y, "component"=> 2, "node_list"=>lx0)
        eclamped3 = FDataDict("displacement"=> displacement_fun_z, "component"=> 3, "node_list"=>lx0)

        modeldata = FDataDict("fens"=>fens, "regions"=>regions, "essential_bcs"=>[eclamped1, eclamped2, eclamped3])
        modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

        modeldata["elementsize"] = t/Refinement
        modeldata["geometricaltolerance"] = tolerance
        push!(modeldatasequence, modeldata)
    end # for refinement

    for (extrap, nodevalmeth) = zip([:extrapmean, :extraptrend, :default], [:averaging, :averaging, :invdistance])
        filebase = "All_EBC_2dir_$(elementtag)_$(extrap)"
        for modeldata = modeldatasequence
            u = modeldata["u"]
            geom = modeldata["geom"]
            modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s", "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3), "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
            modeldata = AlgoDeforLinearModule.exportstress(modeldata)
        end

        println("")
        println("Stress RMS error")
        for md = modeldatasequence
            md["targetfields"] = [e["field"] for e in md["postprocessing"]["exported"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Stress" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        println("")
        println("Displacement RMS error")
        for md = modeldatasequence
            md["targetfields"] = [md["u"] for r in md["regions"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Displ" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        # @async run(`"paraview.exe" $csvFile`)
    end # extrap
    println("Done")

end # All_EBC_2dir_MSH8_conv

function All_EBC_2dir_MSH8_conv_alt()
    elementtag = "MSH8"
    println("""
    Fiber-reinforced block: $(elementtag) Trapezoidal rule
    """)

    modeldatasequence = FDataDict[]
    for Refinement = Refinements

        # This is the material  model
        MR = DeforModelRed3D
        material = MatDeforElastOrtho(MR, 0.0, E1, E2, E3, nu12, nu13, nu23, G12, G13, G23, CTE1, CTE2, CTE3)

        # Select how fine the mesh should be
        nts= Refinement * ones(Int, nLayers);# number of elements per layer
        tnts = sum(nts)
        na, nb = tnts, tnts;

        xs = collect(linearspace(0.0, a, na+1))
        ys = collect(linearspace(0.0, b, nb+1))
        fens,fes = H8layeredplatex(xs, ys, ts, nts)
        println("Mesh: na, nb, nts = $na, $nb, $nts")
        println("count(fens) = $(count(fens))")

        # The volume integrals are evaluated using this rule
        gr = TrapezoidalRule(3)

        # We will create one region per layer
        regions = FDataDict[]
        for layer = 1:nLayers
            rls = selectelem(fens, fes, label =  layer)
            push!(regions, FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegData(subset(fes, rls), gr), CSys(3, 3, updatecs!), material)))
        end

        # The essential boundary conditions: the entire surface
        bfes = meshboundary(fes)
        lx0 = connectednodes(bfes)
        eclamped1 = FDataDict("displacement"=> displacement_fun_x, "component"=> 1, "node_list"=>lx0)
        eclamped2 = FDataDict("displacement"=> displacement_fun_y, "component"=> 2, "node_list"=>lx0)
        eclamped3 = FDataDict("displacement"=> displacement_fun_z, "component"=> 3, "node_list"=>lx0)

        modeldata = FDataDict("fens"=>fens, "regions"=>regions, "essential_bcs"=>[eclamped1, eclamped2, eclamped3])
        modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

        modeldata["elementsize"] = t/Refinement
        modeldata["geometricaltolerance"] = tolerance
        push!(modeldatasequence, modeldata)
    end # for refinement

    for (extrap, nodevalmeth) = zip([:extrapmean, :extraptrend, :default], [:averaging, :averaging, :invdistance])
        filebase = "All_EBC_2dir_$(elementtag)_$(extrap)_trapezoidal"
        for modeldata = modeldatasequence
            u = modeldata["u"]
            geom = modeldata["geom"]
            modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s", "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3), "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
            modeldata = AlgoDeforLinearModule.exportstress(modeldata)
        end

        println("")
        println("Stress RMS error")
        for md = modeldatasequence
            md["targetfields"] = [e["field"] for e in md["postprocessing"]["exported"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Stress" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        println("")
        println("Displacement RMS error")
        for md = modeldatasequence
            md["targetfields"] = [md["u"] for r in md["regions"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Displ" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        # @async run(`"paraview.exe" $csvFile`)
    end # extrap
    println("Done")

end # All_EBC_2dir_MSH8_conv_alt

function All_EBC_2dir_T10_conv()
    elementtag = "T10"
    println("""
    Fiber-reinforced block: $(elementtag)
    """)

    modeldatasequence = FDataDict[]
    for Refinement = Refinements

        # This is the material  model
        MR = DeforModelRed3D
        material = MatDeforElastOrtho(MR, 0.0, E1, E2, E3, nu12, nu13, nu23, G12, G13, G23, CTE1, CTE2, CTE3)

        # Select how fine the mesh should be
        nts= Refinement * ones(Int, nLayers);# number of elements per layer
        tnts = sum(nts)
        na, nb = tnts, tnts;

        xs = collect(linearspace(0.0, a, na+1))
        ys = collect(linearspace(0.0, b, nb+1))
        fens,fes = T10layeredplatex(xs, ys, ts, nts)
        println("Mesh: na, nb, nts = $na, $nb, $nts")
        println("count(fens) = $(count(fens))")

        # The volume integrals are evaluated using this rule
        gr = SimplexRule(3, 4)

        # We will create one region per layer
        regions = FDataDict[]
        for layer = 1:nLayers
            rls = selectelem(fens, fes, label =  layer)
            push!(regions, FDataDict("femm"=>FEMMDeforLinear(MR, IntegData(subset(fes, rls), gr), CSys(3, 3, updatecs!), material)))
        end

        # The essential boundary conditions: the entire surface
        bfes = meshboundary(fes)
        lx0 = connectednodes(bfes)
        eclamped1 = FDataDict("displacement"=> displacement_fun_x, "component"=> 1, "node_list"=>lx0)
        eclamped2 = FDataDict("displacement"=> displacement_fun_y, "component"=> 2, "node_list"=>lx0)
        eclamped3 = FDataDict("displacement"=> displacement_fun_z, "component"=> 3, "node_list"=>lx0)

        modeldata = FDataDict("fens"=>fens, "regions"=>regions, "essential_bcs"=>[eclamped1, eclamped2, eclamped3])
        modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

        modeldata["elementsize"] = t/Refinement
        modeldata["geometricaltolerance"] = tolerance
        push!(modeldatasequence, modeldata)
    end # for refinement

    for (extrap, nodevalmeth) = zip([:default], [:invdistance])
        filebase = "All_EBC_2dir_$(elementtag)_$(extrap)"
        for modeldata = modeldatasequence
            u = modeldata["u"]
            geom = modeldata["geom"]
            modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s", "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3), "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
            modeldata = AlgoDeforLinearModule.exportstress(modeldata)
        end

        println("")
        println("Stress RMS error")
        for md = modeldatasequence
            md["targetfields"] = [e["field"] for e in md["postprocessing"]["exported"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Stress" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        println("")
        println("Displacement RMS error")
        for md = modeldatasequence
            md["targetfields"] = [md["u"] for r in md["regions"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Displ" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        # @async run(`"paraview.exe" $csvFile`)
    end # extrap
    println("Done")

end # All_EBC_2dir_T10_conv

function All_EBC_2dir_H8_conv()
    elementtag = "H8"
    println("""
    Fiber-reinforced block: $(elementtag)
    """)

    modeldatasequence = FDataDict[]
    for Refinement = Refinements

        # This is the material  model
        MR = DeforModelRed3D
        material = MatDeforElastOrtho(MR, 0.0, E1, E2, E3, nu12, nu13, nu23, G12, G13, G23, CTE1, CTE2, CTE3)

        # Select how fine the mesh should be
        nts= Refinement * ones(Int, nLayers);# number of elements per layer
        tnts = sum(nts)
        na, nb = tnts, tnts;

        xs = collect(linearspace(0.0, a, na+1))
        ys = collect(linearspace(0.0, b, nb+1))
        fens,fes = H8layeredplatex(xs, ys, ts, nts)
        println("Mesh: na, nb, nts = $na, $nb, $nts")
        println("count(fens) = $(count(fens))")

        # The volume integrals are evaluated using this rule
        gr = GaussRule(3, 2)

        # We will create one region per layer
        regions = FDataDict[]
        for layer = 1:nLayers
            rls = selectelem(fens, fes, label =  layer)
            push!(regions, FDataDict("femm"=>FEMMDeforLinear(MR, IntegData(subset(fes, rls), gr), CSys(3, 3, updatecs!), material)))
        end

        # The essential boundary conditions: the entire surface
        bfes = meshboundary(fes)
        lx0 = connectednodes(bfes)
        eclamped1 = FDataDict("displacement"=> displacement_fun_x, "component"=> 1, "node_list"=>lx0)
        eclamped2 = FDataDict("displacement"=> displacement_fun_y, "component"=> 2, "node_list"=>lx0)
        eclamped3 = FDataDict("displacement"=> displacement_fun_z, "component"=> 3, "node_list"=>lx0)

        modeldata = FDataDict("fens"=>fens, "regions"=>regions, "essential_bcs"=>[eclamped1, eclamped2, eclamped3])
        modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

        modeldata["elementsize"] = t/Refinement
        modeldata["geometricaltolerance"] = tolerance
        push!(modeldatasequence, modeldata)
    end # for refinement

    for (extrap, nodevalmeth) = zip([:default], [:invdistance])
        filebase = "All_EBC_2dir_$(elementtag)_$(extrap)"
        for modeldata = modeldatasequence
            u = modeldata["u"]
            geom = modeldata["geom"]
            modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s", "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3), "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
            modeldata = AlgoDeforLinearModule.exportstress(modeldata)
        end

        println("")
        println("Stress RMS error")
        for md = modeldatasequence
            md["targetfields"] = [e["field"] for e in md["postprocessing"]["exported"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Stress" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        println("")
        println("Displacement RMS error")
        for md = modeldatasequence
            md["targetfields"] = [md["u"] for r in md["regions"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Displ" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        # @async run(`"paraview.exe" $csvFile`)
    end # extrap
    println("Done")

end # All_EBC_2dir_H8_conv

function allrun()
    println("#####################################################")
    println("# All_EBC_2dir_MSH8_conv_alt ")
    All_EBC_2dir_MSH8_conv_alt()
    println("#####################################################")
    println("# All_EBC_2dir_MST10_conv ")
    All_EBC_2dir_MST10_conv()
    println("#####################################################")
    println("# All_EBC_2dir_MSH8_conv ")
    All_EBC_2dir_MSH8_conv()
    println("#####################################################")
    println("# All_EBC_2dir_T10_conv ")
    All_EBC_2dir_T10_conv()
    println("#####################################################")
    println("# All_EBC_2dir_H8_conv ")
    All_EBC_2dir_H8_conv()
    return true
end # function allrun

end # module Pagano_examples


text:
#!/usr/bin/env julia
# Execute this with the JULIA_NUM_THREADS environment variable set to
# the number of threads on the machine.
#
# Prior to running, ensure all packages have been added. For
# LightGraphs, ensure that LightGraphs master has been added
# since we're using the new shortest paths syntax and a new
# core_number algorithm.
using BenchmarkTools
using LightGraphs, LightGraphs.ShortestPaths, LightGraphs.Centrality
using LightGraphs.Traversals: BreadthFirst, ThreadedBreadthFirst, distances
using LightGraphs.ShortestPaths: shortest_paths, BFS, ThreadedBFS
using LightGraphs.Centrality: centrality, PageRank, ThreadedPageRank
using LightGraphs.Degeneracy: core_number, KabirMadduri
using LightGraphs.Connectivity: connected_components, Kosaraju
using LinearAlgebra.BLAS

using GraphIO
using SortingAlgorithms: RadixSort
f = ARGS[1]
n = parse(Int, ARGS[2])
s = 120

println("using $(Base.Threads.nthreads()) threads")

println("loadgraph")
# we use @time here because we only want a single sample, and we want
# to keep the result.
@time g = squash(loadgraph(f, GraphIO.EdgeList.EdgeListFormat()))
println("loaded $g")
println("\n")

println("shortest paths: with path data")
show(stdout, MIME"text/plain"(), @benchmark b = shortest_paths($g, 1, BFS(RadixSort)) samples = n seconds = s)
println("\n")

println("shortest paths: threaded, with path data")
show(stdout, MIME"text/plain"(), @benchmark b = shortest_paths($g, 1, ThreadedBFS()) samples = n seconds = s)
println("\n")

println("shortest paths: distances only")
show(stdout, MIME"text/plain"(), @benchmark b = distances($g, 1, BreadthFirst(sort_alg=RadixSort)) samples = n seconds = s)
println("\n")

println("shortest paths: threaded, distances only")
show(stdout, MIME"text/plain"(), @benchmark b = distances($g, 1, ThreadedBreadthFirst()) samples = n seconds = s)
println("\n")

println("pagerank")
show(stdout, MIME"text/plain"(), @benchmark centrality($g, PageRank(0.85, 10000000, 1e-3)) samples = n seconds = s)
println("\n")

# note: this might not be as performant because BLAS is inherently multi-threaded.
println("pagerank: threaded")
nb = ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
BLAS.set_num_threads(1)
show(stdout, MIME"text/plain"(), @benchmark centrality($g, ThreadedPageRank(0.85, 10000000, 1e-3)) samples = n seconds = s)
BLAS.set_num_threads(nb)
println("\n")

println("core_number")
show(stdout, MIME"text/plain"(), @benchmark core_number($g) samples = n seconds = s)
println("\n")

println("core_number: threaded")
show(stdout, MIME"text/plain"(), @benchmark core_number($g, KabirMadduri()) samples = n seconds = s)
println("\n")

println("scc: Tarjan")
show(stdout, MIME"text/plain"(), @benchmark connected_components($g) samples = n seconds = s)
println("\n")

ext: jl

max_stars_repo_name: timlrx/graph-benchmarks

max_issues_repo_name: timlrx/graph-benchmarks

max_forks_repo_name: timlrx/graph-benchmarks


text:
using QuOptimalControl

using QuantumInformation
using DelimitedFiles

ψ1 = [1.0 + 0.0im 0.0]
ψt = [0.0 + 0.0im 1.0]

ρ1 = ψ1' * ψ1
ρt = ψt' * ψt

prob = ClosedStateTransfer(
    B = [[sx, sy]],
    A = [0.0 * sz],
    X_init = [ρ1],
    X_target = [ρt],
    duration = 1.0,
    n_timeslices = 10,
    n_controls = 2,
    n_ensemble = 1,
    norm2 = 1.0,
    alg = GRAPE_approx(inplace = true),
    initial_guess = rand(2, 10),
)

sol = solve(prob)

prob = ClosedStateTransfer(
    B = [[sx, sy]],
    A = [0.0 * sz],
    X_init = [ρ1],
    X_target = [ρt],
    duration = 1.0,
    n_timeslices = 10,
    n_controls = 2,
    n_ensemble = 1,
    norm2 = 1,
    alg = GRAPE_AD(),
    initial_guess = rand(2, 10),
)

sol = solve(prob)

prob_dCRAB = ClosedStateTransfer(
    B = [[sx, sy]],
    A = [0.0 * sz],
    X_init = ρ1,
    X_target = ρt,
    duration = 1.0,
    n_timeslices = 10,
    n_controls = 2,
    n_ensemble = 1,
    norm2 = 1.0,
    alg = dCRAB_options(),
    initial_guess = rand(2, 10),
)

sol = solve(prob_dCRAB)

visualise_pulse(sol.optimised_pulses, duration = prob.duration)

"""
Dummy function that simply saves a random number in the result.txt file
"""
function start_exp()
    open("result.txt", "w") do io
        writedlm(io, rand())
    end
end

"""
Example functional for closed loop optimal control using dCRAB
"""
function user_functional_expt(x)
    # we get a 2D array of pulses, with delimited files we can write this to a file
    # writedlm("./")
    open("pulses.txt", "w") do io
        writedlm(io, x')
    end
    # then lets imagine we use PyCall or something to start exp
    start_exp()
    # then read in the result
    infid = readdlm("result.txt")[1]
end


text:
function clean(phone_number)
    # Main machinery is in the regex expression.
    NXX = "[2-9][\\d]{2}"
    country = "(?:\\+?1)?"
    area = "\\(?\\s*($NXX)\\s*\\)?"
    exchange = "($NXX)"
    subscriber = "([\\d]{4})"
    fillers = "\\s*(?:\\.|-)?\\s*"
    r = Regex("^\\s*$country$fillers$area$fillers$exchange$fillers$subscriber\\s*\$")

    result = match(r, phone_number)

    if result != nothing
        result = Base.string(result.captures...)
    end

    return result
end


text:
module ParticleAccelerations

export total_acceleration, total_acceleration!,
    AbstractNBodySystem, AbstractInteraction,
    Pairwise, PairwiseSymmetric

using .Threads
using DiffEqBase

include("api.jl")
include("diffeq.jl")

function total_acceleration(system)
    dv = zero(atype(system))
    total_acceleration!(dv, system)
    return dv
end

@inline function total_acceleration!(dv, system)
    inters = interactions(system)
    for i in inters
        interaction_sum!(dv, summation_alg(i), issparse(i), i, system)
    end
end

@inline interaction_sum!(dv, alg, ::Val{false}, i, sys) = interaction_sum!(dv, alg, i, sys)

interaction_sum!(dv, alg, ::Val{true}, i, sys) = force!(dv, i, sys)

function interaction_sum!(dv, ::PairwiseSymmetric{false}, interaction, system)
    n = number_of_particles(system, interaction)
    @inbounds for i in Base.OneTo(n)
        for j in Base.OneTo(i-1)
            force!(dv, interaction, system, i, j)
        end
    end

    dv ./= particle_mass.(Ref(system), Base.OneTo(n))
end

function interaction_sum!(dv, ::PairwiseSymmetric{true}, interaction, system)
    n = number_of_particles(system, interaction)
    dv_thread = [zero(dv) for _ in nthreads()]
    @threads for i in Base.OneTo(n)
        for j in Base.OneTo(i-1)
            force!(dv_thread[threadid()], interaction, system, i, j)
        end
        dv_thread[threadid()][i] /= particle_mass(system, i)
    end
    sum!(dv, dv_thread)
end

function interaction_sum!(dv, ::PairList{false}, interaction, system)
    n = number_of_particles(system, interaction)
    list = pairlist(system, interaction)
    for ni in Base.OneTo(n)
        i, j = list[ni]
        force!(dv, interaction, system, i, j)
        dv[i] /= particle_mass(system, i)
    end
end

function interaction_sum!(dv, ::PairList{true}, interaction, system)
    n = number_of_particles(system, interaction)
    list = pairlist(system, interaction)
    dv_thread = [zero(dv) for _ in nthreads()]
    @threads for ni in Base.OneTo(n)
        i, j = list[ni]
        force!(dv_thread[threadid()], interaction, system, i, j)
        dv_thread[threadid()][i] /= particle_mass(system, i)
    end
end

end


text:
const Shift = 1
const ReduceL = 2
const ReduceR = 3

type State
    tokens::Vector{Token}
    step::Int
    top::Int
    l::State
    r::Int
    lc::State
    rc::State

    State() = new()
    State(tokens) = new(tokens, 1, 1, State(), 2, State(), State())
    State(tokens, step, top, l, r, lc, rc) = new(tokens, step, top, l, r, lc, rc)
end

left(s::State) = isdefined(s, :tokens) ? s.l : nothing

Base.isvalid(s::State) = isdefined(s, :tokens)
Base.done(s::State) = s.step >= length(s.tokens) * 2
Base.length(s::State) = s.step

function Base.next(st::State, act::Int)
    if act == Shift
        State(st.tokens, st.step+1, st.r, st, st.r+1, State(), State())
    elseif act == ReduceL
        State(st.tokens, st.step+1, st.top, st.l.l, st.r, st.l, st.rc)
    elseif act == ReduceR
        State(st.tokens, st.step+1, st.l.top, st.l.l, st.r, st.l.lc, st)
    else
        throw("Invalid action: $(act)")
    end
end

function Base.next(st::State)
    tokens = st.tokens
    nexts = Tuple{State,Int}[]
    if tokens[1].head < 0
        st.r <= length(st.tokens) && push!(nexts, next(st,Shift))
        isvalid(st.l) && push!(nexts, next(st,ReduceL), next(st,ReduceR))
        nexts
    else
        if isvalid(st.l)
            s0, s1 = tokens[st.top], tokens[st.l.top]
            if s1.head == st.top
                act = ReduceL
            elseif s0.head == st.l.top
                reducable = all(i -> tokens[i].head != st.top, st.r:length(tokens))
                act = reducable ? ReduceR : Shift
            elseif st.r <= length(tokens)
                act = Shift
            else
                throw("Invalid")
            end
        else
            act = Shift
        end
        [next(st,act)]
    end
end


text:
t=Trie{Int}()

t["amy"]=56
t["ann"]=15
t["emma"]=30
t["rob"]=27
t["roger"]=52

@test haskey(t, "roger")
@test get(t,"rob",nothing) == 27
@test sort(keys(t)) == ["amy", "ann", "emma", "rob", "roger"]
@test t["rob"] == 27
@test sort(keys_with_prefix(t,"ro")) == ["rob", "roger"]

# constructors
ks = ["amy", "ann", "emma", "rob", "roger"]
vs = [56, 15, 30, 27, 52]
kvs = collect(zip(ks, vs))
@test isa(Trie(ks, vs), Trie{Int})
@test isa(Trie(kvs), Trie{Int})
@test isa(Trie(Dict(kvs)), Trie{Int})
@test isa(Trie(ks), Trie{Void})

# path iterator
t0 = t
t1 = t0.children['r']
t2 = t1.children['o']
t3 = t2.children['b']
@test collect(path(t, "b")) == [t0]
@test collect(path(t, "rob")) == [t0, t1, t2, t3]
@test collect(path(t, "robb")) == [t0, t1, t2, t3]
@test collect(path(t, "ro")) == [t0, t1, t2]
@test collect(path(t, "roa")) == [t0, t1, t2]


text:
# cell.jl

# This is the main script, to run this script type: include("cell.jl") into Julia
# Set up the following scripts to define your problem and solution:
# (1) genetics.jl -> This is your system, feel free to design or evolve an ideal genetic system
# (2) input.jl -> This defines the problem you are trying to solve in terms of a flow of data into the cell
# (3) output.jl -> This defines the conditions by which you monitor success/failure/end conditions for the system

#Set a random seed to ensure consistent results
srand(0)

###### Initial Set Up ######
println("Creating Cell")

# Set Up Data Types

type Molecule
	# data stored in molecule
	data
end

type Gene
	# tree structure
	parent
	# protein function
	protein
	# length of protein life
	access
	# name of protein it creates
	lifetime
end

type ProteinAbundance
	# Gene for Protein
	gene
	# Raw Number (not normalized)
	abundance
end

include("output.jl")

include("genetics.jl")

include("input.jl")

# (1) Inputter
# This gets molecules for proteins on demand using a specified function from include.jl

pool = Molecule[]

function fromFactory()
	while true
		produce(getMolecule())
	end
end

function fromEnvironment(outside, permeability)
	while true
		r = rand()
		if r < permeability
			produce(consume(outside))
		else
			if length(fetch(pool)) > 0
				println("Processing Molecule from Pool")
				produce(shift!(fetch(pool)))
			else
				produce(consume(outside))
			end
		end
	end
end

outside = @task fromFactory()
environment = @task fromEnvironment(outside, 0.5)

# (2) Consumer

# This assigns molecules to proteins based on distribution
# Using Proteins to do (terrible) things to molecules

function consumer()
	println("Consume a Molecule")
	result = processMolecule(consume(environment), proteinDistribution)
	#println("")
	#println(pool)
	#println("")
	append!(pool, result[1])
	return (result[2], result[3])
end

function processMolecule(molecule, pd)

	# 1 - Normalised Protein Distribution
	normPD = ProteinAbundance[]
	sumPD = 0

	for p in pd
		sumPD = sumPD + p.abundance
	end

	for p in pd
		normPA = ProteinAbundance(p.gene, (p.abundance/sumPD))
		push!(normPD,normPA)
	end

	# 2 - Distribute
	handOfGod = rand()
	sumNPD = 0
	# Bad Default
	assignedProtein = root

	for p in normPD
		sumNPD = sumNPD + p.abundance
		if (sumNPD > handOfGod)
			assignedProtein = p.gene
			break
		end
	end

	# 3 - Run Process

	return (assignedProtein.protein(molecule))

end

###### Main Loop ######
println("Cell is Alive")

complete = false
#outputCount = 0
while complete == false
	# Run it like its hot
	result = consumer()
	outputCount = outputCount + result[2]
	backreg = result[1]
	complete = checkEndCon()
end

println(outputCount)
###### Final Steps ######
println("Cell is Dead")

max_stars_repo_name: voutasaurus/metalbolicalgorithm

max_issues_repo_name: voutasaurus/metalbolicAlgorithm

max_forks_repo_name: voutasaurus/metalbolicAlgorithm


text:
################################################################################
## Loss utility functions

function vector_mse(x, x̂)
    res = zero(eltype(x[1]))
    @inbounds for i in eachindex(x)
        res += sum((x[i] .- x̂[i]).^2)
    end
    # divide per number of time steps and batch size
    res /= size(x, 1) * size(x[1], 2)
    return res
end

# KL divergence
kl(μ, logσ²) = (exp(logσ²) + μ^2 - logσ² - 1) / 2

function vector_kl(μ::T, logσ²::T) where T <: Tuple{Matrix, Matrix}
    P = eltype(μ[1])
    s = zero(P)
    # go through initial conditions and parameters
    @inbounds for i in 1:2
        s1 = zero(P)
        @inbounds for k in eachindex(μ[i])
            s1 += kl(μ[i][k], logσ²[i][k])
        end
        # divide per batch size
        s1 /= size(μ[i], 2)
        s += s1
    end
    return s
end

function vector_kl(μ::T, logσ²::T) where T <: Matrix
    P = eltype(μ)
    s = zero(P)
    # go through initial conditions
    @inbounds for k in eachindex(μ)
        s += kl(μ[k], logσ²[k])
    end
    # divide per batch size
    s /= size(μ, 2)
    return s
end

function vector_kl(μ::T, logσ²::T) where T <: Tuple{CuArray, CuArray}
	return sum( [ mean(sum(kl.(μ[i], logσ²[i]), dims=1)) for i in 1:length(μ) ] )
end

## annealing factor scheduler
# based on https://github.com/haofuml/cyclical_annealing
function frange_cycle_linear(n_iter, start::T=0.0f0, stop::T=1.0f0,  n_cycle=4, ratio=0.5) where T
    L = ones(n_iter) * stop
    period = n_iter/n_cycle
    step = T((stop-start)/(period*ratio)) # linear schedule

    for c in 0:n_cycle-1
        v, i = start, 1
        while (v ≤ stop) & (Int(round(i+c*period)) < n_iter)
            L[Int(round(i+c*period))] = v
            v += step
            i += 1
        end
    end
    return T.(L)
end

################################################################################
## Data pre-processing

function normalize_to_unit_segment(X)
    min_val = minimum(X)
    max_val = maximum(X)

    X̂ = (X .- min_val) ./ (max_val - min_val)
    return X̂, min_val, max_val
end

denormalize_unit_segment(X̂, min_val, max_val) = X̂ .* (max_val .- min_val) .+ min_val

################################################################################
## Training help function

function time_loader(x, full_seq_len, seq_len)

    x_ = Array{Float32, 3}(undef, (size(x,1), size(x,2), seq_len))
    idxs = rand_time(full_seq_len, seq_len)
    for i in 1:size(x,2)
        x_[:,i,:] = x[:, i, idxs]
    end
    return x_
end

function rand_time(full_seq_len, seq_len)
    start_time = rand(1:full_seq_len - seq_len)
    idxs = start_time:start_time+seq_len-1
    return idxs
end

max_issues_repo_name: gabrevaya/LatentDiffEq

max_forks_repo_name: gabrevaya/LatentDiffEq


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_events

using AWS.Compat
using AWS.UUIDs
"""
    CreateDetectorModel()

Creates a detector model.

# Required Parameters
- `detectorModelDefinition`: Information that defines how the detectors operate.
- `detectorModelName`: The name of the detector model.
- `roleArn`: The ARN of the role that grants permission to AWS IoT Events to perform its operations.

# Optional Parameters
- `detectorModelDescription`: A brief description of the detector model.
- `evaluationMethod`: Information about the order in which events are evaluated and how actions are executed.
- `key`: The input attribute key used to identify a device or system to create a detector (an instance of the detector model) and then to route each input received to the appropriate detector (instance). This parameter uses a JSON-path expression in the message payload of each input to specify the attribute-value pair that is used to identify the device associated with the input.
- `tags`: Metadata that can be used to manage the detector model.
"""

create_detector_model(detectorModelDefinition, detectorModelName, roleArn; aws_config::AWSConfig=global_aws_config()) = iot_events("POST", "/detector-models", Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition, "detectorModelName"=>detectorModelName, "roleArn"=>roleArn); aws_config=aws_config)
create_detector_model(detectorModelDefinition, detectorModelName, roleArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("POST", "/detector-models", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition, "detectorModelName"=>detectorModelName, "roleArn"=>roleArn), args)); aws_config=aws_config)

"""
    CreateInput()

Creates an input.

# Required Parameters
- `inputDefinition`: The definition of the input.
- `inputName`: The name you want to give to the input.

# Optional Parameters
- `inputDescription`: A brief description of the input.
- `tags`: Metadata that can be used to manage the input.
"""

create_input(inputDefinition, inputName; aws_config::AWSConfig=global_aws_config()) = iot_events("POST", "/inputs", Dict{String, Any}("inputDefinition"=>inputDefinition, "inputName"=>inputName); aws_config=aws_config)
create_input(inputDefinition, inputName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("POST", "/inputs", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("inputDefinition"=>inputDefinition, "inputName"=>inputName), args)); aws_config=aws_config)

"""
    DeleteDetectorModel()

Deletes a detector model. Any active instances of the detector model are also deleted.

# Required Parameters
- `detectorModelName`: The name of the detector model to be deleted.

"""

delete_detector_model(detectorModelName; aws_config::AWSConfig=global_aws_config()) = iot_events("DELETE", "/detector-models/$(detectorModelName)"; aws_config=aws_config)
delete_detector_model(detectorModelName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("DELETE", "/detector-models/$(detectorModelName)", args; aws_config=aws_config)

"""
    DeleteInput()

Deletes an input.

# Required Parameters
- `inputName`: The name of the input to delete.

"""

delete_input(inputName; aws_config::AWSConfig=global_aws_config()) = iot_events("DELETE", "/inputs/$(inputName)"; aws_config=aws_config)
delete_input(inputName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("DELETE", "/inputs/$(inputName)", args; aws_config=aws_config)

"""
    DescribeDetectorModel()

Describes a detector model. If the version parameter is not specified, information about the latest version is returned.

# Required Parameters
- `detectorModelName`: The name of the detector model.

# Optional Parameters
- `version`: The version of the detector model.
"""

describe_detector_model(detectorModelName; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/detector-models/$(detectorModelName)"; aws_config=aws_config)
describe_detector_model(detectorModelName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/detector-models/$(detectorModelName)", args; aws_config=aws_config)

"""
    DescribeInput()

Describes an input.

# Required Parameters
- `inputName`: The name of the input.

"""

describe_input(inputName; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/inputs/$(inputName)"; aws_config=aws_config)
describe_input(inputName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/inputs/$(inputName)", args; aws_config=aws_config)

"""
    DescribeLoggingOptions()

Retrieves the current settings of the AWS IoT Events logging options.

"""

describe_logging_options(; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/logging"; aws_config=aws_config)
describe_logging_options(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/logging", args; aws_config=aws_config)

"""
    ListDetectorModelVersions()

Lists all the versions of a detector model. Only the metadata associated with each detector model version is returned.

# Required Parameters
- `detectorModelName`: The name of the detector model whose versions are returned.

# Optional Parameters
- `maxResults`: The maximum number of results to return at one time.
- `nextToken`: The token for the next set of results.
"""

list_detector_model_versions(detectorModelName; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/detector-models/$(detectorModelName)/versions"; aws_config=aws_config)
list_detector_model_versions(detectorModelName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/detector-models/$(detectorModelName)/versions", args; aws_config=aws_config)

"""
    ListDetectorModels()

Lists the detector models you have created. Only the metadata associated with each detector model is returned.

# Optional Parameters
- `maxResults`: The maximum number of results to return at one time.
- `nextToken`: The token for the next set of results.
"""

list_detector_models(; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/detector-models"; aws_config=aws_config)
list_detector_models(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/detector-models", args; aws_config=aws_config)

"""
    ListInputs()

Lists the inputs you have created.

# Optional Parameters
- `maxResults`: The maximum number of results to return at one time.
- `nextToken`: The token for the next set of results.
"""

list_inputs(; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/inputs"; aws_config=aws_config)
list_inputs(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/inputs", args; aws_config=aws_config)

"""
    ListTagsForResource()

Lists the tags (metadata) you have assigned to the resource.

# Required Parameters
- `resourceArn`: The ARN of the resource.

"""

list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/tags", Dict{String, Any}("resourceArn"=>resourceArn); aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("GET", "/tags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn), args)); aws_config=aws_config)

"""
    PutLoggingOptions()

Sets or updates the AWS IoT Events logging options. If you update the value of any loggingOptions field, it takes up to one minute for the change to take effect. If you change the policy attached to the role you specified in the roleArn field (for example, to correct an invalid policy), it takes up to five minutes for that change to take effect.

# Required Parameters
- `loggingOptions`: The new values of the AWS IoT Events logging options.

"""

put_logging_options(loggingOptions; aws_config::AWSConfig=global_aws_config()) = iot_events("PUT", "/logging", Dict{String, Any}("loggingOptions"=>loggingOptions); aws_config=aws_config)
put_logging_options(loggingOptions, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("PUT", "/logging", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("loggingOptions"=>loggingOptions), args)); aws_config=aws_config)

"""
    TagResource()

Adds to or modifies the tags of the given resource. Tags are metadata that can be used to manage a resource.

# Required Parameters
- `resourceArn`: The ARN of the resource.
- `tags`: The new or modified tags for the resource.

"""

tag_resource(resourceArn, tags; aws_config::AWSConfig=global_aws_config()) = iot_events("POST", "/tags", Dict{String, Any}("resourceArn"=>resourceArn, "tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("POST", "/tags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes the given tags (metadata) from the resource.

# Required Parameters
- `resourceArn`: The ARN of the resource.
- `tagKeys`: A list of the keys of the tags to be removed from the resource.

"""

untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = iot_events("DELETE", "/tags", Dict{String, Any}("resourceArn"=>resourceArn, "tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("DELETE", "/tags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateDetectorModel()

Updates a detector model. Detectors (instances) spawned by the previous version are deleted and then re-created as new inputs arrive.

# Required Parameters
- `detectorModelDefinition`: Information that defines how a detector operates.
- `detectorModelName`: The name of the detector model that is updated.
- `roleArn`: The ARN of the role that grants permission to AWS IoT Events to perform its operations.

# Optional Parameters
- `detectorModelDescription`: A brief description of the detector model.
- `evaluationMethod`: Information about the order in which events are evaluated and how actions are executed.
"""

update_detector_model(detectorModelDefinition, detectorModelName, roleArn; aws_config::AWSConfig=global_aws_config()) = iot_events("POST", "/detector-models/$(detectorModelName)", Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition, "roleArn"=>roleArn); aws_config=aws_config)
update_detector_model(detectorModelDefinition, detectorModelName, roleArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("POST", "/detector-models/$(detectorModelName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition, "roleArn"=>roleArn), args)); aws_config=aws_config)

"""
    UpdateInput()

Updates an input.

# Required Parameters
- `inputDefinition`: The definition of the input.
- `inputName`: The name of the input you want to update.

# Optional Parameters
- `inputDescription`: A brief description of the input.
"""

update_input(inputDefinition, inputName; aws_config::AWSConfig=global_aws_config()) = iot_events("PUT", "/inputs/$(inputName)", Dict{String, Any}("inputDefinition"=>inputDefinition); aws_config=aws_config)
update_input(inputDefinition, inputName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events("PUT", "/inputs/$(inputName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("inputDefinition"=>inputDefinition), args)); aws_config=aws_config)


text:
println("Aa) Test of the radial integration in CI computations (Randolf).")

function printList(list)
  for element in list
    println(element)
  end
end

function getOrbital(orbitalList, subshell)
  for orbital in orbitalList
    if(subshell.n == orbital.subshell.n && subshell.kappa == orbital.subshell.kappa)
      return orbital
    end
  end

  throw(ErrorException("Orbital not found"))

end

nucleus = Nuclear.Model(10.)

grid = JAC.Radial.Grid("grid: by given parameters", rnt=2.0e-6, h=5.0e-2, NoPoints=390)
setDefaults("standard grid", grid)

# configurations = [Configuration("1s^2 2s^2 2p^6")];      orbitalsIn = JAC.read("orbital list: Grasp92", "Work/Ne-jac/Ne-0+-1-scf.out")
# configurations = [Configuration("1s^2 2s^2 2p^5 3d")];   orbitalsIn = JAC.read("orbital list: Grasp92", "Work/Ne-jac/Ne-0+-2-scf.out")
configurations = [Configuration("1s^2 2s^2 2p^5")];      orbitalsIn = JAC.read("orbital list: Grasp92", "Work/Ne-jac/Ne-0+-3-scf.out")
# configurations = [Configuration("1s^2 2s^2 2p^5 3p")];   orbitalsIn = JAC.read("orbital list: Grasp92", "Work/Ne-jac/Ne-0+-4-scf.out")

# configurations = [Configuration("3s^2 3p^4 3d")]
println("Start from the non-relativistic configurations:")
printList(configurations)

# orbitals1 = JAC.read("orbital list: Grasp92", "Ne-0+-scf.exp.out")

relativisticConfigurationsList = ConfigurationR[]
for  configuration in configurations
    wa = JAC.generate("configuration list: relativistic", configuration)
    append!( relativisticConfigurationsList, wa)
end

println("Have the relativistic configurations:")
printList(relativisticConfigurationsList)

subshellList = JAC.generate("subshells: ordered list for relativistic configurations", relativisticConfigurationsList)
println("The subshells are: $subshellList")

setDefaults("relativistic subshell list", subshellList)

# Generate the relativistic CSF's for the given subshell list
csfList = CsfR[]
for  configuration in relativisticConfigurationsList
    newCsfs = generate("CSF list: from single ConfigurationR", configuration, subshellList)
    append!( csfList, newCsfs)
end
##x for  i = 1:length(csfList)    println("perform-ac: ", csfList[i])    end

println("Have the following $(length(csfList)) CSFs")
printList(csfList)

# Determine the number of electrons and the list of coreSubshells
NoElectrons      = sum( csfList[1].occupation )
coreSubshellList = Subshell[]
for  k in 1:length(subshellList)
    mocc = JAC.subshell_2j(subshellList[k]) + 1;    is_filled = true
    for  csf in csfList
        if  csf.occupation[k] != mocc    is_filled = false;    break   end
    end
    if   is_filled    push!( coreSubshellList, subshellList[k])    end
end
println("The number of electron is $NoElectrons")
println("The core subshells are: $coreSubshellList")

orbitals = Dict{Subshell, Orbital}()
for  subshell in subshellList
  try
    oin = getOrbital(orbitalsIn, subshell)
    orbitals = Base.merge( orbitals, Dict( subshell => oin) )
  catch err
    println("Not all orbitals found in input file, missing: $subshell")
    exit()
  end
end

# println(orbitals)

basis = (true, NoElectrons, subshellList, csfList, coreSubshellList, orbitals)

multiplet = perform("computation: CI", basis, nucleus, grid, MultipletSettings())


text:
using IRKGaussLegendre, DiffEqProblemLibrary, DiffEqDevTools, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_bigfloat2Dlinear

function NbodyODE!(du,u,Gm,t)
     N = length(Gm)
     du[1,:,:] .= 0
     for i in 1:N
        qi = u[2,:,i]
        Gmi = Gm[i]
        du[2,:,i] = u[1,:,i]
        for j in (i+1):N
           qj = u[2,:,j]
           Gmj = Gm[j]
           qij = qi - qj
           auxij = (qij[1]*qij[1]+qij[2]*qij[2]+qij[3]*qij[3])^(-3/2)
           du[1,:,i] -= Gmj*auxij*qij
           du[1,:,j] += Gmi*auxij*qij
        end
     end

    return
end

Gm = [5, 4, 3]
N=length(Gm)
q=[1,-1,0,-2,-1,0,1,3,0]
v=zeros(size(q))
q0 = reshape(q,3,:)
v0 = reshape(v,3,:)
u0 = Array{Float64}(undef,2,3,N)
u0[1,:,:] = v0
u0[2,:,:] = q0
tspan = (0.0,63.0)
prob=ODEProblem(NbodyODE!,u0,tspan,Gm)
sol1=solve(prob,IRKGL16(),reltol=1e-12, abstol=1e-12)

# Analytical tests

sol = solve(prob_ode_2Dlinear,IRKGL16())
@test sol.errors[:l2] < 1e-16

dts = (1//2) .^ (4:-1:2)
sim = test_convergence(dts,prob_ode_bigfloat2Dlinear,IRKGL16())
@test abs(sim.𝒪est[:final]-16) < 0.5


text:
"A grid corresponding to an unstructured collection of points."
struct ScatteredGrid{T} <: AbstractGrid{T,1}
    points     ::  Vector{T}
    domain     ::  Domain
    ScatteredGrid(points::Vector{T}, domain=DomainSets.FullSpace{T}()) where T =
        new{T}(points, domain)
end

name(grid::ScatteredGrid) = "Scattered grid"

size(g::ScatteredGrid) = (length(g.points),)

unsafe_getindex(g::ScatteredGrid, idx) = g.points[idx]

support(g::ScatteredGrid) = g.domain


text:
using Base.Test
using Bebop

const duration = 2

const sampRate = 44100
# the A above middle C
const musicFreq = 440.0

const t = [1:duration * sampRate] / sampRate
const xmono = float32(sin(2 * pi * musicFreq * t))
const xstereo = [xmono'; xmono']

initialize() do
    openStream(sampRate, 1, Float32) do stream
        writeStream(stream, xmono)
    end
    openStream(sampRate, 2, Float32) do stream
        writeStream(stream, xstereo)
    end
end


text:
# JuliaCon 2019 Demonstration
if (isdir("Z:"))
    include("Z:\\Ultimate Datetime Datatype\\Ultimate Datetime Git\\Ultimate-Datetime\\Julia Integration\\UltimateDT.jl")
elseif (isdir("C:"))
    include("C:\\Users\\jay\\Documents\\Ultimate Datetime Datatype\\Ultimate Datetime\\Julia Integration\\UltimateDT.jl")
end
# Create a struct to contain code fragments
struct CodeFragment
    description::String
    expression::Expr
end
# Create a function to describe then demonstrate a code fragment
function demo_code(codeFragment::CodeFragment)
    if (codeFragment.description != "")
        printstyled(stdout, "\n *** ", codeFragment.description, " ***\n", bold=true, color=:light_green)
    end
    # The replace functions removes comment lines inserted by the system
    print(stdout, replace(string(codeFragment.expression), r"#.*#\n" => ""), "  ")
    if eval(codeFragment.expression) != nothing
        printstyled(stdout, eval(codeFragment.expression), color=:light_cyan)
    end
    readline() # Wait for a carraige return
    nothing
end
# Use the Unicode plotting package
using UnicodePlots
# Create a function to create the vector of vectors required to generate a boxplot
#   comparing two UncertainFloat64s.
function boxCompareUF(uf1::UncertainFloat64, uf2::UncertainFloat64)
    val1 = uf1.value - uf1.precision*uf1.uncertainty
    val2 = uf1.value - uf1.precision*uf1.uncertainty
    val3 = uf1.value
    val4 = uf1.value + uf1.precision*uf1.uncertainty
    val5 = uf1.value + uf1.precision*uf1.uncertainty
    val6 = uf2.value - uf2.precision*uf2.uncertainty
    val7 = uf2.value - uf2.precision*uf2.uncertainty
    val8 = uf2.value
    val9 = uf2.value + uf2.precision*uf2.uncertainty
    val10 = uf2.value + uf2.precision*uf2.uncertainty
    return [[val1,val2,val3,val4,val5],[val6,val7,val8,val9,val10]]
end
# Create an array with strings for all of the comparison operators
compops = ["==","!=","==′","!=′","==″","!=″","==‴","!=‴",">","<=″",
           ">′","<=′",">″","<=",">=","<″",">=′","<′",">=″","<"]
# Create an array with descriptions of the operators and print them in boxAndCompare
compopsDesc = ["Must equal", "May not equal", "Has same midpoint as", "Does not have same midpoint as",
               "May equal", "Cannot equal", "Has same range as", "Does not have same range as",
               "Must be greater than", "Can be less than or equal to", "Has larger midpoint than",
               "Has midpoint less than or equal to that of", "Can be greater than", "Must be less than or equal to",
               "Must be greater than or equal to", "Can be less than", "Has midpoint greater than or equal to that of",
               "Has smaller midpoint than", "Can be greater than or equal to", "Must be less than"
               ]
# Create a function to draw a box plot and display the results of the comparisons
# uf1 and uf2 need to be in the global scope so that the eval in boxAndCompare receives the proper values
uf1 = UncertainFloat64(1.,p=1,u=1)
uf2 = uf1
function boxAndCompare(uf1t::UncertainFloat64, uf2t::UncertainFloat64)
    # Need to pass the global variables to eval
    global uf1 = uf1t
    global uf2 = uf2t
    # When used in an expression to be evaluated, boxplot works if it is the last line in
    #  the expression.  Otherwise, must use show(boxplot).
    show(boxplot([string(uf1), string(uf2)], boxCompareUF(uf1,uf2), border=:none))
    print(stdout,"\n")
    for i = 1:20
        code = Meta.parse("uf1" * compops[i] * "uf2")
        print(stdout, "\n", rpad(string(code),11), "   ")
        printstyled(stdout, rpad(compopsDesc[i],45), "  ", color=:light_green)
        if (eval(code) == true)
            printstyled(stdout, eval(code), color=:light_cyan)
        else
            printstyled(stdout, eval(code), color=:light_red)
        end
    end
    return nothing
end
# Create a function to do a simple comparison of two LocalCalCoordsDT
lcc1 = LocalCalCoordsDT(7, 23, 2019, "America/New_York", 2; cal=0, h=10)
lcc2 = lcc1
simplecomps = [:(lcc1 == lcc2), :(lcc1 != lcc2), :(lcc1 > lcc2), :(lcc1 <= lcc2), :(lcc1 < lcc2), :(lcc1 >= lcc2)]
function simpleCompare(lcc1t::LocalCalCoordsDT, lcc2t::LocalCalCoordsDT)
    # Need to pass the global variables to eval
    global lcc1 = lcc1t
    global lcc2 = lcc2t
    print(stdout,"\n")
    for i=1:6
        code = simplecomps[i]
        print(stdout, "\n", rpad(string(code),12), "   ")
        if (eval(code) == true)
            printstyled(stdout, eval(code), color=:light_cyan)
        else
            printstyled(stdout, eval(code), color=:light_red)
        end
    end
end
# Create an array to contain the code fragments
code = Array{CodeFragment,1}()
# Uncertain Floats
# Create an Uncertain Float with default precision and uncertainty
push!(code,CodeFragment("Uncertain float with default precision and uncertainty",:(UncertainFloat64(1005.25))))
# Show change in default precision and uncertainty with magnitude
push!(code,CodeFragment("Default precision and uncertainty change with magnitude",:(UncertainFloat64(5005.25))))
# Create 2 uncertain floats to demonstrate arithmetic functions
push!(code,CodeFragment("Create uf1 with specified precision and uncertainty",:(uf1 = UncertainFloat64(10,p=1,u=2))))
push!(code,CodeFragment("Create uf2 with specified precision and uncertainty",:(uf2 = UncertainFloat64(8,p=1,u=1))))
push!(code,CodeFragment("Addition - add absolute uncertainties, use minimum absolute precision",:(uf1 + uf2)))
push!(code,CodeFragment("Subtraction - add absolute uncertainties, use minimum absolute precision",:(uf1 - uf2)))
push!(code,CodeFragment("Create uf3 with specified precision and uncertainty",:(uf3 = UncertainFloat64(100,p=1,u=2))))
push!(code,CodeFragment("Multiplication - add fractional uncertainties, use minimum fractional precision",:(uf3 * uf2)))
push!(code,CodeFragment("Division - add fractional uncertainties, use minimum fractional precision",:(uf3 / uf2)))
# Exact values have a precision of 1 and uncertainy of 0 (i.e., as Integers)
push!(code,CodeFragment("Create uf1 with specified precision and uncertainty",:(uf1 = UncertainFloat64(1005.25,p=.01,u=2))))
push!(code,CodeFragment("Create uf2 as an exact value",:(uf2 = UncertainFloat64(4,p=1,u=0))))
push!(code,CodeFragment("Multiply the 2 values.  Uf1 is only contributor to p and u.",:(uf1 * uf2)))
# Compare with integber result
push!(code,CodeFragment("Multiply an uncertain value by an integer.  Uf1 is only contributor to p and u.",:(4 * uf1)))
# Zero values require special handling
push!(code,CodeFragment("Zero value with an uncertainty of 1.",:(uf1 = UncertainFloat64(0,p=1,u=1))))
push!(code,CodeFragment("Zero value with an uncertainty of 1 billion",:(uf2 = UncertainFloat64(0,p=1e9,u=1))))
push!(code,CodeFragment("Non-zero value with a 90% uncertainty",:(uf3 = UncertainFloat64(10.,p=1,u=9))))
push!(code,CodeFragment("Multiplication multiplies the uncertainty by the largest possible value of the other factor.",:(uf3 * uf1)))
push!(code,CodeFragment("Division increases the uncertainty by (1 + fractional uncertainty of the denominator).",:(uf2 / uf3)))
# Demonstrate range comparisons
push!(code,CodeFragment("Range based comparison operators with overlapping values",
    :(boxAndCompare(UncertainFloat64(8.,p=1,u=2),UncertainFloat64(9.,p=1,u=2)))))
push!(code,CodeFragment("Range based comparison operators with exact values",
    :(boxAndCompare(UncertainFloat64(8,p=1,u=0),UncertainFloat64(8,p=1,u=0)))))
push!(code,CodeFragment("Range based comparison operators with disjoint values",
    :(boxAndCompare(UncertainFloat64(5,p=1,u=2),UncertainFloat64(10,p=1,u=2)))))
push!(code,CodeFragment("Range based comparison operators with a single overlapping value",
    :(boxAndCompare(UncertainFloat64(5,p=1,u=2),UncertainFloat64(9,p=1,u=2)))))
push!(code,CodeFragment("Range based comparison operators with the same range",
    :(boxAndCompare(UncertainFloat64(9,p=1,u=2),UncertainFloat64(9,p=1,u=2)))))
# Demonstrate that relative time is different from UTCDatetime
push!(code,CodeFragment("Create a UTCDatetime with default precision and uncertainty.",
    :(utc1 = UTCDatetime(y=2019,m=7,d=23,h=11,min=45))))
push!(code,CodeFragment("Create a 2nd UTCDatetime with default precision and uncertainty.",
    :(utc2 = UTCDatetime(y=2019,m=7,d=23,h=10,min=45))))
# UTC Datetime Formatting
push!(code,CodeFragment("Formatting is automatic based on precision.",
    :(for i = -18:15;
        utc = UTCDatetime(gy=99,y=999999999,m=12,d=31,h=23,min=59,s=59,ns=999999999,as=999999999,p=i,u=1);
        printstyled(stdout,"\n", utc,color=:light_cyan);
      end)))
push!(code,CodeFragment("Subtract the 2 UTCDatetimes to produce a relative datetime.",:(rel1 = utc1 - utc2)))
push!(code,CodeFragment("Relative datetimes can be positive or negative.",:(rel2 = utc2 - utc1)))
# Show allowed operations on relative datetimes
push!(code,CodeFragment("Create a UTCDatetime with nanosecond precision and uncertainty.",
    :(utc1 = UTCDatetime(y=2019,m=7,d=23,h=11,min=45,p=-9,u=1))))
push!(code,CodeFragment("Create a 2nd UTCDatetime with nanosecond precision and uncertainty.",
    :(utc2 = UTCDatetime(y=2019,m=7,d=23,h=10,min=45,p=-9,u=1))))
push!(code,CodeFragment("Subtract the 2 UTCDatetimes to produce a relative datetime.",:(rel1 = utc1 - utc2)))
push!(code,CodeFragment("Relative datetimes can be multiplied by integers and floating point numbers.",:(2*rel1 + 1.5*rel1)))
push!(code,CodeFragment("A relative datetime can be divided by another relative datetime.",:( x = rel1 / (2*rel1 + 1.5*rel1))))
push!(code,CodeFragment("The result of this division is an UncertianFloat64 (which is the original reason the type was created.)!", :(typeof(x))))
push!(code,CodeFragment("When the UncertainFloat64 is multiplied by the denominator, the precision is restored with the propogated uncertainty.",
    :(x * (2*rel1 + 1.5*rel1))))
# Demonstrate Julian to Gregorian transition
push!(code,CodeFragment("Create a UTCDatetime just before the earliest transtion to the Gregorian calendar.",
    :(utcj = UTCDatetime(y=1582,m=10,d=4,h=23,min=0,p=0,u=1))))
push!(code,CodeFragment("Attempt to create a UTCDatetime during the skipped period.",
    :(utcs = UTCDatetime(y=1582,m=10,d=5,h=1,min=0,p=0,u=1))))
push!(code,CodeFragment("Create a UTCDatetime just after the earliest transition to the Gregorian calendar.",
    :(utcg = UTCDatetime(y=1582,m=10,d=15,h=1,min=0,p=0,u=1))))
push!(code,CodeFragment("Very little time elapsed between those two dates!",
    :(utcg - utcj)))
# Demonstrate leap seconds
push!(code,CodeFragment("Leap seconds occur at midnight UTC, either on June 30 or December 31.",
    :(lsminus1 = UTCDatetime(y=2015,m=6,d=30,h=23,min=59,s=59))))
push!(code,CodeFragment("This is the actual leap second.  Notice the seconds field.",
    :(ls = UTCDatetime(y=2015,m=6,d=30,h=23,min=59,s=60))))
push!(code,CodeFragment("This is the next second.",
    :(lsplus1 = UTCDatetime(y=2015,m=7,d=1,h=0,min=0,s=0))))
push!(code,CodeFragment("The arithmetic works.",
    :(lsplus1 - lsminus1)))
push!(code,CodeFragment("Only minutes containing leap seconds have a 60th second.",
    :(ls = UTCDatetime(y=2015,m=6,d=30,h=23,min=58,s=60))))
# Explore time zones
push!(code,CodeFragment("Prior to the late 19th century, time zones were based on local mean time, a type of solar time.",
    :(lchi = LocalCalCoordsDT(5, 15, 1870, "America/Chicago", 2; cal=0, h=4))))
push!(code,CodeFragment("The offset from UTC (GMT, at the time) was based solely on the longitude of a location.",
    :(translate(lchi,f=0))))
push!(code,CodeFragment("When this Chicago time is translated to a New York time, the difference is 54:34.4, not 1 hour!",
    :(lny = translate(lchi,tz = "America/New_York"))))
push!(code,CodeFragment("In the 20th century, we get the expected result.",
    :(lchi = LocalCalCoordsDT(5, 15, 1950, "America/Chicago", 2; cal=0, h=4))))
push!(code,CodeFragment("When this Chicago time is translated to a New York time, the difference is 1 hour.",
    :(lny = translate(lchi,tz = "America/New_York"))))
push!(code,CodeFragment("Daylight savings time began during World War I.  The Chicago time above is during daylight savings.",
    :(lny = translate(lchi,f=1))))
push!(code,CodeFragment("The transition to daylight savings time in 1950 took place on April 30th, the last Sunday in the month.",
    :(lchi = LocalCalCoordsDT(4, 30, 1950, "America/Chicago", 2; h=1, min=59))))
push!(code,CodeFragment("1:59 wall time is equivalent to 1:59 standard time.", :(translate(lchi,f=1))))
push!(code,CodeFragment("2:01 wall time does not exist!",
    :(lchi = LocalCalCoordsDT(4, 30, 1950, "America/Chicago", 2; h=2, min=1))))
push!(code,CodeFragment("2:01 standard time does exist...",
    :(lchi = LocalCalCoordsDT(4, 30, 1950, "America/Chicago", 1; h=2, min=1))))
push!(code,CodeFragment("... and is equivalent to 3:01 wall time.",:(translate(lchi,f=2))))
push!(code,CodeFragment("The transition back to standard time in 1950 took place on September 24th, the last Sunday in the month.",
    :(lchi = LocalCalCoordsDT(9, 24, 1950, "America/Chicago", 2; h=0, min=59))))
push!(code,CodeFragment("0:59 wall time is equivalent to 23:59 standard time.", :(translate(lchi,f=1))))
push!(code,CodeFragment("1:01 wall time is ambiguous, as it occurs twice.",
    :(lchi = LocalCalCoordsDT(9, 24, 1950, "America/Chicago", 2; h=1, min=1))))
push!(code,CodeFragment("1:01 standard time is not ambiguous...",
    :(lchi = LocalCalCoordsDT(9, 24, 1950, "America/Chicago", 1; h=1, min=1))))
push!(code,CodeFragment("... and is equivalent to 1:01 wall time, after the transition.",:(translate(lchi,f=2))))
push!(code,CodeFragment("Wall times must be disambiguated during ambiguous periods.",
    :(lchi = LocalCalCoordsDT(9, 24, 1950, "America/Chicago", 2; h=1, min=1, ba=1))))
push!(code,CodeFragment("1:01 before the transition is equivalent to 00:01 standard time.", :(translate(lchi,f=1))))
push!(code,CodeFragment("Daylight savings rules have changed over time.  Prior to 2007, March 15th was standard time in New York.",
    :(lny = LocalCalCoordsDT(3, 15, 2006, "America/New_York", 2; h=12))))
push!(code,CodeFragment("ie, standard and wall times were the same.", :(translate(lny,f=1))))
push!(code,CodeFragment("Beginning in 2007, daylight savings time starts the second Sunday in March, not the first Sunday in April.",
    :(lny = LocalCalCoordsDT(3, 15, 2007, "America/New_York", 2; h=12))))
push!(code,CodeFragment("Standard time is 1 hour less than wall time.", :(translate(lny,f=1))))
# Revisit calendars
push!(code,CodeFragment("Construct a Gregorian date in the proleptic America/New_York time zone.",
    :(lccg = LocalCalCoordsDT(10, 15, 1582, "America/New_York", 2; cal=0, h=10))))
push!(code,CodeFragment("The Julian calendar continued on - the last country to adopt the Gregorian calendar was Greece in 1923.",
    :(lccj = LocalCalCoordsDT(10, 5, 1582, "America/New_York", 2; cal=1, h=10))))
push!(code,CodeFragment("For every Gregorian datetime, there is a corresponding Julian datetime.",
    :(lccg == lccj)))
push!(code,CodeFragment("Simple comparisons work.", :(simpleCompare(lccg, lccj))))
push!(code,CodeFragment("Construct a later Gregorian date.",
    :(lccg2 = LocalCalCoordsDT(10, 16, 1582, "America/New_York", 2; cal=0, h=10))))
push!(code,CodeFragment("Compare.", :(simpleCompare(lccg2, lccj))))
push!(code,CodeFragment("Arithmetic works for this type, too.", :(lccg2 - lccj)))
# Show Swedish calendar
push!(code,CodeFragment("Sweden did something strange in 1712.",
    :(lccsw = LocalCalCoordsDT(2, 30, 1712, "America/New_York", 2; cal=2, h=10))))
push!(code,CodeFragment("The Julian calendar is 1 day ahead, since Sweden skipped 2-29-1700.",
    :(lccj = LocalCalCoordsDT(2, 29, 1712, "America/New_York", 2; cal=1, h=10))))
push!(code,CodeFragment("Arithmetic is calendar aware.", :(lccsw - lccj)))
push!(code,CodeFragment("When last we checked, the Gregorian calendar was 10 days ahead of the Julian calendar.",
    :(lccg = LocalCalCoordsDT(3, 10, 1712, "America/New_York", 2; cal=0, h=10))))
push!(code,CodeFragment("What the heck???", :(lccg - lccj)))
# Investigate the early atomic period
push!(code,CodeFragment("Let's take a deeper look at atomic time.  Create a datetime in the early atomic period.",
    :(lcce1 = LocalCalCoordsDT(6, 15, 1959, "America/New_York", 2; h=5, min=3))))
push!(code,CodeFragment("Create another datetime 1 second later.",
    :(lcce2 = LocalCalCoordsDT(6, 15, 1959, "America/New_York", 2; h=5, min=3, s=1))))
push!(code,CodeFragment("Huh! Have we discovered a bug in local datetimes?",
    :(lcce2 - lcce1)))
push!(code,CodeFragment("Let's try it with UTCDatetims.  Surely they have been better tested.",
    :(utce1 = UTCDatetime(m=6, d=15, y=1959, h=5, min=3, p=-9))))
push!(code,CodeFragment("Create another datetime 1 second later.",
    :(utce2 = UTCDatetime(m=6, d=15, y=1959, h=5, min=3, s=1, p=-9))))
push!(code,CodeFragment("Uh oh!  Wait, maybe it isn't a bug...",
    :(utce2 - utce1)))
push!(code,CodeFragment("Things changed in 1961.",
    :(utce1 = UTCDatetime(m=12, d=31, y=1960, h=23, min=59, s=59, p=-9))))
push!(code,CodeFragment("Create another datetime 1 second later.",
    :(utce2 = UTCDatetime(m=1, d=1, y=1961, h=0, min=0, s=0, p=-9))))
push!(code,CodeFragment("A jump!", :(utce2 - utce1)))
push!(code,CodeFragment("Create another datetime 1 second later.",
    :(utce3 = UTCDatetime(m=1, d=1, y=1961, h=0, min=0, s=1, p=-9))))
push!(code,CodeFragment("A longer second!", :(utce3 - utce2)))
push!(code,CodeFragment("Things changed, again, in 1968.",
    :(utce1 = UTCDatetime(m=12, d=31, y=1967, h=23, min=59, s=59, p=-9))))
push!(code,CodeFragment("Create another datetime 1 second later.",
    :(utce2 = UTCDatetime(m=1, d=1, y=1968, h=0, min=0, s=0, p=-9))))
push!(code,CodeFragment("Another jump!", :(utce2 - utce1)))
push!(code,CodeFragment("Create another datetime 1 second later.",
    :(utce3 = UTCDatetime(m=1, d=1, y=1968, h=0, min=0, s=1, p=-9))))
push!(code,CodeFragment("A slightly longer second!", :(utce3 - utce2)))
push!(code,CodeFragment("Things settled down in 1972.",
    :(utce1 = UTCDatetime(m=12, d=31, y=1971, h=23, min=59, s=59, p=-9))))
push!(code,CodeFragment("Create another datetime 1 second later.",
    :(utce2 = UTCDatetime(m=1, d=1, y=1972, h=0, min=0, s=0, p=-9))))
push!(code,CodeFragment("The last irregular jump.", :(utce2 - utce1)))
push!(code,CodeFragment("Create another datetime 1 second later.",
    :(utce3 = UTCDatetime(m=1, d=1, y=1972, h=0, min=0, s=1, p=-9))))
push!(code,CodeFragment("Seconds only last a second after the introduction of leap seconds.", :(utce3 - utce2)))
map(demo_code, code)
nothing

max_stars_repo_path: Julia Integration/JuliaCon 2019 Demonstration.jl

max_stars_repo_name: JayDweck/Ultimate-Datetime

max_issues_repo_name: JayDweck/Ultimate-Datetime

max_forks_repo_name: JayDweck/Ultimate-Datetime


text:
import Base: size

export dcdp_primitives,
    dcdp_Emax,
    flow,
    dcdp_tmpvars,
    check_EVgrad,
    check_EVgrad!,
    _zspace,
    _ψspace,
    _ψ1clamp,
    _nθt,
    _θt,
    _σv,
    print_summary

struct dcdp_primitives{FF,T<:Real,UP<:AbstractUnitProblem,AM<:AbstractMatrix{T},TT<:Tuple,AV<:AbstractVector{T}}
    β::T
    wp::UP            # structure of endogenous choice vars
    zspace::TT        # z-space (tuple)
    Πz::AM            # transition for z
    ψspace::AV        # ψspace = u + σv*v
    nθt::Int          # Num parameters in flow payoffs MINUS 1 for σv
    anticipate_e::Bool # do we anticipate the ϵ shocks assoc w/ each choice?
end

function dcdp_primitives(FF::Symbol, β::T, wp::UP, zspace::TT, Πz::AM, ψspace::AV, anticipate_e::Bool=true) where {T,TT,AM,AV,UP<:AbstractUnitProblem}
    dcdp_primitives{Val{FF},T,UP,AM,TT,AV}(β, wp, zspace, Πz, ψspace, number_of_model_parms(FF), anticipate_e)
end

flow(prim::dcdp_primitives{FF}) where {FF} = FF

# help us go from big parameter vector for all types to the relevant one
_σv(θ::AbstractVector) = last(θ)

# allows for adding 1 unit for σv
_θt(x::AbstractVector, nθt::Integer,          p1::Integer=0) = view(x, 1:nθt+p1)
_θt(x::AbstractVector, prim::dcdp_primitives, p1::Integer=0) = _θt(x, _nθt(prim), p1)

function print_summary(p::dcdp_primitives)
    print("\n--------------------------------------\n    Primitives\n--------------------------------------\n")
    println("Flow payoffs are $(flow(p)) with nθt+1 = $(_nθt(p))+1 parameters")
    println("Shocks anticipated? $(_anticipate_e(p))")
    println("β = $(p.β)")
    println("wp = $(p.wp)")
    println("zsize = $(length.(_zspace(p))) and zspace = $(_zspace(p))")
    println("nψ = $(_nψ(p)) and ψspace = $(_ψspace(p))")
    println("--------------------------------------")
end

# functions to retrieve elements from dcdp_primitives
_nθt(   prim::dcdp_primitives) = prim.nθt
_nz(    prim::dcdp_primitives) = size(prim.Πz,1)
_nψ(    prim::dcdp_primitives) = length(prim.ψspace) # prim.nψ
_nS(    prim::dcdp_primitives) = _nS(prim.wp)
_nSexp( prim::dcdp_primitives) = _nSexp(prim.wp)
_nd(    prim::dcdp_primitives) = _dmax(prim.wp)+1
_anticipate_e(prim::dcdp_primitives) = prim.anticipate_e

_zspace(prim::dcdp_primitives) = prim.zspace

# just in case we have things lying around...
_ψspace(prim::dcdp_primitives) = prim.ψspace
@inline _ψ1clamp(u::Real,v::Real,ρ::Real, prim::dcdp_primitives) = clamp(_ψ1(u,v,ρ), extrema(_ψspace(prim))...)

sprime_idx(prim::dcdp_primitives, i::Integer) = sprime_idx(prim.wp, i)
wp_info(prim::dcdp_primitives, i::Integer) = wp_info(prim.wp, i)
state(prim::dcdp_primitives, i::Integer) = state(prim.wp, i)
_dmax(prim::dcdp_primitives, i::Integer) = _dmax(prim.wp, i)
_dmax(prim::dcdp_primitives) = _dmax(prim.wp)

size(prim::dcdp_primitives) = _nz(prim), _nψ(prim), _nS(prim)

# --------------------- Emax --------------------------

struct dcdp_Emax{T<:Real,A1<:AbstractArray{T,3},A2<:AbstractArray{T,4}}
    EV::A1
    dEV::A2
    dEVσ::A1
end

dcdp_Emax(EV::AbstractArray3{T}, dEV::AbstractArray4{T}, dEVσ::AbstractArray3{T}) where {T} =  dcdp_Emax{T,typeof(EV),typeof(dEV)}(EV,dEV,dEVσ)

function dcdp_Emax(p::dcdp_primitives{FF,T}) where {FF,T}
    EV   = zeros(T, _nz(p), _nψ(p),          _nS(p))
    dEV  = zeros(T, _nz(p), _nψ(p), _nθt(p), _nS(p))
    dEVσ = zeros(T, _nz(p), _nψ(p),          _nSexp(p))
    dcdp_Emax(EV,dEV,dEVσ) # ,dEVψ)
end

# --------------------- check conformable --------------------------

function check_size(prim::dcdp_primitives, evs::dcdp_Emax)
    nz, nψ, nS = size(prim)
    nθ, nd = _nθt(prim), _nd(prim)
    nSexp = _nSexp(prim)

    prod(length.(_zspace(prim))) == nz ||  throw(error("Πz and zspace not compatible"))
    (nz,nψ,nS)      == size(evs.EV)    ||  throw(error("EV not conformable"))
    (nz,nψ,nθ,nS)   == size(evs.dEV)   ||  throw(error("dEV not conformable"))
    (nz,nψ,nSexp)   == size(evs.dEVσ) ||  throw(error("dEVσ not conformable"))
    # (nz,nψ,nSexp)   == size(evs.dEV_ψ) ||  throw(error("dEV_ψ not conformable"))
    return true
end

# --------------------- tmp vars --------------------------

struct dcdp_tmpvars{AA2<:AbstractMatrix{Float64}, AA3<:AbstractArray3{Float64}, AA4<:AbstractArray4{Float64}}
    ubVfull::AA3
    dubVfull::AA4
    dubV_σ::AA3

    q::AA3
    lse::Matrix{Float64}
    tmp::Matrix{Float64}
    tmp_cart::Matrix{CartesianIndex{3}}
    Πψtmp::Matrix{Float64}
    IminusTEVp::AA2
end

function check_size(prim::dcdp_primitives, t::dcdp_tmpvars)
    nz, nψ, nS = size(prim)
    nθ, nd = _nθt(prim), _nd(prim)
    nSexp =_nSexp(prim)

    # TODO: allow dmaxp1 to vary with regime

    (nz,nψ,nd)    == size(t.q)        || throw(DimensionMismatch())
    (nz,nψ,nθ,nd) == size(t.dubVfull) || throw(DimensionMismatch())
    (nz,nψ,nd)    == size(t.dubV_σ)   || throw(DimensionMismatch())
end

function dcdp_tmpvars(prim::dcdp_primitives)

    T = Float64
    nz, nψ, nS = size(prim)
    nd = _nd(prim)
    nθt = _nθt(prim)
    nSexp = _nSexp(prim)

    nθt > nψ &&  throw(error("Must have more length(ψspace) > length(θt)"))

    # choice-specific value functions
    ubVfull  = zeros(T,nz,nψ,nd)
    dubVfull = zeros(T,nz,nψ,nθt,nd)
    dubV_σ   = zeros(T,nz,nψ,nd)

    # other tempvars
    q        = zeros(T,nz,nψ,nd)
    lse      = zeros(T,nz,nψ)
    tmp      = zeros(T,nz,nψ)
    tmp_cart = Matrix{CartesianIndex{3}}(undef,nz,nψ)

    # transition matrices
    Πψtmp = Matrix{T}(undef,nψ,nψ)
    IminusTEVp = ensure_diagonal(prim.Πz)

    return dcdp_tmpvars(ubVfull, dubVfull, dubV_σ, q, lse, tmp, tmp_cart, Πψtmp, IminusTEVp)
end

function dcdp_tmpvars_view(t::dcdp_tmpvars, idxd::AbstractVector)
    first(idxd) == 1 || throw(DomainError())
    last(idxd) <= size(t.ubVfull,3) || throw(DomainError())

    @views ubV    = t.ubVfull[ :,:,  idxd]
    @views dubV   = t.dubVfull[:,:,:,idxd]
    @views dubV_σ = t.dubV_σ[  :,:,  idxd]
    @views q      = t.q[:,:,idxd]

    return dcdp_tmpvars(ubV, dubV, dubV_σ, q, t.lse, t.tmp, t.tmp_cart, t.Πψtmp, t.IminusTEVp)
end

function zero!(t::dcdp_tmpvars)
    zero!(t.ubVfull )
    zero!(t.dubVfull)
    zero!(t.dubV_σ  )
    zero!(t.q       )
    zero!(t.lse     )
    zero!(t.tmp     )
    zero!(t.Πψtmp   )
end

function zero!(evs::dcdp_Emax)
    zero!(evs.EV)
    zero!(evs.dEV)
    zero!(evs.dEVσ)
end

#


text:
module SSLTestMosaic
using SimpleSDMLayers
using Test
using Statistics

n1 = fill(1.0f0, 5, 5)
s1 = SimpleSDMResponse(n1, 0., 2., 0., 2.)

n2 = fill(2.0, 5, 5)
s2 = SimpleSDMResponse(n2, 1., 3., 1., 3.)

n3 = fill(3.0, 10, 10)
s3 = SimpleSDMResponse(n3, -1.8, 2.2, 1.4, 5.4)

n4 = fill(4.0, 5, 3)
s4 = SimpleSDMResponse(n4, 3.0, 4.2, 3.0, 5.0)

# Test
L = mosaic(mean, [s1, s2, s3, s4])

@test typeof(L) == SimpleSDMResponse{Float32}
@test extrema(L) == (1.0f0, 4.0f0)

end


text:
#Based on the book about Loop Quantum Gravity by Rodolfo Gambini and Jorge Pullin.

#Ch 2.
#1. The pole in the barn paradox. ✓
#Lorentz transformation:
gam=@equ γ=1/(1-v^2)	#c=1
tp=@equ t´=γ*(t-v*x)
xp=@equ x´=γ*(x-v*t)
println("A 3 meters long pole is traveling quickly toward a barn that is shorter than its length.")

println("Measured in the inertial system of the pole with l=length of pole and c=1,")
#println("At time 0:")
#a=[gam;@equs(x=0,t=0,v=0.5)];println("Base of pole: ", xp&a,",\t",tp&a)
#a=[gam;@equs(x=l,t=0,v=0.5)];println("Tip of pole: ", xp&a,",\t",tp&a)
println("at time t:")
a=[gam;@equs(x=0+t*v,v=0.5)];println("Base of pole: ", xp&a,",\t",tp&a)
a=[gam;@equs(x=l+t*v,v=0.5)];println("Tip of pole: ", xp&a,",\t",tp&a)
println("Time at the tip is earlier than at the base, hence when the tip touches the barn wall the base has had time to enter.")

println("The pole is rotated about the time axis of the barn. When one object is rotated relative to another their projections become dwarfed, like when drawing a line at right angle from the minute pointer of a clock to its hour pointer, the projection becomes progressively smaller until it disappears at ninety degrees which in this analogy is light speed. Photons can therefore be complete universes that are rotated at right angles with respect to this universe.")

#2
using LinearAlgebra
M=[:t^2 :t*:x1 :q :q;:x1*:t :x1^2 :q :q;:q :q :x2 :q;:q :q :q :x3]
d=@equ Ten(D,[v,w])=Ten($(M),[v,w]) #D is not mentioned in the texts...
η=@equ Ten(η,[μ,v])=Ten($(Matrix(Diagonal(([-1,1,1,1])))),[μ,v])
eq2_11=@equ Ten(Λ,[μ´,μ])=Ten([Cosh(ϕ) -Sinh(ϕ) 0 0;-Sinh(ϕ) Cosh(ϕ) 0 0;0 0 1 0;0 0 0 1], [μ´,μ])
d´=(eq2_11*(d&@equ(v=μ)))
d2´=(eq2_11*η*d)
#print(d´&@equ w=μ´) #should not be invariant
#print(d2´&@equ w=μ´) #should be invariant
#should make sense...

#3 ✓
a=@equs(Δt=0.91,Δx1=3.5,Δx2=5,Δx3=-1.54,ϕ=1.59,c=1)
eq2_8=@equ Δs^2=-(c*Δt)^2+(Δx1)^2+(Δx2)^2+(Δx3)^2
#2_8 Δs^2 should be invariant under 2_11
dX=@equ ΔX=Ten([Δt, Δx1, Δx2, Δx3], v)
@assert (η*dX)&@equ(μ=2) == (Ten(:η,Any[2,:v])*:ΔX ≖ :Δx1)

dX´=@equ(ΔX´=Ten(Λ,[μ´,μ])*Ten([Δt, Δx1, Δx2, Δx3], μ))&eq2_11
eq2_8´=@equ(Δs^2=-($(dX´.rhs)&@equ(μ´=1))^2+($(dX´.rhs)&@equ(μ´=2))^2+(Δx2)^2+(Δx3)^2)
#@assert isapprox(eq2_8.rhs&a, eq2_8´.rhs&a) #Nay, distance is no longer unchanged!
#eq2_8´=@equ(Δs^2=-(Δt*Cosh(ϕ)-Δx1*Sinh(ϕ))^2+(-Δt*Sinh(ϕ)+Δx1*Cosh(ϕ))^2+(Δx2)^2+(Δx3)^2) #substituting manually gives this

#4 ✓
emM=[0 -:E1 -:E2 -:E3;:E1 0 :B3 -:B2;:E2 -:B3 0 :B1;:E3 :B2 -:B1 0]
eq2_31=@equ Ten(F,[μ,v])=Ten($emM,[μ,v])
eq2_35=@equ Alt([s,m,n,l])*Ten([∂t,∂1,∂2,∂3],m)*Ten($emM,[n,l])=0
@assert (eq2_35/2)&@equ(s=1) == simplify(:∂1*:B1+:∂2*:B2+:∂3*:B3 ≖ 0.0) #Yay, ∂iBi!
@assert (-eq2_35/2)&@equ(s=2) == simplify(:∂t*:B1+:∂2*:E3+-1.0*:∂3*:E2 ≖ 0.0)

#n1=Ten(diagm([-1,1,1,1]),[:j,:j´])
#n2=Ten(diagm([-1,1,1,1]),[:k,:k´])
#simplify(@equ(2*Delta(i,m)=Alt([i,j,k])*($n1*$n2*Alt([j,k,m]))))


text:
using Documenter
using SpectralDistances, ControlSystems, Clustering, JuMP, Convex, SCS, GLPK, DSP, DynamicAxisWarping, SlidingDistancesBase

using Plots
plotly()

@info "makedocs"
makedocs(
    sitename = "SpectralDistances",
    # format = LaTeX(),
    format = Documenter.HTML(prettyurls = haskey(ENV, "CI")),
    modules = [SpectralDistances],
    pages = ["index.md", "ltimodels.md", "distances.md", "time.md", "interpolations.md", "plotting.md", "misc.md", "examples.md", "taskview.md"]
)

# For debugging purposes
# makedocs(
#     sitename = "SpectralDistances",
#     # format = LaTeX(),
#     format = Documenter.HTML(prettyurls = haskey(ENV, "CI")),
#     modules = [SpectralDistances],
#     pages = ["time.md"]
# )

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
@info "deploydocs"
deploydocs(
    repo = "github.com/baggepinnen/SpectralDistances.jl"
)
@info "done"


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "LibAOM"
version = v"1.0.0"
# Collection of sources required to build imagemagick

sources = [
    "http://www.andrews-corner.org/downloads/SBo/libaom-1.0.0.r1804.gcb43f766c.tar.gz" =>
    "365ebf37fdbe6265c4e2a6815c9d19b0caa0a1279a04e5b0eea8fc77e3847420",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mkdir aom_build
cd aom_build
apk install yasm
cmake $WORKSPACE/srcdir/libaom-1.0.0.r1804.gcb43f766c -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DBUILD_SHARED_LIBS=1
make -j${ncore}
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
      LibraryProduct(prefix, "addition", :addition),
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: ianshmean/LibAOMBuilder

max_issues_repo_name: ianshmean/LibAOMBuilder

max_forks_repo_name: ianshmean/LibAOMBuilder


text:
export ConvMap, LaplacianMap, DiffZMap, DiffXMap, GradDivMap
import DSP: hilbert, conv
import LinearMaps: LinearMap

"""
`ConvMap([T,] wavelet, n...) -> W`

Construct a convolution/correlation operator `W` which will act on an `AbstractVector`.

Parameters
----------
* `T` : `Type`, optional

`Type` of argument of `W`. Defaults to `Float64`.

* `wavelet`

One-dimensional wavelet.

* `n...`

Sequence of dimensions of the data.

Usage
-----
* Forward map and Adjoint maps

The forward map `W` convolves the last dimension of the input with `wavelet`.
The adjoint map `W` correlates the last dimension of the input with `wavelet`.

"""
ConvMap(T::Type, wav, n...) = LinearMap{T}(x -> wav_conv(x, wav, n...)[:],
                                           x -> wav_corr(x, wav, n...)[:], prod(n), prod(n))
ConvMap(wav, n...) = ConvMap(Float64, wav, n...)

"""
`LaplacianMap([T,] n...) -> Δ`

Construct a discretized Laplacian operator `Δ` which will act on an `AbstractVector`.

Parameters
----------
* `T` : `Type`, optional

`Type` of argument of `Δ`. Defaults to `Float64`.

* `n...`

Sequence of spatial dimensions of `Δ`.

Usage
-----
* Forward map and Adjoint maps

The forward map `Δ` multiplies a model vector of size `nz × nx × ny × ...` its second order derivative. It is symmetric.

Description
-----------

The forward and adjoint maps computes the following operation
```math
Δm = \\sum_{i} δ_i m
```
where
```math
δ_l m_{i,j,k,...} = - m_{...,l+1,...} + 2m_{...,l,...} - m_{...,l-1,...}
```
"""
LaplacianMap(T::Type, n::Int...) = LinearMap{T}(x -> laplacian(x, n...)[:], prod(n), issymmetric=true)
LaplacianMap(n::Int...) = LaplacianMap(Float64, n...)

"""
`DiffZMap([T,] n...) -> δz`

Construct a discretized z-derivative operator `δz` which will act on an `AbstractVector`.

Parameters
----------
* `T` : `Type`, optional

`Type` of argument of `δz`. Defaults to `Float64`.

* `n...`

Sequence of spatial dimensions of `δz`.

Usage
-----
* Forward map and Adjoint maps

The forward map `δz` multiplies a model vector of size `nz × nx × ny × ...` its first order z-derivative. The adjoint map is minus the forward map.

Description
-----------

The forward map computes the following operation
```math
δ_z m_{i,j,k,...} = (m_{l+1,...} - m_{l-1,...})/2
```
and the adjoint map computes `-δz`.
"""
DiffZMap(T::Type, n::Int...) = LinearMap{T}(x -> copy(deriv_z(x, n...)[:]),
                                            x -> copy(-deriv_z(x, n...)[:]),
                                            prod(n), prod(n))
DiffZMap(n::Int...) = DiffZMap(Float64, n...)

"""
`DiffXMap([T,] n...) -> δx`

Construct a discretized x-derivative operator `δx` which will act on an `AbstractVector`.

Parameters
----------
* `T` : `Type`, optional

`Type` of argument of `δx`. Defaults to `Float64`.

* `n...`

Sequence of spatial dimensions of `δx`.

Usage
-----
* Forward map and Adjoint maps

The forward map `δx` multiplies a model vector of size `nz × nx × ny × ...` its first order x-derivative. The adjoint map is minus the forward map.

Description
-----------

The forward map computes the following operation
```math
δ_x m_{i,j,k,...} = (m_{l+1,...} - m_{l-1,...})/2
```
and the adjoint map computes `-δx`.
"""
DiffXMap(T::Type, n::Int...) = LinearMap{T}(x ->  deriv_x(x, n...)[:],
                                            x -> -deriv_x(x, n...)[:], prod(n), prod(n))
DiffXMap(n::Int...) = DiffXMap(Float64, n...)

"""
`GradDivMap([T,] n...) -> GD`

Construct a discretized gradient operator GD which will act on an `AbstractVector`.

Parameters
----------
* `T` : `Type`, optional

`Type` of argument of `GD`. Defaults to `Float64`.

* `n...`

Sequence of spatial dimensions of `GD`.

Usage
-----
* Forward map

Calculates the discrete gradient of a `nz × nx × ny × ...` using first order central differences.
```math
∇m_{i,j,k,...} = [δ_x m_{i,j,k,...},  ..., δ_z m_{i,j,k,...}]
```

* Adjoint map

Calculates the discrete negative divergence of a `nz × nx × ny × ...` using first order central differences.
```math
∇\\cdot m_{i,j,k,...} = δ_x m_{i,j,k,...} +  ... + δ_z m_{i,j,k,...}
```
"""
GradDivMap(T::Type, nz, nx) = LinearMap{T}(x -> gradient(x, nz, nx)[:],
                                           x -> -divergence(x, nz, nx)[:], 2nz*nx, nz*nx)
GradDivMap(T::Type, nz, nx, ny) = LinearMap{T}(x -> gradient(x, nz, nx, ny)[:],
                                               x -> -divergence(x, nz, nx, ny)[:], 3nz*nx*ny, nz*nx*ny)
GradDivMap(n...) = GradDivMap(Float64, n...)

# Auxiliary functions
function wav_conv(data::AbstractArray, wavelet::AbstractVector, nr, nt)
    kt = argmax(abs.(hilbert(wavelet)))
    return mapslices(x->conv(x, wavelet)[kt:kt+nt-1],
                     reshape(copy(data), nr, nt), dims=[2])
end

function wav_corr(data::AbstractArray, wavelet::AbstractVector, nr, nt)
    kt = argmax(abs.(hilbert(wavelet)))
    return mapslices(x->conv(x, wavelet[end:-1:1])[nt-kt+1:2nt-kt],
                     reshape(copy(data), nr, nt), dims=[2]);
end

function wav_conv(data::AbstractArray, wavelet::AbstractVector, nr, ns, nt)
    kt = argmax(abs.(hilbert(wavelet)))
    return mapslices(x->conv(x, wavelet)[kt:kt+nt-1],
                     reshape(copy(data), nr, ns, nt), dims=[3])
end

function wav_corr(data::AbstractArray, wavelet::AbstractVector, nr, ns, nt)
    kt = argmax(abs.(hilbert(wavelet)))
    return mapslices(x->conv(x, wavelet[end:-1:1])[nt-kt+1:2nt-kt],
                     reshape(copy(data), nr, ns, nt), dims=[3]);
end

laplacian(x::AbstractArray) = laplacian(x, size(x)...)
function laplacian(x::AbstractArray{T}, nz, nx) where T
    x_ = zeros(T, nz+2, nx+2)
    x_[2:nz+1, 2:nx+1] = reshape(x, nz, nx)

    dx = zero(x_)
    dz = zero(x_)
    @fastmath @inbounds @simd for i in 2:nz+1
        dz[i,:] = -x_[i-1,:] + 2x_[i,:] + -x_[i+1,:]
    end
    @fastmath @inbounds @simd for i in 2:nx+1
        dx[:,i] = -x_[:,i-1] + 2x_[:,i] + -x_[:,i+1]
    end
    return dx[2:nz+1, 2:nx+1] + dz[2:nz+1, 2:nx+1]
end

function laplacian(x::AbstractArray{T}, nz, nx, ny) where T
    x_ = zeros(T, nz+2, nx+2, ny+2)
    x_[2:nz+1, 2:nx+1, 2:ny+1] = reshape(x, nz, nx, ny)

    dx = zero(x_)
    dz = zero(x_)
    dy = zero(x_)
    @fastmath @inbounds @simd for i in 2:nz+1
        dz[i,:,:] = -x_[i-1,:,:] + 2x_[i,:,:] + -x_[i+1,:,:]
    end
    @fastmath @inbounds @simd for i in 2:nx+1
        dx[:,i,:] = -x_[:,i-1,:] + 2x_[:,i,:] + -x_[:,i+1,:]
    end
    @fastmath @inbounds @simd for i in 2:ny+1
        dy[:,:,i] = -x_[:,:,i-1] + 2x_[:,:,i] + -x_[:,:,i+1]
    end
    return dx[2:nz+1, 2:nx+1, 2:ny+1] + dy[2:nz+1, 2:nx+1, 2:ny+1] + dz[2:nz+1, 2:nx+1, 2:ny+1]
end

deriv_z(x::AbstractArray) = deriv_z(x, size(x)...)
function deriv_z(x::AbstractArray{T}, n...) where T
    x_ = zeros(T, n[1]+2, n[2:end]...)

    colons = [Colon() for i=2:length(n)]

    x_[2:n[1]+1, colons...] = reshape(x, n[1], n[2:end]...)

    dz = zero(x_)
    @fastmath @inbounds @simd for i in 2:n[1]+1
        dz[i,colons...] = (x_[i+1,colons...] - x_[i-1,colons...])/2.
    end
    return dz[2:n[1]+1, colons...]
end

deriv_x(x::AbstractArray) = deriv_x(x, size(x)...)
function deriv_x(x::AbstractArray{T}, n...) where T
    x_ = zeros(T, n[1], n[2]+2, n[3:end]...)

    colons = [Colon() for i=3:length(n)]

    x_[:, 2:n[2]+1, colons...] = reshape(x, n[1], n[2], n[3:end]...)

    dx = zero(x_)
    @fastmath @inbounds @simd for i in 2:n[2]+1
        dx[:,i,colons...] = (x_[:,i+1,colons...] - x_[:,i-1,colons...])/2.
    end
    return dx[:,2:n[2]+1, colons...]
end

deriv_y(x::AbstractArray) = deriv_y(x, size(x)...)
function deriv_y(x::AbstractArray{T}, n...) where T
    x_ = zeros(T, n[1:2]..., n[3]+2, n[4:end]...)

    colons = [Colon() for i=4:length(n)]

    x_[:, :, 2:n[3]+1, colons...] = reshape(x, n[1:2]..., n[3], n[4:end]...)

    dy = zero(x_)
    @fastmath @inbounds @simd for i in 2:n[3]+1
        dy[:,:,i,colons...] = (x_[:,:,i+1,colons...] - x_[:,:,i-1,colons...])/2.
    end
    return dy[:,:,2:n[3]+1, colons...]
end

gradient(x::AbstractArray) = gradient(x, size(x)...)
function gradient(x::AbstractArray{T}, nz, nx) where T
    grad = zeros(T, nz, nx, 2)
    grad[:,:,1] = deriv_x(x, nz, nx)
    grad[:,:,2] = deriv_z(x, nz, nx)
    return grad
end

function gradient(x::AbstractArray{T}, nz, nx, ny) where T
    grad = zeros(T, nz, nx, ny, 3)
    grad[:,:,:,1] = deriv_x(x, nz, nx, ny)
    grad[:,:,:,2] = deriv_y(x, nz, nx, ny)
    grad[:,:,:,3] = deriv_z(x, nz, nx, ny)
    return grad
end

divergence(x::AbstractArray) = divergence(x, size(x)...)
function divergence(x::AbstractArray, nz, nx)
    x_ = reshape(x, nz, nx, 2)
    return deriv_x(x_[:,:,1], nz, nx) + deriv_z(x_[:,:,2], nz, nx)
end

function divergence(x::AbstractArray, nz, nx, ny)
    x_ = reshape(x, nz, nx, ny, 3)
    return deriv_x(x_[:,:,:,1], nz, nx, ny) + deriv_y(x_[:,:,:,2], nz, nx, ny) + deriv_z(x_[:,:,:,3], nz, nx, ny)
end


text:
"""
Test that `DiscreteIterator` has the same semantics as
`DifferentialEquations.solve`.
"""
module TestDiscreteIterator

using Test
using OrdinaryDiffEq: DiscreteProblem
using LyapunovExponents: DiscreteIterator, keepgoing!, change_tspan,
    solve_discrete

one_step() = DiscreteProblem(
    (u_next, u, p, t) -> u_next .= 2 .* u,
    [1.0],
    (0, 1),
)

two_steps() = change_tspan(one_step(), 2)
ten_steps() = change_tspan(one_step(), 10)

@testset "$make_prob" for make_prob in [
            one_step,
            two_steps,
            ten_steps,
        ]
    prob = make_prob()

    iter = DiscreteIterator(prob)
    keepgoing!(iter)

    sol = solve_discrete(prob)

    @test iter.u0 ≈ sol[end]
end

end


text:
"""
    latexify(t::RootedTree)

Return a LaTeX representation of the rooted tree `t`. This makes use of the
LaTeX package [forest](https://ctan.org/pkg/forest) and assumes that you use
the following LaTeX code in the preamble.

```
% Butcher trees, cf. https://tex.stackexchange.com/questions/283343/butcher-trees-in-tikz
\\usepackage{forest}
\\forestset{
  */.style={
    delay+={append={[]},}
  },
  rooted tree/.style={
    for tree={
      grow'=90,
      parent anchor=center,
      child anchor=center,
      s sep=2.5pt,
      if level=0{
        baseline
      }{},
      delay={
        if content={*}{
          content=,
          append={[]}
        }{}
      }
    },
    before typesetting nodes={
      for tree={
        circle,
        fill,
        minimum width=3pt,
        inner sep=0pt,
        child anchor=center,
      },
    },
    before computing xy={
      for tree={
        l=5pt,
      }
    }
  }
}
\\DeclareDocumentCommand\\rootedtree{o}{\\Forest{rooted tree [#1]}}
```

# Examples

```jldoctest
julia> rootedtree([1, 2, 2]) |> RootedTrees.latexify |> println
\\rootedtree[[][]]

julia> rootedtree([1, 2, 3, 3, 2]) |> RootedTrees.latexify |> println
\\rootedtree[[[][]][]]
```
"""
function latexify(t::RootedTree)
  if isempty(t)
    return "\\varnothing"
  end
  list_representation = butcher_representation(t, false)
  "\\rootedtree" * replace(list_representation, "τ" => "[]")
end

Latexify.@latexrecipe function _(t::RootedTree)
  return Latexify.LaTeXString(RootedTrees.latexify(t))
end