text:
using ProgressMeter
const update! = ProgressMeter.update!

"""
Canonicalize a GitHub repository URL
"""
function canonicalize_source_url(url)
    repo_regex = r"(https:\/\/)?github.com\/([^\/]+)\/([^\/]+)\/?$"

    m = match(repo_regex, url)
    if m !== nothing
        _, user, repo = m.captures
        if !endswith(repo, ".git")
            return "https://github.com/$user/$repo.git"
        end
    end
    url
end

"""
Canonicalize URL to a file within a GitHub repo
"""
function canonicalize_file_url(url)
    blob_regex = r"(https:\/\/)?github.com\/([^\/]+)\/([^\/]+)\/blob\/([^\/]+)\/(.+)"
    m = match(blob_regex, url)
    if m !== nothing
        _, user, repo, ref, filepath = m.captures
        if length(ref) != 40 || !all(c->isnumeric(c) || c in 'a':'f' || c in 'A':'F', ref)
            # Ask github to resolve this ref for us
            ref = GitHub.reference(GitHub.Repo("$user/$repo"), "heads/$ref").object["sha"]
        end
        return "https://raw.githubusercontent.com/$user/$repo/$ref/$(filepath)"
    end
    url
end

struct GitTransferProgress
    total_objects::Cuint
    indexed_objects::Cuint
    received_objects::Cuint
    local_objects::Cuint
    total_deltas::Cuint
    indexed_deltas::Cuint
    received_bytes::Csize_t
end

function transfer_progress(progress::Ptr{GitTransferProgress}, p::Any)
    progress = unsafe_load(progress)
    p.n = progress.total_objects
    if progress.total_deltas != 0
        p.desc = "Resolving Deltas: "
        p.n = progress.total_deltas
        update!(p, Int(max(1, progress.indexed_deltas)))
    else
        update!(p, Int(max(1, progress.received_objects)))
    end
    return Cint(0)
end

"""
    clone(url::String, source_path::String)

Clone a git repository hosted at `url` into `source_path`, with a progress bar
displayed to stdout.
"""
function clone(url::String, source_path::String)
    # Clone with a progress bar
    p = Progress(0, 1, "Cloning: ")
    GC.@preserve p begin
        callbacks = LibGit2.RemoteCallbacks(
            transfer_progress=@cfunction(
                transfer_progress,
                Cint,
                (Ptr{GitTransferProgress}, Any)
            ),
            payload = p
        )
        fetch_opts = LibGit2.FetchOptions(callbacks=callbacks)
        clone_opts = LibGit2.CloneOptions(fetch_opts=fetch_opts, bare = Cint(true))
        return LibGit2.clone(url, source_path, clone_opts)
    end
end

"""
    download_source(state::WizardState)

Ask the user where the source code is coming from, then download and record the
relevant parameters, returning the source `url`, the local `path` it is stored
at after download, and a `hash` identifying the version of the code. In the
case of a `git` source URL, the `hash` will be a git treeish identifying the
exact commit used to build the code, in the case of a tarball, it is the
`sha256` hash of the tarball itself.
"""
function download_source(state::WizardState)
    entered_url = nothing
    while entered_url === nothing #&& !eof(state.ins)
        # First, ask the user where this is all coming from
        msg = replace(strip("""
        Please enter a URL (git repository or compressed archive) containing the
        source code to build:
        """), "\n" => " ")
        new_entered_url = nonempty_line_prompt("URL", msg; ins=state.ins, outs=state.outs)

        # Early-exit for invalid URLs, using HTTP.URIs.parse_uri() to ensure
        # it is a valid URL
        try
            HTTP.URIs.parse_uri(new_entered_url; strict=true)
            entered_url = new_entered_url
        catch e
            printstyled(state.outs, e.msg, color=:red, bold=true)
            println(state.outs)
            println(state.outs)
            continue
        end
    end

    # Did the user exit out with ^D or something else go horribly wrong?
    if entered_url === nothing
        error("Could not obtain source URL")
    end

    url = string(canonicalize_source_url(entered_url))
    if url != entered_url
        print(state.outs, "The entered URL has been canonicalized to\n")
        printstyled(state.outs, url, bold=true)
        println(state.outs)
        println(state.outs)
    end

    # Record the source path and the source hash
    source_path = joinpath(state.workspace, basename(url))
    local source_hash

    if endswith(url, ".git")
        # Clone the URL, record the current gitsha for the given branch
        repo = clone(url, source_path)

        msg = "You have selected a git repository. Please enter a branch, commit or tag to use.\n" *
        "Please note that for reproducability, the exact commit will be recorded, \n" *
        "so updates to the remote resource will not be used automatically; \n" *
        "you will have to manually update the recorded commit."
        #print(state.outs, msg, "\n> ")
        treeish = nonempty_line_prompt("git reference", msg; ins=state.ins, outs=state.outs)

        obj = try
            LibGit2.GitObject(repo, treeish)
        catch
            LibGit2.GitObject(repo, "origin/$treeish")
        end
        source_hash = LibGit2.string(LibGit2.GitHash(obj))

        # Tell the user what we recorded the current commit as
        print(state.outs, "Recorded as ")
        printstyled(state.outs, source_hash, bold=true)
        println(state.outs)
        close(repo)
    else
        # Download the source tarball
        source_path = joinpath(state.workspace, basename(url))

        if isfile(source_path)
            name, ext = splitext(basename(source_path))
            n = 1
            while isfile(joinpath(state.workspace, "$(name)_$n$ext"))
                n += 1
            end
            source_path = joinpath(state.workspace, "$(name)_$n$ext")
        end

        download_cmd = Pkg.PlatformEngines.gen_download_cmd(url, source_path)
        oc = OutputCollector(download_cmd; verbose=true, tee_stream=state.outs)
        try
            if !wait(oc)
                error()
            end
        catch
            error("Could not download $(url) to $(state.workspace)")
        end

        # Save the source hash
        open(source_path) do file
            source_hash = bytes2hex(sha256(file))
        end
    end

    # Spit back the url, local path and source hash
    return url, source_path, source_hash
end

"""
    step1(state::WizardState)

It all starts with a single step, the unabashed ambition to leave your current
stability and engage with the universe on a quest to create something new,
beautiful and unforseen.  It all ends with compiler errors.

This step selects the relevant platform(s) for the built binaries.
"""
function step1(state::WizardState)
    print_wizard_logo(state.outs)

    # Select a platform
    msg = "\t\t\t# Step 1: Select your platforms\n\n"
    printstyled(state.outs, msg, bold=true)
    terminal = TTYTerminal("xterm", state.ins, state.outs, state.outs)
    platform_select = request(terminal,
        "Make a platform selection",
        RadioMenu([
            "All Supported Platforms",
            "Select by Operating System",
            "Fully Custom Platform Choice",
        ])
    )
    println(state.outs)

    # Set `state.platforms` accordingly
    result = nothing
    if platform_select == 1
        state.platforms = supported_platforms()
    elseif platform_select == 2
        oses = sort(unique(map(typeof, supported_platforms())), by = repr)
        while true
            result = request(terminal,
                "Select operating systems",
                MultiSelectMenu(map(repr, oses))
            )
            result = map(x->oses[x], collect(result))
            if isempty(result)
                println("Must select at least one operating system")
            else
                break
            end
        end
        state.platforms = collect(filter(x->typeof(x) in result, supported_platforms()))
    elseif platform_select == 3
        platfs = supported_platforms()
        while true
            result = request(terminal,
                "Select platforms",
                MultiSelectMenu(map(repr, platfs))
            )
            if isempty(result)
                println("Must select at least one platform")
            else
                break
            end
        end
        state.platforms = collect(map(x->platfs[x], collect(result)))
    else
        error("Somehow platform_select was not a valid choice!")
    end

    println(state.outs)
end

function obtain_binary_deps(state::WizardState)
    msg = "\t\t\t# Step 2b: Obtain binary dependencies (if any)\n\n"
    printstyled(state.outs, msg, bold=true)

    q = "Do you require any (binary) dependencies? "
    # Pkg.PackageSpec return different types in different Julia versions so...
    state.dependencies = typeof(Pkg.PackageSpec(name="dummy"))[]
    if yn_prompt(state, q, :n) == :y
        terminal = TTYTerminal("xterm", state.ins, state.outs, state.outs)
        local resolved_deps
        jll_names = []
        while true
            jll_name = nonempty_line_prompt("package name", "Enter JLL package name:"; ins=state.ins, outs=state.outs)
            if !endswith(jll_name, "_jll")
                jll_name *= "_jll"
            end

            # Check to see if this JLL package name can be resolved:
            push!(jll_names, jll_name)
            all_resolved, resolved_deps = resolve_jlls(jll_names, outs=state.outs)

            if !all_resolved
                pop!(jll_names)
                if yn_prompt(state, "Unable to resolve \"$(jll_name)\"; enter a new one?", :y) != :y
                    break
                else
                    continue
                end
            end

            q = "Would you like to provide additional dependencies? "
            if yn_prompt(state, q, :n) != :y
                break
            end
        end
        # jll_names contains the valid names, resolved_deps potentially contains
        # unresolved deps so we filter them out here.
        state.dependencies = filter(x -> x.name ∈ jll_names, resolved_deps)
    end
    println(state.outs)
end

function obtain_source(state::WizardState)
    msg = "\t\t\t# Step 2a: Obtain the source code\n\n"
    printstyled(state.outs, msg, bold=true)

    # Create the workspace that we'll stash everything within
    state.workspace = tempname()
    mkpath(state.workspace)

    # These are the metadata we need to know about all the sources we'll be
    # building over the course of this journey we're on together.
    state.source_urls = String[]
    state.source_files = String[]
    state.source_hashes = String[]

    while true
        url, file, hash = download_source(state)
        push!(state.source_urls, url)
        push!(state.source_files, file)
        push!(state.source_hashes, hash)
        println(state.outs)

        q = "Would you like to download additional sources? "
        if yn_prompt(state, q, :n) != :y
            break
        end
    end
    println(state.outs)
end

function get_name_and_version(state::WizardState)
    yggdrasil_path = get_yggdrasil()

    while state.name === nothing
        msg = "Enter a name for this project.  This will be used for filenames:"
        new_name = nonempty_line_prompt("Name", msg; ins=state.ins, outs=state.outs)

        # Check to see if this project name already exists
        if case_insensitive_file_exists(joinpath(yggdrasil_path, yggdrasil_build_tarballs_path(new_name)))
            println(state.outs, "A build recipe with that name already exists within Yggdrasil.")

            if yn_prompt(state, "Choose a new project name?", :y) == :n
                break
            end
        else
            state.name = new_name
        end
    end

    msg = "Enter a version number for this project:"
    while state.version === nothing
        try
            state.version = VersionNumber(nonempty_line_prompt("Version", msg; ins=state.ins, outs=state.outs))
        catch e
            if isa(e, ArgumentError)
                println(state.outs, e.msg)
                continue
            end
            rethrow(e)
        end
    end
