text:
module Pagano_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshUtilModule
using FinEtools.AlgoBaseModule
using Statistics: mean

function Pagano_3lay_cyl_bend_MST10_conv()
    elementtag = "MST10"
    println("""
    Pagano 3layer cylindrical bending: $(elementtag)
    """)

    # This example provides three-dimensional finite element model for  the
    # transverse shear stress calculations. The problem consists of a one-, two- or
    # three-layer plate subjected to a sinusoidal distributed load, as
    # described by Pagano (1969). The resulting transverse shear and axial
    # stresses through the thickness of the plate are compared to two existing
    # analytical solutions by Pagano (1969). The first solution is derived from
    # classical laminated plate theory (CPT), while the second is an exact
    # solution from linear elasticity theory.

    for (extrap, nodevalmeth) = zip([:extrapmean, :extraptrend, :default], [:averaging, :averaging, :invdistance])
        for Span_to_thickness = [4, 50]

            filebase = "Pagano_3lay_cyl_bend_conv_$(elementtag)_$(Span_to_thickness)_$(extrap)"
            modeldatasequence = FDataDict[]
            for Refinement = [1, 2, 4, 8]

                # Orthotropic material for the 3 layers
                E1 = 25e6*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;
                T = 2.5*phun("in"); # total thickness of the plate
                L = Span_to_thickness*T;
                h = 1.0/Refinement*phun("in");  # depth of the plate
                q0 = 1*phun("PSI")
                CTE1 =  CTE2 =  CTE3 = 0.0

                # Here we define the layout and the thicknesses of the layers.
                angles = vec([0.0 90.0 0.0]);
                nLayers = length(angles)
                ts = T/nLayers * ones(nLayers); # layer thicknesses

                tolerance = 0.0001*T

                # Select how fine the mesh should be
                nL, nh = Refinement * 2 * 6, 1;
                nts= Refinement * 2 * ones(Int, nLayers);# number of elements per layer

                xs = collect(linearspace(0.0, L, nL+1))
                ys = collect(linearspace(0.0, h, nh+1))
                fens,fes = T10layeredplatex(xs, ys, ts, nts)
                display(typeof(fes))
                println("Mesh: nL, nh, nts = $nL, $nh, $nts")
                println("count(fens) = $(count(fens))")

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

                # 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; 0.0; 1.0]);
                end

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

                # We will create 3 regions, one for each of the layers
                regions = FDataDict[]
                for layer = 1:nLayers
                    rls = selectelem(fens, fes, label =  layer)
                    push!(regions, FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial)))
                end

                # The essential boundary conditions are applied to enforce the plane strain constraint.
                ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
                lyh = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
                ey = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly0, lyh))
                # The transverse displacement is fixed at the two ends.
                lz0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
                lzL = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
                ez = FDataDict("displacement"=>  0.0, "component"=> 3, "node_list"=>vcat(lz0, lzL))
                ex = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>[1])

                # The traction boundary condition is applied at the top of the plate.
                bfes = meshboundary(fes)
                function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                    forceout[1] = 0.0
                    forceout[2] = 0.0
                    forceout[3] = -q0*sin(pi*XYZ[1]/L)
                    return forceout
                end
                # From  the entire boundary we select those quadrilaterals that lie on the plane
                # Z = thickness
                tl = selectelem(fens, bfes, box = [-Inf Inf -Inf Inf T T], inflate=tolerance)
                Trac = FDataDict("traction_vector"=>pfun,
                "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), SimplexRule(2, 3))))

                modeldata = FDataDict("fens"=>fens,
                "regions"=>regions,
                "essential_bcs"=>[ex, ey, ez],
                "traction_bcs"=> [Trac]
                )
                modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

                # modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
                # modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

                u = modeldata["u"]
                geom = modeldata["geom"]

                # The results of the displacement and stresses will be reported at
                # nodes located at the appropriate points.
                ntopcenter = selectnode(fens, box=[L/2 L/2 0.0 h T T], inflate=tolerance)
                ncenterline = selectnode(fens, box=[L/2 L/2 0.0 0.0 0.0 T], inflate=tolerance)
                nx0line = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 T], inflate=tolerance)

                zclo = sortperm(vec(geom.values[ncenterline, 3]))
                ncenterline = ncenterline[zclo]
                centerz = geom.values[ncenterline, 3]

                println("Top Center deflection: $(mean(u.values[ntopcenter, 3], dims = 1)/phun("in")) [in]")

                #  Compute  all stresses
                modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s", "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3), "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
                modeldata = AlgoDeforLinearModule.exportstress(modeldata)

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

            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 # Span_to_thickness
    end # extrap
    println("Done")

end # Pagano_3lay_cyl_bend_MST10_conv

function Pagano_3lay_cyl_bend_H8_conv()
    elementtag = "H8"
    println("""
    Pagano 3layer cylindrical bending: $(elementtag)
    """)

    # This example provides three-dimensional finite element model for  the
    # transverse shear stress calculations. The problem consists of a one-, two- or
    # three-layer plate subjected to a sinusoidal distributed load, as
    # described by Pagano (1969). The resulting transverse shear and axial
    # stresses through the thickness of the plate are compared to two existing
    # analytical solutions by Pagano (1969). The first solution is derived from
    # classical laminated plate theory (CPT), while the second is an exact
    # solution from linear elasticity theory.

    for (extrap, nodevalmeth) = zip([:extrapmean, :extraptrend, :default], [:averaging, :averaging, :invdistance])
        for Span_to_thickness = [4, 50]

            filebase = "Pagano_3lay_cyl_bend_conv_$(elementtag)_$(Span_to_thickness)_$(extrap)"
            modeldatasequence = FDataDict[]
            for Refinement = [1, 2, 4, 8]

                # Orthotropic material for the 3 layers
                E1 = 25e6*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;
                T = 2.5*phun("in"); # total thickness of the plate
                L = Span_to_thickness*T;
                h = 1.0/Refinement*phun("in");  # depth of the plate
                q0 = 1*phun("PSI")
                CTE1 =  CTE2 =  CTE3 = 0.0

                # Here we define the layout and the thicknesses of the layers.
                angles = vec([0.0 90.0 0.0]);
                nLayers = length(angles)
                ts = T/nLayers * ones(nLayers); # layer thicknesses

                tolerance = 0.0001*T

                # Select how fine the mesh should be
                nL, nh = Refinement * 2 * 6, 1;
                nts= Refinement * 2 * ones(Int, nLayers);# number of elements per layer

                xs = collect(linearspace(0.0, L, nL+1))
                ys = collect(linearspace(0.0, h, nh+1))
                fens,fes = H8layeredplatex(xs, ys, ts, nts)

                println("Mesh: nL, nh, nts = $nL, $nh, $nts")
                println("count(fens) = $(count(fens))")

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

                # 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; 0.0; 1.0]);
                end

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

                # We will create 3 regions, one for each of the layers
                regions = FDataDict[]
                for layer = 1:nLayers
                    rls = selectelem(fens, fes, label =  layer)
                    push!(regions, FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial)))
                end

                # The essential boundary conditions are applied to enforce the plane strain constraint.
                ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
                lyh = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
                ey = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly0, lyh))
                # The transverse displacement is fixed at the two ends.
                lz0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
                lzL = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
                ez = FDataDict("displacement"=>  0.0, "component"=> 3, "node_list"=>vcat(lz0, lzL))
                ex = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>[1])

                # The traction boundary condition is applied at the top of the plate.
                bfes = meshboundary(fes)
                function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                    forceout[1] = 0.0
                    forceout[2] = 0.0
                    forceout[3] = -q0*sin(pi*XYZ[1]/L)
                    return forceout
                end
                # From  the entire boundary we select those quadrilaterals that lie on the plane
                # Z = thickness
                tl = selectelem(fens, bfes, box = [-Inf Inf -Inf Inf T T], inflate=tolerance)
                Trac = FDataDict("traction_vector"=>pfun,                "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), GaussRule(2, 2))))

                modeldata = FDataDict("fens"=>fens, "regions"=>regions, "essential_bcs"=>[ex, ey, ez], "traction_bcs"=> [Trac])
                modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

                # modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
                # modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

                u = modeldata["u"]
                geom = modeldata["geom"]

                # The results of the displacement and stresses will be reported at
                # nodes located at the appropriate points.
                ntopcenter = selectnode(fens, box=[L/2 L/2 0.0 h T T], inflate=tolerance)
                ncenterline = selectnode(fens, box=[L/2 L/2 0.0 0.0 0.0 T], inflate=tolerance)
                nx0line = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 T], inflate=tolerance)

                zclo = sortperm(vec(geom.values[ncenterline, 3]))
                ncenterline = ncenterline[zclo]
                centerz = geom.values[ncenterline, 3]

                println("Top Center deflection: $(mean(u.values[ntopcenter, 3], dims = 1)/phun("in")) [in]")

                #  Compute  all stresses
                modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s",
                "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3),
                "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
                modeldata = AlgoDeforLinearModule.exportstress(modeldata)

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

            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 # Span_to_thickness
    end # extrap
    println("Done")

end # Pagano_3lay_cyl_bend_H8_conv

