text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using Test
using OpticSim.GlassCat

using Base.Filesystem
using DataFrames: DataFrame
using StaticArrays
using Unitful
using Unitful.DefaultSymbols

@testset "GlassCat" begin
    @testset "Build Tests" begin
        # check that all automatic downloads are working
        # this shouldn't be a test because we cannot guarantee that all downloads will work. Random network issues, changes in webpages, etc. can temporarily prevent downloads.
        # for catname in split("HOYA NIKON OHARA SCHOTT SUMITA")
        #     agffile = joinpath(GlassCat.AGF_DIR, catname * ".agf")
        #     @test isfile(agffile)
        # end

        # check that particularly problematic glasses are parsing correctly
        @test !isnan(NIKON.LLF6.C10)
    end

    @testset "sources.jl" begin
        @testset "add_agf" begin
            tmpdir = mktempdir()
            agfdir = mktempdir(tmpdir)
            sourcefile, _ = mktemp(tmpdir)

            @testset "bad implicit catalog names" begin
                for agffile in split("a 1.agf a.Agf")
                    @test_logs (:error, "invalid implicit catalog name \"$agffile\". Should be purely alphabetical with a .agf/.AGF extension.") add_agf(agffile; agfdir, sourcefile)
                end
            end

            @testset "file not found" begin
                @test_logs(
                    (:info, "Downloading source file from nonexistentfile.agf"),
                    (:error, ArgumentError("missing or unsupported scheme in URL (expected http(s) or ws(s)): nonexistentfile.agf")),
                    (:error, "failed to download from nonexistentfile.agf"),
                    add_agf("nonexistentfile.agf"; agfdir, sourcefile)
                )
            end

            @testset "add to source file" begin
                for name in split("a b")
                    open(joinpath(tmpdir, "$name.agf"), "w") do io
                        write(io, "")
                    end
                end
                empty_sha = "e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855"

                @test isempty(readlines(sourcefile))

                add_agf(joinpath(tmpdir, "a.agf"); agfdir, sourcefile, rebuild=false)
                @test length(readlines(sourcefile)) === 1
                @test readlines(sourcefile)[1] === "A $empty_sha"

                add_agf(joinpath(tmpdir, "b.agf"); agfdir, sourcefile, rebuild=false)
                @test length(readlines(sourcefile)) === 2
                @test readlines(sourcefile)[1] === "A $empty_sha"
                @test readlines(sourcefile)[2] === "B $empty_sha"

                @test_logs(
                    (:error, "adding the catalog name \"A\" would create a duplicate entry in source file $sourcefile"),
                    add_agf(joinpath(tmpdir, "a.agf"); agfdir, sourcefile)
                )
            end

            # TODO rebuild=true
        end

        # integration test
        @testset "verify_sources!" begin
            tmpdir = mktempdir()
            agfdir = mktempdir(tmpdir)

            sources = split.(readlines(GlassCat.SOURCES_PATH))
            GlassCat.verify_sources!(sources, agfdir)

            # this doesn't work if any of the glass catalogs can't be downloaded. Need a different test
            # @test first.(sources) == first.(split.(readlines(GlassCat.SOURCES_PATH)))

            # TODO missing_sources
        end

        # TODO unit tests
    end

    @testset "generate.jl" begin
        CATALOG_NAME = "TEST_CAT"
        SOURCE_DIR = joinpath(@__DIR__, "..", "..", "test")
        SOURCE_FILE = joinpath(SOURCE_DIR, "$(CATALOG_NAME).agf")

        TMP_DIR = mktempdir()
        MAIN_FILE = joinpath(TMP_DIR, "AGF_TEST_CAT.jl")

        TEST_CAT_VALUES = DataFrame(
            name = split("MG463_ MG523 _5MG448 MG452_STAR"),
            raw_name = split("MG463& MG523 5MG448 MG452*"),
            dispform = repeat([2], 4),
            Nd = repeat([1.0], 4),
            Vd = repeat([0.0], 4),
            exclude_sub = [5, 0, 0, 0],
            status = [6, 0, 0, 0],
            meltfreq = repeat([-1], 4),
            TCE = [19.8, 18.3, 30.2, 25.7],
            p = [4.63, 5.23, 4.48, 4.52],
            ΔPgF = repeat([0.0], 4),
            ignore_thermal_exp = [3, 0, 0, 0],
            C1 = [6.781409960E+00, 4.401393677E+00, 8.427948454E+00, 4.656568861E+00],
            C2 = [9.868538020E-02, 2.234722960E-02, 6.325268390E-02, 1.170534894E-01],
            C3 = [3.469529960E-03, 4.562935070E+00, -2.441576539E+00, 1.409322578E+00],
            C4 = [2.277823700E+00, 3.239318260E-01, 1.997453950E-03, 2.445302672E-04],
            C5 = [2.126055280E+00, 6.863473749E-01, 4.975938104E-01, 6.352738622E-01],
            C6 = [3.227490100E+03, 1.625888344E+03, 1.443307391E+03, 1.397305161E+03],
            C7 = repeat([0.0], 4),
            C8 = repeat([0.0], 4),
            C9 = repeat([0.0], 4),
            C10 = repeat([0.0], 4),
            D₀ = [2.226000000E-05, 1.122401430E-04, -1.718689957E-05, -2.545621409E-07],
            D₁ = [5.150100000E-08, 2.868995096E-13, -6.499841076E-13, 1.258087492E-10],
            D₂ = [1.309700000E-10, -4.157332734E-16, -4.153311088E-16, 1.119703940E-11],
            E₀ = [4.382400000E-05, 9.532969159E-05, 8.220254442E-06, 2.184273853E-05],
            E₁ = [4.966000000E-07, -1.045247704E-11, 6.086168578E-12, -2.788674395E-09],
            λₜₖ = [3.728000000E-01, 6.338475915E-01, 3.610058073E-01, -1.948154058E+00],
            temp = [2.002000000E+01, 2.000000000E+01, 2.000000000E+01, 2.000000000E+01],
            relcost = [1.1, -1.0, -1.0, -1.0],
            CR = [1.2, -1.0, -1.0, -1.0],
            FR = [1.3, -1.0, -1.0, -1.0],
            SR = [1.4, -1.0, -1.0, -1.0],
            AR = [1.5, -1.0, -1.0, -1.0],
            PR = [1.6, -1.0, -1.0, -1.0],
            λmin = repeat([2.0], 4),
            λmax = repeat([14.0], 4),
            transmission = [
                [[2.0, 1.0, 3.0], [14.0, 4.0, 5.0]],
                [[2.0, 1.0, 1.0], [14.0, 1.0, 1.0]],
                [[2.0, 1.0, 1.0], [14.0, 1.0, 1.0]],
                [[2.0, 1.0, 1.0], [14.0, 1.0, 1.0]],
            ],
        )
        FIELDS = names(TEST_CAT_VALUES)[2:end]

        @testset "Parsing Tests" begin
            cat = GlassCat.agffile_to_catalog(SOURCE_FILE)

            for glass in eachrow(TEST_CAT_VALUES)
                name = glass["name"]
                @test name ∈ keys(cat)
                for field in FIELDS
                    @test glass[field] == cat[name][field]
                end
            end
        end

        @testset "Module Gen Tests" begin
            OpticSim.GlassCat.generate_jls([CATALOG_NAME], MAIN_FILE, TMP_DIR, SOURCE_DIR, test=true)
            include(MAIN_FILE)

            for row in eachrow(TEST_CAT_VALUES)
                name = row["name"]
                @test "$CATALOG_NAME.$name" ∈ TEST_GLASS_NAMES

                glass = getfield(getfield(Main, Symbol(CATALOG_NAME)), Symbol(name))
                @test glass ∈ TEST_GLASSES

                for field in FIELDS
                    if field === "raw_name"
                    elseif field === "transmission"
                        @test row[field] == getfield(glass, Symbol(field))[1:length(row[field])]
                        @test zero(SVector{100-length(row[field]),SVector{3,Float64}}) == getfield(glass, Symbol(field))[length(row[field])+1:end]
                    else
                        @test row[field] == getfield(glass, Symbol(field))
                    end
                end
            end

            # these used to be in the "Glass Tests" testset, but they rely on the generated AGF_TEST_CAT.jl file
            g = TEST_CAT.MG523
            @test index(g, ((g.λmin + g.λmax) / 2)μm) ≈ 3.1560980389455593 atol=1e-14
            @test_throws ErrorException index(g, (g.λmin - 1)μm)
            @test_throws ErrorException index(g, (g.λmax + 1)μm)
        end
    end

    @testset "Glass Tests" begin
        @test absairindex(500nm) ≈ 1.0002741948670688 atol=1e-14
        @test absairindex(600nm, temperature=35°C, pressure=2.0) ≈ 1.0005179096900811 atol=1e-14
        @test index(Air, 500nm) == 1.0
        @test index(Air, 600nm) == 1.0

        # test against true values
        g = OpticSim.Examples.Examples_N_BK7
        @test index(g, 533nm) ≈ 1.519417351519283 atol=1e-14
        @test index(g, 533nm, temperature=35°C) ≈ 1.519462486258311 atol=1e-14
        @test index(g, 533nm, pressure=2.0) ≈ 1.518994119690216 atol=1e-14
        @test index(g, 533nm, temperature=35°C, pressure=2.0) ≈ 1.519059871499476 atol=1e-14

        # test transmission
        @test absorption(Air, 500nm) == 0.0
        @test absorption(Air, 600nm) == 0.0
        # TODO these are currently taken from this package (i.e. regression tests), ideally we would get true values somehow
        @test absorption(g, 500nm) == 0.0002407228930225164
        @test absorption(g, 600nm) == 0.00022060722752440445
        @test absorption(g, 3000nm) == 0.04086604990127926
        @test absorption(g, 600nm, temperature=35°C, pressure=2.0) == 0.00022075540719494738

        # test that everything is alloc-less
        @test (@wrappedallocs absorption(g, 600nm)) == 0
        @test (@wrappedallocs index(g, 533nm)) == 0
        @test (@allocated OpticSim.Examples.Examples_N_BK7) == 0

        @test glassforid(glassid(OpticSim.Examples.Examples_N_BK7)) == OpticSim.Examples.Examples_N_BK7

        @test isair(Air) == true
        @test isair(OpticSim.Examples.Examples_N_BK7) == false

        # test MIL and model glasses
        fit_acc = 0.001
        bk7 = glassfromMIL(517642)
        @test glassforid(glassid(bk7)) == bk7
        @test index(bk7, 0.5875618) ≈ 1.517 atol=fit_acc
        @test index(bk7, 0.533) ≈ 1.519417351519283 atol=fit_acc
        @test index(bk7, 0.743) ≈ 1.511997032563557 atol=fit_acc
        @test index(bk7, 533nm, temperature=35°C, pressure=2.0) ≈ 1.519059871499476 atol=fit_acc

        t2 = glassfromMIL(1.135635)
        @test index(t2, 0.5875618) ≈ 2.135 atol=fit_acc

        fit_acc = 0.0001
        bk7 = modelglass(1.5168, 64.167336, -0.0009)
        @test glassforid(glassid(bk7)) == bk7
        @test index(bk7, 0.5875618) ≈ 1.5168 atol=fit_acc
        @test index(bk7, 0.533) ≈ 1.519417351519283 atol=fit_acc
        @test index(bk7, 0.743) ≈ 1.511997032563557 atol=fit_acc
        @test index(bk7, 533nm, temperature=35°C, pressure=2.0) ≈ 1.519059871499476 atol=fit_acc

        # test other glass
        @test index(CARGILLE.OG0608, 0.578) == 1.4596475735607324

        # make sure that the other functions work
        plot_indices(OpticSim.Examples.Examples_N_BK7; polyfit = true, fiterror = true)
    end

    @testset "Air.jl" begin
        @test repr(Air) === "Air"
        @test glassname(Air) === "GlassCat.Air"

        info_lines = [
            "GlassCat.Air",
            "Material representing air, RI is always 1.0 at system temperature and pressure, absorption is always 0.0.",
            ""
        ]
        io = IOBuffer()
        info(io, Air)
        @test String(take!(io)) === join(info_lines, '\n')
    end

    @testset "GlassTypes.jl" begin
        @test repr(GlassID(GlassCat.MODEL, 1)) === "MODEL:1"
        @test repr(GlassID(GlassCat.MIL, 1)) === "MIL:1"
        @test repr(GlassID(GlassCat.AGF, 1)) === "AGF:1"
        @test repr(GlassID(GlassCat.OTHER, 1)) === "OTHER:1"
        @test repr(GlassID(GlassCat.AIR, 1)) === "AIR:1" # TODO should this fail?

        empty_args = SVector{36,Union{Int,Nothing}}(repeat([0], 27)..., nothing, repeat([0], 8)...)
        @test repr(GlassCat.Glass(GlassID(GlassCat.MODEL, 1), empty_args...)) === "GlassCat.ModelGlass.1"
        @test repr(GlassCat.Glass(GlassID(GlassCat.MIL, 1), empty_args...)) === "GlassCat.GlassFromMIL.1"
        @test repr(GlassCat.Glass(GlassID(GlassCat.AGF, 1), empty_args...)) === "NIKON.SF9" # fails if sources.txt is changed
        @test repr(GlassCat.Glass(GlassID(GlassCat.OTHER, 1), empty_args...)) === "CARGILLE.OG0607"
        @test repr(GlassCat.Glass(GlassID(GlassCat.AIR, 1), empty_args...)) === "GlassCat.Air" # repeated code

        @test glassforid(GlassID(GlassCat.AIR, 1)) === Air
        @test glassforid(GlassID(GlassCat.OTHER, 1)) === CARGILLE.OG0607

        info_lines = [
            "ID:                                                OTHER:2",
            "Dispersion formula:                                Cauchy (-2)",
            "Dispersion formula coefficients:",
            "     A:                                            1.44503",
            "     B:                                            0.0044096",
            "     C:                                            -2.85878e-5",
            "     D:                                            0.0",
            "     E:                                            0.0",
            "     F:                                            0.0",
            "Valid wavelengths:                                 0.32μm to 1.55μm",
            "Reference temperature:                             25.0°C",
            "Thermal ΔRI coefficients:",
            "     D₀:                                           -0.0009083144750540808",
            "     D₁:                                           0.0",
            "     D₂:                                           0.0",
            "     E₀:                                           0.0",
            "     E₁:                                           0.0",
            "     λₜₖ:                                          0.0",
            "TCE (÷1e-6):                                       700.0",
            "Ignore thermal expansion:                          false",
            "Density (p):                                       0.878g/m³",
            "ΔPgF:                                              0.008",
            "RI at sodium D-Line (587nm):                       1.457587",
            "Abbe Number:                                       57.19833",
            "Cost relative to N_BK7:                            ?",
            "Status:                                            Standard (0)",
            "Melt frequency:                                    ?",
            "Exclude substitution:                              false",
            "Transmission data:",
            "     Wavelength   Transmission      Thickness",
            "         0.32μm           0.15         10.0mm",
            "        0.365μm           0.12        100.0mm",
            "       0.4047μm           0.42        100.0mm",
            "         0.48μm           0.78        100.0mm",
            "       0.4861μm           0.79        100.0mm",
            "       0.5461μm           0.86        100.0mm",
            "       0.5893μm            0.9        100.0mm",
            "       0.6328μm           0.92        100.0mm",
            "       0.6439μm            0.9        100.0mm",
            "       0.6563μm           0.92        100.0mm",
            "       0.6943μm           0.98        100.0mm",
            "         0.84μm           0.99        100.0mm",
            "      0.10648μm           0.61        100.0mm",
            "         0.13μm           0.39        100.0mm",
            "        0.155μm           0.11        100.0mm",
            ""
        ]
        io = IOBuffer()
        info(io, CARGILLE.OG0607)
        @test String(take!(io)) === join(info_lines, '\n')

        # TODO the rest of the string tests
    end

    @testset "search.jl" begin
        #this test set doesn't work as is. It assumes that all the glass catalogs have downloaded correctly which may not happen. Commenting them all out till we figure out a better set of tests.

        # @test glasscatalogs() == [
        #     CARGILLE,
        #     HOYA,
        #     NIKON,
        #     OHARA,
        #     SCHOTT,
        #     SUMITA,
        # ]

        # @test glassnames(CARGILLE) == [
        #     :OG0607,
        #     :OG0608,
        #     :OG081160,
        # ]

        # @test first.(glassnames()) == [
        #     CARGILLE,
        #     HOYA,
        #     NIKON,
        #     OHARA,
        #     SCHOTT,
        #     SUMITA,
        # ]
        # @test length.(last.(glassnames())) == [
        #     3,
        #     210,
        #     379,
        #     160,
        #     160,
        #     178,
        # ]

        # @test findglass(x -> (x.Nd > 2.1 && x.λmin < 0.5 && x.λmax > 0.9)) == [
        #     HOYA.E_FDS3,
        #     SUMITA.K_PSFn214P,
        #     SUMITA.K_PSFn214P_M_,
        # ]

        # TODO _child_modules() unit test
    end
