text:
filename_prefix = "figures\\"

datapoint = fld(T,datat)
#virus related output
    inf_refy = zeros(datapoint+1,k1,k2)
    for t = 1:datapoint
        for i = 1:k1
            for j = 1:k2
                khelp = (i-1)*k2+j
                idex = datahhy[t,:,1] .== khelp
                inf_refy[t,i,j]= sum(datahhy[t,:,2][idex] .== 2)
            end
        end
    end

    plyinf = plot(datat * collect(0:fld(T,datat)),inf_refy[:,1,1])
    for i=1:k1
        for j=1:k2
            if i > 1 || j > 1
                plot!(datat * collect(0:fld(T,datat)),inf_refy[:,i,j])
            end
        end
    end
    savefig(plyinf, "$(filename_prefix)reg_infy.pdf")

    plycas = plot(datat * collect(0:fld(T,datat)),castraj[:,1,1,1])
    for i=1:k1
        for j=1:k2
            if i > 1 || j > 1
                plot!(datat * collect(0:fld(T,datat)),castraj[:,i,j,1])
            end
        end
    end
    savefig(plycas, "$(filename_prefix)reg_casy.pdf")

    inf_refo = zeros(datapoint+1,k1,k2)
    for t = 1:datapoint
        for i = 1:k1
            for j = 1:k2
                khelp = (i-1)*k2+j
                idex = datahho[t,:,1] .== khelp
                inf_refo[t,i,j]= sum(datahho[t,:,2][idex] .== 2)
            end
        end
    end

    ploinf = plot(datat * collect(0:fld(T,datat)),inf_refo[:,1,1])
    for i=1:k1
        for j=1:k2
            if i > 1 || j > 1
                plot!(datat * collect(0:fld(T,datat)),inf_refo[:,i,j])
            end
        end
    end
    savefig(ploinf, "$(filename_prefix)reg_info.pdf")

    plocas = plot(datat * collect(0:fld(T,datat)),castraj[:,1,1,2])
    for i=1:k1
        for j=1:k2
            if i > 1 || j > 1
                plot!(datat * collect(0:fld(T,datat)),castraj[:,i,j,2])
            end
        end
    end
    savefig(plocas, "$(filename_prefix)reg_caso.pdf")

    # plot fracction of infected
    hcapfrac = zeros(datapoint+1)
    for th = 1:datapoint+1
        hcapfrac[th] = hcap / (nhh * icufrac)
    end
    pl1 = plot(datat * collect(0:fld(T,datat)),[sum(datahhy[:,:,2] .== 2, dims=2) / nhhy, sum(datahho[:,:,2] .== 2, dims=2) / nhho, (sum(datahhy[:,:,2] .== 2, dims=2) .+ sum(datahho[:,:,2] .== 2, dims=2)) / nhh, hcapfrac], label = ["y" "o" "tot" "cap"])

    # plot mean of variable
    #pl2 = plot(mean(datahhy[:,:,2], dims=2))

    hhy = sum(castraj[:,:,:,1],dims=2)
    hho = sum(castraj[:,:,:,2],dims=2)
    pl3 = plot(datat * collect(0:fld(T,datat)),[sum(hhy[:,1,:],dims=2) / nhhy, sum(hho[:,1,:],dims=2) / nhho, (sum(hhy[:,1,:],dims=2) .+ sum(hho[:,1,:],dims=2)) / nhh], label = ["y" "o" "tot"])
    plt_all = plot(pl1,pl3, layout = 2)
    savefig(plt_all, "$(filename_prefix)inf_cas.pdf")
    plt_all

    # recovered
    pl1a = plot(datat * collect(0:fld(T,datat)),[sum(datahhy[:,:,2] .== 3, dims=2) / nhhy, sum(datahho[:,:,2] .== 3, dims=2) / nhho, (sum(datahhy[:,:,2] .== 3, dims=2) .+ sum(datahho[:,:,2] .== 3, dims=2)) / nhh], label = ["y" "o" "tot"])
    savefig(pl1a, "$(filename_prefix)recov.pdf")

    # determine RKIR0 factor of infected
    RKIR0traj = zeros(datat*datapoint+1)
    for t= 1:datat*datapoint+1
        if t > 8 && totinftraj[t-8] > 0
            global RKIR0traj[t] = (totinftraj[t]-totinftraj[t-4]) / (totinftraj[t-4]-totinftraj[t-8])
        end
    end
    pl1b = plot(collect(0:datat*datapoint), totinftraj[:])
    savefig(pl1b, "$(filename_prefix)totinf.pdf")

    pl1c = plot(collect(0:datat*datapoint), RKIR0traj[:])
    savefig(pl1c, "$(filename_prefix)RKIR0.pdf")

RKIR0smtraj = zeros(datat*datapoint+1)
for t= 1:datat*datapoint+1
    if t > 6 && RKIR0traj[t-6] > 0
        global RKIR0smtraj[t] = mean(RKIR0traj[t-6:t])
    end
end
pl1d = plot(collect(0:datat*datapoint), RKIR0smtraj[:])
savefig(pl1d, "$(filename_prefix)RKIR0sm.pdf")

pl4 = plot(datat * collect(0:fld(T,datat)),[f_av_workers[:,1],f_av_workers[:,2],f_av_workers[:,3],f_av_workers[:,4]], label = ["man" "ser" "food" "pub"])
savefig(pl4, "$(filename_prefix)av_empl.pdf")

pl4a = plot(datat * collect(0:fld(T,datat)),[f_var_workers[:,1],f_var_workers[:,2],f_var_workers[:,3],f_var_workers[:,4]], label = ["man" "ser" "food" "pub"])
savefig(pl4a, "$(filename_prefix)var_empl.pdf")

pl5 = plot(datat * collect(0:fld(T,datat)),[f_av_stock[:,1],f_av_stock[:,2],f_av_stock[:,3]], label = ["man" "ser" "food"])
savefig(pl5, "$(filename_prefix)av_stock.pdf")

consh = sum(constraj[:,:,:,:],dims=3)
pl6 = plot(datat * collect(1:fld(T,datat)),[sum(consh[2:datapoint+1,1,1,:],dims=2),sum(consh[2:datapoint+1,2,1,:],dims=2),sum(consh[2:datapoint+1,3,1,:],dims=2)], label = ["man" "ser" "food"])
savefig(pl6, "$(filename_prefix)sec_cons.pdf")

pl6a = plot(datat * collect(0:fld(T,datat)), conspercaptraj[:])
savefig(pl6a, "$(filename_prefix)cons_percap.pdf")

pl7 = plot(datat * collect(0:fld(T,datat)),[f_av_demandexp[:,1],f_av_demandexp[:,2],f_av_demandexp[:,3]], label = ["man" "ser" "food"])
savefig(pl7, "$(filename_prefix)av_demexp.pdf")

pl7a = plot(datat * collect(0:fld(T,datat)),[f_var_demandexp[:,1],f_var_demandexp[:,2],f_var_demandexp[:,3]], label = ["man" "ser" "food"])
savefig(pl7a, "$(filename_prefix)var_demexp.pdf")

pl8 = plot(datat * collect(0:fld(T,datat)),pubacctraj[:])
savefig(pl8, "$(filename_prefix)pubaccpercap.pdf")

pl9 = plot(datat * collect(0:fld(T,datat)),unempltraj[:])
savefig(pl9, "$(filename_prefix)unemp.pdf")

pl10 = plot(datat * collect(0:fld(T,datat)),tautraj[:])
savefig(pl10, "$(filename_prefix)tax.pdf")

pl11 = plot(datat * collect(0:fld(T,datat)),[f_av_savings[:,1],f_av_savings[:,2],f_av_savings[:,3]], label = ["man" "ser" "food"])
savefig(pl11, "$(filename_prefix)av_savings.pdf")

pl11a = plot(datat * collect(0:fld(T,datat)),[f_var_savings[:,1],f_var_savings[:,2],f_var_savings[:,3]], label = ["man" "ser" "food"])
savefig(pl11a, "$(filename_prefix)var_savings.pdf")

pl12 = plot(datat * collect(0:fld(T,datat)),[unempsectraj[:,1],unempsectraj[:,2],unempsectraj[:,3]], label = ["man" "ser" "food"])
savefig(pl12, "$(filename_prefix)sec_unempl.pdf")

pl13 = plot(datat * collect(0:fld(T,datat)),R0counttraj[:])
savefig(pl13, "$(filename_prefix)R0.pdf")

pl14 = plot(datat * collect(0:fld(T,datat)),gdppercaptraj[:])
savefig(pl14, "$(filename_prefix)GDPpercap.pdf")

pl14a = plot(datat * collect(0:fld(T,datat)),gdplosstraj[:])
savefig(pl14a, "$(filename_prefix)GDPloss.pdf")

pl15 = plot(datat * collect(0:fld(T,datat)),inactivefirmstraj[:])
savefig(pl15, "$(filename_prefix)inactive_firms.pdf")

pl16 = plot(datat * collect(0:fld(T,datat)),bailouttraj[:])
savefig(pl16, "$(filename_prefix)bailouts.pdf")

pl17 = plot(datat * collect(0:fld(T,datat)),shorttimetraj[:])
savefig(pl17, "$(filename_prefix)shorttime.pdf")

pl18 = plot(datat * collect(0:fld(T,datat)),[shorttimesectraj[:,1],shorttimesectraj[:,2],shorttimesectraj[:,3]], label = ["man" "ser" "food"])
savefig(pl18, "$(filename_prefix)sec_shorttime.pdf")

pl19 = plot(datat * collect(0:fld(T,datat)),totalaccountstraj[:])
savefig(pl19, "$(filename_prefix)totalaccounts.pdf")

pl20 = plot(collect(virustime:virustime+size(contact_count_traj)[1]-1), contact_count_traj[:])
savefig(pl20, "$(filename_prefix)av_contact.pdf")

pl20a = plot(collect(virustime:virustime+size(contact_work_traj)[1]-1), contact_work_traj[:])
savefig(pl20a, "$(filename_prefix)av_contact_work.pdf")

pl20b = plot(collect(virustime:virustime+size(contact_social_traj)[1]-1), contact_social_traj[:])
savefig(pl20b, "$(filename_prefix)av_contact_social.pdf")

pl20c = plot(collect(virustime:virustime+size(contact_shop_traj)[1]-1), contact_shop_traj[:])
savefig(pl20c, "$(filename_prefix)av_contact_shop.pdf")

pl21 = plot(datat * collect(0:fld(T,datat)),totfirmtraj[:])
savefig(pl21, "$(filename_prefix)totfirms.pdf")

pl22 = plot(datat * collect(0:fld(T,datat)),[bankrupttraj[:,1],bankrupttraj[:,2],bankrupttraj[:,3]], label = ["man" "ser" "food"])
savefig(pl22, "$(filename_prefix)bankruptcies.pdf")


text:
using ..Mesh.Grids: Direction, HorizontalDirection, VerticalDirection

advective_courant(m::AtmosLinearModel, a...) = advective_courant(m.atmos, a...)

nondiffusive_courant(m::AtmosLinearModel, a...) =
    nondiffusive_courant(m.atmos, a...)

diffusive_courant(m::AtmosLinearModel, a...) = diffusive_courant(m.atmos, a...)

norm_u(state::Vars, k̂::AbstractVector, ::VerticalDirection) =
    abs(dot(state.ρu, k̂)) / state.ρ
norm_u(state::Vars, k̂::AbstractVector, ::HorizontalDirection) =
    norm((state.ρu .- dot(state.ρu, k̂) * k̂) / state.ρ)
norm_u(state::Vars, k̂::AbstractVector, ::Direction) = norm(state.ρu / state.ρ)

norm_ν(ν::Real, k̂::AbstractVector, ::Direction) = ν
norm_ν(ν::AbstractVector, k̂::AbstractVector, ::VerticalDirection) = dot(ν, k̂)
norm_ν(ν::AbstractVector, k̂::AbstractVector, ::HorizontalDirection) =
    norm(ν - dot(ν, k̂) * k̂)
norm_ν(ν::AbstractVector, k̂::AbstractVector, ::Direction) = norm(ν)

function advective_courant(
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    diffusive::Vars,
    Δx,
    Δt,
    t,
    direction,
)
    k̂ = vertical_unit_vector(m, aux)
    normu = norm_u(state, k̂, direction)
    return Δt * normu / Δx
end

function nondiffusive_courant(
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    diffusive::Vars,
    Δx,
    Δt,
    t,
    direction,
)
    k̂ = vertical_unit_vector(m, aux)
    normu = norm_u(state, k̂, direction)
    return Δt * (normu + soundspeed(m, m.moisture, state, aux)) / Δx
end

function diffusive_courant(
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    diffusive::Vars,
    Δx,
    Δt,
    t,
    direction,
)
    ν, _, _ = turbulence_tensors(m, state, diffusive, aux, t)
    ν = ν isa Real ? ν : diag(ν)
    k̂ = vertical_unit_vector(m, aux)
    normν = norm_ν(ν, k̂, direction)
    return Δt * normν / (Δx * Δx)
end

ext: jl


text:
module TestDiary

using Test, Diary

using FileWatching
using Pkg