function Pagano_3lay_cyl_bend_MSH8_conv()
    elementtag = "MSH8"
    println("""
    Pagano 3layer cylindrical bending: $(elementtag)
    """)

    # This example provides three-dimensional finite element model for  the
    # transverse shear stress calculations. The problem consists of a one-, two- or
    # three-layer plate subjected to a sinusoidal distributed load, as
    # described by Pagano (1969). The resulting transverse shear and axial
    # stresses through the thickness of the plate are compared to two existing
    # analytical solutions by Pagano (1969). The first solution is derived from
    # classical laminated plate theory (CPT), while the second is an exact
    # solution from linear elasticity theory.

    for (extrap, nodevalmeth) = zip([:extrapmean, :extraptrend, :default], [:averaging, :averaging, :invdistance])
        for Span_to_thickness = [4, 50]

            filebase = "Pagano_3lay_cyl_bend_conv_$(elementtag)_$(Span_to_thickness)_$(extrap)"
            modeldatasequence = FDataDict[]
            for Refinement = [1, 2, 4, 8]

                # Orthotropic material for the 3 layers
                E1 = 25e6*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;
                T = 2.5*phun("in"); # total thickness of the plate
                L = Span_to_thickness*T;
                h = 1.0/Refinement*phun("in");  # depth of the plate
                q0 = 1*phun("PSI")
                CTE1 =  CTE2 =  CTE3 = 0.0

                # Here we define the layout and the thicknesses of the layers.
                angles = vec([0.0 90.0 0.0]);
                nLayers = length(angles)
                ts = T/nLayers * ones(nLayers); # layer thicknesses

                tolerance = 0.0001*T

                # Select how fine the mesh should be
                nL, nh = Refinement * 2 * 6, 1;
                nts= Refinement * 2 * ones(Int, nLayers);# number of elements per layer

                xs = collect(linearspace(0.0, L, nL+1))
                ys = collect(linearspace(0.0, h, nh+1))
                fens,fes = H8layeredplatex(xs, ys, ts, nts)

                println("Mesh: nL, nh, nts = $nL, $nh, $nts")
                println("count(fens) = $(count(fens))")

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

                # 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; 0.0; 1.0]);
                end

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

                # We will create 3 regions, one for each of the layers
                regions = FDataDict[]
                for layer = 1:nLayers
                    rls = selectelem(fens, fes, label =  layer)
                    push!(regions, FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial)))
                end

                # The essential boundary conditions are applied to enforce the plane strain constraint.
                ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
                lyh = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
                ey = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly0, lyh))
                # The transverse displacement is fixed at the two ends.
                lz0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
                lzL = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
                ez = FDataDict("displacement"=>  0.0, "component"=> 3, "node_list"=>vcat(lz0, lzL))
                ex = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>[1])

                # The traction boundary condition is applied at the top of the plate.
                bfes = meshboundary(fes)
                function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                    forceout[1] = 0.0
                    forceout[2] = 0.0
                    forceout[3] = -q0*sin(pi*XYZ[1]/L)
                    return forceout
                end
                # From  the entire boundary we select those quadrilaterals that lie on the plane
                # Z = thickness
                tl = selectelem(fens, bfes, box = [-Inf Inf -Inf Inf T T], inflate=tolerance)
                Trac = FDataDict("traction_vector"=>pfun, "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), GaussRule(2, 2))))

                modeldata = FDataDict("fens"=>fens, "regions"=>regions, "essential_bcs"=>[ex, ey, ez], "traction_bcs"=> [Trac])
                modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

                # modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
                # modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

                u = modeldata["u"]
                geom = modeldata["geom"]

                # The results of the displacement and stresses will be reported at
                # nodes located at the appropriate points.
                ntopcenter = selectnode(fens, box=[L/2 L/2 0.0 h T T], inflate=tolerance)
                ncenterline = selectnode(fens, box=[L/2 L/2 0.0 0.0 0.0 T], inflate=tolerance)
                nx0line = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 T], inflate=tolerance)

                zclo = sortperm(vec(geom.values[ncenterline, 3]))
                ncenterline = ncenterline[zclo]
                centerz = geom.values[ncenterline, 3]

                println("Top Center deflection: $(mean(u.values[ntopcenter, 3], dims = 1)/phun("in")) [in]")

                #  Compute  all stresses
                modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s", "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3), "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
                modeldata = AlgoDeforLinearModule.exportstress(modeldata)

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

            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 # Span_to_thickness
    end # extrap
    println("Done")

end # Pagano_3lay_cyl_bend_MSH8_conv

function Pagano_3lay_cyl_bend_T10_conv()
    elementtag = "T10"
    println("""
    Pagano 3layer cylindrical bending: $(elementtag)
    """)

    # This example provides three-dimensional finite element model for  the
    # transverse shear stress calculations. The problem consists of a one-, two- or
    # three-layer plate subjected to a sinusoidal distributed load, as
    # described by Pagano (1969). The resulting transverse shear and axial
    # stresses through the thickness of the plate are compared to two existing
    # analytical solutions by Pagano (1969). The first solution is derived from
    # classical laminated plate theory (CPT), while the second is an exact
    # solution from linear elasticity theory.

    for (extrap, nodevalmeth) = zip([:default], [:invdistance])
        for Span_to_thickness = [4, 50]

            filebase = "Pagano_3lay_cyl_bend_conv_$(elementtag)_$(Span_to_thickness)_$(extrap)"
            modeldatasequence = FDataDict[]
            for Refinement = [1, 2, 4, 8]

                # Orthotropic material for the 3 layers
                E1 = 25e6*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;
                T = 2.5*phun("in"); # total thickness of the plate
                L = Span_to_thickness*T;
                h = 1.0/Refinement*phun("in");  # depth of the plate
                q0 = 1*phun("PSI")
                CTE1 =  CTE2 =  CTE3 = 0.0

                # Here we define the layout and the thicknesses of the layers.
                angles = vec([0.0 90.0 0.0]);
                nLayers = length(angles)
                ts = T/nLayers * ones(nLayers); # layer thicknesses

                tolerance = 0.0001*T

                # Select how fine the mesh should be
                nL, nh = Refinement * 2 * 6, 1;
                nts= Refinement * 2 * ones(Int, nLayers);# number of elements per layer

                xs = collect(linearspace(0.0, L, nL+1))
                ys = collect(linearspace(0.0, h, nh+1))
                fens,fes = T10layeredplatex(xs, ys, ts, nts)

                println("Mesh: nL, nh, nts = $nL, $nh, $nts")
                println("count(fens) = $(count(fens))")

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

                # 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; 0.0; 1.0]);
                end

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

                # We will create 3 regions, one for each of the layers
                regions = FDataDict[]
                for layer = 1:nLayers
                    rls = selectelem(fens, fes, label =  layer)
                    push!(regions, FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial)))
                end

                # The essential boundary conditions are applied to enforce the plane strain constraint.
                ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
                lyh = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
                ey = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly0, lyh))
                # The transverse displacement is fixed at the two ends.
                lz0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
                lzL = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
                ez = FDataDict("displacement"=>  0.0, "component"=> 3, "node_list"=>vcat(lz0, lzL))
                ex = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>[1])

                # The traction boundary condition is applied at the top of the plate.
                bfes = meshboundary(fes)
                function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                    forceout[1] = 0.0
                    forceout[2] = 0.0
                    forceout[3] = -q0*sin(pi*XYZ[1]/L)
                    return forceout
                end
                # From  the entire boundary we select those quadrilaterals that lie on the plane
                # Z = thickness
                tl = selectelem(fens, bfes, box = [-Inf Inf -Inf Inf T T], inflate=tolerance)
                Trac = FDataDict("traction_vector"=>pfun,
                "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), SimplexRule(2, 3))))

                modeldata = FDataDict("fens"=>fens,
                "regions"=>regions,
                "essential_bcs"=>[ex, ey, ez],
                "traction_bcs"=> [Trac]
                )
                modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

                # modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
                # modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

                u = modeldata["u"]
                geom = modeldata["geom"]

                # The results of the displacement and stresses will be reported at
                # nodes located at the appropriate points.
                ntopcenter = selectnode(fens, box=[L/2 L/2 0.0 h T T], inflate=tolerance)
                ncenterline = selectnode(fens, box=[L/2 L/2 0.0 0.0 0.0 T], inflate=tolerance)
                nx0line = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 T], inflate=tolerance)

                zclo = sortperm(vec(geom.values[ncenterline, 3]))
                ncenterline = ncenterline[zclo]
                centerz = geom.values[ncenterline, 3]

                println("Top Center deflection: $(mean(u.values[ntopcenter, 3], dims = 1)/phun("in")) [in]")

                #  Compute  all stresses
                modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s",  "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3), "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
                modeldata = AlgoDeforLinearModule.exportstress(modeldata)

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

            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 # Span_to_thickness
    end # extrap
    println("Done")

end # Pagano_3lay_cyl_bend_T10_conv

function allrun()
    println("#####################################################")
    println("# Pagano_3lay_cyl_bend_MST10_conv ")
    Pagano_3lay_cyl_bend_MST10_conv()
    println("#####################################################")
    println("# Pagano_3lay_cyl_bend_MSH8_conv ")
    Pagano_3lay_cyl_bend_MSH8_conv()
    println("#####################################################")
    println("# Pagano_3lay_cyl_bend_H8_conv ")
    Pagano_3lay_cyl_bend_H8_conv()
    println("#####################################################")
    println("# Pagano_3lay_cyl_bend_MSH8_conv ")
    Pagano_3lay_cyl_bend_MSH8_conv()
    println("#####################################################")
    println("# Pagano_3lay_cyl_bend_T10_conv ")
    Pagano_3lay_cyl_bend_T10_conv()
    return true
end # function allrun

end # module Pagano_examples


text:
__precompile__()
module WRDSFuns

using SASLib
using DataFrames
using Missings
using ShiftedArrays
using Compat.Dates

export compustat2crspjoin, portfolioassignments!
export dropcols,alltomissing!,tomissing!,disallowmissing!

include("crspcompustatfuns.jl")
include("utilities.jl")

end

ext: jl


text:
using Molly

using Base.Threads
using Statistics
using Test

@testset "Lennard-Jones gas energy conservation" begin
    temp = 1.
    timestep = 0.005
    n_steps = 10_000
    box_size = 50.0
    n_atoms = 2000
    m = 0.8 * box_size^3 / n_atoms

    parallel_list = nthreads() > 1 ? (false, true) : (false,)
    lj_potentials = (LennardJones{false}(ShiftedPotentialCutoff(3.), false),
                     LennardJones{true}(ShiftedPotentialCutoff(3.), false),
                     LennardJones{false}(ShiftedForceCutoff(3.), false),
                     LennardJones{true}(ShiftedForceCutoff(3.), false))

    for parallel in parallel_list
        @testset "$lj_potential" for lj_potential in lj_potentials
            s = Simulation(
                simulator=VelocityVerlet(),
                atoms=[Atom(attype="Ar", name="Ar", resnum=i, resname="Ar", charge=0.0,
                            mass=m, σ=0.3, ϵ=0.2) for i in 1:n_atoms],
                general_inters=(lj_potential,),
                coords=placeatoms(n_atoms, box_size, 0.6),
                velocities=[velocity(10.0, temp) .* 0.01 for i in 1:n_atoms],
                temperature=temp,
                box_size=box_size,
                loggers=Dict("coords" => CoordinateLogger(100),
                            "energy" => EnergyLogger(100)),
                timestep=timestep,
                n_steps=n_steps
            )

            E0 = energy(s)
            @time simulate!(s, parallel=parallel)

            ΔE = energy(s) - E0
            @test abs(ΔE) < 2e-2

            Es = s.loggers["energy"].energies
            maxΔE = maximum(abs.(Es .- E0))
            @test maxΔE < 2e-2

            @test abs(Es[end] - Es[1]) < 2e-2
            @test std(Es.-Es[1])/n_atoms < timestep^2

            final_coords = last(s.loggers["coords"].coords)
            @test minimum(minimum.(final_coords)) > 0.0
            @test maximum(maximum.(final_coords)) < box_size
        end
    end
end


text:
function SIR_simple__sir__assign__I__0(I, infected, recovered)
    return I + infected - recovered
end

function SIR_simple__sir__assign__S__0(S, infected)
    return S - infected
end

function SIR_simple__sir__assign__R__0(R, recovered)
    return R + recovered
end

function SIR_simple__sir__assign__infected__0(beta, S, I, R, dt)
    return (-(beta*S*I) / (S + I + R)) * dt
end

function SIR_simple__sir__assign__recovered__0(gamma, I, dt)
    return (gamma*I) * dt
end


text:
using Plots, Utilities, ForwardDiff, AverageShiftedHistograms,JLD
plotly()