end # testset GlassCat


text:
"""
Define inference problem by supplying
    s: observation locations
    X: covariates for observed locations
    z: observed values/labels
    X0: matrix of covariates at prediction location
    s0: prediction location
"""
struct setting{T<:Array{Float64, 2}, S<:Array{Float64, 1}}
    s::T
    s0::T
    X::T
    X0::T
    z::S
end

struct trainingData{T<:Array{Float64, 2}, S<:Array{Float64, 1}}
    s::T
    X::T
    z::S
end

struct testingData{T<:Array{Float64, 2}}
    s0::T
    X0::T
end

ext: jl

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using ITensors
using Zygote
using OptimKit

function ising(n; h)
  ℋ = Sum{Op}()
  for j in 1:(n - 1)
    ℋ -= "Z", j, "Z", j + 1
  end
  for j in 1:n
    ℋ += h, "X", j
  end
  return ℋ
end

n = 4
s = siteinds("S=1/2", n)

β = 1.0
hᶠ = 1.1

ℋᶠ = ising(n; h=hᶠ)
𝒰ᶠ = exp(-β * ℋᶠ; alg=Trotter{1}(5))

Uᶠ = Prod{ITensor}(𝒰ᶠ, s)

ψ = prod(MPS(s, "0"))

# Target state
Uᶠψ = Uᶠ(ψ)

function loss(h)
  ℋ = ising(n; h=h[1])
  𝒰ʰ = exp(-β * ℋ; alg=Trotter{1}(5))
  Uʰ = Prod{ITensor}(𝒰ʰ, s)
  Uʰψ = Uʰ(ψ)
  return -abs(inner(Uᶠψ, Uʰψ))^2 / (norm(Uᶠψ) * norm(Uʰψ))^2
end

h⁰ = [0.0]
@show loss(h⁰)
@show loss'(h⁰)

@show loss(hᶠ)
@show loss'(hᶠ)

loss_∇loss(h) = (loss(h), convert(Vector, loss'(h)))
algorithm = LBFGS(; gradtol=1e-3, verbosity=2)
hᵒᵖᵗ, _ = optimize(loss_∇loss, h⁰, algorithm)

@show loss(hᵒᵖᵗ)
@show loss'(hᵒᵖᵗ)


text:
function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintName,
    c::MOI.ConstraintIndex,
    name::String,
)
    info = _info(model, c)
    info.name = name
    cpo_java_addable_setname(model.inner, info.constraint, name)
    model.name_to_constraint = nothing
    return
end

function MOI.get(model::Optimizer, ::Type{MOI.ConstraintIndex}, name::String)
    if model.name_to_constraint === nothing
        _rebuild_name_to_constraint(model)
    end
    return get(model.name_to_constraint, name, nothing)
end

function MOI.get(
    model::Optimizer,
    C::Type{MOI.ConstraintIndex{F, S}},
    name::String,
) where {F, S}
    index = MOI.get(model, MOI.ConstraintIndex, name)
    if typeof(index) == C
        return index::MOI.ConstraintIndex{F, S}
    end
    return nothing
end

function _rebuild_name_to_constraint_add!(
    model::Optimizer,
    name::String,
    cindex::MOI.ConstraintIndex,
)
    if name == ""
        return
    end

    if haskey(model.name_to_constraint, name)
        error("Duplicate variable name detected: $(name)")
    end
    return model.name_to_constraint[name] = cindex
end

function _rebuild_name_to_constraint(model::Optimizer)
    model.name_to_constraint = Dict{String, MOI.ConstraintIndex}()

    # Variable bounds.
    for (index, info) in model.variable_info
        T = _variabletype_to_type(info.type)

        if _check_bound_compatible(model, index, MOI.EqualTo{T})
            _rebuild_name_to_constraint_add!(
                model,
                info.equalto_name,
                MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{T}}(
                    index.value,
                ),
            )
        end
        if _check_bound_compatible(model, index, MOI.Interval{T})
            _rebuild_name_to_constraint_add!(
                model,
                info.interval_name,
                MOI.ConstraintIndex{MOI.VariableIndex, MOI.Interval{T}}(
                    index.value,
                ),
            )
        end
        if _check_bound_compatible(model, index, MOI.LessThan{T})
            _rebuild_name_to_constraint_add!(
                model,
                info.ub_name,
                MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{T}}(
                    index.value,
                ),
            )
        end
        if _check_bound_compatible(model, index, MOI.GreaterThan{T})
            _rebuild_name_to_constraint_add!(
                model,
                info.lb_name,
                MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{T}}(
                    index.value,
                ),
            )
        end

        if info.type == INTEGER || info.integer !== nothing
            _rebuild_name_to_constraint_add!(model, info.integer_name, index)
        elseif info.type == BINARY || info.binary !== nothing
            _rebuild_name_to_constraint_add!(model, info.binary_name, index)
        end
    end

    # Other constraints.
    for (index, info) in model.constraint_info
        _rebuild_name_to_constraint_add!(
            model,
            info.name,
            MOI.ConstraintIndex{typeof(info.f), typeof(info.set)}(index.value),
        )
    end

    return
end


text:
"""
    struct SkeletonTriangulation{Dc,Dp,B} <: Triangulation{Dc,Dp}
      left::B
      right::B
    end

The inner constructor enforces `B<:BoundaryTriangulation`
"""
struct SkeletonTriangulation{Dc,Dp,B} <: Triangulation{Dc,Dp}
  left::B
  right::B
  function SkeletonTriangulation(left::B,right::B) where B<:Triangulation
    Dc = num_cell_dims(left)
    Dp = num_point_dims(left)
    @assert Dc + 1 == Dp
    new{Dc,Dp,B}(left,right)
  end
end

"""
    SkeletonTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
    SkeletonTriangulation(model::DiscreteModel)
"""
function SkeletonTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
  left_cell_around = 1
  left = BoundaryTriangulation(model,face_to_mask,left_cell_around)
  right_cell_around = 2
  right = BoundaryTriangulation(model,face_to_mask,right_cell_around)
  SkeletonTriangulation(left,right)
end

function SkeletonTriangulation(model::DiscreteModel)
  topo = get_grid_topology(model)
  D = num_cell_dims(model)
  face_to_mask = collect(Bool, .!get_isboundary_face(topo,D-1))
  SkeletonTriangulation(model,face_to_mask)
end

function SkeletonTriangulation(
  model::DiscreteModel,
  reffe::ReferenceFE,
  face_own_dofs::Vector{Vector{Int}})

  ncells = num_cells(model)
  cell_to_reffe = Fill(1,ncells)
  SkeletonTriangulation(model,[reffe],cell_to_reffe,[face_own_dofs])
end

function SkeletonTriangulation(
  model::DiscreteModel,
  reffes::Vector{<:ReferenceFE},
  cell_to_reffe::AbstractVector,
  reffe_to_face_own_dofs::Vector{Vector{Vector{Int}}})

  reffe_to_polytope = map(get_polytope,reffes)
  reffe_to_faces = map(get_faces,reffe_to_polytope)
  reffe_to_range = map(p->get_dimrange(p,num_dims(p)-1),reffe_to_polytope)

  reffe_to_lfacet_to_mask = _compute_reffe_to_facet_to_mask(
    reffe_to_faces,reffe_to_range,reffe_to_face_own_dofs)

  topo = get_grid_topology(model)

  D = num_cell_dims(model)
  nfacets = num_facets(model)
  facet_to_isboundary = get_isboundary_face(topo,D-1)
  cell_to_facets = get_faces(topo,D,D-1)
  facet_to_mask = _compute_disc_facet_mask(
    facet_to_isboundary,cell_to_facets,cell_to_reffe,reffe_to_lfacet_to_mask)

  SkeletonTriangulation(model,facet_to_mask)

end

function _compute_disc_facet_mask(
  facet_to_isboundary, cell_to_facets::Table, cell_to_reffe, reffe_to_lfacet_to_mask)
  nfacets = length(facet_to_isboundary)
  facet_to_mask = Vector{Bool}(undef,nfacets)
  ncells = length(cell_to_reffe)
  for cell in 1:ncells
    pini = cell_to_facets.ptrs[cell]
    pend = cell_to_facets.ptrs[cell+1]-1
    reffe = cell_to_reffe[cell]
    lfacet_to_mask = reffe_to_lfacet_to_mask[reffe]
    for (lfacet,p) in enumerate(pini:pend)
      facet = cell_to_facets.data[p]
      mask = lfacet_to_mask[lfacet]
      isinterior = ! facet_to_isboundary[facet]
      facet_to_mask[facet] = mask && isinterior
    end
  end
  facet_to_mask
end

