text:
using PersistenceDiagrams
using Test

@testset "Constructors" begin
    pi_1 = PersistenceImage((0, 2), (0, 3); size=(10, 15))
    @test length(pi_1.ys) == 11
    @test length(pi_1.xs) == 16
    @test first(pi_1.ys) == 0
    @test last(pi_1.ys) == 2
    @test first(pi_1.xs) == 0
    @test last(pi_1.xs) == 3
    @test pi_1.weight == PersistenceDiagrams.DefaultWeightingFunction(2.0)
    @test pi_1.distribution == PersistenceDiagrams.Binormal(0.375)

    pi_2 = PersistenceImage((1, 3), (2, 4); sigma=2.0, slope_end=0.5)
    @test length(pi_2.ys) == 6
    @test length(pi_2.xs) == 6
    @test first(pi_2.ys) == 1
    @test last(pi_2.ys) == 3
    @test first(pi_2.xs) == 2
    @test last(pi_2.xs) == 4
    @test pi_2.weight == PersistenceDiagrams.DefaultWeightingFunction(1.5)
    @test pi_2.distribution == PersistenceDiagrams.Binormal(2.0)

    @test PersistenceImage((0, 8), (0, 7); distribution=*).distribution == (*)
    @test PersistenceImage((0, 8), (0, 7); weight=*).weight == (*)

    @test_throws ArgumentError PersistenceImage((1, 2), (3, 4); sigma=4, distribution=*)
    @test_throws ArgumentError PersistenceImage((1, 2), (3, 4); slope_end=1, weight=*)
    @test_throws ArgumentError PersistenceImage((1, 2), (3, 4); slope_end=5)
    @test_throws ArgumentError PersistenceImage((1, 2), (3, 4); slope_end=0)
    @test_throws ArgumentError PersistenceImage((1, 2), (3, 4); sigma=0)

    diagram_1 = PersistenceDiagram([(1, 2), (1, 7), (3, 4), (2, 3)])
    diagram_2 = PersistenceDiagram([(0, 2), (-1, Inf), (4, 5)])

    pi_3 = PersistenceImage([diagram_1, diagram_2]; margin=0, zero_start=false)
    @test pi_3.ys == range(1, 6; length=6)
    @test pi_3.xs == range(0, 4; length=6)
    @test pi_3.weight == PersistenceDiagrams.DefaultWeightingFunction(6.0)

    pi_4 = PersistenceImage([diagram_1, diagram_2]; margin=0)
    @test pi_4.ys == range(0, 6; length=6)
    @test pi_4.xs == range(0, 4; length=6)
    @test pi_4.weight == PersistenceDiagrams.DefaultWeightingFunction(6.0)

    pi_5 = PersistenceImage([diagram_1, diagram_2]; margin=0.5)
    @test pi_5.ys == range(0, 9; length=6)
    @test pi_5.xs == range(-2, 6; length=6)
    @test pi_5.weight == PersistenceDiagrams.DefaultWeightingFunction(9.0)

    @test_throws ArgumentError PersistenceImage([diagram_1, diagram_2]; margin=-1)

    @test sprint(show, pi_1) == "10×15 PersistenceImage"
    @test sprint(show, pi_2) == "5×5 PersistenceImage"
    @test sprint(show, pi_3) == "5×5 PersistenceImage"

    @test sprint((io, x) -> show(io, MIME"text/plain"(), x), pi_1) ==
        "10×15 PersistenceImage(\n" *
          "  distribution = PersistenceDiagrams.Binormal(0.375),\n" *
          "  weight = PersistenceDiagrams.DefaultWeightingFunction(2.0),\n)"
end

@testset "Transform" begin
    diagram_1 = PersistenceDiagram([(1, 2), (1, 6), (3, 4), (2, 3)])
    diagram_2 = PersistenceDiagram([(0, 2), (0, Inf), (5, 5)])

    image_1 = PersistenceImage((0, 2), (0, 3); size=20, weight=*)
    image_2 = PersistenceImage((0, 2), (0, 3); size=(10, 15), slope_end=0.5)
    image_3 = PersistenceImage([diagram_1, diagram_2])

    @test size(image_1(diagram_1)) == (20, 20)
    @test size(image_2(diagram_1)) == (10, 15)
    @test size(image_3(diagram_1)) == (5, 5)
    @test size(image_1(diagram_2)) == (20, 20)
    @test size(image_2(diagram_2)) == (10, 15)
    @test size(image_3(diagram_2)) == (5, 5)

    diagram_3 = PersistenceDiagram([(1, 1), (2, 2), (3, 3)])
    @test all(iszero, image_1(diagram_3))
    @test all(iszero, image_2(diagram_3))
    @test all(iszero, image_3(diagram_3))
end


text:
using Clang.Generators
using SCIP_jll

cd(@__DIR__)

const (HEADER_BASE, is_default) = if haskey(ENV, "SCIPDIR")
    (joinpath(ENV["SCIPDIR"], "src"), false)
else
    (joinpath(SCIP_jll.artifact_dir, "include"), true)
end

const SCIP_TYPES_H = filter(readdir(joinpath(HEADER_BASE, "scip"), join=true)) do f
    occursin("type", f) && endswith(f, ".h")
end
push!(SCIP_TYPES_H, joinpath(HEADER_BASE, "scip/def.h"))
push!(SCIP_TYPES_H, joinpath(HEADER_BASE, "scip/nlpi.h"))
push!(SCIP_TYPES_H, joinpath(HEADER_BASE, "scip/scipdefplugins.h"))

# add the generated config.h file, either from JLL or build dir
if is_default
    push!(SCIP_TYPES_H, joinpath(HEADER_BASE, "scip/config.h"))
end

const SCIP_PUB_H = filter(readdir(joinpath(HEADER_BASE, "scip"), join=true)) do f
    occursin("pub", f) && endswith(f, ".h")
end

const SCIP_MEM_H = filter(readdir(joinpath(HEADER_BASE, "blockmemshell"), join=true)) do f
    endswith(f, ".h")
end

const SCIP_LPI_H = filter(readdir(joinpath(HEADER_BASE, "lpi"), join=true)) do f
    endswith(f, ".h")
end

headers = append!(
    SCIP_TYPES_H,
    SCIP_PUB_H,
    SCIP_MEM_H,
    SCIP_LPI_H,
)

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$HEADER_BASE")
if !is_default
    push!(args, "-I$SCIP_BUILD_DIR/scip")
end

ctx = create_context(headers, args, options)

build!(ctx)

ext: jl


text:
struct NotImplemetedError <: Exception
    msg :: String
end

Base.showerror(io::IO, e::NotImplemetedError) = print(io, e.msg)


text:
include("../functions/Functions.jl")
include("../evaluate/Evaluate.jl")
include("../opt/pgd.jl")
include("../data/data.jl")
include("../plot/Plotting.jl")

#using Plotting

X, y = getOldData()

x0 = zeros(size(X,2),1)

lambda = 0

iter = 1000
step = 0.01

f = f_logreg(X,y,x0,lambda)
g = g_logreg(X,y,x0,0)

#f = f_svm(X,y,x0,lambda)
#g = sub_g_svm()

prox_operator = SqrNormL2(lambda/2)

@time (theta,state,val_array,stop_array) = pgd(X,y,x0,f,g,1e-8,lambda,iter,step,size(X,1),prox_operator)

#plotArray(val_array, "sgd1.svg")

grad = g(theta)
grad_norm = norm(grad, Inf)
@printf("Inf-Norm of gradient is %f.\n",grad_norm)
println(evaluate(X,y,theta,"logreg",false))

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
const DEFAULT_JACOBIAN_ϵ = 8sqrt(eps())

abstract type JacobianParameters{T} end

compute_jacobian!(j::AbstractMatrix, x::AbstractVector, jacobian::JacobianParameters) = jacobian(j,x)

function check_jacobian(J::AbstractMatrix)
    println("Condition Number of Jacobian: ", cond(J))
    println("Determinant of Jacobian:      ", det(J))
    println("minimum(|Jacobian|):          ", minimum(abs.(J)))
    println("maximum(|Jacobian|):          ", maximum(abs.(J)))
    println()
end

function print_jacobian(J::AbstractMatrix)
    display(J)
    println()
end

struct JacobianParametersUser{T, JT} <: JacobianParameters{T}
    J!::JT
end

JacobianParametersUser(J!::JT, ::AbstractArray{T}) where {T,JT} = JacobianParametersUser{T,JT}(J!)

function (jac::JacobianParametersUser{T})(J::AbstractMatrix{T}, x::AbstractVector{T}) where {T}
    jac.J!(J, x)
end

struct JacobianParametersAD{T, FT, JT <: ForwardDiff.JacobianConfig, YT <: AbstractVector{T}} <: JacobianParameters{T}
    F!::FT
    Jconfig::JT
    ty::YT

    function JacobianParametersAD{T}(F!::FT, Jconfig::JT, y::YT) where {T, FT, JT, YT}
        new{T, FT, JT, YT}(F!, Jconfig, y)
    end

end

function JacobianParametersAD(F!::FT, x::AbstractVector{T}, y::AbstractVector{T}) where {T, FT}
    Jconfig = ForwardDiff.JacobianConfig(F!, y, x)
    JacobianParametersAD{T}(F!, Jconfig, zero(y))
end

function JacobianParametersAD{T}(F!::FT, nx::Int, ny::Int) where {T, FT}
    tx = zeros(T, nx)
    ty = zeros(T, ny)
    JacobianParametersAD(F!, tx, ty)
end

JacobianParametersAD{T}(F!, n) where {T} = JacobianParametersAD{T}(F!, n, n)

function (jac::JacobianParametersAD{T})(J::AbstractMatrix{T}, x::AbstractVector{T}) where {T}
    ForwardDiff.jacobian!(J, jac.F!, jac.ty, x, jac.Jconfig)
end

function compute_jacobian_ad!(J::AbstractMatrix{T}, x::AbstractVector{T}, F!::FT) where {T, FT}
    jac = JacobianParametersAD{T}(F!, length(x))
    jac(J,x)
end

struct JacobianParametersFD{T, FT} <: JacobianParameters{T}
    ϵ::T
    F!::FT
    f1::Vector{T}
    f2::Vector{T}
    e::Vector{T}
    tx::Vector{T}
end

function JacobianParametersFD{T}(F!::FT, nx::Int, ny::Int; ϵ=DEFAULT_JACOBIAN_ϵ) where {T, FT}
    f1 = zeros(T, ny)
    f2 = zeros(T, ny)
    e  = zeros(T, nx)
    tx = zeros(T, nx)
    JacobianParametersFD{T,FT}(ϵ, F!, f1, f2, e, tx)
end

JacobianParametersFD{T}(F!, n; kwargs...) where {T} = JacobianParametersFD{T}(F!, n, n; kwargs...)

function (jac::JacobianParametersFD{T})(J::AbstractMatrix{T}, x::AbstractVector{T}) where {T}
    local ϵⱼ::T

    for j in eachindex(x)
        ϵⱼ = jac.ϵ * x[j] + jac.ϵ
        fill!(jac.e, 0)
        jac.e[j] = 1
        jac.tx .= x .- ϵⱼ .* jac.e
        jac.F!(jac.f1, jac.tx)
        jac.tx .= x .+ ϵⱼ .* jac.e
        jac.F!(jac.f2, jac.tx)
        for i in eachindex(x)
            J[i,j] = (jac.f2[i] - jac.f1[i]) / (2ϵⱼ)
        end
    end
end

function compute_jacobian_fd!(J::AbstractMatrix{T}, x::AbstractVector{T}, F!::FT; kwargs...) where {T, FT}
    jac = JacobianParametersFD{T}(F!, length(x); kwargs...)
    jac(J,x)
end

function JacobianParameters{T}(ForJ, nx::Int, ny::Int; mode = :autodiff, diff_type = :forward, kwargs...) where {T}
    if mode == :autodiff
        if diff_type == :forward
            Jparams = JacobianParametersAD{T}(ForJ, nx, ny)
        else
            Jparams = JacobianParametersFD{T}(ForJ, nx, ny; kwargs...)
        end
    else
        Jparams = JacobianParametersUser{T, typeof(ForJ)}(ForJ)
    end
    return Jparams
end

JacobianParameters{T}(ForJ, n::Int; kwargs...) where {T} = JacobianParameters{T}(ForJ, n, n; kwargs...)

JacobianParameters{T}(J!, F!, nx, ny; kwargs...) where {T} = JacobianParameters{T}(J!, nx, ny; mode = :user, kwargs...)

JacobianParameters{T}(J!::Nothing, F!, nx, ny; kwargs...) where {T} = JacobianParameters{T}(F!, nx, ny;  mode = :autodiff, kwargs...)

JacobianParameters{T}(J!, F!, n; kwargs...) where {T} = JacobianParameters{T}(J!, F!, n, n; kwargs...)

function compute_jacobian!(j::AbstractMatrix{T}, x::AbstractVector{T}, ForJ; kwargs...) where {T}
    jacobian = JacobianParameters{T}(ForJ, size(j,1), size(j,2); kwargs...)
    jacobian(j,x)
end


text:
func(x) = cos.(x)
Inputs = 1.1
fp_simple   = fixed_point(func, Inputs; Algorithm = Simple)
fp_simple.Convergence_ < 1e-10
fp_anderson = fixed_point(func, Inputs; Algorithm = Anderson)
fp_anderson.Convergence_ < 1e-10
fp_aitken   = fixed_point(func, Inputs; Algorithm = Aitken)
fp_aitken.Convergence_ < 1e-10
fp_newton   = fixed_point(func, Inputs; Algorithm = Newton)
fp_newton.Convergence_ < 1e-10
fp_VEA      = fixed_point(func, Inputs; Algorithm = VEA)
fp_VEA.Convergence_ < 1e-10
fp_SEA      = fixed_point(func, Inputs; Algorithm = SEA)
fp_SEA.Convergence_ < 1e-10
fp_MPE      = fixed_point(func, Inputs; Algorithm = MPE)
fp_MPE.Convergence_ < 1e-10
fp_RRE      = fixed_point(func, Inputs; Algorithm = RRE)
fp_RRE.Convergence_ < 1e-10