Model = fooSdErgmDirBin1
T = 300
N=10
deltaN = 3
GBA = 1
GW = N # N#N # 20
UMdgp,dgpdegs = linSpacedPar(Model,N; NgroupsW = GW,deltaN = deltaN)
groupsIndAB_I = Utilities.distributeAinVecN(Array{Int64}(1:GBA),GW)
groupsIndAB_O = groupsIndAB_I
splitVec(Vector(1:10))
groupsInds = [Utilities.distributeAinVecN(Array{Int64}(1:GW),N), groupsIndAB_I, groupsIndAB_O ]
Bdgp = 0.9*ones(GBA)
Adgp = 0.1*ones(GBA)
Wdgp = UMdgp.*(1-Bdgp)
scoreScal ="FISHER-DIAG"# ""#"FISHER-EWMA"#
# Test sample
Mod = SdErgmDirBin1(ones(2,2),[Wdgp , Bdgp, Adgp],groupsInds,scoreScal)
Mod,Y_T,Fitness_T = sampl(Mod,T)

## Load Data
save_fold =   "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
file_nameStart = "Weekly_eMid_Estimates"
file_nameEnd = "_from_$(Dates.format(startDate,"yyyy_mm_dd"))"*
        "_to_" *"$(Dates.format(endDate,"yyyy_mm_dd"))"*".jld"
load_path = save_fold*file_nameStart* "_Bin1_" *file_nameEnd#
@load(save_path,estSS,estGas,estGasTarg,estSS_1GW,estGas_1GW,AeMidWeekly_T, YeMidWeekly_T ,weekInd,degsIO_T)

N = round(Int,length(degsIO_T[:,1])/2)
eMidMod = DynNets.SdErgmDirBin1(degsIO_T,estGas_1GW[1], [Int.(ones(2N)),Int.(ones(1)),Int.(ones(1))],"")
score_driven_filter_or_dgp(eMidMod,[estGas_1GW[1][1];estGas_1GW[1][2];estGas_1GW[1][3]])
estSS =  DynNets.estimateSnapSeq(eMidMod1Snap)'
## test fisher scoring on data
#time = @elapsed  out = estimate(eMidMod)

println((T,N,GW ,round(time),out[1][2],out[1][3]))
Fitness_T_estTarg, tmp = score_driven_filter_or_dgp( Mod, [out[1][1];out[1][2];out[1][3]])
out
#test the pox of autodiff loglikelihood
function loglikeRE(BAvecRe::Vector)
    B = BAvecRe[1]
    A = BAvecRe[2]
    ~, ll = score_driven_filter_or_dgp( eMidMod, [SSMean;B;A])
    return ll
end
loglikeUN(BAvecUn::Vector) = ( ReB = exp(BAvecUn[1])./(1+exp(BAvecUn[1])); ReA = exp(BAvecUn[2]) ; loglikeRE([ReB;ReA]) )

score = x -> ForwardDiff.gradient(loglikeUN,x)
hess = x -> ForwardDiff.hessian(loglikeUN,x)
x0Re = [0.05;0.5]
    α = 1
    x0Un = [log(x0Re[1]/(1-x0Re[1]));log(x0Re[2])]
    x = x0Un
    llVal = loglikeUN(x0Un)
    maxIt = 1000;it= 0
    minScoreNorm = 10*eps(); scoreNorm = 10
    minObjImp = 10*eps();objImp = 10
    @elapsed while (it<maxIt)& (scoreNorm>minScoreNorm)&(objImp>minObjImp)
        it+=1
        score_x = score(x)
        scoreNorm = sum(score_x.^2)
        x = x - α* hess(x)\score_x
        llValNew = loglikeUN(x)
        objImp = abs(llVal - llValNew)
        llVal = llValNew
        println((it,exp(x[1])/(1+exp(x[1])),exp(x[2]),llValNew ,objImp,scoreNorm))
    end

plot(score_driven_filter_or_dgp( eMidMod, [SSMean;x0Re])[1][:,1+N:end])
##
Fitness_T = Float64.(Fitness_T)
corrMat=  squeeze(StatsBase.crosscor( Fitness_T,Fitness_T,[0]),1)
err = corrMat[isfinite(corrMat)]
tmpAsh = ash(err)
plot(tmpAsh)

## Test Target Estimate

estOutTarg = estimateTargeting(Mod)
Mod_estTarg = SdErgm1(Mod.obsT,estOutTarg[1],groupsInd,scoreScal)

vresParTarg = [Mod_estTarg.Par[1];Mod_estTarg.Par[2];Mod_estTarg.Par[3]]
Fitness_T_estTarg, tmp = score_driven_filter_or_dgp( Mod, vresParTarg)

ppar = plot(Fitness_T_estTarg,title = "$(round.(squeeze(mean(Fitness_T_estTarg,1),1),1))  $(round.(UMdgp,1))")
pdegs = plot(degs_T,title = "$(dgpdegs)  $(round.(squeeze(mean(degs_T,1),1),1))")
plot(ppar,pdegs,layout = (2,1),legend=:none,size = (1200,600))

## Test estimate

estOut = estimate(Mod)
Mod_est = SdErgm1(Mod.obsT,estOut[1],groupsInds,scoreScal)

vresPar = [Mod_est.Par[1];Mod_est.Par[2];Mod_est.Par[3]]
Fitness_T_est, tmp = score_driven_filter_or_dgp( Mod, vresPar)

ppar = plot(Fitness_T_est,title = "$(round.(squeeze(mean(Fitness_T_est,1),1),1))  $(round.(UMdgp,1))")
pdegs = plot(degs_T,title = "$(dgpdegs)  $(round.(squeeze(mean(degs_T,1),1),1))")
plot(ppar,pdegs,layout = (2,1),legend=:none,size = (1200,600))

#

##
using JLD2

N_est = 300
Tvals =  [1000]
N_ind = length(Tvals)[1]

n  =1
N = 10 ;
GBAdgp =1; GWdgp = N  ;
GBAest =1; GWest = GWdgp; useStartVal = false

useStartVal? startValFlag = "T" : startValFlag = "F";# false #false #
# save
save_fold = "./data/estimatesTest/"
file_name = "EstTest_$(N)_$(GWdgp)_$(GBAdgp)_$(Tvals)_$(N_est)_$(useStartVal)_Miss_$(GWest)_$(GBAest)_rescaling.jld"
save_path = save_fold*file_name#
estDataT = @load(save_path,groupsIndsEst,groupsIndsDgp,allObs,
        simPar_1,simPar_2,simPar_3,estPar_1,estPar_2,estPar_3,est_conv_flag,est_times)

vestPar = [estPar_1[1,n,:]; estPar_2[n];estPar_3[n]]
Fitness_Tfil ,like= score_driven_filter_or_dgp( Mod,vestPar)

    pparfil = plot(Fitness_Tfil)
    pparfil = plot(Fitness_Tfil)
    plot(ppar,pparfil,layout = (2,1),legend=:none,size = (1200,600))
##
vResGasPar = [Mod.Par[1];Mod.Par[2];Mod.Par[3];]
Fitness_Tfil ,like= score_driven_filter_or_dgp( Mod,vResGasPar)

    pparfil = plot(Fitness_Tfil)
    pparfil = plot(Fitness_Tfil)
    plot(ppar,pparfil,layout = (2,1),legend=:none,size = (1200,600))

## test

#

# using Plots, AverageShiftedHistograms
# plotly()
#
# T = 500
# N=10
# GBA = 1
# GW = 2N#20
#
# UMdgp = identify!(fooSdErgmDBin1, [linSpacedFitnesses(N,mode = "UNDIRECTED");linSpacedFitnesses(N;mode = "UNDIRECTED")] )
#
# Bdgp = 0.9*ones(GBA)
# Adgp = 0.2*ones(GBA)
# Wdgp = UMdgp.*(1-Bdgp)
# groupsInd = [distributeAinVecN(Array{Int64}(1:GW),2N), distributeAinVecN(Array{Int64}(1:GBA),GW)]
#
#
#
# ##
#
# Mod= SdErgmDBin1(ones(2,2),[Wdgp , Bdgp, Adgp],groupsInd)
# Mod,Y_T,Fitness_T = sampl(Mod,T)
# degs_T = Mod.obsT
#
# N = round(Int,length(degs_T[1,:])/2)
#     pparI = plot(Fitness_T[:,1:N])
#     pparO = plot(Fitness_T[:,N+1:2N])
#     pdegsI = plot(degs_T[:,1:N])
#     pdegsO = plot(degs_T[:,N+1:2N]);
#     plot(pparI,pdegsI,pparO,pdegsO,layout = (4,1),legend=:none,size = (1200,600))
#
#
# ##
# err =degs_T[:,1]
#
# rng_space = 0:.001:6
# tmpAsh = ash(err;rng = rng_space)
# plot(tmpAsh)
# ##
# vResGasPar = [Mod.Par[1];Mod.Par[2];Mod.Par[3];]
# Fitness_Tfil ,like= score_driven_filter_or_dgpAndLikeliood( Mod,vResGasPar)
#
#
#     pparIfil = plot(Fitness_Tfil[:,1:N])
#     pparOfil = plot(Fitness_Tfil[:,N+1:2N])
#     plot(pparI,pparIfil,pparO,pparOfil,layout = (4,1),legend=:none,size = (1200,600))
#
#
# ## Test estimate
# estimate(Mod)
#
#
# # Mod2 = SdErgmDBin1(degs_T[1:2,:])
# # Fitness_T_snap = estSnapSeq(Mod2;degsIOT = repmat(Array{Int,1}(1:10)',2,1))
# #
# # plot(Fitness_T_snap')
#
# sampleDynNetsDBin1(T,"CONST";groupsInds=groupsInd, sinMean = UMdgp)
# ##
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
# #


text:
# Autogenerated wrapper script for GLU_jll for armv7l-linux-musleabihf
export libGLU