end

@enum Compilers C=1 Go Rust
function get_compilers(state::WizardState)
    while state.compilers === nothing
        compiler_descriptions = Dict(C => "C/C++/Fortran", Go => "Go", Rust => "Rust")
        compiler_symbols = Dict(Int(C) => :c, Int(Go) => :go, Int(Rust) => :rust)
        terminal = TTYTerminal("xterm", state.ins, state.outs, state.outs)
        result = nothing
        while true
            select_menu = MultiSelectMenu([compiler_descriptions[i] for i in instances(Compilers)])
            select_menu.selected = Set([Int(C)])
            result = request(terminal,
                             "Select compilers for the project",
                             select_menu
                             )
            if isempty(result)
                println("Must select at least one platform")
            else
                break
            end
        end
        state.compilers = map(c -> compiler_symbols[c], collect(result))
    end
end

function get_preferred_version(state::WizardState, compiler::AbstractString,
                               available_versions=Vector{Integer})
    terminal = TTYTerminal("xterm", state.ins, state.outs, state.outs)
    version_selected = request(terminal, "Select the preferred $(compiler) version",
                               RadioMenu(string.(available_versions)))
    if compiler == "GCC"
        state.preferred_gcc_version = available_versions[version_selected]
    elseif compiler == "LLVM"
        state.preferred_llvm_version = available_versions[version_selected]
    end
end

"""
    step2(state::WizardState)

This step obtains the source code to be built and required binary dependencies.
"""
function step2(state::WizardState)
    obtain_source(state)
    obtain_binary_deps(state)
    get_name_and_version(state)
    if yn_prompt(state, "Do you want to customize the set of compilers?", :n) == :y
        get_compilers(state)
        get_preferred_version(state, "GCC", getversion.(available_gcc_builds))
        get_preferred_version(state, "LLVM", getversion.(available_llvm_builds))
    else
        state.compilers = [:c]
        # Default GCC version is the oldest one
        state.preferred_gcc_version = getversion(available_gcc_builds[1])
        # Default LLVM version is the latest one
        state.preferred_llvm_version = getversion(available_llvm_builds[end])
    end
end


text:
module RulesProbitInTest

using Test
using ReactiveMP
using Random
using Distributions
using StatsFuns: normcdf, normccdf, normlogcdf, normlogccdf, normlogpdf, normpdf

import ReactiveMP: @test_rules

@testset "rules:Probit:in" begin

    @testset "Belief Propagation: (m_out::PointMass, )" begin

        @test_rules [ with_float_conversions = true ] Probit(:in, Marginalisation) [
            (input = (m_out = PointMass(1.0), ), output = ContinuousUnivariateLogPdf((z) -> log(normcdf(z)))),
            (input = (m_out = PointMass(0.8), ), output = ContinuousUnivariateLogPdf((z) -> log(0.2 + 0.6*normcdf(z)))),
            (input = (m_out = PointMass(0.5), ), output = ContinuousUnivariateLogPdf((z) -> log(0.5))),
            (input = (m_out = PointMass(0.0), ), output = ContinuousUnivariateLogPdf((z) -> log(1 - normcdf(z)))),
        ]

    end

    @testset "Expectation Propagation: (m_out::PointMass, m_in::UnivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] Probit(:in, MomentMatching) [
            (input = (m_out = PointMass(1.0), m_in = NormalMeanVariance(1.0, 0.5) ), output = NormalWeightedMeanPrecision(0.6723616582693972, 0.32950039939606945) ),
            (input = (m_out = PointMass(1.0), m_in = NormalMeanPrecision(1.0, 2.0) ), output = NormalWeightedMeanPrecision(0.6723616582693972, 0.32950039939606945) ),
            (input = (m_out = PointMass(1.0), m_in = NormalWeightedMeanPrecision(2.0, 2.0) ), output = NormalWeightedMeanPrecision(0.6723616582693972, 0.32950039939606945) ),
            (input = (m_out = PointMass(0.0), m_in = NormalMeanVariance(1.0, 0.5) ), output = NormalWeightedMeanPrecision(-0.821224653874111, 0.7003447377360019) ),
            (input = (m_out = PointMass(0.0), m_in = NormalMeanPrecision(1.0, 2.0) ), output = NormalWeightedMeanPrecision(-0.821224653874111, 0.7003447377360019) ),
            (input = (m_out = PointMass(0.0), m_in = NormalWeightedMeanPrecision(2.0, 2.0) ), output = NormalWeightedMeanPrecision(-0.821224653874111, 0.7003447377360019) )
        ]

    end

    @testset "Expectation Propagation: (m_out::Bernoulli, m_in::UnivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] Probit(:in, MomentMatching) [
            (input = (m_out = Bernoulli(1.0), m_in = NormalMeanVariance(1.0, 0.5) ), output = NormalWeightedMeanPrecision(0.6723616582693972, 0.32950039939606945) ),
            (input = (m_out = Bernoulli(1.0), m_in = NormalMeanVariance(1.0, 0.5) ), output = NormalWeightedMeanPrecision(0.6723616582693972, 0.32950039939606945) ),
            (input = (m_out = Bernoulli(1.0), m_in = NormalMeanVariance(1.0, 0.5) ), output = NormalWeightedMeanPrecision(0.6723616582693972, 0.32950039939606945) ),
            (input = (m_out = Bernoulli(0.8), m_in = NormalMeanVariance(1.0, 0.5) ), output = NormalWeightedMeanPrecision(0.427017495944859, 0.199141999223396) ),
            (input = (m_out = Bernoulli(0.8), m_in = NormalMeanPrecision(1.0, 2.0) ), output = NormalWeightedMeanPrecision(0.427017495944859, 0.199141999223396) ),
            (input = (m_out = Bernoulli(0.8), m_in = NormalWeightedMeanPrecision(2.0, 2.0) ), output = NormalWeightedMeanPrecision(0.427017495944859, 0.199141999223396) ),
            (input = (m_out = Bernoulli(0.0), m_in = NormalMeanVariance(1.0, 0.5) ), output = NormalWeightedMeanPrecision(-0.821224653874111, 0.7003447377360019) ),
            (input = (m_out = Bernoulli(0.0), m_in = NormalMeanPrecision(1.0, 2.0) ), output = NormalWeightedMeanPrecision(-0.821224653874111, 0.7003447377360019) ),
            (input = (m_out = Bernoulli(0.0), m_in = NormalWeightedMeanPrecision(2.0, 2.0) ), output = NormalWeightedMeanPrecision(-0.821224653874111, 0.7003447377360019) )
        ]

        # Test against an extreme case with m_out = Bernoulli(0.5)

        @test_rules [ with_float_conversions = false, atol = 1e-13 ] Probit(:in, MomentMatching) [
            (input = (m_out = Bernoulli(0.5), m_in = NormalMeanVariance(1.0, 0.5) ), output = NormalWeightedMeanPrecision(0.0, 1.0 * tiny) ),
            (input = (m_out = Bernoulli(0.5), m_in = NormalMeanPrecision(1.0, 2.0) ), output = NormalWeightedMeanPrecision(0.0, 1.0 * tiny) ),
            (input = (m_out = Bernoulli(0.5), m_in = NormalWeightedMeanPrecision(2.0, 2.0) ), output = NormalWeightedMeanPrecision(0.0, 1.0 * tiny) ),
        ]

        @test_rules [ with_float_conversions = false, atol = 1e-7 ] Probit(:in, MomentMatching) [
            (input = (m_out = Bernoulli(0.5f0), m_in = NormalMeanVariance(1.0f0, 0.5f0) ), output = NormalWeightedMeanPrecision(0f0, 1f0 * tiny) ),
            (input = (m_out = Bernoulli(0.5f0), m_in = NormalMeanPrecision(1.0f0, 2.0f0) ), output = NormalWeightedMeanPrecision(0f0, 1f0 * tiny) ),
            (input = (m_out = Bernoulli(0.5f0), m_in = NormalWeightedMeanPrecision(2.0f0, 2.0f0) ), output = NormalWeightedMeanPrecision(0f0, 1f0 * tiny) ),
        ]

        @test_rules [ with_float_conversions = false, atol = 1e-25 ] Probit(:in, MomentMatching) [
            (input = (m_out = Bernoulli(big"0.5"), m_in = NormalMeanVariance(big"1.0", big"0.5") ), output = NormalWeightedMeanPrecision(big"0.0", big"1.0" * tiny) ),
            (input = (m_out = Bernoulli(big"0.5"), m_in = NormalMeanPrecision(big"1.0", big"2.0") ), output = NormalWeightedMeanPrecision(big"0.0", big"1.0" * tiny) ),
            (input = (m_out = Bernoulli(big"0.5"), m_in = NormalWeightedMeanPrecision(big"2.0", big"2.0") ), output = NormalWeightedMeanPrecision(big"0.0", big"1.0" * tiny) ),
        ]

    end

end

end

ext: jl


text:
module Cambrian

using Random
using Logging
using Statistics
using YAML
using JSON
import Formatting
import Dates

include("config.jl")
include("logger.jl")
include("individual.jl")
include("geography.jl")
include("evolution.jl")
include("selection.jl")
include("evaluation.jl")
include("step.jl")
include("oneplus.jl")
include("GA.jl")

##
# Giving these methods bang names, since they are expected to mutate the
# AbstractEvolution instance.
##

@deprecate populate(e) Cambrian.populate!(e) true
@deprecate evaluate(e) Cambrian.evaluate!(e) true
@deprecate generation(e) Cambrian.generation!(e) true

end


text:
tests = ["gauss"]

for t in tests
    fpath = "$t.jl"
    println("running $fpath ...")
    include(fpath)
end


text:
function next(side::HorizonSide)
    return HorizonSide(mod(Int(side) + 1, 4))
end

function f()

    # Тут решение.
    side = Nord
    a = 1
    while !ismarker(robot)
        for i = 1:a
            if !ismarker(robot)
                move!(robot, side)
            end
        end
        side = next(side)
        for i = 1:a
            if !ismarker(robot)
                move!(robot, side)
            end
        end
        side = next(side)
        a += 1
    end

end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
using TimeSeriesClustering
using Test
using JLD2
using Random
using Cbc
using StatsBase

include("test_utils.jl")

include("run_clust.jl")
include("extreme_vals.jl")
include("datastructs.jl")
include("load_data.jl")
include("utils.jl")


text:
"""
    XFloats.jl

See: [`XFloat16`](@ref), [`XFloat32`](@ref)
"""
module XFloats

export XFloat16, XFloat32

import LinearAlgebra: (*), (/), (\), dot, cross, factorize, inv, det, logdet, logabsdet,
                      size, eigvals, eigvecs, svdvals
using LinearAlgebra
using Random

# import SpecialFunctions
# using SpecialFunctions

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

"""
    XFloat16 is a floating point type

A more accurate alternative to Float16 that holds a Float32 and that performs as a Float32.
"""
primitive type XFloat16 <: AbstractFloat 32 end