@testset "History Parsing" begin
    history_lines = [
        "# time: ***",
        "# mode: julia",
	"\ta = rand(100);",
    ]
    @test Diary.parse_history(history_lines) == ["a = rand(100)"]

    history_lines = [
        "# time: ***",
        "# mode: julia",
	"\tstruct MyStruct",
	"\t    field::Any",
	"\tend",
    ]
    @test Diary.parse_history(history_lines) == ["struct MyStruct", "    field::Any", "end"]

    history_lines = [
        "# time: ***",
        "# mode: pkg",
	"\tstatus",
    ]
    @test Diary.parse_history(history_lines) == []

    history_lines = [
        "# time: ***",
        "# mode: julia",
	"\t# diary: commit",
    ]
    @test Diary.parse_history(history_lines) == []
end

@testset "Command Parsing" begin
    cmd = "commit"
    @test Diary.parse_command(cmd) === 0
    cmd = "commit 1"
    @test Diary.parse_command(cmd) === 0
    cmd = "commit 1 2"
    @test_logs (:error, "Diary.jl: too many arguments to `commit`: 2, expected 1") Diary.parse_command(cmd)
    cmd = "erase"
    @test Diary.parse_command(cmd) === true
    cmd = "erase diary"
    @test_logs (:error, "Diary.jl: too many arguments to `erase`: 1, expected 0") Diary.parse_command(cmd)
    cmd = "unsupported command"
    @test_logs (:error, "Diary.jl: could not parse command: $cmd") Diary.parse_command(cmd)
end

@testset "Finding configuration" begin
    haskey(ENV, "JULIA_DIARY_CONFIG") && return

    Pkg.activate()
    default_configuration = Dict{String,Any}(
        "author" => "",
        "autocommit" => true,
        "blacklist" => [
            joinpath(DEPOT_PATH[1], "environments"),
        ],
        "create_if_missing" => true,
        "date_format" => "E U d HH:MM",
        "diary_name" => "diary.jl",
        "directory_mode" => false,
        "file_polling" => false,
        "file_polling_interval" => 1.0,
        "persistent_history" => true,
    )

    ENV["JULIA_DIARY_CONFIG"] = tempname()
    @test Diary.default_configuration() == default_configuration
    @test Diary.read_configuration()["author"] == default_configuration["author"]
    delete!(ENV, "JULIA_DIARY_CONFIG")

    configuration_file = joinpath(
        dirname(Pkg.project().path),
        "Diary.toml",
    )
    touch(configuration_file)
    @test Diary.find_configuration_file() == configuration_file

    open(configuration_file, write=true) do io
        println(io, "author = \"Test\"")
    end
    @test Diary.read_configuration()["author"] == "Test"

    open(configuration_file, write=true) do io
        println(io, "12 34")
    end

    @test_logs(
        (
            :warn,
            string(
                "Diary.jl: an error occured while parsing configuration file: ",
                configuration_file,
            )
        ),
        Diary.read_configuration(),
    )

    # Clean-up
    rm(configuration_file)
end

@testset "Finding Diary" begin
    @testset "Basic functionality" begin
        Pkg.activate()
        default_path = joinpath(
            dirname(Pkg.project().path),
            "diary.jl",
        )
        @test Diary.find_diary() == default_path
        @test isfile(default_path)
    end

    @testset "Environment variable" begin
        env_path = joinpath(dirname(Pkg.project().path), "env_diary.jl")
        ENV["JULIA_DIARY"] = env_path
        @test Diary.find_diary() == env_path
        @test isfile(env_path)
        delete!(ENV, "JULIA_DIARY")
    end

    @testset "Blacklist" begin
        # Test that the default environment is blacklisted by default.
        Pkg.activate("v$(VERSION.major).$(VERSION.minor)"; shared=true)
        @test isnothing(Diary.find_diary(; configuration=Diary.default_configuration()))
    end

    @testset "Create missing" begin
        Pkg.activate()

        configuration = merge!(
            Diary.default_configuration(),
            Dict{String,Any}(
                "create_if_missing" => false,
            ),
        )

        diary_file_path = Diary.find_diary_path(; configuration=configuration)
        rm(diary_file_path, force=true)
        @test isnothing(Diary.find_diary(; configuration=configuration))
    end

    @testset "Directory mode" begin
        Pkg.activate()

        configuration = merge!(
            Diary.default_configuration(),
            Dict{String,Any}(
                "directory_mode" => true,
            ),
        )

        directory = mktempdir(dirname(Pkg.project().path))
        expected_dir = abspath(joinpath(directory, "diary.jl"))
        cd(directory) do
            @test Diary.find_diary(; configuration=configuration) == expected_dir
        end
    end
end

@testset "Header" begin
    Pkg.activate()

    configuration = merge!(
        Diary.default_configuration(),
        Dict{String,Any}(
            "author" => "Test",
            "date_format" => "",
        ),
    )

    diary_file = Diary.find_diary()
    open(diary_file, write=true) do io
        Diary.write_header(io; configuration=configuration)
    end

    @test readline(diary_file) == "# Test: "

    # Clean-up
    rm(diary_file)
end

@testset "Committing" begin
    Pkg.activate()

    configuration = merge!(
        Diary.default_configuration(),
        Dict{String,Any}(
            "author" => "Test",
            "date_format" => "",
        ),
    )

    push!(Diary.GLOBAL_SEGMENT_BUFFER, ["a = rand(100)"])

    diary_file = Diary.find_diary(; configuration=configuration)
    Diary.commit(; diary_file=diary_file, configuration=configuration)

    @test readlines(Diary.find_diary(; configuration=configuration)) == ["# Test: ", "", "a = rand(100)"]

    push!(Diary.GLOBAL_SEGMENT_BUFFER, ["function f(x)", "    return x^2", "end"])

    diary_file = Diary.find_diary(; configuration=configuration)
    Diary.commit(; diary_file=diary_file, configuration=configuration, with_header=false)

    expected_output = [
        "# Test: ",
        "",
        "a = rand(100)",
        "",
        "function f(x)",
        "    return x^2",
        "end",
        "",
    ]
    @test readlines(Diary.find_diary(; configuration=configuration)) == expected_output

    @testset "Erasing" begin
        Diary.erase_diary(; diary_file=diary_file, configuration=configuration)
        @test readlines(Diary.find_diary(; configuration=configuration)) == []
    end

    # Clean-up
    rm(diary_file)
end

@testset "Interactivity" begin
    Pkg.activate()
    # Re-initialise Diary, enabling the watcher.
    ENV["JULIA_HISTORY"] = tempname()
    touch(ENV["JULIA_HISTORY"])
    repl_history_file = ENV["JULIA_HISTORY"]
    Diary.__init__(enabled=true)
    # Locate relevant files.
    history_file = ENV["JULIA_HISTORY"]
    configuration_file = joinpath(dirname(Pkg.project().path), "Diary.toml")
    # Write the configuration.
    open(configuration_file, write=true) do io
        join(io, ["author = \"Test\"", "date_format = \"\""], "\n")
        print(io, "\n")
    end
    diary_file = Diary.find_diary()

    @test Diary.read_configuration()["author"] == "Test"
    @test Diary.read_configuration()["date_format"] == ""

    # Simulate user interaction by writing lines to the history file.
    history_lines = [
        "# time: ***",
        "# mode: julia",
	"\ta = rand(100);",
    ]
    open(history_file, read=true, write=true) do io
        seekend(io)
        join(io, history_lines, "\n")
        print(io, "\n")
    end
    # Allow the `watch_task` to update the diary file.
    file_event = FileWatching.watch_file(diary_file, 5)
    @test file_event.changed || file_event.timedout
    @test readlines(diary_file) == ["# Test: ", "", "a = rand(100)"]
    @test readlines(repl_history_file) == history_lines

    # Add line with incorrect syntax and check that it does not update the diary file.
    history_lines = [
        "# time: ***",
        "# mode: julia",
	"\ta b c d e",
    ]
    open(history_file, read=true, write=true) do io
        seekend(io)
        join(io, history_lines, "\n")
        print(io, "\n")
    end
    # Allow the `watch_task` to update the diary file.
    file_event = FileWatching.watch_file(diary_file, 5)
    @test file_event.timedout
    @test readlines(diary_file) == ["# Test: ", "", "a = rand(100)"]

    @testset "Non-persistent history" begin
        previous_history_lines = readlines(repl_history_file)

        open(configuration_file, write=true) do io
            configuration = [
                "author = \"Test\"",
                "date_format = \"\"",
                "persistent_history = false",
            ]
            join(io, configuration, "\n")
            print(io, "\n")
        end
        # Simulate user interaction by writing lines to the history file.
        open(history_file, read=true, write=true) do io
            seekend(io)
            join(io, ["# time: ***", "# mode: julia", "\tb = 42"], "\n")
            print(io, "\n")
        end
        # Allow the `watch_task` to update the diary file.
        file_event = FileWatching.watch_file(diary_file, 5)
        @test file_event.changed || file_event.timedout
        @test readlines(diary_file) == ["# Test: ", "", "a = rand(100)", "b = 42"]
        # Test that the original history file hasn't changed.
        @test readlines(repl_history_file) == previous_history_lines
    end

    @testset "File polling" begin
        open(configuration_file, write=true) do io
            configuration = [
                "author = \"Test\"",
                "date_format = \"\"",
                "file_polling = true",
            ]
            join(io, configuration, "\n")
            print(io, "\n")
        end
        # Trigger a configuration refresh to make the task switch to polling, and erase the
        # diary.
        open(history_file, read=true, write=true) do io
            seekend(io)
            join(io, ["# time: ***", "# mode: julia", "\t# diary: erase"], "\n")
            print(io, "\n")
        end
        # Allow the diary to synchronise before wiping it.
        sleep(5)
        # Write the actual history lines we want to check up against.
        open(history_file, read=true, write=true) do io
            seekend(io)
            join(io, ["# time: ***", "# mode: julia", "\tc = 13"], "\n")
            print(io, "\n")
        end
        file_event = FileWatching.watch_file(diary_file, 5)
        @test file_event.changed || file_event.timedout
        @test readlines(diary_file) == ["c = 13"]
    end

    # Clean-up
    rm(configuration_file)
    rm(diary_file)
end

end # module


text:
function gsr(
    equation::Union{String, Array{String}, Array{Symbol}},
    data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}, Tuple, DataFrame, Nothing};
    datanames::Union{Array, Array{Symbol, 1}, Nothing}=nothing,
    method::Union{Symbol, String}=Preprocessing.METHOD_DEFAULT,
    intercept::Bool=Preprocessing.INTERCEPT_DEFAULT,
    panel::Union{Symbol, String, Nothing}=Preprocessing.PANEL_DEFAULT,
    time::Union{Symbol, String, Nothing}=Preprocessing.TIME_DEFAULT,
    seasonaladjustment::Union{Dict, Array, Nothing}=Preprocessing.SEASONALADJUSTMENT_DEFAULT,
    removeoutliers::Bool=Preprocessing.REMOVEOUTLIERS_DEFAULT,
    fe_sqr::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_log::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_inv::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_lag::Union{Nothing, Array}=nothing,
    interaction::Union{Nothing, Array}=nothing,
    preliminaryselection::Union{Nothing, Symbol}=nothing,
    fixedvariables::Union{Nothing, Array}=AllSubsetRegression.FIXEDVARIABLES_DEFAULT,
    outsample::Union{Nothing, Int, Array}=AllSubsetRegression.OUTSAMPLE_DEFAULT,
    criteria::Array=AllSubsetRegression.CRITERIA_DEFAULT,
    ttest::Bool=AllSubsetRegression.TTEST_DEFAULT,
    modelavg::Bool=AllSubsetRegression.MODELAVG_DEFAULT,
    residualtest::Bool=AllSubsetRegression.RESIDUALTEST_DEFAULT,
    orderresults::Bool=AllSubsetRegression.ORDERRESULTS_DEFAULT,
    kfoldcrossvalidation::Bool=KFOLDCROSSVALIDATION_DEFAULT,
    numfolds::Int=NUMFOLDS_DEFAULT,
    testsetshare::Union{Float32, Float64}=TESTSETSHARE_DEFAULT,
    exportcsv::Union{Nothing, String}=EXPORTCSV_DEFAULT,
    exportsummary::Union{Nothing, String}=EXPORTSUMMARY_DEFAULT,
    exportlatex::Union{Nothing, String}=EXPORTLATEX_DEFAULT
    )

    gsr(
        equation,
        data=data,
        datanames=datanames,
        method=method,
        intercept=intercept,
        panel=panel,
        time=time,
        seasonaladjustment=seasonaladjustment,
        removeoutliers=removeoutliers,
        fe_sqr=fe_sqr,
        fe_log=fe_log,
        fe_inv=fe_inv,
        fe_lag=fe_lag,
        interaction=interaction,
        preliminaryselection=preliminaryselection,
        fixedvariables=fixedvariables,
        outsample=outsample,
        criteria=criteria,
        ttest=ttest,
        modelavg=modelavg,
        residualtest=residualtest,
        orderresults=orderresults,
        kfoldcrossvalidation=kfoldcrossvalidation,
        numfolds=numfolds,
        testsetshare=testsetshare,
        exportcsv=exportcsv,
        exportsummary=exportsummary,
        exportlatex=exportlatex
    )