using Libglvnd_jll
JLLWrappers.@generate_wrapper_header("GLU")
JLLWrappers.@declare_library_product(libGLU, "libGLU.so.1")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll)
    JLLWrappers.@init_library_product(
        libGLU,
        "lib/libGLU.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct KernelParams
    # Size of the array for each thread.
    array_size::Int
    sampletime::Float64
    mode::String
    # The number of times to run each kernel. Larger kernels need to run for fewer
    # iterations.
    inner_iterations::Int
    # Callback for measurements.
    measurements::Vector{Any}
    vector_sizes::Vector{Int}
    temporals::Vector{Bool}
    array_type::String
end

function runkernels(params, f = (sz) -> Vector{Float32}(undef, sz))
    # Allocate and instantiate the array
    sz = Threads.nthreads() * (2 ^ params.array_size)
    @show sz
    A = f(sz)
    @show sizeof(A)
    threadme(vector_write, A, Val{16}())

    # Setup the counter
    pipe = connect(PIPEPATH)
    println(pipe, "sampletime $(params.sampletime)")

    # Do the benchmarking!
    benchmark(A, pipe, params)
end

function benchmark(A::AbstractArray, pipe, params::KernelParams)
    # Set up a bunch of nested loops to run through all the variations of tests we want
    # to execute.
    fns = [
        vector_sum,
        vector_write,
        vector_increment
    ]

    for (f, sz, _nt) in Iterators.product(fns, params.vector_sizes, params.temporals)
        # Create a NamedTuple of parameters to forward to the function.
        nt = (
            vector = Val{sz}(),
            aligned = Val{true}(),
            nontemporal = Val{_nt}(),
        )

        run(f, A, pipe, params, nt)
    end

    # # Prepare for the random access benchmarks.
    # for vector_size in params.vector_sizes
    #     sz = params.array_size

    #     # Get the number of unique indices to cycle through.
    #     lfsr_size = div(2 ^ sz, vector_size)
    #     for _nt in params.temporals
    #         nt = (
    #             vector = Val{vector_size}(),
    #             lfsr = LFSR(lfsr_size),
    #             nontemporal = Val{_nt}(),
    #         )

    #         run(hop_sum, A, pipe, params, nt)
    #         run(hop_write, A, pipe, params, nt)
    #         run(hop_increment, A, pipe, params, nt)
    #     end
    # end
    return nothing
end

function run(f, A, pipe, params, nt)
    # warm up the function.
    threadme(f, A, nt...; prepare = true, iterations = 1)
    iters = params.inner_iterations

    # Run all the desired measurements
    first = true
    for m in params.measurements
        # Create a named tuple of the all the run parameters.
        run_params = (
            sz = div(sizeof(A), 10^9),
            mode = params.mode,
            threads = Threads.nthreads(),
            array = params.array_type,
        )

        run_params = merge(run_params, nt)
        final_params = make_params(f, run_params)

        # Always store data in the same spot I guess ...
        filepath = joinpath(DATADIR, "array_data.jls")

        # Transfer the measurements and parameters to the `server` process.
        transfer(pipe, m)
        paramtransfer(pipe, final_params)

        println(pipe, "filepath $filepath")
        println(pipe, "start")
        sleep(2)
        threadme(f, A, nt...; iterations = iters)
        sleep(2)
        println(pipe, "stop")
    end
    return nothing
end

####
#### Main
####

# Hack to get Julia threads to stop being moved by the linux scheduler
function assign_affinities(cpustart = 24)
    # Get Julia's PID
    pid = getpid()

    # Get the threads from Julia. The first thread belonds to Julia.
    # The next 8 threads (by defualt) belong to BLAS.
    # The rest of the threads are extra Julia threads.
    threads = split(chomp(read(pipeline(`ps -e -T`, `grep $pid`), String)), "\n")

    julia_thread_strings = vcat(threads[1], threads[10:end])
    @assert length(julia_thread_strings) == Threads.nthreads()

    # Parse out the thread PID
    _parse(x) = parse(Int, split(x)[2])
    julia_threads = _parse.(julia_thread_strings)

    # Start assigning affinity.
    for (index, pid) in enumerate(julia_threads)
        Base.run(`taskset -cp $(cpustart - 1 + index) $pid`)
    end
    return nothing
end

#####
##### Kernel for testing if the DRAM cache is inclusive.
#####

function cache_smasher(size_per_thread, f = (sz) -> Vector{Float32}(undef, sz))
    # Determine how big our array should be and instantiate it.
    sz = Threads.nthreads() * (2 ^ size_per_thread)
    A = f(sz)
    # Make sure the OS instantiates the whole array.
    threadme(vector_write, A, Val{16}())

    # Now, make sure everything in the cache is clean.
    threadme(vector_sum, A, Val{16}())
    return _cache_smasher(A, f)
end

# Allow an entry point for if the array is already allocated.
function _cache_smasher(A, f)
    # Allocate another array that figs in the L3 cache (33 MB)
    #
    # Here, we make an array that is 4 MB.
    B = f(2^17)
    @assert eltype(B) == Float32

    # Pass in another vector to record timings.
    timings = Float64[]
    threadme(_smasher, A, B, timings)
    return timings
end

function _smasher(x, small, timings)
    # We want one thread to iterate over the small array.
    #
    # All other threads blow out the cache
    if Threads.threadid() == 1
        _cycle(small, timings)
    else
        # Sleep so the cycler can get some inital data
        sleep(5)
        for _ in 1:2
            vector_sum(x, Val{16}())
        end
    end
end

function _cycle(x, timings)
    # Time for 20 seconds.
    start = now()
    while now() < start + Second(60)
        t = @elapsed for _ in 1:1000
            vector_sum(x, Val{16}())
        end
        push!(timings, t)
    end
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using MesocosmMetabolism, CSV, Turing, DataFrames, MCMCChains, StatsPlots, Interpolations, ReverseDiff
times = collect(4.0:0.1:75.0)
nDays = 3
dayStart = 4
days = [convert(Int64, floor(times[i]/24.0 - (dayStart/24.0)) + 1) for i in 1:length(times)]
PARData= [if (7<= (times[i] - 24*(days[i]-1)) < 19) 1 else 0 end for i in 1:length(times)]
OSatData = repeat([8.2], length(times))
OSat = Interpolations.interpolate((times,), OSatData, Interpolations.Gridded(Interpolations.Linear()))
GPP_daily_mu = 1.72 # g/(L*hr) - only applies to 12 hours when lights are on
GPP_daily_sigma = 0.01
ER_daily_mu = -2.0 # g/(L*hr)
ER_daily_sigma = 0.01
k600_meanlog = -0.2 # hr^-1
k600_sdlog = 0.01
Q = 1.794
dayStart = 4
K02_conv = 1.0
err_obs_iid_sigma_scale = 0.015
err_proc_iid_sigma_scale = 0.01
u0 = [6.0]

problem = initialize_process_model(
    times = times,
    PAR = Interpolations.ConstantInterpolation((times,), PARData),
    OSat = Interpolations.interpolate((times,), OSatData, Interpolations.Gridded(Interpolations.Linear())),
    day = Interpolations.ConstantInterpolation((times,), days),
    GPP_daily = fill(GPP_daily_mu, nDays),
    ER_daily = fill(ER_daily_mu, nDays),
    kGas = exp(k600_meanlog),
    Q = Q,
    dayStart = dayStart,
    u0 = u0)

synthDataGenerator = simulate_data(
    DOdata = missing,
    ODE = problem,
    times = times,
    nDays = nDays,
    GPP_daily_mu = GPP_daily_mu,
    GPP_daily_sigma = GPP_daily_sigma,
    ER_daily_mu = ER_daily_mu,
    ER_daily_sigma = ER_daily_sigma,
    k600_meanlog = k600_meanlog,
    k600_sdlog = k600_sdlog,
    K02_conv = K02_conv,
    err_obs_iid_sigma_scale = err_obs_iid_sigma_scale,
    err_proc_iid_sigma_scale= err_proc_iid_sigma_scale)

synthData = synthDataGenerator()
plot(times, synthData, seriestype = :scatter, legend = false)

model = statisticalModel(
    nDays,
    GPP_daily_mu,
    GPP_daily_sigma,
    ER_daily_mu,
    ER_daily_sigma,
    k600_meanlog,
    k600_sdlog,
    K02_conv,
    err_obs_iid_sigma_scale,
    err_proc_iid_sigma_scale)

sampleMe = model(DOdata = synthData, ODE = problem, times = times)
setadbackend(:reversediff)
chains = sample(sampleMe, NUTS(0.65), MCMCThreads(), 100, 3)

# Evaluate chain mixing
gelmandiag(chains)

# Trace plot for GPP on day 2
df = DataFrame(chains)
@df df plot(
    :iteration,
    cols(4),
    group = (:chain),
    lw = 3,
    xguidefontsize=15,
    xlabel = "Iteration",
    ylabel = "GPP: Day 2",
    yguidefontsize = 15,
    legendtitle = "Chain",
    background_color = :black,
    color_palette = palette([colorant"#3D8D59", colorant"#77BDC6",colorant"#18537E"]),
    legend = :false)

# Plot posterior and prior for GPP on day 2
@df df density(
    cols(4),
    group = (:chain),
    lw = 3,
    xguidefontsize=15,
    xlabel = "GPP: Day 2",
    ylabel = "Probability Density",
    yguidefontsize = 15,
    legendtitle = "Chain",
    background_color = :black,
    color_palette = palette([colorant"#3D8D59", colorant"#77BDC6",colorant"#18537E"]),
    legend = :false)

max_stars_repo_name: elizabethjmohr/mesocosm-metabolism

max_issues_repo_name: elizabethjmohr/mesocosm-metabolism

max_forks_repo_name: elizabethjmohr/mesocosm-metabolism


text:
function formatinfo(s::Spec, x::AbstractString)
    size = ncodeunits(x) * sizeof(codeunit(x))
    width = length(x)
    return size + paddingsize(s, width), width
end

function formatfield(data::Vector{UInt8}, p::Int, s::Spec, x::AbstractString, width::Int)
    align = default(s.align, ALIGN_LEFT)
    pw = paddingwidth(s, width)
    p = padleft(data, p, s.fill, align, pw)
    n = ncodeunits(x)
    if isspecified(s.precision)
        n = min(nextind(x, 1, s.precision) - 1, n)
    end
    copyto!(data, p, codeunits(x), 1, n)
    p += n
    p = padright(data, p, s.fill, align, pw)
    return p
end


text:
module UniversalCredit
#
# This module implements the Universal Credit system,
# largely as its defined in CPAG 19/20 and 20/21.
#
using Base: Bool
using Dates: TimeType, now

using ScottishTaxBenefitModel

# FIXME prune these imports to things actually used.

using .ModelHousehold:
    BenefitUnit,
    Household,
    Person,
    empl_status_in,
    get_benefit_units,
    get_head,
    get_spouse,
    is_head,
    is_spouse,
    pers_is_disabled,
    pers_is_carer

using .STBParameters:
    UniversalCreditSys,
    AgeLimits,
    ChildLimits,
    HoursLimits,
    HousingRestrictions,
    MinimumWage,
    get_minimum_wage

using .Intermediate:
    MTIntermediate,
    HHIntermed,
    born_before,
    has_limited_capactity_for_work_activity,
    has_limited_capactity_for_work,
    make_recipient,
    reached_state_pension_age

using .Results:
    BenefitUnitResult,
    HouseholdResult,
    IndividualResult,
    UCResults,
    has_any,
    to_string,
    total

using .LocalLevelCalculations:
    apply_rent_restrictions

using .LegacyMeansTestedBenefits:
    num_qualifying_for_severe_disability,
    tariff_income

using .STBIncomes
using .Definitions

export
    basic_conditions_satisfied,
    calc_elements!,
    calc_standard_allowance,
    calc_tariff_income!,
    calc_uc_child_costs!,
    calc_uc_income!,
    calc_universal_credit!,
    disqualified_on_capital,
    make_min_se,
    qualifiying_16_17_yo

function pid_of_youngest_adult( bu :: BenefitUnit ) :: BigInt
    y = 99999
    yp :: BigInt = -1
    for pid in bu.adults
        if bu.people[pid].age < y
            y = bu.people[pid].age
            yp = pid
        end
    end
    return yp
end

"""
The test of basic eligibility from p37- 2020/1 CPAG
"""
function basic_conditions_satisfied(
    benefit_unit :: BenefitUnit,
    intermed     :: MTIntermediate,
    uc           :: UniversalCreditSys, # FIXME uc->ucsys everywhere
    age_limits   :: AgeLimits ) :: Bool
    bu = benefit_unit # shortcut
    if intermed.all_pension_age
        return false
    elseif intermed.age_oldest_adult < 18 # FIXME parameterise this
        q1617 = false
        for pid in benefit_unit.adults
            if qualifiying_16_17_yo( bu, bu.people[pid], intermed, uc )
                q1617 = true
                break
            end
        end
        return q1617
    else
        all_in_educ = true
        for pid in bu.adults
            # FIXME need a better test than this
            in_educ :: Bool =
                (bu.people[pid].employment_status == Student) &&
                (intermed.num_children == 0) &&
                (! pers_is_disabled( bu.people[pid]))
            if ! in_educ
                all_in_educ = false
                break
            end
        end
        return ! all_in_educ
    end
    @assert false "should never get to end of basic_conditions_satisfied"
end

function disqualified_on_capital(
    benefit_unit :: BenefitUnit,
    uc           :: UniversalCreditSys ) :: Bool
    cap = 0.0
    bu = benefit_unit # shortcut
    # FIXME we're doing this twice
    for pid in bu.adults
        if bu.people[pid].over_20_k_saving
            return true
        else
            for (at,val) in bu.people[pid].assets
                cap += val
            end
        end
    end
    return cap > uc.capital_max
end

function qualifiying_16_17_yo(
    benefit_unit   :: BenefitUnit,
    pers :: Person,
    intermed :: MTIntermediate,
    uc :: UniversalCreditSys ) :: Bool
    bu = benefit_unit # shortcut
    if pers_is_carer( pers )
        return true
    elseif has_limited_capactity_for_work( pers )
        return true
    elseif (is_head( bu, pers ) || is_spouse( bu, pers )) && ( intermed.num_children > 0)
        return true
        ## And some others; see CPAG 2020/1 p 37
    end
    return false
end

function calc_standard_allowance(
    benefit_unit     :: BenefitUnit,
    intermed         :: MTIntermediate,
    uc               :: UniversalCreditSys )
    sa = 0.0
    bu = benefit_unit # shortcut
    yp :: BigInt = pid_of_youngest_adult( bu )
    if intermed.num_adults == 1
        if intermed.age_oldest_adult < 18
            # under 18 single FIXME we should only every need 1 call to this, at the beginning
            if qualifiying_16_17_yo( bu, bu.people[yp], intermed, uc )
                sa = uc.age_18_24
            end
        elseif intermed.age_oldest_adult < 25
            # u 26 single
            sa = uc.age_18_24
        else
            # 25+ single
            sa = uc.age_25_and_over
        end
    else
        if  intermed.age_oldest_adult < 18 # FIXME parameterise these numbers, at least
            # under 18 couple
            np = 0
            for pid in bu.adults
                if qualifiying_16_17_yo( bu, bu.people[pid], intermed, uc )
                    np += 1
                end
            end
            if np == 1
                sa = uc.age_18_24
            elseif np == 2
                sa = uc.couple_both_under_25
            end
        elseif intermed.age_oldest_adult < 25
            # under 25 couple(probably - might be u25 single if 2nd adult under 18)
            if  intermed.age_youngest_adult >= 18
                sa = uc.couple_both_under_25
            else
                if qualifiying_16_17_yo( bu, bu.people[yp], intermed, uc )
                    sa = uc.couple_both_under_25
                else
                    sa = uc.age_18_24 # single person, ...
                end
            end
        else
            # 25+ couple (probably - might be 25+ single if 2nd adult under 18)
            if intermed.age_youngest_adult >= 18
                sa = uc.couple_oldest_25_plus
            else
                ## kinda TODO
                if qualifiying_16_17_yo( bu, bu.people[yp], intermed, uc )
                    sa = uc.couple_oldest_25_plus
                else
                    sa = uc.age_25_and_over # single person, ...
                end
            end
        end # over 25 couple
    end # 2 adult bu
    return sa
end # standard allowance

function calc_elements!(
    ucr          :: UCResults,
    benefit_unit :: BenefitUnit,
    intermed     :: MTIntermediate,
    uc           :: UniversalCreditSys,
    hours_limits :: HoursLimits,
    child_limits :: ChildLimits )
    bu = benefit_unit # shortcut

    # child elements

    if intermed.num_allowed_children > 0
        if born_before(
            intermed.age_oldest_child,
            child_limits.policy_start,
            now()) # FIXME we need to parameterise model_run_date somehow
            # this is the abolished 1st child thing see 20/21 p66
            ucr.child_element = uc.first_child
        else
            ucr.child_element = uc.subsequent_child
        end
        ucr.child_element += (intermed.num_allowed_children-1)*uc.subsequent_child
    end
    # limited capacity for work-related Activity
    # 1 per BU; see cpag 20/1 p 71-73
    lcwa = false
    lcwa_pid :: BigInt = -1
    for pid in bu.adults
        if has_limited_capactity_for_work_activity( bu.people[pid])
            lcwa = true
            lcwa_pid = pid
            break
        end
    end
    if lcwa
        ucr.limited_capacity_for_work_activity_element = uc.limited_capcacity_for_work_activity
    end
    carer = false
    for pid in bu.adults
        if pid != lcwa_pid # can't also give carers to same person, but could give to spouse
            if pers_is_carer( bu.people[pid] )
                carer = true
                break
            end # is carer
        end # pid
    end
    if carer
        ucr.carer_element = uc.carer
    end
end

"""

FIXME: Near dup of WTC calculation
"""
function calc_uc_child_costs!(
    ucr                 :: UCResults,
    benefit_unit        :: BenefitUnit,
    intermed            :: MTIntermediate,
    uc                  :: UniversalCreditSys )
    bu = benefit_unit # shortcut
    cost_of_childcare = 0.0
    for pid in bu.children
        cost_of_childcare += bu.people[pid].cost_of_childcare
    end
    cost_of_childcare *= uc.childcare_proportion
    if intermed.num_children > 1
        cost_of_childcare = min( uc.childcare_max_2_plus_children, cost_of_childcare )
    else
        cost_of_childcare = min( uc.childcare_max_1_child, cost_of_childcare )
    end
    ucr.childcare_costs = cost_of_childcare
end

function make_min_se(
    seinc    :: T,
    age      :: Integer,
    uc       :: UniversalCreditSys,
    minwage  :: MinimumWage ) :: T where T
    mw :: T = get_minimum_wage( minwage, age )
    min_se :: T = mw*uc.minimum_income_floor_hours
    return max( min_se, seinc )
end

"""
Implements CPAG 19/20 ch 7
"""
function calc_uc_income!(
    benefit_unit_result :: BenefitUnitResult,
    benefit_unit        :: BenefitUnit,
    intermed            :: MTIntermediate,
    uc                  :: UniversalCreditSys,
    minwage             :: MinimumWage )
    bu = benefit_unit # shortcut
    bur = benefit_unit_result # shortcut
    inc = 0.0
    earn = 0.0
    for pid in bu.adults
        inc += isum( bur.pers[pid].income, uc.other_income )
        earn += isum( bur.pers[pid].income, uc.earned_income )
        seinc = bur.pers[pid].income[SELF_EMPLOYMENT_INCOME]
        # FIXME any self employed? what about losses?
        if bu.people[pid].employment_status == Full_time_Self_Employed
            seinc = make_min_se( seinc, bu.people[pid].age, uc, minwage )
        end
        earn += seinc
    end
    bur.uc.earnings_before_allowances = earn # the PiP calculator uses this for CTR for UC recipients
    if( intermed.num_children > 0 ) || intermed.limited_capacity_for_work
        bur.uc.work_allowance = bur.uc.housing_element > 0 ?
            uc.work_allowance_w_housing : uc.work_allowance_no_housing
        earn = max( 0.0, earn-bur.uc.work_allowance)
    end
    bur.uc.other_income = inc
    bur.uc.earned_income = earn*uc.taper
end

## FIXME we need the extra capital var here benunit.Totsav
function calc_tariff_income!(
    benefit_unit_result :: BenefitUnitResult,
    benefit_unit :: BenefitUnit,
    uc           :: UniversalCreditSys )
    bu = benefit_unit # shortcut
    ucr = benefit_unit_result.uc # shortcut
    cap = 0.0
    # FIXME we're doing this twice!
    for pid in bu.adults
        for (at,val) in bu.people[pid].assets
            cap += val
        end
    end
    ucr.assets = cap
    ucr.tariff_income = tariff_income( cap, uc.capital_min, uc.capital_tariff )
 end

#
# FIXME CPAG 19/20 p 122 shows a earned/unearned
# split calculation with earnings
#
function calc_universal_credit!(
    benefit_unit_result :: BenefitUnitResult,
    benefit_unit     :: BenefitUnit,
    intermed         :: MTIntermediate,
    uc               :: UniversalCreditSys,
    age_limits       :: AgeLimits,
    hours_limits     :: HoursLimits,
    child_limits     :: ChildLimits,
    hr               :: HousingRestrictions,
    minwage          :: MinimumWage )
    bu = benefit_unit # shortcut
    bur = benefit_unit_result # shortcut

    if ! basic_conditions_satisfied(
        bu,
        intermed,
        uc,
        age_limits )
        bur.uc.basic_conditions_satisfied = false
        return
    end
    bur.uc.basic_conditions_satisfied = true
    if disqualified_on_capital( bu, uc )
        bur.uc.disqualified_on_capital = true
        return
    end
    bur.uc.disqualified_on_capital = false
    bur.uc.standard_allowance = calc_standard_allowance( benefit_unit, intermed, uc )
    calc_elements!( bur.uc, bu, intermed, uc, hours_limits, child_limits )
    if( intermed.num_working_ft > 0 ) || ( intermed.num_working_pt > 0 )
        calc_uc_child_costs!( bur.uc, bu, intermed, uc )
    end
    calc_uc_income!( bur, bu, intermed, uc, minwage )
    calc_tariff_income!( bur, bu, uc )
    bur.uc.maximum =
        bur.uc.standard_allowance +
        bur.uc.limited_capacity_for_work_activity_element +
        bur.uc.child_element +
        bur.uc.housing_element +
        bur.uc.carer_element +
        bur.uc.childcare_costs
    bur.uc.total_income =
        bur.uc.earned_income -
        bur.uc.other_income -
        bur.uc.tariff_income
    uce = max( 0.0, bur.uc.maximum - bur.uc.total_income )

    # Make the recipient the bu head if the head isn't
    # retired.
    head = get_head( bu )
    target_pid = head.pid
    if reached_state_pension_age( age_limits, head.age, head.sex )
        target_pid = get_spouse( bu ).pid # fail if spouse is nothing - but this must be the case
    end
    bur.uc.recipient= target_pid # save a record of who gets for e.g. CT calculation
    bur.pers[target_pid].income[UNIVERSAL_CREDIT] = uce
end

function calc_uc_housing_element!(
    household_result :: HouseholdResult,
    household        :: Household,
    intermed         :: HHIntermed,
    uc               :: UniversalCreditSys,
    hr               :: HousingRestrictions
)
    hhr = household_result # shortcut
    hh = household # shortcut
    eligible_amount = 0.0
    if renter( hh.tenure )
        eligible_amount = household_result.housing.allowed_rent
    elseif ! has_any( hhr.bus[1], WAGES,SELF_EMPLOYMENT_INCOME )
        eligible_amount = household.other_housing_charges  # FIXME go through the list
    end
    bus = get_benefit_units(household)
    nbus = size(bus)[1]
    ndds = 0.0
    if nbus > 1
        if num_qualifying_for_severe_disability(
            bus[1],
            hhr.bus[1],
            1 ) == 0
            for buno in 2:nbus
                # children under 5
                buint = intermed.buint[buno] # shortcut
                if (buint.num_children == 0) || (buint.age_youngest_child >= 5)
                    for adno in bus[buno].adults
                        pers = bus[buno].people[adno]
                        pr = hhr.bus[buno].pers[adno] # fixme make this `people` to match the bu
                        exempt =  # FIXME parameterise these
                            pers.age < 21 ||
                            any_positive( pr.income,
                                [ATTENDANCE_ALLOWANCE,
                                PENSION_CREDIT,
                                PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING,
                                DLA_SELF_CARE,
                                CARERS_ALLOWANCE,
                                SCOTTISH_CARERS_SUPPLEMENT,
                                SCOTTISH_DISABILITY_ASSISTANCE_CHILDREN_DAILY_LIVING,
                                SCOTTISH_DISABILITY_ASSISTANCE_OLDER_PEOPLE,
                                SCOTTISH_DISABILITY_ASSISTANCE_WORKING_AGE_DAILY_LIVING] )
                        if ! exempt
                            ndds += uc.ndd
                        end
                    end # adult loop
                end # not responsible for u5 child
            end # bu loop
        end # bu 1 isn't severe disabled
    end # > 1 bu
    hhr.bus[1].uc.housing_element =
        max(0.0, eligible_amount - ndds )
end

function calc_universal_credit!(
    household_result :: HouseholdResult,
    household        :: Household,
    intermed         :: HHIntermed,
    uc               :: UniversalCreditSys,
    age_limits       :: AgeLimits,
    hours_limits     :: HoursLimits,
    child_limits     :: ChildLimits,
    hr               :: HousingRestrictions,
    minwage          :: MinimumWage )
    # fixme duped with LMTBens
    hhr = household_result
    hh = household # shortcuts
    if uc.abolished
        return
    end
    household_result.housing = apply_rent_restrictions(
        household, intermed.hhint, hr )
    bus = get_benefit_units(household)
    nbus = size( bus )[1]
    #
    # Do housing 1st since it needs a seperate
    # loop round the bus, but we can do everything else
    # bu by bu. housing element is assigned always
    # to the first benefit unit; everyone else
    # contributes NDDs.
    #
    calc_uc_housing_element!(
        hhr,
        hh,
        intermed,
        uc,
        hr
    )
    for buno in eachindex(bus)
        calc_universal_credit!(
            hhr.bus[buno],
            bus[buno],
            intermed.buint[buno],
            uc,
            age_limits,
            hours_limits,
            child_limits,
            hr,
            minwage
        )
    end

    # hack - council tax rebates for UC Retail_trade_except_of_motor_vehicles_and_motorcycles
    # 20% of any income above Maximum Universal Credit is deducted from maximum CT support
    # income seems to *include* universal credit, and full wages
    # so grossed back up from the tapered calculated wages above, but
    # otherwise following the same rules. Something like that, anyway.
    #
    recipient = hhr.bus[1].uc.recipient
    if recipient > 0 # some *something* done with UC
        ucrec = hhr.bus[1].pers[recipient].income[UNIVERSAL_CREDIT]
        if ucrec > 0
            bur = household_result.bus[1]
            ct = total(household_result, LOCAL_TAXES )
            # grossed_up_earn = (bur.uc.earned_income/uc.taper)
            ucincome =
                bur.uc.earnings_before_allowances + # gross earned income back up
                bur.uc.other_income +
                bur.uc.tariff_income +
                ucrec
            excess = max(0.0, ucincome - bur.uc.maximum)
            if excess > 0
                ct = max( 0.0, ct - excess*uc.ctr_taper )
            end
            bur.uc.ctr = ct
        end
    end
end

end # module


text:
using Test, Printf, ForwardDiff

using SingleColumnModels.FiniteDifferenceGrids

@testset "Grid interface" begin
    for FT in (Float32, Float64)
        for n_ghost in (1, 2)
            n_elems = 12
            n_elems_real = n_elems - 2 * n_ghost
            elem_indexes = 1:n_elems
            elem_indexes_real = elem_indexes[(1 + n_ghost):(end - n_ghost)]
            Δz = FT(1 / n_elems_real)
            grid = UniformGrid(FT(0.0), FT(1.0), n_elems_real, n_ghost)
            @test n_hat(Zmin()) == -1
            @test n_hat(Zmax()) == 1
            @test binary(Zmin()) == 0
            @test binary(Zmax()) == 1
            @test ghost_dual(Zmin()) == [1, 0]
            @test ghost_dual(Zmax()) == [0, 1]
            @test eltype(grid) == FT
            @test all(over_elems(grid) .== elem_indexes)
            @test all(over_elems_real(grid) .== elem_indexes_real)
            @test length(over_elems(grid)) == n_elems
            @test length(over_elems_real(grid)) == n_elems_real
            @test first_interior(grid, Zmax()) == n_elems - n_ghost
            @test first_interior(grid, Zmin()) == 1 + n_ghost
            @test boundary(grid, Zmin()) == 1 + n_ghost
            @test boundary(grid, Zmax()) == n_elems - n_ghost
            @test over_elems_ghost(grid) ==
                  [(1:n_ghost)..., ((n_elems + 1 - n_ghost):n_elems)...]
            @test grid.zn_min ≈ FT(0.0)
            @test grid.zn_max ≈ FT(1.0)
            sprint(show, grid)
        end
    end
end

function ∇(f, x::T) where {T}
    tag = typeof(ForwardDiff.Tag(f, T))
    y = f(ForwardDiff.Dual{tag}(x, one(x)))
    ForwardDiff.partials(tag, y, 1)
end

@testset "Grid operators" begin
    for FT in (Float32, Float64)
        n_ghost = 1
        n_elems = 12
        n_elems_real = n_elems - 2 * n_ghost
        elem_indexes = 1:n_elems
        elem_indexes_real = elem_indexes[(1 + n_ghost):(end - n_ghost)]
        Δz = FT(1 / n_elems_real)
        grid = UniformGrid(FT(0.0), FT(1.0), n_elems_real, n_ghost)
        f = FT[1, 5, 3]
        w = FT[2, 3, 8]
        @test grad(f[1:2], grid) ≈ FT(40)
        @test grad(f, grid) ≈ FT(10)
        @test_throws AssertionError grad(f[1], grid)
        @test_throws AssertionError advect(f[1:2], w[1:2], grid)
        @test advect(f, w, grid) ≈ FT(40)
        @test advect(f, -w, grid) ≈ FT(-20)
        @test ∇_up(f, grid) ≈ FT(40)
        @test ∇_dn(f, grid) ≈ FT(-20)
    end
end

@testset "Grid operators convergence" begin
    z_min = 0.0
    z_max = 2.0 * π
    tol = 1 / 2
    ϕ_(x) = 5 + sin(x)
    ∇ϕ_(x) = cos(x) # computed manually, since ForwardDiff doesn't like gradient functions as inputs.
    ψ_(x) = 5 + cos(x)
    w_(x) = 5 + sin(x)
    K_(x) = 5 + sin(x)
    ϕψ_(x) = ψ_(x) * ϕ_(x)
    wϕ_(x) = w_(x) * ϕ_(x)
    K∇ϕ_(x) = K_(x) * ∇ϕ_(x)

    for n_elems_real in [2^k for k in 3:8]
        grid = UniformGrid(z_min, z_max, n_elems_real)
        z = grid.zc

        ϕ = ϕ_.(z)
        ψ = ψ_.(z)
        ϕψ = ϕψ_.(z)
        w = w_.(z)
        K = K_.(z)

        ∇ϕ = ∇.(ϕ_, z)
        ∇ψ = ∇.(ψ_, z)
        ∇ϕψ = ∇.(ϕψ_, z)
        ∇w = ∇.(w_, z)
        ∇K∇ϕ = ∇.(K∇ϕ_, z)
        ∇ϕ = ∇.(ϕ_, z)
        ∇wϕ = ∇.(wϕ_, z)

        del_err = [
            abs(
                ∇ϕ[k] -
                ∇_z_upwind(ϕ[(k - 1):(k + 1)], w[(k - 1):(k + 1)], grid),
            ) for k in over_elems_real(grid)
        ]
        @test all([e < tol * grid.Δz for e in del_err])

        del_err = [
            abs((∇ϕ[k - 1] + ∇ϕ[k]) / 2 - ∇_z_flux(ϕ[(k - 1):k], grid))
            for k in over_elems_real(grid)
        ]
        @test all([e < tol * grid.Δz^2 for e in del_err])
        del_err = [
            abs((∇ϕ[k] + ∇ϕ[k + 1]) / 2 - ∇_z_flux(ϕ[k:(k + 1)], grid))
            for k in over_elems_real(grid)
        ]
        @test all([e < tol * grid.Δz^2 for e in del_err])

        del_err = [
            abs(∇ϕ[k] - ∇_z_centered(ϕ[(k - 1):(k + 1)], grid))
            for k in over_elems_real(grid)
        ]
        @test all([e < tol * grid.Δz^2 for e in del_err])

        del_err = [
            abs(
                (∇ϕ[k - 1] + ∇ϕ[k]) / 2 - ∇_z_dual(ϕ[(k - 1):(k + 1)], grid)[1],
            ) for k in over_elems_real(grid)
        ]
        @test all([e < tol * grid.Δz^2 for e in del_err])
        del_err = [
            abs(
                (∇ϕ[k] + ∇ϕ[k + 1]) / 2 - ∇_z_dual(ϕ[(k - 1):(k + 1)], grid)[2],
            ) for k in over_elems_real(grid)
        ]
        @test all([e < tol * grid.Δz^2 for e in del_err])

        Lap_err = [
            abs(∇ψ[k] - Δ_z(ϕ[(k - 1):(k + 1)], grid))
            for k in over_elems_real(grid)
        ]
        @test all([e < tol * grid.Δz^2 for e in Lap_err])

        Lap_err = [
            abs(∇K∇ϕ[k] - Δ_z(ϕ[(k - 1):(k + 1)], grid, K[(k - 1):(k + 1)])) for k in over_elems_real(grid)
        ]
        @test all([e < 2 * tol * grid.Δz^2 for e in Lap_err])

        Lap_err = [
            abs(
                ∇K∇ϕ[k] - Δ_z_dual(
                    ϕ[(k - 1):(k + 1)],
                    grid,
                    0.5 * (K[(k - 1):k] .+ K[k:(k + 1)]),
                ),
            ) for k in over_elems_real(grid)
        ]
        @test all([e < 2 * tol * grid.Δz^2 for e in Lap_err])
        Δt = 1.0

        # ∇•(wϕ)
        for conservatove_scheme in
            (UpwindCollocated(), OneSidedUp(), OneSidedDn(), CenteredUnstable())
            adv_err = [
                abs(
                    ∇wϕ[k] - advect_old(
                        ϕ[(k - 1):(k + 1)],
                        0.5 * (ϕ[(k - 1):k] .+ ϕ[k:(k + 1)]),
                        w[(k - 1):(k + 1)],
                        0.5 * (w[(k - 1):k] .+ w[k:(k + 1)]),
                        grid,
                        conservatove_scheme,
                        Δt,
                    ),
                ) for k in over_elems_real(grid)
            ]
            @test all([e < 12 * tol * grid.Δz for e in adv_err])
        end
        # w•∇ϕ
        for non_conservative_scheme in (UpwindAdvective(),)
            adv_err = [
                abs(
                    w[k] * ∇ϕ[k] - advect_old(
                        ϕ[(k - 1):(k + 1)],
                        0.5 * (ϕ[(k - 1):k] .+ ϕ[k:(k + 1)]),
                        w[(k - 1):(k + 1)],
                        0.5 * (w[(k - 1):k] .+ w[k:(k + 1)]),
                        grid,
                        non_conservative_scheme,
                        Δt,
                    ),
                ) for k in over_elems_real(grid)
            ]
            @test all([e < 6 * tol * grid.Δz for e in adv_err])
        end

    end
end


text:
s = """
throw(MethodError("blah"))
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
MethodError("blah")
"""
msgs = lintstr(s)
@test msgs[1].code == :W448
@test contains(msgs[1].message, "MethodError is an Exception but it is not enclosed in a throw")


text:
function _gram(kernel::RK,
               knots::Vector{T}) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
    n_1 = length(knots)
    mat = Matrix{T}(undef, n_1, n_1)
    @inbounds for j = 1:n_1
        for i = 1:j
            mat[i,j] = _rk(kernel, knots[i], knots[j])
            mat[j,i] = mat[i,j]
        end
    end
    return mat
end

function _gram_d1(kernel::RK,
                  knots::Vector{T},
                  d1_knots::Vector{T}) where {T <: AbstractFloat, RK <: ReproducingKernel_2}
    n_1 = length(knots)
    n_2 = length(d1_knots)
    n = n_1 + n_2
    mat = Matrix{T}(undef, n, n)
    @inbounds for j = 1:n_1
        for i = 1:j
            mat[i,j] = _rk(kernel, knots[i], knots[j])
            mat[j,i] = mat[i,j]
        end
    end
    n_1_p1 = n_1 + 1
    @inbounds for j = n_1_p1:n
        for i = 1:n_1
            mat[i,j] = _∂rk_∂ξ(kernel, knots[i], d1_knots[j-n_1])
            mat[j,i] = mat[i,j]
        end
    end
    @inbounds for j = n_1_p1:n
        for i = j:n
            mat[i,j] = _∂²rk_∂η∂ξ(kernel, d1_knots[i-n_1], d1_knots[j-n_1])
            mat[j,i] = mat[i,j]
        end
    end
    return mat
end

function _gram_d2(kernel::RK,
                  knots::Vector{T},
                  d2_knots::Vector{T}) where {T <: AbstractFloat, RK <: ReproducingKernel_3}
    n_1 = length(knots)
    n_2 = length(d2_knots)
    n = n_1 + n_2
    mat = Matrix{T}(undef, n, n)
    @inbounds for i = 1:n_1
        for j = 1:i
            mat[i,j] = _rk(kernel, knots[i], knots[j])
            mat[j,i] = mat[i,j]
        end
    end
    n_1_p1 = n_1 + 1
    @inbounds for i = 1:n_1
        for j = n_1_p1:n
            mat[i,j] = _∂²rk_∂ξ²(kernel, knots[i], d2_knots[j-n_1])
            mat[j,i] = mat[i,j]
        end
    end
    @inbounds for i = n_1_p1:n
        for j = i:n
            mat[i,j] = _∂⁴rk_∂η²∂ξ²(kernel, d2_knots[i-n_1], d2_knots[j-n_1])
            mat[j,i] = mat[i,j]
        end
    end
    return mat
end

function _gram12(kernel::RK,
               knots::Vector{T},
               d1_knots::Vector{T},
               d2_knots::Vector{T}) where {T <: AbstractFloat, RK <: ReproducingKernel_3}
    n_1 = length(knots)
    n_2 = length(d1_knots)
    n_3 = length(d2_knots)
    n = n_1 + n_2 + n_3
    mat = Matrix{T}(undef, n, n)
    @inbounds for j = 1:n_1
        for i = 1:j
            mat[i,j] = _rk(kernel, knots[i], knots[j])
            mat[j,i] = mat[i,j]
        end
    end
    n_1_p1 = n_1 + 1
    ne = n_1 + n_2
    @inbounds for j = n_1_p1:ne
        for i = 1:n_1
            mat[i,j] = _∂rk_∂ξ(kernel, knots[i], d1_knots[j-n_1])
            mat[j,i] = mat[i,j]
        end
    end
    ns = n_1 + n_2 + 1
    @inbounds for j = ns:n
        for i = 1:n_1
            mat[i,j] = _∂²rk_∂ξ²(kernel, knots[i], d2_knots[j-ns+1])
            mat[j,i] = mat[i,j]
        end
    end
    ne = n_1 + n_2
    @inbounds for j = n_1_p1:ne
        for i = j:ne
            mat[i,j] = _∂²rk_∂η∂ξ(kernel, d1_knots[i-n_1], d1_knots[j-n_1])
            mat[j,i] = mat[i,j]
        end
    end
    ns = n_1 + n_2 + 1
    ne = n_1 + n_2
    @inbounds for j = ns:n
        for i = n_1_p1:ne
            mat[i,j] = _∂³rk_∂η∂ξ²(kernel, d1_knots[i-n_1], d2_knots[j-ns+1])
            mat[j,i] = mat[i,j]
        end
    end
    ns = n_1 + n_2 + 1
    @inbounds for j = ns:n
        for i = j:n
            mat[i,j] = _∂⁴rk_∂η²∂ξ²(kernel, d2_knots[i-ns+1], d2_knots[j-ns+1])
            mat[j,i] = mat[i,j]
        end
    end
    return mat
end


text:
# =========================================================================================
"""
    readpqr{T <: AbstractFloat}(
        stream::IOStream,
              ::Type{T}=Float64
    )

Reads a charge model from the given PQR file.

# Specification
<https://pdb2pqr.readthedocs.io/en/latest/formats/pqr.html>

# Return type
`Vector{Charge{T}}`

# Alias

    readpqr{T}(fname::String, ::Type{T}=Float64)

Reads the charge model using a file name rather than a `IOStream` object.
"""
function readpqr(
        stream::IOStream,
              ::Type{T}=Float64
    ) where T <: AbstractFloat
    charges = Charge{T}[]
    for line in eachline(stream)
        startswith(line, "ATOM") || continue # skip remarks and water molecules
        vals = map(e -> parse(T, e), split(line)[end-4:end-1])
        vals[4] == 0 && continue             # skip zero charges
        push!(charges, Charge(vals...))
    end
    charges
end

function readpqr(
        fname::String,
             ::Type{T}=Float64
    ) where T
    open(fh -> readpqr(fh, T), fname)
end


text: default(titlefontsize = 16, legendfontsize = 12, guidefontsize = 16, tickfontsize = 14, linewidth=2, fontfamily = "serif")


text:
using MetidaFreq
using Test
using DataFrames, CSV

path     = dirname(@__FILE__)
io       = IOBuffer();
freqdat  = CSV.File(path*"/csv/freqdat.csv") |> DataFrame

MetidaFreq.contab(freqdat, :row, :col)

freqtable(freqdat, :row, :col)

v = vcat(freqdat, freqdat)
for i = 1:10
    v = vcat(v, freqdat)
end

using Profile
using PProf
using BenchmarkTools

pprof()
PProf.kill()
Profile.clear()

@profile  for i = 1:10000  MetidaFreq.contab(v, :row, :col) end
@profile  for i = 1:10000  freqtable(v, :row, :col) end
@benchmark MetidaFreq.contab($v, :row, :col)
@benchmark freqtable($v, :row, :col)

freqdat = append!(freqdat, freqdat)
a = Vector{String}(undef, 86)
fill!(view(a, 1:40), "a")
fill!(view(a, 41:86), "b")
freqdat.s1 = a
freqdat = append!(freqdat, freqdat)
a = Vector{String}(undef, 172)
fill!(view(a, 1:100), "c")
fill!(view(a, 101:172), "d")
freqdat.s2 = a
MetidaFreq.contab(freqdat, :row, :col; sort = :s1)
MetidaFreq.contab(freqdat, :row, :col)

freqdat.s2 = Symbol.(freqdat.s2)


text:
module Persist

# TODO: Implement Arrays
# TODO: Handle Circular References
# TODO: Cleanup loading and storing code
#   - Review existing constructs - there may be a simpler implementation
#   - Make generic persisting and storing type stable, either through code tricks or by
#       using @generated.
#
# TODO: write tests!

export Persistent, persist, retrieve

import Base: getproperty, setproperty!, size, getindex, setindex!, IndexStyle, sizeof

using ..Lib
using ..Transaction

abstract type AbstractClass end

struct Fallback <: AbstractClass end
struct IsBits <: AbstractClass end

class(::Type{T}) where {T} = isbitstype(T) ? IsBits() : Fallback()
class(x::T) where {T} = class(T)

"""
Persistent wrapper for (mostly) arbitrary types. The layout of a `Persistent{T}` is similar
similar to the layout of standard Julia object with the following modifications

* If a field element is `isbits`, the element is stored inline
* If a field element is not `isbits`, then we store a `PersistentOID` to that object. This
    means that fields offsets below the non `isbits` object may have to be shifted down to
    accomodate the 16 bytes required for a `PersistentOID`.
"""
struct Persistent{T}
    __oid::PersistentOID{T}
end
getoid(P::Persistent) = getfield(P, :__oid)
getpool(P::Persistent) = Lib.pool(getoid(P))

#####
##### persist
#####

persist(pool, x::T) where {T} = transaction(() -> store(pool, x), pool)

#####
##### Store
#####

# Special behavior for strings and arrays
store(pool, x::String) = Persistent(Lib.tx_strdup(x, String))
store(pool, x::Symbol) = Persistent(Lib.tx_strdup(String(x), Symbol))

# arrays
store(pool, x::Array{T,N}) where {T,N} = _store(pool, x, Val{isbitstype(T)}())

# Fallback - check if isbits
store(pool, x::T) where {T} = store(pool, x::T, class(T))
function store(pool, x::T, ::IsBits) where {T}
    oid = Lib.tx_alloc(sizeof(x), T)
    ptr = Lib.direct(oid)
    unsafe_store!(ptr, x)

    return Persistent(oid)
end

function store(pool, x::T, ::Fallback) where {T}
    # Create OIDs for every field of `x`
    oid = Lib.tx_alloc(nfields(x) * sizeof(Persistent{T}), T)
    ptr = Lib.direct(oid)

    # Store all members of `x`
    for (index, field) in enumerate(fieldnames(T))
        object = store(pool, getfield(x, field))
        thisptr = convert(Ptr{typeof(object)}, ptr) + (index - 1) * sizeof(Persistent{T})
        unsafe_store!(thisptr, object)
    end
    return Persistent{T}(oid)
end

#####
##### retrieve
#####

retrieve(x::Persistent) = retrieve(Lib.direct(getoid(x)))

# Specializtions
retrieve(ptr::Ptr{String}) = unsafe_string(convert(Ptr{UInt8}, ptr))
retrieve(ptr::Ptr{Symbol}) = Symbol(unsafe_string(convert(Ptr{UInt8}, ptr)))

# Fallback
retrieve(ptr::Ptr{T}) where {T} = retrieve(ptr, class(T))
retrieve(ptr::Ptr{T}, ::IsBits) where {T} = unsafe_load(ptr)

function retrieve(ptr::Ptr{T}, ::Fallback) where {T}
    fields = []
    for (index, field) in enumerate(fieldnames(T))
        # Need to make an OID from the 16 bytes here to get the actual object.
        ptr_offset = ptr + (index - 1) * sizeof(Persistent{T})
        handle = unsafe_load(convert(Ptr{Persistent{fieldtype(T, field)}}, ptr_offset))

        push!(fields, retrieve(handle))
    end
    return T(fields...)
end

#####
##### destroy
#####

destroy(x::Persistent) = destroy(getoid(x))

# Specializations
destroy(oid::PersistentOID{String}) = Lib.free(oid)
destroy(oid::PersistentOID{Symbol}) = Lib.free(oid)

# Fallback
destroy(oid::PersistentOID{T}) where {T} = destroy(oid, class(T))
destroy(oid::PersistentOID{T}, ::IsBits) where {T} = Lib.free(oid)

end


text:
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using Printf

"""
    ShowExampleAppAutoResize(p_open::Ref{Bool})
Create a window which gets auto-resized according to its content.
"""
function ShowExampleAppAutoResize(p_open::Ref{Bool})
    CImGui.Begin("Example: Auto-resizing window", p_open, CImGui.ImGuiWindowFlags_AlwaysAutoResize) || (CImGui.End(); return)

    @cstatic lines=Cint(10) begin
        CImGui.Text("Window will resize every-frame to the size of its content.\nNote that you probably don't want to query the window size to\noutput your content because that would create a feedback loop.")
        @c CImGui.SliderInt("Number of lines", &lines, 1, 20)
        for i = 1:lines
            # pad with space to extend size horizontally
            CImGui.Text("$(" "^4i)This is line $i")
        end
        CImGui.End()
    end # @cstatic
end


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "flamingo", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# get hopper model
s_sim = get_simulation("flamingo", "slope_smooth_2D_lc", "slope")
s = get_simulation("flamingo", "flat_2D_lc", "flat")
model = s.model
env = s.env
const ContactImplicitMPC = Main
ref_traj_ = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/flamingo/gaits/gait_forward_36_4.jld2"),
    load_type = :split_traj_alt))