function _compute_reffe_to_facet_to_mask(
  reffe_to_faces,reffe_to_range,reffe_to_face_own_dofs)

  reffe_to_facet_to_mask = Vector{Bool}[]
  nreffes = length(reffe_to_faces)
  for reffe in 1:nreffes
    range = reffe_to_range[reffe]
    facet_to_faces = reffe_to_faces[reffe][range]
    facet_to_own_dofs = reffe_to_face_own_dofs[reffe][range]
    face_to_own_dofs = reffe_to_face_own_dofs[reffe]
    nfacets = length(facet_to_faces)
    facet_to_mask = fill(false,nfacets)
    for facet in 1:nfacets
      n = length(facet_to_own_dofs[facet])
      faces = facet_to_faces[facet]
      for face in faces
        n += length(face_to_own_dofs[face])
      end
      facet_to_mask[facet] = n == 0
    end
    push!(reffe_to_facet_to_mask,facet_to_mask)
  end

  reffe_to_facet_to_mask
end

const IN = -1
const OUT = 1

"""
"""
function InterfaceTriangulation(model::DiscreteModel,cell_to_is_in::Vector{Bool})
  cell_to_inout = fill(Int8(OUT),length(cell_to_is_in))
  cell_to_inout[cell_to_is_in] .= IN
  InterfaceTriangulation(model,cell_to_inout)
end

function InterfaceTriangulation(model::DiscreteModel,cells_in,cells_out)
  cell_to_inout = fill(Int8(OUT),num_cells(model))
  cell_to_inout[cells_in] .= IN
  cell_to_inout[cells_out] .= OUT
  InterfaceTriangulation(model,cell_to_inout)
end

function InterfaceTriangulation(model::DiscreteModel,cell_to_inout::AbstractVector{<:Integer})

  D = num_cell_dims(model)
  facet_grid = Grid(ReferenceFE{D-1},model)
  cell_grid = Grid(ReferenceFE{D},model)

  topo = get_grid_topology(model)
  facet_to_cells = Table(get_faces(topo,D-1,D))
  ifacet_to_facet, facet_to_lcell_left, facet_to_lcell_right = _find_interface_facets(
    cell_to_inout, facet_to_cells)

  ifacet_trian = TriangulationPortion(facet_grid,ifacet_to_facet)

  left = GenericBoundaryTriangulation(ifacet_trian,cell_grid,topo,ifacet_to_facet,facet_to_lcell_left)
  right = GenericBoundaryTriangulation(ifacet_trian,cell_grid,topo,ifacet_to_facet,facet_to_lcell_right)

  SkeletonTriangulation(left,right)
end

function _find_interface_facets( cell_to_inout, facet_to_cells::Table)

  nifacets = 0
  for facet in 1:length(facet_to_cells)
    a = facet_to_cells.ptrs[facet]
    b = facet_to_cells.ptrs[facet+1]
    if b-a == 2
      cell1 = facet_to_cells.data[a]
      cell2 = facet_to_cells.data[a+1]
      inout_1 = cell_to_inout[cell1]
      inout_2 = cell_to_inout[cell2]
      if (inout_1== IN && inout_2==OUT) || (inout_1== OUT && inout_2==IN)
        nifacets += 1
      end
    end
  end

  T = eltype(eltype(facet_to_cells))
  ifacet_to_facet = zeros(T,nifacets)
  nfacets = length(facet_to_cells)
  facet_to_lcell_left = fill(Int8(1),nfacets)
  facet_to_lcell_right = fill(Int8(2),nfacets)

  nifacets = 0
  for facet in 1:length(facet_to_cells)
    a = facet_to_cells.ptrs[facet]
    b = facet_to_cells.ptrs[facet+1]
    if b-a == 2
      cell1 = facet_to_cells.data[a]
      cell2 = facet_to_cells.data[a+1]
      inout_1 = cell_to_inout[cell1]
      inout_2 = cell_to_inout[cell2]
      if (inout_1== IN && inout_2==OUT) || (inout_1== OUT && inout_2==IN)
        nifacets += 1
        ifacet_to_facet[nifacets] = facet
        if inout_1 == IN
          facet_to_lcell_left[facet] = 1
          facet_to_lcell_right[facet] = 2
        else
          facet_to_lcell_left[facet] = 2
          facet_to_lcell_right[facet] = 1
        end
      end
    end
  end

   ifacet_to_facet, facet_to_lcell_left, facet_to_lcell_right
 end

# Triangulation interface

function get_cell_coordinates(trian::SkeletonTriangulation)
  get_cell_coordinates(trian.left)
end

function get_reffes(trian::SkeletonTriangulation)
    get_reffes(trian.left)
end

function get_cell_type(trian::SkeletonTriangulation)
  get_cell_type(trian.left)
end

function get_cell_id(trian::SkeletonTriangulation)
  left = get_cell_id(trian.left)
  right = get_cell_id(trian.right)
  SkeletonPair(left,right)
end

function restrict(f::AbstractArray, trian::SkeletonTriangulation)
  left = restrict(f,trian.left)
  right = restrict(f,trian.right)
  SkeletonPair(left,right)
end

# Delegating into the left side

"""
    get_volume_triangulation(trian::SkeletonTriangulation)
"""
function get_volume_triangulation(trian::SkeletonTriangulation)
  get_volume_triangulation(trian.left)
end

"""
    get_normal_vector(trian::SkeletonTriangulation)
"""
function get_normal_vector(trian::SkeletonTriangulation)
  get_normal_vector(trian.left)
end

function TriangulationPortion(oldtrian::SkeletonTriangulation,cell_to_oldcell::Vector{Int})
  left = TriangulationPortion(oldtrian.left,cell_to_oldcell)
  right = TriangulationPortion(oldtrian.right,cell_to_oldcell)
  SkeletonTriangulation(left,right)
end

# Specific API

"""
"""
get_left_boundary(trian::SkeletonTriangulation) = trian.left

"""
"""
get_right_boundary(trian::SkeletonTriangulation) = trian.right


text:
module JavaTranspile

# @enum AccessLevelModifiers public private protected
# @enum Modifiers static final abstract synchronized volatile

# Java primitive types
# immutable byte value::Int8 end
# immutable short value::Int16 end
# immutable int value::Int32 end
# immutable long value::Int64 end
# immutable float value::Float32 end
# immutable double value::Float64 end
# immutable boolean value::Bool end
# immutable char value::Char end

# immutable String value::Base.String end

type Variable
    name::String
    annotations::Vector{String}
    modifiers::Vector{String}
    typ::String
    init::String
end

Variable() = Variable("", String[], String[], "", "")

function Base.show(io::IO, m::Variable)
    println(io, "Java.Variable: $(m.name)")
    println(io, "  annotations: $(join(m.annotations, ", "))")
    println(io, "   definition: $(join(m.modifiers, ' ')) $(m.typ) $(m.name)")
end

function Base.showcompact(io::IO, m::Variable)
    print(io, "$(join(m.modifiers, ' ')) $(m.typ) $(m.name)")
end

type Method
    name::String
    annotations::Vector{String}
    modifiers::Vector{String}
    returntype::String
    args::Vector{String}
    body::String
end

Method() = Method("", String[], String[], "", String[], "")

function Base.show(io::IO, m::Method)
    println(io, "Java.Method: $(m.name)")
    println(io, "  annotations: $(join(m.annotations, ", "))")
    println(io, "   definition: $(join(m.modifiers, ' ')) $(m.returntype == m.name ? "" : string(m.returntype, " "))$(m.name)($(join(m.args, ", ")))")
end

function Base.showcompact(io::IO, m::Method)
    print(io, "$(join(m.modifiers, ' ')) $(m.returntype == m.name ? "" : string(m.returntype, " "))$(m.name)($(join(m.args, ", ")))")
end

type Class
    name::String
    classtype::String
    package::String
    imports::Vector{String}
    annotations::Vector{String}
    classdef::String
    enums::Vector{String}
    variables::Vector{Variable}
    methods::Vector{Method}
    innerclasses::Vector{Class}
end

Class() = Class("", "", "", String[], String[], "", String[], Variable[], Method[], Class[])

function Base.show(io::IO, c::Class)
    println(io, "Java.Class: $(c.name)")
    println(io, "        package: $(chop(c.package))")
    println(io, "      # imports: $(length(c.imports))")
    println(io, "    annotations: $(join(c.annotations, ", "))")
    println(io, "     definition: $(chop(c.classdef))")
    if c.classtype == "enum"
    println(io, "          enums: ")
    for e in c.enums
        println(io, (" " ^ 16) * e)
    end
    end
    println(io, "         fields: ")
    buf = IOBuffer()
    for f in c.variables
        showcompact(buf, f)
        println(io, (" " ^ 16) * takebuf_string(buf))
    end
    println(io, "        methods: ")
    for m in c.methods
        showcompact(buf, m)
        println(io, (" " ^ 16) * takebuf_string(buf))
    end
    println(io, "  inner classes: ")
    for f in c.innerclasses
        println(io, (" " ^ 16) * f.name)
    end
end

function parseline(f)
    line = strip(readline(f))
    while !eof(f)
        if line == ""
            # if empty line, skip
            line = strip(readline(f))
        elseif startswith(line, "//")
            # if comment, skip line
            line = strip(readline(f))
        elseif startswith(line, "/*")
            # if multi-line comment, keep reading till end
            while !contains(line, "*/")
                line = strip(readline(f))
            end
            line = strip(readline(f))
        else
            break
        end
    end
    return line
end

function parsepackage(f, line, class)
    if startswith(line, "package")
        class.package = line
        line = parseline(f)
    end
    return line
end

function parseimports(f, line, class)
    while !eof(f)
        if startswith(line, "import")
            push!(class.imports, line)
            line = parseline(f)
        else
            break
        end
    end
    return line
end

function parseclassannotations(f, line, class)
    while !eof(f)
        if startswith(line, "@")
            if contains(line, "(") && !contains(line, ")")
                while !contains(line, ")")
                    line *= parseline(f)
                end
            end
            push!(class.annotations, line)
            line = parseline(f)
        else
            break
        end
    end
    return line
end

function parseclassdefinition(f, line, class)
    while !contains(line, "{")
        line *= string(" ", parseline(f))
    end
    class.classdef = line
    if contains(line, "class")
        # abstract or regular class
        if contains(line, "abstract")
            class.classtype = "abstract"
        else
            class.classtype = "class"
        end
    elseif contains(line, "interface")
        # interface class
        class.classtype = "interface"
    elseif contains(line, "enum")
        # enum class
        class.classtype = "enum"
        # parse comma-separated enum list until semicolon
        class.enums = String[strip(string(x)) for x in split(chop(parseuntil(f, ';')), ',')]
    else
        error("unsupported class type detected: $line")
    end
    return
end

function skipblock(f, bracecount)
    line = parseline(f)
    while true
        bracecount += count(c -> c == '{', line)
        bracecount -= count(c -> c == '}', line)
        bracecount <= 0 && break
        line = parseline(f)
    end
end

function nextline(f)
    ch = string(Char(Base.peek(f)))
    while ismatch(r"\s", ch)
        read(f, Char)
        eof(f) && break
        ch = string(Char(Base.peek(f)))
    end
end

function parseclassbody(f, line, class)
    annotations = String[]
    while !eof(f)
        nextline(f)
        token = parseuntil(f, 0, '<', '>', ' ', '(')
        if startswith(token, "//")
            parseline(f)
            continue
        elseif startswith(token, "/*")
            while !contains(token, "*/")
                token = parseline(f)
            end
            continue
        end
        if token == "}" || token == "};" || token == ""
            # end of class/file
            break
        elseif token[1] == '@'
            # method (or variable?) annotations
            token *= token[end] == '(' ? parseuntil(f, 1, '(', ')', ')') : ""
            push!(annotations, token)
        elseif token == "{"
            # initialization block
            skipblock(f, 1)
        elseif token == "static" && Char(Base.peek(f)) == '{'
            # static initialization block
            skipblock(f, 0)
        else
            # we're parsing an inner class, member variable, or method
            modifiers = String[]
            while token in ("public", "private", "protected", "static", "final", "abstract", "synchronized", "volatile")
                push!(modifiers, token)
                token = parseuntil(f, 0, '<', '>', ' ', '(')
            end
            # we're now positioned at the type (or method name for constructors)
            if token in ("class", "enum", "interface")
                # inner class
                ch = ' '
                while ch != '{'
                    ch = read(f, Char)
                end
                skipblock(f, 1)
                continue
            elseif token[end] == '('
                # class constructor
                name = token[1:end-1]
                returntype = name
            else
                # return type
                returntype = token
                # method or variable name
                token = parseuntil(f, 0, '<', '>', ' ', '(', ';')
                name = chop(token)
            end
            if token[end] == '('
                # we're parsing a method
                args = parseargs(chop(parseuntil(f, 1, '(', ')', ')')))
                body = string(read(f, Char))
                if body == ";"
                    body = ""
                else
                    body *= parseuntil(f, Int(body == "{"), '{', '}', '}')
                end
                push!(class.methods, Method(name, splice!(annotations, 1:length(annotations)), modifiers, returntype, args, body))
            else
                # we're parsing a variable
                init = token[end] == ';' ? "" : parseuntil(f, ';')
                push!(class.variables, Variable(name, splice!(annotations, 1:length(annotations)), modifiers, returntype, init))
            end
        end
    end