"""
    XFloat32 is a floating point type

A more accurate alternative to Float32 that holds a Float64 and that performs as a Float64.
"""
primitive type XFloat32 <: AbstractFloat 64 end

XFloat16(x::XFloat16) = x
XFloat16(x::XFloat32) = reinterpret(XFloat16, Float32(reinterpret(Float64, x)))
XFloat32(x::XFloat16) = reinterpret(XFloat32, Float64(reinterpret(Float32, x)))
XFloat32(x::XFloat32) = x

XFloat16(x::Float32) = reinterpret(XFloat16, x)
Float32(x::XFloat16) = reinterpret(Float16, x)
XFloat32(x::Float64) = reinterpret(XFloat32, x)
Float64(x::XFloat32) = reinterpret(Float64, x)

include("type/construct.jl")
include("type/promote_convert.jl")
include("type/representations.jl")

include("imports.jl")
include("dispatch/specialize.jl")
include("dispatch/linearalgebra.jl")

end # module XFloats


text:
# Copyright 2017-2019 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

export LabelEncoder, transform, inverse_transform

import Base: broadcastable

struct LabelEncoder{LabelType}
    imap::Dict{LabelType,Int}
    labels::Vector{LabelType}
    freqs::Vector{Int}

end

function fit(LabelEncoder, y::AbstractVector{LabelType}) where LabelType
    L = Dict{LabelType,Int}()
    for c in y
        L[c] = get(L, c, 0) + 1
    end

    labels = collect(keys(L))
    sort!(labels)
    freqs = [L[c] for c in labels]
    imap = Dict(c => i for (i, c) in enumerate(labels))
    LabelEncoder{LabelType}(imap, labels, freqs)
end

function transform(le::LabelEncoder{LabelType}, y::LabelType)::Int where LabelType
    le.imap[y]
end

function inverse_transform(le::LabelEncoder{LabelType}, y::Int)::LabelType where LabelType
    le.labels[y]
end

function broadcastable(le::LabelEncoder)
    (le,)
end


text:
using Revise, ProgressMeter, OffsetArrays, Kinetic

cd(@__DIR__)
D = read_dict("shear.txt")
for key in keys(D)
    s = Symbol(key)
    @eval $s = $(D[key])
end

begin
    γ = heat_capacity_ratio(inK, 2)
    set = Setup(matter, case, space, flux, collision, nSpecies, interpOrder, limiter, boundary, cfl, maxTime)
    pSpace = PSpace1D(x0, x1, nx, nxg)
    vSpace = VSpace2D(umin, umax, nu, vmin, vmax, nv, vMeshType, nug, nvg)

    μᵣ = ref_vhs_vis(knudsen, alphaRef, omegaRef)
    gas = Gas(
        knudsen,
        mach,
        prandtl,
        inK,
        γ,
        omega,
        alphaRef,
        omegaRef,
        μᵣ,
    )

    primL = [1.0, 0.0, 1.0, 1.0]
    primR = [1.0, 0.0, -1.0, 2.0]

    wL = prim_conserve(primL, γ)
    wR = prim_conserve(primR, γ)

    HL = maxwellian(vSpace.u, vSpace.v, primL)
    HR = maxwellian(vSpace.u, vSpace.v, primR)

    BL = HL .* inK ./ (2.0 * primL[end])
    BR = HR .* inK ./ (2.0 * primR[end])

    bc = zeros(4)

    ib = IB2F(
            wL,
            primL,
            HL,
            BL,
            bc,
            wR,
            primR,
            HR,
            BR,
            bc,
        )

    outputFolder = pwd()

    ks = SolverSet(set, pSpace, vSpace, gas, ib, outputFolder)
    KS = ks
end

begin
    ctr = OffsetArray{ControlVolume1D2F}(undef, axes(KS.pSpace.x, 1))
    face = Array{Interface1D2F}(undef, KS.pSpace.nx + 1)

    idx0 = (eachindex(pSpace.x) |> collect)[1]
    idx1 = (eachindex(pSpace.x) |> collect)[end]

    for i in eachindex(ctr)
        if i <= KS.pSpace.nx ÷ 2
            ctr[i] = ControlVolume1D2F( KS.pSpace.x[i], KS.pSpace.dx[i], ks.ib.wL, ks.ib.primL,
            ks.ib.hL, ks.ib.bL )
        else
            ctr[i] = ControlVolume1D2F( KS.pSpace.x[i], KS.pSpace.dx[i], ks.ib.wR, ks.ib.primR,
            ks.ib.hR, ks.ib.bR )
        end
    end

    face = Array{Interface1D2F}(undef, KS.pSpace.nx+1)
    for i=1:KS.pSpace.nx+1
        face[i] = Interface1D2F(ks.ib.wL, ks.ib.hL)
    end
end

begin
    iter = 0
    res = zeros(4)
    simTime = 0.0
    dt = Kinetic.timestep(KS, ctr, simTime)
    maxTime = vhs_collision_time(ks.ib.primL, μᵣ, omega)
    nt = Int(floor(maxTime / dt))
end

# There're no default solver for 1D simulation with 2D setting
# Let's do it manually
@showprogress for iter in 1:nt
    #Kinetic.reconstruct!(KS, ctr)

    @inbounds Threads.@threads for i in eachindex(face)
        flux_kfvs!(
            face[i].fw,
            face[i].fh,
            face[i].fb,
            ctr[i-1].h,
            ctr[i-1].b,
            ctr[i].h,
            ctr[i].b,
            ks.vSpace.u,
            ks.vSpace.v,
            ks.vSpace.weights,
            dt,
            1.0
        )
    end

    @inbounds Threads.@threads for i in 1:ks.pSpace.nx
        #--- store W^n and calculate shakhov term ---#
        w_old = deepcopy(ctr[i].w)

        #--- update W^{n+1} ---#
        @. ctr[i].w += (face[i].fw - face[i+1].fw) / ctr[i].dx
        ctr[i].prim .= conserve_prim(ctr[i].w, ks.gas.γ)

        #--- calculate M^{n+1} and tau^{n+1} ---#
        MH = maxwellian(ks.vSpace.u, ks.vSpace.v, ctr[i].prim)
        MB = MH .* ks.gas.K ./ (2.0 * ctr[i].prim[end])
        τ = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)

        #--- update distribution function ---#
        for q in axes(MH,2), p in axes(MH,1)
            ctr[i].h[p,q] = (ctr[i].h[p,q] + (face[i].fh[p,q] - face[i+1].fh[p,q]) / ctr[i].dx + dt / τ * MH[p,q]) / (1.0 + dt / τ)
            ctr[i].b[p,q] = (ctr[i].b[p,q] + (face[i].fb[p,q] - face[i+1].fb[p,q]) / ctr[i].dx + dt / τ * MB[p,q]) / (1.0 + dt / τ)
        end
    end
end

sol = zeros(ks.pSpace.nx, 10)
for i in 1:ks.pSpace.nx
    sol[i, 1:3] = ctr[i].prim[1:3]
    sol[i, 4] = 1. / ctr[i].prim[4]
end

using Plots
plot(ks.pSpace.x[1:ks.pSpace.nx], sol[:,1])
plot(ks.pSpace.x[1:ks.pSpace.nx], sol[:,2])
plot(ks.pSpace.x[1:ks.pSpace.nx], sol[:,3])
plot(ks.pSpace.x[1:ks.pSpace.nx], sol[:,4])

u1d = VSpace1D(umin, umax, nu, vMeshType)
f = zeros(ks.vSpace.nv)
for j in axes(ctr[1].h, 2), i in axes(ctr[1].h, 1)
    f[j] = 0.5 * (sum(@. u1d.weights * ctr[ks.pSpace.nx÷2].h[:,j]) +
        sum(@. u1d.weights * ctr[ks.pSpace.nx÷2+1].h[:,j]))
end
plot(ks.vSpace.v[end÷2,:,1], f)


text:
function _population(pop::Population, ids)
  x = Float64[]
  y = Float64[]
  for i in ids
    push!(x, pop.risks[i, :x])
    push!(y, pop.risks[i, :y])
  end
  return x, y
end

function _ids_by_state(events::Events{T},
                       state::DiseaseState,
                       time::Float64) where T <: EpidemicModel
  if time < 0.0
    @error "Time must be ≥ 0.0"
  elseif state ∉ _state_progressions[T]
    @error "Invalid state specified"
  end
  ids = Int64[]
  if state == _state_progressions[T][1] # S
    nextstate = advance(state, T) # Either E or I
    append!(ids, findall(events[nextstate] .> Ref(time))) # E/I after `time`
    append!(ids, findall(isnan.(events[nextstate]))) # Never E/I
  elseif state in _state_progressions[T][2:end-1]
    nextstate = advance(state, T) # Either I or R
    append!(ids, findall((events[state] .<= Ref(time)) .& (events[nextstate] .> Ref(time)))) # E/I at or before time and I/R after time
    append!(ids, findall((events[state] .<= Ref(time)) .& isnan.(events[nextstate]))) # E/I at or before time and never I/R
  elseif state == _state_progressions[T][end] # I or R
    append!(ids, findall(events[state] .<= Ref(time))) # I/R at or before time
  end
  @debug "Individual(s) in state $state at t = $time" Individuals = ids
  return ids
end

@recipe function f(pop::Population) where T <: EpidemicModel
  xguide --> ""
  yguide --> ""
  legend --> :none
  xlims --> extrema(pop.risks[!, :x]) .+ (sum(extrema(pop.risks[!, :x]).*(-1,1)) .* (-0.05, 0.05))
  ylims --> extrema(pop.risks[!, :y]) .+ (sum(extrema(pop.risks[!, :y]).*(-1,1)) .* (-0.05, 0.05))
  aspect_ratio := :equal
  seriestype := :scatter
  markerstrokecolor --> :black
  markercolor --> :black
  markersize --> 2.75
  axis --> nothing
  titlefontcolor --> :black
  framestyle --> :none
  _population(pop, 1:pop.individuals)
end

@recipe function f(pop::Population,
                   events::Events{T},
                   time::Float64) where T <: EpidemicModel
  xguide --> ""
  yguide --> ""
  legend --> :topright
  xlims --> extrema(pop.risks[!, :x]) .+ (sum(extrema(pop.risks[!, :x]).*(-1,1)) .* (-0.05, 0.05))
  ylims --> extrema(pop.risks[!, :y]) .+ (sum(extrema(pop.risks[!, :y]).*(-1,1)) .* (-0.05, 0.05))
  aspect_ratio := :equal
  markerstrokecolor --> :black
  markersize --> 2.75
  axis --> nothing
  titlefontcolor --> :black
  framestyle --> :none
  @series begin
    ids_susceptible = _ids_by_state(events, State_S, time)
    x, y = _population(pop, ids_susceptible)
    seriestype := :scatter
    seriescolor --> :purple
    label --> "S"
    x, y
  end
  if T in [SEIR; SEI]
    @series begin
      ids_exposed = _ids_by_state(events, State_E, time)
      x, y = _population(pop, ids_exposed)
      seriestype := :scatter
      seriescolor --> :lightblue4
      label --> "E"
      x, y
    end
  end
  @series begin
    ids_infected = _ids_by_state(events, State_I, time)
    x, y = _population(pop, ids_infected)
    seriestype := :scatter
    seriescolor --> :lightgreen
    label --> "I"
    x, y
  end
  if T in [SEIR; SIR]
    @series begin
      ids_removed = _ids_by_state(events, State_R, time)
      x, y = _population(pop, ids_removed)
      seriestype := :scatter
      seriescolor --> :yellow
      label --> "R"
      x, y
    end
  end
