text:
@testset "alldifferent.jl" begin
    @testset "AllDifferent(::Vector{AbstractIntVar}, ::Trailer)" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(1, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        z = SeaPearl.IntVar(2, 3, "Z", trailer)
        vec = Vector{SeaPearl.IntVar}([x, y, z])

        constraint = SeaPearl.AllDifferent(vec, trailer)

        @test constraint.active.value
        @test !constraint.initialized.value
        @test constraint.nodesMin == 1
        @test constraint.numberOfVals == 3
        @test constraint in x.onDomainChange
        @test constraint in y.onDomainChange
        @test constraint in z.onDomainChange
    end
    @testset "orderEdge(::Edge)::Edge" begin
        @test SeaPearl.orderEdge(Edge(1, 2)) == Edge(1, 2)
        @test SeaPearl.orderEdge(Edge(2, 1)) == Edge(1, 2)
    end
    @testset "initializeGraphs!(::AllDifferent)::Pair{Graph{Int}, DiGraph{Int}}" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(1, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        z = SeaPearl.IntVar(2, 3, "Z", trailer)
        vec = Vector{SeaPearl.IntVar}([x, y, z])

        constraint = SeaPearl.AllDifferent(vec, trailer)
        SeaPearl.setValue!(constraint.edgesState[Edge(1, 6)], SeaPearl.removed)
        graph, digraph = SeaPearl.initializeGraphs!(constraint)

        @test Edge(1, 4) in edges(graph)
        @test Edge(1, 5) in edges(graph)
        @test Edge(2, 5) in edges(graph)
        @test Edge(2, 6) in edges(graph)
        @test Edge(3, 5) in edges(graph)
        @test Edge(3, 6) in edges(graph)
        @test ne(graph) == 6
    end
    @testset "getAllEdges(::DiGraph, ::Vector{Int})" begin
        bipartite = DiGraph(7)
        add_edge!(bipartite, 4, 1)
        add_edge!(bipartite, 5, 1)
        add_edge!(bipartite, 1, 6)
        add_edge!(bipartite, 2, 5)
        add_edge!(bipartite, 6, 2)
        add_edge!(bipartite, 3, 7)
        parents = bfs_parents(bipartite, 4; dir=:out)
        edgeset = SeaPearl.getAllEdges(bipartite, parents)

        @test length(edgeset) == 5
        @test Edge(1, 4) in edgeset
        @test Edge(1, 6) in edgeset
        @test Edge(2, 5) in edgeset
        @test Edge(2, 6) in edgeset
        @test Edge(1, 5) in edgeset
    end
    @testset "getAllEdges(::DiGraph, ::Vector{Int}, ::Vector{Int})" begin
        bipartite = DiGraph(7)
        add_edge!(bipartite, 4, 1)
        add_edge!(bipartite, 5, 1)
        add_edge!(bipartite, 1, 6)
        add_edge!(bipartite, 2, 5)
        add_edge!(bipartite, 6, 2)
        add_edge!(bipartite, 3, 7)
        edgeset = SeaPearl.getAllEdges(bipartite, [1, 2], [5, 6])

        @test length(edgeset) == 4
        @test Edge(1, 5) in edgeset
        @test Edge(1, 6) in edgeset
        @test Edge(2, 5) in edgeset
        @test Edge(2, 6) in edgeset
    end
    @testset "removeEdges!(::AllDifferent, ::Vector{Vector{Int}}, ::Graph, ::DiGraph)" begin

        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(1, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        z = SeaPearl.IntVar(1, 3, "z", trailer)
        SeaPearl.remove!(z.domain, 2)
        a = SeaPearl.IntVar(2, 4, "a", trailer)
        SeaPearl.remove!(a.domain, 3)
        b = SeaPearl.IntVar(3, 6, "b", trailer)
        c = SeaPearl.IntVar(6, 7, "c", trailer)
        vars = [x, y, z, a, b, c]
        constraint = SeaPearl.AllDifferent(vars, trailer)

        graph, digraph = SeaPearl.initializeGraphs!(constraint)
        matching = SeaPearl.Matching(6, [Pair(1, 7), Pair(2, 8), Pair(3, 9), Pair(4, 10), Pair(5, 11), Pair(6, 12)])
        for (idx, match) in enumerate(matching.matches)
            constraint.matching[idx] = SeaPearl.StateObject{Pair{Int, Int}}(match, trailer)
        end
        SeaPearl.setValue!(constraint.initialized, true)
        SeaPearl.buildDigraph!(digraph, graph, matching)
        prunedValues = Vector{Vector{Int}}(undef, constraint.numberOfVars)
        for i = 1:constraint.numberOfVars
            prunedValues[i] = Int[]
        end
        SeaPearl.removeEdges!(constraint, prunedValues, graph, digraph)

        @test isempty(prunedValues[1])
        @test isempty(prunedValues[2])
        @test isempty(prunedValues[3])
        @test prunedValues[4] == [2]
        @test Set(prunedValues[5]) == Set([3, 4])
        @test isempty(prunedValues[6])

        @test !(Edge(8, 4) in edges(digraph))
        @test !(Edge(9, 5) in edges(digraph))
        @test !(Edge(10, 5) in edges(digraph))

        @test !(Edge(8, 4) in edges(graph))
        @test !(Edge(9, 5) in edges(graph))
        @test !(Edge(10, 5) in edges(graph))
    end
    @testset "propagate!(::AllDifferent, ::Set{Constraint}, ::CPModification)" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(1, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        z = SeaPearl.IntVar(1, 3, "z", trailer)
        SeaPearl.remove!(z.domain, 2)
        a = SeaPearl.IntVar(2, 4, "a", trailer)
        SeaPearl.remove!(a.domain, 3)
        b = SeaPearl.IntVar(3, 6, "b", trailer)
        c = SeaPearl.IntVar(6, 7, "c", trailer)
        vars = Matrix{SeaPearl.AbstractIntVar}([x y z; a b c])
        constraint = SeaPearl.AllDifferent(vars, trailer)

        toPropagate = Set{SeaPearl.Constraint}([constraint])
        modif = SeaPearl.CPModification(Dict("z" => [1]))
        SeaPearl.remove!(z.domain, 1)

        res = SeaPearl.propagate!(constraint, toPropagate, modif)

        @test res
        @test !(constraint in toPropagate)
        @test modif["x"] == [2]
        @test modif["y"] == [3]
        @test modif["z"] == [1]
        @test modif["a"] == [2]
        @test 3 in modif["b"] && 4 in modif["b"]
        @test !("c" in keys(modif))

        @test length(x.domain) == 1
        @test length(y.domain) == 1
        @test length(z.domain) == 1
        @test length(a.domain) == 1
        @test length(b.domain) == 2
        @test length(c.domain) == 2

    end
    @testset "3 queens dummy" begin
        trailer = SeaPearl.Trailer()
        n = 3
        rows = Vector{SeaPearl.AbstractIntVar}(undef, n)
        for i = 1:n
            rows[i] = SeaPearl.IntVar(1, n, "row_"*string(i), trailer)
        end

        rows_plus = Vector{SeaPearl.AbstractIntVar}(undef, n)
        for i = 1:n
            rows_plus[i] = SeaPearl.IntVarViewOffset(rows[i], i, rows[i].id*"+"*string(i))
        end

        rows_minus = Vector{SeaPearl.AbstractIntVar}(undef, n)
        for i = 1:n
            rows_minus[i] = SeaPearl.IntVarViewOffset(rows[i], -i, rows[i].id*"-"*string(i))
        end

        con1 = SeaPearl.AllDifferent(rows, trailer)
        con2 = SeaPearl.AllDifferent(rows_plus, trailer)
        con3 = SeaPearl.AllDifferent(rows_minus, trailer)

        modif = SeaPearl.CPModification()
        toPropagate = Set{SeaPearl.Constraint}([con1, con2, con3])

        @test SeaPearl.propagate!(con1, toPropagate, modif)
        @test SeaPearl.propagate!(con3, toPropagate, modif)
        @test SeaPearl.propagate!(con2, toPropagate, modif)

        SeaPearl.withNewState!(trailer) do
            SeaPearl.assign!(rows[1].domain, 1)
            toPropagate = Set{SeaPearl.Constraint}([con1, con2, con3])

            @test SeaPearl.propagate!(con1, toPropagate, modif)
            @test SeaPearl.propagate!(con3, toPropagate, modif)
            @test !SeaPearl.propagate!(con2, toPropagate, modif)
        end
        modif = SeaPearl.CPModification()
        SeaPearl.withNewState!(trailer) do
            SeaPearl.assign!(rows[1].domain, 2)
            toPropagate = Set{SeaPearl.Constraint}([con1, con2, con3])

            @test SeaPearl.propagate!(con1, toPropagate, modif)
            @test SeaPearl.propagate!(con3, toPropagate, modif)
            @test !SeaPearl.propagate!(con2, toPropagate, modif)
        end
    end
    @testset "5 queens full" begin
        n=5
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        rows = Vector{SeaPearl.IntVar}(undef, n)
        for i = 1:n
            rows[i] = SeaPearl.IntVar(1, n, "row_"*string(i), trailer)
            SeaPearl.addVariable!(model, rows[i]; branchable=true)
        end

        rows_plus = Vector{SeaPearl.IntVarView}(undef, n)
        for i = 1:n
            rows_plus[i] = SeaPearl.IntVarViewOffset(rows[i], i, rows[i].id*"+"*string(i))
            #SeaPearl.addVariable!(model, rows_plus[i]; branchable=false)
        end

        rows_minus = Vector{SeaPearl.IntVarView}(undef, n)
        for i = 1:n
            rows_minus[i] = SeaPearl.IntVarViewOffset(rows[i], -i, rows[i].id*"-"*string(i))
            #SeaPearl.addVariable!(model, rows_minus[i]; branchable=false)
        end

        SeaPearl.addConstraint!(model, SeaPearl.AllDifferent(rows, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.AllDifferent(rows_plus, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.AllDifferent(rows_minus, trailer))

        variableSelection = SeaPearl.MinDomainVariableSelection{false}()
        status = SeaPearl.solve!(model; variableHeuristic=variableSelection)

        @test status == :Optimal
        @test model.statistics.numberOfSolutions == 10
    end
    @testset "7 queens full" begin
        n=7
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        rows = Vector{SeaPearl.AbstractIntVar}(undef, n)
        for i = 1:n
            rows[i] = SeaPearl.IntVar(1, n, "row_"*string(i), trailer)
            SeaPearl.addVariable!(model, rows[i]; branchable=true)
        end

        rows_plus = Vector{SeaPearl.AbstractIntVar}(undef, n)
        for i = 1:n
            rows_plus[i] = SeaPearl.IntVarViewOffset(rows[i], i, rows[i].id*"+"*string(i))
            #SeaPearl.addVariable!(model, rows_plus[i]; branchable=false)
        end

        rows_minus = Vector{SeaPearl.AbstractIntVar}(undef, n)
        for i = 1:n
            rows_minus[i] = SeaPearl.IntVarViewOffset(rows[i], -i, rows[i].id*"-"*string(i))
            #SeaPearl.addVariable!(model, rows_minus[i]; branchable=false)
        end

        SeaPearl.addConstraint!(model, SeaPearl.AllDifferent(rows, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.AllDifferent(rows_plus, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.AllDifferent(rows_minus, trailer))

        variableSelection = SeaPearl.MinDomainVariableSelection{false}()
        status = SeaPearl.solve!(model; variableHeuristic=variableSelection)

        @test status == :Optimal
        @test model.statistics.numberOfSolutions == 40
    end

    @testset "variablesArray()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(1, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        z = SeaPearl.IntVar(2, 3, "Z", trailer)
        vec = Vector{SeaPearl.IntVar}([x, y, z])

        constraint = SeaPearl.AllDifferent(vec, trailer)
        @test length(SeaPearl.variablesArray(constraint)) == 3
    end
end


text:
"""
A Julia Package for handling Spherical Geometry. Spherical geometry is the
geometry of the two-dimensional surface of a sphere. This package only
handles geometries generated using great circle sections (arcs).

Angles are specified in [deg].

Types:
* Point: specified by latitude ϕ [deg] and longitude λ [deg]
* Arc: specified as the shortest great circle between two points.
* Line: is a great circle line specified by a point and azimuth [deg]
* Arcs: is a string of continuous line sections defined by a set of points.
* Polygon: is a spherical polygon defined by a set of points.

It includes the calculation of:
* The angular distance to a point, line, arc, multi-arc, or polygon border
* The along line angular distance between a point and a line.
* The intersection points between lines, arcs, multi-arcs, and polygon borders.
* The self intersection points of multi-arcs and polygon borders
* The bounding box of a given polygon or set of arcs.
* The convexhull of a given polygon.
* The normalized point
* The (final) azimuth [deg] between two points
* The spherical angle [deg] and spherical excess [deg] between three points
* The midpoint between two points, of an arc, of an arcs
* The intermediate point at a given fraction between two points, of an arc, of an arcs
* The destination point given a start point, a direction and distance.
* The intersection points of arcs, and polygons
* The self intersection points of arcs or a polygon
* The highest/lowest latitude (point) of a great circle
* The area of a polygon/spherical triangle given a radius

And testing if:
* a point, arc, arcs or polygon is inside a polygon
* a point is on a line, arc, arcs or polygon border (within tolerance)
* a polygon/arcs is self-isselfintersecting
* a polygon is simple, complex, convex or concave
"""
module SphericalGeometry

import Base:(-)
import Base:(*)
import Base.Math.rad2deg, Base.Math.deg2rad
import LinearAlgebra.cross, LinearAlgebra.dot, LinearAlgebra.norm

include("types.jl")
include("utility.jl")
include("angles.jl")
include("distances.jl")
include("points.jl")
include("operations.jl")
include("area.jl")
end

#TODO Check https://spacetelescope.github.io/spherical_geometry/index.html for additional functions
#TODO Intersection ∩ of polygons
#TODO Union ∪ of polygons
#TODO Substraction - of polygons
#TODO isclockwise
#TODO iscounterclockwise

ext: jl


text:
#!/usr/bin/env julia

module YAMLTests

import YAML
import Base.Filesystem
using Test

const tests = [
    "spec-02-01",
    "spec-02-02",
    "spec-02-03",
    "spec-02-04",
    "spec-02-05",
    "spec-02-06",
    "spec-02-07",
    "spec-02-08",
    "spec-02-09",
    "spec-02-10",
    "spec-02-11",
    "spec-02-12",
    "spec-02-13",
    "spec-02-14",
    "spec-02-15",
    "spec-02-16",
    "spec-02-17",
    "spec-02-18",
    "spec-02-19",
    "spec-02-20",
    "spec-02-21",
    "spec-02-22",
    "spec-02-23",
    "empty_scalar",
    "no_trailing_newline",
    "windows_newlines",
    "escape_sequences",
    "issue15",
    "issue30",
    "issue36",
    "issue39",
    "cartesian",
    "ar1",
    "ar1_cartesian",
    "multi-constructor",
    "merge-01"
]

# ignore some test cases in write_and_load testing
const test_write_ignored = [
    "spec-02-17",
    "escape_sequences",
    "cartesian",
    "ar1",
    "ar1_cartesian",
    "multi-constructor"
]

function equivalent(xs::Dict, ys::Dict)
    if Set(collect(keys(xs))) != Set(collect(keys(ys)))
        return false
    end

    for k in keys(xs)
        if !equivalent(xs[k], ys[k])
            return false
        end
    end

    true
end

function equivalent(xs::AbstractArray, ys::AbstractArray)
    if length(xs) != length(ys)
        return false
    end

    for (x, y) in zip(xs, ys)
        if !equivalent(x, y)
            return false
        end
    end

    true
end

function equivalent(x::Float64, y::Float64)
    isnan(x) && isnan(y) ? true : x == y
end

function equivalent(x::AbstractString, y::AbstractString)
    while endswith(x, "\n")
        x = x[1:end-1] # trailing newline characters are ambiguous
    end
    while endswith(y, "\n")
        y = y[1:end-1]
    end
    x == y
end

function equivalent(x, y)
    x == y
end

# test custom tags
function construct_type_map(t::Symbol, constructor::YAML.Constructor,
                            node::YAML.Node)
    mapping = YAML.construct_mapping(constructor, node)
    mapping[:tag] = t
    mapping
end

function TestConstructor()
    pairs = [("!Cartesian", :Cartesian),
             ("!AR1", :AR1)]
    ret = YAML.SafeConstructor()
    for (t,s) in pairs
        YAML.add_constructor!(ret, t) do c, n
            construct_type_map(s, c, n)
        end
    end

    YAML.add_multi_constructor!(ret, "!addtag:") do constructor::YAML.Constructor, tag, node
        construct_type_map(Symbol(tag), constructor, node)
    end

    ret
end

const more_constructors = let
    pairs = [("!Cartesian", :Cartesian),
             ("!AR1", :AR1)]
    Dict{String,Function}([(t, (c, n) -> construct_type_map(s, c, n))
                           for (t, s) in pairs])
end

const multi_constructors = Dict{String, Function}(
    "!addtag:" => (c, t, n) -> construct_type_map(Symbol(t), c, n)
)
# write a file, then load its contents to be tested again
function write_and_load(data::Any)
    path = Filesystem.tempname() * ".yml" # path to a temporary file
    try
        YAML.write_file(path, data)
        return YAML.load_file(path, more_constructors)
    finally
        Filesystem.rm(path, force=true)
    end
end

const testdir = dirname(@__FILE__)
@testset for test in tests
    yamlString = open(joinpath(testdir, string(test, ".data"))) do f
        read(f, String)
    end

    expected = evalfile(joinpath(testdir, string(test, ".expected")))

    # Test Loading File with Constructor Object
    data = YAML.load_file(
        joinpath(testdir, string(test, ".data")),
        TestConstructor()
    )
    @test equivalent(data, expected)
    dictData = YAML.load_file(
        joinpath(testdir, string(test, ".data")),
        more_constructors, multi_constructors
    )
    @test equivalent(dictData, expected)

    stringData = YAML.load(yamlString, TestConstructor())
    @test equivalent(stringData, expected)
    dictStringData = YAML.load(yamlString, more_constructors, multi_constructors)
    @test equivalent(dictStringData, expected)

    @test equivalent(first(YAML.load_all_file(joinpath(testdir, string(test, ".data")), TestConstructor())), expected)
    @test equivalent(first(YAML.load_all_file(joinpath(testdir, string(test, ".data")), more_constructors, multi_constructors)), expected)

    @test equivalent(first(YAML.load_all(yamlString, TestConstructor())), expected)
    @test equivalent(first(YAML.load_all(yamlString, more_constructors, multi_constructors)), expected)
    if !in(test, test_write_ignored)
        @test equivalent(write_and_load(data), expected)
    else
        println("WARNING: I do not test the writing of $test")
    end
end

const test_errors = [
    "invalid-tag"
]

@testset "YAML Errors" "error test = $test" for test in test_errors
    @test_throws YAML.ConstructorError YAML.load_file(
        joinpath(testdir, string(test, ".data")),
        TestConstructor()
    )
end

@testset "Custom Constructor" begin

    function MySafeConstructor()
        yaml_constructors = copy(YAML.default_yaml_constructors)
        delete!(yaml_constructors, nothing)
        YAML.Constructor(yaml_constructors)
    end

    function MyReallySafeConstructor()
        yaml_constructors = copy(YAML.default_yaml_constructors)
        delete!(yaml_constructors, nothing)
        ret = YAML.Constructor(yaml_constructors)
        YAML.add_multi_constructor!(ret, nothing) do constructor::YAML.Constructor, tag, node
            throw(YAML.ConstructorError(nothing, nothing,
                "could not determine a constructor for the tag '$(tag)'",
                node.start_mark))
        end
        ret
    end

    yamlString = """
    Test: !test
        test1: !test data
        test2: !test2
            - test1
            - test2
    """

    expected = Dict{Any,Any}("Test" => Dict{Any,Any}("test2"=>["test1", "test2"],"test1"=>"data"))

    @test equivalent(YAML.load(yamlString, MySafeConstructor()), expected)
    @test_throws YAML.ConstructorError YAML.load(
        yamlString,
        MyReallySafeConstructor()
    )
end

end  # module


text:
include("./a.jl")

# [function](https://docs.julialang.org/en/v1/manual/functions/#man-functions-1) to calculate the volume of a sphere
function sphere_vol(r)
    # julia allows [Unicode names](https://docs.julialang.org/en/v1/manual/unicode-input/#Unicode-Input-1) (in UTF-8 encoding)
    # so either "pi" or the symbol π can be used
    return 4/3*pi*r^3
end

# functions can also be defined more succinctly
quadratic(a, sqr_term, b) = (-b + sqr_term) / 2a

vol = sphere_vol(3)
# @printf allows number formatting but does not automatically append the \n to statements, see below
using Printf
@printf "volume = %0.3f\n" vol
#> volume = 113.097

quad1, quad2 = quadratic2(2.0, -2.0, -12.0)
println("result 1: ", quad1)
#> result 1: 3.0
println("result 2: ", quad2)
#> result 2: -2.0

max_stars_repo_name: timm/misc

max_issues_repo_name: timm/misc

max_forks_repo_name: timm/misc


text:
default_partitioner(::Type{<:SparseMatrix1DVBC{W, Tv, Ti}}) where {W, Tv, Ti} =
    DynamicTotalChunker(model_SparseMatrix1DVBC_memory(Tv, Ti), W)

function SparseMatrix1DVBC{W}(A::SparseMatrixCSC{Tv, Ti}, method=default_partitioner(SparseMatrix1DVBC{W, Tv, Ti})) where {W, Tv, Ti}
    Φ = pack_stripe(A, method)
    return SparseMatrix1DVBC{W}(A, Φ)
end

function SparseMatrix1DVBC{W}(A::SparseMatrixCSC{Tv, Ti}, Φ::SplitPartition{Ti}) where {W, Tv, Ti}
    @inbounds begin
        (m, n) = size(A)
        hst = zeros(Ti, m + 1)
        A_pos = A.colptr
        A_idx = A.rowval
        A_val = A.nzval
        L = length(Φ)
        Φ_spl = Φ.spl
        pos = undefs(Ti, L + 1)
        ofs = undefs(Ti, L + 1)
        pos[1] = 1
        ofs[1] = 1
        for l = 1:L
            pos[l + 1] = pos[l]
            j = Φ_spl[l]
            j′ = Φ_spl[l + 1]
            for q = A_pos[j]:A_pos[j′] - 1
                i = A_idx[q]
                pos[l + 1] += (hst[i] < l)
                hst[i] = l
            end
            ofs[l + 1] = ofs[l] + (pos[l + 1] - pos[l]) * (j′ - j)
        end

        idx = Vector{Ti}(undef, pos[end] - 1)
        Δw = fld(CpuId.simdbytes(), sizeof(Tv)) #TODO need to define an upper limit on this one
        val = Vector{Tv}(undef, ofs[end] - 1 + (Δw * cld(W, Δw)))
        for q = ofs[end] : ofs[end]  - 1 + (Δw * cld(W, Δw)) #extra stuff to keep vector access in bounds
            val[q] = zero(Tv)
        end

        A_q = ones(Int, W)

        for l = 1:L
            j = Φ_spl[l]
            w = Φ_spl[l + 1] - j
            @assert w <= W
            if w == 1
                Q = pos[l]
                q = ofs[l]
                for A_q_1 = A_pos[j]:(A_pos[j + w] - 1)
                    idx[Q] = A_idx[A_q_1]
                    val[q] = A_val[A_q_1]
                    Q += 1
                    q += 1
                end
            else
                i = m + 1
                for Δj = 1:w
                    A_q[Δj] = A_pos[j + Δj - 1]
                    if A_q[Δj] < A_pos[j + Δj]
                        i = min(i, A_idx[A_q[Δj]])
                    end
                end
                Q = pos[l]
                q = ofs[l]
                while i != m + 1
                    i′ = m + 1
                    for Δj = 1:w
                        if A_q[Δj] < A_pos[j + Δj]
                            if A_idx[A_q[Δj]] == i
                                val[q] = A_val[A_q[Δj]]
                                A_q[Δj] += 1
                            else
                                val[q] = zero(Tv)
                            end
                            if A_q[Δj] < A_pos[j + Δj]
                                i′ = min(i′, A_idx[A_q[Δj]])
                            end
                        else
                            val[q] = zero(Tv)
                        end
                        q += 1
                    end
                    idx[Q] = i
                    Q += 1
                    i = i′
                end
            end
        end
        return SparseMatrix1DVBC{W, Tv, Ti}(m, n, Φ, pos, idx, ofs, val)
    end
end

function SparseMatrix1DVBC{W}(A::SparseMatrixCSC{Tv, Ti}, method::StrictChunker) where {W, Tv, Ti}
    @inbounds begin
        # matrix notation...
        # i = 1:m rows, j = 1:n columns
        m, n = size(A)

        Φ = pack_stripe(A, method)

        A_pos = A.colptr
        A_idx = A.rowval
        A_val = A.nzval

        L = length(Φ)
        Φ_spl = Φ.spl
        pos = undefs(Ti, L + 1)
        ofs = undefs(Ti, L + 1)
        pos[1] = 1
        ofs[1] = 1
        for l = 1:L
            j = Φ_spl[l]
            j′ = Φ_spl[l + 1]
            pos[l + 1] = pos[l] + A_pos[min(j + 1, j′)] - A_pos[j]
            ofs[l + 1] = A_pos[j′]
        end
        idx = Vector{Ti}(undef, pos[end] - 1)

        Δw = fld(CpuId.simdbytes(), sizeof(Tv))
        val = Vector{Tv}(undef, ofs[end] - 1 + (Δw * cld(W, Δw)))
        for q = ofs[end] : ofs[end]  - 1 + (Δw * cld(W, Δw)) #extra stuff to keep vector access in bounds
            val[q] = zero(Tv)
        end

        A_q = ones(Int, W)

        for l = 1:L
            j = Φ_spl[l]
            w = Φ_spl[l + 1] - j
            @assert w <= W
            for Q = 0 : A_pos[j + 1] - A_pos[j] - 1
                idx[pos[l] + Q] = A_idx[A_pos[j] + Q]
            end
            for j′ = Φ_spl[l] : (Φ_spl[l + 1] - 1)
                for Q = 0 : A_pos[j + 1] - A_pos[j] - 1
                    val[ofs[l] + Q * w + j′ - j] = A_val[A_pos[j′] + Q]
                end
            end
        end
        return SparseMatrix1DVBC{W, Tv, Ti}(m, n, Φ, pos, idx, ofs, val)
    end
end


text:
using Test

using DataFrames

using ScottishTaxBenefitModel
using .Randoms
using .UCTransition: route_to_uc_or_legacy, route_to_uc_or_legacy!
using .ModelHousehold
using .RunSettings
using .Results
using .Intermediate
using .STBIncomes
using .ExampleHelpers

sys = get_system( scotland=true )
settings = Settings()
settings.means_tested_routing = modelled_phase_in

@testset "Transition Tests on Example HHlds" begin
    examples = get_all_examples()
    for (hht,hh) in examples
        println( "on hhld '$hht'")
        lhh = deepcopy( hh )
        bus = get_benefit_units( lhh )
        nbus = size(bus)[1]
        for bno in eachindex(bus)
            @test typeof(bus[bno]) <: BenefitUnit
            intermed = make_intermediate(
                bno,
                bus[bno],
                sys.lmt.hours_limits,
                sys.age_limits,
                sys.child_limits,
                nbus )
            route = route_to_uc_or_legacy(
                settings,
                bus[bno],
                intermed )
        end
    end
end # example tests

@testset "Live Data Transitions" begin

    if IS_LOCAL
        rc = @timed begin
            num_households,total_num_people,nhh2 = FRSHouseholdGetter.initialise( DEFAULT_SETTINGS )
        end
        r = 0
        df = DataFrame( typ=["Job Seekers", "Disabled", "W/Children", "W/Housing", "Other", "Total"], legacy=zeros(6), uc=zeros(6))
        println( "num_households=$num_households, num_people=$(total_num_people)")
        @time for hhno in 1:num_households
            hh = FRSHouseholdGetter.get_household( hhno )
            bus = get_benefit_units(hh)
            r += 1
            intermed = make_intermediate(
                hh,
                sys.hours_limits,
                sys.age_limits,
                sys.child_limits )
            im = intermed.hhint
            if on_mt_benefits( hh ) && ! (im.someone_pension_age)
                # on actual data !! FIXME the HoC thing is actually Benefit Units
                route = route_to_uc_or_legacy(
                    settings,
                    bus[1],
                    im )
                col = route == uc_bens ? :uc : :legacy
                n = 1
                if im.num_job_seekers > 0
                    row = "Job Seekers"
                    n = im.num_adults
                elseif im.limited_capacity_for_work
                    row = "Disabled"
                elseif im.num_children > 0
                    row = "W/Children"
                elseif im.benefit_unit_number == 1
                    row = "W/Housing"
                else
                    row = "Other"
                end
                k = ( df.typ .== row )
                println( "adding |$row| |$col| weight $(hh.weight)" )
                # println( col )
                df[ k, col ] .+= hh.weight*n
                df[ (df.typ.=="Total"), col ] .+= hh.weight
            end
        end
        df[:, :sum] = df.uc.+df.legacy
        df[:, :pctuc] .= 100.0 .* df.uc./(df.uc.+df.legacy)
        print( df )
    end

end

@testset "Single HH Test - UC Route checks" begin
    #
    # chasing a bug when using UC route but contrib benefits
    # are being set to zero
    #
    examples = get_all_examples()
    hh = deepcopy( examples[cpl_w_2_children_hh])
    hhres = init_household_result( hh )
    intermed = make_intermediate(
        hh,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits
    )
    settings.means_tested_routing = uc_full
    bus = get_benefit_units( hh )
    bres = hhres.bus[1]
    tozero!( bres, ALL_INCOMES ...)
    head = get_head( bus[1])
    hpid = head.pid
    for i in BENEFITS
        bres.pers[hpid].income[i] = 1.0
    end

    route_to_uc_or_legacy!( hhres, settings, hh, intermed )
    println( inctostr( bres.pers[hpid].income ))
    for i in LEGACY_MTBS
        println( "on $(iname(i))" )
        @test bres.pers[hpid].income[i] == 0
    end
    for i in [CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE, CONTRIB_JOBSEEKERS_ALLOWANCE]
        println( "on $(iname(i)) testing for zero " )
        @test bres.pers[hpid].income[i] == 1.0
    end
end


text:
@testset "IP" begin

@testset "decode" begin
    ep = decode_ethernet(dns_packet)
    ipp = decode_ipv4(ep.payload)
    @test ipp.header.header_length == 20
    @test ipp.header.dscp == 0
    @test ipp.header.total_length == 68
    @test ipp.header.id == 0xad0b
    @test ipp.header.flags == 0
    @test ipp.header.fragment_offset == 0
    @test ipp.header.ttl == 64
    @test ipp.header.protocol == IPPROTOCOL_UDP
    @test ipp.header.src_ip == ip"172.20.2.253"
    @test ipp.header.dst_ip == ip"172.20.0.6"
    @test length(ipp.payload) == 48
end

end

max_stars_repo_name: ancapdev/NetworkProtocols

max_issues_repo_name: ancapdev/NetworkProtocols


text:
struct RedisError <: Exception
    err_type::AbstractString
    err_msg::AbstractString
end

function Base.showerror(io::IO, ex::RedisError; backtrace=true)
    printstyled(io, "$(typeof(ex)): $(ex.err_type)\n\n" * ex.err_msg * "\n", color=Base.error_color())
end


text:
module ArrayConstructors

import Compat.uninitialized

import Base: Array, Vector, Matrix

# constructors that take a value and dims
Array(x::T, dims::Dims{N}) where {T,N} = Array{T,N}(x, dims)
Array{T}(x, dims::Dims{N}) where {T,N} = Array{T,N}(x, dims)
function Array{T,N}(x, dims::Dims{N}) where {T,N}
    A = Array{T,N}(uninitialized, dims)
    for i in eachindex(A)
        A[i] = copy(x)
    end
    return A
end
# Vector and Matrix need special treatment
Vector(x::T, dims::Dims{1}) where {T} = Vector{T}(x, dims)
Matrix(x::T, dims::Dims{2}) where {T} = Matrix{T}(x, dims)

# constructors that take an iterable with size
Array(x::AbstractArray{T,N}) where {T,N} = Array{T,N}(x)
Array{T}(x::AbstractArray{S,N}) where {T,S,N} = Array{T,N}(x)
Array{T,N}(x::AbstractArray{S,N}) where {T,S,N} = Array{T,N}(uninitialized, size(x)) .= x
# Vector and Matrix need special treatment
Vector(x::AbstractVector{T}) where {T} = Vector{T}(x)
Matrix(x::AbstractMatrix{T}) where {T} = Matrix{T}(x)

# ambiguities with Base
Vector(x::Vector{T}) where {T} = Vector{T}(x)
Matrix(x::Matrix{T}) where {T} = Matrix{T}(x)

end # module


text:
export Gaussian,
        Disk,
        MRing,
        Crescent,
        ConcordanceCrescent,
        ExtendedRing,
        Ring,
        ParabolicSegment,
        Wisp

"""
$(TYPEDEF)
A type that defines it is a geometric model. These are usually
primitive models, and are usually analytic in Fourier and the image domain.
"""
abstract type GeometricModel <: AbstractModel end
@inline flux(::GeometricModel) = 1.0

@inline isprimitive(::Type{<:GeometricModel}) = IsPrimitive()

@inline visanalytic(::Type{<:GeometricModel}) = IsAnalytic()
@inline imanalytic(::Type{<:GeometricModel}) = IsAnalytic()

"""
    $(TYPEDEF)
Gaussian geometrical model.
This is a Gaussian with unit flux and standard deviation.

## Notes
To change the Gaussian flux, and shape please use the modifier functions
"""
struct Gaussian{T} <: GeometricModel end
Gaussian() = Gaussian{Float64}()
radialextent(::Gaussian) = 5.0

@inline function intensity_point(::Gaussian, x,y)
    return exp(-(x^2+y^2)/2)/2π
end

@inline function visibility_point(::Gaussian{T}, u, v, args...) where {T}
    return exp(-2π^2*(u^2 + v^2)) + zero(T)im
end

raw"""
    $(TYPEDEF)
Tophat disk geometrical model. The model is given by
```math
    I(x,y) = \begin{cases} \pi^{-1} & x^2+y^2 < 1 \\ 0 & x^2+y^2 \geq 0 \end{cases}
```
"""
struct Disk{T} <: GeometricModel end
Disk() = Disk{Float64}()

@inline function intensity_point(::Disk{T}, x, y, args...) where {T}
    r = x^2 + y^2
    return r < 1 ?  one(T)/(π) : zero(T)
end

@inline function visibility_point(::Disk{T}, x, y, args...) where {T}
    ur = 2π*(hypot(x,y) + eps(T))
    return 2*besselj1(ur)/(ur) + zero(T)im
end

radialextent(::Disk) = 3.0

"""
    $(TYPEDEF)
m-ring geometric model. This corresponds to a delta ring with a fourier expansion
in θ. The m in m-ring refers to the order of the Fourier expansion. The radius is unity.
"""
struct Ring{T} <: GeometricModel end
Ring() = Ring{Float64}()
radialextent(::Ring) = 1.5

@inline function intensity_point(::Ring{T}, x::Number, y::Number) where {T}
    r = hypot(x,y)
    θ = atan(x,y)
    dr = 1e-2
    if (abs(r-1) < dr/2)
        acc = one(T)
        return acc/(2π*dr)
    else
        return zero(T)
    end
end

@inline function visibility_point(::Ring{T}, u, v, args...) where {T}
    k = 2π*sqrt(u^2 + v^2) + eps(T)
    vis = besselj0(k) + zero(T)*im
    return vis
end

"""
    $(TYPEDEF)
m-ring geometric model. This corresponds to a delta ring with a fourier expansion
in θ. The m in m-ring refers to the order of the Fourier expansion. The radius is unity.
"""
struct MRing{T,N} <: GeometricModel
    """
    Real Fourier mode coefficients
    """
    α::NTuple{N,T}
    """
    Imaginary Fourier mode coefficients
    """
    β::NTuple{N,T}
end

function MRing{N}(α::T, β::T) where {N,T<:AbstractVector}
    (N != length(α)) && throw("N must be the length of the vector")
    S = promote_type(eltype(α), eltype(β))
    return MRing{S, N}(NTuple{N,S}(α), NTuple{N,S}(β))
end

function MRing(c::NTuple{N, <:Complex}) where {N}
    α = real.(c)
    β = imag.(c)
    return MRing(α, β)
end

radialextent(::MRing) = 1.5

@inline function intensity_point(m::MRing{T,N}, x::Number, y::Number) where {T,N}
    r = hypot(x,y)
    θ = atan(x,y)
    dr = 0.025
    if (abs(r-1) < dr/2)
        acc = one(T)
        for n in 1:N
            s,c = sincos(n*θ)
            acc += m.α[n]*c - m.β[n]*s
        end
        return acc/(2π*dr)
    else
        return zero(T)
    end
end

@inline function visibility_point(m::MRing{T,N}, u, v, args...) where {T,N}
    k = 2π*sqrt(u^2 + v^2) + eps(T)
    vis = besselj0(k) + zero(T)*im
    θ = atan(u, v)
    @inbounds for n in 1:N
        s,c = sincos(n*θ)
        vis += 2*(m.α[n]*c - m.β[n]*s)*(1im)^n*besselj(n, k)
    end
    return vis
end

"""
    $(TYPEDEF)
Creates a [Kamruddin and Dexter](https://academic.oup.com/mnras/article/434/1/765/1005984)
crescent model. This works by composing two disk models together.

# Arguments
- router: The radius of the outer disk
- rinner: The radius of the inner disk
- shift: How much the inner disk radius is shifted (positive is to the right)
- floor: The floor of the inner disk 0 means the inner intensity is zero and 1 means it is a large disk.
"""
function Crescent(router, rinner, shift, floor)
    m = stretched(Disk(), router, router)*(π*router^2) - shifted(stretched(Disk(), rinner, rinner)*((1-floor)*π*rinner^2), shift, zero(typeof(shift)))
    return m/flux(m)
end

"""
    $(TYPEDEF)
Creates the ConcordanceCrescent model, i.e. a flat-top crescent
with a displacment and a slash. Note this creates a crescent with
unit flux. If you want a different flux please use the `renomed`
modifier.

## Fields
## Notes
Unlike the Gaussian and Disk models this does not create the
unit version. In fact, this model could have been created using
the `Disk` and primitives by using Comrade.jl's model composition
functionality.
"""
struct ConcordanceCrescent{T} <: GeometricModel
    """
    Outer radius of the crescent
    """
    router::T
    """
    Inner radius of the crescent
    (i.e. inside this radius there is a hole)
    """
    rinner::T
    """
    Displacment of the inner disk radius
    """
    shift::T
    """
    Strength of the linear slash. Note that
    s∈[0.0,1.0] to ensure positivity in the image.
    """
    slash::T
end

radialextent(m::ConcordanceCrescent) = m.router*1.5

# Crescent normalization to ensure the
function _crescentnorm(m::ConcordanceCrescent)
    f = (1+m.slash)*(m.router^2 - m.rinner^2) -
        (1-m.slash)*m.shift*m.rinner*m.rinner/m.router
    return 2/π/f
end

function intensity_point(m::ConcordanceCrescent{T}, x, y, args...) where {T}
    r2 = x^2 + y ^2
    norm = _crescentnorm(m)
    if (r2 < m.router^2 && (x-m.shift)^2 + y^2 > m.rinner^2 )
        return norm/2*((1+x/m.router) + m.slash*(1-x/m.router))
    else
        return zero(T)
    end
end

function visibility_point(m::ConcordanceCrescent{T}, u, v, args...) where {T}
    k = 2π*sqrt(u^2 + v^2) + eps(T)
    norm = π*_crescentnorm(m)/k
    phaseshift = cispi(2*m.shift*u)
    b0outer,b0inner = besselj0(k*m.router), besselj0(k*m.rinner)
    b1outer,b1inner = besselj1(k*m.router), besselj1(k*m.rinner)
    b2outer,b2inner = besselj(2,k*m.router), besselj(2, k*m.rinner)

    v1 = (1+m.slash)*m.router*b1outer
    v2 = ((1+m.slash) + (1-m.slash)*m.shift/m.router)*
            phaseshift*m.rinner*b1inner
    v3 = -2im*π*u*(1-m.slash)*(m.router*b0outer -
                           m.router*b2outer -
                           2*b1outer/k
                          )/(2*k)
    v4 = -2im*π*u*(1-m.slash)*(m.rinner*b0inner -
                          m.rinner*b2inner -
                          2*b1inner/k
                         )/(2*k)*(m.rinner/m.router)*phaseshift
    return norm*(v1-v2+v3-v4)
end

"""
    $(TYPEDEF)
A symmetric extended ring whose radial profile follows an inverse
gamma distributions.

# Note
We mainly use this as an example of a non-analytic Fourier transform
(although it has a complicated expression)

# Fields

$(FIELDS)

"""
struct ExtendedRing{F} <: GeometricModel
    """shape of the radial distribution"""
    shape::F
end
visanalytic(::Type{<:ExtendedRing}) = NotAnalytic()

radialextent(m::ExtendedRing) = 6.0

function intensity_point(m::ExtendedRing, x, y)
    r = hypot(x, y) + eps()
    β = (m.shape + 1)
    α = m.shape
    β^α*r^(-α-2)*exp(-β/r)/gamma(α)/(2*π)
end

"""
    $(TYPEDEF)
A delta parabolic segment in the image domain.
The segment is centered at zero, with roots ±1 and a yintercept of 1.
"""
struct ParabolicSegment{F} <: GeometricModel end
ParabolicSegment() = ParabolicSegment{Float64}()
radialextent(::ParabolicSegment{T}) where {T} = one(T)*sqrt(2)

"""
    ParabolicSegment(a, h)
A parabolic segment with x-intercepts `±a`` and a yintercept of `h``.

# Note
This is just a convenience function for `stretched(ParabolicSegment(), a, h)`
"""
@inline function ParabolicSegment(a, h)
    # Define stretched model from unital model
    stretched(ParabolicSegment(), a, h)
end

function intensity_point(::ParabolicSegment{T}, x, y) where {T}
    #length = (√5 + asinh(2)/2)
    yw = (1-x^2)
    if abs(y - yw) < 0.01/2 && abs(x) < 1
        return 1/(2*0.01)
    else
        return zero(T)
    end
end

function visibility_point(::ParabolicSegment{T}, u, v) where {T}
    ϵ = sqrt(eps(T))
    vϵ = v + ϵ + 0im
    phase = cispi(3/4 + 2*vϵ + u^2/(2vϵ))
    #length = (√5 + asinh(2)/2)
    Δ1 = erf(√(π/(2vϵ))*cispi(1/4)*(u-2vϵ))
    Δ2 = erf(√(π/(2vϵ))*cispi(1/4)*(u+2vϵ))
    return phase/(√(2vϵ))*(Δ1-Δ2)/4
end


text:
# Test script to test Stan on Windows

try
  global CMDSTAN_HOME = @static Sys.iswindows() ? "C:\\cmdstan" : CMDSTAN_HOME
  Pkg.rm("Stan")
  Pkg.rm("Stan")
  Pkg.clone("Stan")
catch e
  println(e)
end

# Go to the Stan.jl package directory

ProjDir = joinpath(dirname(@__FILE__), "..", "Examples", "NoMamba", "Bernoulli")
cd(ProjDir) do

isdir("tmp") &&
  rm("tmp", recursive=true);

try
  include(joinpath(ProjDir, "bernoulli.jl"))
catch e
  println(e)
end

println("After include('bernoulli.jl'), CMDSTAN_HOME is set to: $(CMDSTAN_HOME)")
cd(stanmodel.tmpdir)
run(`ls`)
println()
println(stanmodel.command)
println()

end # cd


text:
const A = 𝐹(0)
const B = 𝐹(7)

oncurve(𝑥::FieldElement,𝑦::FieldElement) = 𝑦^2 == 𝑥^3 + A*𝑥 + B

xField = Union{FieldElement, Infinity}

"""
    Point{T<:xField}

`Point{T<:xField}(𝑥::T, 𝑦::T)` represents a point in an scep256k1 field.
`xField` is equal to `Union{FieldElement, Infinity}`
"""
struct Point{T<:xField}
    𝑥::T
    𝑦::T
end

Point(::Infinity,::Infinity) = Point{Infinity}(∞,∞)
Point(𝑥::FieldElement,𝑦::FieldElement) = !oncurve(𝑥,𝑦) ? throw(NotOnCurve()) : Point{FieldElement}(𝑥,𝑦)
Point(𝑥::Integer,𝑦::Integer) = !oncurve(𝐹(𝑥),𝐹(𝑦)) ? throw(NotOnCurve()) : Point{FieldElement}(𝐹(𝑥),𝐹(𝑦))

const G = Point(big"0x79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798",
                big"0x483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8")

"Formats Point{FieldElement} as `(𝑥, 𝑦)` in hexadecimal format"
function show(io::IO, z::Point{FieldElement})
    x, y = z.𝑥.𝑛, z.𝑦.𝑛
    print(io, "scep256k1 Point(𝑥,𝑦):\n", string(x, base = 16), ",\n", string(y, base = 16))
end

"Formats Point{Infinity} as `(∞, ∞)`"
function show(io::IO, z::Point{Infinity})
    print(io, "scep256k1 Point(∞, ∞)")
end

"Compares two Point, returns true if coordinates are equal"
==(x::Point, y::Point) = x.𝑥 == y.𝑥 && x.𝑦 == y.𝑦

"""
Returns the point resulting from the intersection of the curve and the
straight line defined by the points P and Q
"""
function +(𝑃::Point,𝑄::Point)
    if 𝑃.𝑥 == ∞
        return 𝑄
    elseif 𝑄.𝑥 == ∞
        return 𝑃
    elseif 𝑃.𝑥 == 𝑄.𝑥 && 𝑃.𝑦 != 𝑄.𝑦
        Point{Infinity}(∞, ∞)

    # Case 1
    elseif 𝑃.𝑥 != 𝑄.𝑥
        λ = (𝑄.𝑦 - 𝑃.𝑦) ÷ (𝑄.𝑥 - 𝑃.𝑥)
        𝑥 = λ^2 - 𝑃.𝑥 - 𝑄.𝑥
    # Case 2
    else
        λ = (3 * 𝑃.𝑥^2 + A) ÷ (2 * 𝑃.𝑦)
        𝑥 = λ^2 - 2 * 𝑃.𝑥
    end
    𝑦 = λ * (𝑃.𝑥 - 𝑥) - 𝑃.𝑦
    return Point{FieldElement}(𝑥, 𝑦)
end

"Scalar multiplication of an Point"
function *(λ::Integer,𝑃::Point)
    𝑅 = Point(∞, ∞)
    λ =  mod(λ, N)
    while λ > 0
        if λ & 1 != 0
            𝑅 += 𝑃
        end
        𝑃 += 𝑃
        λ >>= 1
    end
    return 𝑅
end

"""
    serialize(P::Point; compressed::Bool) -> Vector{UInt8}

Serialize a `Point` to its SEC format
- `compressed=true` by default
"""
function serialize(P::Point; compressed::Bool=true)
    xbin = bytes(P.𝑥.𝑛)
    if length(xbin) < 32
        prepend!(xbin, UInt8.(zeros(32 - length(xbin))))
    end
    if compressed
        if mod(P.𝑦.𝑛, 2) == 0
            prefix = 0x02
        else
            prefix = 0x03
        end
        return pushfirst!(xbin,prefix)
    else
        pushfirst!(xbin, 0x04)
        ybin = bytes(P.𝑦.𝑛)
        if length(ybin) < 32
            prepend!(ybin, UInt8.(zeros(32 - length(ybin))))
        end
        return append!(xbin, ybin)
    end
end

"""
    Point(io::IOBuffer) -> Point

Parse a SEC binary to an `Point`
"""
function Point(io::IOBuffer)
    prefix = read(io, 1)[1]
    if prefix == 4
        𝑥 = to_int(read(io, 32))
        𝑦 = to_int(read(io, 32))
        return Point(𝑥, 𝑦)
    end
    is_even = prefix == 2
    𝑥 = 𝐹(to_int(read(io, 32)))
    α = 𝑥^3 + 𝐹(B)
    β = sqrt(α)
    if mod(β.𝑛, 2) == 0
        evenβ = β
        oddβ = 𝐹(P - β.𝑛)
    else
        evenβ = 𝐹(P - β.𝑛)
        oddβ = β
    end
    if is_even
        return Point(𝑥, evenβ)
    else
        return Point(𝑥, oddβ)
    end
end

Point(sec::Vector{UInt8}) = Point(IOBuffer(sec))


text:
# Download Pfam
# =============

"""
Download a gzipped stockholm full alignment for the `pfamcode`.
The extension of the downloaded file is `.stockholm.gz` by default.
The `filename` can be changed, but the `.gz` at the end is mandatory.
"""
function downloadpfam(pfamcode::String; filename::String="$pfamcode.stockholm.gz", kargs...)
    @assert endswith(filename,".gz") "filename must end with the .gz extension"
    if occursin(r"^PF\d{5}$"i, pfamcode)
        number = pfamcode[3:end]
        download_file("http://pfam.xfam.org/family/PF$(number)/alignment/full/gzipped",
                      filename; kargs...)
    else
        throw(ErrorException("$pfamcode is not a correct Pfam code"))
    end
end


text:
#=
    Heater modelling
 =#

mutable struct Heater <: AbstractConverter
     # Parameters
     η_E_H::Float64
     lifetime::Int64
     bounds::NamedTuple{(:lb, :ub), Tuple{Float64, Float64}}
     # Initial conditions
     powerMax_ini::Float64
     soh_ini::Float64
     # Variables
     powerMax::AbstractArray{Float64,2}
     carrier::Vector{Any}
     # Eco
     cost::AbstractArray{Float64,2}
     # Inner constructor
     Heater(; η_E_H = 1.,
             lifetime = 25,
             bounds = (lb = 30., ub = 30.),
             powerMax_ini = 30.,
             soh_ini = 1.) =
             new(η_E_H, lifetime, bounds, powerMax_ini)
end

### Preallocation
function preallocate!(heater::Heater, nh::Int64, ny::Int64, ns::Int64)
     heater.powerMax = convert(SharedArray,zeros(ny+1, ns)) ; heater.powerMax[1,:] .= heater.powerMax_ini
     heater.carrier = [Electricity(), Heat()]
     heater.carrier[1].power = convert(SharedArray,zeros(nh, ny, ns))
     heater.carrier[2].power = convert(SharedArray,zeros(nh, ny, ns))
     heater.cost = convert(SharedArray,zeros(ny, ns))
     return heater
end

 ### Operation dynamic
function compute_operation_dynamics!(h::Int64, y::Int64, s::Int64, heater::Heater, decision::Float64, Δh::Int64)
    heater.carrier[1].power[h,y,s], heater.carrier[2].power[h,y,s] = compute_operation_dynamics(heater, (powerMax = heater.powerMax[y,s],), decision, Δh)
end

function compute_operation_dynamics(heater::Heater, state::NamedTuple{(:powerMax,), Tuple{Float64}}, decision::Float64, Δh::Int64)
    # Power constraint and correction
    power_E = min(max(decision, -state.powerMax), 0.)
    # Power computation
    power_H = - heater.η_E_H * power_E
    return power_E, power_H
end

### Investment dynamic
function compute_investment_dynamics!(y::Int64, s::Int64, heater::Heater, decision::Union{Float64, Int64})
    heater.powerMax[y+1,s] = compute_investment_dynamics(heater, (powerMax = heater.powerMax[y,s],), decision)
end

function compute_investment_dynamics(heater::Heater, state::NamedTuple{(:powerMax,), Tuple{Float64}}, decision::Union{Float64, Int64})
    if decision > 1e-2
        powerMax_next = decision
    else
        powerMax_next = state.powerMax
    end
    return powerMax_next
end


text:
immutable Transformation
    self :: Symbol
    deriv :: Symbol
end

function Base.show(io :: IO, m :: Transformation)
    print(io, string("Transformation(", m.self, ")"))
end

import Core.eval

function eval(transformation :: Transformation)
    eval(transformation.self)
end

function deriv(transformation :: Transformation)
    eval(transformation.deriv)
end

function sigmoid(x)
        2. / (1. + exp(-x)) - 1.
end

function dsigmoid(x)
    2. * exp(-x) / (1. + exp(-x)) ^ 2
end

const Sigmoid = Transformation(:sigmoid, :dsigmoid)

const Identity = Transformation(:identity, :one)

max_stars_repo_name: Non-Contradiction/NextGeneration

max_issues_repo_name: Non-Contradiction/NextGeneration

max_forks_repo_name: Non-Contradiction/NextGeneration


text:
# branch and bound

function remove_arc_in_routes!(routes, travel_time)
    routes_to_be_removed = []
    for r_idx in eachindex(routes)
        r = routes[r_idx]
        if calculate_path_cost(r, travel_time) == Inf
            push!(routes_to_be_removed, r_idx)
        end
    end
    deleteat!(routes, routes_to_be_removed)
end

function update_travel_time!(new_vrptw, arc::Tuple{Int, Int}, value::Int)
    i, j = arc
    n_nodes = size(new_vrptw.travel_time, 1)
    depot0 = n_nodes - 1
    depot_dummy = n_nodes

    if value == 0
        new_vrptw.travel_time[i, j] = Inf
    elseif value == 1
        for k in 1:n_nodes
            if k != i && k != depot0 && j != depot_dummy
                new_vrptw.travel_time[k, j] = Inf
            end
            if k != j && k != depot_dummy && i != depot0
                new_vrptw.travel_time[i, k] = Inf
            end
        end
    end
end

function generate_branch_instance(arc::Tuple{Int, Int}, value::Int, current_branch)
    i, j = arc

    new_vrptw = deepcopy(current_branch.vrptw_instance)
    new_routes = deepcopy(current_branch.root_routes)

    n_nodes = size(new_vrptw.travel_time, 1)

    update_travel_time!(new_vrptw, arc, value)

    remove_arc_in_routes!(new_routes, new_vrptw.travel_time)

    return new_routes, new_vrptw
end

function select_new_arc(branching_history, branch_priority, vrptw::VRPTW_Instance)
    n_nodes = size(vrptw.travel_time, 1)

    if isempty(branch_priority)
        for i in 1:n_nodes, j in 1:n_nodes
            arc = (i, j)
            if !in(arc, keys(branching_history))
                return arc, branch_priority
            end
        end
    else
        new_branch_priority = copy(branch_priority)
        arc = popfirst!(new_branch_priority)[1]
        while in(arc, keys(branching_history))
            arc = popfirst!(new_branch_priority)[1]
        end
        return arc, new_branch_priority
    end
end

function show_current(sol_y, sol_routes)
    for n in 1:length(sol_y)
        if sol_y[n] > 0.01
            @show n, sol_y[n], sol_routes[n]
        end
    end
end

function  branch_decision!(next_branches, best_sol, sol_y, sol_routes, sol_obj, new_branch)
    if isempty(sol_y)
        # stop. Infeasible

    elseif is_binary(sol_y)
        if sol_obj < best_sol.objective
            best_sol.y = sol_y
            best_sol.routes = sol_routes
            best_sol.objective = sol_obj
        end
        # stop. a binary feasible solution found.

    elseif sol_obj < best_sol.objective
        # Fractional LP Solution. Still alive. Move on
        push!(next_branches, new_branch)

    else
        # stop. Fathomed.

    end
end

function solve_BnB!(best_sol::BestIncumbent, branches; pricing_method="pulse")
    next_branches = Branch[]

    while !isempty(branches)

        # Sorting branches with the lower bound
        # The last element has the lowest bound
        sort!(branches, by= x->x.lower_bound, rev=true)
        current_branch = pop!(branches)

        # current_branch.branching_priority is updated here
        arc, new_branch_priority = select_new_arc(current_branch.history, current_branch.branch_priority, current_branch.vrptw_instance)

        sol_y, sol_routes, sol_obj = [], [], Inf
        for arc_value in 0:1

            # branching history is updated here
            new_history = copy(current_branch.history)
            new_history[arc] = arc_value

            new_routes, new_vrptw = generate_branch_instance(arc, arc_value, current_branch)

            veh_cond = current_branch.veh_cond

            if isempty(new_routes)
                sol_y, sol_routes, sol_obj = [], [], Inf
            else

                sol_y, sol_routes, sol_obj = solve_cg_rmp(new_vrptw, initial_routes=new_routes, veh_cond=veh_cond, tw_reduce=false, pricing_method=pricing_method)

                new_branch = Branch(new_history, new_vrptw, current_branch.root_routes, new_branch_priority, sol_obj, veh_cond)

                branch_decision!(next_branches, best_sol, sol_y, sol_routes, sol_obj, new_branch)

            end
        end

        if isempty(branches)
            branches = deepcopy(next_branches)
            next_branches = Branch[]
        end

    end

end

function generate_branch_priority(root_y, root_routes, vrptw)
    branch_score_dict = Dict()
    for n in eachindex(root_y)
        if root_y[n] > 0.01
            r = root_routes[n]
            for i in 1:length(r)-1
                arc = (r[i], r[i+1])
                if in(arc, keys(branch_score_dict))
                    branch_score_dict[arc] += max(1-root_y[n], root_y[n]) * vrptw.travel_time[arc[1], arc[2]]
                else
                    branch_score_dict[arc] = max(1-root_y[n], root_y[n]) * vrptw.travel_time[arc[1], arc[2]]
                end
            end
        end
    end
    return sort(collect(branch_score_dict), by= x->x[2])
end

function initial_BnB!(best_sol, vrptw, root_y, root_routes, root_obj; pricing_method="pulse")
    branch_priority = generate_branch_priority(root_y, root_routes, vrptw)
    new_vrptw = deepcopy(vrptw)
    history = OrderedDict()

    for n in eachindex(root_y)
        if root_y[n] > 1.0 - 1e-6
            r = root_routes[n]
            for i in 1:length(r)-1
                arc = (r[i], r[i+1])
                history[arc] = 1
                update_travel_time!(new_vrptw, arc, 1)
            end
        end
    end

    branches = []
    if round(sum(root_y)) == sum(root_y)
        veh_cond = ("<=", -1)
        push!(branches, Branch(history, new_vrptw, root_routes, branch_priority, root_obj, veh_cond))
    else
        new_routes = copy(root_routes)
        remove_arc_in_routes!(new_routes, new_vrptw.travel_time)
        n_vehicles = floor(sum(root_y)) |> Int
        vehicle_conditions = [("<=", n_vehicles), (">=", n_vehicles+1)]
        for veh_cond in vehicle_conditions
            sol_y, sol_routes, sol_obj = solve_cg_rmp(new_vrptw, initial_routes=new_routes, veh_cond=veh_cond, tw_reduce=false, pricing_method=pricing_method)

            new_branch = Branch(history, new_vrptw, root_routes, branch_priority, sol_obj, veh_cond)

            branch_decision!(branches, best_sol, sol_y, sol_routes, sol_obj, new_branch)
        end

    end

    solve_BnB!(best_sol, branches)

end

function complete_BnB!(best_sol, vrptw, root_y, root_routes, root_obj; pricing_method="pulse")
    # Most Priority is at the end.
    branch_priority = generate_branch_priority(root_y, root_routes, vrptw)

    branches = []
    if round(sum(root_y)) == sum(root_y)
        veh_cond = ("<=", -1)
        history = OrderedDict()
        root_branch = Branch(history, vrptw, root_routes, branch_priority, root_obj, veh_cond)
        push!(branches, root_branch)
    else
        n_vehicles = floor(sum(root_y)) |> Int
        vehicle_conditions = [("<=", n_vehicles), (">=", n_vehicles+1)]
        for veh_cond in vehicle_conditions
            sol_y, sol_routes, sol_obj = solve_cg_rmp(vrptw, initial_routes=root_routes, veh_cond=veh_cond, tw_reduce=false, pricing_method=pricing_method)

            history = OrderedDict()
            new_branch = Branch(history, vrptw, root_routes, branch_priority, sol_obj, veh_cond)

            branch_decision!(branches, best_sol, sol_y, sol_routes, sol_obj, new_branch)
        end

    end

    solve_BnB!(best_sol, branches)

end

function solveVRP(vrptw::VRPTW_Instance; tw_reduce=true, pricing_method="pulse")
    println("Pricing method = $pricing_method")

    start_time = time()

    # Time Windows Reduction
    if tw_reduce
        time_window_reduction!(vrptw)
        println("Time Windows Reduced. Cumulative Time: ", time() - start_time)
    end

    # initialization
    best_sol = BestIncumbent([], [], Inf)

    # solve root node
    @info("Solving the root LP relaxation with column generation...")
    root_y, root_routes, root_obj = solve_cg_rmp(vrptw, initial_routes=[], tw_reduce=false, pricing_method=pricing_method)
    println("Root node solved. Cumulative Time: ", time() - start_time)
    @show root_obj

    @info("Root LP relaxation solution:")
    show_current(root_y, root_routes)
    println("Is the root solution binary? ", is_binary(root_y))

    if is_binary(root_y)
        best_sol.y = root_y
        best_sol.routes = root_routes
        best_sol.objective = root_obj
    else

        # # Initial branching
        # # First branching on the number of vehicles
        # # Then on arc flow, first fixing x_ij = 1
        # # then branch-and-bound for the rest
        # initial_BnB!(best_sol, vrptw, root_y, root_routes, root_obj)

        # @info("Initial BnB is done.")
        # println("Initial BnB. Cumulative Time: ", time() - start_time)
        # @show best_sol.objective
        # show_current(best_sol.y, best_sol.routes)

        # Complete branching
        # First branching on the number of vehicles
        # Then on arc flow, from the beginning
        complete_BnB!(best_sol, vrptw, root_y, root_routes, root_obj)

        @info("Complete BnB is done.")
        println("Complete BnB. Cumulative Time: ", time() - start_time)
        @show best_sol.objective
        show_current(best_sol.y, best_sol.routes)

    end

    depot = size(vrptw.travel_time, 1) - 1
    dummy = size(vrptw.travel_time, 1)
    final_routes = []
    for i in eachindex(best_sol.y)
        if best_sol.y[i] > 1 - EPS
            route = best_sol.routes[i]
            if route[1] == depot
                route[1] = 0
            end
            if route[end] == dummy
                route[end] = 0
            end
            push!(final_routes, route)
        end
    end
    return final_routes, best_sol.objective

    # return best_sol.y, best_sol.routes, best_sol.objective
end

function solveVRP(solomon::Solomon, dists::Matrix{Float64}; tw_reduce=true, pricing_method="pulse", digits=1)
    dists_OA = OffsetArray(dists, 0:size(dists,1)-1, 0:size(dists,2)-1)
    return solveVRP(solomon, dists_OA, tw_reduce=tw_reduce, pricing_method=pricing_method, digits=digits)
end

# function solveVRPpy(solomon::Solomon; tw_reduce=true, pricing_method="pulse", digits=1)
#     return solveVRP(solomon, tw_reduce=tw_reduce, pricing_method=pricing_method, digits=digits)
# end

function solveVRP(solomon::Solomon, dists::OffsetArray{Float64,2,Array{Float64,2}}; tw_reduce=true, pricing_method="pulse", digits=1)
    vrptw = generate_solomon_vrptw_instance(solomon, dists=dists, digits=digits)
    return solveVRP(vrptw, tw_reduce=tw_reduce, pricing_method=pricing_method)
end

function solveVRP(solomon::Solomon; tw_reduce=true, pricing_method="pulse", digits=1)
    vrptw = generate_solomon_vrptw_instance(solomon, digits=digits)
    return solveVRP(vrptw, tw_reduce=tw_reduce, pricing_method=pricing_method)
end


text:
#
# """
# `split_smooth(cf::Vararg{Term}) -> (smooth, nonsmooth)`
#
# Splits cost function into `SmoothFunction` and `NonSmoothFunction` terms.
# """
# split_smooth(cf::Vararg{Term}) = cf[findall(is_smooth(cf))],cf[findall((!).(is_smooth(cf)))]
# split_smooth{N}(cf::NTuple{N,Term}) = split_smooth(cf...)
#
# """
# `split_AAc_diagonal(cf::Vararg{Term}) -> (proximable, non_proximable)`
#
# Splits cost function into terms with L'*L diagonal operator.
# """
# split_AAc_diagonal(cf::Vararg{Term}) = cf[findall(is_AAc_diagonal(cf))],cf[findall((!).(is_AAc_diagonal(cf)))]
# split_AAc_diagonal{N}(cf::NTuple{N,Term}) = split_AAc_diagonal(cf...)
#
# #""" TODO
# #`split_Quadratic(cf::Vararg{Term}) -> (quadratic, non_quadratic)`
# #
# #Splits cost function into `QuadraticFunction` and non `QuadraticFunction` terms.
# #"""

split_smooth(terms::Tuple) =
    terms[findall(is_smooth.(terms))], terms[findall((!).(is_smooth.(terms)))]

split_quadratic(terms::Tuple) =
    terms[findall(is_quadratic.(terms))], terms[findall((!).(is_quadratic.(terms)))]

split_AAc_diagonal(terms::Tuple) =
    terms[findall(is_AAc_diagonal.(terms))], terms[findall((!).(is_AAc_diagonal.(terms)))]


text:
using NeuralPDE
using Test, Flux
using StochasticDiffEq
using LinearAlgebra

d = 1
m = Chain(Dense(3 , 16, elu),Dense(16, 16, elu) , Dense(16 , 5 , elu) , Dense(5 , 1))
ensemblealg = EnsembleThreads()
γ_mu_prototype = nothing
γ_sigma_prototype = zeros(d , d , 1)
γ_phi_prototype = nothing

sdealg = EM()
tspan = (0.00 , 1.00)
trajectories = 10000
function phi(x , y_phi)
  x.^2
end
g(x , γ_sigma) = γ_sigma[: , : , 1]
f(x , γ_mu_1 ,γ_mu_2 ) = [0.00]
xspan = (0.00 , 3.00)
y_domain = KolmogorovParamDomain((0.00 , 2.00) , (0.00 , 0.00) , (0.00 , 0.00) )
dt = 0.01
dx = 0.01
dy = 0.01
opt = Flux.ADAM(1e-2)
prob = ParamKolmogorovPDEProblem(f , g , phi , xspan , tspan , d  , y_domain  ; Y_sigma_prototype = γ_sigma_prototype)
sol = solve(prob, NNParamKolmogorov(m,opt , sdealg,ensemblealg) , verbose = true, dt = 0.01,
            abstol=1e-10, dx = 0.01 , trajectories = trajectories ,  maxiters = 150 , use_gpu = false)

x_test = rand(xspan[1]:dx:xspan[2] , d , 1  , 1000)
t_test = rand(tspan[1]:dt:tspan[2] ,  1 , 1000)
γ_sigma_test = rand(0.3:dy:0.3 , d , d,  1, 1000)

function getTestData(x_test, t_test, γ_sigma_test)
  X_test  = reshape([], 3 ,0)
  for i in 1:length(t_test)
    K = vcat(t_test[i] , x_test[: , : , i] , reshape(γ_sigma_test[: , : , : , i] , d^2*(size(γ_sigma_prototype))[3] , 1))
    X_test = hcat(X_test , K)
  end
  return X_test
end

X_test = getTestData(x_test, t_test, γ_sigma_test)

function analytical(x , t , y)
  #return sum(x.^2) + t*tr(y*y) for multidimensional
  return x^2 + t*(y*y)
end
function y_analytical(X)
  y_test  = reshape([], 1 ,0)
  for i in 1:length(t_test)
    y_test = hcat(y_test , analytical(X[: , i][2] , X[: , i][1] , γ_sigma_test[: , : , 1, i][1] ) )
  end
  return y_test
end

y_test = y_analytical(X_test)

m(X_test)

@test Flux.mse(m(X_test), y_test) < 0.1


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: PHE.mso 250 2007-04-27 16:32:02Z bicca $
#*------------------------------------------------------------------
type PHE_HeatTransfer
	PHE_HeatTransfer()=begin
		new(
			positive(Dict{Symbol,Any}(
				:Brief=>"Reynolds Number",
				:Default=>100,
				:Lower=>1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Prandtl Number",
				:Default=>0.5,
				:Lower=>1e-8
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Number of Units Transference",
				:Default=>0.05,
				:Lower=>1E-10
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Stream Heat Capacity",
				:Lower=>1E-3,
				:Default=>1E3,
				:Unit=>"W/K"
			)),
			heat_trans_coeff(Dict{Symbol,Any}(
				:Brief=>"Film Coefficient",
				:Default=>1,
				:Lower=>1E-12,
				:Upper=>1E6
			)),
			flux_mass(Dict{Symbol,Any}(
				:Brief=>"Channel Mass Flux",
				:Default=>1,
				:Lower=>1E-6,
				:Symbol=>"G^{channel}"
			)),
			flux_mass(Dict{Symbol,Any}(
				:Brief=>"Ports Mass Flux",
				:Default=>1,
				:Lower=>1E-6,
				:Symbol=>"G^{ports}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Viscosity Correction",
				:Default=>1,
				:Lower=>1E-6,
				:Symbol=>"\\phi"
			)),
			[:Re,:PR,:NTU,:WCp,:hcoeff,:Gchannel,:Gports,:Phi,]
		)
	end
	Re::positive
	PR::positive
	NTU::positive
	WCp::positive
	hcoeff::heat_trans_coeff
	Gchannel::flux_mass
	Gports::flux_mass
	Phi::positive
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export PHE_HeatTransfer
function atributes(in::PHE_HeatTransfer,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
PHE_HeatTransfer(_::Dict{Symbol,Any})=begin
	newModel=PHE_HeatTransfer()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module MPI

using Libdl, Serialization
using Requires
using DocStringExtensions

export mpiexec

function serialize(x)
    s = IOBuffer()
    Serialization.serialize(s, x)
    take!(s)
end

function deserialize(x)
    s = IOBuffer(x)
    Serialization.deserialize(s)
end

primitive type SentinelPtr Sys.WORD_SIZE
end

primitive type MPIPtr Sys.WORD_SIZE
end
Base.cconvert(::Type{MPIPtr}, x::SentinelPtr) = x
Base.unsafe_convert(::Type{MPIPtr}, x::SentinelPtr) = reinterpret(MPIPtr, x)

function _doc_external(fname)
"""
- `$fname` man page: [OpenMPI](https://www.open-mpi.org/doc/current/man3/$fname.3.php), [MPICH](https://www.mpich.org/static/docs/latest/www3/$fname.html)
"""
end

try
    include(joinpath(dirname(@__DIR__), "deps","deps.jl"))
catch e
    error("MPI.jl not properly configured, please run `Pkg.build(\"MPI\")`.")
end
include("implementations.jl")
include("error.jl")
include("handle.jl")
include("info.jl")
include("comm.jl")
include("environment.jl")
include("datatypes.jl")
include("buffers.jl")
include("operators.jl")
include("pointtopoint.jl")
include("collective.jl")
include("topology.jl")
include("onesided.jl")
include("io.jl")
include("mpiexec_wrapper.jl")

include("deprecated.jl")

function __init__()

    @static if Sys.isunix()
        # need to open libmpi with RTLD_GLOBAL flag for Linux, before
        # any ccall cannot use RTLD_DEEPBIND; this leads to segfaults
        # at least on Ubuntu 15.10
        Libdl.dlopen(libmpi, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
    end

    __init__deps()

    # disable UCX memory cache, since it doesn't work correctly
    # https://github.com/openucx/ucx/issues/5061
    if !haskey(ENV, "UCX_MEMTYPE_CACHE")
        ENV["UCX_MEMTYPE_CACHE"] = "no"
    end

    # Julia multithreading uses SIGSEGV to sync thread
    # https://docs.julialang.org/en/v1/devdocs/debuggingtips/#Dealing-with-signals-1
    # By default, UCX will error if this occurs (issue #337)
    if !haskey(ENV, "UCX_ERROR_SIGNALS")
        # default is "SIGILL,SIGSEGV,SIGBUS,SIGFPE"
        ENV["UCX_ERROR_SIGNALS"] = "SIGILL,SIGBUS,SIGFPE"
    end

    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("cuda.jl")
end

end


text:
using Documenter
using GridPotentialFlow

ENV["GKSwstype"] = "nul"

makedocs(
    modules = [GridPotentialFlow],
    sitename = "GridPotentialFlow.jl",
    doctest = true,
    clean = true,
    pages = [
        "Home" => "index.md",
        "Manual" => ["manual/1.-Basic-potential-flow-problem.md",
                     "manual/2.-Potential-flow-with-an-impenetrable-surface.md",
                     "manual/3.-Non-uniqueness-and-discrete-circulation.md",
                     "manual/4.-The-Kutta-condition.md",
                     "manual/5.-Generalized-edge-conditions.md",
                     "manual/6.-Force-and-the-added-mass.md"
                     ],
        "Library" => ["public.md",
                      "private.md"]
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = MathJax2(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict()
            )
        ))
    )
)

deploydocs(
    repo = "github.com/JuliaIBPM/GridPotentialFlow.jl.git",
    target = "build"
)


text:
# record.jl
# elementary audio recording utility function

using PortAudio: PortAudioStream, PortAudioDevice, read, devices
import PortAudio as PA # terminate, initialize, get_default_input_index

export record

"""
    get_default_input_device()
Determine current system-wide default input audio device.
"""
function get_default_input_device()
    PA.terminate()
    PA.initialize()
    isempty(devices()) && throw("No audio devices.")
    return devices()[1 + PA.get_default_input_index()]
end

"""
    data, S = record(time::Real = 5; input_device, args=(1,0), chat::Bool=true)

Record `time` seconds of audio data
using `input_device` (typically defaults to the built-in microphone).

# Input
* `time` : duration; 5 seconds by default

# Option
* `input_device::PortAudioDevice = get_default_input_device()` system default
* `args` : arguments to PortAudioStream;
  `(1,0)` for single-channel input by default
* `chat` : show begin/end message? `true` by default

# Output
* `data` : Vector of length `time * S`
* `S` : `sample_rate` of input stream
"""
function record(
    time::Real = 5 ;
    input_device::PortAudioDevice = get_default_input_device(),
    args = (1,0),
    chat::Bool=true,
    print::Function = (x) -> printstyled(x * "\n"; color=:blue),
)

    input_device.input_bounds.max_channels > 0 ||
        throw("no input channels for $input_device")

    PortAudioStream(input_device, args...) do stream
        S = stream.sample_rate
        N = round(Int, time * S)
        if chat
             print("Begin $time sec recording at rate $S for $N samples from $input_device.")
        end
        buf = read(stream, N)
        data = buf.data
        if size(data) == (size(data,1), 1) # N × 1
            data = vec(data)
        end
        @assert buf.samplerate == S
        if chat
            ext = round.(extrema(data), sigdigits=2)
            print("Recording finished; extrema: $ext.")
        end
        return data, S
    end
end


text:
function randa(T, K)
    data = valloc(T, div(32, sizeof(T)), K)
    data .= rand(T, K)
    data
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot
@gsp "unset key" "set multi layout 2,2 title 'Multiplot title'" :-
@gsp :- 1 x y Z "w linespoints pt 4 ps 2" "set view 55, 65" :-
@gsp :- 2 x y Z "w points pt 3 ps 3" "set view 55, 65" :-
@gsp :- 3 x y Z "w l palette" "set view 55, 65" :-
@gsp :- 4 x y Z "w pm3d" "set view 55, 65" "set key off" :-
@gsp :- x y Z "w l lc 'white'"
save(term="png size 1600,1200", output="./code/plt3d_ex2.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
abstract type AbstractControlSystem <: AbstractTimeDependentSystem end

function namespace_controls(sys::AbstractControlSystem)
    [rename(x,renamespace(nameof(sys),nameof(x))) for x in controls(sys)]
end

function controls(sys::AbstractControlSystem,args...)
    name = last(args)
    extra_names = reduce(Symbol,[Symbol(:₊,nameof(x)) for x in args[1:end-1]])
    newname = renamespace(extra_names,name)
    rename(x,renamespace(nameof(sys),newname))(get_iv(sys))
end

function controls(sys::AbstractControlSystem,name::Symbol)
    x = get_controls(sys)[findfirst(x->nameof(x)==name,sys.ps)]
    rename(x,renamespace(nameof(sys),nameof(x)))
end

controls(sys::AbstractControlSystem) = isempty(get_systems(sys)) ? get_controls(sys) : [get_controls(sys);reduce(vcat,namespace_controls.(get_systems(sys)))]

"""
$(TYPEDEF)

A system describing an optimal control problem. This contains a loss function
and ordinary differential equations with control variables that describe the
dynamics.

# Fields
$(FIELDS)

# Example

```julia
using ModelingToolkit

@variables t x(t) v(t) u(t)
D = Differential(t)

loss = (4-x)^2 + 2v^2 + u^2
eqs = [
    D(x) ~ v
    D(v) ~ u^3
]

sys = ControlSystem(loss,eqs,t,[x,v],[u],[])
```
"""
struct ControlSystem <: AbstractControlSystem
    """The Loss function"""
    loss::Any
    """The ODEs defining the system."""
    eqs::Vector{Equation}
    """Independent variable."""
    iv::Sym
    """Dependent (state) variables. Must not contain the independent variable."""
    states::Vector
    """Control variables."""
    controls::Vector
    """Parameter variables. Must not contain the independent variable."""
    ps::Vector
    observed::Vector{Equation}
    """
    Name: the name of the system. These are required to have unique names.
    """
    name::Symbol
    """
    systems: The internal systems
    """
    systems::Vector{ControlSystem}
    """
    defaults: The default values to use when initial conditions and/or
    parameters are not supplied in `ODEProblem`.
    """
    defaults::Dict
    function ControlSystem(loss, deqs, iv, dvs, controls, ps, observed, name, systems, defaults)
        check_variables(dvs, iv)
        check_parameters(ps, iv)
        check_equations(deqs, iv)
        check_equations(observed, iv)
        new(loss, deqs, iv, dvs, controls, ps, observed, name, systems, defaults)
    end
end

function ControlSystem(loss, deqs::AbstractVector{<:Equation}, iv, dvs, controls, ps;
                       observed = [],
                       systems = ODESystem[],
                       default_u0=Dict(),
                       default_p=Dict(),
                       defaults=_merge(Dict(default_u0), Dict(default_p)),
                       name=nothing)
    name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
    if !(isempty(default_u0) && isempty(default_p))
        Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", :ControlSystem, force=true)
    end
    sysnames = nameof.(systems)
    if length(unique(sysnames)) != length(sysnames)
        throw(ArgumentError("System names must be unique."))
    end
    iv′ = value(iv)
    dvs′ = value.(dvs)
    controls′ = value.(controls)
    ps′ = value.(ps)
    defaults = todict(defaults)
    defaults = Dict(value(k) => value(v) for (k, v) in pairs(defaults))
    collect_defaults!(defaults, dvs′)
    collect_defaults!(defaults, ps′)
    ControlSystem(value(loss), deqs, iv′, dvs′, controls′,
                  ps′, observed, name, systems, defaults)
end

struct ControlToExpr
    sys::AbstractControlSystem
    states::Vector
    controls::Vector
end
ControlToExpr(@nospecialize(sys)) = ControlToExpr(sys,states(sys),controls(sys))
function (f::ControlToExpr)(O)
    !istree(O) && return O
    res = if isa(operation(O), Sym)
        # normal variables and control variables
        (any(isequal(O), f.states) || any(isequal(O), controls(f))) && return tosymbol(O)
        build_expr(:call, Any[operation(O).name; f.(arguments(O))])
    else
        build_expr(:call, Any[Symbol(operation(O)); f.(arguments(O))])
    end
end
(f::ControlToExpr)(x::Sym) = nameof(x)

function constructRadauIIA5(T::Type = Float64)
  sq6 = sqrt(convert(T, 6))
  A = [11//45-7sq6/360 37//225-169sq6/1800 -2//225+sq6/75
       37//225+169sq6/1800 11//45+7sq6/360 -2//225-sq6/75
       4//9-sq6/36 4//9+sq6/36 1//9]
  c = [2//5-sq6/10;2/5+sq6/10;1]
  α = [4//9-sq6/36;4//9+sq6/36;1//9]
  A = map(T,A)
  α = map(T,α)
  c = map(T,c)
  return DiffEqBase.ImplicitRKTableau(A,c,α,5)
end

"""
```julia
runge_kutta_discretize(sys::ControlSystem,dt,tspan;
                       tab = ModelingToolkit.constructRadauIIA5())
```

Transforms a nonlinear optimal control problem into a constrained
`OptimizationProblem` according to a Runge-Kutta tableau that describes
a collocation method. Requires a fixed `dt` over a given `timespan`.
Defaults to using the 5th order RadauIIA tableau, and altnerative tableaus
can be specified using the SciML tableau style.
"""
function runge_kutta_discretize(sys::ControlSystem,dt,tspan;
                       tab = ModelingToolkit.constructRadauIIA5())
    n = length(tspan[1]:dt:tspan[2]) - 1
    m = length(tab.α)

    sts = states(sys)
    ctr = controls(sys)
    ps = parameters(sys)
    lo = get_loss(sys)
    iv = get_iv(sys)
    f = @RuntimeGeneratedFunction(build_function([x.rhs for x in equations(sys)],sts,ctr,ps,iv,conv = ModelingToolkit.ControlToExpr(sys))[1])
    L = @RuntimeGeneratedFunction(build_function(lo,sts,ctr,ps,iv,conv = ModelingToolkit.ControlToExpr(sys)))

    var(n, i...) = var(nameof(n), i...)
    var(n::Symbol, i...) = variable(n, i..., T=FnType)
    # Expand out all of the variables in time and by stages
    timed_vars = [[var(operation(x),i)(iv) for i in 1:n+1] for x in states(sys)]
    k_vars = [[var(Symbol(:ᵏ,nameof(operation(x))),i,j)(iv) for i in 1:m, j in 1:n] for x in states(sys)]
    states_timeseries = [getindex.(timed_vars,j) for j in 1:n+1]
    k_timeseries = [[Num.(getindex.(k_vars,i,j)) for i in 1:m] for j in 1:n]
    control_timeseries = [[[var(operation(x),i,j)(iv) for x in sts] for i in 1:m] for j in 1:n]
    ps = parameters(sys)
    iv = iv

    # Calculate all of the update and stage equations
    mult = [tab.A * k_timeseries[i] for i in 1:n]
    tmps = [[states_timeseries[i] .+ mult[i][j] for j in 1:m] for i in 1:n]

    bs = [states_timeseries[i] .+ dt .* reduce(+, tab.α .* k_timeseries[i],dims=1)[1] for i in 1:n]
    updates = reduce(vcat,[states_timeseries[i+1] .~ bs[i] for i in 1:n])

    df = [[dt .* Base.invokelatest(f,tmps[j][i],control_timeseries[j][i],ps,iv) for i in 1:m] for j in 1:n]
    stages = reduce(vcat,[k_timeseries[i][j] .~ df[i][j] for i in 1:n for j in 1:m])

    # Enforce equalities in the controls
    control_equality = reduce(vcat,[control_timeseries[i][end] .~ control_timeseries[i+1][1] for i in 1:n-1])

    # Create the loss function
    losses = [Base.invokelatest(L,states_timeseries[i],control_timeseries[i][1],ps,iv) for i in 1:n]
    losses = vcat(losses,[Base.invokelatest(L,states_timeseries[n+1],control_timeseries[n][end],ps,iv)])

    # Calculate final pieces
    equalities = vcat(stages,updates,control_equality)
    opt_states = vcat(reduce(vcat,reduce(vcat,states_timeseries)),reduce(vcat,reduce(vcat,k_timeseries)),reduce(vcat,reduce(vcat,control_timeseries)))

    OptimizationSystem(reduce(+,losses, init=0),opt_states,ps,equality_constraints = equalities, name=nameof(sys))
end


text:
# script to run the forward simulations

using PyPlot
using KernelDensity
using Distributions
using DataFrames

include(pwd()*"/DWPSDE model/generate_data.jl")
include(pwd()*"/DWPSDE model/set_up.jl")
include(pwd()*"/DWPSDE model/mcmc.jl")

remove_missing_values(x) = reshape(collect(skipmissing(x)),7,:)

problem = "real data"
problem = "sim data small problem"

algorithm = "MCWM"
#algorithm = "DA"
algorithm = "ADA"

if problem == "real data" && algorithm == "MCWM"
    jobname = "gp_training_7_par_lunarc_new_data_4_coresnew_data"
elseif problem == "real data" && algorithm == "DA"
    jobname = "_dagpest7new_datada_gp_mcmc"
elseif problem == "real data" && algorithm == "ADA"
    jobname = "_dagpest7new_dataada_gp_mcmc_dt"
elseif problem == "sim data small problem" && algorithm == "MCWM"
    jobname = "gp_training_7_par_lunarc_simdata_4_coressimdata"
elseif problem == "sim data small problem" && algorithm == "DA"
    jobname = "_dagpest7simdatada_gp_mcmc"
else #problem == "sim data scaled up problem" && algorithm = "ADA"
    jobname = "_dagpest7simdataada_gp_mcmc_dt"
end

# set-up parmaters for real data or sim data
if problem == "real data"

    dt = 0.35
    dt_U = 1

    B = 1.
    A = 0.01
    A_sign = 1
    f = 0.
    g = 0.03
    theta_known =  [A;A_sign;B;f;g]

else

    dt = 0.035
    dt_U = 1

    B = 1.
    A = 0.0025
    A_sign = -1
    f = 0.
    g = 0.
    theta_known =  [A;A_sign;B;f;g]

end

problem

# load data
data_res = convert(Array,readtable("DWPSDE model/analyses/"*problem*"/output_res"*jobname*".csv"))

M, N = size(data_res)

data_param = convert(Array,readtable("DWPSDE model/analyses/"*problem*"/output_param"*jobname*".csv"))

theta_true = data_param[1:N-2]
burn_in = Int64(data_param[N-2+1])

Z = convert(Array,readtable("DWPSDE model/analyses/"*problem*"/data_used"*jobname*".csv"))
Z = Z[:,1]

if algorithm == "DA" || algorithm == "ADA"
  burn_in = 1
end

# compute paramter estimations
Theta = remove_missing_values(data_res[:,1:N-2]') # stor data in column-major order
Theta = Theta[:,burn_in+1:end]
posterior_param_est = round.(mean(Theta,2),2)

# set paramters for forward sim
N_simulations = 3
nbr_sim_steps = length(Z)
start_val = Z[1]

# high density region of posterior
dist_from_posterior_mean = zeros(size(Theta,2))

for i = 1:size(Theta,2)
    dist_from_posterior_mean[i] = norm(Theta[:,i]-posterior_param_est)
end

idx = find(x -> x < quantile(dist_from_posterior_mean, 0.25),dist_from_posterior_mean)
posterior_high_dens = Theta[:,idx]
posterior_dist = Categorical(1/size(posterior_high_dens,2)*ones(size(posterior_high_dens,2)))

# pre-allocate data martix
forward_sim = zeros(N_simulations, nbr_sim_steps+1)

# run forward simulations
for i = 1:N_simulations
    idx = rand(posterior_dist)
    theta = posterior_high_dens[:,idx]
    forward_sim[i,:] = generate_data(theta, theta_known, 1., dt, dt_U, nbr_sim_steps, start_val)[1]
end

# plot forward simulations
text_size = 25
label_size = 20

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.subplot(221)
PyPlot.plot(forward_sim[1,:])
PyPlot.subplot(222)
PyPlot.plot(forward_sim[2,:])
PyPlot.subplot(223)
PyPlot.plot(forward_sim[3,:])
PyPlot.subplot(224)
PyPlot.plot(Z, "k")
#PyPlot.xlabel("Index",fontsize=text_size)
ax[:tick_params]("both",labelsize = label_size)

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.subplot(221)
PyPlot.plt[:hist](forward_sim[1,:],50)
PyPlot.subplot(222)
PyPlot.plt[:hist](forward_sim[2,:],50)
PyPlot.subplot(223)
PyPlot.plt[:hist](forward_sim[3,:],50)
PyPlot.subplot(224)
PyPlot.plt[:hist](Z,50, color = "k")
ax[:tick_params]("both",labelsize = label_size)

max_stars_repo_path: DWPSDE model/analyses/forwardsim.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: DWPSDE model/analyses/forwardsim.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: DWPSDE model/analyses/forwardsim.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


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

mutable struct ApplicationGatewayBackendHttpSettingsPropertiesFormat <: SwaggerModel
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    cookieBasedAffinity::Any # spec type: Union{ Nothing, String } # spec name: cookieBasedAffinity
    requestTimeout::Any # spec type: Union{ Nothing, Int32 } # spec name: requestTimeout
    probe::Any # spec type: Union{ Nothing, SubResource } # spec name: probe
    authenticationCertificates::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: authenticationCertificates
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    connectionDraining::Any # spec type: Union{ Nothing, ApplicationGatewayConnectionDraining } # spec name: connectionDraining

    function ApplicationGatewayBackendHttpSettingsPropertiesFormat(;port=nothing, protocol=nothing, cookieBasedAffinity=nothing, requestTimeout=nothing, probe=nothing, authenticationCertificates=nothing, provisioningState=nothing, connectionDraining=nothing)
        o = new()
        validate_property(ApplicationGatewayBackendHttpSettingsPropertiesFormat, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(ApplicationGatewayBackendHttpSettingsPropertiesFormat, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(ApplicationGatewayBackendHttpSettingsPropertiesFormat, Symbol("cookieBasedAffinity"), cookieBasedAffinity)
        setfield!(o, Symbol("cookieBasedAffinity"), cookieBasedAffinity)
        validate_property(ApplicationGatewayBackendHttpSettingsPropertiesFormat, Symbol("requestTimeout"), requestTimeout)
        setfield!(o, Symbol("requestTimeout"), requestTimeout)
        validate_property(ApplicationGatewayBackendHttpSettingsPropertiesFormat, Symbol("probe"), probe)
        setfield!(o, Symbol("probe"), probe)
        validate_property(ApplicationGatewayBackendHttpSettingsPropertiesFormat, Symbol("authenticationCertificates"), authenticationCertificates)
        setfield!(o, Symbol("authenticationCertificates"), authenticationCertificates)
        validate_property(ApplicationGatewayBackendHttpSettingsPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ApplicationGatewayBackendHttpSettingsPropertiesFormat, Symbol("connectionDraining"), connectionDraining)
        setfield!(o, Symbol("connectionDraining"), connectionDraining)
        o
    end
end # type ApplicationGatewayBackendHttpSettingsPropertiesFormat

const _property_map_ApplicationGatewayBackendHttpSettingsPropertiesFormat = Dict{Symbol,Symbol}(Symbol("port")=>Symbol("port"), Symbol("protocol")=>Symbol("protocol"), Symbol("cookieBasedAffinity")=>Symbol("cookieBasedAffinity"), Symbol("requestTimeout")=>Symbol("requestTimeout"), Symbol("probe")=>Symbol("probe"), Symbol("authenticationCertificates")=>Symbol("authenticationCertificates"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("connectionDraining")=>Symbol("connectionDraining"))
const _property_types_ApplicationGatewayBackendHttpSettingsPropertiesFormat = Dict{Symbol,String}(Symbol("port")=>"Int32", Symbol("protocol")=>"String", Symbol("cookieBasedAffinity")=>"String", Symbol("requestTimeout")=>"Int32", Symbol("probe")=>"SubResource", Symbol("authenticationCertificates")=>"Vector{SubResource}", Symbol("provisioningState")=>"String", Symbol("connectionDraining")=>"ApplicationGatewayConnectionDraining")
Base.propertynames(::Type{ ApplicationGatewayBackendHttpSettingsPropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayBackendHttpSettingsPropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayBackendHttpSettingsPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewayBackendHttpSettingsPropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendHttpSettingsPropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendHttpSettingsPropertiesFormat[property_name]

const _allowed_ApplicationGatewayBackendHttpSettingsPropertiesFormat_protocol = ["Http", "Https"]

const _allowed_ApplicationGatewayBackendHttpSettingsPropertiesFormat_cookieBasedAffinity = ["Enabled", "Disabled"]

function check_required(o::ApplicationGatewayBackendHttpSettingsPropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayBackendHttpSettingsPropertiesFormat }, name::Symbol, val)
    if name === Symbol("protocol")
        Swagger.validate_param(name, "ApplicationGatewayBackendHttpSettingsPropertiesFormat", :enum, val, _allowed_ApplicationGatewayBackendHttpSettingsPropertiesFormat_protocol)
    end
    if name === Symbol("cookieBasedAffinity")
        Swagger.validate_param(name, "ApplicationGatewayBackendHttpSettingsPropertiesFormat", :enum, val, _allowed_ApplicationGatewayBackendHttpSettingsPropertiesFormat_cookieBasedAffinity)
    end
end


text:
using Plots

function level_set(fun; level_center=[0,0], rangeX=[-1,1], rangeY=[-1,1], number_of_levels=50, color_level=true)
    gr()

    centered_RangeX = [level_center[1] + rangeX[1], level_center[1] + rangeX[2]]
    centered_RangeY = [level_center[2] + rangeY[1], level_center[2] + rangeY[2]]

    lengthX = (centered_RangeX[2] - centered_RangeX[1]) / number_of_levels
    x = centered_RangeX[1]:lengthX:centered_RangeX[2]
    lengthY = (centered_RangeY[2] - centered_RangeY[1]) / number_of_levels
    y = centered_RangeY[1]:lengthY:centered_RangeY[2]

    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))
    Z = map(fun, X, Y)
    p1 = contour(x, y, fun, fill=color_level, c=:starrynight)

    return p1
end

function add_vector(x, y)
    return quiver!([x[1]],[x[2]], quiver=([y[1]],[y[2]]))
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

__precompile__(false)
module Gaston

export closefigure, closeall, figure,
       plot, plot!, histogram, imagesc, surf,
       printfigure, set

import Base.show

using Random
using DelimitedFiles

# before doing anything else, verify gnuplot is present on this system
try
  success(`gnuplot --version`)
catch
  error("Gaston cannot be loaded: gnuplot is not available on this system.")
end

# load files
include("gaston_types.jl")
include("gaston_aux.jl")
include("gaston_llplot.jl")
include("gaston_hilvl.jl")
include("gaston_config.jl")

# determine if running in an IJulia notebook
isjupyter = false
if isdefined(Main, :IJulia) && Main.IJulia.inited
    isjupyter = true
end

# initialize internal state
gnuplot_state = GnuplotState()

# when gnuplot_state goes out of scope, exit gnuplot
finalizer(gnuplot_exit, gnuplot_state)

# initialize default configuration
gaston_config = GastonConfig()

# load initialization file
include("gaston_init.jl")

end


text:
using Test, Yao, YaoExtensions, Yao.AD
using BitBasis, Random

@testset "kernels" begin
    @test rbf_kernel(0.5)(3, 4) == rbf_kernel(0.5)(BitStr64{5}(3), BitStr64{5}(4))
    @test brbf_kernel(0.5)(3, 4) == brbf_kernel(0.5)(BitStr64{5}(3), BitStr64{5}(4))
    for sf in [rbf_functional(0.5), StatFunctional{1}(sin∘bfloat)]
        wv = YaoExtensions.witness_vec(sf, randn(4))
        @test ndims(wv) == 1
        wv = YaoExtensions.witness_vec(sf, randn(4, 10))
        @test ndims(wv) == 2
    end
end

@testset "stats basic 2d" begin
    Random.seed!(2)
    for V in [rbf_functional(0.5), StatFunctional{1}(sin∘bfloat)]
        res = expect(V, zero_state(2) |> put(2,1=>H) |> as_weights)
        @test res isa Float64
        res2 = expect(V, measure(zero_state(2) |> put(2,1=>H); nshots=1000))
        @test res2 isa Float64
        @test isapprox(res, res2, atol=1e-2)

        res3 = expect(V, zero_state(2, nbatch=10) |> put(2,1=>H) |> as_weights)
        @test length(res3) == 10 && all(res3 .≈ res)

        res4 = expect(V, measure(zero_state(2, nbatch=10) |> put(2,1=>H); nshots=1000))
        @test length(res3) == 10 && all(isapprox.(res3, res, atol=1e-2))
    end
end

@testset "stat diff 2d" begin
    Random.seed!(2)
    nbit = 3
    # 2D
    x = 0:1<<nbit-1
    V = rbf_functional(0.5)
    c = variational_circuit(nbit)
    dispatch!(c, :random)
    dbs = collect_blocks(Diff,c)

    for nbatch=[1, 10]
        p0 = zero_state(nbit; nbatch=nbatch) |> c |> probs
        sample0 = measure(zero_state(nbit) |> c; nshots=5000)
        loss0 = expect(V, p0 |> as_weights)
        gradsn = numdiff(c->expect(V, zero_state(nbit; nbatch=nbatch) |> c |> as_weights), c)
        gradse = faithful_grad(V, zero_state(nbit; nbatch=nbatch) => c)
        gradsa = expect'(V, zero_state(nbit; nbatch=nbatch) => c)[2]
        @test isapprox(gradse, gradsn, atol=1e-4)
        @test isapprox(sum.(gradse), gradsa, atol=1e-4)
    end
end

@testset "stat diff 1d" begin
    Random.seed!(2)
    nbit = 3
    # 1D
    V = StatFunctional{1}(sin ∘ bfloat)
    c = variational_circuit(nbit)
    dispatch!(c, :random)

    for nbatch=[1, 100]
        p0 = zero_state(nbit; nbatch=nbatch) |> c |> probs
        loss0 = expect(V, p0 |> as_weights)
        gradsn = numdiff(c->expect(V, zero_state(nbit; nbatch=nbatch) |> c |> as_weights), c)
        gradse = faithful_grad(V, zero_state(nbit; nbatch=nbatch) => c)
        gradsa = expect'(V, zero_state(nbit; nbatch=nbatch) => c)[2]
        @test isapprox(gradse, gradsn, atol=1e-4)
        @test isapprox(sum.(gradse), gradsa, atol=1e-4)
    end
end

@testset "mmd loss" begin
    Random.seed!(4)
    nbit = 3
    # 2D
    x = 0:1<<nbit-1
    V = rbf_mmd_loss(0.5, normalize!(rand(1<<nbit)))
    c = variational_circuit(nbit)
    dispatch!(c, :random)
    dbs = collect_blocks(Diff,c)

    for nbatch=[1, 10]
        p0 = zero_state(nbit; nbatch=nbatch) |> c |> probs
        sample0 = measure(zero_state(nbit) |> c; nshots=5000)
        loss0 = expect(V, p0 |> as_weights)
        gradsn = numdiff(c->expect(V, zero_state(nbit; nbatch=nbatch) |> c |> as_weights), c)
        gradse = faithful_grad(V, zero_state(nbit; nbatch=nbatch) => c)
        gradsa = expect'(V, zero_state(nbit; nbatch=nbatch) => c)[2]
        @test isapprox(gradse, gradsn, atol=1e-3)
        @test isapprox(sum.(gradse), gradsa, atol=1e-3)
    end
end


text:
using Random
using LinearAlgebra

function uniformOnSphere(rng = Random.GLOBAL_RNG)
    phi = 2.0 * pi * rand(rng)
    z = 2.0 * rand(rng) - 1.0;
    r = sqrt(1.0 - z * z)
    return (r * cos(phi), r * sin(phi), z)
end

function exchangeEnergy(s1, M::InteractionMatrix, s2)::Float64
    return s1[1] * (M.m11 * s2[1] + M.m12 * s2[2] + M.m13 * s2[3]) + s1[2] * (M.m21 * s2[1] + M.m22 * s2[2] + M.m23 * s2[3]) + s1[3] * (M.m31 * s2[1] + M.m32 * s2[2] + M.m33 * s2[3])
end

function getEnergy(lattice::Lattice{D,N})::Float64 where {D,N}
    energy = 0.0

    for site in 1:length(lattice)
        s0 = getSpin(lattice, site)

        #two-spin interactions
        interactionSites = getInteractionSites(lattice, site)
        interactionMatrices = getInteractionMatrices(lattice, site)
        for i in 1:length(interactionSites)
            if site > interactionSites[i]
                energy += exchangeEnergy(s0, interactionMatrices[i], getSpin(lattice, interactionSites[i]))
            end
        end

        #onsite interaction
        energy += exchangeEnergy(s0, getInteractionOnsite(lattice, site), s0)

        #field interaction
        energy += dot(s0, getInteractionField(lattice, site))
    end

    return energy
end

function getEnergyDifference(lattice::Lattice{D,N}, site::Int, newState::Tuple{Float64,Float64,Float64})::Float64 where {D,N}
    dE = 0.0
    oldState = getSpin(lattice, site)
    ds = newState .- oldState

    #two-spin interactions
    interactionSites = getInteractionSites(lattice, site)
    interactionMatrices = getInteractionMatrices(lattice, site)
    for i in 1:length(interactionSites)
        dE += exchangeEnergy(ds, interactionMatrices[i], getSpin(lattice, interactionSites[i]))
    end

    #onsite interaction
    interactionOnsite = getInteractionOnsite(lattice, site)
    dE += exchangeEnergy(newState, interactionOnsite, newState) - exchangeEnergy(oldState, interactionOnsite, oldState)

    #field interaction
    dE += dot(ds, getInteractionField(lattice, site))

    return dE
end

function getMagnetization(lattice::Lattice{D,N}) where {D,N}
    mx, my, mz = 0.0, 0.0, 0.0
    for i in 1:length(lattice)
        spin = getSpin(lattice, i)
        mx += spin[1]
        my += spin[2]
        mz += spin[3]
    end
    return [mx, my, mz] / length(lattice)
end

function getCorrelation(lattice::Lattice{D,N}, spin::Int = 1) where {D,N}
    corr = zeros(length(lattice))
    s0 = getSpin(lattice, spin)
    for i in 1:length(lattice)
        corr[i] = dot(s0, getSpin(lattice, i))
    end
    return corr
end


text:
using Test

input = joinpath(@__DIR__, "input")
thing = map(x -> parse(Int, x), readlines(input))

preamble = thing[1:25]
numbers = thing[26:end]

function calculate(n, preamble)
    for m in preamble
        if (n - m) in preamble
            return true
        end
    end
    return false
end

function calibrate(numbers, preamble)
    for n in numbers
        if !calculate(n, preamble)
            return n
        end
        push!(preamble, n)
        popfirst!(preamble)
    end
end

p1 = calibrate(numbers, preamble)

function crack(target, seq)
    for i in 1:length(seq)
        for j in i:length(seq)
            slice = seq[i:j]
            s = sum(slice)
            if s == target
                return minimum(slice) + maximum(slice)
            elseif s > target
                break
            end
        end
    end
end

p2 = crack(p1, thing)

@assert(p1 == 23278925)
@assert(p2 == 4011064)

sample = [
    35,
    20,
    15,
    25,
    47,
    40,
    62,
    55,
    65,
    95,
    102,
    117,
    150,
    182,
    127,
    219,
    299,
    277,
    309,
    576
]

@test calibrate(sample[6:end], sample[1:5]) == 127
@test crack(127, sample) == 62

println("-----------------------------------------------------------------------")
println("encoding error -- part one :: $p1")
println("encoding error -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
export IterativeSampler

mutable struct IterativeSampler{N,P,IC,EC,S,SC,V} <: AbstractIterativeSampler
    net::N
    problem::P
    itercache::IC
    sampled_values::EC

    sampler::S
    sampler_cache::SC
    v::V
end

"""
    IterativeSampler(net, sampler, algorithm, problem)

Create a single-thread iterative sampler for the quantities defined
by algorithm.
"""
function IterativeSampler(net,
                          sampl,
                          prob,
                          algo=prob)

    evaluated_vals = EvaluatedNetwork(algo, net)
    itercache      = SamplingCache(algo, prob, net)
    v              = state(prob, net)
    sampler_cache  = init_sampler!(sampl, net, v)

    IterativeSampler(net, prob, itercache, evaluated_vals, sampl, sampler_cache, v)
end

"""
    sample!(is::IterativeSampler)

Samples the quantities accordingly, returning the sampled values and sampling history.
"""
function sample!(is::IterativeSampler)
    init_sampler!(is.sampler, is.net, is.v, is.sampler_cache)
    zero!(is.itercache);
    while true
        if !(is.sampler isa FullSumSampler)
            sample_network!(is.itercache, is.problem, is.net, is.v)
        else
            sample_network_wholespace!(is.itercache, is.problem, is.net, is.v)
        end
        !samplenext!(is.v, is.sampler, is.net, is.sampler_cache) && break
    end

    evaluation_post_sampling!(is.sampled_values, is.itercache)
    return is.sampled_values
end

Base.show(io::IO, is::IterativeSampler) = print(io,
    "IterativeSampler for :"*
    "\n\tnet\t\t: $(is.net)"*
    "\n\tproblem\t: $(is.problem)"*
    "\n\tsampler\t: $(is.sampler)")