end

function gsr(
    equation::Union{String, Array{String}, Array{Symbol}};
    data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}, Tuple, DataFrame, Nothing},
    datanames::Union{Array, Array{Symbol, 1}, Nothing}=nothing,
    method::Union{Symbol, String}=Preprocessing.METHOD_DEFAULT,
    intercept::Bool=Preprocessing.INTERCEPT_DEFAULT,
    panel::Union{Symbol, String, Nothing}=Preprocessing.PANEL_DEFAULT,
    time::Union{Symbol, String, Nothing}=Preprocessing.TIME_DEFAULT,
    seasonaladjustment::Union{Dict, Array, Nothing}=Preprocessing.SEASONALADJUSTMENT_DEFAULT,
    removeoutliers::Bool=Preprocessing.REMOVEOUTLIERS_DEFAULT,
    fe_sqr::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_log::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_inv::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_lag::Union{Nothing, Array}=nothing,
    interaction::Union{Nothing, Array}=nothing,
    preliminaryselection::Union{Nothing, Symbol}=nothing,
    fixedvariables::Union{Nothing, Array}=AllSubsetRegression.FIXEDVARIABLES_DEFAULT,
    outsample::Union{Nothing, Int, Array}=AllSubsetRegression.OUTSAMPLE_DEFAULT,
    criteria::Array=AllSubsetRegression.CRITERIA_DEFAULT,
    ttest::Bool=AllSubsetRegression.TTEST_DEFAULT,
    modelavg::Bool=AllSubsetRegression.MODELAVG_DEFAULT,
    residualtest::Bool=AllSubsetRegression.RESIDUALTEST_DEFAULT,
    orderresults::Bool=AllSubsetRegression.ORDERRESULTS_DEFAULT,
    kfoldcrossvalidation::Bool=KFOLDCROSSVALIDATION_DEFAULT,
    numfolds::Int=NUMFOLDS_DEFAULT,
    testsetshare::Union{Float32, Float64}=TESTSETSHARE_DEFAULT,
    exportcsv::Union{Nothing, String}=EXPORTCSV_DEFAULT,
    exportsummary::Union{Nothing, String}=EXPORTSUMMARY_DEFAULT,
    exportlatex::Union{Nothing, String}=EXPORTLATEX_DEFAULT
)

    if fe_lag == nothing
        removemissings = true
    else
        removemissings = false
    end

    data = Preprocessing.input(
        equation,
        data=data,
        datanames=datanames,
        method=method,
        intercept=intercept,
        panel=panel,
        time=time,
        seasonaladjustment=seasonaladjustment,
        removeoutliers=removeoutliers,
        removemissings=removemissings
    )

    if featureextraction_enabled(fe_sqr, fe_log, fe_inv, fe_lag, interaction)
        data = FeatureExtraction.featureextraction!(data, fe_sqr=fe_sqr, fe_lag=fe_lag, fe_log=fe_log, fe_inv=fe_inv, interaction=interaction, removemissings=true)
    end

    original_data = copy_data(data)

    if preliminaryselection_enabled(preliminaryselection)
        preliminaryselection = Symbol(preliminaryselection)
        if !validate_preliminaryselection(preliminaryselection)
            error(INVALID_PRELIMINARYSELECTION)
        end
        data = PreliminarySelection._lasso!(data)
        original_data.extras = data.extras
    end

    AllSubsetRegression.ols!(
        data,
        fixedvariables=fixedvariables,
        outsample=outsample,
        criteria=criteria,
        ttest=ttest,
        modelavg=modelavg,
        residualtest=residualtest,
        orderresults=orderresults
    )

    original_data.extras = data.extras

    if kfoldcrossvalidation
        CrossValidation.kfoldcrossvalidation!(data, original_data, numfolds, testsetshare)
    end

    if exportcsv != nothing
        ModelSelection.Output.csv(data, filename=exportcsv)
    end

    if exportlatex != nothing
        ModelSelection.Output.latex(data, original_data, path=exportlatex)
    end

    println(ModelSelection.Output.summary(data, filename=exportsummary))

    return data
end


text:
using MixedIntegerExperiments
using Polyhedra
using Base.Test

import MixedIntegerExperiments: perp

@testset "util" begin
    p1 = SimpleHRepresentation(rand(2, 2), rand(2))
    p2 = SimpleHRepresentation(rand(2, 2), rand(2))
    pcart = p1 * p2
    @test fulldim(pcart) == fulldim(p1) + fulldim(p2)
    vr1, vr2, vrcart = map(p -> vrep(polyhedron(p)), (p1, p2, pcart))
    for v1 in points(vr1), v2 in points(vr2)
        @test vcat(v1, v2) ∈ points(vrcart)
    end
end

@testset "region descriptions" begin
    # FIXME:
    # for i = 1 :  10
    #     p1 = randn(2)
    #     p2 = randn(2)
    #     μ = rand()
    #     fmax = 20.
    #     region = contact_region(p1, p2, μ, fmax)
    #     n = normalize(perp(p2 - p1))
    #     @test all(region.force.A * fmax * (1 - eps()) * n .>= region.force.b)
    #
    #     for j = 1 : 100
    #         f = randn(2)
    #         infrictioncone = norm(f) <= sqrt(1 + μ^2) * dot(f, n)
    #         @test all(region.force.A * f .>= region.force.b) == infrictioncone
    #     end
    # end
    #
    # for i = 1 : 10
    #     p1 = rand(2)
    #     p2 = rand(2)
    #     xaxis = [1.; 0]
    #     region = axis_aligned_free_box_region(p1, p2, xaxis)
    #
    #     for j = 1 : 100
    #         λ1 = rand()
    #         λ = [λ1; 1 - λ1]
    #         p = [p1 p2] * λ
    #         @test all(region.position.A * p .>= region.position.b)
    #     end
    #     for j = 1 : 100
    #         p = randn(2)
    #         λ1 = randn()
    #         p = λ1 * (p1 - p2) + p2
    #         @test (zero(λ1) <= λ1 <= one(λ1)) == all(region.position.A * p .>= region.position.b)
    #     end
    # end
end


text:
using DSGE, HDF5, JLD
using Base.Test

path = dirname(@__FILE__)

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:saveroot, tempdir())
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:date_conditional_end, quartertodate("2015-Q4"))
m <= Setting(:forecast_uncertainty_override, Nullable(false))
m <= Setting(:use_population_forecast, true)
m <= Setting(:compute_shockdec_bands, true)

estroot = normpath(joinpath(dirname(@__FILE__), "..", "reference"))
overrides = forecast_input_file_overrides(m)
overrides[:mode] = joinpath(estroot, "optimize.h5")
overrides[:full] = joinpath(estroot, "metropolis_hastings.h5")

output_vars = add_requisite_output_vars([:histpseudo, :histobs,
                                         :hist4qpseudo, :hist4qobs,
                                         :forecastpseudo, :forecastobs,
                                         :forecast4qpseudo, :forecast4qobs,
                                         :shockdecpseudo, :shockdecobs,
                                         :irfpseudo, :irfobs])

# Read expected output
exp_modal_means, exp_modal_bands, exp_full_means, exp_full_bands =
    jldopen("$path/../reference/means_bands_out.jld", "r") do file
        read(file, "exp_modal_means"), read(file, "exp_modal_bands"),
        read(file, "exp_full_means"),  read(file, "exp_full_bands")
    end

# Modal
forecast_one(m, :mode, :none, output_vars, verbose = :none)
compute_meansbands(m, :mode, :none, output_vars; compute_shockdec_bands = true, verbose = :none)
meansbands_to_matrix(m, :mode, :none, output_vars; verbose = :none)

@testset "Check modal meansbands computation" begin
    for var in output_vars
        filename = get_forecast_filename(m, :mode, :none, Symbol("mb_matrix_", var),
                                         pathfcn = workpath, fileformat = :h5)
        @test @test_matrix_approx_eq exp_modal_means[var] h5read(filename, "means")
        @test @test_matrix_approx_eq exp_modal_bands[var] h5read(filename, "bands")
    end
end

# Full-distribution
@everywhere using DSGE
m <= Setting(:forecast_block_size, 5)
forecast_one(m, :full, :none, output_vars, verbose = :none)
compute_meansbands(m, :full, :none, output_vars; compute_shockdec_bands = true, verbose = :none)
meansbands_to_matrix(m, :full, :none, output_vars; verbose = :none)

@testset "Check full meansbands computation" begin
    for var in output_vars
        filename = get_forecast_filename(m, :full, :none, Symbol("mb_matrix_", var),
                                         pathfcn = workpath, fileformat = :h5)
        @test @test_matrix_approx_eq exp_full_means[var] h5read(filename, "means")
        @test @test_matrix_approx_eq exp_full_bands[var] h5read(filename, "bands")
    end
end

nothing


text:
using MPI, PETSc2

  PetscInitialize()
  b = PetscVec(10, VECMPI, MPI.COMM_WORLD)

  idx = PetscInt[1, 2, 3]
  vals = PetscScalar[1.0, 2.0, 3.0]

  set_values1!(b, idx, vals, ADD_VALUES)

  # don't assemble the vector
  PetscView(b)

  PetscDestroy(b)


text:
mutable struct TimeDerivatives{N,MD,ED}
    accelerations::Vector{SVectorF{N}}
    mass_derivatives::MD
    energy_derivatives::ED

    function TimeDerivatives(
        velocities::Velocities{N},
        mass_component::MassComponent{MD,<:Kernel{N}},
        energy_component::EnergyComponent{ED},
    ) where {N,MD,ED}
        new{N,MD,ED}(
            initderivatives(velocities),
            initderivatives(mass_component),
            initderivatives(energy_component),
        )
    end
end

"""
    initderivatives!(
        derivatives::TimeDerivatives{N,MD,ED},
        velocities::Velocities{N},
        mass_component::MassComponent{MD,<:Kernel{N}},
        energy_component::EnergyComponent{ED},
    )

Resets the given fluid derivatives.
"""
function initderivatives!(
    derivatives::TimeDerivatives{N,MD,ED},
    velocities::Velocities{N},
    mass_component::MassComponent{MD,<:Kernel{N}},
    energy_component::EnergyComponent{ED},
) where {N,MD,ED}
    initderivatives!(derivatives.accelerations, velocities)
    initderivatives!(derivatives.mass_derivatives, mass_component)
    initderivatives!(derivatives.energy_derivatives, energy_component)
end

"""
    updatefluid!(
        time_derivatives::TimeDerivatives{N,MD,ED},
        mass_component::MassComponent{MD,<:Kernel{N}},
        energy_component::EnergyComponent{ED},
        diffusion_component::DiffusionComponent,
        pressures::Pressures,
        velocities::Velocities{N},
        positions::Vector{SVectorF{N}},
        eos::EquationOfState,
        boundary_component::BoundaryComponent,
    )

Update the state of the fluid variables and derivatives according to the
current particle positions.
"""
function updatefluid!(
    time_derivatives::TimeDerivatives{N,MD,ED},
    mass_component::MassComponent{MD,<:Kernel{N}},
    energy_component::EnergyComponent{ED},
    diffusion_component::DiffusionComponent,
    pressures::Pressures,
    velocities::Velocities{N},
    positions::Vector{SVectorF{N}},
    eos::EquationOfState,
    boundary_component::BoundaryComponent{N},
) where {N,MD,ED}
    updatemasses!(mass_component, positions)
    updatepressures!(pressures, mass_component, energy_component, eos)
    updatediffusion!(diffusion_component, energy_component, eos)
    updateboundaries!(
        boundary_component,
        positions,
        velocities,
        mass_component,
        pressures,
    )

    initderivatives!(
        time_derivatives,
        velocities,
        mass_component,
        energy_component,
    )

    n::UInt = length(positions)
    for i = 1:n-1, j = i:n
        rᵢⱼ = positions[j] - positions[i]
        r = norm(rᵢⱼ)
        hᵢ_hⱼ = get_kernel_widths(mass_component, (i, j))
        qᵢ_qⱼ = normalize_distance(mass_component, r, hᵢ_hⱼ)

        if r > 0 && nonzeroat(get_kernel(mass_component), qᵢ_qⱼ)
            r̂ᵢⱼ = rᵢⱼ / r
            updatederivatives!(
                time_derivatives,
                (i, j),
                (rᵢⱼ, r̂ᵢⱼ, r),
                hᵢ_hⱼ,
                qᵢ_qⱼ,
                velocities,
                pressures,
                mass_component,
                energy_component,
                diffusion_component,
            )
        end
    end

    boundary_positions = get_positions(boundary_component)
    nb::UInt = length(boundary_positions)
    for i = 1:n, j = 1:nb
        rᵢⱼ = boundary_positions[j] - positions[i]
        r = norm(rᵢⱼ)
        hᵢ = get_kernel_width(mass_component, i)
        qᵢ = normalize_distance(mass_component, r, hᵢ)

        if r > 0 && nonzeroat(get_kernel(mass_component), qᵢ)
            r̂ᵢⱼ = rᵢⱼ / r
            updatederivatives_boundary!(
                time_derivatives,
                (i, j),
                (rᵢⱼ, r̂ᵢⱼ, r),
                hᵢ,
                qᵢ,
                velocities,
                pressures,
                mass_component,
                energy_component,
                diffusion_component,
                boundary_component,
            )
        end
    end