text:
module GaugeAction_module
    import ..Abstractsmearing_module:CovNeuralnet
    import ..AbstractGaugefields_module:AbstractGaugefields,evaluate_gaugelinks!,add_U!,clear_U!,set_wing_U!,getvalue,
                            evaluate_gaugelinks_eachsite!

    import Wilsonloop:Wilsonline,make_staple
    using LinearAlgebra
    using InteractiveUtils

    struct GaugeAction_dataset{Dim}
        #β::Float64
        β::ComplexF64
        closedloops::Vector{Wilsonline{Dim}}
        staples::Vector{Vector{Wilsonline{Dim}}}
    end

    struct GaugeAction{Dim,T}
        hascovnet::Bool
        covneuralnet::Union{Nothing,CovNeuralnet{Dim}}
        dataset::Vector{GaugeAction_dataset{Dim}}
        _temp_U::Vector{T}
    end

    function GaugeAction_dataset(β,closedloops::Vector{Wilsonline{Dim}}) where Dim
        allstaples = Array{Vector{Wilsonline{Dim}},1}(undef,Dim)
        numloops = length(closedloops)
        for μ=1:Dim
            staples = Wilsonline{Dim}[]
            for i=1:numloops
                staples_i = make_staple(closedloops[i],μ)
                for j=1:length(staples_i)
                    push!(staples,staples_i[j])
                end
            end
            allstaples[μ] = staples
        end
        return GaugeAction_dataset{Dim}(β,closedloops,allstaples)
    end

    function get_temporary_gaugefields(S::GaugeAction)
        return S._temp_U
    end

    function calc_dSdUμ(S,μ,U::Vector{<: AbstractGaugefields{NC,Dim}}) where {Dim,NC}
        dSdUμ = similar(U[1])
        calc_dSdUμ!(dSdUμ,S,μ,U)
        return dSdUμ
    end

    function calc_dSdUμ!(dSdUμ,S,μ,U::Vector{<: AbstractGaugefields{NC,Dim}}) where {Dim,NC}
        temp1 = S._temp_U[1]
        temp2 = S._temp_U[2]
        temp3 = S._temp_U[3]
        numterm = length(S.dataset)

        clear_U!(dSdUμ)
        for i=1:numterm
            dataset = S.dataset[i]
            β =  dataset.β
            staples_μ =  dataset.staples[μ]
            evaluate_gaugelinks!(temp3,staples_μ,U,S._temp_U)

            #println("temp3 in dSdUμ! ",getvalue(temp3,1,1,1,1,1,1))
            add_U!(dSdUμ,β,temp3)
            #println("dSdUμ! ",getvalue(dSdUμ,1,1,1,1,1,1))
        end
        set_wing_U!(dSdUμ)

    end

    function evaluate_GaugeAction(S::GaugeAction,U::Vector{<: AbstractGaugefields{NC,Dim}}) where {Dim,NC}
        temp1 = S._temp_U[4]
        evaluate_GaugeAction_untraced!(temp1,S,U)
        value = tr(temp1)
        return value
    end

    function evaluate_GaugeAction_untraced(S::GaugeAction,U::Vector{<: AbstractGaugefields{NC,Dim}}) where {Dim,NC}
        uout = similar(U[1])
        clear_U!(uout)

        evaluate_GaugeAction_untraced!(uout,S,U)

        return uout
    end

    function evaluate_staple_eachindex!(mat_U,μ,S::GaugeAction,U::Vector{<: AbstractGaugefields{NC,Dim}},mat_temps,indices...) where {Dim,NC,T}
        temp3 = mat_temps[5]
        numterm = length(S.dataset)
        mat_U .= 0
        for i=1:numterm
            dataset = S.dataset[i]
            β = dataset.β
            staples_μ =  dataset.staples[μ]
            evaluate_gaugelinks_eachsite!(temp3,staples_μ,U,view(mat_temps,1:4),indices...)
            mat_U .+= β*temp3
        end
    end

    function evaluate_GaugeAction_untraced!(uout,S::GaugeAction,U::Vector{<: AbstractGaugefields{NC,Dim}}) where {Dim,NC,T}
        numterm = length(S.dataset)
        temp1 = S._temp_U[1]
        temp2 = S._temp_U[2]
        temp3 = S._temp_U[3]
        clear_U!(uout)

        for i=1:numterm
            dataset = S.dataset[i]
            β = dataset.β
            w = dataset.closedloops
            evaluate_gaugelinks!(temp3,w,U,[temp1,temp2])
            add_U!(uout,β,temp3)
        end
        set_wing_U!(uout)

        return
    end

    function GaugeAction(U::Vector{<: AbstractGaugefields{NC,Dim}};hascovnet = false) where {NC,Dim}
        if hascovnet
            covneuralnet = CovNeuralnet(Dim=Dim)
        else
            covneuralnet = nothing
        end
        dataset = GaugeAction_dataset{Dim}[]
        num = 4
        _temp_U = Array{eltype(U)}(undef,num)
        for i=1:num
            _temp_U[i] = similar(U[1])
        end

        return GaugeAction{Dim,eltype(U)}(hascovnet,covneuralnet,dataset,_temp_U)
    end

    function Base.push!(S::GaugeAction{Dim,T1},β::T,closedloops::Vector{Wilsonline{Dim}}) where {Dim,T <: Number,T1}
        dataset = GaugeAction_dataset(β,closedloops)
        push!(S.dataset,dataset)
    end

    function Base.show(s::GaugeAction{Dim,T}) where {Dim,T}
        println("----------------------------------------------")
        println("Structure of the actions for Gaugefields")
        println("num. of terms: ", length(s.dataset))
        for i=1:length(s.dataset)
            if i==1
                string = "st"
            elseif i==2
                string = "nd"
            elseif i==3
                string = "rd"
            else
                string = "th"
            end
            println("-------------------------------")
            println("      $i-$string term: ")
            println("          coefficient: ",s.dataset[i].β)
            println("      -------------------------")
            show(s.dataset[i].closedloops)
            println("      -------------------------")
        end
        println("----------------------------------------------")
    end

end


text:
using CSV, DataFrames, wgregseq, CairoMakie

##

# Open Gene information
dir = @__DIR__
home_dir = joinpath(split(dir, '/')[1:end-2]...)
# Import gene list to generate sequences for
gene_table = CSV.read("/$home_dir/data/100_genes.csv", DataFrames.DataFrame)
##
# Import promoter data
df_DB_prom = CSV.read(
    "/$home_dir/data/regulonDB/promoter.txt",
    DataFrames.DataFrame,
    comment="#",
    header=[
        "PROMOTER_ID",
        "PROMOTER_NAME",
        "PROMOTER_STRAND",
        "POS_1",
        "SIGMA_FACTOR",
        "BASAL_TRANS_VAL",
        "EQUILIBRIUM_CONST",
        "KINETIC_CONST",
        "STRENGTH_SEQ",
        "PROMOTER_SEQUENCE",
        "KEY_ID_ORG",
        "PROMOTER_NOTE",
        "PROMOTER_INTERNAL_COMMENT",
    ]
)

# Remove import artifacts
df_DB_prom = df_DB_prom[.~ ismissing.(df_DB_prom.PROMOTER_NAME), :]

# Change direction to be compatible with other database
direction_dict = Dict("forward" => "+", "reverse" => "-")
df_DB_prom.PROMOTER_STRAND = [direction_dict[x] for x in df_DB_prom.PROMOTER_STRAND]

# Fix type
df_DB_prom.PROMOTER_NAME = convert(Vector{String}, df_DB_prom.PROMOTER_NAME)
df_DB_prom

##
df_DB_genes = CSV.read(
    "/$home_dir/data/regulonDB/gene.txt",
    DataFrames.DataFrame,
    comment="#",
    header=[
        "GENE_ID",
        "GENE_NAME",
        "GENE_POSLEFT",
        "GENE_POSRIGHT",
        "GENE_STRAND",
        "GENE_SEQUENCE",
        "GC_CONTENT",
        "CRI_SCORE",
        "GENE_NOTE",
        "GENE_INTERNAL_COMMENT",
        "KEY_ID_ORG",
        "GENE_TYPE",
    ]
)

df_DB_genes = df_DB_genes[.~ ismissing.(df_DB_genes.GENE_NAME), :]
gene_pos = []

# Gene Strand information is missing for some genes, turns out those are split by transposons
df_trans = df_DB_genes[ismissing.(df_DB_genes.GENE_STRAND), :]
df_DB_genes = df_DB_genes[.~ ismissing.(df_DB_genes.GENE_STRAND), :]

for i in 1:nrow(df_DB_genes)
    if df_DB_genes[i, "GENE_STRAND"] == "forward"
        push!(gene_pos, df_DB_genes[i, "GENE_POSLEFT"])
    else
        push!(gene_pos, df_DB_genes[i, "GENE_POSRIGHT"])
    end
end
insertcols!(df_DB_genes, 13, :gene_position => gene_pos)
## Transcription units

df_DB_TU = CSV.read(
    "/$home_dir/data/regulonDB/transcription_unit.txt",
    DataFrames.DataFrame,
    comment="#",
    header=[
        "TRANSCRIPTION_UNIT_ID",
        "PROMOTER_ID",
        "TRANSCRIPTION_UNIT_NAME",
        "OPERON_ID",
        "KEY_ID_ORG",
        "TRANSCRIPTION_UNIT_NOTE",
        "TU_INTERNAL_COMMENT",
    ],
    delim='\t'
)

# Remove import artifacts
df_DB_TU = df_DB_TU[.~ ismissing.(df_DB_TU.TRANSCRIPTION_UNIT_ID), :]
# Remove transcription units without promoters
df_DB_TU = df_DB_TU[.~ ismissing.(df_DB_TU.PROMOTER_ID), :]

# Fix type
df_DB_TU.TRANSCRIPTION_UNIT_ID = convert(Vector{String}, df_DB_TU.TRANSCRIPTION_UNIT_ID)
df_DB_TU.PROMOTER_ID = convert(Vector{String}, df_DB_TU.PROMOTER_ID)
df_DB_TU

## Link fof tr

df_DB_link = CSV.read(
    "/$home_dir/data/regulonDB/tu_gene_link.txt",
    DataFrames.DataFrame,
    comment="#",
    header=[
        "TRANSCRIPTION_UNIT_ID",
        "GENE_ID",
    ],
    delim='\t'
)
df_DB_link
##
df_DB_tss = DataFrames.DataFrame()
for i in 1:nrow(df_DB_TU)
    # Get TU ID and promoter ID
    TU_ID, promoter_ID= df_DB_TU[i, ["TRANSCRIPTION_UNIT_ID", "PROMOTER_ID"]]

    # Find gene IDs of genes in TU
    gene_IDs = df_DB_link[map(x -> x == TU_ID, df_DB_link.TRANSCRIPTION_UNIT_ID), "GENE_ID"]

    # Get gene names and locations
    genes = convert(Vector{String}, df_DB_genes[map(x -> x in gene_IDs, df_DB_genes.GENE_ID), "GENE_NAME"])
    gene_positions = convert(Vector{Int64}, df_DB_genes[map(x -> x in gene_IDs, df_DB_genes.GENE_ID), "gene_position"])

    # get promoter data
    promoter_df = df_DB_prom[df_DB_prom.PROMOTER_ID .== promoter_ID, ["PROMOTER_NAME", "PROMOTER_STRAND", "POS_1",]]
    insertcols!(promoter_df, 2, :genes =>fill(genes, nrow(promoter_df)))
    insertcols!(promoter_df, 3, :gene_position =>fill(gene_positions, nrow(promoter_df)))
    rename!(
        promoter_df,
        Dict(
            "PROMOTER_NAME"=>"promoter",
            "PROMOTER_STRAND"=>"direction",
            "POS_1"=>"tss"
        )
    )
    append!(df_DB_tss, promoter_df)
end
df_DB_tss
CSV.write("/$home_dir/data/promoter_list_regulon_DB.csv", df_DB_tss)
##

df_DB_prom[df_DB_prom.PROMOTER_NAME .== "TSS_1868", :]
filter(x -> "dicA" in x, df_DB_tss.genes)

_d = dropmissing(df_DB_TU, "TRANSCRIPTION_UNIT_NAME")
_d[_d.TRANSCRIPTION_UNIT_NAME .== "dicA", :]

df_DB_tss


text:
# https://discourse.julialang.org/t/union-splitting-vs-c/61772/20
module Lines

using BenchmarkTools
using Test
using SingleDispatchArrays

abstract type AbstractLine end

for T in [:Line1, :Line2, :Line3, :Line4, :Line5]
    @eval begin
        struct $T <:AbstractLine
            length::Float64
        end
        paint(l::$T) = l.length
    end
end

struct Picture{T<:AbstractLine}
       lines::Vector{T}
end

# Dynamical dispatch at runtime

function paint1(p)
  s = 0.
  for l in p.lines
    s += paint(l)
  end
  s
end

# Union splitting

function paint2(p)
  s = 0.
  for l in p.lines
    if l isa Line1
      s += paint(l)
    elseif l isa Line2
      s += paint(l)
    elseif l isa Line3
      s += paint(l)
    elseif l isa Line4
      s += paint(l)
    elseif l isa Line5
      s += paint(l)
    end
  end
  s
end

function paint3(a::SingleDispatchVector)
    # Will be lowered to `Core.Box` if s is not constant
    s = Ref(0.)
    foreach(a) do l
        s[] += paint(l)
    end
    s[]
end

function run(n)

  line_types = [ Line1, Line2, Line3, Line4, Line5 ]
  p = Picture([line_types[rand(1:5)](rand()) for i in 1:n])

  SingleDispatchArrays.add_subtype.(AbstractLine, line_types)
  a = SingleDispatchVector{AbstractLine}(p.lines)

  @test paint1(p) ≈ paint3(a)
  println("result = ", paint1(p))

  print(" with runtime dispatch: "); @btime paint1($p)
  print(" with splitting: "); @btime paint2($p)
  print(" with SingleDispatchArrays.jl: "); @btime paint3($a)
end

end

# running
Lines.run(1_000_000)


text:
using LinearAlgebraicRepresentation
LAR = LinearAlgebraicRepresentation