end

function parseuntil(io, bracecount, ob, cb, chars...)
    buf = IOBuffer()
    stripleadingwhitespace = true
    while !eof(io)
        ch = read(io, Char)
        if stripleadingwhitespace
            while ch == ' '
                ch = read(io, Char)
            end
            stripleadingwhitespace = false
        end
        write(buf, ch)
        if ch == ob
            bracecount += 1
        elseif ch == cb
            bracecount -= 1
        end
        detectclosing = false
        for c in chars
            ch == c && (detectclosing = true)
        end
        detectclosing && bracecount == 0 && break
    end
    return strip(takebuf_string(buf))
end
parseuntil(io, c) = parseuntil(io, 0, '<', '>', c)

function parseargs(args)
    length(args) <= 1 && return String[]
    return String[string(x) for x in split(args, ',')]
end

function parse(file)
    isfile(file) || throw(ArgumentError("'$file' is not a valid file"))
    class = Class()
    class.name = basename(file)[1:end-5]
    open(file) do f
        line = parseline(f)
        line = parsepackage(f, line, class)
        line = parseimports(f, line, class)
        line = parseclassannotations(f, line, class)
        parseclassdefinition(f, line, class)
        parseclassbody(f, line, class)
    end
    return class
end

function generateclient(class::Class, output_dir)

end

function generatetype(class::Class)

end

#TODO
 # JavaParser:
   # parse inner classes
   # setup proper test suite
   # return actual types in classes/methods?
   # don't show constructor type twice in method show?
 # write generateclient method
 # write generatetype method
 # rename this package to "JavaParser"
 # create "Java" submodule to hold type definitions, methods; export this submodule

# need:
 # transpile a regular class or enum to corresponding Julia type
 # transpile a "client" interface to corresponding Julia module + methods
 # Client generation
  # create a module for client interface w/ using Requests, JSON
  # create a method for each interface method
  # for each non primitive type:
   # find path from imports
   # parse custom class path
   # generate corresponding Julia type

end # module


text:
module VAEs

# """ Init VAE as struct consisting the network for the variable reduction (encoder), the mu and logsigma layer and the the decoder network """
using Flux
using Flux: throttle, params
import Random
using LinearAlgebra

struct VAE
	encoder
	mu
	logsigma
	decoder
end

function initvae(; nvariables, nhidden, nz)
   VAE(Flux.Chain(Flux.Dense(nvariables,nhidden,tanh)),
      Flux.Dense(nhidden,nz),Flux.Dense(nhidden,nz),
      Flux.Chain(Flux.Dense(nz,nhidden,tanh),
      Flux.Dense(nhidden,nvariables,Flux.NNlib.σ)))
end

function normalz(mu,logsigma)
	return mu  + exp(logsigma) * randn()
end

function decoderpen(vae)
	return 0.01 * sum(x->sum(x.^2),Flux.params(vae.decoder))
end

function kullbackleiblerqp(mu,logsigma)
	return 0.5 * sum(exp.(2 .* logsigma) + mu.^2 .- 1 .- 2 .* logsigma)
end

function reconstructionloss(xhat,x)
	ce(x,y) =  y * log(x + eps()) + (1 - y) * log(1 - x + eps())
	# logpdf(b::Bernoulli,y::Bool) = y * log(b.p + eps()) + (1-y) * log(1-b.p + eps())
	# return sum(Distributions.logpdf.(Bernoulli.(xhat),x))
	return sum(ce.(xhat,x))
end

""" The functions to sample from the vae conditional or unconditional of the data """
function getzfordata(vae,data)
	transformeddata = vae.encoder(data)

	mu,logsigma = (vae.mu(transformeddata).data,vae.logsigma(transformeddata).data)
	z = normalz.(mu,logsigma)
	z
end

function getsamplefordata(vae,data,prob=true)
	z = getzfordata(vae,data)
	samplegivenz(vae,z,prob)
end

function samplegivenz(vae,z,prob=prob)
	s = vae.decoder(z).data
	if prob == false
		s = Float64.(s .> rand(length(s)))
	end
	s
end

function getlatentvarsfordata(vae,data,prob=true,varlevel=1)
	z = getzfordata(vae,data)
	getlatentvarsforz(vae,z,prob=prob,varlevel=varlevel)
end

function getlatentvarsforz(vae,z,prob=true,varlevel=1)
	latent = vae.decoder[1:varlevel](z).data
	if prob == false
		latent = splitvars(latent)
	end
	latent
end

function splitvars(x)
	return Float64.(x .>= 0)
end

function samplefromtanh(x)
	ifelse(x < (2 * rand() - 1),0.0,1.0)
end

function samples(vae::VAE, nsamples)
   dz = size(vae.mu.W, 1)
	z = randn(dz, nsamples)
	s = vae.decoder(z).data
	copy(Transpose(float.(s .> 0.5)))
end