end

"""
    updatederivatives!(
        time_derivatives::TimeDerivatives{N,MD,ED},
        i_j::Tuple{Unsigned,Unsigned},
        (rᵢⱼ, r̂ᵢⱼ, r)::Tuple{SVectorF{N},SVectorF{N},Number},
        hᵢ_hⱼ::Union{Number,Tuple{Number,Number}},
        qᵢ_qⱼ::Union{Number,Tuple{Number,Number}},
        velocities::Velocities{N},
        pressures::Pressures,
        mass_component::MassComponent{MD,<:Kernel{N}},
        energy_component::EnergyComponent{ED},
        diffusion_component::DiffusionComponent,
    )

Update time derivatives for particles `i` and `j`.
"""
function updatederivatives!(
    time_derivatives::TimeDerivatives{N,MD,ED},
    i_j::Tuple{Unsigned,Unsigned},
    (rᵢⱼ, r̂ᵢⱼ, r)::Tuple{SVectorF{N},SVectorF{N},Number},
    hᵢ_hⱼ::Union{Number,Tuple{Number,Number}},
    qᵢ_qⱼ::Union{Number,Tuple{Number,Number}},
    velocities::Velocities{N},
    pressures::Pressures,
    mass_component::MassComponent{MD,<:Kernel{N}},
    energy_component::EnergyComponent{ED},
    diffusion_component::DiffusionComponent,
) where {N,MD,ED}
    vᵢ, vⱼ = get_velocities(velocities, i_j)
    vᵢⱼ = vⱼ - vᵢ
    vᵣ = vᵢⱼ ⋅ r̂ᵢⱼ

    ρᵢ_ρⱼ = get_mass_densities(mass_component, i_j)
    Pᵢ_Pⱼ = get_pressures(pressures, i_j)

    ddr_W_r_hᵢ_hⱼ = ddr(get_kernel(mass_component), qᵢ_qⱼ, hᵢ_hⱼ)

    m_ddr_W_r_hᵢ_hⱼ = compute_m_ddr_W_r_h(i_j, ddr_W_r_hᵢ_hⱼ, mass_component)

    Δaᵢ_pressure, Δaⱼ_pressure =
        compute_Δa_pressure(ρᵢ_ρⱼ, Pᵢ_Pⱼ, m_ddr_W_r_hᵢ_hⱼ)

    Δa_diffusion = compute_Δa_diffusion(
        r,
        vᵣ,
        hᵢ_hⱼ,
        ρᵢ_ρⱼ,
        ddr_W_r_hᵢ_hⱼ,
        mass_component,
        diffusion_component,
    )

    updatederivatives!(
        time_derivatives.mass_derivatives,
        i_j,
        vᵣ,
        m_ddr_W_r_hᵢ_hⱼ,
        mass_component,
    )
    updatederivatives!(
        time_derivatives.energy_derivatives,
        i_j,
        vᵣ,
        (Δaᵢ_pressure, Δaⱼ_pressure, Δa_diffusion),
        energy_component,
    )
    updatederivatives!(
        time_derivatives.accelerations,
        i_j,
        r̂ᵢⱼ,
        (Δaᵢ_pressure, Δaⱼ_pressure, Δa_diffusion),
        velocities,
    )
end

"""
    updatederivatives_boundary!(
        time_derivatives::TimeDerivatives{N,MD,ED},
        (i, j)::Tuple{Unsigned,Unsigned},
        (rᵢⱼ, r̂ᵢⱼ, r)::Tuple{SVectorF{N},SVectorF{N},Number},
        hᵢ::Number,
        qᵢ::Number,
        velocities::Velocities{N},
        pressures::Pressures,
        mass_component::MassComponent{MD,<:Kernel{N}},
        energy_component::EnergyComponent{ED},
        diffusion_component::DiffusionComponent,
        boundary_component::BoundaryComponent{N},
    )

Update time derivatives for particle `i` due to boundary particle `j`.
"""
function updatederivatives_boundary!(
    time_derivatives::TimeDerivatives{N,MD,ED},
    (i, j)::Tuple{Unsigned,Unsigned},
    (rᵢⱼ, r̂ᵢⱼ, r)::Tuple{SVectorF{N},SVectorF{N},Number},
    hᵢ::Number,
    qᵢ::Number,
    velocities::Velocities{N},
    pressures::Pressures,
    mass_component::MassComponent{MD,<:Kernel{N}},
    energy_component::EnergyComponent{ED},
    diffusion_component::DiffusionComponent,
    boundary_component::BoundaryComponent{N},
) where {N,MD,ED}
    vᵢ = get_velocity(velocities, i)
    vⱼ = get_velocity(boundary_component, j)
    vᵢⱼ = vⱼ - vᵢ
    vᵣ = vᵢⱼ ⋅ r̂ᵢⱼ

    ρᵢ = get_mass_density(mass_component, i)
    ρⱼ = compute_mass_density(boundary_component, energy_component, eos, i, j)

    Pᵢ = get_pressure(pressures, i)
    Pⱼ = get_pressure(boundary_component, j)

    ddr_W_r_hᵢ = ddr(get_kernel(mass_component), qᵢ, hᵢ)

    m_ddr_W_r_hᵢ = compute_m_ddr_W_r_h(i, ddr_W_r_hᵢ, mass_component)

    Δaᵢ_pressure, Δaⱼ_pressure =
        compute_Δa_pressure((ρᵢ, ρⱼ), (Pᵢ, Pⱼ), m_ddr_W_r_hᵢ)

    Δa_diffusion = compute_Δa_diffusion(
        r,
        vᵣ,
        hᵢ,
        (ρᵢ, ρⱼ),
        ddr_W_r_hᵢ,
        mass_component,
        diffusion_component,
    )

    updatederivative!(
        time_derivatives.mass_derivatives,
        i,
        vᵣ,
        ddr_W_r_hᵢ,
        mass_component,
    )
    updatederivative!(
        time_derivatives.energy_derivatives,
        i,
        vᵣ,
        (Δaᵢ_pressure, Δaⱼ_pressure, Δa_diffusion),
        energy_component,
    )
    updatederivative!(
        time_derivatives.accelerations,
        i,
        r̂ᵢⱼ,
        (Δaᵢ_pressure, Δaⱼ_pressure, Δa_diffusion),
        velocities,
    )
end

"""
    evolvevariables!(
        velocities::Velocities{N},
        mass_component::MassComponent{MD},
        energy_component::EnergyComponent{ED},
        time_derivatives::TimeDerivatives{N,MD,ED},
        Δt::Number,
    )

Evolve the components with the given derivatives and time step.
"""
function evolvevariables!(
    velocities::Velocities{N},
    mass_component::MassComponent{MD},
    energy_component::EnergyComponent{ED},
    time_derivatives::TimeDerivatives{N,MD,ED},
    Δt::Number,
) where {N,MD,ED}
    evolvevariables!(velocities, time_derivatives.accelerations, Δt)
    evolvevariables!(mass_component, time_derivatives.mass_derivatives, Δt)
    evolvevariables!(energy_component, time_derivatives.energy_derivatives, Δt)
end

"""
    evolvevariables!(
        velocities::Velocities{N},
        mass_component::MassComponent{MD},
        energy_component::EnergyComponent{ED},
        Δt::Number,
        time_derivatives_and_weights::Tuple{TimeDerivatives{N,MD,ED},Float}...,
    )

Evolve the components with the given time step and weighed sets of derivatives.
"""
function evolvevariables!(
    velocities::Velocities{N},
    mass_component::MassComponent{MD},
    energy_component::EnergyComponent{ED},
    Δt::Number,
    time_derivatives_and_weights::Tuple{TimeDerivatives{N,MD,ED},Float}...,
) where {N,MD,ED}
    for (time_derivatives, weight) in time_derivatives_and_weights
        weighted_Δt = weight * Δt
        evolvevariables!(
            velocities,
            mass_component,
            energy_component,
            time_derivatives,
            weighted_Δt,
        )
    end
end

"""
    evolvepositions!(
        positions::Vector{SVectorF{N}},
        Δt::Number,
        velocities_and_weights::Tuple{Velocities{N},Number}...,
    )

Evolves the positions with the given time step weighed sets of velocities.
"""
function evolvepositions!(
    positions::Vector{SVectorF{N}},
    Δt::Number,
    velocities_and_weights::Tuple{Velocities{N},Number}...,
) where {N}
    for (velocities, weight) in velocities_and_weights
        weighted_Δt = weight * Δt
        evolvepositions!(positions, velocities, weighted_Δt)
    end
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/EAGO.jl
# The main file for EAGO.
#############################################################################

__precompile__()

module EAGO

    import MathOptInterface

    using Reexport, Requires, Cassette, IntervalArithmetic, DocStringExtensions,
          FastRounding, SpecialFunctions, Ipopt, Cbc, Printf, PrettyTables

    using JuMP
    import JuMP
    import JuMP._Derivatives: operators, NodeData
    using JuMP._Derivatives: univariate_operators,
                             univariate_operator_to_id
    import JuMP: _SubexpressionStorage
    import JuMP._Derivatives: NodeType, UserOperatorRegistry
    const JuMPOpReg = JuMP._Derivatives.UserOperatorRegistry

    using DataStructures: OrderedDict, BinaryMinMaxHeap, popmin!, popmax!, top
    using SparseArrays: SparseMatrixCSC, spzeros, rowvals, nzrange, nonzeros, sparse, findnz
    using LinearAlgebra: eigmin, norm
    using Base: @propagate_inbounds
    import Base: setindex!, + , *, -, ^, /, zero, one, inv, log, log10, exp, exp10, isempty

    @reexport using McCormick
    @reexport using SpecialFunctions
    #@reexport using ReverseMcCormick

    #using IntervalContractors
    using IntervalContractors
    #=
    using IntervalContractors: plus_rev, minus_rev, inv_rev,
       mul_rev, div_rev, power_rev,
       max_rev, min_rev,
        sqr_rev, sqrt_rev, abs_rev,
        exp_rev, exp2_rev, exp10_rev, expm1_rev,
        log_rev, log2_rev, log10_rev, log1p_rev,
        sin_rev, cos_rev, tan_rev,
        asin_rev, acos_rev, atan_rev,
        sinh_rev, cosh_rev, tanh_rev,
        asinh_rev, acosh_rev, atanh_rev
        =#

    const MOI = MathOptInterface
    const MOIU = MOI.Utilities
    const MOIB = MOI.Bridges

    const SAF = MOI.ScalarAffineFunction{Float64}
    const SQF = MOI.ScalarQuadraticFunction{Float64}
    const VECOFVAR = MOI.VectorOfVariables

    const SQT = MOI.ScalarQuadraticTerm{Float64}
    const SAT = MOI.ScalarAffineTerm{Float64}

    const LT = MOI.LessThan{Float64}
    const GT = MOI.GreaterThan{Float64}
    const ET = MOI.EqualTo{Float64}
    const IT = MOI.Interval{Float64}
    const ZO = MOI.ZeroOne
    const SOC = MOI.SecondOrderCone

    const VI = MOI.VariableIndex
    const CI = MOI.ConstraintIndex

    const SCoefC = MOI.ScalarCoefficientChange
    const SConsC = MOI.ScalarConstantChange

    const LT_ZERO = LT(0.0)

    # Add storage types for EAGO optimizers
    export NodeBB, get_history, get_lower_bound, get_upper_bound, get_lower_time,
           get_upper_time, get_preprocess_time, get_postprocess_time, get_lower_bound, get_solution_time,
           get_iteration_number, get_node_count, get_absolute_gap, get_relative_gap, SubSolvers,
           EAGOModel, AuxiliaryVariableRef
    export auxillary_variable, @auxillary_variable, add_mimo_expression, @add_mimo_expression

    export register_eago_operators!

    # map/reduce nonallocating no bounds checking map-reduce like utilities
    include(joinpath(@__DIR__, "eago_optimizer", "utilities.jl"))

    # creates a context that removes domain violations when constructing bounds
    #include("eago_optimizer/guarded_context.jl")

    include(joinpath(@__DIR__, "eago_optimizer", "types", "log.jl"))
    include(joinpath(@__DIR__, "eago_optimizer", "types", "variable_info.jl"))
    include(joinpath(@__DIR__, "eago_optimizer", "types", "node_bb.jl"))
    include(joinpath(@__DIR__, "eago_optimizer", "types", "extension.jl"))
    include(joinpath(@__DIR__, "eago_optimizer", "types", "incremental.jl"))
    include(joinpath(@__DIR__, "eago_optimizer", "types", "subsolver_block.jl"))

    # load internal storage functions
    include(joinpath(@__DIR__, "eago_optimizer", "functions", "functions.jl"))

    include(joinpath(@__DIR__, "eago_optimizer", "types", "global_optimizer.jl"))

    # defines the optimizer structures
    include(joinpath(@__DIR__, "eago_optimizer", "optimizer.jl"))

    # defines routines to add variable, saf, sqf, and nlp block constraints
    include(joinpath(@__DIR__, "eago_optimizer", "moi_wrapper.jl"))

    #
    include(joinpath(@__DIR__, "eago_optimizer", "optimize", "nonconvex", "relax.jl"))
    include(joinpath(@__DIR__, "eago_optimizer", "optimize", "nonconvex", "bound.jl"))

    #
    include(joinpath(@__DIR__, "eago_optimizer", "domain_reduction.jl"))

    #
    include(joinpath(@__DIR__, "eago_optimizer", "parse.jl"))

    #
    include(joinpath(@__DIR__, "eago_optimizer", "optimize", "optimize.jl"))

    # import the script solving utilities
    include(joinpath(@__DIR__, "eago_script", "script.jl"))

    # routines for solving SIPs
    export SIPResult, SIPProblem, SIPCallback, SIPSubResult,
           sip_solve, SIPRes, SIPResRev, SIPHybrid,
           build_model, set_tolerance_inner!, set_tolerance!, get_disc_set,
           sip_llp!, sip_bnd!, sip_res!, get_sip_optimizer, check_convergence,
           LowerLevel1, LowerLevel2, LowerLevel3, LowerProblem, UpperProblem,
           ResProblem, AbstractSIPAlgo, AbstractSubproblemType
    include(joinpath(@__DIR__, "eago_semiinfinite", "semiinfinite.jl"))

    include(joinpath(@__DIR__, "subsolvers", "cbc.jl"))
    include(joinpath(@__DIR__, "subsolvers", "ipopt.jl"))
    function __init__()
        #@require Cbc="9961bab8-2fa3-5c5a-9d89-47fab24efd76"        include(joinpath(@__DIR__, "subsolvers", "cbc.jl"))
        @require Clp="e2554f3b-3117-50c0-817c-e040a3ddf72d"        include(joinpath(@__DIR__, "subsolvers", "clp.jl"))
        @require CPLEX="a076750e-1247-5638-91d2-ce28b192dca0"      include(joinpath(@__DIR__, "subsolvers", "cplex.jl"))
        @require GLPK="60bf3e95-4087-53dc-ae20-288a0d20c6a6"       include(joinpath(@__DIR__, "subsolvers", "glpk.jl"))
        @require Gurobi="2e9cd046-0924-5485-92f1-d5272153d98b"     include(joinpath(@__DIR__, "subsolvers", "gurobi.jl"))
        @require Hypatia="b99e6be6-89ff-11e8-14f8-45c827f4f8f2"    include(joinpath(@__DIR__, "subsolvers", "hypatia.jl"))
        @require KNITRO="67920dd8-b58e-52a8-8622-53c4cffbe346"     include(joinpath(@__DIR__, "subsolvers", "knitro.jl"))
        @require MosekTools="1ec41992-ff65-5c91-ac43-2df89e9693a4" include(joinpath(@__DIR__, "subsolvers", "mosek.jl"))
        @require Xpress="9e70acf3-d6c9-5be6-b5bd-4e2c73e3e054"     include(joinpath(@__DIR__, "subsolvers", "xpress.jl"))
    end

    #include("precompile.jl")
    #_precompile_()
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function RandomDataGenerator()
    return RandomDataGenerator(())