function triangulate2D(model::Lar.Model)::Array{Lar.Cells}
    return LAR.triangulate2D(convert(Lar.Points, model.G'), [model.T[1], model.T[2]])
end

function FV2EVs(model::Lar.Model)::Array{Lar.Cells}
    return LAR.FV2EVs(model.T[1], model.T[2])
end


text:
import Pkg

VLC = try; success(`vlc --version`); catch; false; end

MAKIE = try
            Pkg.add("GLFW")
            using GLFW
            true
            Pkg.rm("GLFW")
        catch
            false
        end

VLC   && @info "VLC backend is available"
MAKIE && @info "Makie backend is available"

open(joinpath(@__DIR__, "includes.jl"), "w") do file

    VLC && println(file, "include(\"vlc.jl\")")

    MAKIE && println(file, "include(\"makie.jl\")")

end


text:
function Curlcurl(lattice::Lattice{1},α::Array{ComplexF64},α_half::Array{ComplexF64},nnm,nnp,indices,interior,surface,domain_wall)

	dx = lattice.dx

	α_half⁻¹ = 1 ./α_half

	N = length(α)

	∂ₓα⁻¹∂ₓ_bulk = _bulk_laplacian(lattice,N,interior,domain_wall,surface,nnm,nnp,α_half⁻¹)
	∂ₓα⁻¹∂ₓ_domain_wall = _domain_wall_laplacian(lattice,N,interior,domain_wall,nnm,nnp,α_half⁻¹,indices)
	∂ₓα⁻¹∂ₓ_surface = _surface_laplacian(lattice,N,interior,surface,nnm,nnp,α_half⁻¹,indices)
	∂ₓα⁻¹∂ₓ = ∂ₓα⁻¹∂ₓ_bulk + ∂ₓα⁻¹∂ₓ_domain_wall + ∂ₓα⁻¹∂ₓ_surface

	∂ₓ_bulk = _bulk_derivative(lattice,N,interior,domain_wall,surface,nnm,nnp)
	∂ₓ_domain_wall = _domain_wall_derivative(lattice,N,interior,domain_wall,nnm,nnp,α_half⁻¹,indices)
	∂ₓ_surface = _surface_derivative(lattice,N,interior,surface,nnm,nnp,indices)
	∂ₓ = ∂ₓ_bulk + ∂ₓ_domain_wall + ∂ₓ_surface
	i∂ₓ = 1im*∂ₓ

	return Curlcurl{1}(∂ₓα⁻¹∂ₓ,i∂ₓ,spdiagm(0=>α))
end

function (cc::Curlcurl{1})(ky::Real,kz::Real)

	ky² = ky^2
	kz² = kz^2
	kykz = ky*kz

	cc0 = kron(sparse([2,3],[2,3],[-1,-1],3,3),cc.cc0)
	cc1 = kron(sparse([2,3,1,1],[1,1,2,2],[ky,kz,ky,kz],3,3),cc.cc1)
	cc2 = kron(sparse([1,2,3,2,3],[1,2,2,3,3],[ky²+kz²,kz²,kykz,kykz,ky²],3,3),cc.cc2)

	return cc0+cc1+cc2
end

function _bulk_laplacian(lattice::Lattice{1},N,interior,domain_wall,surface,nnm,nnp,α_half⁻¹)
	rows = findall(interior .& .!(domain_wall .| surface))
	colsm = nnm[1][interior .& .!(domain_wall .| surface)]
	valsm = α_half⁻¹[colsm]
	colsp = nnp[1][interior .& .!(domain_wall .| surface)]
	valsp = α_half⁻¹[colsp]
	cols0 = rows
	vals0 = - α_half⁻¹[colsm] - α_half⁻¹[colsp]
	return sparse(vcat(rows,rows,rows),vcat(colsm,cols0,colsp),vcat(valsm,vals0,valsp)/lattice.dx^2,N,N)
end
function _domain_wall_laplacian(lattice::Lattice{1},N,interior,domain_wall,nnm,nnp,α_half⁻¹,indices)
	rows_dw = findall(domain_wall .& interior)
	cols_dwm = map(i->findfirst(isequal(indices[i]-CartesianIndex(1,)),indices),rows_dw)
	vals_dwm = α_half⁻¹[cols_dwm]
	cols_dwp = map(i->findfirst(isequal(indices[i]+CartesianIndex(1,)),indices),rows_dw)
	vals_dwp = α_half⁻¹[cols_dwp]
	cols_dw0 = rows_dw
	vals_dw0 = - α_half⁻¹[cols_dwm] - α_half⁻¹[cols_dwp]
	return sparse(vcat(rows_dw,rows_dw,rows_dw),vcat(cols_dwm,cols_dw0,cols_dwp),vcat(vals_dwm,vals_dw0,vals_dwp)/lattice.dx^2,N,N)
end
function _surface_laplacian(lattice,N,interior,surface,nnm,nnp,α_half⁻¹,indices)
	rows_sf = findall(surface .& interior)
	cols_sfm = map(i->findfirst(isequal(indices[i]-CartesianIndex(1,)),indices),rows_sf)
	rows_sfm = [rows_sf[findfirst(.!isnothing.(cols_sfm))]]
	cols_sfm = [cols_sfm[findfirst(.!isnothing.(cols_sfm))]]
	vals_sfm = α_half⁻¹[cols_sfm]
	cols_sfp = map(i->findfirst(isequal(indices[i]+CartesianIndex(1,)),indices),rows_sf)
	rows_sfp = [rows_sf[findfirst(.!isnothing.(cols_sfp))]]
	cols_sfp = [cols_sfp[findfirst(.!isnothing.(cols_sfp))]]
	vals_sfp = α_half⁻¹[cols_sfp]
	cols_sf0 = rows_sf
	vals_sf0 = - α_half⁻¹[[0,1] .+ rows_sf] - α_half⁻¹[[0,1] .+ rows_sf]
	return sparse(vcat(rows_sfm,rows_sf,rows_sfp),vcat(cols_sfm,cols_sf0,cols_sfp),vcat(vals_sfm,vals_sf0,vals_sfp)/lattice.dx^2,N,N)
end

function _bulk_derivative(lattice::Lattice{1},N,interior,domain_wall,surface,nnm,nnp)
	rows = findall(interior .& .!(domain_wall .| surface))
	colsm = nnm[1][interior .& .!(domain_wall .| surface)]
	valsm = -ones(size(colsm))
	colsp = nnp[1][interior .& .!(domain_wall .| surface)]
	valsp = ones(size(colsm))
	return sparse(vcat(rows,rows),vcat(colsm,colsp),vcat(valsm,valsp)/lattice.dx,N,N)
end
function _domain_wall_derivative(lattice::Lattice{1},N,interior,domain_wall,nnm,nnp,α_half⁻¹,indices)
	rows_dw = findall(domain_wall .& interior)
	cols_dwm = map(i->findfirst(isequal(indices[i]-CartesianIndex(1,)),indices),rows_dw)
	vals_dwm = -ones(size(cols_dwm))
	cols_dwp = map(i->findfirst(isequal(indices[i]+CartesianIndex(1,)),indices),rows_dw)
	vals_dwp = ones(size(cols_dwp))
	return sparse(vcat(rows_dw,rows_dw),vcat(cols_dwm,cols_dwp),vcat(vals_dwm,vals_dwp)/lattice.dx,N,N)
end
function _surface_derivative(lattice::Lattice{1},N,interior,surface,nnm,nnp,indices)
	rows_sf = findall(surface .& interior)
	cols_sfm = map(i->findfirst(isequal(indices[i]-CartesianIndex(1,)),indices),rows_sf)
	rows_sfm = [rows_sf[findfirst(.!isnothing.(cols_sfm))]]
	cols_sfm = [cols_sfm[findfirst(.!isnothing.(cols_sfm))]]
	vals_sfm = -ones(size(rows_sfm))
	cols_sfp = map(i->findfirst(isequal(indices[i]+CartesianIndex(1,)),indices),rows_sf)
	rows_sfp = [rows_sf[findfirst(.!isnothing.(cols_sfp))]]
	cols_sfp = [cols_sfp[findfirst(.!isnothing.(cols_sfp))]]
	vals_sfp = ones(size(rows_sfm))
	return sparse(vcat(rows_sfm,rows_sfp),vcat(cols_sfm,cols_sfp),vcat(vals_sfm,vals_sfp)/lattice.dx,N,N)
end


text:
using MPI
using ClimateMachine
using Logging
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
import ClimateMachine.DGMethods.FVReconstructions: FVConstant, FVLinear
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using LinearAlgebra
using Printf
using Test
using Dates
using ClimateMachine.GenericCallbacks:
    EveryXWallTimeSeconds, EveryXSimulationSteps
using ClimateMachine.VTK: writevtk, writepvtu

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

const output = parse(Bool, lowercase(get(ENV, "JULIA_CLIMA_OUTPUT", "false")))

include("advection_diffusion_model.jl")

struct Pseudo1D{u, v, ν} <: AdvectionDiffusionProblem end

# This test has two 2D equations : the first one is an advection equation with
# the Gaussian initial condition the second one is an advection diffusion
# equation with the Gaussian initial condition
function init_velocity_diffusion!(
    ::Pseudo1D{u, v, ν},
    aux::Vars,
    geom::LocalGeometry,
) where {u, v, ν}
    # Advection velocity of the flow is [u, v]
    uvec = SVector(u, v, 0)
    aux.advection.u = hcat(uvec, uvec)

    # Diffusion of the flow is νI (isentropic diffusivity)
    I3 = @SMatrix [1 0 0; 0 1 0; 0 0 1]
    aux.diffusion.D = hcat(0 * I3, ν * I3)

end

function initial_condition!(
    ::Pseudo1D{u, v, ν},
    state,
    aux,
    localgeo,
    t,
) where {u, v, ν}
    FT = typeof(u)
    # The computational domain is [-1.5 1.5]×[-1.5 1.5]
    Lx, Ly = 3, 3
    x, y, _ = localgeo.coord

    μx, μz, σ = 0, 0, Lx / 10
    ρ1 = exp.(-(((x - μx) / σ)^2 + ((y - μz) / σ)^2) / 2) / (σ * sqrt(2 * pi))
    ρ2 = exp.(-(((x - μx) / σ)^2 + ((y - μz) / σ)^2) / 2) / (σ * sqrt(2 * pi))

    state.ρ = (ρ1, ρ2)
end

inhomogenous_data!(::Val{0}, P::Pseudo1D, x...) = initial_condition!(P, x...)

function do_output(mpicomm, vtkdir, vtkstep, dgfvm, Q, Qe, model, testname)
    ## Name of the file that this MPI rank will write
    filename = @sprintf(
        "%s/%s_mpirank%04d_step%04d",
        vtkdir,
        testname,
        MPI.Comm_rank(mpicomm),
        vtkstep
    )

    statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
    exactnames = statenames .* "_exact"

    writevtk(filename, Q, dgfvm, statenames, Qe, exactnames)

    ## Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        ## Name of the pvtu file
        pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

        ## Name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
        end

        writepvtu(
            pvtuprefix,
            prefixes,
            (statenames..., exactnames...),
            eltype(Q),
        )

        @info "Done writing VTK: $pvtuprefix"
    end
end

function test_run(
    mpicomm,
    vtkdir,
    fvmethod,
    polynomialorders,
    level,
    ArrayType,
    FT,
)

    dim = 2
    Lx, Ly = FT(3), FT(3)
    u, v = FT(1), FT(1)
    ν = FT(1 // 100)

    # Grid/topology information
    base_num_elem = 4
    Ne = 2^(level - 1) * base_num_elem

    # Match number of points
    N_dg_point, N_fvm_point = Ne + 1, Ne * polynomialorders[1] + 1

    brickrange = (
        range(-Lx / 2; length = N_dg_point, stop = Lx / 2),
        range(-Ly / 2; length = N_fvm_point, stop = Ly / 2),
    )

    periodicity = ntuple(j -> true, dim)
    bc = ntuple(j -> (1, 2), dim)
    connectivity = dim == 3 ? :full : :face
    topl = StackedBrickTopology(
        mpicomm,
        brickrange;
        periodicity = periodicity,
        boundary = bc,
        connectivity = connectivity,
    )

    # One period
    timeend = Lx / u
    # dt ≤ CFL (Δx / Np²)/u   u = √2  CFL = 1/√2
    dt = Lx / (Ne * polynomialorders[1]^2)

    Nt = (Ne * polynomialorders[1]^2)

    outputtime = Ne * dt

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorders,
    )

    # Model being tested
    model = AdvectionDiffusion{dim}(Pseudo1D{u, v, ν}(), num_equations = 2)

    # Main DG discretization
    dgfvm = DGFVModel(
        model,
        grid,
        fvmethod,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        direction = EveryDirection(),
    )

    # Initialize all relevant state arrays and create solvers
    Q = init_ode_state(dgfvm, FT(0))

    solver = LSRK54CarpenterKennedy(dgfvm, Q; dt = dt, t0 = 0)

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            @info @sprintf(
                """Update
                simtime = %.16e
                runtime = %s
                norm(Q) = %.16e""",
                gettime(solver),
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )
        end
    end
    callbacks = (cbinfo,)
    if ~isnothing(vtkdir)
        # Create vtk dir
        mkpath(vtkdir)

        vtkstep = 0
        # Output initial step
        do_output(
            mpicomm,
            vtkdir,
            vtkstep,
            dgfvm,
            Q,
            Q,
            model,
            "advection_diffusion_periodic",
        )

        # Setup the output callback
        cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
            vtkstep += 1
            Qe = init_ode_state(dgfvm, gettime(solver))
            do_output(
                mpicomm,
                vtkdir,
                vtkstep,
                dgfvm,
                Q,
                Qe,
                model,
                "advection_diffusion_periodic",
            )
        end
        callbacks = (callbacks..., cbvtk)
    end

    numberofsteps = convert(Int64, cld(timeend, dt))
    dt = timeend / numberofsteps

    @info "time step" dt numberofsteps dt * numberofsteps timeend

    solve!(Q, solver; timeend = timeend, callbacks = callbacks)

    engf = norm(Q, dims = (1, 3))

    # Reference solution
    Q_ref = init_ode_state(dgfvm, FT(0))
    engfe = norm(Q_ref, dims = (1, 3))

    errf = norm(Q_ref .- Q, dims = (1, 3))

    metrics = [engf; engfe; errf]

    @info @sprintf """Finished
    Advection equation:
    norm(Q)                 = %.16e
    norm(Qe)                = %.16e
    norm(Q - Qe)            = %.16e
    Advection diffusion equation:
    norm(Q)                 = %.16e
    norm(Qe)                = %.16e
    norm(Q - Qe)            = %.16e
    """ metrics[:]...

    return errf
end

# Run this test problem
function main()

    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()
    mpicomm = MPI.COMM_WORLD

    # Dictionary keys: dim, level, polynomial order, FT, and direction
    expected_result = Dict()

    # Dim 2, degree 4 in the horizontal, FV order 1, refinement level, Float64, equation number
    expected_result[2, 4, FVConstant(), 1, Float64, 1] = 4.5196354911392578e-01
    expected_result[2, 4, FVConstant(), 1, Float64, 2] = 4.8622171647472179e-01
    expected_result[2, 4, FVConstant(), 2, Float64, 1] = 3.5051983693457450e-01
    expected_result[2, 4, FVConstant(), 2, Float64, 2] = 4.1168975732563706e-01
    expected_result[2, 4, FVConstant(), 3, Float64, 1] = 2.5130141068332995e-01
    expected_result[2, 4, FVConstant(), 3, Float64, 2] = 3.4635415661628755e-01
    expected_result[2, 4, FVConstant(), 4, Float64, 1] = 1.6320856055402777e-01
    expected_result[2, 4, FVConstant(), 4, Float64, 2] = 2.9647989774687805e-01
    expected_result[2, 4, FVConstant(), 5, Float64, 1] = 9.6641259241910610e-02
    expected_result[2, 4, FVConstant(), 5, Float64, 2] = 2.6372336617960646e-01

    expected_result[2, 4, FVConstant(), 1, Float32, 1] = 4.5196333527565002e-01
    expected_result[2, 4, FVConstant(), 1, Float32, 2] = 4.8622152209281921e-01
    expected_result[2, 4, FVConstant(), 2, Float32, 1] = 3.5051953792572021e-01
    expected_result[2, 4, FVConstant(), 2, Float32, 2] = 4.1168937087059021e-01
    expected_result[2, 4, FVConstant(), 3, Float32, 1] = 2.5130018591880798e-01
    expected_result[2, 4, FVConstant(), 3, Float32, 2] = 3.4635293483734131e-01
    expected_result[2, 4, FVConstant(), 4, Float32, 1] = 1.6320574283599854e-01
    expected_result[2, 4, FVConstant(), 4, Float32, 2] = 2.9647585749626160e-01
    expected_result[2, 4, FVConstant(), 5, Float32, 1] = 9.6632070839405060e-02
    expected_result[2, 4, FVConstant(), 5, Float32, 2] = 2.6370745897293091e-01

    # Dim 2, degree 4 in the horizontal, FV order 1, refinement level, Float64, equation number

    expected_result[2, 4, FVLinear(), 1, Float64, 1] = 2.2783152269907422e-01
    expected_result[2, 4, FVLinear(), 1, Float64, 2] = 3.1823753821941969e-01
    expected_result[2, 4, FVLinear(), 2, Float64, 1] = 9.2628269590376469e-02
    expected_result[2, 4, FVLinear(), 2, Float64, 2] = 2.4517823755458742e-01
    expected_result[2, 4, FVLinear(), 3, Float64, 1] = 3.1401247771425542e-02
    expected_result[2, 4, FVLinear(), 3, Float64, 2] = 2.2608977452273774e-01
    expected_result[2, 4, FVLinear(), 4, Float64, 1] = 9.8710350648653390e-03
    expected_result[2, 4, FVLinear(), 4, Float64, 2] = 2.2377315397364847e-01
    expected_result[2, 4, FVLinear(), 5, Float64, 1] = 2.9619222883137744e-03
    expected_result[2, 4, FVLinear(), 5, Float64, 2] = 2.2359698753564772e-01

    expected_result[2, 4, FVLinear(), 1, Float32, 1] = 2.2783152269907422e-01
    expected_result[2, 4, FVLinear(), 1, Float32, 2] = 3.1823753821941969e-01
    expected_result[2, 4, FVLinear(), 2, Float32, 1] = 9.2628269590376469e-02
    expected_result[2, 4, FVLinear(), 2, Float32, 2] = 2.4517823755458742e-01
    expected_result[2, 4, FVLinear(), 3, Float32, 1] = 3.1401247771425542e-02
    expected_result[2, 4, FVLinear(), 3, Float32, 2] = 2.2608977452273774e-01
    expected_result[2, 4, FVLinear(), 4, Float32, 1] = 9.8710350648653390e-03
    expected_result[2, 4, FVLinear(), 4, Float32, 2] = 2.2377315397364847e-01
    expected_result[2, 4, FVLinear(), 5, Float32, 1] = 2.9614968225359917e-03
    expected_result[2, 4, FVLinear(), 5, Float32, 2] = 2.2358775138854980e-01

    polynomialorders = (4, 0)
    numlevels = integration_testing ? 5 : 1

    # Dictionary keys: dim, level, polynomial order, FT, and direction
    @testset "$(@__FILE__)" begin
        for FT in (Float64, Float32)
            result = Dict()
            for fvmethod in (FVConstant(), FVLinear(), FVLinear{3}())
                @info @sprintf """Test parameters:
                ArrayType                   = %s
                FloatType                   = %s
                FV Reconstruction           = %s
                Dimension                   = %s
                Horizontal polynomial order = %s
                Vertical polynomial order   = %s
                """ ArrayType FT fvmethod 2 polynomialorders[1] polynomialorders[end]
                for level in 1:numlevels
                    result[level] = test_run(
                        mpicomm,
                        output ?
                        "vtk_advection_diffusion_2d" *
                        "_poly$(polynomialorders)" *
                        "_$(ArrayType)_$(FT)" *
                        "_level$(level)" :
                        nothing,
                        fvmethod,
                        polynomialorders,
                        level,
                        ArrayType,
                        FT,
                    )

                    fv_key = fvmethod isa FVLinear ? FVLinear() : fvmethod
                    @test result[level][1] ≈ FT(expected_result[
                        2,
                        polynomialorders[1],
                        fv_key,
                        level,
                        FT,
                        1,
                    ])
                    @test result[level][2] ≈ FT(expected_result[
                        2,
                        polynomialorders[1],
                        fv_key,
                        level,
                        FT,
                        2,
                    ])
                end

                @info begin
                    msg = "advection equation"
                    for l in 1:(numlevels - 1)
                        rate = log2(result[l][1]) - log2(result[l + 1][1])
                        msg *= @sprintf("\n  rate for level %d = %e", l, rate)
                    end
                    msg
                end
                # Not an analytic solution, convergence doesn't make sense
                # @info begin
                #     msg = "advection-diffusion equation"
                #     for l in 1:(numlevels - 1)
                #         rate = log2(result[l][2]) - log2(result[l + 1][2])
                #         msg *= @sprintf("\n  rate for level %d = %e", l, rate)
                #     end
                #     msg
                # end
            end
        end
    end
end

main()


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct MvDistDensity{D<:Distribution{Multivariate,Continuous}} <: AbstractDensity
    d::D
end

export MvDistDensity

function unsafe_density_logval(
    density::MvDistDensity,
    params::AbstractVector{<:Real},
    exec_context::ExecContext = ExecContext()
)
    Distributions.logpdf(density.d, params)
end

Base.parent(density::MvDistDensity) = density.d

nparams(density::MvDistDensity) = length(density.d)

# Assume that implementations of logpdf are thread-safe and remote-safe:
exec_capabilities(::typeof(unsafe_density_logval), density::MvDistDensity, params::AbstractVector{<:Real}) =
    ExecCapabilities(0, true, 0, true)

function unsafe_density_logval!(
    r::AbstractArray{<:Real},
    density::MvDistDensity,
    params::AbstractMatrix{<:Real},
    exec_context::ExecContext = ExecContext()
)
    # TODO: Parallel execution, depending on exec_context
    Distributions.logpdf!(r, density.d, params)
end

# Assume that implementations of logpdf! are thread-safe and remote-safe:
exec_capabilities(::typeof(unsafe_density_logval!), r::AbstractArray{<:Real}, density::MvDistDensity, params::AbstractMatrix{<:Real}) =
    ExecCapabilities(0, true, 0, true) # Change when implementation of density_logval! for MvDistDensity becomes multithreaded.

Distributions.sampler(density::MvDistDensity) = bat_sampler(parent(density))


text:
module HartGEMs

    import ...Chemostat_Rath2017
    import ..Human1
    import MAT
    import Chemostat
    import Chemostat.MetNets
    using ProjAssistant
    @gen_sub_proj

    include("const.jl")
    include("metnet.jl")

    function __init__()
        @create_proj_dirs
    end
end


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

# TYPE

mutable struct Tan <: TwoStageModel

    first_stage::OneStageModel
    second_stage::IV
    pscore::Vector{Float64}
    eweights::ProbabilityWeights{Float64, Float64, Vector{Float64}}

    Tan() = new()
end

#==========================================================================================#

# FIRST STAGE

function first_stage(::Type{Tan}, MM::Type{<:OneStageModel}, MD::Microdata; kwargs...)

    FSD                    = Microdata(MD)
    FSD.mapping[:response] = MD.mapping[:instrument]

    pop!(FSD.mapping, :treatment)
    pop!(FSD.mapping, :instrument)

    fit(MM, FSD; kwargs...)
end

#==========================================================================================#

# ESTIMATION

function fit(
        ::Type{Tan}, MM::Type{<:OneStageModel}, MD::Microdata; novar::Bool = false, kwargs...
    )

    m = first_stage(Tan, MM, MD; novar = novar)
    fit(Tan, m, MD; novar = novar, kwargs...)
end

function fit(
        ::Type{Tan},
        MM::OneStageModel,
        MD::Microdata;
        novar::Bool = false,
        trim::AbstractFloat = 0.0,
    )

    w = getweights(MD)
    z = getvector(MD, :instrument)
    p = mean(z, getweights(MD))
    π = predict(MM)
    v = [(1.0 - zi) * (1.0 - p) / (1.0 - πi) + zi * p / πi for (zi, πi) in zip(z, π)]

    v[((trim .> π) .| (1.0 - trim .< π))] .= 0.0

    SSD                   = Microdata(MD)
    SSD.mapping[:control] = asgn(MD.model, InterceptTerm{true}())

    obj              = Tan()
    obj.first_stage  = MM
    obj.second_stage = IV(SSD)
    obj.pscore       = π
    obj.eweights     = ProbabilityWeights(v)

    _fit!(second_stage(obj), reweight(w, obj.eweights))
    novar || _vcov!(obj, getcorr(obj), w)

    return obj
end

#==========================================================================================#

# SCORE (MOMENT CONDITIONS)

score(obj::Tan) = lmul!(Diagonal(obj.eweights), score(second_stage(obj)))

# EXPECTED JACOBIAN OF SCORE × NUMBER OF OBSERVATIONS

jacobian(obj::Tan, ::UnitWeights) = jacobian(second_stage(obj), obj.eweights)

function jacobian(obj::Tan, w::AbstractWeights)
    jacobian(second_stage(obj), reweight(w, obj.eweights))
end

# EXPECTED JACOBIAN OF SCORE W.R.T. FIRST-STAGE PARAMETERS × NUMBER OF OBSERVATIONS

function crossjacobian(obj::Tan, ::UnitWeights)

    z = getvector(obj, :instrument)
    p = mean(z)
    π = obj.pscore
    D = [(1.0 - zi) * (1.0 - p) / abs2(1.0 - πi) - zi * p / abs2(πi)
         for (zi, πi) in zip(z, π)]

    D[iszero.(obj.eweights)] .= 0.0

    g₁ = jacobexp(obj.first_stage)
    g₂ = score(second_stage(obj))

    return g₂' * lmul!(Diagonal(D), g₁)
end

function crossjacobian(obj::Tan, w::AbstractWeights)

    z = getvector(obj, :instrument)
    p = mean(z, w)
    π = obj.pscore
    D = [wi * ((1.0 - zi) * (1.0 - p) / abs2(1.0 - πi) - zi * p / abs2(πi))
         for (zi, πi, wi) in zip(z, π, w)]

    D[iszero.(obj.eweights)] .= 0.0

    g₁ = jacobexp(obj.first_stage)
    g₂ = score(second_stage(obj))

    return g₂' * lmul!(Diagonal(D), g₁)
end

#==========================================================================================#

# FITTED VALUES

predict(obj::Tan, MD::Microdata) = predict(second_stage(obj), MD)

# DERIVATIVE OF FITTED VALUES

jacobexp(obj::Tan) = jacobexp(second_stage(obj))

#==========================================================================================#

# UTILITIES

coefnames(obj::Tan)   = coefnames(second_stage(obj))
model_title(obj::Tan) = "Tan (2006)"

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
using CSDP

X_sol = Blockmatrix(
   [2 1
    1 2],
        [3 0 1
         0 2 0
         1 0 3],
      Diagonal([0,
                  0]))

Z_sol = Blockmatrix(
   [4 1
    1 2],
        [3 0 1
         0 2 0
         1 0 3],
      Diagonal([0,
                  0]))

y_origin = [5.5, 6.2]
y = Cvector{Cdouble}(CSDP.fptr(y_origin))

k = 2
n = 7

Xs = CSDP.blockmatrix(X_sol)
println("********** Alright **********")

ys = y.e
Zs = CSDP.blockmatrix(Z_sol)
@printf "********** Alright %s\n"  @__FILE__

CSDP.write_sol("prob.dat-s",
               n,
               k,
               X_sol,
               y,
               Z_sol)


text:
using MLDatasets, LightGraphs, Plots, LaTeXStrings, NGWP, MTSG
# Load local module
push!(LOAD_PATH, @__DIR__)
using pSGWT

example, label = MNIST.traindata(Float64, 28)
digit_img = example[4:26, 26:-1:5]
# heatmap(digit_img', ratio = 1, c=:viridis, frame = :none, xlim = [1, 22])

Nx, Ny = size(digit_img)
G = LightGraphs.grid([Nx, Ny]); N = nv(G);
W = Matrix(adjacency_matrix(G))
L = Matrix(laplacian_matrix(G))
Q = incidence_matrix(G; oriented = true)
𝛌, 𝚽 = eigen(L); 𝚽 = 𝚽 .* sign.(𝚽[1, :])';

D = natural_eigdist(𝚽, 𝛌, Q; distance = :DAG)

function grid_eigenvector_plot(l)
    heatmap(reshape(𝚽[:, l], (Nx, Ny))', c = :viridis, ratio = 1, frame = :none,
        xlim = [1, Nx], size = (500, 400), cbar = true)
end

function grid_NGWFvector_plot(l, x, D, 𝚽; σ = 0.2 * maximum(D))
    heatmap(reshape(ngwf_vector(D, l, x, 𝚽; σ = σ)', (Nx, Ny))',
        c = :viridis, ratio = 1, frame = :none, xlim = [1, Nx], size = (500, 400))
end

function plot_edge!(A, B; style = :solid, subplot = 1)
    plot!([A[1], B[1], NaN], [A[2], B[2], NaN], c = :red, legend = false,
        width = 10, style = style, subplot = subplot)
end

function plot_square!(Nx, Ny; subplot = 1)
    plot_edge!([-0.33, 0], [Nx + 1.24, 0]; subplot = subplot)
    plot_edge!([0, 0], [0, Ny + 1]; subplot = subplot)
    plot_edge!([-0.33, Ny + 1], [Nx + 1.24, Ny + 1]; subplot = subplot)
    plot_edge!([Nx + 0.92, 0], [Nx + 0.92, Ny + 1]; subplot = subplot)
end

function grid_vector_plot!(l, i, VECs)
    v = deepcopy(VECs[:, l])
    v ./= norm(v, 1)
    heatmap!(reshape(v, (Nx, Ny))', c = :viridis, ratio = 1, frame = :none,
        xlim = [-0.5, Nx+1.5], cbar = false, subplot = i)
end

function grid_vec_heatmap(VEC, Nx, Ny; l = 1)
    v = VEC[:, l]
    heatmap(reshape(v, (Nx, Ny))', c = :viridis, ratio = 1, frame = :none,
        xlim = [1, Nx], size = (500, 400), cbar = true)
end


text:
"""
    Dᶜ(N, Δ)

Returns a discrete 1D derivative operator for taking the derivative of a face-centered field with `N+1` grid points and `Δ` grid spacing and producing a cell-centered field with `N` grid points.
"""
function Dᶜ(N, Δ)
    D = zeros(N, N+1)
    for k in 1:N
        D[k, k]   = -1.0
        D[k, k+1] =  1.0
    end
    D = 1/Δ * D
    return D
end

"""
    Dᶠ(N, Δ)

Returns a discrete 1D derivative operator for taking the derivative of a cell-centered field with `N` grid points and `Δ` grid spacing and producing a face-centered field with `N+1` grid points.
"""
function Dᶠ(N, Δ)
    D = zeros(N+1, N)
    for k in 2:N
        D[k, k-1] = -1.0
        D[k, k]   =  1.0
    end
    D = 1/Δ * D
    return D
end

function cell_to_cell_derivative(D, data)
    face_data = D * data
    cell_data = 0.5 .* (@view(face_data[1:end-1]) + @view(face_data[2:end]))
    return cell_data
end


text:
push!(LOAD_PATH,"../src/")

using Documenter
using Raytracing

Documenter.makedocs(
	format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
	modules = [Raytracing],
	sitename = "Raytracer.jl",
	pages = [
			"Introduction" => "index.md",
			"Demo" => "demo.md",
			"The Render Function" => "render.md",
			"Interpreter for the Scene File" => "interpreter.md",
			"Base Structs" => "base_structs.md",
			"Reading and Writing PFM files" => "readingwritingpfm.md",
			"BRDFs and Pigments" => "brdfs_and_pigments.md",
			"Cameras" => "cameras.md",
			"Tone Mapping" => "tone_mapping.md",
			"Renderers" => "renderers.md",
			"Transformations" => "transformations.md",
			"Avaiable Shapes" => "shapes.md",
			"Range Tester Functions" => "range_testers.md",
			],
)

deploydocs(repo = "github.com/cosmofico97/Raytracing.git")


text:
```Write function that translates a text to Pig Latin and back. English is
translated to Pig Latin by taking the first letter of every word, moving it to
the end of the word and adding ‘ay’. “The quick brown fox” becomes “Hetay
uickqay rownbay oxfay”.```

st = "the quick brown fox"

function en_pig_lat(s::String)
    sp = split(s)
    for i in sp
        f = i[1]
        println(i[2:end], f, "ay")
    end

end

en_pig_lat(st)

st = "hetay uickqay rownbay oxfay"

function pig_lat_en(s::String)
    sp = split(s)
    for i in sp
        f = i[end-2]
        println(f,i[1:(end-3)])
    end
end

pig_lat_en(st)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
# slt_load.jl
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.6.0
# 09.12.2016
# Last Edit: 28.04.2021

"""# slt_load

Load a range of SLT files and convert the analog channels if conversion values are present in the corresponding CFG files

IMPORTANT NOTE: Gill sonic anemometers have a 30° offset from North, when a Gill instrument is listed on the "Sonic:" line of a configuration file then the 30° offset will be corrected.

`slt_load(dr,mindate,maxdate;average=false,verbose=true,angle_offset=0.0)`\n
* **dr**::String = SLT directory
* **mindate**::DateTime = Minimum of the time range to load
* **maxdate**::DateTime = Maximum of the time range to load
* **average**::Bool (optional) = Half hour average the data when TRUE, FALSE is default
* **verbose**::Bool (optional) = Display information of what is going on when TRUE, TRUE is default
* **angle_offset**::Real (optional) = Angle to add to the wind direction it will be added to the "Sonic Alignment" in the configuration file, 0.0 is default
"""
function slt_load(dr::String,mindate::DateTime,maxdate::DateTime;average::Bool=false,verbose::Bool=true,angle_offset::Real=0.0)
	###############
	## Constants ##
	###############
	GillSonics = 30.0 # Gill sonics have an offset of 30° towards West from North, see page 43 of the Gill R3 manual

	############
	## Checks ##
	############
	if !isdir(dr)
		error("First input must be a directory.")
	end

	#########################
	## List and Sort Files ##
	#########################
	verbose ? println("Listing SLT Files") : nothing
	(files,folders) = dirlist(dr,regex=r"\w\d{8}\.(slt|cfg)$")

	# Parse File Dates
	verbose ? println("Sorting " * string(length(files)) * " Files") : nothing
	fdates = Array{DateTime}(undef,length(files))
	i = 0
	for i=1:1:length(files)
		temp = basename(files[i])
		try
			yearstr = Meta.parse(temp[2:5])
			daystr = Meta.parse(temp[6:8])
			hourstr = Meta.parse(temp[9:10])
			minutestr = Meta.parse(temp[11:12])
			fdates[i] = DateTime(yearstr) + Dates.Day(daystr) - Dates.Day(1) + Dates.Hour(hourstr) + Dates.Minute(minutestr)
			catch
				println(files[i])
				error("Failed")
			end
	end

	# Save Config files before dates are filtered
	f = [occursin(r"\.cfg$",i) for i in files]
	cfgfiles = files[f]
	cfgdates = fdates[f]

	# Sort CFG Files by Date
	f = sortperm(cfgdates)
	cfgdates = cfgdates[f]
	cfgfiles = cfgfiles[f]

	# Save SLT files before dates are filtered
	f = [occursin(r"\.slt$",i) for i in files]
	files = files[f]
	fdates = fdates[f]

	# Remove files beyond the given time bounds
	f = findall(mindate .<= fdates .< maxdate)
	fdates = fdates[f]
	files = files[f]

	# Sort the Files By Date
	f = sortperm(fdates)
	fdates = fdates[f]
	files = files[f]

	# If no files are found return nothing
	if isempty(files)
		if average == false
			return DataFrame()
		else
			return DateTime[],[],[],[],[]
		end
	end

	###############################
	## Prepare CFG and SLT Info  ##
	###############################
	verbose ? println("Loading Settings for " * string(length(files)) * " files") : nothing

	# Remove Unnecessary Config Files
	f = findall(cfgdates .<= mindate)
	if isempty(f)
		error("No corresponding configuration files found, cannot continue.")
	end
	f = f[end]
	cfgdates = cfgdates[f:end]
	cfgfiles = cfgfiles[f:end]

	f = findall(cfgdates .< maxdate)
	cfgdates = cfgdates[f]
	cfgfiles = cfgfiles[f]

	# Load CFG info
	configs = slt_config(cfgfiles)
	configs = [configs;deepcopy(configs[1])] # Setup a fake final row with a time well beyond any real file name
	configs[end]["Time"] = DateTime(9999)

	# Check for column changes
	if length(unique(get.(configs,"Analog_Inputs",[-9999.0.*ones(Int,6)]))) > 1
		error("Analog inputs are not all the same, dimension mismatch.")
	end

	# Load SLT Info
	sltinfo = DataFrame()
	cfg = Array{Int8}(undef,length(files)) # List of which CFG file each
	offset = 1 # File processing list offset
	times = get.(configs,"Time",Dates.DateTime(0))
	while offset < length(files)
		cfgf = findall(times .<= fdates[offset])[end] # Find latest config
		config1 = configs[cfgf] # Current config and lower bound on file list
		config2 = configs[cfgf+1] # Upper bound on file list

		f = findall(get(config1,"Time",Dates.DateTime(1)) .<= fdates .< get(config2,"Time",Dates.DateTime(1)))
		tempfiles = files[f]
		analog_count = get(config1,"Analog_Count",0)
		sample_frequency = Meta.parse(get(config1,"Frequency","0"))
		slope = get(config1,"Slope",zeros(Float64,6))
		for i=1:1:length(tempfiles)
			cfg[offset] = cfgf
			temp = slt_header(tempfiles[i],analog_count,sample_frequency)
			temp[!,:Analog_Count] .= analog_count
			temp[!,:Sample_Frequency] .= sample_frequency
			temp[!,:Slope] = collect(Array[slope])
			if isempty(sltinfo)
				sltinfo = temp
			else
				sltinfo = [sltinfo;temp]
			end
			offset += 1
		end
	end

	# Rename DataFrame columns
	h = ["Time";"u";"v";"w";"sonic_temp";"speed_of_sound";"wind_direction";get(configs[1],"Analog_Names","")]
	h_unique = unique(h);
	if length(h) !== length(unique(h_unique))
		for i = 1:1:length(h_unique)
			n = 1
			for j = 1:1:length(h)
				if h[j] == h_unique[i]
					if n > 1
						h[j] = h_unique[i] * "_" * string(n)
					end
					n += 1
				end
			end
		end
	end

	# Preallocate Final Arrays
	verbose ? println("Preallocating Final Arrays (" * string(Int(sum(sltinfo.Line_Count))) * "," * string(4+get(configs[1],"Analog_Count",0)) * ")") : nothing
	l = Int(sum(sltinfo.Line_Count))
	col_types = fill!(Array{DataType}(undef,length(h)),Float64)
	col_types[1] = DateTime
	D = DataFrame([Array{i}(undef,l) for i in col_types],Symbol[Symbol(i) for i in h])
	D[:,2:end] .= NaN
	D[:,1] .= DateTime(0)

	###################
	## Load the Data ##
	###################
	verbose ? println("Loading Data (" * string(length(files)) * " files)") : nothing
	instrument_offset = 0.0 # If instrument in use by its design has an offset
	offset = 0
	u = 0.0
	v = 0.0
	filefolder = ""
	for i=1:1:length(files)
		if splitdir(files[i])[1] != filefolder
			filefolder = splitdir(files[i])[1]
			println("  " * filefolder)
		end
		println("\t" * string(i) * ": " * basename(files[i]))

		# Instrument Angle Offset
		if occursin("Gill",get.(configs,"Sonic","")[cfg[i]])
			instrument_offset = GillSonics
		else
			instrument_offset = 0.0
		end

		fid = open(files[i],"r")
		try
			seek(fid,sltinfo.Start_Pos[i])

			for j=1:1:Int(sltinfo.Line_Count[i])
				ms = Int(floor((j-1)*(1/sltinfo.Sample_Frequency[i])*1000)) # milliseconds
				D[j+offset,1] = sltinfo.T0[i] + Dates.Millisecond(ms) # Years, Months, Days, Hours, Minutes, Seconds, Milliseconds
				u = Float64(read!(fid,Array{Int16}(undef,1))[1])/100 # u
				v = Float64(read!(fid,Array{Int16}(undef,1))[1])/100 # v
				D[j+offset,2] = u
				D[j+offset,3] = v
				D[j+offset,4] = Float64(read!(fid,Array{Int16}(undef,1))[1])/100 # w
				temp = Float64(read!(fid,Array{Int16}(undef,1))[1])
				D[j+offset,5] = ((temp/50)^2)/403 - 273.16 # Tc
				D[j+offset,6] = temp/50 # c

				# Wind Direction
				if u == v == 0
					D[j+offset,7] = 0 # Otherwise u = v = 0 → NaN which messes up other calculations
				else
					if v < 0
						D[j+offset,7] = acosd(u/sqrt(u^2 + v^2)) - get.(configs,"Sonic_Alignment",[0])[cfg[i]] + instrument_offset + angle_offset # Wind Direction, A·B = |A||B|cos(Θ) → Θ = acos(Au/|A|) if B = [0 1 0]m/s (positive N wind)
					else
						D[j+offset,7] = 360 - acosd(u/sqrt(u^2 + v^2)) - get.(configs,"Sonic_Alignment",[0])[cfg[i]] + instrument_offset + angle_offset # Wind Direction, A·B = |A||B|cos(Θ) → Θ = acos(Au/|A|) if B = [0 1 0]m/s (positive N wind)
					end
				end

				for k=1:1:sltinfo.Analog_Count[i]
					V = Float64(read!(fid,Array{Int16}(undef,1))[1]) # V1
					if sltinfo.Bit_Mask[i][k] & 2^(1-1) > 0 # If the first bit is high
						# If the bit mask is high, use the following formula to convert the binary value to mV
						V = (V[1] + 25000)/10 # Binary to mV
					end
					b = Float64(get.(configs,"Analog_Lower",[0])[cfg[i]][k]) # Y-intercept
					m = get.(configs,"Slope",[0])[cfg[i]][k] # Slope
					D[j+offset,k+7] = m*V + b # Convert to a value
				end
			end
			offset += Int(sltinfo.Line_Count[i])
		catch e
			println("Error loading file: " * files[i])
			println(e)
		end
		close(fid)
	end

	##############################
	##  Correct Wind Direction  ##
	##############################
	f = findall(D[!,7] .> 360.0)
	D[f,7] = D[f,7] .- 360.0

	f = findall(D[!,7] .< 0.0)
	D[f,7] = D[f,7] .+ 360.0

	####################
	##  Average Data  ##
	####################
	tmean = DateTime[]
	Dmean = DataFrame[]
	Dstd = DataFrame[]
	Dmin = DataFrame[]
	Dmax = DataFrame[]
	if average
		verbose ? println("\tAveraging SLT Data") : nothing

		targets = collect(minimum(D.Time) - Dates.Millisecond(minimum(D.Time)):Dates.Minute(30):maximum(D.Time))
		actual = findnewton(D.Time,targets)
		f = [actual;length(D.Time) + 1]

		# Preallocate
		tmean = fill!(Array{DateTime}(undef,length(f) - 1),DateTime(0))
		l = length(f) - 1
		Dcols = [:u,:v,:w,:sonic_temp,:speed_of_sound,:wind_direction,:Analog1,:Analog2,:Analog3,:Analog4,:Analog5,:Analog6]
		Dtypes = [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64]

		Dmean = DataFrame([Array{i}(undef,l) for i in Dtypes[1:6+get.(configs,"Analog_Count",[0])[1]]],
				Dcols[1:6+get.(configs,"Analog_Count",[0])[1]])
		Dstd = DataFrame([Array{i}(undef,l) for i in Dtypes[1:6+get.(configs,"Analog_Count",[0])[1]]],
				Dcols[1:6+get.(configs,"Analog_Count",[0])[1]])
		Dmin = DataFrame([Array{i}(undef,l) for i in Dtypes[1:6+get.(configs,"Analog_Count",[0])[1]]],
				Dcols[1:6+get.(configs,"Analog_Count",[0])[1]])
		Dmax = DataFrame([Array{i}(undef,l) for i in Dtypes[1:6+get.(configs,"Analog_Count",[0])[1]]],
				Dcols[1:6+get.(configs,"Analog_Count",[0])[1]])
		for i = length(Dcols[1:4+get.(configs,"Analog_Count",[0])[1]])
			for j = l
				Dmean[j,i] = NaN
				Dstd[j,i] = NaN
				Dmin[j,i] = NaN
				Dmax[j,i] = NaN
			end
		end
		#names!(Dmean,names(D)[2:end])
		rename!(Dmean,names(D)[2:end])
		#names!(Dstd,names(D)[2:end])
		rename!(Dmean,names(D)[2:end])
		#names!(Dmin,names(D)[2:end])
		rename!(Dmean,names(D)[2:end])
		#names!(Dmax,names(D)[2:end])
		rename!(Dmean,names(D)[2:end])

		# Process the Data
		temp_Dmean = []
		temp_Dstd = []
		temp_Dmin = []
		temp_Dmax = []
		for j=1:1:length(f)-1
			try
				tmean[j] = D[f[j],1]
				temp_Dmean = mean(Array(D[f[j]:f[j+1]-1,2:end]),dims=1)
				temp_Dstd = std(Array(D[f[j]:f[j+1]-1,2:end]),dims=1)
				temp_Dmin = minimum(Array(D[f[j]:f[j+1]-1,2:end]),dims=1)
				temp_Dmax = maximum(Array(D[f[j]:f[j+1]-1,2:end]),dims=1)

				for k=1:1:length(temp_Dmean)
					Dmean[j,k] = temp_Dmean[k]
					Dstd[j,k] = temp_Dstd[k]
					Dmin[j,k] = temp_Dmin[k]
					Dmax[j,k] = temp_Dmax[k]
				end
			catch e
				println(e)
				println("\n\tsize(D) = " * string(size(D)))
				println("\tnames(D) = " * string(names(D)))
				println("\tsize(Dmean) = " * string(size(Dmean)))
				println("\tsize(temp_Dmean) = " * string(size(temp_Dmean)))
				println("\tminimum(t) = " * string(minimum(D.Time)))
				println("\tmaximum(t) = " * string(maximum(D.Time)))
				println("\tf[j] = " * string(f[j]))
				println("\tj = " * string(j))
			end
		end
	end

	verbose ? println("Complete") : nothing
	if average
		return tmean, Dmean, Dstd, Dmin, Dmax
	else
		return D
	end
end # slt_load(dr::String,mindate::DateTime,maxdate::DateTime)

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
using SafeTestsets

@safetestset "Day08" begin
    using AdventOfCode.TestUtils
    using AdventOfCode.Day08
    using AdventOfCode.Day08: checksum, decode

    @test checksum("123456789012", 3*2) == 1
    @test decode("0222112222120000", 2, 2) == [0 1; 1 0]

    btest("8.1") do
        part1() == 1320
    end

    btest("8.2") do
        part2() == ("OOO   OO  O   OO  O OOO  \n"*
                    "O  O O  O O   OO O  O  O \n"*
                    "O  O O     O O OO   O  O \n"*
                    "OOO  O      O  O O  OOO  \n"*
                    "O O  O  O   O  O O  O O  \n"*
                    "O  O  OO    O  O  O O  O \n")
    end
end


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

### Default Speed Limits in Kilometers Per Hour ###
const SPEED_ROADS_URBAN = Dict(
    1 => 95,    # Motorway
    2 => 72,    # Trunk
    3 => 48,    # Primary
    4 => 32,    # Secondary
    5 => 22,    # Tertiary
    6 => 12,    # Residential/Unclassified
    7 => 8,     # Service
    8 => 5)     # Living street

const SPEED_ROADS_RURAL = Dict(
    1 => 110,
    2 => 90,
    3 => 80,
    4 => 72,
    5 => 55,
    6 => 40,
    7 => 15,
    8 => 15)

### The following dictionaries are used for grouping ways    ###
### into limited, discrete classes for routing and plotting. ###

# Ordered by typical significance
const ROAD_CLASSES = Dict(
    "motorway" => 1,
    "trunk" => 2,
    "primary" => 3,
    "secondary" => 4,
    "tertiary" => 5,
    "unclassified" => 6,
    "residential" => 6,
    "service" => 7,
    "motorway_link" => 1,
    "trunk_link" => 2,
    "primary_link" => 3,
    "secondary_link" => 4,
    "tertiary_link" => 5,
    "living_street" => 8,
    "pedestrian" => 8,
    "road" => 6)

# Level 1: Cycleways, walking paths, and pedestrian streets
# Level 2: Sidewalks
# Level 3: Pedestrians typically allowed but unspecified
# Level 4: Agricultural or horse paths, etc.
const PED_CLASSES = Dict(
    "cycleway" => 1,
    "pedestrian" => 1,
    "living_street" => 1,
    "footway" => 1,
    "sidewalk" => 2,
    "sidewalk:yes" => 2,
    "sidewalk:both" => 2,
    "sidewalk:left" => 2,
    "sidewalk:right" => 2,
    "steps" => 2,
    "path" => 3,
    "residential" => 3,
    "service" => 3,
    "secondary" => 4,
    "tertiary" => 4,
    "primary" => 4,
    "track" => 4,
    "bridleway" => 4,
    "unclassified" => 4)

# Level 1: Bike paths
# Level 2: Separated bike lanes (tracks)
# Level 3: Bike lanes
# Level 4: Bikes typically allowed but not specified
const CYCLE_CLASSES = Dict(
    "cycleway" => 1,
    "cycleway:track" => 2,
    "cycleway:opposite_track" => 2,
    "bicycle:use_sidepath" => 2,
    "bicycle:designated" => 2,
    "cycleway:lane" => 3,
    "cycleway:opposite" => 3,
    "cycleway:opposite_lane" => 3,
    "cycleway:shared" => 3,
    "cycleway:share_busway" => 3,
    "cycleway:shared_lane" => 3,
    "bicycle:permissive" => 3,
    "bicycle:yes" => 3,
    "bicycle:dismount" => 4,
    "residential" => 4,
    "pedestrian" => 4,
    "living_street" => 4,
    "service" => 4,
    "unclassified" => 4)

const FEATURE_CLASSES = Dict(
    "aerialway" => 1,
    "aeroway" => 2,
    "amenity" => 3,
    "barrier" => 4,
    "boundary" => 5,
    "building" => 6,
    "craft" => 7,
    "emergency" => 8,
    "geological" => 9,
    "highway" => 10,
    "historic" => 11,
    "landuse" => 12,
    "leisure" => 13,
    "man_made" => 14,
    "military" => 15,
    "natural" => 16,
    "office" => 17,
    "place" => 18,
    "power" => 19,
    "public_transport" => 20,
    "railway" => 21,
    "route" => 22,
    "shop" => 23,
    "sport" => 24,
    "tourism" => 25,
    "waterway" => 26)

# Class 1: Residential/Accomodation
# Class 2: Commercial
# Class 3: Civic/Amenity
# Class 4: Other
# Class 5: Unclassified
const BUILDING_CLASSES = Dict(
    "accomodation" => 1,
    "apartments" => 1,
    "bungalow" => 1,
    "detached" => 1,
    "dormitory" => 1,
    "dwelling_house" => 1,
    "farm" => 1,
    "ger" => 1,
    "hotel" => 1,
    "house" => 1,
    "houseboat" => 1,
    "residential" => 1,
    "Residential" => 1,
    "semidetached_house" => 1,
    "static_caravan" => 1,
    "terrace" => 1,

    "commercial" => 2,
    "factory" => 2,
    "industrial" => 2,
    "manufacture" =>2,
    "office" => 2,
    "retail" => 2,
    "supermarket" => 2,
    "warehouse" => 2,

    "administrative" => 3,
    "cathedral" => 3,
    "chapel" => 3,
    "church" => 3,
    "civic" => 3,
    "kindergarten" => 3,
    "hospital" => 3,
    "mosque" => 3,
    "pavilion" => 3,
    "public" => 3,
    "school" => 3,
    "shrine" => 3,
    "stadium" => 3,
    "synagogue" => 3,
    "temple" => 3,
    "train_station" => 3,
    "transportation" => 3,
    "university" => 3,

    "barn" => 4,
    "bridge" => 4,
    "bunker" => 4,
    "cabin" => 4,
    "collapsed" => 4,
    "construction" => 4,
    "cowshed" => 4,
    "damaged" => 4,
    "farm_auxiliary" => 4,
    "garage" => 4,
    "garages" => 4,
    "greenhouse" => 4,
    "hangar" => 4,
    "hut" => 4,
    "kiosk" => 4,
    "roof" => 4,
    "semi" => 4,
    "service" => 4,
    "shed" => 4,
    "stable" => 4,
    "storage_tank" => 4,
    "sty" => 4,
    "tank" => 4,
    "transformer_tower" => 4,
    "ruins" => 4,

    "yes" => 5,
    "Yes" => 5)

### Standard map display "layers." ###
# const LAYER_STANDARD = Dict(
#     1 => Style(0x7BB6EF, 4), # Soft blue
#     2 => Style(0x66C266, 3), # Soft green
#     3 => Style(0xE68080, 3), # Soft red
#     4 => Style(0xFF9900, 3), # Soft orange
#     5 => Style(0xDADA47, 3), # Dark yellow
#     6 => Style(0x999999, 2), # Dark gray
#     7 => Style(0xE0E0E0, 2), # Light gray
#     8 => Style(0x999999, 1)) # Dark gray

# const LAYER_CYCLE = Dict(
#     1 => Style(0x006600, 3), # Green
#     2 => Style(0x5C85FF, 3), # Blue
#     3 => Style(0x5C85FF, 2), # Blue
#     4 => Style(0x999999, 2)) # Dark gray

# const LAYER_PED = Dict(
#     1 => Style(0x999999, 3), # Dark gray
#     2 => Style(0x999999, 3), # Dark gray
#     3 => Style(0x999999, 2), # Dark gray
#     4 => Style(0xE0E0E0, 2)) # Light gray

# const LAYER_FEATURES = Dict(
#     1 => Style(0x9966FF, 1.5, "."), # Lavender
#     2 => Style(0xFF0000, 1.5, "."), # Red
#     3 => Style(0x000000, 1.5, "."), # Black
#     4 => Style(0xFF66FF, 1.5, "."), # Pink
#     5 => Style(0x996633, 1.5, "."), # Brown
#     6 => Style(0xFF9900, 2.0, "."), # Orange
#     7 => Style(0xCC00CC, 1.5, ".")) # Brown

# const LAYER_BUILDINGS = Dict(
#     1 => Style(0xE1E1EB, 1, "-"), # Lighter gray
#     2 => Style(0xB8DBFF, 1, "-"), # Light blue
#     3 => Style(0xB5B5CE, 1, "-"), # Light gray
#     4 => Style(0xFFFF99, 1, "-"), # Pale yellow
#     5 => Style(0x006600, 1, "-")) # Green


text:
# This file was generated, do not modify it. # hide
using Gnuplot, ColorSchemes, Random
Random.seed!(123)
n = 15
x, y, color = rand(n), rand(n), rand(n)
cmaps = [:cool, :viridis, :plasma, :inferno, :thermal,
    :leonardo, :winter, :spring, :ice]
@gp "set multiplot layout 3,3; set key off" :-
for i in 1:9
    @gp :- i title = "$(cmaps[i]), pt $(i)" "set size square" :-
    @gp :- x y color "w p pt $(i) ps 3 lc palette" palette(cmaps[i]) :-
    @gp :- "set cbtics out nomirror" "set xtics 0,1" "set ytics 0,1" :-
    @gp :- "set cbtics 0.0,0.2,1.0" xrange= (-0.1,1.1) yrange= (-0.1,1.1) :-
end
@gp
save(term="pngcairo size 900,800", output="./code/plt_ex12.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
struct BatchClosure{F, A, B}
    f::F
end
function (b::BatchClosure{F,A,B})(p::Ptr{UInt}) where {F,A,B}
    (offset, args) = ThreadingUtilities.load(p, A, 2*sizeof(UInt))
    (offset, start) = ThreadingUtilities.load(p, UInt, offset)
    (offset, stop ) = ThreadingUtilities.load(p, UInt, offset)
    b.f(args, (start+one(UInt))%Int, stop%Int)
    B && free_local_threads!()
    nothing
end

@inline function batch_closure(f::F, args::A, ::Val{B}) where {F,A,B}
    bc = BatchClosure{F,A,B}(f)
    @cfunction($bc, Cvoid, (Ptr{UInt},))
end

@inline function setup_batch!(p::Ptr{UInt}, fptr::Ptr{Cvoid}, argtup, start::UInt, stop::UInt)
    offset = ThreadingUtilities.store!(p, fptr, sizeof(UInt))
    offset = ThreadingUtilities.store!(p, argtup, offset)
    offset = ThreadingUtilities.store!(p, start, offset)
    offset = ThreadingUtilities.store!(p, stop, offset)
    nothing
end
@inline function launch_batched_thread!(cfunc, tid, argtup, start, stop)
    fptr = Base.unsafe_convert(Ptr{Cvoid}, cfunc)
    ThreadingUtilities.launch(tid, fptr, argtup, start, stop) do p, fptr, argtup, start, stop
        setup_batch!(p, fptr, argtup, start, stop)
    end
end
_extract_params(::Type{T}) where {T<:Tuple} = T.parameters
_extract_params(::Type{NamedTuple{S,T}}) where {S,T<:Tuple} = T.parameters
push_tup!(x, ::Type{T}, t) where {T<:Tuple} = push!(x,t)
push_tup!(x, ::Type{NamedTuple{S,T}}, t) where {S,T<:Tuple} = push!(x, Expr(:call, Expr(:curly, :NamedTuple, S), t))

function add_var!(q, argtup, gcpres, ::Type{T}, argtupname, gcpresname, k) where {T}
    parg_k = Symbol(argtupname, :_, k)
    garg_k = Symbol(gcpresname, :_, k)
    if (T <: Tuple) || (T <: NamedTuple)
        push!(q.args, Expr(:(=), parg_k, Expr(:ref, argtupname, k)))
        t = Expr(:tuple)
        for (j,p) ∈ enumerate(_extract_params(T))
            add_var!(q, t, gcpres, p, parg_k, garg_k, j)
        end
        push_tup!(argtup.args, T, t)
    else
        push!(q.args, Expr(:(=), Expr(:tuple, parg_k, garg_k), Expr(:call, :object_and_preserve, Expr(:ref, argtupname, k))))
        push!(argtup.args, parg_k)
        push!(gcpres.args, garg_k)
    end
end

@generated function _batch_no_reserve(
    f!::F, threadmask, nthread, torelease, Nr, Nd, ulen, args::Vararg{Any,K}
    ) where {F,K}
    1+2
    q = quote
        $(Expr(:meta,:inline))
        threads = UnsignedIteratorEarlyStop(threadmask, nthread)
        Ndp = Nd + one(Nd)
    end
    block = quote
        start = zero(UInt)
        i = 0x00000000
        tid = 0x00000000
        tm = mask(threads)
        while true
            VectorizationBase.assume(tm ≠ zero(tm))
            tz = trailing_zeros(tm) % UInt32
            stop = start + ifelse(i < Nr, Ndp, Nd)
            i += 0x00000001
            tz += 0x00000001
            tid += tz
            tm >>>= tz
            launch_batched_thread!(cfunc, tid, argtup, start, stop)
            start = stop
            i == nthread && break
        end
        f!(map(dereference, argtup), (start+one(UInt)) % Int, ulen % Int)
        tm = mask(threads)
        tid = 0x00000000
        while true
            VectorizationBase.assume(tm ≠ zero(tm))
            tz = trailing_zeros(tm) % UInt32
            tz += 0x00000001
            tm >>>= tz
            tid += tz
            ThreadingUtilities.wait(tid)
            iszero(tm) && break
        end
        free_threads!(torelease)
        nothing
    end
    gcpr = Expr(:gc_preserve, block, :cfunc)
    argt = Expr(:tuple)
    for k ∈ 1:K
        add_var!(q, argt, gcpr, args[k], :args, :gcp, k)
    end
    push!(q.args, :(argtup = $argt), :(cfunc = batch_closure(f!, argtup, Val{false}())), gcpr)
    push!(q.args, nothing)
    q
end
@generated function _batch_reserve(
    f!::F, threadmask, nthread, unused_threads, torelease, Nr, Nd, ulen, args::Vararg{Any,K}
) where {F,K}
    q = quote
        nbatch = nthread + one(nthread)
        threads = UnsignedIteratorEarlyStop(threadmask, nthread)
        Ndp = Nd + one(Nd)
        nres_per = Base.udiv_int(unused_threads, nbatch)
        nres_rem = unused_threads - nres_per * nbatch
        nres_prr = nres_per + one(nres_per)
    end
    block = quote
        start = zero(UInt)
        i = zero(nres_rem)
        tid = 0x00000000
        tm = mask(threads)
        wait_mask = zero(worker_type())
        while true
            VectorizationBase.assume(tm ≠ zero(tm))
            tz = trailing_zeros(tm) % UInt32
            reserve = ifelse(i < nres_rem, nres_prr, nres_per)
            tz += 0x00000001
            stop = start + ifelse(i < Nr, Ndp, Nd)
            tid += tz
            tid_to_launch = tid
            wait_mask |= (one(wait_mask) << (tid - one(tid)))
            tm >>>= tz
            reserved_threads = zero(worker_type())
            for _ ∈ 1:reserve
                VectorizationBase.assume(tm ≠ zero(tm))
                tz = trailing_zeros(tm) % UInt32
                tz += 0x00000001
                tid += tz
                tm >>>= tz
                reserved_threads |= (one(reserve) << (tid - one(tid)))
            end
            reserve_threads!(tid_to_launch, reserved_threads)
            launch_batched_thread!(cfunc, tid_to_launch, argtup, start, stop)
            i += one(i)
            start = stop
            i == nthread && break
        end
        f!(map(dereference, argtup), (start+one(UInt)) % Int, ulen % Int)
        tid = 0x00000000
        while true
            VectorizationBase.assume(wait_mask ≠ zero(wait_mask))
            tz = (trailing_zeros(wait_mask) % UInt32) + 0x00000001
            wait_mask >>>= tz
            tid += tz
            ThreadingUtilities.wait(tid)
            iszero(wait_mask) && break
        end
        nothing
    end
    gcpr = Expr(:gc_preserve, block, :cfunc)
    argt = Expr(:tuple)
    for k ∈ 1:K
        add_var!(q, argt, gcpr, args[k], :args, :gcp, k)
    end
    push!(q.args, :(argtup = $argt), :(cfunc = batch_closure(f!, argtup, Val{true}())), gcpr)
    push!(q.args, nothing)
    q
end

@inline function batch(
    f!::F, (len, nbatches)::Tuple{Vararg{Integer,2}}, args::Vararg{Any,K}
) where {F,K}
    threads, torelease = request_threads(Base.Threads.threadid(), nbatches - one(nbatches))
    nthread = length(threads)
    ulen = len % UInt
    if nthread % Int32 ≤ zero(Int32)
        f!(args, one(Int), ulen % Int)
        return
    end
    nbatch = nthread + one(nthread)
    Nd = Base.udiv_int(ulen, nbatch % UInt) # reasonable for `ulen` to be ≥ 2^32
    Nr = ulen - Nd * nbatch

    _batch_no_reserve(f!, mask(threads), nthread, torelease, Nr, Nd, ulen, args...)
end
function batch(
    f!::F, (len, nbatches, reserve_per_worker)::Tuple{Vararg{Integer,3}}, args::Vararg{Any,K}
) where {F,K}
    ulen = len % UInt
    nbatches ≤ 1 && @goto NOTHREADS
    requested_threads = reserve_per_worker*nbatches
    threads, torelease = request_threads(Base.Threads.threadid(), requested_threads - one(nbatches))
    nthread = length(threads)
    iszero(nthread) && @goto NOTHREADS
    total_threads = nthread + one(nthread)
    nbatch = min(total_threads, nbatches % UInt32)

    Nd = Base.udiv_int(ulen, nbatch % UInt)
    Nr = ulen - Nd * nbatch

    unused_threads = total_threads - nbatch
    if iszero(unused_threads)
        _batch_no_reserve(f!, mask(threads), nthread, torelease, Nr, Nd, ulen, args...)
    else
        _batch_reserve(f!, mask(threads), nbatch-one(nbatch), unused_threads, torelease, Nr, Nd, ulen, args...)
    end
    return nothing
    @label NOTHREADS
    f!(args, one(Int), ulen%Int)
    return nothing
end


text:
using FractionalTimeDG
using PyPlot

α = 3/4
N = 1000
r = 4
M = 2r
version = 1

store = Store(α, r, M)
H = coef_H_uniform!(r, N, M, store, version)

ℓ = collect(1:N-1)

function max_antidiags(H)
    r = size(H[0], 1)
    N = length(H)
    γ = zeros(N-1, 2r-1)
    for m = 2:2r
        for ℓbar = 1:N-1
            if m ≤ r
                antidiag = [ H[ℓbar][i,m-i] for i = 1:m-1 ]
            else
                antidiag = [ H[ℓbar][i,m-i] for i = m-r:r ]
            end
#            println(antidiag)
            largest = maximum(abs.(antidiag))
            if largest ≥ eps(Float64)
                γ[ℓbar, m-1] = largest
            else
                γ[ℓbar, m-1] = NaN
            end
        end
    end
    return γ
end

γ = max_antidiags(H)

figure(2)
fmt = [":", "--", "-.", "-", ":", "--", "-."]
for m = 2:2r
    loglog(ℓ, γ[:,m-1], fmt[m-1])
end
legend([latexstring("i+j=$m") for m = 2:2r])
grid(true)
xlabel(L"$\bar\ell$", fontsize=12)
savefig("fig2.eps")


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

mutable struct VirtualMachineScaleSetInstanceView <: SwaggerModel
    virtualMachine::Any # spec type: Union{ Nothing, VirtualMachineScaleSetInstanceViewStatusesSummary } # spec name: virtualMachine
    extensions::Any # spec type: Union{ Nothing, Vector{VirtualMachineScaleSetVMExtensionsSummary} } # spec name: extensions
    statuses::Any # spec type: Union{ Nothing, Vector{InstanceViewStatus} } # spec name: statuses

    function VirtualMachineScaleSetInstanceView(;virtualMachine=nothing, extensions=nothing, statuses=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetInstanceView, Symbol("virtualMachine"), virtualMachine)
        setfield!(o, Symbol("virtualMachine"), virtualMachine)
        validate_property(VirtualMachineScaleSetInstanceView, Symbol("extensions"), extensions)
        setfield!(o, Symbol("extensions"), extensions)
        validate_property(VirtualMachineScaleSetInstanceView, Symbol("statuses"), statuses)
        setfield!(o, Symbol("statuses"), statuses)
        o
    end
end # type VirtualMachineScaleSetInstanceView

const _property_map_VirtualMachineScaleSetInstanceView = Dict{Symbol,Symbol}(Symbol("virtualMachine")=>Symbol("virtualMachine"), Symbol("extensions")=>Symbol("extensions"), Symbol("statuses")=>Symbol("statuses"))
const _property_types_VirtualMachineScaleSetInstanceView = Dict{Symbol,String}(Symbol("virtualMachine")=>"VirtualMachineScaleSetInstanceViewStatusesSummary", Symbol("extensions")=>"Vector{VirtualMachineScaleSetVMExtensionsSummary}", Symbol("statuses")=>"Vector{InstanceViewStatus}")
Base.propertynames(::Type{ VirtualMachineScaleSetInstanceView }) = collect(keys(_property_map_VirtualMachineScaleSetInstanceView))
Swagger.property_type(::Type{ VirtualMachineScaleSetInstanceView }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetInstanceView[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetInstanceView }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetInstanceView[property_name]

function check_required(o::VirtualMachineScaleSetInstanceView)
    true
end

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


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

module GeoSimulation

using Meshes
using GeoStatsBase
using Variography
using KrigingEstimators

using Distributions
using LinearAlgebra
using Statistics
using FFTW
using CpuId

import GeoStatsBase: preprocess, solve, solvesingle

include("lu.jl")
include("fft.jl")
include("seq.jl")
include("sgs.jl")
include("cookie.jl")

export
  # generic solvers
  SeqSim,

  # concrete solvers
  LUGS, FFTGS, SGS,

  # meta solvers
  CookieCutter

end


text: Dict( "hello" => "Hallo Latn de")


text:
"""
    BSplines

Module defining B-spline bases and B-spline related functions.
"""
module BSplines

export
    AbstractBSplineBasis,
    BasisFunction,
    BSplineBasis,
    BSplineOrder,
    boundaries,
    order,
    knots,
    basis,
    evaluate,
    evaluate!,
    evaluate_all,
    nonzero_in_segment,
    support

using ..DifferentialOps

"""
    AbstractBSplineBasis{k,T}

Abstract type defining a B-spline basis, or more generally, a functional basis
defined from B-splines.

The basis is represented by a B-spline order `k` and a knot element type `T`.

---

    (B::AbstractBSplineBasis)(
        x::Real, [op = Derivative(0)], [T = float(typeof(x))];
        [ileft = nothing],
    ) -> (i, bs)

Evaluates all basis functions which are non-zero at `x`.

This is a convenience alias for `evaluate_all`.
See [`evaluate_all`](@ref) for details on optional arguments and on the returned values.

# Examples

```jldoctest
julia> B = BSplineBasis(BSplineOrder(4), -1:0.1:1)
23-element BSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4  …  0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0]

julia> i, bs = B(0.42)
(18, (0.0013333333333333268, 0.28266666666666657, 0.6306666666666667, 0.08533333333333339))

julia> sum(bs)
1.0

julia> bs[1] - B[i](0.42)
0.0

julia> bs[2] - B[i - 1](0.42)
-5.551115123125783e-17

julia> B(0.44; ileft = i)
(18, (0.01066666666666666, 0.4146666666666667, 0.5386666666666665, 0.03599999999999999))

julia> B(0.42, Float32)
(18, (0.0013333336f0, 0.28266668f0, 0.6306667f0, 0.085333325f0))

julia> B(0.42, Derivative(1))
(18, (0.19999999999999937, 6.4, -3.3999999999999977, -3.200000000000001))
```
"""
abstract type AbstractBSplineBasis{k,T} end

@inline (B::AbstractBSplineBasis)(args...; kws...) =
    evaluate_all(B, args...; kws...)

"""
    BSplineOrder(k::Integer)

Specifies the B-spline order `k`.
"""
struct BSplineOrder{k} end

@inline BSplineOrder(k::Integer) = BSplineOrder{Int(k)}()

include("basis.jl")
include("basis_function.jl")
include("knots.jl")
include("evaluate_all.jl")

end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
using GaussianBasis.Libcint

export ∇FD_overlap, ∇FD_kinetic, ∇FD_nuclear, ∇FD_ERI_2e4c, ∇FD_sparseERI_2e4c
export ∇overlap, ∇kinetic, ∇nuclear
export ∇ERI_2e4c, ∇sparseERI_2e4c

include("Gradients/FiniteDifferences.jl")
include("Gradients/OneElectronGrad.jl")
include("Gradients/TwoElectronGrad.jl")


text:
"""
    cov2cor(mat::AbstractMatrix)

Converts a matrix (representing a covariance matrix) into a `Hermitian` correlation matrix and a vector of standard deviations.
### Inputs
* `cor` - A matrix.
### Returns
* A `Hermitian`.
* A `Vector` of standard deviations (not volatilities).
"""
function cov2cor(mat::AbstractMatrix)
    sdevs = sqrt.(diag(mat))
    cor = mat ./ (sdevs * transpose(sdevs))
    cor[diagind(cor)] .= 1
    return Hermitian(cor), sdevs
end

"""
    cov2cor_and_vol(mat::AbstractMatrix, duration_of_covariance_matrix::Dates.Period, duration_for_desired_vols::Dates.Period)
    cov2cor_and_vol(mat::AbstractMatrix, duration_of_covariance_matrix::Real, duration_for_desired_vols::Real)

Converts a matrix (representing a covariance matrix) into a `Hermitian` correlation matrix and a vector of volatilities.
### Inputs
* `cor` - A correlation matrix.
* `duration_of_covariance_matrix` - The duration of the covariance matrix. If these are input as reals they must have the same units.
* `duration_for_desired_vols` - The duration you want a volatility for. If these are input as reals they must have the same units.
### Returns
* A `Hermitian`.
* A `Vector` of volatilities.

    cov2cor_and_vol(mat::AbstractMatrix, duration_of_covariance_matrix_in_nanoseconds::Real, duration_for_desired_vols::Dates.Period)

### Inputs
* `cor` - A correlation matrix.
* `duration_of_covariance_matrix_in_natural_units` - The duration of the covariance matrix. It duration must be input in units that you know of (for instance the `time_period_per_unit` of a `SortedDataFrame`).
### Returns
* A `Hermitian`.
* A `Vector` of volatilities.
"""
function cov2cor_and_vol(mat::AbstractMatrix, duration_of_covariance_matrix::Dates.Period, duration_for_desired_vols::Dates.Period)
    cor, sdevs = cov2cor(mat)
    vols = sdevs / sqrt(  time_period_ratio(duration_of_covariance_matrix, duration_for_desired_vols)  )
    return Hermitian(cor), vols
end
function cov2cor_and_vol(mat::AbstractMatrix, duration_of_covariance_matrix::Real, duration_for_desired_vols::Real)
    cor, sdevs = cov2cor(mat)
    vols = sdevs / sqrt(  duration_of_covariance_matrix / duration_for_desired_vols  )
    return Hermitian(cor), vols
end
function cov2cor_and_vol(mat::AbstractMatrix, duration_of_covariance_matrix_in_natural_units::Real)
    cor, sdevs = cov2cor(mat)
    vols = sdevs / sqrt( duration_of_covariance_matrix_in_natural_units  )
    return Hermitian(cor), vols
end

"""
    cor2cov(cor::AbstractMatrix,sdevs::Vector{<:Real})

Converts a correlation matrix and some standard deviations into a `Hermitian` covariance matrix.
### Inputs
* `cor` - A correlation matrix.
* `sdevs` - A vector of standard deviations (not volatilities).
### Returns
* A `Hermitian`.
"""
function cor2cov(cor::AbstractMatrix, sdevs::Vector{<:Real})
    mat = cor .* (sdevs * transpose(sdevs))
    return Hermitian(mat)
end

"""
    covariance(cm::CovarianceMatrix, period::Dates.Period = cm.time_period_per_unit, assets::Vector{Symbol} = cm.labels)

This makes a `Hermitian` matrix for the covariance matrix over some duration.
### Inputs
* `cm` - A `CovarianceMatrix` struct.
* `period` - A duration for which you want a covariance matrix. This should be in a Dates.Period.
* `assets` - What assets in include in the covariance matrix.
### Returns
* A `Hermitian`. The labelling of assets for each row/column is as per the input `assets` vector.
"""
function covariance(cm::CovarianceMatrix, period::Dates.Period = cm.time_period_per_unit, assets::Vector{Symbol} = cm.labels)
    cm2 = rearrange(cm, assets)
    sds = convert_vol(cm2.volatility, cm2.time_period_per_unit, period)
    return cor2cov(cm2.correlation, sds)
end

"""
    construct_matrix_from_eigen(eigenvalues::Vector{<:Real},
                                eigenvectors::Matrix{<:Real})

Constructs a matrix from its eigenvalue decomposition.
### Inputs
* `eigenvalues` - A vector of eigenvalues.
* `eigenvectors` - A matrix of eigenvectors. The i'th column corresponds to the i'th eigenvalue.
### Returns
* A `Matrix`.
"""
function construct_matrix_from_eigen(eigenvalues::Vector{<:Real}, eigenvectors::Matrix{<:Real})
    value_matrix = zeros(length(eigenvalues), length(eigenvalues))
    value_matrix[diagind(value_matrix)] = eigenvalues
    final_mat = eigenvectors * value_matrix * transpose(eigenvectors)
    return Matrix(final_mat)
end

"""
    simple_differencing(new::Vector,old::Vector)

Does simple differencing of two vectors.
### Inputs
* `new` - A `Vector` of `Real`s.
* `old` - A `Vector` of `Real`s
### Returns
* A `Vector` of `Real`s
"""
simple_differencing(new::Vector,old::Vector) =  (new .- old)

"""
    get_returns(dd::DataFrame; rescale_for_duration::Bool = false)

Converts a long format `DataFrame` of prices into a `DataFrame` of returns.
### Inputs
* `dd` - A `DataFrame` with a column called :Time and all other columns being asset prices in each period.
* `rescale_for_duration` - Should returns be rescaled to reflect a common time interval.
### Returns
* A `DataFrame` of returns.
"""
function get_returns(dd::DataFrame; rescale_for_duration::Bool = false)
    N = nrow(dd)
    assets = setdiff(Symbol.(collect(names(dd))), [:Time])
    dd_mat = Array{Float64,2}(dd[1:N,assets])
    diffs  = reduce(hcat, map(i -> simple_differencing(dd_mat[2:end,i] , dd_mat[1:(end-1),i]) , 1:length(assets)) )
    dd_returns = diffs
    if rescale_for_duration
        time_diffs = dd[2:N,:Time] - dd[1:(N-1),:Time]
        dd_returns = (1 ./ sqrt.(time_diffs)) .*  diffs
    end
    dd2 = DataFrame(dd_returns, assets)
    return dd2
end

function weighted_mean(x::Vector, w::Vector)
    if length(x) == 0 return NaN end
    if length(x) == 1 return x[1] end
    if sum(abs.(w)) < 100*eps() return mean(x) end
    return sum(x .* w) / sum(w)
end

is_missing_nan_inf(x) = (ismissing(x) | isnan(x)) | isinf(x)

"""
    combine_covariance_matrices(vect::Vector{CovarianceMatrix{T}},
                                cor_weights::Vector{<:Real} = repeat([1.0], length(vect)),
                                vol_weights::Vector{<:Real} = cor_weights,
                                time_period_per_unit::Union{Missing,Dates.Period} = vect[1].time_period_per_unit) where T<:Real

Combines a vector of `CovarianceMatrix` structs into one `CovarianceMatrix` struct.
### Inputs
* `vect` - A vector of `CovarianceMatrix` structs.
* `cor_weights` - A vector for how much to weight the correlations from each covariance matrix (by default they will be equalweighted).
* `vol_weights` - A vector for how much to weight the volatilities from each covariance matrix (by default they will be equalweighted).
* `time_period_per_unit` - What time period should the volatilities be scaled to.
### Returns
* A matrix and a vector of labels for each row/column of the matrix.
"""
function combine_covariance_matrices(vect::Vector{CovarianceMatrix{T}}, cor_weights::Vector{<:Real} = repeat([1.0], length(vect)), vol_weights::Vector{<:Real} = cor_weights,
                                     time_period_per_unit::Union{Missing,Dates.Period} = vect[1].time_period_per_unit) where T<:Real
    if length(vect) < 1
        error("An empty vector of covariance matrices was input. So not possible to combine.")
    end
    all_labels = union(map(x -> x.labels, vect)...)
    dims = length(all_labels)
    R = promote_type(map(x -> eltype(x.correlation), vect)...)
    new_mat  = Array{R,2}(undef, dims, dims)
    new_vols = Array{R,1}(undef, dims)
    for row in 1:dims
        row_label = all_labels[row]
        for col in (row+1):dims
            col_label = all_labels[col]
            correls = map(i -> get_correlation(vect[i], row_label, col_label),  1:length(vect))
            valid_entries = setdiff(1:length(correls), findall(is_missing_nan_inf.(correls)))
            new_mat[row,col] = weighted_mean(correls[valid_entries], cor_weights[valid_entries])
        end
        vols = map(i -> convert_vol(get_volatility(vect[i], row_label), vect[i].time_period_per_unit, time_period_per_unit )             ,  1:length(vect))

        valid_entries = setdiff(1:length(vols), findall(is_missing_nan_inf.(vols)))
        new_vols[row] = weighted_mean(vols[valid_entries], vol_weights[valid_entries])
    end
    new_mat[diagind(new_mat)] .= 1
    hermitian_new_mat = Hermitian(new_mat)
    return CovarianceMatrix(hermitian_new_mat, new_vols, all_labels, time_period_per_unit)
end

"""
    rearrange(cm::CovarianceMatrix, labels::Vector{Symbol})

Rearrange the order of labels in a `CovarianceMatrix`.
### Takes
* `cm` - A `CovarianceMatrix`.
* `labels` - A `Vector` of labels.
* `time_period_per_unit` - The time period you want for the resultant Covariance Matrix
### Returns
* A `CovarianceMatrix`.
"""
function rearrange(cm::CovarianceMatrix, labels::Vector{Symbol},
                   time_period_per_unit::Union{Missing,Dates.Period} = cm.time_period_per_unit)
  if length(setdiff(labels, cm.labels)) > 0 error("You put in labels that are not in the CovarianceMatrix") end
  same_assets = (length(cm.labels) == length(labels)) && (all(cm.labels .== labels))
  same_vol = (time_period_per_unit == cm.time_period_per_unit)
  if (same_assets & same_vol) return cm end
  reordering = map(x -> findfirst(x .== cm.labels)[1], labels)
  Acor = same_assets ? cm.correlation[reordering,reordering] : Hermitian(cm.correlation[reordering,reordering])
  Avol = same_vol ? cm.volatility[reordering] : convert_vol(cm.volatility[reordering], cm.time_period_per_unit, time_period_per_unit)
  return CovarianceMatrix(Acor, Avol, labels, time_period_per_unit)
end

"""
    squared_frobenius_distance(x1::AbstractMatrix, x2::AbstractMatrix = x1)

Returns the squared frobenius distance between two matrices. This is a real.
### Inputs
* `x1` The first matrix.
* `x2` The second matrix.
### Returns
* A Scalar.
"""
function squared_frobenius_distance(x1::AbstractMatrix, x2::AbstractMatrix = x1)
    return squared_frobenius(x1 .- x2)
end
"""
    squared_frobenius(x1::AbstractMatrix)

Returns the squared frobenius norm of a matrix. This is a real.
### Inputs
* `x1` The matrix that you want the squared frobenius norm for.
### Returns
* A Scalar.
"""
function squared_frobenius(x1::AbstractMatrix)
    p = size(x1)[1]
    return tr(x1 * transpose(x1))/p
end

"""
    time_between_refreshes(ts::SortedDataFrame;
                           assets::Vector{Symbol} = get_assets(ts))

Get a `DataFrame` showing how many time is between each refresh and how many ticks in total.
### Inputs
* `ts` - Tick data.
* `assets` - A `Vector` of labels.
### Returns
* A `DataFrame` summarising the average number of time between ticks for each asset.
"""
function time_between_refreshes(ts::SortedDataFrame; assets::Vector{Symbol} = get_assets(ts))
    total_secs = maximum(ts.df[:,ts.time]) - minimum(ts.df[:,ts.time])
    R = eltype(ts.df[:,ts.value])
    N = nrow(ts.df)
    dd = DataFrame(Asset = Array{Symbol,1}([]), time_between_ticks = Array{R,1}([]), number_of_ticks = Array{Int64,1}([]))
    for a in assets
        asset_rows = ts.groupingrows[a]
        dff = ts.df[asset_rows,:]
        number_of_ticks = length(asset_rows)
        push!(dd, Dict(:Asset => a, :time_between_ticks => number_of_ticks/total_secs, :number_of_ticks => number_of_ticks))
    end
    sort!(dd, :time_between_ticks)
    return dd
end


text:
"""
	struct ExtractScalar{T} <: AbstractExtractor
		c::T
		s::T
		uniontypes::Bool
	end

Extracts a numerical value, centred by subtracting `c` and scaled by multiplying by `s`.
Strings are converted to numbers.

The extractor returns `ArrayNode{Matrix{Union{Missing, Int64}},Nothing}` or it subtypes.
If passed `missing`, it extracts missing values which Mill understands and can work with.

The `uniontypes` field determines whether extractor may or may not accept `missing`.
If `uniontypes` is false, it does not accept missing values. If `uniontypes` is true, it accepts missing values,
and always returns Mill structure of type Union{Missing, T} due to type stability reasons.

It can be created also using `extractscalar(Float32, 5, 2)`

# Example
```jldoctest
julia> ExtractScalar(Float32, 2, 3, true)(1)
1×1 Mill.ArrayNode{Matrix{Union{Missing, Float32}}, Nothing}:
 -3.0f0

julia> ExtractScalar(Float32, 2, 3, true)(missing)
1×1 Mill.ArrayNode{Matrix{Union{Missing, Float32}}, Nothing}:
 missing

julia> ExtractScalar(Float32, 2, 3, false)(1)
1×1 Mill.ArrayNode{Matrix{Float32}, Nothing}:
 -3.0

```
"""
struct ExtractScalar{T} <: AbstractExtractor
	c::T
	s::T
	uniontypes::Bool
end

ExtractScalar(T = Float32) = ExtractScalar(zero(T), one(T), true)
ExtractScalar(T, c, s, uniontypes = true) = ExtractScalar(T(c), T(s), uniontypes)

extractscalar(::Type{T}, m = zero(T), s = one(T), uniontypes = true) where {T<:Number} = ExtractScalar(T, m, s, uniontypes)
extractscalar(::Type{T}, uniontypes::Bool) where {T<:Number} = ExtractScalar(T, zero(T), one(T), uniontypes)
function extractscalar(::Type{T}, e::Entry, uniontypes = true) where {T<:Number}
	if unify_types(e) <: AbstractString
		values = parse.(T, keys(e.counts))
	else
		values = keys(e.counts)
	end
	min_val = minimum(values)
	max_val = maximum(values)
	c = min_val
	s = max_val == min_val ? 1 : 1 / (max_val - min_val)
	ExtractScalar(Float32(c), Float32(s), uniontypes)
end

_fill_and_normalize(s::ExtractScalar{T}, v::T) where {T} = s.s .* (fill(v,1,1) .- s.c)
stabilize_types_scalar(s::ExtractScalar{T}, x) where {T} = s.uniontypes ? Matrix{Union{Missing, T}}(x) : x
make_missing_scalar(s::ExtractScalar{T}, v, store_input) where {T} =
	s.uniontypes ?
	_make_array_node(Matrix{Union{Missing, T}}(fill(missing,1,1)), fill(v,1,1), store_input) :
	error("This extractor does not support missing values")
make_empty_scalar(s::ExtractScalar{T}, store_input) where {T} =
	_make_array_node(stabilize_types_scalar(s, fill(zero(T),1,0)), fill(undef,1,0), store_input)

(s::ExtractScalar)(v; store_input=false) = make_missing_scalar(s, v, store_input)
(s::ExtractScalar{T})(v::MissingOrNothing; store_input=false) where {T} = make_missing_scalar(s, v, store_input)
(s::ExtractScalar{T})(v::ExtractEmpty; store_input=false) where {T} = make_empty_scalar(s, store_input)
(s::ExtractScalar{T})(v::Number; store_input=false) where {T} =
	_make_array_node(stabilize_types_scalar(s, _fill_and_normalize(s, T(v))), fill(v,1,1), store_input)
function (s::ExtractScalar{T})(v::AbstractString; store_input=false) where {T}
	w = tryparse(T,v)
	isnothing(w) && return make_missing_scalar(s, v, store_input)
	_make_array_node(stabilize_types_scalar(s, _fill_and_normalize(s, T(w))), fill(v,1,1), store_input)
end

Base.length(e::ExtractScalar) = 1

# data type has different hashes for each patch version of julia
# see https://discourse.julialang.org/t/datatype-hash-differs-per-patch-version/48827
Base.hash(e::ExtractScalar{T}, h::UInt) where {T} = hash((e.c, e.s, e.uniontypes), h)
Base.:(==)(e1::ExtractScalar{T}, e2::ExtractScalar{T}) where {T} = e1.c === e2.c && e1.s === e2.s && e1.uniontypes === e2.uniontypes