ref_traj = deepcopy(ref_traj_)

H = ref_traj.H
h = ref_traj.h
N_sample = 5
H_mpc = 15
h_sim = h / N_sample
H_sim = 10000 #35000

# barrier parameter
κ_mpc = 2.0e-4

obj = TrackingVelocityObjective(model, env, H_mpc,
    v = [Diagonal(1e-3 * [1e0,1,1e4,1,1,1,1,1e4,1e4]) for t = 1:H_mpc],
    q = [Diagonal(1e-1 * [3e2, 1e-6, 3e2, 1, 1, 1, 1, 0.1, 0.1]) for t = 1:H_mpc],
    u = [Diagonal(3e-1 * [0.1; 0.1; 0.3; 0.3; ones(model.dim.u-6); 2; 2]) for t = 1:H_mpc],
    γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
    b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

p = linearized_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
	ip_type = :interior_point,
	mode = :configuration,
	ip_opts = InteriorPointOptions(
					undercut = 5.0,
					κ_tol = κ_mpc,
					r_tol = 1.0e-8,
					diff_sol = true,
					solver = :empty_solver,
					max_time = 1e5,),
    n_opts = NewtonOptions(
        r_tol = 3e-4,
		# verbose = true,
        max_iter = 5),
    mpc_opts = LinearizedMPCOptions(
        # live_plotting=true,
        altitude_update = true,
        altitude_impact_threshold = 0.02,
        # altitude_verbose = true,
        )
    )