end

function RandomDataGenerator(arg0::jlong)
    return RandomDataGenerator((jlong,), arg0)
end

function next_beta(obj::RandomDataGenerator, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "nextBeta", jdouble, (jdouble, jdouble), arg0, arg1)
end

function next_deviate(obj::RandomDataGenerator, arg0::IntegerDistribution)
    return jcall(obj, "nextDeviate", jint, (IntegerDistribution,), arg0)
end

function next_deviate(obj::RandomDataGenerator, arg0::RealDistribution)
    return jcall(obj, "nextDeviate", jdouble, (RealDistribution,), arg0)
end

function next_deviates(obj::RandomDataGenerator, arg0::IntegerDistribution, arg1::jint)
    return jcall(obj, "nextDeviates", Vector{jint}, (IntegerDistribution, jint), arg0, arg1)
end

function next_deviates(obj::RandomDataGenerator, arg0::RealDistribution, arg1::jint)
    return jcall(obj, "nextDeviates", Vector{jdouble}, (RealDistribution, jint), arg0, arg1)
end

function next_exponential(obj::RandomDataGenerator, arg0::jdouble)
    return jcall(obj, "nextExponential", jdouble, (jdouble,), arg0)
end

function next_gamma(obj::RandomDataGenerator, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "nextGamma", jdouble, (jdouble, jdouble), arg0, arg1)
end

function next_hex_string(obj::RandomDataGenerator, arg0::jint)
    return jcall(obj, "nextHexString", JString, (jint,), arg0)
end

function next_int(obj::RandomDataGenerator, arg0::jint, arg1::jint)
    return jcall(obj, "nextInt", jint, (jint, jint), arg0, arg1)
end

function next_log_normal(obj::RandomDataGenerator, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "nextLogNormal", jdouble, (jdouble, jdouble), arg0, arg1)
end

function next_long(obj::RandomDataGenerator, arg0::jlong, arg1::jlong)
    return jcall(obj, "nextLong", jlong, (jlong, jlong), arg0, arg1)
end

function next_normal(obj::RandomDataGenerator, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "nextNormal", jdouble, (jdouble, jdouble), arg0, arg1)
end

function next_permutation(obj::RandomDataGenerator, arg0::jint, arg1::jint)
    return jcall(obj, "nextPermutation", Vector{jint}, (jint, jint), arg0, arg1)
end

function next_poisson(obj::RandomDataGenerator, arg0::jdouble)
    return jcall(obj, "nextPoisson", jint, (jdouble,), arg0)
end

function next_sample(obj::RandomDataGenerator, arg0::Collection, arg1::jint)
    return jcall(obj, "nextSample", Vector{Object}, (Collection, jint), arg0, arg1)
end

function next_sample(obj::RandomDataGenerator, arg0::Vector{jdouble}, arg1::jint)
    return jcall(obj, "nextSample", Vector{jdouble}, (Vector{jdouble}, jint), arg0, arg1)
end

function next_sample_with_replacement(obj::RandomDataGenerator, arg0::jint, arg1::Vector{jdouble})
    return jcall(obj, "nextSampleWithReplacement", Vector{jint}, (jint, Vector{jdouble}), arg0, arg1)
end

function next_uniform(obj::RandomDataGenerator, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "nextUniform", jdouble, (jdouble, jdouble), arg0, arg1)
end

function next_zipf(obj::RandomDataGenerator, arg0::jint, arg1::jdouble)
    return jcall(obj, "nextZipf", jint, (jint, jdouble), arg0, arg1)
end

function of(::Type{RandomDataGenerator}, arg0::RandomGenerator)
    return jcall(RandomDataGenerator, "of", RandomDataGenerator, (RandomGenerator,), arg0)
end


text:
function generateBaseNodeList( mpSim::MPsim, nodes::Vector{String} )

    isBaseNode = haskey.( Ref( mpSim.baseNodeList ), nodes )
    nodeList = nodes[isBaseNode]

    for node in nodes[.!isBaseNode]
        append!( nodeList, mpSim.compoundNodeList[node].baseNodeList )
    end  # for node in nodes[.!isBaseNode]

    return unique( nodeList )

end  # generateBaseNodeList( mpSim, nodes )

function generateCompositionReport( node::CompoundNode,
    baseNodeReport::DataFrame )

    result = baseNodeReport[:, vcat( :timePoint,
        Symbol.( node.baseNodeList ) )]
    nRows = size( result, 1 )
    counts = zeros( Int, nRows )

    for baseNode in node.baseNodeList
        counts += result[:, Symbol( baseNode )]
    end  # for baseNode in node.baseNodeList

    insertcols!( result, 2, Symbol( node.name ) => counts )

    return result

end  # generateCompositionReport( node, baseNodeReport )


text:
using NExOS
using ProximalOperators
using Test
using Random, LinearAlgebra

@testset "sparse_regression" begin

    M = 1.0

    r = 3

    C = SparseSet(M, r)

    A = [0.0026751207469390887 -1.7874010558441433 1.9570015234024314 0.48297821744005304 0.5499471860787609; -1.1926192117844148 -0.2962711608418759 0.34270253049651844 -0.9366799738457191 2.6135543879953094; -2.3956100254067567 1.5957363630260428 -0.5408329087542932 -1.595958552841769 -1.6414598361443185; 0.6791192541307551 2.895535427621422 -0.6676549207144326 -0.49424052539586016 -0.4336822064078515; 2.201382399088591 -0.578227952607833 -0.17602060199064307 -0.46261612338452723 0.24784285430931077; -0.49522601422578916 0.44636764792131967 1.2906418721153354 0.8301596930345838 0.5837176074011505]

    b = [-0.7765194270735608, 0.28815810383363427, -0.07926006593887291, 0.5659412043036721, 1.3718921913877151, 0.8834692513807884]

    m, n = size(A)

    z0 = [0.7363657041121843, 1.7271074489333478, -1.7929456862340958, 0.9359690211495514, 1.771070535933835]

    f = LeastSquares(A, b, iterative = true)

    settings = NExOS.Settings(μ_max = 2, μ_min = 1e-8, μ_mult_fact = 0.85, verbose = true, freq = 250, γ_updt_rule = :adaptive)

    problem = NExOS.Problem(f, C, settings.β, z0)

    state_final = NExOS.solve!(problem, settings)

    @test log10(state_final.fxd_pnt_gap) <= -4

    @test log10(state_final.fsblt_gap) <= -4

    @test abs(f(state_final.x) - 0.934) <= 1e-3

end

@testset "affine_rank_minimization" begin

    m = 10

    n = 2*m

    M = 1.0

    k = convert(Int64, m*n/20) # k is the number of components in the affine operator A: R^m×n → R^k

    r = convert(Int64,round(m*.35))  # r corresponds to the rank of the matrix rank(X) <= r

    mcA = randn(k, m*n)

    b = randn(k)

    Z0 = zeros(m,n)

    f = LeastSquaresOverMatrix(mcA, b, 1.0, iterative = true);

    D = NExOS.RankSet(M, r)

    settings = NExOS.Settings(μ_max = 2, μ_min = 1e-8, μ_mult_fact = 0.5, verbose = true, freq = 250, γ_updt_rule = :safe)

    problem = NExOS.Problem(f, D, settings.β, Z0)

    state_final = NExOS.solve!(problem, settings)

    f(state_final.x)

    @test log10(state_final.fxd_pnt_gap) <= -4

    @test log10(state_final.fsblt_gap) <= -4

end

@testset "matrix_completion" begin

    Random.seed!(1234)

    # Construct a random m-by-n matrix matrix of rank k
    m,n,k = 5,5,2

    Zfull = randn(m,k)*randn(k,n) # ground truth data

    M = opnorm(Zfull,2)

    n_obs = 13

    Zobs = fill(NaN,(m,n))

    obs = randperm(m*n)[1:n_obs]

    Zobs[obs] = Zfull[obs] # partially observed matrix

    # let us find the indices of all the elements that are available

    f = SquaredLossMatrixCompletion(Zobs, iterative = true)

    r = k

    Z0 = zeros(size(Zobs))

    C = RankSet(M, r)

    settings = NExOS.Settings(μ_max = 5, μ_min = 1e-8, μ_mult_fact = 0.5, n_iter_min = 1000, n_iter_max = 1000, verbose = true, freq = 250, tol = 1e-4, γ_updt_rule = :safe)

    problem = NExOS.Problem(f, C, settings.β, Z0)

    state_final = NExOS.solve!(problem, settings)

    @test log10(state_final.fxd_pnt_gap) <= -4

    @test log10(state_final.fsblt_gap) <= -4

end


text:
#### Test BlockMap with SerialComm ####

function SerialMapTests(map::BlockMap{Int, Int, Int}, map2::BlockMap{Int, Int, Int}, diffMap::BlockMap{Int, Int, Int})
#    quote
        mapCopy = BlockMap{Int, Int, Int}(map.data)

        @test uniqueGIDs(map)

        for i = 1:5
            @test myLID(map, i)
            @test myGID(map, i)
            @test i == lid(map, i)
            @test i == gid(map, i)
        end

        @test !myLID(map, -1)
        @test !myLID(map, 0)
        @test !myLID(map, 6)
        @test !myLID(map, 30)
        @test !myGID(map, -1)
        @test !myGID(map, 0)
        @test !myGID(map, 6)
        @test !myGID(map, 30)

        @test 0 == lid(map, -1)
        @test 0 == lid(map, 0)
        @test 0 == lid(map, 6)
        @test 0 == lid(map, 30)
        @test 0 == gid(map, -1)
        @test 0 == gid(map, 0)
        @test 0 == gid(map, 6)
        @test 0 == gid(map, 30)

        @test !distributedGlobal(map)

        @test 5 == numGlobalElements(map)
        @test 5 == numMyElements(map)

        @test 1 == minMyGID(map)
        @test 5 == maxMyGID(map)
        @test 1 == minAllGID(map)
        @test 5 == maxAllGID(map)
        @test 1 == minLID(map)
        @test 5 == maxLID(map)

        @test ([1, 1, 1, 1, 1], [1, 2, 3, 4, 5]) == remoteIDList(map, [1, 2, 3, 4, 5])

        @test [1, 2, 3, 4, 5] == myGlobalElements(map)

        @test sameBlockMapDataAs(map, mapCopy)
        @test sameBlockMapDataAs(mapCopy, map)
        @test !sameBlockMapDataAs(map, map2)
        @test !sameBlockMapDataAs(map2, map)
        @test !sameBlockMapDataAs(map, diffMap)
        @test !sameBlockMapDataAs(diffMap, map)

        @test sameAs(map, mapCopy)
        @test sameAs(mapCopy, map)
        @test sameAs(map, map2)
        @test sameAs(map2, map)
        @test !sameAs(map, diffMap)
        @test !sameAs(diffMap, map)
        @test !sameAs(map2, diffMap)
        @test !sameAs(diffMap, map2)

        @test linearMap(map)

        @test [1, 2, 3, 4, 5] == myGlobalElementIDs(map)

        @test commVal == getComm(map)