end


text:
using Test
using MatrixLM

using DataFrames
using Random
using LinearAlgebra
using GLM

@testset "GLM compatible" begin

    # Tolerance for tests
    tol = 10.0^(-7)

    # Dimensions of matrices
    n = 100
    m = 200
    p = 10
    q = 20

    # Generate some matrices.
    Random.seed!(4)
    X = rand(n,p)
    Z = rand(m,q)
    B = rand(1:20,p,q)
    E = randn(n,m)
    Y = X*B*transpose(Z)+E

    # Data frame to be passed into lm
    GLMData = DataFrame(hcat(vec(Y), kron(Z,X)), :auto)
    # lm estimate
    GLMEst = lm(Matrix(GLMData[:,2:end]), Vector(GLMData[:,1]))

    # Put together RawData object for MLM
    MLMData = RawData(Response(Y), Predictors(X, Z))
    # mlm estimate
    # MLMEst = mlm(MLMData, hasXIntercept = false, hasZIntercept = false)
    MLMEst = mlm(MLMData, hasXIntercept = false, hasZIntercept = false)

    @test isapprox(GLM.coef(GLMEst), vec(MatrixLM.coef(MLMEst)), atol=tol)
    @test isapprox(GLM.predict(GLMEst), vec(MatrixLM.predict(MLMEst).Y), atol=tol)
    @test LinearAlgebra.issymmetric(round.(MLMEst.sigma, digits=10)) # and also positive semi-definite?

end

# Tests for variance shrinkage, WLS
# Contrasts
# Permutations


text:
using LAP_julia: inds_to_points, assess_flow_quality

## testing

flow_size = size(flow)

## method
img, imgw, flow = gen_init()
timer=TimerOutput("reg alg: sparse lap");
method_kwargs = Dict(:timer => timer, :point_count => 200, :spacing => 15)
flow_est, source_reg, timer, results, (estim_at_inds, inds) = test_registration_alg(sparse_lap, img, imgw, flow, [25, (51,51)], method_kwargs, timer=timer)

showflow(flow)
showflow(flow_est)
##

estim_flow = interpolate_flow_quad(estim_at_inds, inds, size(flow))
estim_flow_rbf = interpolate_flow(estim_at_inds, inds, size(flow))

assess_flow_quality(flow, estim_flow)

assess_flow_quality(flow, estim_flow_rbf)

showflow(flow, figtitle="original")
showflow(estim_flow, figtitle="quad estimate")
showflow(estim_flow_rbf, figtitle="rbf estimate")

estim_flow = interpolate_flow_quad(estim_at_inds, inds, size(flow))
estim_flow_rbf = interpolate_flow(estim_at_inds, inds, size(flow))

assess_flow_quality(flow, estim_flow)
assess_flow_quality(flow, estim_flow_rbf)

showflow(estim_flow)
showflow(flow)
showflow(flow_truth .- estim_flow)

showflow(rbf_flow)
assess_flow_quality(flow, rbf_flow)
showflow(flow_truth .- rbf_flow)

@benchmark interpolate_flow_quad(estim_at_inds, inds, size(flow))

## benchmark both

rand_inds = rand(CartesianIndices(size(flow)), 15)
lsdkf = interpolate_flow_quad(map(x -> flow[x], rand_inds), rand_inds, size(flow))
flow
showflow(lsdkf)
showflow(flow)

# function make_flow_at_ind(ind, a, b, c, d, e, f)
#     return a +
#            b * ind[1] +
#            c * ind[2] +
#            d * ind[1]^2 +
#            e * ind[2]^2 +
#            f * ind[1]*ind[2]
# end

# function make_flow(coeffs, flow_size)
#     return reduce((coeff, first) -> coeff .* ones(flow_size) .+ first, coeffs; init=0)
# end


text:
## This is a crude replication of Carvalho and Vilela (2015)
## By Daniel Coutinho and Rafaela

####### WARNING############
## NOT DEBUGGED AS MUCH AS IT SHOULD!
#################

### TODO

## 1. Pi and Psi matrices
## 2. More comments
## 3. Coefficients as strings for latexify

#######

### ORDER
## c_t
##c_{t-1}
##i_{t-1}
##i_{t-2}
## pi_t
## gamma_t
## y_t
## tot_{t-1}
## y^{star}_t
## q_{t-1}
## s_t
## s_{t-1}
## pi^{star}_t
## pi_{D,t}
## pi_{D,t-1}
## pi_{I,t}
## pi_{I,t-1}
## a_t
## z_t
## i^{star}_t
## varepsilon_{cp,t}
## phi_t

#Creating Matrices

G0 = Array{Any}(nothing,23,23)
G1 = Array{Any}(nothing,23,23)
Pi = zeros(23,6)
Psi = zeros(23,8)

#First equation

G0[1,1] = 1
G0[1,5] = "(1-h)/sigma()"
G0[1,6] = "(1-h)/sigma()"
G0[1,3] = "-(1-h)/sigma()"

G1[1,1] = "(1-h)"
G1[1,2] = "h"
G1[1,6] = "(1-h)/sigma()"

#Second Equation

G1[2,7] = -1
G1[2,1] = "1-alpha()"
G1[2,8] = "alpha()*eta()*(2-alpha())"
G1[2,9] = "alpha()*eta()"
G1[2,10] = "alpha()"

#Third Equation

G0[3,8] = 1

G1[3,8] = 1
G1[3,15] = -1
G1[3,17] = 1

#Fourth Equation

G1[4,11] = -1
G1[4,9] = 1
G1[4,8] = "1 -alpha()"

#Fifth equation

G0[5,11] = -1
G0[5,14] = -1

G1[5,11] = 1
G1[5,12] = 1
G1[5,13] = -1
G1[5,5] = -1

#Equation number 6

G0[6,15] = "beta()"

G1[6,15] = "1+delta()_d"
G1[6,16] = "beta()delta()_d"
G1[6,8] = "-kappa()_d*alpha()"
G1[6,7] = "-kappa()_d*phi()"
G1[6,1] = "-kappa()_d*sigma()/(1-h)"
G1[6,2] = "kappa()_d*sigma()/(1-h)"
G1[6,19] = "kappa()_d*(1+phi())"

#Equation number 7

G0[7,17] = "beta()"

G1[7,17] = "(1+beta()*delta()_i)"
G1[7,18] = "delta()_i"
G1[7,9] = "kappa()_i"

#Equation number 8

G0[8,8] = "-alpha()"

G1[8,5] = -1
G1[8,15] = 1
G1[8,8] = "alpha()"

#Equation number 9

G0[9,20] = 1
G0[9,8] = "alpha()"

G1[9,20] = "1/beta()"
G1[9,9] = "-alpha()"
G1[9,7] = 1
G1[9,1] = -1

# Equation number 10

G0[10,12] = 1
G0[10,21] = 1
G0[10,3] = -1

G1[10,20] = "chi()"

# Equation number 11: Taylor rule

G0[11,3] = 1
G1[11,3] = "rho()_it1"
G1[11,4] = "rho()_it2"
G1[11,5] = "(1-rho()_it1-rho()_it2)*lambda()_pi"
G1[11,7] = "(1-rho()_it1-rho()_it2)*lambda()_y"
G1[11,12] = "(1-rho()_it1-rho()_it2)*lambda()_s"
G1[11,13] = "-(1-rho()_it1-rho()_it2)*lambda()_s" ## its for the delta tot

###################
## Foreign Block ##
###################

# Equation 12

G0[12,10] = 1

G1[12,10] = "a_1yy"
G1[12,14] = "a_1y*pi"
G1[12,21] = "a_1yi"

#Equation 13

G0[13,14] = 1
G0[13,10] = "a_0piy"

G1[13,10] = "a_1piy"
G1[13,14] = "a_1pipi"
G1[13,21] = "a_1pii"

#Equation 14

G0[14,21] = 1
G0[14,14] = "a_0ipi"
G0[14,10] = "a_0iy"

G1[14,10] = "a_1iy"
G1[14,14] = "a_1ipi"
G1[14,21] = "a_1ii"

#################
## Shocks block #
#################

# Equation 15

G0[15,19] = 1

G1[15,19] = "rho()_a"

# Equation 16

G0[16,6] = 1

G1[16,6] = "rho()_gamma"

# Equation 17

G0[17,22] = 1

G1[17,22] = "rho()_varepsilon_cp"

# Equation 18

G0[18,23] = 1

G1[18,23] = "rho()_phi"

####################
## Identity block ##
####################

# Equation 19

G0[19,1] = 1
G1[19,2] = 1

# Eqquation 20

G0[20,3] = 1
G1[20,4] = 1

# Equation 21

G0[21,12] = 1
G1[21,13] = 1

# Equation 22

G0[22,15] = 1
G1[22,16] = 1

# Equation 23

G0[23,17] = 1
G1[23,18] = 1

G0[findall(G0 .== nothing)] .= 0
G1[findall(G1 .== nothing)] .= 0

latexify(G0)
latexify(G1)

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


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

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