q1_ref = copy(ref_traj.q[2])
q0_ref = copy(ref_traj.q[1])
q1_sim = SVector{model.dim.q}(q1_ref)
q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

sim = simulator(s_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
	ip_opts = InteriorPointOptions(
		undercut = Inf,
		γ_reg = 0.0,
		# verbose = true,
        r_tol = 1.0e-8,
        κ_tol = 1.0e-8),
    sim_opts = SimulatorOptions(warmstart = true),
    )

@time status = simulate!(sim, verbose = true)
l = 9
lu = 1
plt = plot(layout=(3,1), legend=false)
plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[1,1], hcat(Vector.([q[l:l] for q in sim.traj.q])...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][lu:lu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[3,1], hcat(Vector.(vcat([fill(ref_traj.γ[i][1:nc], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[lu:lu] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
# plot!(plt[3,1], hcat(Vector.([γ[1:nc] for γ in sim.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)
# plot!(plt[3,1], hcat(Vector.([b[1:nb] for b in sim.traj.b]*N_sample)...)', color=:red, linewidth=1.0)

plot_surface!(vis, s_sim.env, xlims=[-0.5, 10.5], ylims = [-0.5, 0.5])
plot_lines!(vis, s_sim.model, sim.traj.q)
# anim = visualize_robot!(vis, model, sim.traj, sample=10)
anim = visualize_meshrobot!(vis, s_sim.model, sim.traj, sample=10)
anim = visualize_force!(vis, s_sim.model, s_sim.env, sim.traj, anim=anim, h=h_sim, sample=10)

(s_sim.env.surf(5.0) - s_sim.env.surf(4.0)) / (5.0 - 4.0) * 360 / (2π)

# filename = "flamingo_40_steps_sine"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/$filename.mp4", overwrite=true)
#
# convert_video_to_gif(
#     "/home/simon/Documents/$filename.mp4",
#     "/home/simon/Documents/$filename.gif", overwrite=true)
#
# flamingo_ghost!(vis, sim, model_sim.env.surf)
# flamingo_animation!(vis, sim, model_sim.env.surf)


text:
"""Julia program to implement Painting Fence algorithm
Given a fence with n posts and k colors
Find out the number of ways of painting the fence
such that at most 2 adjacent posts have the same color.

"""

function painting_fence(n, k)
    if (k < (n ÷ 2) || n < 2)
        return -1
    end
    tot = [0 for i in 1:(n - 1)]
    same = [0 for i in 1:(n - 1)]
    diff = [0 for i in 1:(n - 1)]

    same[1] = k
    diff[1] = k * (k - 1)
    tot[1] = k * k
    for i in 2:(n-1)
        same[i] = diff[i-1]
        diff[i] = tot[i-1] * (k-1)
        tot[i] = same[i] + diff[i]
    end
    return tot[n-1]
end

print("Enter the number of posts: ")
n = readline()
n = parse(Int, n)
print("Enter the number of colours available: ")
k = readline()
k = parse(Int, k)
ans = painting_fence(n, k)
if( ans == -1)
    println("The given posts cannot be painted with the given colours satisfying the given conditions")
else
    println("The number of ways to paint the given posts is $ans")
end

"""
Time Complexity: O(n), where 'n' is the number of posts
Space Complexity: O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE 1

Enter the number of posts: 10
Enter the number of colours available: 4
The given posts cannot be painted with the given colours satisfying the given conditions

SAMPLE 2

Enter the number of posts: 3
Enter the number of colours available: 5
The number of ways to paint the given posts is 120
"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
# Thompson sampling.
# Based on https://arxiv.org/abs/1111.1797.

mutable struct ThompsonSampling <: Policy end

struct ThompsonSamplingDetails <: PolicyDetails
  solver_time::Float64
end

function choose_action(instance::CombinatorialInstance{T}, policy::ThompsonSampling, state::State{T}; with_trace::Bool=false) where T
  # Determine the weights for each arm and use them to solve the combinatorial problem.
  # println("Counts: $(state.arm_counts); rewards: $(state.arm_average_reward)")
  weights = Dict(arm => rand(Beta(state.arm_counts[arm] * state.arm_average_reward[arm] + 1,
                                  state.arm_counts[arm] * (1 - state.arm_average_reward[arm]) + 1)) for arm in keys(state.arm_counts))

  t0 = time_ns()
  sol = solve_linear(instance, weights)
  t1 = time_ns()

  if with_trace
    return sol, ThompsonSamplingDetails((t1 - t0) / 1_000_000_000)
  else
    return sol
  end
end


text:
module TestPkgLH

#using ConfigLH

export print_version

pkg_path() = normpath(joinpath(@__DIR__, ".."));

function print_version()
	verNo = 43;
	println("TestPkg version $verNo")
	return verNo
end

end # module


text:
# Energy spectrum of oscillations at a fixed point.

using FFTW, JLD2, CurveFit, PyPlot
using Vlasiator: RE

file = "satellites_uniform_sampled.jld2"

data = JLD2.load(file)

nSatellite = length(data["t"])
nI, nJ = size(data["rho"])[2:3]

t = data["t"]
# Select spatial point
i, j = 5, 5
var = data["rho"][:,i,j]

dt = t[2] - t[1] # uniform sample interval [s]
Fs = 1 / dt      # sample frequency, [Hz]
Fn = Fs / 2      # Nyquist frequency, [Hz]

## Frequency calculation

nPoints = length(var)
nFFT = nPoints
df = Fs / nFFT
freq_fullrange = -Fn:df:Fn

freq = freq_fullrange[(nPoints ÷ 2 + 1):end-1]

var_freq = fft(var)
var_power = abs.(fftshift(var_freq))[(nPoints ÷ 2 + 1):end]

# k is the exponential coefficient
a, k = @views power_fit(freq[10:end], var_power[10:end])

figure(figsize=(6,8))
loglog(freq, var_power)
axis("scaled")

min_power, max_power = extrema(@view var_power[10:end])
xlim(freq[8], Fs)
ylim(min_power * 0.75, max_power * 2.0)

xlabel("Frequency [Hz]"; fontsize=14)
ylabel("Power Density "; fontsize=14)
title(string(round.(data["locations"][i,j]./RE, digits=1))*"Re"; fontsize=14)

loglog(freq[10:end], a.*freq[10:end].^k, label="k = $(round(k, digits=1))")

legend()


text:
module Neurodevo

using Distances
using Random

include("config.jl")
include("controller.jl")
include("controllers/random.jl")
include("controllers/const.jl")
include("controllers/static.jl")
include("controllers/snn.jl")
include("conn.jl")
include("cell.jl")
include("model.jl")
include("step.jl")
include("inits.jl")

end


text:
x, y, z = activate_github("arnavs/InstantiationTest")
@test z == Dict{AbstractString,Any}("Expectations"=>"2fe49d83-0758-5602-8f54-1f90ad0d522b","Distances"=>"b4f34e82-e78d-54a5-968a-f98e89d6e8f7","Suppressor"=>"fd094767-a336-5f1f-9728-57cf17d0bbfb","StatsFuns"=>"4c63d2b9-4356-54db-8cca-17b64c39e42c")