#    end
end

commVal = SerialComm{Int, Int, Int}()

## constructor 1 ##
@test_throws InvalidArgumentError BlockMap(-8, commVal)
@test_throws InvalidArgumentError BlockMap(-1, commVal)

BlockMap(0, commVal)
BlockMap(1, commVal)

map = BlockMap(5, commVal)
map2 = BlockMap(5, commVal)
diffMap = BlockMap(6, commVal)
#@SerialMapTests
SerialMapTests(map, map2, diffMap)

## constructor 2 ##
@test_throws InvalidArgumentError BlockMap(-8, 4, commVal)
@test_throws InvalidArgumentError BlockMap(-2, 4, commVal)
@test_throws InvalidArgumentError BlockMap(5, -6, commVal)
@test_throws InvalidArgumentError BlockMap(4, -1, commVal)

BlockMap(0, 0, commVal)
BlockMap(1, 1, commVal)

map = BlockMap(5, 5, commVal)
map2 = BlockMap(5, 5, commVal)
diffMap = BlockMap(6, 6, commVal)
#@SerialMapTests
SerialMapTests(map, map2, diffMap)

map = BlockMap(-1, 5, commVal)
map2 = BlockMap(-1, 5, commVal)
diffMap = BlockMap(-1, 6, commVal)
#@SerialMapTests
SerialMapTests(map, map2, diffMap)

## constructor 3 ##
BlockMap(0, Int[], commVal)
BlockMap(1, [1], commVal)

map = BlockMap(5, [1, 2, 3, 4, 5], commVal)
map2 = BlockMap(5, [1, 2, 3, 4, 5], commVal)
diffMap = BlockMap(6, [1, 2, 3, 4, 5, 6], commVal)
#@SerialMapTests
SerialMapTests(map, map2, diffMap)

## constructor 4 ##
@test_throws InvalidArgumentError BlockMap(-8, 4, [1, 2, 3, 4], false, 1, 4, commVal)
@test_throws InvalidArgumentError BlockMap(-2, 4, [1, 2, 3, 4], false, 1, 4, commVal)
@test_throws InvalidArgumentError BlockMap(5, -6, [1, 2, 3, 4, 5], false, 1, 5, commVal)
@test_throws InvalidArgumentError BlockMap(4, -1, [1, 2, 3, 4], false, 1, 4, commVal)

BlockMap(0, 0, Int[], false, 1, 0, commVal)
BlockMap(1, 1, [1], false, 1, 1, commVal)

map = BlockMap(5, 5, [1, 2, 3, 4, 5], false, 1, 5, commVal)
map2 = BlockMap(5, 5, [1, 2, 3, 4, 5], false, 1, 5, commVal)
diffMap = BlockMap(6, 6, [1, 2, 3, 4, 5, 6], false, 1, 6, commVal)
#@SerialMapTests
SerialMapTests(map, map2, diffMap)

#stability tests
for (GID, PID, LID) in Base.product(stableGIDs, stablePIDs, stableLIDs)
        @test is_stable(check_method(gid, (BlockMap{GID, PID, LID}, LID)))
        @test is_stable(check_method(lid, (BlockMap{GID, PID, LID}, GID)))
        @test is_stable(check_method(myGID, (BlockMap{GID, PID, LID}, GID)))
        @test is_stable(check_method(myLID, (BlockMap{GID, PID, LID}, LID)))
        @test is_stable(check_method(numGlobalElements, (BlockMap{GID, PID, LID},)))
        @test is_stable(check_method(numMyElements, (BlockMap{GID, PID, LID},)))
end


text:
"""
    MDCProblem(cost, p0, dp0, momentum, tspan)
Creates an MDCProblem, that can then generate a minimally disruptive curve using evolve(c::MDCProblem, ...; ...)
"""
struct MDCProblem{A,B,C,D,E} <: CurveProblem
    cost::A
    p0::B
    dp0::C
    momentum::D
    tspan::E
    ## reverse initial direction and signflip curve span if the latter is nonpositive
    function MDCProblem(a::A, b::B, c::C, d::D, e::E) where A where B where C where D where E
        if max(e...) <= 0.
            e = map(x -> -x |> abs, e) |> reverse
            c = -c
        end
        new{A,B,C,D,E}(a, b, c, d, e)
    end
end
isjumped(c::MDCProblem) = ZeroStart()
whatdynamics(c::MDCProblem) = MDCDynamics()

num_params(c::CurveProblem) = length(c.p0)
param_template(c::CurveProblem) = deepcopy(c.p0)
initial_params(c::CurveProblem) = c.p0

"""
DEPRECATED. use MDCProblem
"""
function curveProblem(a, b, c, d, e)
    @warn("curveProblem and specify_curve are DEPRECATED. please use MDCProblem (with the same arguments) instead")
    return MDCProblem(a, b, c, d, e)
end

"""
DEPRECATED. use MDCproblem
"""
specify_curve(cost, p0, dp0, momentum, tspan) = curveProblem(cost, p0, dp0, momentum, tspan)
specify_curve(;cost=nothing, p0=nothing, dp0=nothing,momentum=nothing,tspan=nothing) = curveProblem(cost, p0, dp0, momentum, tspan)

"""
    Callback to readjust momentum in the case that the numerical residual from the identity dHdu = 0 crosses a user-specified threshold
"""
(m::MomentumReadjustment)(c::CurveProblem) = readjustment(c, ResidualCondition(), CostateAffect(), m.tol, m.verbose)
"""
    Callback to readjust state in the case that the numerical residual from the identity dHdu = 0 crosses a user-specified threshold. EXPERIMENTAL AND WILL PROBABLY BREAK
"""
(m::StateReadjustment)(c::CurveProblem) = readjustment(c, ResidualCondition(), StateAffect(), m.tol, m.verbose)

"""
    (c::MDCProblem)()
returns a tuple of ODEProblems specificed by the MDCProblem. Usually a single ODEProblem. Two are provided if the curve crosses zero, so that one can run two curves in parallel going backwards/forwards from zero
"""
function (c::MDCProblem)()
    spans = make_spans(c, c.tspan)
    cs = map(spans) do span
        mult = sign(span[end])
        return MDCProblem(c.cost, c.p0, mult * c.dp0, c.momentum, abs.(span))
    end
    spans = map(x -> abs.(x), spans)
    u0s = initial_conditions.(cs)
    u0 = map(span -> initial_conditions(c), spans)
    fs = dynamics.(cs)
    return ODEProblem.(fs, u0s, spans)
    # return map(sp -> ODEProblem(f, u0, sp), spans)  # make two problems for 2-sided tspan
end

"""
    make_spans(c::MDCProblem, span)
- makes sure span of curve is increasing.
- if the span crosses zero, then returns two separate spans. evolve then runs two curves in parallel, going backwards/forwards from zero.
"""
function make_spans(c::CurveProblem, span, ::ZeroStart)
    (span[1] > span[2]) && error("make your curve span monotone increasing")
    if (span[2] > 0) && (span[1] < 0)
        spans = ((0., span[1]), (0., span[2]))
    else
        spans = (span,)
    end
    return spans
end

"""
    initial_costate(c::MDCProblem)
solves for the initial costate required to evolve a MD curve.
"""
function initial_costate(c::MDCProblem)
    μ₂ = (-c.momentum + c.cost(c.p0)) / 2.
    λ₀ = -2. * μ₂ * c.dp0
    return λ₀
end

"""
Generate initial conditions of an MDCProblem
"""
function initial_conditions(c::MDCProblem)
    λ₀ = initial_costate(c)
    return cat(c.p0, λ₀, dims=1)
end