"""

Creates or updates an authorization in the specified express route circuit.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: authorizationName::String (required)
Param: authorizationParameters::ExpressRouteCircuitAuthorization (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ExpressRouteCircuitAuthorization
"""
function _swaggerinternal_expressRouteCircuitAuthorizationsCreateOrUpdate(_api::ExpressRouteCircuitAuthorizationsApi, resourceGroupName::String, circuitName::String, authorizationName::String, authorizationParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", ExpressRouteCircuitAuthorization, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/authorizations/{authorizationName}", ["azure_auth"], authorizationParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "authorizationName", authorizationName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function expressRouteCircuitAuthorizationsCreateOrUpdate(_api::ExpressRouteCircuitAuthorizationsApi, resourceGroupName::String, circuitName::String, authorizationName::String, authorizationParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitAuthorizationsCreateOrUpdate(_api, resourceGroupName, circuitName, authorizationName, authorizationParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function expressRouteCircuitAuthorizationsCreateOrUpdate(_api::ExpressRouteCircuitAuthorizationsApi, response_stream::Channel, resourceGroupName::String, circuitName::String, authorizationName::String, authorizationParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitAuthorizationsCreateOrUpdate(_api, resourceGroupName, circuitName, authorizationName, authorizationParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified authorization from the specified express route circuit.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: authorizationName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_expressRouteCircuitAuthorizationsDelete(_api::ExpressRouteCircuitAuthorizationsApi, resourceGroupName::String, circuitName::String, authorizationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/authorizations/{authorizationName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "authorizationName", authorizationName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function expressRouteCircuitAuthorizationsDelete(_api::ExpressRouteCircuitAuthorizationsApi, resourceGroupName::String, circuitName::String, authorizationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitAuthorizationsDelete(_api, resourceGroupName, circuitName, authorizationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function expressRouteCircuitAuthorizationsDelete(_api::ExpressRouteCircuitAuthorizationsApi, response_stream::Channel, resourceGroupName::String, circuitName::String, authorizationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitAuthorizationsDelete(_api, resourceGroupName, circuitName, authorizationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified authorization from the specified express route circuit.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: authorizationName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ExpressRouteCircuitAuthorization
"""
function _swaggerinternal_expressRouteCircuitAuthorizationsGet(_api::ExpressRouteCircuitAuthorizationsApi, resourceGroupName::String, circuitName::String, authorizationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ExpressRouteCircuitAuthorization, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/authorizations/{authorizationName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "authorizationName", authorizationName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function expressRouteCircuitAuthorizationsGet(_api::ExpressRouteCircuitAuthorizationsApi, resourceGroupName::String, circuitName::String, authorizationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitAuthorizationsGet(_api, resourceGroupName, circuitName, authorizationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function expressRouteCircuitAuthorizationsGet(_api::ExpressRouteCircuitAuthorizationsApi, response_stream::Channel, resourceGroupName::String, circuitName::String, authorizationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitAuthorizationsGet(_api, resourceGroupName, circuitName, authorizationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all authorizations in an express route circuit.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: AuthorizationListResult
"""
function _swaggerinternal_expressRouteCircuitAuthorizationsList(_api::ExpressRouteCircuitAuthorizationsApi, resourceGroupName::String, circuitName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", AuthorizationListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/authorizations", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function expressRouteCircuitAuthorizationsList(_api::ExpressRouteCircuitAuthorizationsApi, resourceGroupName::String, circuitName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitAuthorizationsList(_api, resourceGroupName, circuitName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function expressRouteCircuitAuthorizationsList(_api::ExpressRouteCircuitAuthorizationsApi, response_stream::Channel, resourceGroupName::String, circuitName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitAuthorizationsList(_api, resourceGroupName, circuitName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export expressRouteCircuitAuthorizationsCreateOrUpdate, expressRouteCircuitAuthorizationsDelete, expressRouteCircuitAuthorizationsGet, expressRouteCircuitAuthorizationsList


text:
# deprecation of specialized version of constructors
const deprecated_constructors =
    Dict(:symbol => :Symbol)

const not_constructible = Set([:Union, :Tuple, :Type])

function lintfuncargtype(ex, ctx::LintContext)
    lintexpr(ex, ctx)
    if isexpr(ex, :curly)
        st = 2
        en = 1
        if ex.args[1] == :Array
            en = 2
        elseif ex.args[1] == :Dict
            en = 3
        end
        for i in st:en
            if in(ex.args[i], [:Number])
                msg(ctx, :E533, ex, "type parameters are invariant, try f{T<:Number}(x::T)...")
            end
        end
    end
end

function istype(ctx::LintContext, x)
    obj = abstract_eval(ctx, x)
    obj !== nothing && isa(obj, Type)
end

# if ctorType isn't symbol("") then we are in the context of
# a constructor for a type. We would check
# * if the function name matches the type name
function lintfunction(ex::Expr, ctx::LintContext; ctorType = Symbol(""), isstaged=false)
    if length(ex.args) == 1 && isa(ex.args[1], Symbol)
        # generic function without methods
        set!(ctx.current, ex.args[1], VarInfo(location(ctx), Function))
        return
    end

    if !isa(ex.args[1], Expr)
        msg(ctx, :E121, ex.args[1], "Lint does not understand the expression")
        return
    end

    fname = Symbol("")
    if ex.args[1].head == :tuple # anonymous
        # do nothing
    elseif isexpr(ex.args[1].args[1], :(.))
        fname = ex.args[1].args[1]
    elseif isa(ex.args[1].args[1], Symbol)
        fname = ex.args[1].args[1]
    elseif !isa(ex.args[1].args[1], Expr)
        msg(ctx, :E121, ex.args[1].args[1], "Lint does not understand the expression")
        return
    elseif ex.args[1].args[1].head == :curly
        fname = ex.args[1].args[1].args[1]
    end
    if isa(fname, Symbol)
        # TODO: warn if it's a using'd thing
        finfo = lookup(ctx.current, fname)
        if finfo == nothing
            set!(ctx.current, fname, VarInfo(location(ctx), Function))
        else
            # TODO: warn if it's something bad
        end
    end

    ctx.scope = string(fname)
    if fname != Symbol("") && !occursin(ctx.file, "deprecate")
        isDeprecated = functionIsDeprecated(ex.args[1])
        if isDeprecated != nothing && !pragmaexists("Ignore deprecated $fname", ctx.current)
            msg(ctx, :E211, ex.args[1], "$(isDeprecated.message); See: " *
                "deprecated.jl $(isDeprecated.line)")
        end
    end

    defer!(ctx.current, MethodInfo(location(ctx), ex, isstaged))
end

function lintfunctionbody(ctx::LintContext, mi::MethodInfo)
    ex = mi.body
    isstaged = mi.isstaged
    oldloc = location(ctx)
    location!(ctx, location(mi))
    temporaryTypes = Any[]
    if isexpr(ex, [:(=), :function]) && isexpr(ex.args[1], :call) &&
       isexpr(ex.args[1].args[1], :curly)
        for i in 2:length(ex.args[1].args[1].args)
            adt = ex.args[1].args[1].args[i]
            if isa(adt, Symbol)
                if istype(ctx, adt)
                    msg(ctx, :E534, adt, "introducing a new name for an implicit " *
                        "argument to the function, use {T<:$(adt)}")
                else
                    push!(temporaryTypes, adt)
                end
            elseif isexpr(adt, :(<:))
                temptype = adt.args[1]
                typeconstraint = adt.args[2]
                if istype(ctx, temptype)
                    msg(ctx, :E536, temptype, "use {T<:...} instead of a known type")
                end
                if istype(ctx, typeconstraint)
                    dt = parsetype(ctx, typeconstraint)
                    if isconcretetype(dt)
                        msg(ctx, :E513, adt, "leaf type as a type constraint makes no sense")
                    end
                end
                push!(temporaryTypes, adt.args[1])
            end
        end
    end

    withcontext(ctx, LocalContext(ctx.current)) do
        # temporaryTypes are the type parameters in curly brackets, make them legal
        # in the current scope
        for t in temporaryTypes
            # TODO: infer when t is a type
            localset!(ctx.current, t, VarInfo(location(ctx)))
        end

        argsSeen = Set{Symbol}()
        optionalposition = 0
        typeRHShints = Dict{Symbol, Any}() # x = 1
        assertions = Dict{Symbol, Any}() # e.g. x::Int

        function resolveArguments(sube::Symbol, position)
            if sube in argsSeen
                msg(ctx, :E331, sube, "duplicate argument")
            elseif sube in keys(ctx.current.localvars)
                msg(ctx, :E331, sube,
                    "function argument duplicates static parameter name")
            end
            if position != 0 && optionalposition != 0
                msg(ctx, :E411, sube, "non-default argument following default arguments")
            end
            push!(argsSeen, sube)
            localset!(ctx.current, sube, VarInfo(location(ctx)))
            if isstaged
                assertions[sube] = Type
            end
            return sube
        end
        function resolveArguments(sube, position)
            # zero position means it's not called at the top level
            if isexpr(sube, :parameters)
                for (j,kw) in enumerate(sube.args)
                    if isexpr(kw, :(...))
                        if j != length(sube.args)
                            msg(ctx, :E412, kw, "named ellipsis ... can only be the last argument")
                            return
                        end
                        sym = resolveArguments(kw, 0)
                        if isa(sym, Symbol)
                            if isstaged
                                assertions[sym] = Type
                            else
                                # This may change to Array{(Symbol,Any), 1} in the future
                                assertions[sym] = Array{Any,1}
                            end
                        end
                        return
                    elseif !isexpr(kw, [:(=), :kw])
                        msg(ctx, :E423, kw, "named keyword argument must have a default")
                        return
                    else
                        resolveArguments(kw, 0)
                    end
                end
            elseif isexpr(sube, [:(=), :kw])
                if position != 0
                    optionalposition = position
                end
                lintexpr(sube.args[2], ctx)
                sym = resolveArguments(sube.args[1], 0)
                if !isstaged
                    rhstype = guesstype(sube.args[2], ctx)
                    if isa(sym, Symbol)
                        typeRHShints[sym] = rhstype
                    end
                end
            elseif isexpr(sube, :(::))
                if length(sube.args) > 1
                    sym = resolveArguments(sube.args[1], 0)
                    if !isstaged
                        if isa(sym, Symbol)
                            dt = parsetype(ctx, sube.args[2])
                            assertions[sym] = dt
                        end
                    end
                    lintfuncargtype(sube.args[2], ctx)
                    return sym
                else
                    lintfuncargtype(sube.args[1], ctx)
                end
            elseif isexpr(sube, :(...))
                if position != 0 && position != length(ex.args[1].args)
                    msg(ctx, :E413, sube, "positional ellipsis ... can only be the last argument")
                end
                sym = resolveArguments(sube.args[1], 0)
                if isa(sym, Symbol)
                    if isstaged
                        assertions[sym] = Tuple{Vararg{Type}}
                    elseif haskey(assertions, sym)
                        assertions[sym] = Tuple{Vararg{assertions[sym]}}
                    else
                        assertions[sym] = Tuple{Vararg{Any}}
                    end
                end
            elseif isexpr(sube, :($))
                lintexpr(sube.args[1], ctx)
            else
                msg(ctx, :E131, sube, "Lint does not understand argument #$(position)")
            end
            return
        end

        params = nothing
        for i = (isexpr(ex.args[1], :call) ? 2 : 1):length(ex.args[1].args)
            if isexpr(ex.args[1].args[i], :parameters)
                params = ex.args[1].args[i]
                continue
            end
            resolveArguments(ex.args[1].args[i], i)
        end
        if params != nothing
            resolveArguments(params, 1)
        end

        for s in argsSeen
            try
                vi = ctx.current.localvars[s]
                if haskey(assertions, s)
                    dt = parsetype(ctx, assertions[s])
                    vi.typeactual = dt
                    if dt != Any && haskey(typeRHShints, s) && typeRHShints[s] != Any &&
                        !(typeRHShints[s] <: dt)
                        msg(ctx, :E516, s, "type assertion and default seem inconsistent")
                    end
                elseif haskey(typeRHShints, s)
                    vi.typeactual = typeRHShints[s]
                end
            catch
            end
        end

        # TODO: deal with staged functions
        lintexpr(ex.args[2], ctx)
    end

    location!(ctx, oldloc)
    # TODO check cyclomatic complexity?
end

function lintlambda(ex::Expr, ctx::LintContext)
    # TODO: do not duplicate this code in function
    withcontext(ctx, LocalContext(ctx.current)) do
        function resolveArguments(sube)
            if isa(sube, Symbol)
                localset!(ctx.current, sube, VarInfo(location(ctx)))
            elseif sube.head == :parameters
                for kw in sube.args
                    resolveArguments(kw)
                end
            elseif isexpr(sube, Symbol[:(=), :(kw), :(::), :(...)])
                if sube.head == :(=) || sube.head == :kw
                    resolveArguments(sube.args[1])
                elseif sube.head == :(::)
                    if length(sube.args) > 1
                        resolveArguments(sube.args[1])
                    end
                elseif sube.head == :(...)
                    resolveArguments(sube.args[1])
                end
            else
                msg(ctx, :E132, sube, "Lint does not understand argument")
            end
        end

        if isa(ex.args[1], Symbol)
            resolveArguments(ex.args[1])
        elseif isexpr(ex.args[1], :tuple)
            for i = 1:length(ex.args[1].args)
                resolveArguments(ex.args[1].args[i])
            end
        else
            resolveArguments(ex.args[1])
        end
        lintexpr(ex.args[2], ctx)
    end
end

function lintfunctioncall(ex::Expr, ctx::LintContext; inthrow::Bool=false)
    if ex.args[1] == :ccall
        return  # TODO: lint ccall arguments too?
    elseif ex.args[1] == :include && ctx.quoteLvl == 0
        if isa(ex.args[2], AbstractString)
            inclfile = String(ex.args[2])
        else
            # Avoid ERROR level warnings about dynamic includes.
            msg(ctx, :I372, ex.args[2], "unable to follow non-literal include file")
            return
        end

        inclfile = joinpath(ctx.path, inclfile)

        if !ispath(inclfile)
            msg(ctx, :E311, inclfile, "cannot find include file")
            return
        else
            lintinclude(ctx, inclfile)
        end
    else
        if withincurly(ex.args[1]) == :new
            # TODO: lint uses of new?
        else
            lintexpr(ex.args[1], ctx)
        end
        func = abstract_eval(ctx, ex.args[1])

        if func !== nothing && isa(func, Type) && func <: Dict
            lintdict(ex, ctx)
        end

        if haskey(deprecated_constructors, ex.args[1])
            repl = string(deprecated_constructors[ex.args[1]])
            msg(ctx, :I481, ex.args[1], "replace $(ex.args[1])() with $(repl)()")
        end
        if ex.args[1] in not_constructible
            msg(ctx, :W441, "type $(ex.args[1]) is not constructible like this")
        elseif ex.args[1] == :(+)
            lintplus(ex, ctx)
        end

        skiplist = Int[]

        #splice! allows empty range such as 3:2, it means inserting an array
        # between position 2 and 3, without taking out any value.
        if ex.args[1] == Symbol("splice!") && Meta.isexpr(ex.args[3], :(:)) &&
            length(ex.args[3].args) == 2 && isa(ex.args[3].args[1], Real) &&
            isa(ex.args[3].args[2], Real) && ex.args[3].args[2] < ex.args[3].args[1]
            push!(skiplist, 3)
        end

        st = 2
        if ex.args[1] == :ifelse
            lintboolean(ex.args[2], ctx)
            st = 3
        end

        num_args = length(ex.args)

        if isexpr(ex.args[1], :(.))
            lintexpr(ex.args[1], ctx)
        end

        if ex.args[1] in COMPARISON_OPS
            if num_args == 3
                #reuse lintcomparison by synthetically construct the expr
                lintcomparison(Expr(:comparison, ex.args[2], ex.args[1], ex.args[3]), ctx)
            elseif !isexpr(ex.args[2], :(...))
                msg(ctx, :E132, ex.args[2], "Lint does not understand argument")
            end
        end

        if !inthrow && isa(ex.args[1], Symbol)
            s = lowercase(string(ex.args[1]))
            if occursin(s,"error") || occursin(s,"exception") || occursin(s,"mismatch") || occursin(s,"fault")
                try
                    dt = parsetype(ctx, ex.args[1])
                    if dt <: Exception && !pragmaexists( "Ignore unthrown " * string(ex.args[1]), ctx.current)
                        msg(ctx, :W448, string(ex.args[1]) * " is an Exception but it is not enclosed in a throw()")
                    end
                catch
                end
            end
        end

        if !inthrow && ex.args[1] == :throw && isexpr(ex.args[2], :call)
            lintfunctioncall(ex.args[2], ctx, inthrow=true)
        else
            for i in st:num_args
                if in(i, skiplist)
                    continue
                elseif isexpr(ex.args[i], :parameters)
                    for kw in ex.args[i].args
                        if isexpr(kw, :(...))
                            lintexpr(kw.args[1], ctx)
                        elseif ispairexpr(kw)
                            lintexpr(lexicalfirst(kw), ctx)
                            lintexpr(lexicallast(kw), ctx)
                        elseif isa(kw, Expr) && length(kw.args) == 2
                            lintexpr(kw.args[2], ctx)
                        elseif isa(kw, Symbol)
                            lintexpr(kw, ctx)
                        else
                            msg(ctx, :E133, kw, "unknown keyword pattern")
                        end
                    end
                elseif isexpr(ex.args[i], :kw)
                    lintexpr(ex.args[i].args[2], ctx)
                else
                    lintexpr(ex.args[i], ctx)
                end
            end
        end
    end
end

function lintplus(ex::Expr, ctx::LintContext)
    for i in 2:length(ex.args)
        if guesstype(ex.args[i], ctx) <: AbstractString
            msg(ctx, :E422, "string uses * to concatenate")
            break
        end
    end
end


text:
function feedforward!(net::Network,input::Vector{Float64})::Vector{Float64}
    net.layers[1] = input

    # a^(l+1) = σ(w^l*a^l + b^1)
    for i=1:(length(net.layers)-1)
        net.z[i] = net.weights[i]*net.neurons[i] + net.biases[i]
        net.neurons[i+1] = net.z[i].net.act()
    end

    return net.neurons[end]
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
using IndirectTrajOpt
using BenchmarkTools
using StaticArrays
using DifferentialEquations

# Initialize parameters
ps = initCR3BPIndirectWithSTMParams("Low Thrust CR3BP")
ps.ϵ = 0.0

# Initial condition vector
y0 = [-0.019488511458668, -0.016033479812051, 0.0,
       8.918881923678198, -4.081793688818725, 0.0, 1.0,
       15.616017, 32.875896, -0.094522,
       -0.101606, 0.044791, -0.000150, 0.133266]
y0static = SVector{14}(y0)

# Initial condition vector with STM
z0 = [y0; zeros(196)]
for i in 1:14
    z0[14 + i + (i-1)*14] = 1.0
end

# Indirect EOMs without STM
@benchmark IndirectTrajOpt.cr3bpEomIndirect($y0, $ps, 0.0)
dy0 = zeros(14)
@benchmark IndirectTrajOpt.cr3bpEomIndirect!($dy0, $y0, $ps, 0.0)

# Indirect EOMs with STM
dz0 = zeros(210)
@benchmark IndirectTrajOpt.cr3bpEomIndirectWithSTM!($dz0, $z0, $ps, 0.0)

# Integration benchmarking
tspan = (0.0, 8.6404*24*3600/ps.crp.TU)

# Without STM
@benchmark cr3bpOptIntegrate($y0static, $tspan, $ps)
@benchmark cr3bpOptIntegrate($y0static, $tspan, $ps; termCallbacks = true)
@benchmark cr3bpOptIntegrate($y0, $tspan, $ps; inPlace = true)
@benchmark cr3bpOptIntegrate($y0, $tspan, $ps; inPlace = true, termCallbacks = true)

# With STM
@benchmark cr3bpOptWithSTMIntegrate($z0, $tspan, $ps)
prob = IndirectTrajOpt.createCR3BPODEWithSTMProb(z0, tspan, ps)
@benchmark solve($prob, Vern9(); reltol = 1e-14, abstol = 1e-14, save_everystep = false,
    save_start = false, initialize_save = false, maxiters = 1e6)


text:
"""
    landcover(layers::Vector{T}; full::Bool=false, path::AbstractString="assets") where {T <: Integer}

Download and prepare the EarthEnv consensus landcover data, and returns them as
an array of `SimpleSDMPredictor`s. Layers are called by their number, from 1 to
19. The list of available layers is given in a table below. The raw data come
from https://www.earthenv.org/landcover.

THe `full` keyword indicates whether the *DISCover* information must be
included. Quoting from the reference website:

> Although DISCover is based on older remote sensing imagery (1992-1993), it
> contains some complementary information which is useful for capturing
> sub-pixel land cover heterogeneity (please see the associated article for
> details). Therefore, it is recommended to use the full version of the
> consensus land cover dataset for most applications. However, the reduced
> version may provide an alternative for applications in regions with large land
> cover change in the past two decades.

It is recommended to *keep* the content of the `path` folder, as it will
eliminate the need to download and/or extract the tiff files. For example,
calling `landcover(1:12)` will download and extract everything, and future calls
will be much faster. Please keep in mind that the layers can be quite large, so
keeping the models stored is particularly important.

| Variable |             Description              |
| -------- | ------------------------------------ |
| 1        | Evergreen/Deciduous Needleleaf Trees |
| 2        | Evergreen Broadleaf Trees            |
| 3        | Deciduous Broadleaf Trees            |
| 4        | Mixed/Other Trees                    |
| 5        | Shrubs                               |
| 6        | Herbaceous Vegetation                |
| 7        | Cultivated and Managed Vegetation    |
| 8        | Regularly Flooded Vegetation         |
| 9        | Urban/Built-up                       |
| 10       | Snow/Ice                             |
| 11       | Barren                               |
| 12       | Open Water                           |

These data are released under a CC-BY-NC license to Tuanmu & Jetz.
"""
function landcover(layer::Integer; full::Bool=false, left=nothing, right=nothing, bottom=nothing, top=nothing)
    return raster(SimpleSDMPredictor, EarthEnv(full), layer=layer, left=left, right=right, bottom=bottom, top=top)
end

landcover(layers::Vector{T}; args...) where {T <: Integer} = [landcover(l; args...) for l in layers]
landcover(layers::UnitRange{T}; args...) where {T <: Integer} = [landcover(l; args...) for l in layers]


text:
using VarianceComponentsHDFE, CSV, DataFrames, Parameters

# Precompile the compute_whole method

y = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.464866 ; 0.730622 ; 0.619239; 0.753429; 0.0863208; 0.372084 ;  0.958089]
id = [1; 1; 2;2 ; 3; 4;4 ;5;5 ;6;6]
firmid = [1;2;1;1;1;1;2;2;2;2;3]
obs_id = collect(1:11)

controls = nothing

settings_exact = VCHDFESettings(leverage_algorithm = ExactAlgorithm(), first_id_effects=true, cov_effects=true)
settings_JLA = VCHDFESettings(leverage_algorithm = JLAAlgorithm(num_simulations=5), first_id_effects=true, cov_effects=true)

# compute_whole(y,id,firmid,controls,settings_exact)
@unpack obs,  y  , first_id , second_id, controls = get_leave_one_out_set(y, id, firmid, settings_exact, controls)
@unpack θ_first, θ_second, θCOV, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = leave_out_estimation(y,first_id,second_id,controls,settings_exact)
# compute_whole(y,id,firmid,controls,settings_JLA)
@unpack θ_first, θ_second, θCOV, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = leave_out_estimation(y,first_id,second_id,controls,settings_JLA)

# Precompile reading and writing CSV files

output = DataFrame(y = y, id = first_id)
CSV.write("csv_file.csv",output)
data = DataFrame!(CSV.File("csv_file.csv"))


text:
module IsingMCMC

function genFunctor( Jd, Bd)
    Jd = Jd
    Bd = Bd

    # E(S) = - 1/2\sum_{i,j, i not j} J_{i,j} si sj - \sum_i B_i si
    function E(S, N)
        local Et = 0
        for i in 1:N
            for (j,Jv) in Jd[i]
                @inbounds Et += -S[i]*Jv*S[j]/2
            end
            @inbounds Et += -S[i]*Bd[i]
        end
        Et
    end

    # differential energy
    function dE(S, k, N)
        local dEt = 0
        for (j,Jv) in Jd[k]
            @inbounds dEt += S[k]*Jv*S[j]*2
        end
        @inbounds dEt += S[k]*Bd[k]*2
        dEt
    end

    # probability for dE
    PrE(dE, T) = exp(-dE/T)

    # flip spin at x
    function flipx!(list, x)
        list[x] *= -1
    end

    # MCMC
    function MCMC( T, N, trial)
        # initialize
        simE = zeros(Int, div(trial,1000)+1)
        S = ones(Int8, N)
        for i in 1:N
            if rand() < 0.5
                flipx!(S,i) # random flip at first
            end
        end

        # MCMC trial
        Ec = E(S,N) # current energy
        simE[1] = Ec
        local k, de
        @inbounds for t in 1:trial
            k  = rand(1:N) # Gibbs sampling position
            de = dE(S, k, N)
            if rand() < PrE(de, T) # MH criteria
                flipx!(S, k) # change k
                Ec += de     # change E
            end
            (t % 1000 == 0) && (simE[div(t,1000)+1] = Ec)
        end

        (simE, S)
    end

    return (E, MCMC)
end # of functor

end # of module


text:
using ForwardDiff

# Densities (and potentials) are represented by arrays
# ρ[ix,iy,iz,iσ] in real space, where iσ ∈ [1:n_spin_components]

function _check_positive(ρ)
    minimum(ρ) < 0 && @warn("Negative ρ detected", min_ρ=minimum(ρ))
end
function _check_total_charge(dvol, ρ, N)
    n_electrons = sum(ρ) * dvol
    if abs(n_electrons - N) > sqrt(eps(eltype(ρ)))
        @warn("Mismatch in number of electrons", sum_ρ=n_electrons, N=N)
    end
end

"""
    compute_density(basis::PlaneWaveBasis, ψ::AbstractVector, occupation::AbstractVector)

Compute the density for a wave function `ψ` discretized on the plane-wave
grid `basis`, where the individual k-points are occupied according to `occupation`.
`ψ` should be one coefficient matrix per ``k``-point.
"""
@views @timing function compute_density(basis::PlaneWaveBasis, ψ, occupation)
    T = promote_type(eltype(basis), real(eltype(ψ[1])))

    # we split the total iteration range (ik, n) in chunks, and parallelize over them
    ik_n = [(ik, n) for ik = 1:length(basis.kpoints) for n = 1:size(ψ[ik], 2)]
    chunk_length = cld(length(ik_n), Threads.nthreads())

    # chunk-local variables
    ρ_chunklocal = [zeros(T, basis.fft_size..., basis.model.n_spin_components)
                    for _ = 1:Threads.nthreads()]
    ψnk_real_chunklocal = [zeros(complex(T), basis.fft_size) for _ = 1:Threads.nthreads()]

    @sync for (ichunk, chunk) in enumerate(Iterators.partition(ik_n, chunk_length))
        Threads.@spawn for (ik, n) in chunk  # spawn a task per chunk
            ψnk_real = ψnk_real_chunklocal[ichunk]
            ρ_loc = ρ_chunklocal[ichunk]

            kpt = basis.kpoints[ik]
            G_to_r!(ψnk_real, basis, kpt, ψ[ik][:, n])
            ρ_loc[:, :, :, kpt.spin] .+= occupation[ik][n] .* basis.kweights[ik] .* abs2.(ψnk_real)
        end
    end

    ρ = sum(ρ_chunklocal)
    mpi_sum!(ρ, basis.comm_kpts)
    ρ = symmetrize_ρ(basis, ρ)

    _check_positive(ρ)
    n_elec_check = weighted_ksum(basis, sum.(occupation))
    _check_total_charge(basis.dvol, ρ, n_elec_check)

    ρ
end

# Variation in density corresponding to a variation in the orbitals and occupations.
@views @timing function compute_δρ(basis::PlaneWaveBasis{T}, ψ, δψ,
                                   occupation, δoccupation=zero.(occupation)) where T
    ForwardDiff.derivative(zero(T)) do ε
        ψ_ε   = [ψk   .+ ε .* δψk   for (ψk,   δψk)   in zip(ψ, δψ)]
        occ_ε = [occk .+ ε .* δocck for (occk, δocck) in zip(occupation, δoccupation)]
        compute_density(basis, ψ_ε, occ_ε)
    end
end

@views @timing function compute_kinetic_energy_density(basis::PlaneWaveBasis, ψ, occupation)
    T = promote_type(eltype(basis), real(eltype(ψ[1])))
    τ = similar(ψ[1], T, (basis.fft_size..., basis.model.n_spin_components))
    τ .= 0
    dαψnk_real = zeros(complex(eltype(basis)), basis.fft_size)
    for (ik, kpt) in enumerate(basis.kpoints)
        G_plus_k = [[Gk[α] for Gk in Gplusk_vectors_cart(basis, kpt)] for α in 1:3]
        for n = 1:size(ψ[ik], 2), α = 1:3
            G_to_r!(dαψnk_real, basis, kpt, im .* G_plus_k[α] .* ψ[ik][:, n])
            @. τ[:, :, :, kpt.spin] += occupation[ik][n] * basis.kweights[ik] / 2 * abs2(dαψnk_real)
        end
    end
    mpi_sum!(τ, basis.comm_kpts)
    symmetrize_ρ(basis, τ)
end

total_density(ρ) = dropdims(sum(ρ; dims=4); dims=4)
@views function spin_density(ρ)
    if size(ρ, 4) == 2
        ρ[:, :, :, 1] - ρ[:, :, :, 2]
    else
        zero(ρ[:, :, :])
    end
end

function ρ_from_total_and_spin(ρtot, ρspin=nothing)
    if ρspin === nothing
        # Val used to ensure inferability
        cat(ρtot; dims=Val(4))  # copy for consistency with other case
    else
        cat((ρtot .+ ρspin) ./ 2,
            (ρtot .- ρspin) ./ 2; dims=Val(4))
    end
end


text:
function departarrive_eccentricity_objective(e, a, Ω, i, r̄ₒ, rₛₚ, μ, v̄rel::AbstractVector{<:Real}, c)
    # for a given parking position, and velocity at SoI, and eccentricity, define a departure/arrival orbit
    rₒ = norm(r̄ₒ)
    θₒ = c * orbital_angle(rₒ, a, e)
    θₛₚ = c * orbital_angle(rₛₚ, a, e)
    ω = wrap_angle(angle_in_plane(r̄ₒ, basis_rotation(Ω, i, 0.))-θₒ)

    v̄ₛₚ = orbital_velocity(θₛₚ,a,e,μ,i,Ω,ω)

    return 1-dot(v̄ₛₚ./norm(v̄ₛₚ), v̄rel./norm(v̄rel)), ω
end

function departarrive_true_anomaly(r̄ₒ::AbstractVector{<:Real}, prim::Union{CelestialBody,Star}, v̄rel::AbstractVector{<:Real}, tₛₚ, c)
    # optimize the eccentricity of the departure/arrival orbit against the departure/arrival vector
    μ = prim.μ
    rₒ = norm(r̄ₒ)
    rₛₚ = prim.SoI
    ĥ = normalize(cross(r̄ₒ,v̄rel))
    n̂ = normalize(cross(@SVector([0.,0.,1.]), ĥ))
    if isnan(norm(n̂))
        n̂ = @SVector([1.,0.,0.])
    end
    a =  1 / (2/rₛₚ - sum(abs2,v̄rel)/μ)
    i = wrap_acos(ĥ[3])
    Ω = wrap_angle(copysign(1,n̂[2]) * acos(n̂[1]))

    if a > 0
        lb = 0.
        ub = 1. - eps(Float64)
    else
        lb = 1. + eps(Float64)
        ub = max(1-rₒ/a, rₒ/a-1)
        ub = ub == 1 ? 1. + 2*eps(Float64) : ub
    end

    objectivefun(e) = departarrive_eccentricity_objective(e, a, Ω, i, r̄ₒ, rₛₚ, μ, v̄rel, c)[1]
    res = optimize(objectivefun, lb, ub, Brent())
    e = Optim.minimizer(res)

    eccobj, ω = departarrive_eccentricity_objective(e, a, Ω, i, r̄ₒ, rₛₚ, μ, v̄rel, c)

    θₒ = c * orbital_angle(rₒ, a, e)
    θₛₚ = c * orbital_angle(rₛₚ, a, e)
    T = period(a,μ)
    Mₒ = true_to_mean(θₒ, e)
    Δt = true_to_time(θₛₚ, e, T) - true_to_time(θₒ, e, T)
    tₒ = tₛₚ - Δt

    return Orbit(a,e,i,Ω,ω,Mₒ,tₒ,prim), eccobj, 0.
end

function departarrive_true_anomaly(θ, pkorb::Orbit{<:CelestialBody}, v̄rel::AbstractVector{<:Real}, tₛₚ, c; match_pkorb_Mo=true)
    r̄ₒ = orbital_position(θ, pkorb)
    daorb, eccobj = departarrive_true_anomaly(r̄ₒ, pkorb.primary, v̄rel::AbstractVector{<:Real}, tₛₚ, c)[1:2]

    a, e, i, Ω, ω, Mₒ = daorb.a, daorb.e, daorb.i, daorb.Ω, daorb.ω, daorb.Mo

    tₒ = match_pkorb_Mo ? true_to_time(θ, pkorb, daorb.epoch - pkorb.period/2) : daorb.epoch

    v̄ₒ = time_orbital_velocity(daorb.epoch, daorb)
    v̄ₚₖ = orbital_velocity(θ, pkorb)
    Δv̄ = v̄ₒ .- v̄ₚₖ

    return Orbit(a,e,i,Ω,ω,Mₒ,tₒ,pkorb.primary), eccobj, Δv̄
end

function departarrive_objectivefun(θ, pkorb, v̄rel, tₛₚ, c)
    err, Δv̄ = departarrive_true_anomaly(θ, pkorb, v̄rel, tₛₚ, c)[2:3]
    return norm(Δv̄) * exp(1000*err)        # penalize mismatched direction
end

function departarrive_orbit(pkorb::Orbit{<:CelestialBody}, v̄rel::AbstractVector{<:Real}, tₛₚ; out=true, match_pkorb_Mo=true)
    c = out ? 1 : -1
    # optimize dv against true anomaly, with each true anomaly being optimized for eccentricity
    if pkorb.e < 1           # elliptical case: any true anomaly
        lb = -2π
        ub = 2π
    else                    # hyperbolic case: true anomalies occuring within the SoI
        ub = orbital_angle(pkorb.primary.SoI, pkorb.a, pkorb.e) - eps(Float64)
        lb = -ub
    end

    objectivefun(θ) = departarrive_objectivefun(θ, pkorb, v̄rel, tₛₚ, c)

    res = optimize(objectivefun, lb, ub, Brent())
    θₚₖ = Optim.minimizer(res)

    daorb, eccobj, Δv̄ = departarrive_true_anomaly(θₚₖ, pkorb, v̄rel, tₛₚ, c; match_pkorb_Mo=match_pkorb_Mo)

    if eccobj > 1e-3
        @warn "departure/arrival velocity mismatch from target direction: $(pkorb.primary.name), $eccobj"
    end
    return daorb, eccobj, Δv̄
end

function fast_departarrive_orbit(pkorb::Orbit{<:CelestialBody}, v̄rel::AbstractVector{<:Real}, tₛₚ; out=true,  match_pkorb_Mo=true, tol=0.1, maxit=50)
    c = out ? 1 : -1
    μ = pkorb.primary.μ
    rₛₚ = pkorb.primary.SoI
    vₛₚ = norm(v̄rel)
    err = tol+1
    it = 0
    rₒnext = pkorb.a
    a = 1 / (2/rₛₚ - vₛₚ^2/μ)

    e, θₛₚ = 0., 0.
    r̄ₒ = @SVector zeros(Float64, 3)
    v̄ₒ = @SVector zeros(Float64, 3)
    while abs(err) > tol && it < maxit
        it = it+1
        rₒ = rₒnext
        vₒ = √(vₛₚ^2 + 2(μ/rₒ - μ/rₛₚ))
        e = √(1 + 2((vₒ^2)/2 - μ/rₒ) * rₒ^2 * vₒ^2 / μ^2)
        θₛₚ = c * orbital_angle(rₛₚ, a, e)
        ϕₛₚ = orbital_direction(θₛₚ, e)
        v̄ₛₚ = vₛₚ * @SVector([cos(ϕₛₚ),
                            sin(ϕₛₚ),
                            0])

        r̄ₒ = @SVector([rₒ, 0, 0])
        v̄ₒ = @SVector([0, vₒ, 0])
        v̄relplane = inertial_to_perifocal_bases(v̄rel, pkorb)
        # R = align_vectors(v̄ₛₚ, v̄relplane)               # align velocity at SoI with input direction,
        #                                                # valid only for undefined parking orbit
        ψ₁ = atan(v̄relplane[3], √abs(vₛₚ^2 - v̄ₛₚ[1]^2 - v̄relplane[3]^2))
        # R₁ = MRP([1        0        0;                   # rotate around x-axis to match i
        #           0  cos(ψ₁) -sin(ψ₁);
        #           0  sin(ψ₁)  cos(ψ₁)])
        R₁ = MRP((1.,0.,0.,0.,cos(ψ₁),sin(ψ₁),0.,-sin(ψ₁),cos(ψ₁)))
        R1v̄ₛₚ = R₁*v̄ₛₚ
        ψ₂ = atan(v̄relplane[2], v̄relplane[1]) - atan(R1v̄ₛₚ[2], R1v̄ₛₚ[1])
        # R₂ = MRP([cos(ψ₂) -sin(ψ₂)  0;                   # rotate around z-axis to match ω
        #           sin(ψ₂)  cos(ψ₂)  0;
        #           0        0        1]))
        R₂ = MRP((cos(ψ₂),sin(ψ₂),0.,-sin(ψ₂),cos(ψ₂),0.,0.,0.,1.))

        v̄ₒ = perifocal_to_inertial_bases(R₂*R₁*v̄ₒ, pkorb)
        r̄ₒ = perifocal_to_inertial_bases(R₂*R₁*r̄ₒ, pkorb)

        δθ = angle_in_plane(r̄ₒ, pkorb)
        r̄ₒpkorb = state_vector(δθ, pkorb)[1]

        errprev = err
        err = norm(r̄ₒ) - norm(r̄ₒpkorb)
        if abs(err/errprev) > 0.9
            rₒnext = (norm(r̄ₒpkorb) + rₒ)/2
        else
            rₒnext = norm(r̄ₒpkorb)
        end
    end
    Δt = true_to_time(θₛₚ, e, period(a, μ))
    θₚₖ = angle_in_plane(r̄ₒ, pkorb)
    tₚₖ =  match_pkorb_Mo ? true_to_time(θₚₖ, pkorb, tₛₚ-Δt-pkorb.period/2) : tₛₚ - Δt
    v̄ₚₖ = orbital_velocity(θₚₖ, pkorb)
    Δv̄ = c*(v̄ₒ .- v̄ₚₖ)
    return Orbit(tₚₖ, r̄ₒ, v̄ₒ, pkorb.primary), Δv̄
end

"""
    orb = departure_orbit(pkorb, v̄rel, tₛₚ)

Computes an optimized departure orbit from the parking orbit `pkorb`, the relative velocity
at the patch `v̄rel`, the time at the patch `tₛₚ`.
"""
departure_orbit(pkorb::Orbit{<:CelestialBody}, v̄rel::AbstractVector{<:Real}, tₛₚ; kwargs...) = departarrive_orbit(pkorb, v̄rel, tₛₚ; out=true,  kwargs...)

"""
    orb = arrival_orbit(pkorb, v̄rel, tₛₚ)

Computes an optimized arrival orbit from the parking orbit `pkorb`, the relative velocity
at the patch `v̄rel`, the time at the patch `tₛₚ`.
"""
arrival_orbit(pkorb::Orbit{<:CelestialBody}, v̄rel::AbstractVector{<:Real}, tₛₚ; kwargs...)   = departarrive_orbit(pkorb, v̄rel, tₛₚ; out=false, kwargs...)

"""
    orb = fast_departure_orbit(pkorb, v̄rel, tₛₚ)

Computes a departure orbit from the parking orbit `pkorb`, the relative velocity
at the patch `v̄rel`, the time at the patch `tₛₚ`.

This function is faster than `departure_orbit`, but may return sub-optimal trajectories
for highly elliptical parking orbits.
"""
fast_departure_orbit(pkorb::Orbit{<:CelestialBody}, v̄rel::AbstractVector{<:Real}, tₛₚ; kwargs...) = fast_departarrive_orbit(pkorb, v̄rel, tₛₚ; out=true, kwargs...)

"""
    orb = fast_arrival_orbit(pkorb, v̄rel, tₛₚ)

Computes a arrival orbit from the parking orbit `pkorb`, the relative velocity
at the patch `v̄rel`, the time at the patch `tₛₚ`.

This function is faster than `arrival_orbit`, but may return sub-optimal trajectories
for highly elliptical parking orbits.
"""
fast_arrival_orbit(pkorb::Orbit{<:CelestialBody}, v̄rel::AbstractVector{<:Real}, tₛₚ; kwargs...)   = fast_departarrive_orbit(pkorb, v̄rel, tₛₚ; out=false, kwargs...)

# true anomaly/time at the SoI
patch_angle(daorb::Orbit, c) = c*orbital_angle(daorb.primary.SoI, daorb)
patch_time(daorb::Orbit, c) = true_to_time(patch_angle(daorb, c), daorb)

"""
    θ = ejection_angle(dorb)

Computes the true anomaly at SoI ejection for a departure orbit `dorb`.
"""
ejection_angle(dorb::Orbit)  = patch_angle(dorb,  1)

"""
    t = ejection_time(dorb)

Computes the time at SoI ejection for a departure orbit `dorb`.
"""
ejection_time(dorb::Orbit)   = patch_time(dorb,  1)

"""
    θ = insertion_angle(aorb)

Computes the true anomaly at SoI insertion for a arrival orbit `aorb`.
"""
insertion_angle(aorb::Orbit) = patch_angle(aorb, -1)

"""
    t = insertion_angle(aorb)

Computes the time at SoI insertion for a arrival orbit `aorb`.
"""
insertion_time(aorb::Orbit)  = patch_time(aorb, -1)

patch_position(daorb::Orbit, c) = orbital_position(patch_angle(daorb, c), daorb)
patch_velocity(daorb::Orbit, c) = orbital_velocity(patch_angle(daorb, c), daorb)
patch_state_vector(daorb::Orbit, c) = state_vector(patch_angle(daorb, c), daorb)

"""
    r̄ = ejection_position(dorb)

Computes the position vector at SoI ejection for a departure orbit `dorb`.
"""
ejection_position(dorb::Orbit) = orbital_position(ejection_angle(dorb), dorb)

"""
    v̄ = ejection_velocity(dorb)

Computes the velocity vector at SoI ejection for a departure orbit `dorb`.
"""
ejection_velocity(dorb::Orbit) = orbital_velocity(ejection_angle(dorb), dorb)

"""
    r̄, v̄ = ejection_state_vector(dorb)

Computes the state vector at SoI ejection for a departure orbit `dorb`.
"""
ejection_state_vector(dorb::Orbit) = state_vector(ejection_angle(dorb), dorb)

"""
    r̄ = insertion_position(aorb)

Computes the position vector at SoI insertion for a arrival orbit `aorb`.
"""
insertion_position(aorb::Orbit) = orbital_position(insertion_angle(aorb), aorb)

"""
    v̄= insertion_velocity(aorb)

Computes the velocity vector at SoI insertion for a arrival orbit `aorb`.
"""
insertion_velocity(aorb::Orbit) = orbital_velocity(insertion_angle(aorb), aorb)

"""
    r̄, v̄ = insertion_state_vector(aorb)

Computes the state vector at SoI insertion for a arrival orbit `aorb`.
"""
insertion_state_vector(aorb::Orbit) = state_vector(insertion_angle(aorb), aorb)