function trainvae!(vae, x; epochs = 1,learningrate = 0.001,batchsize = 10)

   g(x) = (h = vae.encoder(x); (vae.mu(h), vae.logsigma(h)))
   loss(x) = ((mu,logsigma) = g(x);
         (reconstructionloss(vae.decoder(normalz.(mu,logsigma)),x) -
               kullbackleiblerqp(mu,logsigma)) / batchsize)
   penloss(x) = -loss(x) + decoderpen(vae)

   n,p =size(x)
   optimizer = ADAM(learningrate)
   vaeparams = Flux.params(vae.encoder, vae.mu, vae.logsigma, vae.decoder)
   for i = 1:epochs
       data = [x[i,:]' for i in Iterators.partition(Random.shuffle(1:n), batchsize)]
       Flux.train!(penloss,vaeparams,zip(data),optimizer)
   end
   vae
end

end # module VAEs

max_stars_repo_name: stefan-m-lenz/dist-gen-comp

max_issues_repo_name: stefan-m-lenz/dist-gen-comp

max_forks_repo_name: stefan-m-lenz/dist-gen-comp


text:
# ==============================================================================
# Xavier Gandibleux - November 2021
#   Implemented in Julia 1.6

# ==============================================================================
# Parser of bi-01BKP instances generated by O. Perederieieva and X. Gandibleux
#   All coefficients P, W, and  ω are integer, all variables are binary

# ==============================================================================
# Datastructure of a multi-objective multi-dimensionnal KP with 0/1 variables
struct _MOMKP
    P  :: Matrix{Int} # profit of items for the objectives, k=1..p, j=1..n
    W  :: Matrix{Int} # weight of items for the constraints, i=1..m, j=1..n
    ω  :: Vector{Int} # capacity of knapsacks, i=1..m
end

# ==============================================================================
function readInstanceMOMKPformatPG(verbose::Bool, fname::String)

    f=open(fname)

    n = parse(Int, readline(f) )
    m = 2  # 2 constraints
    p = 2  # 2 objectives

    momkp = _MOMKP(zeros(Int,p,n),zeros(Int,m,n),zeros(Int,m))
    !verbose ? nothing : println("Instance of bi-01BKP : ",n," variables")

    momkp.ω[1] = parse(Int, readline(f) )
    momkp.ω[2] = parse(Int, readline(f) )
    !verbose ? nothing : println("ω[1]=",momkp.ω[1])
    !verbose ? nothing : println("ω[2]=",momkp.ω[2])

    # Extract the data of the MKP from the following lines ---------------------
    for j = 1:n
        momkp.P[1,j],
        momkp.P[2,j],
        momkp.W[1,j],
        momkp.W[2,j] = parse.(Int, split(readline(f)) )

        !verbose ? nothing :
          println("j=",j, "  |  ",
                  "p[1,",j,"]=",momkp.P[1,j], "   ",
                  "p[2,",j,"]=",momkp.P[2,j], "  |  ",
                  "w[1,",j,"]=",momkp.W[1,j], "   ",
                  "w[2,",j,"]=",momkp.W[2,j]
                 )
    end
    close(f)

    return momkp
end

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

# Example on how to use it :
# fname="set1/ZL28.DAT"; verbose = true
# instancePG = readInstanceMOMKPformatPG(verbose,fname)


text:
## Allows for lazy calculation (when requested)
## Also will have observer mixin in here

type LazyMixin
  calculated::Bool
  frozen::Bool
  observe::ObserverMixin

  LazyMixin() = new(false, false, ObserverMixin())
end

is_calculated{L <: LazyObject}(lazy::L) = lazy.lazyMixin.calculated
is_frozen{L <: LazyObject}(lazy::L) = lazy.lazyMixin.frozen

calculated!{L <: LazyObject}(lazy::L, setting::Bool = true) = lazy.lazyMixin.calculated = setting

function calculate!{L <: LazyObject}(lazy::L)
  if !is_calculated(lazy)
    calculated!(lazy)
    perform_calculations!(lazy)
  end

  return lazy
end

function recalculate!{L <: LazyObject}(lazy::L)
  calculated!(lazy, false)
  calculate!(lazy)

  return lazy
end

get_observers(lazy::LazyObject) = lazy.lazyMixin.observe.observers


text:
export UnionFind, addnode!, unify!, clusterize!, clusterid

@doc doc"""
Union-find algorithm.
"""
mutable struct UnionFind
    parents :: Vector{Int}
    weights :: Vector{Int}
    ids :: Vector{Int}
    nnodes :: Int
    nclusters :: Int
    UnionFind(n::Int=0) = new(collect(1:n),ones(Int,n),zeros(Int,n),n,n)
end

numnodes(u::UnionFind) = u.nnodes
numclusters(u::UnionFind) = u.nclusters

isroot(u::UnionFind, n::Integer) = u.parents[n] == n

@doc doc"""
    root!(u::UnionFind, n::Integer)

Returns the root node of the cluster where `n` belongs.
This may changes graph connection by "Path halving" method.
"""
function root!(u::UnionFind, n::Integer)
    @inbounds while !isroot(u,n)
        p = u.parents[n]
        u.parents[n] = u.parents[p]
        n = p
    end
    return n
end

@doc doc"""
    unify!(u, n1, n2)

Connects `n1` and `n2` nodes and returns the root.
"""
function unify!(u::UnionFind, n1::Integer, n2::Integer)
    r1 = root!(u,n1)
    r2 = root!(u,n2)
    @inbounds if r1 != r2
        u.nclusters -= 1
        if u.weights[r1] < u.weights[r2]
            r1,r2 = r2,r1
        end
        u.parents[r2] = r1
        u.weights[r2] += u.weights[r1]
    end
    return r1
end

@doc doc"""
    addnode!(u::UnionFind)

Adds a new node into `u` and returns the number of nodes including the added node.
"""
function addnode!(u::UnionFind)
    push!(u.parents,length(u.parents)+1)
    push!(u.weights,1)
    push!(u.ids, 0)
    u.nclusters += 1
    u.nnodes += 1
    return u.nnodes
end

@doc doc"""
    clusterize!(u::UnionFind)

Assigns cluster ID to each node and returns the number of clusters.
"""
function clusterize!(u::UnionFind)
    u.nclusters = 0
    @inbounds for i in 1:length(u.parents)
        if isroot(u,i)
            u.nclusters += 1
            u.ids[i] = u.nclusters
        end
    end
    @inbounds for i in 1:length(u.parents)
        u.ids[i] = u.ids[root!(u,i)]
    end
    return u.nclusters
end

@doc doc"""
    clusterid(u::UnionFind, i::Integer)

Returns the index of the cluster where `i` node belongs.
"""
clusterid(u::UnionFind, i::Integer) = u.ids[i]


text:
using CalibrationTests
using CalibrationErrors, Distributions
using Random, Test, Statistics

Random.seed!(1234)

# sample data
function generate_binary_data(nsamples)
    # generate predictions
    predictions = rand(Dirichlet(2, 1), nsamples)

    # generate labels
    labels_consistent = [rand() < predictions[1, i] ? 1 : 2 for i in 1:nsamples]
    labels_only_two = fill(2, nsamples)

    (predictions, labels_consistent), (predictions, labels_only_two)
end
const data_consistent, data_only_two = generate_binary_data(500)

# define matrix-valued kernel (using the mean total variation distance as bandwidth)
const kernel = UniformScalingKernel(ExponentialKernel(3))

@testset "Consistency test" begin
    println("consistency tests")

    # define estimators
    estimators = (ECE(UniformBinning(10)), BiasedSKCE(kernel),
                  LinearUnbiasedSKCE(kernel), QuadraticUnbiasedSKCE(kernel))

    for estimator in estimators
        test_consistent = ConsistencyTest(estimator, data_consistent)
        @test pvalue(test_consistent) > 0.7

        test_only_two = ConsistencyTest(estimator, data_only_two)
        @test pvalue(test_only_two) < 1e-6
    end
end

@testset "Distribution-free tests" begin
    println("distribution-free tests")

    # define estimators
    estimators = (BiasedSKCE(kernel), LinearUnbiasedSKCE(kernel),
                  QuadraticUnbiasedSKCE(kernel))

    for estimator in estimators
        test_consistent = DistributionFreeTest(estimator, data_consistent)
        @test pvalue(test_consistent) > 0.7

        test_only_two = ConsistencyTest(estimator, data_only_two)
        @test pvalue(test_only_two) < 1e-6
    end
end

@testset "Asymptotic linear test" begin
    println("asymptotic linear test")

    test_consistent = AsymptoticLinearTest(kernel, data_consistent)
    @test pvalue(test_consistent) > 0.7

    test_only_two = AsymptoticLinearTest(kernel, data_only_two)
    @test pvalue(test_only_two) < 1e-6
end

@testset "Asymptotic quadratic test" begin
    println("asymptotic quadratic test")

    test_consistent = AsymptoticQuadraticTest(kernel, data_consistent)
    @test pvalue(test_consistent) > 0.7

    test_only_two = AsymptoticQuadraticTest(kernel, data_only_two)
    @test pvalue(test_only_two) < 1e-6
end


text:
export relative_simple_extension, isprimitive, isprimitive_over

################################################################################
#
#  Relative extension
#
################################################################################

@doc Markdown.doc"""
    relative_simple_extension(K::NumField, k::NumField) -> NfRel

Given two fields $K\supset k$, it returns $K$ as a simple relative
extension $L$ of $k$ and an isomorphism $L \to K$.
"""
function relative_simple_extension(K::NumField, k::NumField)
  if _issubfield_in_tower(k, K)
    return relative_ext_in_tower(K, k)
  end
  error("Not yet implemented")
end

function relative_ext_in_tower(K::NumField, k::NumField)
  el = primitive_element(K, k)
  f = minpoly(el, k)
  L = number_field(f, cached = false, check = false)[1]
  mL = hom(L, K, el)
  return L, mL
end

function relative_simple_extension(K::AnticNumberField, k::AnticNumberField)
  fl, mp = issubfield(k, K)
  @assert fl
  return relative_simple_extension(mp)
end

function relative_simple_extension(m::NfToNfMor)
  k = domain(m)
  K = codomain(m)
  lf = factor(K.pol, k)
  rel_deg = divexact(degree(K), degree(k))
  pols = [f for (f, v) in lf if degree(f) == rel_deg]
  p = pols[1]
  if length(pols) > 1
    i = 2
    while !iszero(map_coefficients(m, p)(gen(K)))
      p = pols[i]
      i += 1
    end
  end
  L, b = number_field(p, cached = false, check = false)
  mp = hom(K, L, b, inverse = (image_primitive_element(m), gen(K)))
  return L, inv(mp)
end

################################################################################
#
#  Relative minpoly
#
################################################################################

function minpoly(a::NumFieldElem, K::NumField)
  if !_issubfield_in_tower(K, parent(a))
    error("Not yet implemented!")
  end
  f = minpoly(a)
  g = norm(f, K)
  if issquarefree(g)
    return g
  else
    h = gcd(g, derivative(g))
    return divexact(g, h)
  end
end

function minpoly(a::NumFieldElem, ::FlintRationalField)
  return absolute_minpoly(a)
end

################################################################################
#
#  IsSubfieldTower
#
################################################################################

function _issubfield_in_tower(k::NumField, K::NumField)
  if absolute_degree(k) == 1
    return true
  end
  found = true
  while k != base_field(K)
    K = base_field(K)
    if absolute_degree(k) > absolute_degree(K)
      found = false
      break
    end
  end
  return found
end

################################################################################
#
#  Relative primitive element
#
################################################################################

function primitive_element(K::NumField, k::NumField; check::Bool = true)
  if _issubfield_in_tower(k, K)
    return primitive_element_tower(K, k)
  else
    if check
      @assert issubfield(k, K)[1]
    end
    return primitive_element_generic(K, k)
  end
end

function primitive_element_generic(K::NumField, k::NumField)
  return absolute_primitive_element(K)
end

function primitive_element_tower(K::NumField, k::NumField)
  pK = primitive_element(K)
  if base_field(K) == k
    return pK
  end
  if isprimitive_over(pK, k)
    return pK
  end
  L = base_field(K)
  gL = primitive_element(L, k)
  pK += gL
  while !isprimitive_over(pK, k)
    pK += gL
  end
  return pK
end

function isprimitive_over(a::NumFieldElem, k::NumField)
  return divexact(absolute_degree(parent(a)), absolute_degree(k)) == degree(minpoly(a, k))
end


text:
import LightGraphs as LG

export min_fill, min_width, GraphModification
# Helper functions

struct GraphModification
    nb::Vector{Int}
    addedEdges::Vector{Tuple{Int, Int}}
end

# joins vertices together if they're already neighbors
@inline function joinVerts!(g::LG.AbstractGraph, v::Int)::GraphModification
    nb = copy(LG.neighbors(g, v))
    numberOfNeighbors = length(nb)
    addedEdges = []
    for i = 1:numberOfNeighbors-1
        for j = i+1:numberOfNeighbors
            if LG.add_edge!(g, nb[i], nb[j])
                push!(addedEdges, (nb[i], nb[j]))
            end
        end
    end
    GraphModification(nb, addedEdges)
end

# num of edges added if vertex removed
@inline function countEliminatedEdges(g::LG.AbstractGraph, v::Int)
    e = 0
    nb = LG.neighbors(g, v)
    numberOfNeighbors = length(nb)
    for i = 1:(numberOfNeighbors-1)
        for j = (i+1):numberOfNeighbors
            if !LG.has_edge(g, nb[i], nb[j])
                e+=1
            end
        end
    end
    e
end

@inline genNumberList(n::Int) = [i for i = 1:n]

# Heuristics Below

function min_fill(g::LG.AbstractGraph)
    nVertices = LG.nv(g)
    sorted = zeros(Int, nVertices)
    nums = genNumberList(nVertices)
    treeWidth = 0

    for i = 1:nVertices
        removal = argmin(map(v -> countEliminatedEdges(g,v), LG.vertices(g)))
        sorted[i] = nums[removal]

        x = nVertices + 1 - i
        if removal < x
            nums[removal] = x
        end

        treeWidth = max(treeWidth, LG.degree(g, removal))
        joinVerts!(g, removal)
        LG.rem_vertex!(g, removal)
    end

    (sorted, treeWidth)
end

function min_width(g::LG.AbstractGraph)
    nVertices = LG.nv(g)
    sorted = zeros(Int, nVertices)
    nums = genNumberList(nVertices)
    treeWidth = 0

    for i = 1:nVertices
        removal = argmin(LG.degree(g))
        sorted[i] = nums[removal]

        x = nVertices + 1 - i
        if removal < x
            nums[removal] = x
        end

        treeWidth = max(treeWidth, LG.degree(g, removal))
        joinVerts!(g, removal)
        LG.rem_vertex!(g, removal)
    end

    (sorted, treeWidth)
end

max_stars_repo_name: OleeL/treewidth-heuristics

max_issues_repo_name: OleeL/treewidth-heuristics

max_forks_repo_name: OleeL/treewidth-heuristics


text:
using Random
import StatsBase

"""
An empirical distribution, defined by a histogram which assigns probabilities to
a set of axis-coordinates in a (possibly high-dimensional) grid.
"""
struct EmpiricalDistribution{Dims, T}
    prob::Array{Float64, Dims}
    axes::SVector{Dims, Vector{T}}
end

Base.eltype(::Type{<:EmpiricalDistribution{Dims, T}}) where {Dims, T} = Vector{T}

function empirical_dist(prob::AbstractArray{Float64,N}, axes::Vararg{<:AbstractVector,N}) where {N}
    EmpiricalDistribution{N, eltype(axes[1])}(prob ./ sum(prob), (axes...,))
end

function Random.rand(rng::AbstractRNG, d::Random.SamplerTrivial{<:EmpiricalDistribution})
    dist = d[]
	x = StatsBase.sample(rng, StatsBase.weights(vec(dist.prob)))
	index = CartesianIndices(size(dist.prob))
	i = index[x]
	map(j->dist.axes[j][i[j]], 1:length(i))
end

function Distributions.logpdf(dist::EmpiricalDistribution, v::AbstractVector)
	indices = [findfirst(==(v[i]), dist.axes[i]) for i=eachindex(dist.axes)]
	if nothing ∈ indices
		-Inf64
	else
		log(dist.prob[indices...])
	end
end
Distributions.logpdf(dist::EmpiricalDistribution, v) = logpdf(dist, [v])

function Distributions.pdf(dist::EmpiricalDistribution, v::AbstractVector)
	indices = [findfirst(==(v[i]), dist.axes[i]) for i=eachindex(dist.axes)]
	if nothing ∈ indices
		0.0
	else
		dist.prob[indices...]
	end
end
Distributions.pdf(dist::EmpiricalDistribution, v) = pdf(dist, [v])

function marginalize(dist::EmpiricalDistribution, axis_index::Integer)
	new_prob = sum(dist.prob, dims=axis_index)
	new_prob = dropdims(new_prob, dims=axis_index)
	new_axes = deleteat(dist.axes, axis_index)
	EmpiricalDistribution(new_prob, new_axes)
end


text:
function stack_by(df, col, val, input, output)
    cols = unique(vcat(propertynames(df), output))

    df[!,:start] = (1:size(df)[1]) .+ 1
    df_split = df[df[:,col] .== val, :]
    df_split = edit_with(df_split, Rename.(input, output))

    sort!(unique!(df_split), :start)
    df_split[!, :stop] .= vcat(df_split[2:end, :start] .- 2, [size(df)[1]])

    [df[!, col] .= "" for col in output]
    [df[row[:start]:row[:stop], col] .= row[col] for row in eachrow(df_split) for col in output]

    df = df[df[:,output[1]] .!= "", :]
    return df[:, cols]
end

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# SCALE NAICS
filename = files_map[occursin.(joinpath("scale", "naics"), files_map)]

if length(filename) > 0
    println(string("Standardizing ", filename[1]))
    y = read_file(filename[1])
    df = read_file(y["PathIn"], y["CSVInput"])

    col = :naics_level
    input = [:naics_code, :naics_desc]
    vals = ["sector", "subsector", "industry group", "naics industry", "national industry"]
    outputs = [Symbol.(replace(val, " " => "_"), ["_code", "_desc"]) for val in vals]

    [global df = stack_by(df, col, val, input, output)
        for (val, output) in collect(zip(vals, outputs))[1:end-1]]
    df = edit_with(df, Rename.(input, outputs[end]))

    global df = edit_with(df, y["Order"])
    CSV.write(joinpath(SLIDE_DIR, y["PathOut"]...), df)
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


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

struct HostsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_apiV1FleetHostsCountGet(_api::HostsApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Any, "/api/v1/fleet/hosts/count", ["bearerAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? [] : [_mediaType])
    return _ctx
end

@doc raw"""Hosts Count

Params:

Return: Any
"""
function apiV1FleetHostsCountGet(_api::HostsApi; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsCountGet(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function apiV1FleetHostsCountGet(_api::HostsApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsCountGet(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_apiV1FleetHostsDeletePost(_api::HostsApi; body=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/api/v1/fleet/hosts/delete", ["bearerAuth"], body)
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

@doc raw"""Bulk Delete Hosts

Params:
- body::Any

Return: Nothing
"""
function apiV1FleetHostsDeletePost(_api::HostsApi; body=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsDeletePost(_api; body=body, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function apiV1FleetHostsDeletePost(_api::HostsApi, response_stream::Channel; body=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsDeletePost(_api; body=body, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_apiV1FleetHostsGet(_api::HostsApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Any, "/api/v1/fleet/hosts", ["bearerAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? [] : [_mediaType])
    return _ctx
end

@doc raw"""List Hosts

Params:

Return: Any
"""
function apiV1FleetHostsGet(_api::HostsApi; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsGet(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function apiV1FleetHostsGet(_api::HostsApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsGet(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_apiV1FleetHostsHostIdGet(_api::HostsApi, host_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Any, "/api/v1/fleet/hosts/{host_id}", ["bearerAuth"])
    Swagger.set_param(_ctx.path, "host_id", host_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? [] : [_mediaType])
    return _ctx
end

@doc raw"""Get Host By ID

Params:
- host_id::String (required)

Return: Any
"""
function apiV1FleetHostsHostIdGet(_api::HostsApi, host_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsHostIdGet(_api, host_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function apiV1FleetHostsHostIdGet(_api::HostsApi, response_stream::Channel, host_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsHostIdGet(_api, host_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_apiV1FleetHostsDeviceMappingGet(_api::HostsApi, host_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Any, "/api/v1/fleet/hosts/{host_id}/device_mapping", ["bearerAuth"])
    Swagger.set_param(_ctx.path, "host_id", host_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? [] : [_mediaType])
    return _ctx
end

@doc raw"""Get Host Device Mapping

Params:
- host_id::String (required)

Return: Any
"""
function apiV1FleetHostsDeviceMappingGet(_api::HostsApi, host_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsDeviceMappingGet(_api, host_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function apiV1FleetHostsDeviceMappingGet(_api::HostsApi, response_stream::Channel, host_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetHostsDeviceMappingGet(_api, host_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export apiV1FleetHostsCountGet, apiV1FleetHostsDeletePost, apiV1FleetHostsGet, apiV1FleetHostsHostIdGet, apiV1FleetHostsDeviceMappingGet


text:
import Base.length
import Base.show
import Base.size

flatten(arr) = reshape(arr, length(arr))

typealias Volume{T} Array{T, 3}

type Vol
    w::Volume{Float64}
    dw::Volume{Float64}

    function Vol(dimensions, value::Float64)
        w = Array(Float64, dimensions)
        dw = Array(Float64, dimensions)
        new(ones(dimensions) * value, zeros(dimensions))
    end

    Vol(w::Array{Float64, 3}, dw::Array{Float64, 3}) = new(w, dw)

    function Vol(dimensions)
        (x, y, depth) = dimensions
        scale = sqrt(1.0 / (x * y * depth))
        new(randn(x, y, depth) * scale, zeros(x, y, depth))
    end
end

function show(io::IO, v::Vol)
    println("Vol")
    println(io, "w: ", v.w)
    println(io, "dw: ", v.dw)
end

function length(v::Vol)
    @assert size(v.w) == size(v.dw)
    length(v.w)
end

function size(v::Vol)
    @assert size(v.w) == size(v.dw)
    size(v.w)
end


text:
include("thermalizemetropolis.jl")
include("finishmetropolis.jl")


text:
using .APIUtils

using CEnum: @cenum

using Printf

import Libdl

using LazyArtifacts

# low-level wrappers
const CUdeviceptr = CuPtr{Cvoid}
const CUarray = CuArrayPtr{Cvoid}
const GLuint = Cuint    # FIXME: get these from somewhere
const GLenum = Cuint
include("libcuda_common.jl")
include("error.jl")
include("libcuda.jl")
include("libcuda_deprecated.jl")

@noinline function find_libcuda()
    # NOTE: we can't use function wrappers from CUDAdrv.jl here,
    #       as those depend on `libcuda()` which is not set yet.

    # find the local system library
    system_driver = if Sys.iswindows()
        Libdl.find_library("nvcuda")
    else
        Libdl.find_library(["libcuda.so.1", "libcuda.so"])
    end
    if system_driver == ""
        if Sys.iswindows()
            error("""Could not find the CUDA driver library. Please make sure you have installed the NVIDIA driver for your GPU.
                     If you're sure it's installed, look for `nvcuda.dll` in your system and make sure it's discoverable by the linker.
                     Typically, that involves adding an entry to PATH.""")
        else
            error("""Could not find the CUDA driver library. Please make sure you have installed the NVIDIA driver for your GPU.
                     If you're sure it's installed, look for `libcuda.so` in your system and make sure it's discoverable by the linker.
                     Typically, that involves an entry in '/etc/ld.so.conf', or setting LD_LIBRARY_PATH.""")
        end
    end

    # query the system driver version
    function get_version(driver)
        library_handle = Libdl.dlopen(driver)
        try
            function_handle = Libdl.dlsym(library_handle, "cuDriverGetVersion")
            version_ref = Ref{Cint}()
            @check ccall(function_handle, CUresult, (Ptr{Cint},), version_ref)
            major, ver = divrem(version_ref[], 1000)
            minor, patch = divrem(ver, 10)
            return VersionNumber(major, minor, patch)
        finally
            Libdl.dlclose(library_handle)
        end
    end

    # XXX: apparently cuDriverGetVersion can be used before cuInit,
    #      despite the docs stating "any function [...] will return
    #      CUDA_ERROR_NOT_INITIALIZED"; is this a recent change?
    _system_version[] = get_version(system_driver)

    # check if we managed to unload the system driver.
    # if we didn't, we can't consider a forward compatible library because that would
    # risk having multiple copies of libcuda.so loaded (also see NVIDIA bug #3418723)
    system_driver_loaded = Libdl.dlopen(system_driver, Libdl.RTLD_NOLOAD; throw_error=false) !== nothing
    if system_driver_loaded
        @debug "Could not unload the system CUDA library; this will prevent use of the forward-compatible package"
        return system_driver
    end

    # check if this process is hooked by CUDA's injection libraries, which prevents
    # unloading libcuda after dlopening. this is problematic, because we might want to
    # after loading a forwards-compatible libcuda and realizing we can't use it. without
    # being able to unload the library, we'd run into issues (see NVIDIA bug #3418723)
    hooked = haskey(ENV, "CUDA_INJECTION64_PATH")
    if hooked
        @debug "Running under CUDA injection tools; this will prevent use of the forward-compatible package"
        return system_driver
    end

    # if we're using an older driver; consider using forward compatibility
    function do_init(driver)
        library_handle = Libdl.dlopen(driver)
        try
            function_handle = Libdl.dlsym(library_handle, "cuInit")
            @check ccall(function_handle, CUresult, (UInt32,), 0)
        finally
            Libdl.dlclose(library_handle)
        end
        return
    end
    if getenv("JULIA_CUDA_USE_COMPAT", !hooked && !system_driver_loaded) && _system_version[] < v"11.6"
        artifact = try
            # work around artifact"" eagerness on unsupported platforms by passing a variable
            f = id -> @artifact_str(id)
            f("CUDA_compat")
        catch ex
            @debug "Could not download forward compatibility package" exception=(ex,catch_backtrace())
            nothing
        end

        if artifact !== nothing
            # TODO: do we need to dlopen the JIT compiler library for it to be discoverable?
            #       doesn't that clash with a system one if compat cuInit fails? or should
            #       we load it _after_ the compat driver initialization succeeds?
            #compat_compiler = joinpath(artifact, "lib", "libnvidia-ptxjitcompiler.so")
            #Libdl.dlopen(compat_compiler)

            compat_driver = joinpath(artifact, "lib", "libcuda.so")
            try
                # just calling cuDriverGetVersion doesn't trigger errors,
                # so actually perform driver intialization here already.
                # XXX: reuse __init_driver__, which detects e.g. running under rr
                do_init(compat_driver)
                return compat_driver
            catch ex
                @debug "Could not use forward compatibility package" exception=(ex,catch_backtrace())

                # see comment above about unloading the system driver
                compat_driver_loaded = Libdl.dlopen(compat_driver, Libdl.RTLD_NOLOAD; throw_error=false) !== nothing
                if compat_driver_loaded
                    error("""Could not unload the forward compatible CUDA driver library.

                                This is probably caused by running Julia under a tool that hooks CUDA API calls.
                                In that case, prevent CUDA.jl from loading multiple drivers by setting JULIA_CUDA_USE_COMPAT=false in your environment.""")
                end
            end
        end
    end

    return system_driver
end

const _libcuda = Ref{Union{String,Exception}}()
function libcuda()
    # libcuda discovery has not been triggered yet
    if !isassigned(_libcuda)
        try
            _libcuda[] = find_libcuda()
            __init_driver__()
        catch err
            _libcuda[] = err
            rethrow()
        end
    end

    # libcuda was not found, or initialization failed.
    # since we may have cached a value of libcuda (through ccalls from __init_driver__)
    # we cannot re-initialize in this session, so re-throw the original error again.
    if _libcuda[] isa Exception
        throw(_libcuda[])
    end

    _libcuda[]::String
end

# high-level wrappers
include("types.jl")
include("version.jl")
include("devices.jl")
include("context.jl")
include("stream.jl")
include("pool.jl")
include("memory.jl")
include("module.jl")
include("events.jl")
include("execution.jl")
include("profile.jl")
include("occupancy.jl")
include("graph.jl")

# global state (CUDA.jl's driver wrappers behave like CUDA's runtime library)
include("state.jl")


text:
function func3(entity_ptr, r_::Ptr{T}, h_::Ptr{T}, m_ptr, u_::Ptr{T}) where T
# an anisotropic function
# populates h with the desired mesh size in all three dimension

#  println("entered func3")
  # load arrays from C pointers
  r = pointer_to_array(r_, (3,3))  # REMEMBER: r is transposed when it is passed back to C
#  println("r = \n", r)
  h = pointer_to_array(h_, 3)
#  println("h = \n", h)
  u = pointer_to_array(u_, 9)
#  println("u = \n", u)

  # get vertex coords
  coords = zeros(3,1)
  getVertCoords(entity_ptr, coords, 3, 1)
  x = coords[1]
  y = coords[2]
  z = coords[3]

  # use xyz frame
  r[1,1] = 1.0
  r[2,2] = 1.0
  r[3,3] = 1.0

#  println("in julia, r = ", r)

  # calculate derivative of rho here

  drho_dx = smoothHeavisideder(x)
  h[1] = abs(0.5/drho_dx)  # make mesh size proportional to 1/drho/dx (larger gradient -> smaller mesh)
     # 2 is an emperical things with weird units

  ubound = 0.50
  lbound = 0.001
  # upper and lower bounds
  if h[1] < lbound
    h[1] = lbound
  elseif h[1] > ubound
    h[1] = ubound
  end

#  h[1] = 0.5
  h[2] = h[1]
  h[3] = 2.0

  println("x = ", x, " h = ", h, " drho_dx = ", drho_dx)
return nothing
end

# smooth heaviside function
function smoothHeavisideder(x)
# calculate the value of the smooth heaviside function at a location x
# x0 is specified within this function

  x0 = 0
  L = 5
  k = 5

#  return (L/(1 + e^(-k*(x-x0))))
  return L*(2*k*e^(-2*k*x))/(e^(-2*k*x) +1 )^2
end

function shockRefine(entity_ptr, r_::Ptr{T}, h_::Ptr{T}, m_ptr, f_ptr) where T
# an anisotropic function
# populates h with the desired mesh size in all three dimension
# f_ptr is a pointer to a solution field (apf::Field)

#  println("entered func3")
  # load arrays from C pointers
  r = pointer_to_array(r_, (3,3))  # REMEMBER: r is transposed when it is passed back to C
#  println("r = \n", r)
  h = pointer_to_array(h_, 3)
#  println("h = \n", h)

  # get vertex coords
  coords = zeros(3,1)
  getVertCoords(entity_ptr, coords, 3, 1)
  x = coords[1]
  y = coords[2]
  z = coords[3]

  # use xyz frame
  r[1,1] = 1.0
  r[2,2] = 1.0
  r[3,3] = 1.0

#  println("in julia, r = ", r)

  u_node = zeros(4)
  retrieveNodeSolution(f_ptr, entity_ptr, u_node)
  # calculate derivative of rho here

  drho_dx = u_node[1]  # rho  value
#  drho_dx = smoothHeavisideder(x)

  h[1] = abs(0.5/drho_dx)  # make mesh size proportional to 1/drho/dx (larger gradient -> smaller mesh)
     # 2 is an emperical things with weird units

  ubound = 0.50
  lbound = 0.001
  # upper and lower bounds
  if h[1] < lbound
    h[1] = lbound
  elseif h[1] > ubound
    h[1] = ubound
  end

#  h[1] = 0.5
  h[2] = h[1]
  h[3] = 2.0

  println("x = ", x, " h = ", h, " drho_dx = ", drho_dx)
return nothing
end

function shockRefine2(entity_ptr, r_::Ptr{T}, h_::Ptr{T}, m_ptr, f_ptr, operator_ptr) where T
# an anisotropic function
# populates h with the desired mesh size in all three dimension
# entity_ptr is a pointer to a vertex
# r_ptr is a pointer to an array of doubles that specifiy the orthonormal
# coordinate system to use
# h_ptr is the desired mesh size in each direction of r
# m_ptr is the mesh pointer
# f_ptr is a pointer to a solution field (apf::Field)
# operator_ptr = pointer to SBP operator

#  println("entered func3")

  # load arrays from C pointers
  r = pointer_to_array(r_, (3,3))  # REMEMBER: r is transposed when it is passed back to C
#  println("r = \n", r)

  h = pointer_to_array(h_, 3)
#  println("h = \n", h)

#  println("operator_ptr = ", operator_ptr)
  sbp = unsafe_pointer_to_objref(operator_ptr)

  # get vertex coords
  coords = zeros(3,1)
  getVertCoords(entity_ptr, coords, 3, 1)
  x = coords[1]
  y = coords[2]
  z = coords[3]

  # use xyz frame
  r[1,1] = 1.0
  r[2,2] = 1.0
  r[3,3] = 1.0

#  println("in julia, r = ", r)

  # get the elements this vertex is part of
  num_elements = countAdjacent(m_ptr, entity_ptr, 2)  # get adjacnet elements
  elements = getAdjacent(num_elements)

  # choose first element (arbitrary)
  el_ptr = elements[1]
  verts, tmp = getDownward(m_ptr, el_ptr, 0)

  # get coordinates
  coords = zeros(3,3, 1)
  sub_coords = sub(coords, :, :, 1)
  getFaceCoords(el_ptr, sub_coords, 3, 3)

  # get solution at vertices, find out the local index of the vertex
  u_vals = zeros(4,3, 1)
  vert_index = 0
  for j=1:3  # get solution value of each node, figure out which vertex is original
    subarray = sub(u_vals, :, j, 1)
    retrieveNodeSolution(f_ptr, verts[j], subarray)

    if (entity_ptr == verts[j])
      vert_index = j
    end
  end

  # get mapping jacobian
  dxi_dx = zeros(2,2,3,1)
  jac = zeros(3,1)
  mappingjacobian!(sbp, coords[1:2,:,:], dxi_dx, jac)

  # differentiate (mapping jacobian is constant for an element)
  u_dxidx = u_vals*dxi_dx[1,1,1,1]
  u_detadx = u_vals*dxi_dx[2,1,1,1]
  res = zeros(size(u_vals))

  differentiate!(sbp, 1, u_dxidx, res)
  differentiate!(sbp, 2, u_detadx, res)

  drho_dx = res[1, vert_index, 1]*jac[vert_index,1]
#  println("drho_dx = ", drho_dx)

  # make mesh size proportional to 1/ drhodx
  h[1] = abs(0.25/drho_dx)

  ubound = 0.50
  lbound = 0.0001
  # upper and lower bounds
  if h[1] < lbound
    h[1] = lbound
  elseif h[1] > ubound
    h[1] = ubound
  end

#  h[1] = 0.5
  h[2] = h[1]  # make mesh square
  h[3] = 2.0

#  println("x = ", x, " h = ", h, " drho_dx = ", drho_dx)
#  print("\n")
return nothing
end


text:
#
# Methods for orbit states and state vectors.
#

keplerian(state::CartesianStateVector, μ) = keplerian(position(state), velocity(state), μ)
keplerian(orbit::CartesianR2BPOrbit) = keplerian(state(orbit), massparameter(system(orbit)))
keplerian(orbit::KeplerianR2BPOrbit) = orbit

cartesian(orbit::CartesianR2BPOrbit) = orbit
cartesian(orbit::KeplerianR2BPOrbit) = cartesian(state(orbit), massparameter(system(orbit)))
cartesian(state::KeplerianState, μ) = cartesian(
    eccentricity(state),
    semimajor_axis(state),
    inclination(state),
    RAAN(state),
    argument_of_periapsis(state),
    true_anomoly(state),
    μ
)

Base.position(state::CartesianStateVector) = OrbitalStates.get_r(state) * lengthunit(state)
Base.position(orbit::KeplerianR2BPOrbit) = position(CartesianState(cartesian(state(orbit), massparameter(state(orbit)))))
Base.position(orbit::Orbit) = position(state(orbit))

distance(state::CartesianStateVector) = norm(position(state))
distance(orbit::Orbit) = norm(position(orbit))

velocity(state::CartesianStateVector) = OrbitalStates.get_v(state) * velocityunit(state)
velocity(orbit::Orbit) = velocity(state(orbit))

speed(state::CartesianStateVector) = norm(velocity(state))
speed(orbit::KeplerianR2BPOrbit) = speed(semi_parameter(semimajor_axis(state(orbit)), eccentricity(state(orbit))), semimajor_axis(state(orbit)), massparameter(system(orbit)))
speed(orbit::Orbit) = norm(velocity(state(orbit)))

eccentricity(state::KeplerianState) = OrbitalStates.get_e(state)
eccentricity(orbit::KeplerianR2BPOrbit) = eccentricity(state(orbit))
eccentricity(orbit::CartesianR2BPOrbit) = eccentricity(position(state(orbit)), velocity(state(orbit)), massparameter(system(orbit)))

semimajor_axis(state::KeplerianState) = OrbitalStates.get_a(state) * lengthunit(state)
semimajor_axis(orbit::CartesianR2BPOrbit) = semimajor_axis(distance(state(orbit)), speed(state(orbit)), massparameter(system(orbit)))
semimajor_axis(orbit::KeplerianR2BPOrbit) = semimajor_axis(state(orbit))

inclination(state::KeplerianState) = OrbitalStates.get_i(state) * angularunit(state)
inclination(orbit::KeplerianR2BPOrbit) = inclination(state(orbit))

RAAN(state::KeplerianState) = OrbitalStates.get_Ω(state) * angularunit(state)
RAAN(orbit::KeplerianR2BPOrbit) = RAAN(state(orbit))

argument_of_periapsis(state::KeplerianState) = OrbitalStates.get_ω(state) * angularunit(state)
argument_of_periapsis(orbit::KeplerianR2BPOrbit) = argument_of_periapsis(state(orbit))

true_anomoly(state::KeplerianState) = OrbitalStates.get_ν(state) * angularunit(state)
true_anomoly(orbit::KeplerianR2BPOrbit) = true_anomoly(state(orbit))
true_anomoly(orbit::CartesianR2BPOrbit) = true_anomoly(distance(orbit), specific_angular_momentum(orbit), eccentricity(orbit), massparameter(system(orbit)))

conic(orbit::R2BPOrbit) = conic(eccentricity(orbit))

specific_angular_momentum_vector(orbit::CartesianR2BPOrbit) = specific_angular_momentum_vector(position(state(orbit)), velocity(state(orbit)))

specific_angular_momentum(orbit::CartesianR2BPOrbit) = specific_angular_momentum(position(state(orbit)), velocity(state(orbit)))

specific_energy(orbit::CartesianR2BPOrbit) = specific_energy(distance(state(orbit)), speed(state(orbit)), massparameter(system(orbit)))
specific_energy(orbit::KeplerianR2BPOrbit) = specific_energy(semimajor_axis(orbit), massparameter(system(orbit)))

C3(orbit::R2BPOrbit) = C3(distance(orbit), speed(orbit), massparameter(system(orbit)))

v_infinity(orbit::R2BPOrbit) = v_infinity(distance(orbit), speed(orbit), massparameter(system(orbit)))

eccentricity_vector(orbit::CartesianR2BPOrbit) = eccentricity_vector(position(state(orbit)), velocity(state(orbit)), massparameter(system(orbit)))

semi_parameter(orbit::R2BPOrbit) = semi_parameter(semimajor_axis(orbit), eccentricity(orbit))

periapsis_radius(orbit::R2BPOrbit) = periapsis_radius(semimajor_axis(orbit), eccentricity(orbit))

apoapsis_radius(orbit::R2BPOrbit) = apoapsis_radius(semimajor_axis(orbit), eccentricity(orbit))

periapsis_velocity(orbit::R2BPOrbit) = speed(periapsis_radius(orbit), semimajor_axis(orbit), massparameter(system(orbit)))

apoapsis_velocity(orbit::R2BPOrbit) = speed(apoapsis_radius(orbit), semimajor_axis(orbit), massparameter(system(orbit)))

period(orbit::R2BPOrbit) = conic(orbit) ∈ (Circular, Elliptical) ? period(semimajor_axis(orbit), massparameter(system(orbit))) : conic(orbit) == Parabolic ? Inf * timeunit(orbit) : NaN * timeunit(orbit)

mean_motion(orbit::R2BPOrbit) = mean_motion(semimajor_axis(orbit), massparameter(system(orbit)))

function mean_motion_vector(orbit::R2BPOrbit)
#   î = SVector{3, Float64}([1, 0, 0])
#   ĵ = SVector{3, Float64}([0, 1, 0])
    k̂ = SVector{3}(0, 0, 1)
    return k̂ × specific_angular_momentum_vector(orbit)
end

time_since_periapsis(orbit::R2BPOrbit) = time_since_periapsis(mean_motion(orbit), eccentricity(orbit), eccentric_anomoly(orbit))

specific_potential_energy(orbit::CartesianR2BPOrbit) = specific_potential_energy(distance(orbit.state), massparameter(system(orbit)))


text:
using Test, YaoBlocks, YaoArrayRegister, LuxurySparse

@testset "test constructor" for T in [Float16, Float32, Float64]
    @test PhaseGate(0.1) isa PrimitiveBlock{1}
    @test_throws TypeError PhaseGate{Complex{T}} # will not accept non-real type
    @test phase(T(0.1)) isa PrimitiveBlock{1}
    @test phase(1) isa PhaseGate{Float64} # default we convert to float64
end

@testset "test copy" begin
    g = PhaseGate(0.1)
    cg = copy(g)
    cg.theta = 0.2
    @test g.theta == 0.1
end

@testset "test matrix" begin
    g = PhaseGate{Float64}(pi)
    @test mat(g) ≈ exp(im * pi) * IMatrix(2)
end

@testset "test operations" begin
    g = PhaseGate{Float64}(pi)
    reg = rand_state(1)
    @test mat(g) * state(reg) ≈ state(apply!(reg, g))

    @test (PhaseGate(2.0) == PhaseGate(2.0)) == true
end

@testset "test dispatch & parameters" begin
    @test nparameters(phase(0.1)) == 1
    @test adjoint(phase(0.1)) == phase(-0.1)
    @test dispatch!(phase(0.1), 0.3) == phase(0.3)

    @testset "test $op" for op in [+, -, *, /]
        @test dispatch!(op, phase(0.1), π) == phase(op(0.1, π))
    end

    @test_throws AssertionError dispatch!(phase(0.1), (0.2, 0.3))
end

@testset "properties" begin
    g = PhaseGate{Float64}(0.1)
    @test nqubits(g) == 1
    @test isreflexive(g) == false
    @test isunitary(g) == true
    @test ishermitian(g) == false
end


text:
module GadflyPlots

import Gadfly
import Weave
try
    import Cairo
catch
    @warn("Cairo.jl is required to be installed to generate raster images")
end

Gadfly.set_default_plot_format(:svg)

function Base.display(report::Weave.Report, m::MIME"image/png", p::Gadfly.Plot)
    display(report, MIME("image/svg+xml"), p)
end

#Gadfly doesn't call the default display methods, this catches
#all Gadfly plots
function Base.display(report::Weave.Report, m::MIME"image/svg+xml", p::Gadfly.Plot)

    chunk = report.cur_chunk

    w = chunk.options[:fig_width]Gadfly.inch
    h = chunk.options[:fig_height]Gadfly.inch
    format = chunk.options[:fig_ext]
    dpi = chunk.options[:dpi]

    full_name, rel_name = Weave.get_figname(report, chunk, ext = format)

    push!(report.figures, rel_name)
    report.fignum += 1

    if format == ".svg"
        Gadfly.draw(Gadfly.SVG(full_name, w, h), p)
    elseif format == ".js.svg"
        Gadfly.draw(Gadfly.SVGJS(full_name, w, h), p)
    elseif format == ".png"
        Gadfly.draw(Gadfly.PNG(full_name, w, h, dpi=dpi), p)
    elseif format == ".pdf"
        Gadfly.draw(Gadfly.PDF(full_name, w, h), p)
    elseif format == ".ps"
        Gadfly.draw(Gadfly.PS(full_name, w, h), p)
    elseif format == ".tex"
        Gadfly.draw(Gadfly.PGF(full_name, w, h, true ), p)
    else
        @warn("Can't save figure. Unsupported format, $format")
    end
end
end


text:
module LintExample2

#submodule
module LintExSubMod1
subfoo(x) = println(x)
end

module LintExSubMod2
using ..LintExSubMod1
subfar(x) = subfoo(x)
end

using LintExample

export foo, bar

@fancyfuncgen(foo)
@fancyfuncgen(bar)

end


text:
#=

Implements iteration on the Bellman equations to solve the McCall growth model

=#

using Distributions

# A default utility function

function u(c, sigma)
    if c > 0
        return (c^(1 - sigma) - 1) / (1 - sigma)
    else
        return -10e6
    end
end

# default wage vector with probabilities

const n = 60                                   # n possible outcomes for wage
const default_w_vec = linspace(10, 20, n)   # wages between 10 and 20
const a, b = 600, 400                          # shape parameters
const dist = BetaBinomial(n-1, a, b)
const default_p_vec = pdf(dist)

type McCallModel
    alpha::Float64        # Job separation rate
    beta::Float64         # Discount rate
    gamma::Float64        # Job offer rate
    c::Float64            # Unemployment compensation
    sigma::Float64        # Utility parameter
    w_vec::Vector{Float64} # Possible wage values
    p_vec::Vector{Float64} # Probabilities over w_vec

    function McCallModel(alpha=0.2,
                         beta=0.98,
                         gamma=0.7,
                         c=6.0,
                         sigma=2.0,
                         w_vec=default_w_vec,
                         p_vec=default_p_vec)

        return new(alpha, beta, gamma, c, sigma, w_vec, p_vec)
    end
end

"""
A function to update the Bellman equations.  Note that V_new is modified in
place (i.e, modified by this function).  The new value of U is returned.

"""
function update_bellman!(mcm, V, V_new, U)
    # Simplify notation
    alpha, beta, sigma, c, gamma = mcm.alpha, mcm.beta, mcm.sigma, mcm.c, mcm.gamma

    for (w_idx, w) in enumerate(mcm.w_vec)
        # w_idx indexes the vector of possible wages
        V_new[w_idx] = u(w, sigma) + beta * ((1 - alpha) * V[w_idx] + alpha * U)
    end

    U_new = u(c, sigma) + beta * (1 - gamma) * U +
                    beta * gamma * sum(max(U, V) .* mcm.p_vec)

    return U_new
end

function solve_mccall_model(mcm; tol::Float64=1e-5, max_iter::Int=2000)

    V = ones(length(mcm.w_vec))  # Initial guess of V
    V_new = similar(V)           # To store updates to V
    U = 1.0                        # Initial guess of U
    i = 0
    error = tol + 1

    while error > tol && i < max_iter
        U_new = update_bellman!(mcm, V, V_new, U)
        error_1 = maximum(abs(V_new - V))
        error_2 = abs(U_new - U)
        error = max(error_1, error_2)
        V[:] = V_new
        U = U_new
        i += 1
    end

    return V, U
end


text:
using Oceananigans.Fields: OneField
using Oceananigans.Grids: architecture
using Oceananigans.Architectures: arch_array
import Oceananigans.Fields: condition_operand, conditional_length, set!

# For conditional reductions such as mean(u * v, condition = u .> 0))

struct ConditionalOperation{LX, LY, LZ, O, F, G, C, M, T} <: AbstractOperation{LX, LY, LZ, G, T}
    operand :: O
    func :: F
    grid :: G
    condition :: C
    mask :: M

     function ConditionalOperation{LX, LY, LZ}(operand::O, func::F, grid::G, condition::C, mask::M) where {LX, LY, LZ, O, F, G, C, M}
         T = eltype(operand)
         return new{LX, LY, LZ, O, F, G, C, M, T}(operand, func, grid, condition, mask)
     end
end

"""
    ConditionalOperation{LX, LY, LZ}(operand, func, grid, condition, mask)

Returns an abstract representation of a masking procedure applied when `condition` is satisfied on a field
described by `func(operand)`.

Positional arguments
====================

- `operand`: The `AbstractField` to be masked (it must have a `grid` property!)

Keyword arguments
=================

- `func`: A unary transformation applied element-wise to the field `operand` at locations where `condition == true`. Default is `identity`

- `condition`: either a function of `(i, j, k, grid, operand)` returning a Boolean,
               or a 3-dimensional Boolean `AbstractArray`. At locations where `condition == false`,
               operand will be masked by `mask`

- `mask`: the scalar mask

`condition_operand` is a conveniece function used to construct a `ConditionalOperation`

`condition_operand(func::Function, operand::AbstractField, condition, mask) = ConditionalOperation(operand; func, condition, mask)`

Example
=======

```jldoctest
julia> using Oceananigans

julia> using Oceananigans.Fields: condition_operand

julia> c = CenterField(RectilinearGrid(size=(2, 1, 1), extent=(1, 1, 1)));

julia> f(i, j, k, grid, c) = i < 2; d = condition_operand(cos, c, f, 10)
ConditionalOperation at (Center, Center, Center)
├── operand: 2×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── func: typeof(cos)
├── condition: typeof(f)
└── mask: 10

julia> d[1, 1, 1]
1.0

julia> d[2, 1, 1]
10
```
"""
function ConditionalOperation(operand::AbstractField; func = identity, condition = nothing, mask = 0)
    grid = operand.grid
    LX, LY, LZ = location(operand)
    return ConditionalOperation{LX, LY, LZ}(operand, func, grid, condition, mask)
end

function ConditionalOperation(c::ConditionalOperation; func = nothing, condition = nothing, mask = nothing)
    LX, LY, LZ = location(c)
    return ConditionalOperation{LX, LY, LZ}(c.operand, c.func, c.grid, c.condition, c.mask)
end

@inline condition_operand(func::Function, operand::AbstractField, condition, mask) = ConditionalOperation(operand; func, condition, mask)
@inline condition_operand(func::Function, operand::AbstractField, ::Nothing, mask) = ConditionalOperation(operand; func, condition = truefunc, mask)
@inline function condition_operand(func::Function, operand::AbstractField, condition::AbstractArray, mask)
    condition = arch_array(architecture(operand.grid), condition)
    return ConditionalOperation(operand; func, condition, mask)
end

@inline condition_operand(func::typeof(identity), c::ConditionalOperation, ::Nothing, mask) = ConditionalOperation(c; mask)
@inline condition_operand(func::Function, c::ConditionalOperation, ::Nothing, mask) = ConditionalOperation(c; func, mask)

@inline truefunc(args...) = true

@inline condition_onefield(c::ConditionalOperation{LX, LY, LZ}, mask) where {LX, LY, LZ} =
                              ConditionalOperation{LX, LY, LZ}(OneField(Int), identity, c.grid, c.condition, mask)

@inline conditional_length(c::ConditionalOperation)       = sum(condition_onefield(c, 0))
@inline conditional_length(c::ConditionalOperation, dims) = sum(condition_onefield(c, 0); dims = dims)

Adapt.adapt_structure(to, c::ConditionalOperation{LX, LY, LZ}) where {LX, LY, LZ} =
            ConditionalOperation{LX, LY, LZ}(adapt(to, c.operand),
                                     adapt(to, c.func),
                                     adapt(to, c.grid),
                                     adapt(to, c.condition),
                                     adapt(to, c.mask))

@inline function Base.getindex(c::ConditionalOperation, i, j, k)
    return ifelse(get_condition(c.condition, i, j, k, c.grid, c),
                  c.func(getindex(c.operand, i, j, k)),
                  c.mask)
end

@inline concretize_condition!(c::ConditionalOperation) = set!(c.operand, c)

function concretize_condition(c::ConditionalOperation)
    f = similar(c.operand)
    set!(f, c)
    return f
end

@inline get_condition(condition, i, j, k, grid, args...)                = condition(i, j, k, grid, args...)
@inline get_condition(condition::AbstractArray, i, j, k, grid, args...) = @inbounds condition[i, j, k]

Base.summary(c::ConditionalOperation) = string("ConditionalOperation of ", summary(c.operand), " with condition ", summary(c.condition))

Base.show(io::IO, operation::ConditionalOperation) =
    print(io,
          "ConditionalOperation at $(location(operation))", '\n',
          "├── operand: ", summary(operation.operand), '\n',
          "├── grid: ", summary(operation.grid), '\n',
          "├── func: ", summary(operation.func), '\n',
          "├── condition: ", summary(operation.condition), '\n',
          "└── mask: ", operation.mask)


text:
module DoublePendulumWithDampers

using ModiaLang

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

import Modia3D
using  Modia3D.ModiaInterface

Bar = Model(
    Lx = 0.1,
    Ly = Par(value=:(0.2*Lx)),
    Lz = Par(value=:(0.2*Lx)),
    vmat1 = VisualMaterial(color="SpringGreen4", transparency=0.5),
    vmat2 = VisualMaterial(color="Red"),
    frame0 = Object3D(feature=Solid(shape=Beam(axis=1, length=:Lx, width=:Ly, thickness=:Lz),
                                    solidMaterial="Aluminium",
                                    visualMaterial=:vmat1)),
	frame1 = Object3D(parent=:frame0,
                      feature=Visual(shape=Cylinder(axis=3, diameter=:(Ly/2), length=:(1.2*Lz)),
                                     visualMaterial=:(vmat2)),
                      translation=:[-Lx/2, 0.0, 0.0]),
	frame2 = Object3D(parent=:frame0,
                      translation=:[ Lx/2, 0.0, 0.0])
)

DoublePendulum = Model(
    Lx = 1.0,
    world = Object3D(feature=Scene()),
    worldFrame = Object3D(parent=:world,
                          feature=Visual(shape=CoordinateSystem(length=:(Lx/2)))),
    bar1 = Bar | Map(Lx=:Lx),
    bar2 = Bar | Map(Lx=:Lx),
    rev1 = RevoluteWithFlange(obj1=:world, obj2=:(bar1.frame1), phi=Var(init=1.0)),
    rev2 = RevoluteWithFlange(obj1=:(bar1.frame2), obj2=:(bar2.frame1)),

    damper1 = Damper | Map(d=20.0),
    fixed1  = Fixed,
    damper2 = Damper | Map(d=20.0),
    fixed2  = Fixed,
    connect = :[(damper1.flange_b, rev1.flange),
                (damper1.flange_a, fixed1.flange),
                (damper2.flange_b, rev2.flange),
                (damper2.flange_a, fixed2.flange)]
)

doublePendulum = @instantiateModel(buildModia3D(DoublePendulum), unitless=true)

#@showModel doublePendulum.parameterExpressions
#@showModel doublePendulum.parameters
#@show doublePendulum.equationInfo

stopTime = 10.0
tolerance = 1.0e-8
requiredFinalStates = [-0.9608496178685947, 1.0258202272580021, -6.0828664801544345, -0.11676963250128454]
simulate!(doublePendulum, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(doublePendulum, [("rev1.phi", "rev2.phi"), ("rev1.w", "rev2.w")])

end


text:
__precompile__(true)
module UnionUtils

import Base.|
|(::Type{A}, ::Type{B}) where {A, B} = Union{A, B}
have(::Nothing) = false
have(x) = true
export have

end # module


text:
using SafeTestsets

@safetestset "Purebas de f." begin
    include("test_f.jl")
end


text:
@info "Starting notebook kernel server"

pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "packages"))
using VSCodeServer
popfirst!(LOAD_PATH)

@info "Core notebook support loaded"

using InteractiveUtils

let
    args = [popfirst!(Base.ARGS) for _ in 1:2]

    conn_pipeline, telemetry_pipeline = args[1:2]

    ccall(:jl_exit_on_sigint, Nothing, (Cint,), false)

    outputchannel_logger = Base.CoreLogging.SimpleLogger(Base.stderr)

    @info "Handing things off to VSCodeServer.serve_notebook"

    VSCodeServer.serve_notebook(conn_pipeline, outputchannel_logger, crashreporting_pipename=telemetry_pipeline)
end

max_stars_repo_name: mfkiwl/julia-vscode

max_issues_repo_name: mfkiwl/julia-vscode

max_forks_repo_name: mfkiwl/julia-vscode


text:
# load packages
using CSV;
using DataFrames;

# setup filepaths
path_source = string(@__DIR__,"\\..\\source");
path_dev = string(@__DIR__,"\\..\\dev");

# read in VAERS files for 2019
df_data19 = CSV.read(joinpath(path_source,"2019VAERSDATA.csv"), DataFrame);
df_symp19 = CSV.read(joinpath(path_source,"2019VAERSSYMPTOMS.csv"), DataFrame);
df_vax19 = CSV.read(joinpath(path_source,"2019VAERSVAX.csv"), DataFrame);
insertcols!(df_data19, 2, :YEAR => "2019");

# read in VAERS files for 2020
df_data20 = CSV.read(joinpath(path_source,"2020VAERSDATA.csv"), DataFrame);
df_symp20 = CSV.read(joinpath(path_source,"2020VAERSSYMPTOMS.csv"), DataFrame);
df_vax20 = CSV.read(joinpath(path_source,"2020VAERSVAX.csv"), DataFrame);
insertcols!(df_data20, 2, :YEAR => "2020");

# read in VAERS files for 2021
df_data21 = CSV.read(joinpath(path_source,"2021VAERSDATA.csv"), DataFrame);
df_symp21 = CSV.read(joinpath(path_source,"2021VAERSSYMPTOMS.csv"), DataFrame);
df_vax21 = CSV.read(joinpath(path_source,"2021VAERSVAX.csv"), DataFrame);
insertcols!(df_data21, 2, :YEAR => "2021");

# append VAERS data files from 2019-2021
print(size(df_data19));
print(size(df_data20));
print(size(df_data21));
df_data_appended = reduce(vcat, [df_data19, df_data20, df_data21]);
print(size(df_data_appended));

# write to dev
CSV.write(joinpath(path_dev,"19-21VAERSDATA.csv"), df_data_appended);

# append VAERS symptoms files from 2019-2021
print(size(df_symp19));
print(size(df_symp20));
print(size(df_symp21));
df_symp_appended = reduce(vcat, [df_symp19, df_symp20, df_symp21]);
print(size(df_symp_appended));

# write to dev
CSV.write(joinpath(path_dev,"19-21VAERSSYMPTOMS.csv"), df_symp_appended);

# append VAERS vaccine files from 2019-2021
print(size(df_vax19));
print(size(df_vax20));
print(size(df_vax21));
df_vax_appended = reduce(vcat, [df_vax19, df_vax20, df_vax21]);
print(size(df_vax_appended));

# write to dev
CSV.write(joinpath(path_dev,"19-21VAERSVAX.csv"), df_vax_appended);