"""
Generate vector field for MD curve, as specified by MDCProblem
"""
function dynamics(c::CurveProblem, ::MDCDynamics)
    cost = c.cost
    ∇C = param_template(c)
    N = num_params(c)
    H = c.momentum
    θ₀ = initial_params(c)
    function upd(du, u, p, t)
        θ = u[1:N] # current parameter vector
        λ = u[N + 1:end] # current costate vector
        dist = sum((θ - θ₀).^2) # which should = t so investigate cost/benefits of using t instead of dist
        C = cost(θ, ∇C) # also updates ∇C as a mutable
        μ2 = (C - H) / 2
        μ1 = dist > 1e-3 ?  (λ' * λ - 4 * μ2^2 ) / (λ' * (θ - θ₀)) : 0.
        # if mu1 < -1e-4 warn of numerical issue
        # if mu1 > 1e-3 and dist > 1e-3 then set mu1 = 0
        du[1:N] = @. (-λ + μ1 * (θ - θ₀)) / (2 * μ2) # ie dθ
        du[1:N] /= (sqrt(sum((du[1:N]).^2)))
        damping_constant = (λ' * du[1:N]) / (H - C)  # theoretically = 1 but not numerically
        du[N + 1:end] = @. (μ1 * du[1:N] - ∇C) * damping_constant # ie dλ
res = λ  + 2 * μ2 * du[1:N]
        return nothing
    end
    return upd
    end

"""
Callback to stop MD Curve evolving if cost > momentum
"""
function (t::TerminalCond)(c::CurveProblem)
    cost = c.cost
    H = c.momentum
    N = num_params(c)
    function condition(u, t, integrator)
        return (cost(u[1:N]) > H)
    end
    return DiscreteCallback(condition, terminate!)
end

function readjustment(c::CurveProblem, cnd::ConditionType, aff::AffectType, momentum_tol, verbose::Bool)
    if isnan(momentum_tol)
        return nothing
    end
    cond = build_cond(c, cnd, momentum_tol)
    affect! = build_affect(c, aff)
    return DiscreteCallback(cond, affect!)
end

    function build_cond(c::CurveProblem, ::ResidualCondition, tol, ::MDCDynamics)
        function rescond(u, t, integ)
            absres = dHdu_residual(c, u, t, integ)
            absres > tol ? begin
            # @info "applying readjustment at t=$t, |res| = $absres"
                return true
            end : return false
        end
    return rescond
end

function build_cond(c::CurveProblem, ::CostCondition, tol)
    N = num_params(c)
    function costcond(u, t, integ)
         (c.cost(u[1:N]) > tol) ? (return true) : (return false)
     end
     return costcond
end

"""
    For dHdu_residual and build_affect(::MDCProblem, ::CostateAffect): there is an unnecessary allocation in the line `dθ = ...`. I initially used `dθ[:] = ....`, but this produced unreliable output (the MDCurve changed on each solution). I found that this was because temporary arrays like this are not safe in callbacks, for some reason. The solution is to use SciMLBase.get_tmp_cache. Don't have time to figure out how to do this right now. Do at some point.
"""

"""
    dHdu_residual(c::MDCProblem, u, t, dθ)
Checks dHdu residual (u deriv of Hamiltonian). Returns true if abs(residual) is greater than some tolerance (it should be zero)
"""
function dHdu_residual(c::CurveProblem, u, t, integ, ::MDCDynamics)
    N = num_params(c)
    H = c.momentum
    θ₀ = initial_params(c)

    θ = u[1:N]
    λ = u[N + 1:end]
    μ2 = (c.cost(θ) - H) / 2.
    μ1 = t > 1e-3 ?  (λ' * λ - 4 * μ2^2 ) / (λ' * (θ - θ₀)) : 0.

    # dθ = SciMLBase.get_tmp_cache(integ)[1][1:N]
    dθ = (-λ + μ1 * (θ - θ₀)) / (2 * μ2)
    dθ /= (sqrt(sum((dθ).^2)))
    return sum(abs.(λ + 2 * μ2 * dθ))
end

"""
    build_affect(c::MDCProblem, ::CostateAffect)
Resets costate to undo effect of cumulative numerical error. Specifically, finds costate so that dHdu = 0, where H is the Hamiltonian.

*I wanted to put dθ[:] = ... here instead of dθ = ... . Somehow the output of the MDC changes each time if I do that, there is a dirty state being transmitted. But I don't at all see how from the code. Figure out.*
"""
function build_affect(c::CurveProblem, ::CostateAffect, ::MDCDynamics)
    N = num_params(c)
    H = c.momentum
    θ₀ = initial_params(c)
    dp = param_template(c)
    function reset_costate!(integ, dθ)
        θ = integ.u[1:N] # current parameter vector
        λ = integ.u[N + 1:end] # current costate vector
        μ2 = (c.cost(θ) - H) / 2
        μ1 = integ.t > 1e-3 ?  (λ' * λ - 4 * μ2^2 ) / (λ' * (θ - θ₀)) : 0.
        # dθ = SciMLBase.get_tmp_cache(integ)[1][1:N]
        dθ = (-λ + μ1 * (θ - θ₀)) / (2 * μ2)
        dθ /= (sqrt(sum((dθ).^2)))
        integ.u[N + 1:end] =  -2 * μ2 * dθ
        return integ
    end
    return integ -> reset_costate!(integ, dp)
end

"""
    build_affect(c::MDCProblem, ::StateAffect)
resets state so that residual is zero. also resets costate necessarily. NOT YET FULLY IMPLEMENTED
min C(θ) such that norm(θ - θ₀)^2 = K where K is current distance
we will do this with unconstrained optimisation and lagrange multipliers
ideally would have an inequality constraint >=K. But Optim.jl doesn't support this
"""
function build_affect(c::CurveProblem, ::StateAffect, ::MDCDynamics)
    N = num_params(c)
    H = c.momentum
    cost = c.cost
    θ₀ = initial_params(c)
    dp = param_template(c)
    _reset_costate! = build_affect(c, CostateAffect())
    function reset_state!(integ, dθ)
            K = sum((integ.u[1:N] - θ₀).^2)
        println(K)
        function constr(x) # constraint func: g = 0
            return K - sum((x - θ₀).^2)
        end

        function L(x)
            θ, λ = x[1:end - 1], x[end]
            return cost(θ) + λ * constr(θ)
        end
        gc = deepcopy(θ₀)
        function L(x, g)
            θ, λ = x[1:end - 1], x[end]
            C = cost(θ, dθ) # dθ is just an arbitrary pre-allocation
            cstr = constr(θ)
            g[1:end - 1] = gc + 2 * λ * (θ - θ₀)
            g[end] = cstr
            return C + λ * cstr
        end
        g = zeros(N + 1)
    opt = optimize(L, cat(integ.u[1:N], 0., dims=1), LBFGS())
        C0 = cost(θ₀)
        @info "cost after readjustment is $(opt.minimum). cost before readjustment was $C0"
        (opt.ls_success == true) && (integ.u[1:N] = opt.minimizer[1:N])
        integ = _reset_costate!(integ, dθ)
        return integ
end
return integ -> reset_state!(integ, dp)
end

function build_callbacks(c::CurveProblem, callbacks::SciMLBase.DECallback)
    # DECallback supertype includes CallbackSet
    return CallbackSet(callback)
    end

build_callbacks(c::CurveProblem, n::Nothing) = nothing

function build_callbacks(c::CurveProblem, mdc_callbacks::Vector{T}, mtol::Number) where T <: CallbackCallable

    if !any(x -> typeof(x) <: MomentumReadjustment, mdc_callbacks)
        push!(mdc_callbacks, MomentumReadjustment(mtol))
    end
    push!(mdc_callbacks, TerminalCond())
    actual_callbacks = map(mdc_callbacks) do cb
        a = cb(c)
    end |> x -> vcat(x...)
    return actual_callbacks
end


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

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

function _swaggerinternal_appsPermissionsUsersList(_api::AppsPermissionsUsersApi; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/apps.permissions.users.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Returns list of user grants and corresponding scopes this app has on a team.
Param: cursor::String
Param: limit::Int32
Return: Dict{String, Any}
"""
function appsPermissionsUsersList(_api::AppsPermissionsUsersApi; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsUsersList(_api; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function appsPermissionsUsersList(_api::AppsPermissionsUsersApi, response_stream::Channel; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsUsersList(_api; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_appsPermissionsUsersRequest(_api::AppsPermissionsUsersApi, scopes::String, trigger_id::String, user::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/apps.permissions.users.request", ["slackAuth"])
    Swagger.set_param(_ctx.query, "scopes", scopes)  # type String
    Swagger.set_param(_ctx.query, "trigger_id", trigger_id)  # type String
    Swagger.set_param(_ctx.query, "user", user)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Enables an app to trigger a permissions modal to grant an app access to a user access scope.
Param: scopes::String (required)
Param: trigger_id::String (required)
Param: user::String (required)
Return: Dict{String, Any}
"""
function appsPermissionsUsersRequest(_api::AppsPermissionsUsersApi, scopes::String, trigger_id::String, user::String; _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsUsersRequest(_api, scopes, trigger_id, user; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function appsPermissionsUsersRequest(_api::AppsPermissionsUsersApi, response_stream::Channel, scopes::String, trigger_id::String, user::String; _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsUsersRequest(_api, scopes, trigger_id, user; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export appsPermissionsUsersList, appsPermissionsUsersRequest


text:
# Stricker Klaus 2020 - License is MIT: http://julialang.org/license
# ReadWriteDlm2 - rwd2_writeutils.jl - https://github.com/strickek/ReadWriteDlm2.jl

"""

    floatformat(a, decimal::AbstractChar)

Convert Int or Float64 numbers to String and change decimal mark.
"""
function floatformat(a, decimal)
    a = string(a)
    (decimal != '.') && (a = replace(a, '.' => decimal))
    return a
end

"""

    timeformat(a, decimal::AbstractChar)

Convert Time to String, optional with change of decimal mark for seconds.
"""
function timeformat(a, decimal)
    a = string(a)
    (decimal != '.') && (a = replace(a, '.' => decimal))
    return a
end

"""

    Complexformat(a, decimal::AbstractChar, imsuffix::AbstractString)

Convert Complex number to String, optional change of decimal and/or imsuffix.
"""
function complexformat(a, decimal, imsuffix)
    a = string(a)
    a = replace(a, " " => "" )  #"1 + 3im" => "1+3im"
    (imsuffix != "im") && (a = string(split(a, "im")[1], imsuffix))
    (decimal != '.') && (a = replace(a, '.' => decimal))
    return a
end


text:
#=
Offset

A point in an array is located by its index in the frame of the array.

An offset is a vector, in another, amorphous frame
(whose origin is the other point in the other frame.)

This hides that indices may be linear or cartesian.
Whereas an offset may be of the opposite variety.
Hides in the sense: I went back and forth choosing types
in the algorithm when testing 1D, 2D,...

In Resynthesizer, offsetting points is done:
1) when choosing neighbors (ScatterPatch is a set of nearest offset points.)
2) in the heuristic search, inverse offsetting back into the target.

TODO when sampling random points in the corpus,
we could choose offsets at random and use this.
=#

#=
Return a point (index) offset from given point.

Given point is in bounds of array.

Returned point may be out of bounds.
Returned point is in frame of array.

IOW this is coordinate arithmetic.
=#

#=
Specific to case where index and offset are same type, usually linear index i.e. Int64,
or both CartesianIndex{N}
TODO cruft?
specific to point is-a Int64 and offset is-a CartesianIndex???
=#
function offsetPatchPoint(
        array::Array{ValueType, DimensionCount},    # not used
        point, # ::CartesianIndex,
        offset::CartesianIndex{DimensionCount}
        ) where{ValueType, DimensionCount}
    return point + offset
end

#=
Specific to case where index is linear index i.e. Int64
but offset is CartesianIndex
TODO maybe we should generate offsets that are always linear indices???
=#
function offsetPatchPoint(
        array::Array{ValueType, DimensionCount},
        point::Int64,   # linear index
        offset::CartesianIndex{DimensionCount}) where{ValueType, DimensionCount}
    #=
    Conversion from linear index to Cartesian is expensive division.
    =#
    return CartesianIndices(array)[point] + offset
end

#= CRUFT
Formerly called from ScatterPatch

function pointAtOffset(pointInImage, offset)
    #= OLD, framed
    return PointInMaskedImage(
        pointInImage.point + offset,    # unchecked coordinate arithmetic
        pointInImage.maskedImage
    )
    =#
    return pointInImage + offset # unchecked coordinate arithmetic

    #result = copy(pointInImage)
    #result.point = pointInImage.point + offset
    #return result
end
=#

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
"""
    NexRTC(nexSerial, name; pid=pid, cid=cid)

A Nextion control to manage Real Time Clock (RTC).

!!! note

    Only for enhanced model.
"""
struct NexRTC <: AbstractNexObject
    _nid::NexID

    function NexRTC(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid)
    end
end


text:
module Foos

import SearchLight: AbstractModel, DbId
import Base: @kwdef

export Foo

@kwdef mutable struct Foo <: AbstractModel
  id::DbId = DbId()
end

end


text:
const MAX_UINT = 2^32 - 1

struct ScryptParameters
    r::Int  # element length multiplier
    N::Int  # processing cost
    p::Int  # parallelization

    function ScryptParameters(r, N, p)
        r > 0 || ArgumentError("r Must be > 0.") |> throw
        N > 0 || ArgumentError("N Must be > 0.") |> throw
        p > 0 || ArgumentError("p Must be > 0.") |> throw

        parameters = new(UInt(r), UInt(N), UInt(p))

        p ≤ MAX_UINT * hashlength(parameters) / elementlength(parameters) ||
            ArgumentError("p and r must satisfy the relationship p ≤ (2^32 - 1) * hashlength / elementlength)") |> throw

        r * N * elementunitlength(parameters) ≤ Sys.total_memory() ||
            ArgumentError("r and N must satisfy the relationship r * N * elementunitlength ≤ Sys.total_memory") |> throw

        parameters
    end
end

bytes(x) = x ÷ 8

hashbitslength(::ScryptParameters) = 256
hashlength(x::ScryptParameters) = hashbitslength(x) |> bytes

elementunitbitslength(::ScryptParameters) = 2 * 8 * sizeof(Salsa512)
elementunitlength(x::ScryptParameters) = elementunitbitslength(x) |> bytes
elementlength(x::ScryptParameters) = elementunitlength(x) * x.r
bufferlength(x::ScryptParameters) = elementlength(x) * x.p
elementblockcount(x::ScryptParameters) = elementlength(x) ÷ sizeof(Salsa512)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Libdl, PCRE2_jll

@testset "PCRE2_jll" begin
    vstr = zeros(UInt8, 32)
    @test ccall((:pcre2_config_8, libpcre2_8), Cint, (UInt32, Ref{UInt8}), 11, vstr) > 0
    vn = VersionNumber(split(unsafe_string(pointer(vstr)), " ")[1])
    @test vn == v"10.35.0"
end

max_stars_repo_name: dmageeLANL/julia

max_issues_repo_name: dmageeLANL/julia

max_forks_repo_name: dmageeLANL/julia


text:
using Printf

function string_examples()
    println("--------------------------------------")
    println("String examples")

    println("Strings have to be surrounded by double-quotes")
    s1 = "Just some random words\n"
    println(s1)

    println("Size of string")
    println(length(s1))

    println("The 1st index is at 1 (NOT ZERO INDEXED!)")
    println(s1[1])

    println("Get last character")
    println(s1[end])

    println("Get a range")
    println(s1[1:4])

    println("Concatenation")
    s2 = string("Yukiteru", " Amano")
    println(s2)

    println("You can also concatenate strings with the * symbol")
    println("Yuno" * " Gasai")

    println("Interpolation (like using {0} is C#)")
    i3 = 2
    i4 = 3
    println("$i3 + $i4 = $(i3 + i4)")

    println("Multiline strings are surrounded by triple-quotes. Note the output includes the new-lines!")
    s3 = """I
    have
    many
    lines"""
    println(s3)

    println("String comparisons with == > < !=")
    println("Takao" > "Hiyama")

    println("Find index for character (note the result is one-indexed, so index of 'i' in 'Keigo' will be 3 not 2)")
    println(findfirst(isequal('i'), "Keigo"))

    println("Find a substring")
    println(occursin("key", "monkey"))

end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
using Test

const SBP_TEST = get(ENV, "SBP_TEST", "all")

@time @testset "SummationByPartsOperators.jl tests" begin
    @time if SBP_TEST == "all" || SBP_TEST == "part1"
        @time @testset "Periodic Operators" begin include("periodic_operators_test.jl") end
        @time @testset "Non-Periodic Operators" begin include("SBP_operators_test.jl") end
        @time @testset "Dissipation Operators" begin include("dissipation_operators_test.jl") end
        @time @testset "Variable Coefficient Operators" begin include("var_coef_operators_test.jl") end
        @time @testset "Banded Matrices" begin include("banded_matrices_test.jl") end
        @time @testset "Fourier Operators" begin include("fourier_operators_test.jl") end
        @time @testset "Fourier Operators 2D" begin include("fourier_operators_2d_test.jl") end
        @time @testset "Legendre Operators" begin include("legendre_operators_test.jl") end
        @time @testset "Sum of Operators" begin include("sum_of_operators_test.jl") end
        @time @testset "Upwind Operators" begin include("upwind_operators_test.jl") end
        @time @testset "Special Matrix Types" begin include("special_matrix_types.jl") end
    end

    @time if SBP_TEST == "all" || SBP_TEST == "part2"
        @time @testset "Coupling" begin include("coupling_test.jl") end
        @time @testset "Conservation Laws" begin
            include("conservation_laws/burgers_test.jl")
            include("conservation_laws/cubic_test.jl")
            include("conservation_laws/variable_linear_advection_test.jl")
            include("conservation_laws/quartic_nonconvex_test.jl")
        end
        @time @testset "Second Order Equations" begin
            include("second_order_eqs/wave_eq_test.jl")
        end
    end
end


text:
export scatter1d

# depth-first scattering tree traversal
"""
    scatter1d_layer(S, i, j, nDone, x_hat, Transform, FilterBanks)

Apply one layer of the wavelet scattering transform, calling future layers recursively.

# Arguments
 - `S::Scattered`: Scattered object for holding scattering coefficients.
 - `i::Int64`: Index of the current scattering layer.
 - `j::Int64`: Max subsampling accepted for the current node
 - `x_hat::Vector{Complex{Float64}}`: Fourier transformed input data vector.
 - `Transform::ScatteringTransform1d`: ScatteringTransform1d object containing information about the total transform.
 - `FilterBanks::Vector{FilterBank1d}`: Filter banks for each layer of the scattering transform.
"""
function scatter1d_layer(S::Scattered, i::Int64, j::Int64, nDone::Vector{Int64}, x_hat::Vector{Complex{Float64}}, Transform::ScatteringTransform1d, FilterBanks::Vector{FilterBank1d}; subsample::Bool=true)
    for ψ in FilterBanks[i].Λ
        if ψ.j > j || subsample==false
            # convolve signal with morlet wavelet, transform back to the time domain, and take the modulus
            U_hat = fft(abs.(ifft(x_hat .* ψ.ψ)))

            # stop scattering when we reach the last layer
            if i < Transform.D
                scatter1d_layer(S, i+1, ψ.j, nDone, U_hat, Transform, FilterBanks, subsample=subsample)
            end

            # apply lowpass filter
            # result should be real-valued, so cast to real
            S.Coeff[i+1][:, nDone[i]] .= real.(ifft(U_hat .* FilterBanks[i].ϕ.ϕ))
            nDone[i]+=1
        end
    end
end

"""
    scatter1d(x, Transform, FilterBanks)

Calculate the 1d wavelet scattering transform.

# Arguments
 - `x::Vector{Float64}`: Time series.
 - `Transform::ScatteringTransform1d`: Scattering transform object.
 - `FilterBanks::Vector{FilterBank1d}`: Filter banks for each layer of the scattering transform.
"""
function scatter1d(x::Vector{Float64}, Transform::ScatteringTransform1d, FilterBanks::Vector{FilterBank1d}; subsample::Bool=true)
    @assert length(FilterBanks)==Transform.D "The number of filter banks should equal the scattering transform depth.
                                              Got $(length(FilterBanks)) filter banks and $(Transform.D) layers."

    # Transform the signal to the frequency domain
    x_hat = FFTW.fft(x)
    # create Scattered object to hold the scattering coefficients
    S = Scattered1d(FilterBanks, subsample=subsample)

    # zeroth order scattering coefficients are a simple averaging of the
    S.Coeff[1] .= real.(ifft(x_hat .* FilterBanks[1].ϕ.ϕ))
    nDone=ones(Int, Transform.D)
    # higher order scattering coefficients are computed recursively
    # seed j=0 to allow all filters for the first layer
    scatter1d_layer(S, 1, -1, nDone, x_hat, Transform, FilterBanks, subsample=subsample)

    return S
end

# breadth-first scattering tree traversal
# TODO: allow any number of layers with this approach
"""
    scatter1d(x, Transform, FilterBanks)

Calculate the 1d wavelet scattering transform.

# Arguments
 - `x::Vector{Float64}`: Time series.
 - `Transform::ScatteringTransform1d`: Scattering transform object.
 - `FilterBanks::Vector{FilterBank1d}`: Filter banks for each layer of the scattering transform.
"""
function scatter1d(x::Vector{Float64}, Transform::ScatteringTransform1d, FilterBanks::Vector{FilterBank1dBlock})
    @assert length(FilterBanks)==Transform.D "The number of filter banks should equal the scattering transform depth.
                                              Got $(length(FilterBanks)) filter banks and $(Transform.D) layers."
    @assert Transform.D <=2 "The number of layers should be 2 or fewer. Got $(Transform.D)"

    # Transform the signal to the frequency domain
    x_hat = FFTW.fft(x)
    # create Scattered object to hold the scattering coefficients
    S = Scattered1d(FilterBanks)

    # zeroth order scattering coefficients are a simple averaging of the
    S.Coeff[1] .= real.(ifft(x_hat .* FilterBanks[1].ϕ))

    # first order scattering coefficients
    U1 = Array{Float64,2}(undef, size(FilterBanks[1].Λ))
    U1 .= abs.(cwt(x, FilterBanks[1]))
    S.Coeff[2] .= real.(ifft(fft(U1, 1) .* FilterBanks[1].ϕ, 1))

    # second order scattering coefficients
    U2i = Array{Float64,2}(undef, size(FilterBanks[2].Λ))
    nDone=0
    for i=1:size(U1,2)
        U2i .= abs.(cwt(U1[:,i], FilterBanks[2]))
        S.Coeff[3][:, nDone+1:nDone+size(U2i,2)] .= real.(ifft(fft(U2i,1) .* FilterBanks[2].ϕ, 1))
        nDone+=size(U2i,2)
    end

    return S
end


text:
"""
Core genie configuration / settings functionality.
"""
module Configuration

"""
  const GENIE_VERSION

References the current Genie version number.
"""
const GENIE_VERSION = v"1.8"

import Logging
import Genie
import MbedTLS

export isdev, isprod, istest, env
export Settings, DEV, PROD, TEST

# app environments
const DEV   = "dev"
const PROD  = "prod"
const TEST  = "test"

haskey(ENV, "GENIE_ENV") || (ENV["GENIE_ENV"] = DEV)
haskey(ENV, "HOST") || (ENV["HOST"] = "127.0.0.1")

"""
    isdev()  :: Bool

Set of utility functions that return whether or not the current environment is development, production or testing.

# Examples
```julia
julia> Configuration.isdev()
true

julia> Configuration.isprod()
false
```
"""
isdev() :: Bool  = (Genie.config.app_env == DEV)

"""
    isprod() :: Bool

Set of utility functions that return whether or not the current environment is development, production or testing.

# Examples
```julia
julia> Configuration.isdev()
true

julia> Configuration.isprod()
false
```
"""
isprod():: Bool = (Genie.config.app_env == PROD)

"""
    istest() :: Bool

Set of utility functions that return whether or not the current environment is development, production or testing.

# Examples
```julia
julia> Configuration.isdev()
true

julia> Configuration.isprod()
false
```
"""
istest():: Bool = (Genie.config.app_env == TEST)

"""
    env() :: String

Returns the current Genie environment.

# Examples
```julia
julia> Configuration.env()
"dev"
```
"""
env() :: String = Genie.config.app_env

"""
    buildpath() :: String

Constructs the temp dir where Genie's view files are built.
"""
buildpath() :: String = Base.Filesystem.mktempdir(prefix = "jl_genie_build_")

"""
    mutable struct Settings

App configuration - sets up the app's defaults. Individual options are overwritten in the corresponding environment file.

# Arguments
- `server_port::Int`: the port for running the web server (default 8000)
- `server_host::String`: the host for running the web server (default "127.0.0.1")
- `server_document_root::String`: path to the document root (default "public/")
- `server_handle_static_files::Bool`: if `true`, Genie will also serve static files. In production, it is recommended to serve static files with a web server like Nginx.
- `server_signature::String`: Genie's signature used for tagging the HTTP responses. If empty, it will not be added.
- `app_env::String`: the environment in which the app is running (dev, test, or prod)
- `cors_headers::Dict{String,String}`: default `Access-Control-*` CORS settings
- `cors_allowed_origins::Vector{String}`: allowed origin hosts for CORS settings
- `cache_duraction::Int`: cache expiration time in seconds
- `log_level::Logging.LogLevel`: logging severity level
- `log_to_file::Bool`: if true, information will be logged to file besides REPL
- `assets_fingerprinted::Bool`: if true, asset fingerprinting is used in the asset pipeline
- `session_key_name::String`: the name of the session cookie
- `session_storage::Symbol`: the backend adapter for session storage (default File)
- `inflector_irregulars::Vector{Tuple{String,String}}`: additional irregular singular-plural forms to be used by the Inflector
- `run_as_server::Bool`: when true the server thread is launched synchronously to avoid that the script exits
- `websockets_server::Bool`: if true, the websocket server is also started together with the web server
- `html_parser_close_tag::String`: default " /". Can be changed to an empty string "" so the single tags would not be closed.
- `ssl_enabled::Bool`: default false. Server runs over SSL/HTTPS in development.
- `ssl_config::MbedTLS.SSLConfig`: default `nothing`. If not `nothing` and `ssl_enabled`, it will use the config to start the server over HTTPS.
"""
mutable struct Settings
  server_port::Int
  server_host::String
  server_document_root::String

  server_handle_static_files::Bool
  server_signature::String

  app_env::String

  cors_headers::Dict{String,String}
  cors_allowed_origins::Vector{String}

  cache_duration::Int
  cache_storage::Symbol

  log_level::Logging.LogLevel
  log_to_file::Bool

  assets_fingerprinted::Bool

  session_key_name::String
  session_storage::Symbol

  inflector_irregulars::Vector{Tuple{String,String}}

  run_as_server::Bool

  websockets_server::Bool
  websockets_port::Int

  initializers_folder::String

  path_config::String
  path_env::String
  path_app::String
  path_resources::String
  path_lib::String
  path_helpers::String
  path_log::String
  path_tasks::String
  path_build::String
  path_plugins::String
  path_cache::String
  path_initializers::String
  path_db::String
  path_bin::String
  path_src::String

  webchannels_default_route::String
  webchannels_js_file::String
  webchannels_subscribe_channel::String
  webchannels_unsubscribe_channel::String
  webchannels_autosubscribe::Bool

  html_parser_close_tag::String
  html_parser_char_at::String
  html_parser_char_dot::String
  html_parser_char_column::String
  html_parser_char_dash::String

  ssl_enabled::Bool
  ssl_config::Union{MbedTLS.SSLConfig,Nothing}

  Settings(;
            server_port                 = (haskey(ENV, "PORT") ? parse(Int, ENV["PORT"]) : 8000), # default port for binding the web server
            server_host                 = ENV["HOST"],
            server_document_root        = "public",
            server_handle_static_files  = true,
            server_signature            = "Genie/$GENIE_VERSION/Julia/$VERSION",

            app_env                     = ENV["GENIE_ENV"],

            cors_headers  = Dict{String,String}(
              "Access-Control-Allow-Origin"       => "", # ex: "*" or "http://mozilla.org"
              "Access-Control-Expose-Headers"     => "", # ex: "X-My-Custom-Header, X-Another-Custom-Header"
              "Access-Control-Max-Age"            => "86400", # 24 hours
              "Access-Control-Allow-Credentials"  => "", # "true" or "false"
              "Access-Control-Allow-Methods"      => "", # ex: "POST, GET"
              "Access-Control-Allow-Headers"      => "", # ex: "X-PINGOTHER, Content-Type"
            ),
            cors_allowed_origins = String[],

            cache_duration    = 0,
            cache_storage     = :File,

            log_level     = Logging.Debug,
            log_to_file   = false,

            assets_fingerprinted  = false,

            session_key_name    = "__geniesid",
            session_storage     = :File,

            inflector_irregulars = Tuple{String,String}[],

            run_as_server = false,

            websockets_server = false,
            websockets_port   = server_port,

            initializers_folder = "initializers",

            path_config         = "config",
            path_env            = joinpath(path_config, "env"),
            path_app            = "app",
            path_resources      = joinpath(path_app, "resources"),
            path_lib            = "lib",
            path_helpers        = joinpath(path_app, "helpers"),
            path_log            = "log",
            path_tasks          = "tasks",
            path_build          = buildpath(),
            path_plugins        = "plugins",
            path_cache          = "cache",
            path_initializers   = joinpath(path_config, initializers_folder),
            path_db             = "db",
            path_bin            = "bin",
            path_src            = "src",

            webchannels_default_route       = "__",
            webchannels_js_file             = "channels.js",
            webchannels_subscribe_channel   = "subscribe",
            webchannels_unsubscribe_channel = "unsubscribe",
            webchannels_autosubscribe       = true,

            html_parser_close_tag = " /",
            html_parser_char_at = "!!",
            html_parser_char_dot = "!",
            html_parser_char_column = "!",
            html_parser_char_dash = "__",

            ssl_enabled = false,
            ssl_config = nothing
        ) =
              new(
                  server_port, server_host,
                  server_document_root, server_handle_static_files, server_signature,
                  app_env,
                  cors_headers, cors_allowed_origins,
                  cache_duration, cache_storage,
                  log_level, log_to_file,
                  assets_fingerprinted,
                  session_key_name, session_storage,
                  inflector_irregulars,
                  run_as_server,
                  websockets_server, websockets_port,
                  initializers_folder,
                  path_config, path_env, path_app, path_resources, path_lib, path_helpers, path_log, path_tasks, path_build,
                  path_plugins, path_cache, path_initializers, path_db, path_bin, path_src,
                  webchannels_default_route, webchannels_js_file, webchannels_subscribe_channel, webchannels_unsubscribe_channel, webchannels_autosubscribe,
                  html_parser_close_tag, html_parser_char_at, html_parser_char_dot, html_parser_char_column, html_parser_char_dash,
                  ssl_enabled, ssl_config
                )
end

end