text:
# Include everything common between our LLVM versions (That's a lot)
include("../llvm_common.jl")

# Build the tarballs, then upload to Yggdrasil releases

name = "LLVMBootstrap"
version = v"9.0.1"
sources, script, products, dependencies = llvm_build_args(;version=version)
build_info = build_tarballs(ARGS, name, version, sources, script, [host_platform], products, dependencies; skip_audit=true)
upload_and_insert_shards("JuliaPackaging/Yggdrasil", name, version, build_info)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using HDF5
function simpleRead(filePath,exIdx)
    fid = h5open(filePath,"r")
    example = fid["$(exIdx)"]
    images = read(example["images"])
    utterences = read(example["utterences"])
    objects = read(example["objects"])
    coords = read(example["coords"])
    scene_id = read(example["scene_id"])
    utt_tokenized = [split(t) for t in utterences]
    lengths = [length(t) for t in utt_tokenized]
    (coords,images,objects,scene_id,utterences,utt_tokenized,lengths)
end
coords,images,objects,scene_id,utterences,utt_tokenized,lengths = simpleRead("codraw_val.h5",150);

f = open("glove_codraw_iclevr.txt");
lines = readlines(f)
glove = Dict()
for line in lines
    splitline = split(line)
    word = splitline[1]
    embedding = [parse(Float32,val) for val in splitline[2:end]]
    glove[word] = embedding
end

turns_word_embeddings = zeros((length(utterences), maximum(lengths), 300))

for (i, turn) in enumerate(utt_tokenized)
    for (j, w) in enumerate(turn)
        turns_word_embeddings[i, j,:] = glove[w]
    end
end

ext: jl


text:
function test_parse(input, args... ; broken=false)
    if broken
        @test_broken texparse(input) == manual_texexpr((:group, args...))
    else
        @test texparse(input) == manual_texexpr((:group, args...))
    end
end

@testset "Parser" begin
    @testset "Accent" begin
        test_parse(
            raw"\vec{a}",
            (:accent, raw"\vec", 'a'),
            broken=true
        )
        test_parse(
            raw"\dot{\vec{x}}",
            (:accent, raw"\dot", (:accent, raw"\vec", 'x')),
            broken=true
        )
    end

    @testset "Delimiter" begin
        test_parse(
            raw"\left( x \right)",
            (:delimited, '(', 'x', ')')
        )

        test_parse(
            raw"\left[ y \right)",
            (:delimited, '[', 'y', ')')
        )

        test_parse(
            raw"\left( a + b \right)",
            (:delimited,
                '(',
                (:group, 'a', (:spaced, '+'), 'b'),
                ')'
            )
        )

        @test_throws TeXParseError texparse(raw"\left( x")
        @test_throws TeXParseError texparse(raw"x \right)")

        # TODO Recursive delim
    end

    @testset "Fraction" begin
        test_parse(raw"\frac{1}{2}", (:frac, '1', '2'))
    end

    @testset "Function" begin
        # Check that braces are not added to the function name
        test_parse(raw"\sin{x}", (:function, "sin"), 'x')
        test_parse(
            raw"\exp(x)",
            (:function, "exp"), '(', 'x', ')'
        )
        test_parse(
            raw"\tan\alpha",
            (:function, "tan"), (:symbol, 'α', raw"\alpha")
        )
    end

    @testset "Group" begin
        test_parse(raw"{x}", 'x')
        test_parse(raw"{fgh}", (:group, 'f', 'g', 'h'))
        test_parse(
            raw"{a{b{cd}}}",
            (:group, 'a', (:group, 'b', (:group, 'c', 'd')))
        )

        @test_throws TeXParseError texparse(raw"{v")
        @test_throws TeXParseError texparse(raw"w}")
    end

    @testset "Integral" begin
        test_parse(
            raw"\int",
            (:integral, (:symbol, '∫', raw"\int"), nothing, nothing)
        )
        test_parse(
            raw"\int_a^b",
            (:overunder, (:symbol, '∫', raw"\int"), 'a', 'b')
        )
    end

    @testset "Overunder" begin
        test_parse(
            raw"\sum",
            (:underover, (:symbol, '∑', raw"\sum"), nothing, nothing)
        )
        test_parse(
            raw"\sum_{k=0}^n",
            (:underover,
                (:symbol, '∑', raw"\sum"),
                (:group, 'k', (:spaced, '='), '0'),
                'n'
            )
        )
        test_parse(
            raw"\lim_x",
            (:underover,
                (:function, "lim"),
                'x',
                nothing
            )
        )
    end

    @testset "Square root" begin
        test_parse(raw"\sqrt{x}", (:sqrt, 'x'))
        test_parse(
            raw"\sqrt{abc}",
            (:sqrt, (:group, 'a', 'b', 'c'))
        )
    end

    @testset "Spaced symbol" begin
        test_parse(raw"=", (:spaced, '='))
        test_parse(
            raw"\rightarrow",
            (:spaced, (:symbol, '→', raw"\rightarrow"))
        )
    end

    @testset "Subscript and superscript" begin
        @test texparse(raw"a^2_3") == texparse("a_3^2")
    end

    @testset "Symbol" begin
        for (char, sym) in zip(split("ϕ φ Φ"), split(raw"\phi \varphi \Phi"))
            test_parse(sym, (:symbol, first(char), sym))
            @test texparse(char) == texparse(sym)
        end

        # Check interaction with decoration
        test_parse(
            raw"ω_k",
            (:decorated, (:symbol, 'ω', "\\omega"), 'k', nothing)
        )
    end
end


text:
"""
    JSON3.@pretty json_str

Pretty print a JSON string or an object as JSON.
"""
macro pretty(json)
    return esc(:(JSON3.pretty($json)))
end

"""
    JSON3.pretty(x; kw...)
    JSON3.pretty(io, x; kw...)

Pretty print a JSON string.

## Args

* `x`: A JSON string, or an object to write to JSON then pretty print.
* `io`: The `IO` object to write the pretty printed string to. [default `stdout`]

## Keyword Args

See [`JSON3.write`](@ref) and [`JSON3.read`](@ref).
"""
pretty(str; kw...) = pretty(stdout, str; kw...)
pretty(out::IO, x; kw...) = pretty(out, JSON3.write(x; kw...); kw...)
function pretty(out::IO, str::String, indent=0, offset=0; kw...)
    buf = codeunits(str)
    len = length(buf)
    if len == 0
        return
    end
    pos = 1
    b = getbyte(buf, pos)
    @wh
    # printing object?
    if b == UInt8('{')
        Base.write(out, "{\n")

        obj = JSON3.read(str; kw...)

        if length(obj) == 0
            Base.write(out, "}")
            return
        end

        ks = collect(keys(obj))
        maxlen = maximum(map(sizeof, ks)) + 5
        indent += 1

        i = 1
        for (key, value) in obj
            Base.write(out, "  "^indent)
            Base.write(out, lpad("\"$(key)\"" * ": ", maxlen + offset, ' '))
            pretty(out, JSON3.write(value; kw...), indent, maxlen + offset; kw...)
            if i == length(obj)
                indent -= 1
                Base.write(out, "\n" * ("  "^indent * " "^offset) * "}")
            else
                Base.write(out, ",\n")
                i += 1
            end
        end
    # printing array?
    elseif b == UInt8('[')
        Base.write(out, "[\n")

        arr = JSON3.read(str; kw...)

        if length(arr) == 0
            Base.write(out, "]")
            return
        end

        indent += 1

        for (i, val) in enumerate(arr)
            Base.write(out, "  "^indent * " "^offset)
            pretty(out, JSON3.write(val; kw...), indent, offset; kw...)
            if i == length(arr)
                indent -= 1
                Base.write(out, "\n" * ("  "^indent * " "^offset) * "]")
            else
                Base.write(out, ",\n")
            end
        end
    # printing constant?
    else
        Base.write(out, str)
    end
    return
@label invalid
    Base.error("error pretty-fying json: $error")
end


text:
myParse(x) = parse.(Int, x)
caveMap = Matrix(hcat(myParse.(collect.(readlines("input.txt")))...)')

using DataStructures

function lowestRisk(start, goal, grid)

    # Utility Functions

    nRows, nCols = size(grid)

    function coordfromIndex(index)
        row = ((index - 1) % nRows) + 1
        col = ((index - 1) ÷ nRows) + 1
        return row, col
    end

    function indexfromCoord(row, col)
        return nRows * (col - 1) + row
    end

    function neighbors(index)
        row, col = coordfromIndex(index)
        arr = [(row - 1, col), (row, col - 1), (row + 1, col), (row, col + 1)]
        return [indexfromCoord(x, y) for (x, y) in arr if x in 1:nRows && y in 1:nCols]
    end

    # End: Utility Functions

    function heuristic(indexA, indexB)

        # Manhattan Distance

        A = coordfromIndex(indexA)
        B = coordfromIndex(indexB)

        return abs(A[2] - B[2]) + abs(A[1] - B[1])
    end

    frontier = PriorityQueue()
    enqueue!(frontier, start => 0)

    prev = Dict()
    costTill = Dict()

    prev[start] = nothing
    costTill[start] = 0

    while !isempty(frontier)
        current = dequeue!(frontier)

        # Early exit

        if current == goal
            break
        end

        for neighbor in neighbors(current)
            newCost = costTill[current] + grid[neighbor]

            # Update pathway if no costs or better cost

            if !haskey(costTill, neighbor) || newCost < costTill[neighbor]
                costTill[neighbor] = newCost
                priority = newCost + heuristic(neighbor, goal)
                frontier[neighbor] = priority
                prev[neighbor] = current
            end
        end
    end

    # Calculate risk
    cost = 0
    current = goal

    while current != start
        cost += grid[current]
        current = prev[current]
    end
    return cost
end

#! Part 1

lowestRisk(1, length(caveMap), caveMap)

#! Part 2

function addVal(x, val)
    sum = x + val
    x = sum != 9 ? sum % 9 : 9
end

newMap = vcat([hcat([addVal.(caveMap, Ref(i)) for i in j:j+4]...) for j in 0:4]...)

lowestRisk(1, length(newMap), newMap)

max_stars_repo_path: Day 15/main.jl

max_issues_repo_path: Day 15/main.jl

max_forks_repo_path: Day 15/main.jl


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 794a5356-6e05-44b4-ad9d-375a1fa3613b
begin
	using PlutoUI
	using DataFrames
	using BenchmarkTools
	using Statistics
	using Test
end

# ╔═╡ b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
md"""
## Introduction DataFrames / 5 - Manipulating columns

ref. [Introduction to DataFrames, Bogumił Kamiński](https://github.com/bkamins/Julia-DataFrames-Tutorial/blob/master/05_columns.ipynb)

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ d6993a5b-8158-4ce3-accb-8e8c514229b3
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ ce0a87f7-8f0c-4c87-8fd1-8198edbec86c
md"""
### Renaming columns

Let us start with a DataFrame of Bools that has default column names.
"""

# ╔═╡ e7020e48-1e43-4f07-93a5-07ea1f1a63e7
x₀ = DataFrame(rand(Bool, 3, 4), :auto)

# ╔═╡ 9528b5ce-da9e-449f-a042-bb9033f8bffa
begin
	using Random
	Random.seed!(42)

	x₀[:, shuffle(names(x₀))]
end

# ╔═╡ 9a500d40-75c6-4914-86af-2c255255c1ad
md"""
With `rename`, we create new DataFrame; here we rename the column `:x1` to `:A`. (`rename` also accepts collections of Pairs.)
"""

# ╔═╡ 33bdbcf6-893e-46ce-9370-917b87045b72
rename(x₀, :x1 => :A)

# ╔═╡ 98379291-e308-4642-a71e-b8cb5cf6f376
md"""
With `rename!` we do an in place transformation.

This time we have applied a function to every column name (note that the function gets a column names as a string).
"""

# ╔═╡ 2f70c7eb-ab73-4447-9529-354995484739
rename!(c -> c^2, x₀)

# ╔═╡ 549aa400-44b0-4bfb-8cb4-81bffdb6cd41
md"""
We can also change the name of a particular column without knowing the original.

Here we change the name of the third column, creating a new DataFrame.
"""

# ╔═╡ 31edf126-68b7-435f-b286-0e762e8affaf
rename(x₀, 3 => :third)

# ╔═╡ bac788c3-fa41-4daa-a34c-4d5e3fd8c285
md"""
If we pass a vector of names to `rename!`, we can change the names of all variables.
"""

# ╔═╡ 3925258e-c5f2-42c0-9af2-8b60cf3d51d2
rename!(x₀, [:a, :b, :c, :d])

# ╔═╡ aedfb361-edc8-46a3-9497-75542db83642
md"""
In all the above examples you could have used strings instead of symbols, e.g.
"""

# ╔═╡ 45f24ad5-ff3e-48b6-8917-36ea094388f2
rename!(x₀, string.('a':'d'))

# ╔═╡ 6b133d04-6912-4c5b-9784-bf7eea64bc55
md"""
`rename!` allows for circular renaming of columns, e.g.:
"""

# ╔═╡ ae7e4869-a3cc-4012-8436-47d3af7d28c0
rename!(x₀, "a"=>"d", "d"=>"a")

# ╔═╡ 96958f4d-df91-44f1-89fd-088fb0c6682d
md"""
We get an error when we try to provide duplicate names
"""

# ╔═╡ bd4ece67-9755-448c-93a8-30b391178c0d
@test_throws ArgumentError rename(x₀, fill(:a, 4))

# ╔═╡ 5b7b3282-a516-4b21-980b-ef552e6b0c0f
rename(x₀, fill(:a, 4), makeunique=true)

# ╔═╡ 4f72dc4a-31b8-40a4-81d9-a727ed2c22b8
md"""
### Reordering columns

We can reorder the names(x) vector as needed, creating a new DataFrame.
"""

# ╔═╡ 41bac74c-b6f9-4896-9b8d-fcb4d6ce7d47
md"""
Also `select!` can be used to achieve this in place (or select to perform a copy):
"""

# ╔═╡ b4e3e3cb-27a6-4002-993b-b5a87587db06
select!(x₀, 4:-1:1); x₀

# ╔═╡ 15106b6a-4f96-4701-bf95-4f0ed2d988e7
md"""
### Merging/Adding columns
"""

# ╔═╡ 44b6e6a2-b22a-485e-b161-6a79513486e5
x₁ = DataFrame([(i, j) for i ∈ 1:3, j ∈ 1:4], :auto)

# ╔═╡ dd2c281c-34a6-4aa1-a339-6d2cc5baa455
md"""
With `hcat` we can merge two DataFrames. Also `[x y]` syntax is supported but only when DataFrames have unique column names.
"""

# ╔═╡ 97acfc11-79a5-48eb-b8e8-5b040bd92b8f
hcat(x₁, x₁, makeunique=true)

# ╔═╡ b1be6837-38dd-4313-8860-e73ac130b30d
md"""
We can also use `hcat` to add a new column; a default name `:x1` will be used for this column, so `makeunique=true` is needed in our case.
"""

# ╔═╡ 72bab76f-91bb-4f78-9dfb-c865ad271cfb
y₁ = hcat(x₁, [1, 2, 3], makeunique=true)

# ╔═╡ 4d2d3d2b-cf7d-4ca6-9da0-67d28f269db9
md"""
We can also prepend a vector with `hcat`.
"""

# ╔═╡ 5881f7e3-f6d9-4dfc-ba1d-dfa0fea2c07c
hcat([1, 2, 3], x₁, makeunique=true)

# ╔═╡ 5d103d4b-9ca8-4274-b556-acbce4aea876
md"""
Alternatively you could append a vector with the following syntax. This is a bit more verbose but cleaner.
"""

# ╔═╡ eb08d57b-69b8-4543-88ba-c56793210681
y₂ = [x₁ DataFrame(A=[1, 2, 3])]

# ╔═╡ a33dc13b-732c-4ecf-9295-9fa0f5bc417c
md"""
And here we do the same but add column `:A` to the front.
"""

# ╔═╡ e1d534a0-dc8e-4fd6-a35a-314985d6eba8
y₃ = [DataFrame(A=[1, 2, 3]) x₁]

# ╔═╡ 129c6825-9c01-400b-9318-74a3bde8c163
md"""
A column can also be added in the middle. Here a brute-force method is used and a new DataFrame is created.
"""

# ╔═╡ 33bb1d78-0bcd-4c04-8ef8-d9e62c9f16c6
with_terminal() do
	@btime [$x₁[!, 1:2] DataFrame(A=[1,2,3]) $x₁[!, 3:4]]
end

# ╔═╡ 483bfe0d-6e54-4355-8688-d34d9b68fad0
md"""
We could also do this with a specialized in place method `insertcols!`. Let's add `:newcol` to the DataFrame y₂.
"""

# ╔═╡ e14dfcec-5710-4051-895e-b2f0fdcd6c78
insertcols!(y₂, 2, "newcol" => [1, 2, 3])

# ╔═╡ c968b160-53cb-4ecc-ad11-78e4b09a73b8
md"""
If you want to insert the same column name several times `makeunique=true` is needed as usual.
"""

# ╔═╡ 016e5a00-49d6-4ce0-a6be-62287d5c928b
insertcols!(y₂, 2, :newcol => [1, 2, 3], makeunique=true)

# ╔═╡ f312fcc0-e8f4-4e1a-9a7a-b2f6675caac5
md"""
We can see how much faster it is to insert a column with `insertcols!` than with `hcat` using `@btime` (note that here we use a Pair notation as an example).
"""

# ╔═╡ c64a8915-c2e8-47fe-96e3-2ab7a066764f
with_terminal() do
	@btime insertcols!(copy($x₁), 3, :A => [1, 2, 3])
end

# ╔═╡ 52d33d1a-0e2c-4f32-a1fd-2380b656460d
md"""
Let's use `insertcols!` to append a column in place (note that we dropped the index at which we insert the column).
"""

# ╔═╡ 33dfd57c-9252-495d-b39b-caf96ab8008f
insertcols!(x₁, :A => 1:3)

# ╔═╡ 520123eb-9673-4c8f-af83-f6f36a621602
md"""
and to in place `prepend` a column.
"""

# ╔═╡ d8d90255-2393-4362-a239-798976e50a92
insertcols!(x₁, 1, :B => 1:3)

# ╔═╡ 76498ebd-e783-4716-8c9a-f0ef05616a8f
md"""
Note that `insertcols!` can be used to insert several columns to a dataframe at once and that it performs broadcasting if needed:
"""

# ╔═╡ 91f780ee-6e10-4124-9d80-d8a0164155fe
df₀ = DataFrame(a = [1, 2, 3])

# ╔═╡ d58ac235-cf27-4540-aaab-413b72d07418
insertcols!(df₀, :b => "x", :c => 'a':'c', :d => Ref([1, 2, 3]))

# ╔═╡ 4a8e2e67-8e97-4e09-bffc-8866ff566648
md"""
Interestingly we can emulate hcat mutating the data frame in-place using `insertcols!`:
"""

# ╔═╡ b693f3ac-bcf8-43fd-9060-753267d7bef9
df₁ = DataFrame(a=[1,2])

# ╔═╡ 60cf6f09-55b8-4a70-a208-9e13f43d64ce
df₂ = DataFrame(b=[2,3], c=[3,4])

# ╔═╡ 686cf84f-8bc1-4a74-90f2-cfc52d3963cf
hcat(df₁, df₂)  ## df₁ is not touched

# ╔═╡ 3510270e-7edc-4edc-bb9b-e1b274479300
insertcols!(df₁, pairs(eachcol(df₂))...)

# ╔═╡ ddce5237-d203-43cf-879c-a6dd91adbf55
df₁   ## now we have changed df1

# ╔═╡ 51011425-7eff-4742-9753-d67405de81d5
md"""
### Subsetting/Removing columns

Let us create a new DataFrame `sdf` and show a few ways to create DataFrames with a subset of `x`'s columns.
"""

# ╔═╡ 8a99674d-1e8a-4991-9390-9dad31b3d998
sdf = DataFrame([(i,j) for i in 1:3, j in 1:5], :auto)

# ╔═╡ fa280e7c-cdeb-484d-9678-6d1359ebe5f2
md"""
First we could do this by index:
"""

# ╔═╡ b3a8f83d-8150-456e-9682-2195a2ea4aea
sdf[:, [1, 2, 4, 5]]   ## use ! instead of : for non-copying operation

# ╔═╡ c8f0a8e0-e433-4607-af1a-dfe54df641d4
md"""
or by column name:
"""

# ╔═╡ 07ae5e0b-df47-4d7e-871d-a0ea8b4bb839
sdf[:, [:x1, :x4]]

# ╔═╡ cbfe95e6-dc80-4b62-9f68-6fd9d72196e9
md"""
We can also choose to keep or exclude columns by Bool (we need a vector whose length is the number of columns in the original DataFrame).
"""

# ╔═╡ 08c6994f-52c2-43be-b547-8b10dbfb0bcd
sdf[:, [true, false, true, false, true]]

# ╔═╡ ab2465b5-e837-44a6-983f-92c7489abb05
md"""
Here we create a single column DataFrame
"""

# ╔═╡ 198df989-3013-429f-94d5-a5af95e1e8fd
sdf[:, [:x1]]

# ╔═╡ 592d59cf-a253-4103-9610-3a7c232d9bbd
md"""
and here we access the vector contained in column `:x1`
"""

# ╔═╡ 8659fa70-cc8f-4a15-bd5d-d9b0cc13c52c
sdf[!, :x1]   ## use : instead of ! to copy

# ╔═╡ a525f4b8-eb1c-4e75-9829-1fe31c2dc620
sdf.x1        ## the same

# ╔═╡ b346fbd0-a135-4200-9f9d-756afe0221ce
md"""
We could grab the same vector by column number
"""

# ╔═╡ b934b67d-8d8e-4854-a038-3e1cf71bf15b
sdf[!, 1]

# ╔═╡ 18487733-491e-4a95-a44d-f5ef109695c7
md"""
**Note** that getting a single column returns it without copying while creating a new DataFrame performs a copy of the column:
"""

# ╔═╡ 21fe70fd-0799-43e1-92ce-340c7c523c3e
@assert sdf[!, 1] !== sdf[!, [1]]

# ╔═╡ 04abd76d-5b51-45e4-8e23-d9b568387fba
md"""
We can also use `Regex`, `All`, `Between` and `Not` from `InvertedIndies.jl` for column selection:
"""

# ╔═╡ 96f3cd04-f903-4f1c-be07-94860d381bd5
sdf[!, r"[12]"]

# ╔═╡ 6dd8de75-fabf-4659-aea5-33b40e219463
sdf[!, Not(1)]

# ╔═╡ 00f4efee-e0da-42a1-8f0a-b6d7a3f4a035
sdf[!, Between(:x2, :x4)]

# ╔═╡ 1fc27c5d-092b-4e54-9eca-6582bc00a241
sdf[!, Cols(:x1, Between(:x3, :x5))]

# ╔═╡ 303882ce-bff9-4295-80b4-873f96de1345
select(sdf, :x1, Between(:x3, :x5), copycols=false)  ## the same as above

# ╔═╡ 6b42aeef-6eda-47c6-9aa9-4d06a9738bce
md"""
We can use `select` and `select!` functions to select a subset of columns from a data frame. `select` creates a new data frame and `select!` operates in place
"""

# ╔═╡ f70a4b29-5ec6-44df-b231-7a429de639d0
ndf = copy(sdf)

# ╔═╡ 3b379422-e957-47de-96f8-62101bf391a8
ndf₂ = select(ndf, [1, 2])

# ╔═╡ 8c7e1e4e-98de-4f50-9999-234c0c7ef309
select(ndf, Not([1, 2]))

# ╔═╡ e41dfcc4-c92c-4548-985e-8cb8eb32c858
md"""
by default `select` copies columns
"""

# ╔═╡ a5acfa9e-6af8-47b1-b85a-95dab43b7fa8
@assert ndf₂[!, 1] !== ndf[!, 1]

# ╔═╡ 74ccb3d6-9ecf-4118-b913-afe49849bd11
md"""
This can be avoided by using `copycols=false` keyword argument
"""

# ╔═╡ c7dc731b-3e97-4439-805f-dab44d1f43de
ndf₃ = select(ndf, [1, 2], copycols=false)

# ╔═╡ 452a8e80-b98a-40df-b942-a9359e3546a1
@assert ndf₃[!, 1] === ndf[!, 1]

# ╔═╡ 98789796-33e4-415f-a7a6-76a55078c92f
md"""
Using `select!` will modify the source dataframe
"""

# ╔═╡ 7e51ffcf-86b6-4832-9246-57e69a8b1d85
select!(ndf, [1, 2])

# ╔═╡ 5b4e71a2-bbb4-42dc-9674-2823e0aafeef
@assert ndf == ndf₂

# ╔═╡ 795c4656-e789-48de-a38b-9b5c70c2b86f
md"""
Here we create a copy of `sdf` and delete the 3rd column from the copy with `select!` and `Not`.
"""

# ╔═╡ ed1fee04-a51b-4909-9dd1-972a917335df
begin
	dfz = copy(sdf)
	select!(dfz, Not(3))
end

# ╔═╡ 386f5384-fb9a-48c9-8a91-c8a812ab3e1a
md"""
Alternatively we can achieve the same by using the `select` function
"""

# ╔═╡ 108b3bdb-8275-4e6b-8b30-6c6a8f81f950
select(sdf, Not(3))

# ╔═╡ 6de67316-b33d-456b-8754-88cc8e43819e
sdf   ## sdf stays unchanged

# ╔═╡ 998d6542-ff01-43f9-aa21-9f48f0cefdd6
md"""
#### Views

Note, that we can also create a view of a DataFrame when we want a subset of its columns:
"""

# ╔═╡ a84445cf-1598-4558-a978-ca1c912e1303
with_terminal() do
	@btime sdf[:, [1,3,5]]
end

# ╔═╡ 2ec79047-6e56-4303-ac57-dce06620ad45
with_terminal() do
	@btime @view sdf[:, [1, 3, 5]]
end

# ╔═╡ ecb25704-6a75-4bd2-a598-145ff9982dc4
md"""
### Modify column by name
"""

# ╔═╡ 2cbc8d22-6e0d-4a39-8669-7e3ab5f11521
xdf = DataFrame([(i,j) for i in 1:3, j in 1:5], :auto)

# ╔═╡ e6386760-0bf2-4dcc-a4ae-849dc89cc220
md"""
With the following syntax, the existing column is *modified without performing any copying* (**this is discouraged as it creates column alias**).
"""

# ╔═╡ 6b10c9af-934e-4f37-8013-120b918626d0

# ╔═╡ ef253e41-17b4-436b-9ecf-fb425814ca35
md"""
This syntax is safer
"""

# ╔═╡ e6072582-80d7-402c-9869-55441f9bafcf
xdf[!, :x1] = xdf[:, :x2]; xdf

# ╔═╡ 74051d06-5b6a-4fec-b76c-db71fe767110
md"""
We can also use the following syntax to add a new column at the end of a DataFrame.
"""

# ╔═╡ 71a75892-2724-4ee1-b20f-02030959b49b
xdf[!, :A] = [1, 2, 3]; xdf

# ╔═╡ e2ebb352-df12-4abc-9c06-fde8175d5bce
md"""
A new column name will be added to our DataFrame with the following syntax as well:
"""

# ╔═╡ 2c436f30-26cd-496d-8f50-a5aecd4ea419
xdf.B = 11:13; xdf

# ╔═╡ 6950fc3b-6436-4894-a7ee-cdf299ae39dc
md"""
### Find column name
"""

# ╔═╡ e1c027a2-2547-4612-ad80-d6b815977ebb
dfx = DataFrame([(i,j) for i in 1:3, j in 1:5], :auto)

# ╔═╡ 3c4f5e66-4903-4a2d-a554-3ddec91d3bd2
md"""
We can check if a column with a given name exists via `hasproperty`:
"""

# ╔═╡ 9cf7e784-0b0d-46d9-b32e-23eee66d1451
@assert hasproperty(dfx, :x1)

# ╔═╡ 3cf08ee6-9555-4ed7-8356-985633984030
md"""
and determine its index via `columnindex`:
"""

# ╔═╡ fcc0b1b4-6e99-40a0-97ae-d1d451946510
columnindex(dfx, :x2)

# ╔═╡ a9a58863-587a-4ee8-91af-438b5049d835
md"""
### Advanced ways of column selection

These are most useful for non-standard column names (e.g. containing spaces)
"""

# ╔═╡ 9c2d17a2-3758-4ff9-a099-d8c50ed9289c
begin
	df = DataFrame()
	df.x1 = 1:3
	df[!, "column 2"] = 4:6
	df
end

# ╔═╡ b8125aba-a8da-436e-b658-de277f992836
df."column 2"

# ╔═╡ bdffb7fd-bcd0-463a-b7a2-f9dc066d191e
df[:, "column 2"]

# ╔═╡ 48d6dcea-c814-4568-a47f-9d92fcb51b3f
md"""
Or we can interpolate column name using `:()` syntax
"""

# ╔═╡ 0a6b0a70-b6e3-4b06-8b52-dd5f9714bacc
with_terminal() do
	for n in names(df)
   		println(n, "\n", df.:($n), "\n")
	end
end

# ╔═╡ 9de09d8d-3e55-402d-8907-1131a2287ae3
md"""
#### Working on a collection of columns

When using `eachcol` of a data frame the resulting object retains reference to its parent and e.g. can be queried with `getproperty`:
"""

# ╔═╡ f8a050a2-dfc6-4473-a54d-045435f2dd0a
dfₓ = DataFrame(reshape(1:12, 3, 4), :auto)

# ╔═╡ baed2eb3-925f-485f-ae18-bc28e413bbe2
ec_df = eachcol(dfₓ)

# ╔═╡ 55f5c1bc-1c05-4b87-8673-aa6156ae2ccf
@assert ec_df.x1 == ec_df[1]

# ╔═╡ bbff0c85-b56f-4168-a035-84816b87b38c
md"""
### Transforming columns

We will get to this subject later in `10_transforms.jl` notebook, but here let us just note that `select`, `select!`, `transform`, `transform!` and `combine` functions allow to generate new columns based on the old columns of a data frame.

The general rules are the following:

  - `select` and `transform` always return the number of rows equal to the source dataframe, while `combine` returns any number of rows (`combine` is allowed to combine rows of the source dataframe)
  - `transform` retains columns from the old data frame
  - `select!` and `transform!` are in-place versions of `select` and `transform`
"""

# ╔═╡ 248fb9c3-caed-4298-bdb2-3d8426fa8514
tdf = DataFrame(reshape(1:12, 3, 4), :auto)

# ╔═╡ b0acc38f-49fe-405e-beb0-3932c59e0680
md"""
Here we add a new column `:res` that is a sum of columns `:x1` and `:x2`. A general syntax of transformations of this kind is:

	source_columns => function_to_apply => target_column_name

then `function_to_apply` gets columns selected by `source_columns` as positional arguments.
"""

# ╔═╡ afc96124-55d6-4b42-8ea1-de5d79647da8
transform(tdf, [:x1, :x2] => (+) => :res)

# ╔═╡ a866a4c5-c8a3-410f-bd26-b240859f28bb
md"""
One can omit passing `target_column_name` in which case it is automatically generated:
"""

# ╔═╡ dedc8234-630c-4fa2-8353-6b260620c3ea
combine(tdf, [:x1, :x2] => cor)

# ╔═╡ 03e05f07-3221-4cf7-a5c2-b81b71e0b4b3
md"""
Note that `combine` allowed the number of columns in the resulting dataframe to be changed. If we used `select` instead it would automatically broadcast the return value to match the number of rows of the source:
"""

# ╔═╡ 62d7860d-d52f-401c-b446-dd63638fc628
select(tdf, [:x1, :x2] => cor)

# ╔═╡ d4cbc501-7d7f-49ac-a296-971e0ada0edd
md"""
If we want to apply some function on each row of the source wrap it in `ByRow`:
"""

# ╔═╡ 28b31f33-b76d-4910-ba04-2d9d8f01ffd6
select(tdf, :x1, :x2, [:x1, :x2] => ByRow(string))

# ╔═╡ ef51c633-f8be-4884-ad62-c7969cb94bf8
md"""
Also if we want columns to be passed as a NamedTuple to a funcion (instead of being positional arguments) wrap them in `AsTable`:
"""

# ╔═╡ 35426049-4cc9-4fdb-b856-9e73fa3d7343
select(tdf, :x1, :x2, AsTable([:x1, :x2]) => x -> x.x1 + x.x2)

# ╔═╡ e26f822a-1649-43e1-b98e-93ed3d8d9a3b
md"""
For simplicity (as this functionality is often needed) there is a special treatement of nrow function that can be just passed as a transformation (without specifying of column selector):
"""

# ╔═╡ 7c6ca147-007d-46db-8dff-51a0e97dfd57
select(tdf, :x1, nrow => "number_of_rows")

# ╔═╡ 2166ab63-1628-4726-8d01-b6fbf25f966b
md"""
*Note* that in `select` the number of rows is automatically broadcasted to match the number of rows of the source data frame.

Finally we can conveniently create multiple columns with one function, e.g.:
"""

# ╔═╡ 4d86e251-c8f7-4281-a9ce-ea17d016f871
select(tdf, :x1, :x1 => ByRow(x -> [x ^ 2, x ^ 3]) => ["x1²", "x1³"])

# ╔═╡ f78fc90b-b3a5-4be7-90c2-25dc3dfdbe1c
md"""
or with the following which produces the same result (as above):
"""

# ╔═╡ f8ee1b70-5588-4328-a036-cf3499bdb1d0
select(tdf, :x1, :x1 => (x -> DataFrame("x1²" => x .^ 2, "x1³" => x .^ 3)) => AsTable)

# ╔═╡ c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e
html"""

"""

# ╔═╡ Cell order:
# ╟─b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
# ╠═794a5356-6e05-44b4-ad9d-375a1fa3613b
# ╟─d6993a5b-8158-4ce3-accb-8e8c514229b3
# ╟─ce0a87f7-8f0c-4c87-8fd1-8198edbec86c
# ╠═e7020e48-1e43-4f07-93a5-07ea1f1a63e7
# ╟─9a500d40-75c6-4914-86af-2c255255c1ad
# ╠═33bdbcf6-893e-46ce-9370-917b87045b72
# ╟─98379291-e308-4642-a71e-b8cb5cf6f376
# ╠═2f70c7eb-ab73-4447-9529-354995484739
# ╟─549aa400-44b0-4bfb-8cb4-81bffdb6cd41
# ╠═31edf126-68b7-435f-b286-0e762e8affaf
# ╟─bac788c3-fa41-4daa-a34c-4d5e3fd8c285
# ╠═3925258e-c5f2-42c0-9af2-8b60cf3d51d2
# ╟─aedfb361-edc8-46a3-9497-75542db83642
# ╠═45f24ad5-ff3e-48b6-8917-36ea094388f2
# ╟─6b133d04-6912-4c5b-9784-bf7eea64bc55
# ╠═ae7e4869-a3cc-4012-8436-47d3af7d28c0
# ╟─96958f4d-df91-44f1-89fd-088fb0c6682d
# ╠═bd4ece67-9755-448c-93a8-30b391178c0d
# ╠═5b7b3282-a516-4b21-980b-ef552e6b0c0f
# ╟─4f72dc4a-31b8-40a4-81d9-a727ed2c22b8
# ╠═9528b5ce-da9e-449f-a042-bb9033f8bffa
# ╟─41bac74c-b6f9-4896-9b8d-fcb4d6ce7d47
# ╠═b4e3e3cb-27a6-4002-993b-b5a87587db06
# ╟─15106b6a-4f96-4701-bf95-4f0ed2d988e7
# ╠═44b6e6a2-b22a-485e-b161-6a79513486e5
# ╟─dd2c281c-34a6-4aa1-a339-6d2cc5baa455
# ╠═97acfc11-79a5-48eb-b8e8-5b040bd92b8f
# ╟─b1be6837-38dd-4313-8860-e73ac130b30d
# ╠═72bab76f-91bb-4f78-9dfb-c865ad271cfb
# ╟─4d2d3d2b-cf7d-4ca6-9da0-67d28f269db9
# ╠═5881f7e3-f6d9-4dfc-ba1d-dfa0fea2c07c
# ╟─5d103d4b-9ca8-4274-b556-acbce4aea876
# ╠═eb08d57b-69b8-4543-88ba-c56793210681
# ╟─a33dc13b-732c-4ecf-9295-9fa0f5bc417c
# ╠═e1d534a0-dc8e-4fd6-a35a-314985d6eba8
# ╟─129c6825-9c01-400b-9318-74a3bde8c163
# ╠═33bb1d78-0bcd-4c04-8ef8-d9e62c9f16c6
# ╟─483bfe0d-6e54-4355-8688-d34d9b68fad0
# ╠═e14dfcec-5710-4051-895e-b2f0fdcd6c78
# ╟─c968b160-53cb-4ecc-ad11-78e4b09a73b8
# ╠═016e5a00-49d6-4ce0-a6be-62287d5c928b
# ╟─f312fcc0-e8f4-4e1a-9a7a-b2f6675caac5
# ╠═c64a8915-c2e8-47fe-96e3-2ab7a066764f
# ╟─52d33d1a-0e2c-4f32-a1fd-2380b656460d
# ╠═33dfd57c-9252-495d-b39b-caf96ab8008f
# ╟─520123eb-9673-4c8f-af83-f6f36a621602
# ╠═d8d90255-2393-4362-a239-798976e50a92
# ╟─76498ebd-e783-4716-8c9a-f0ef05616a8f
# ╠═91f780ee-6e10-4124-9d80-d8a0164155fe
# ╠═d58ac235-cf27-4540-aaab-413b72d07418
# ╟─4a8e2e67-8e97-4e09-bffc-8866ff566648
# ╠═b693f3ac-bcf8-43fd-9060-753267d7bef9
# ╠═60cf6f09-55b8-4a70-a208-9e13f43d64ce
# ╠═686cf84f-8bc1-4a74-90f2-cfc52d3963cf
# ╠═3510270e-7edc-4edc-bb9b-e1b274479300
# ╠═ddce5237-d203-43cf-879c-a6dd91adbf55
# ╟─51011425-7eff-4742-9753-d67405de81d5
# ╠═8a99674d-1e8a-4991-9390-9dad31b3d998
# ╟─fa280e7c-cdeb-484d-9678-6d1359ebe5f2
# ╠═b3a8f83d-8150-456e-9682-2195a2ea4aea
# ╟─c8f0a8e0-e433-4607-af1a-dfe54df641d4
# ╠═07ae5e0b-df47-4d7e-871d-a0ea8b4bb839
# ╟─cbfe95e6-dc80-4b62-9f68-6fd9d72196e9
# ╠═08c6994f-52c2-43be-b547-8b10dbfb0bcd
# ╟─ab2465b5-e837-44a6-983f-92c7489abb05
# ╠═198df989-3013-429f-94d5-a5af95e1e8fd
# ╟─592d59cf-a253-4103-9610-3a7c232d9bbd
# ╠═8659fa70-cc8f-4a15-bd5d-d9b0cc13c52c
# ╠═a525f4b8-eb1c-4e75-9829-1fe31c2dc620
# ╟─b346fbd0-a135-4200-9f9d-756afe0221ce
# ╠═b934b67d-8d8e-4854-a038-3e1cf71bf15b
# ╟─18487733-491e-4a95-a44d-f5ef109695c7
# ╠═21fe70fd-0799-43e1-92ce-340c7c523c3e
# ╟─04abd76d-5b51-45e4-8e23-d9b568387fba
# ╠═96f3cd04-f903-4f1c-be07-94860d381bd5
# ╠═6dd8de75-fabf-4659-aea5-33b40e219463
# ╠═00f4efee-e0da-42a1-8f0a-b6d7a3f4a035
# ╠═1fc27c5d-092b-4e54-9eca-6582bc00a241
# ╠═303882ce-bff9-4295-80b4-873f96de1345
# ╟─6b42aeef-6eda-47c6-9aa9-4d06a9738bce
# ╠═f70a4b29-5ec6-44df-b231-7a429de639d0
# ╠═3b379422-e957-47de-96f8-62101bf391a8
# ╠═8c7e1e4e-98de-4f50-9999-234c0c7ef309
# ╟─e41dfcc4-c92c-4548-985e-8cb8eb32c858
# ╠═a5acfa9e-6af8-47b1-b85a-95dab43b7fa8
# ╟─74ccb3d6-9ecf-4118-b913-afe49849bd11
# ╠═c7dc731b-3e97-4439-805f-dab44d1f43de
# ╠═452a8e80-b98a-40df-b942-a9359e3546a1
# ╟─98789796-33e4-415f-a7a6-76a55078c92f
# ╠═7e51ffcf-86b6-4832-9246-57e69a8b1d85
# ╠═5b4e71a2-bbb4-42dc-9674-2823e0aafeef
# ╟─795c4656-e789-48de-a38b-9b5c70c2b86f
# ╠═ed1fee04-a51b-4909-9dd1-972a917335df
# ╟─386f5384-fb9a-48c9-8a91-c8a812ab3e1a
# ╠═108b3bdb-8275-4e6b-8b30-6c6a8f81f950
# ╠═6de67316-b33d-456b-8754-88cc8e43819e
# ╟─998d6542-ff01-43f9-aa21-9f48f0cefdd6
# ╠═a84445cf-1598-4558-a978-ca1c912e1303
# ╠═2ec79047-6e56-4303-ac57-dce06620ad45
# ╟─ecb25704-6a75-4bd2-a598-145ff9982dc4
# ╠═2cbc8d22-6e0d-4a39-8669-7e3ab5f11521
# ╟─e6386760-0bf2-4dcc-a4ae-849dc89cc220
# ╠═6b10c9af-934e-4f37-8013-120b918626d0
# ╟─ef253e41-17b4-436b-9ecf-fb425814ca35
# ╠═e6072582-80d7-402c-9869-55441f9bafcf
# ╟─74051d06-5b6a-4fec-b76c-db71fe767110
# ╠═71a75892-2724-4ee1-b20f-02030959b49b
# ╟─e2ebb352-df12-4abc-9c06-fde8175d5bce
# ╠═2c436f30-26cd-496d-8f50-a5aecd4ea419
# ╟─6950fc3b-6436-4894-a7ee-cdf299ae39dc
# ╠═e1c027a2-2547-4612-ad80-d6b815977ebb
# ╟─3c4f5e66-4903-4a2d-a554-3ddec91d3bd2
# ╠═9cf7e784-0b0d-46d9-b32e-23eee66d1451
# ╟─3cf08ee6-9555-4ed7-8356-985633984030
# ╠═fcc0b1b4-6e99-40a0-97ae-d1d451946510
# ╟─a9a58863-587a-4ee8-91af-438b5049d835
# ╠═9c2d17a2-3758-4ff9-a099-d8c50ed9289c
# ╠═b8125aba-a8da-436e-b658-de277f992836
# ╠═bdffb7fd-bcd0-463a-b7a2-f9dc066d191e
# ╟─48d6dcea-c814-4568-a47f-9d92fcb51b3f
# ╠═0a6b0a70-b6e3-4b06-8b52-dd5f9714bacc
# ╟─9de09d8d-3e55-402d-8907-1131a2287ae3
# ╠═f8a050a2-dfc6-4473-a54d-045435f2dd0a
# ╠═baed2eb3-925f-485f-ae18-bc28e413bbe2
# ╠═55f5c1bc-1c05-4b87-8673-aa6156ae2ccf
# ╟─bbff0c85-b56f-4168-a035-84816b87b38c
# ╠═248fb9c3-caed-4298-bdb2-3d8426fa8514
# ╟─b0acc38f-49fe-405e-beb0-3932c59e0680
# ╠═afc96124-55d6-4b42-8ea1-de5d79647da8
# ╟─a866a4c5-c8a3-410f-bd26-b240859f28bb
# ╠═dedc8234-630c-4fa2-8353-6b260620c3ea
# ╟─03e05f07-3221-4cf7-a5c2-b81b71e0b4b3
# ╠═62d7860d-d52f-401c-b446-dd63638fc628
# ╟─d4cbc501-7d7f-49ac-a296-971e0ada0edd
# ╠═28b31f33-b76d-4910-ba04-2d9d8f01ffd6
# ╟─ef51c633-f8be-4884-ad62-c7969cb94bf8
# ╠═35426049-4cc9-4fdb-b856-9e73fa3d7343
# ╟─e26f822a-1649-43e1-b98e-93ed3d8d9a3b
# ╠═7c6ca147-007d-46db-8dff-51a0e97dfd57
# ╟─2166ab63-1628-4726-8d01-b6fbf25f966b
# ╠═4d86e251-c8f7-4281-a9ce-ea17d016f871
# ╟─f78fc90b-b3a5-4be7-90c2-25dc3dfdbe1c
# ╠═f8ee1b70-5588-4328-a036-cf3499bdb1d0
# ╟─c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
##############################################################################
##
## Non Linear solver using Pseudo-Transient Continuation Method
##
##############################################################################

# Implicit time step
function implicit_timestep(F!, ypost, Δ; is_algebraic = fill(false, size(ypost)...), iterations = 100, verbose = true, method = :newton, autodiff = :forward, maxdist = 1e-9, J0c = (nothing, nothing))
    F_helper!(ydot, y) = (F!(ydot, y) ; ydot .+= .!is_algebraic .* (ypost .- y) ./ Δ)
    J0, colorvec = J0c
    if J0 == nothing
        result = nlsolve(F_helper!, ypost; iterations = iterations, show_trace = verbose, ftol = maxdist, method = method, autodiff = autodiff)
    else
        if autodiff == :forward
            jac_cache = ForwardColorJacCache(F_helper!, deepcopy(ypost); colorvec = colorvec, sparsity = J0)
            j_helper! = (J, y) -> forwarddiff_color_jacobian!(J, F_helper!, y, jac_cache)
        else
            j_helper! = (J, y) -> DiffEqDiffTools.finite_difference_jacobian!(J, F_helper!, y; colorvec = colorvec)
        end
        result = nlsolve(OnceDifferentiable(F_helper!, j_helper!, deepcopy(ypost), deepcopy(ypost), J0), ypost; iterations = iterations, show_trace = verbose, ftol = maxdist, method = method)
    end
    return result.zero, result.residual_norm
end

# Solve for steady state
function finiteschemesolve(F!, y0; Δ = 1.0, is_algebraic = fill(false, size(y0)...), iterations = 100, inner_iterations = 10, verbose = true, inner_verbose = false, method = :newton, autodiff = :forward, maxdist = 1e-9, scale = 10.0, J0c = (nothing, nothing))
    ypost = y0
    ydot = zero(y0)
    F!(ydot, ypost)
    distance = norm(ydot) / length(ydot)
    isnan(distance) && throw("F! returns NaN with the initial value")
    if Δ == Inf
        ypost, distance = implicit_timestep(F!, y0, Δ; is_algebraic = is_algebraic, verbose = verbose, iterations = iterations,  method = method, autodiff = autodiff, maxdist = maxdist, J0c = J0c)
    else
        coef = 1.0
        olddistance = distance
        iter = 0
        while (iter < iterations) & (Δ >= 1e-12) & (distance > maxdist)
            iter += 1
            y, nldistance = implicit_timestep(F!, ypost, Δ; is_algebraic = is_algebraic, verbose = inner_verbose, iterations = inner_iterations, method = method, autodiff = autodiff, maxdist = maxdist, J0c = J0c)
            F!(ydot, y)
            distance, olddistance = norm(ydot) / length(ydot), distance
            distance = isnan(distance) ? Inf : distance
            if nldistance <= maxdist
                # if the implicit time step is correctly solved
                verbose && @show iter, Δ, distance
                coef = (distance <= olddistance) ? scale * coef : 1.0
                Δ = Δ * coef * olddistance / distance
                ypost, y = y, ypost
            else
                verbose && @show iter, Δ, NaN
                # if the implict time step is not solved
                # revert and diminish the time step
                coef = 1.0
                Δ = Δ / 10
                distance = olddistance
            end
        end
    end
    verbose && (distance > maxdist) && @warn "Iteration did not converge"
    return ypost, distance
end


text:
function GRGSFormatReader(arg0::JString, arg1::jboolean)
    return GRGSFormatReader((JString, jboolean), arg0, arg1)
end

function GRGSFormatReader(arg0::JString, arg1::jboolean, arg2::TimeScale)
    return GRGSFormatReader((JString, jboolean, TimeScale), arg0, arg1, arg2)
end

function get_provider(obj::GRGSFormatReader, arg0::jboolean, arg1::jint, arg2::jint)
    return jcall(obj, "getProvider", RawSphericalHarmonicsProvider, (jboolean, jint, jint), arg0, arg1, arg2)
end

function load_data(obj::GRGSFormatReader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end


text:
const DIR = dirname(@__FILE__)
const GENDIR = joinpath(DIR, "..")
const TEMPCFGFILE = joinpath(GENDIR, "config.json")
const SRCDIR = joinpath(GENDIR, "src")
const JULIA = joinpath(Sys.BINDIR, "julia")
const SWAGGERDIR=read(`$(JULIA) -e 'import Swagger; print(normpath(joinpath(dirname(pathof(Swagger)), "..")))'`, String)
const SWAGGERGEN = joinpath(SWAGGERDIR, "plugin", "generate.sh")

function gen(pkg, swg)
    open(TEMPCFGFILE, "w") do f
        println(f, """{ "packageName": "$pkg" } """)
    end
    outpath = joinpath(SRCDIR, pkg*"tmp")
    pkgpath = joinpath(SRCDIR, pkg)
    mkpath(outpath)
    run(`$SWAGGERGEN -i $swg -o $outpath -c $TEMPCFGFILE`)
    mv(joinpath(outpath, "src"), pkgpath; force=true)
    rm(outpath; force=true, recursive=true)
    rm(TEMPCFGFILE; force=true)
    nothing
end

const USAGE = """Usage: julia generate.jl <specs>
Specifications: https://github.com/firecracker-microvm/firecracker/blob/master/api_server/swagger/firecracker.yaml"""

if isempty(ARGS)
    println(USAGE)
else
    gen("FirecrackerApis", ARGS[1])
end


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: HairpinIncr.mso								$
#*------------------------------------------------------------------
type UpperPipe_basic
	UpperPipe_basic()=begin
		new(
			HairpinIncr_basic(),
			[
				:(_base_1.Details.Q(1:N) = _base_1.InletOuter.F*(_base_1.Outer.HeatTransfer.Enth(2:_base_1.Npoints) - _base_1.Outer.HeatTransfer.Enth(1:N))),
				:(_base_1.Details.Q(1:N) = -_base_1.InletInner.F*(_base_1.Inner.HeatTransfer.Enth(1:N) - _base_1.Inner.HeatTransfer.Enth(2:_base_1.Npoints))),
				:(_base_1.Details.Q = _base_1.Details.Ud*_base_1.Pi*_base_1.DoInner*(_base_1.Lpipe/N)*(_base_1.Outer.Properties.Average.T - _base_1.Inner.Properties.Average.T)),
				:(_base_1.Details.Q(1:N) = _base_1.InletInner.F*(_base_1.Inner.HeatTransfer.Enth(1:N)-_base_1.Inner.HeatTransfer.Enth(2:_base_1.Npoints))),
				:(_base_1.Details.Q(1:N) = -_base_1.InletOuter.F*(_base_1.Outer.HeatTransfer.Enth(2:_base_1.Npoints) - _base_1.Outer.HeatTransfer.Enth(1:N))),
				:(_base_1.Details.Q = _base_1.Details.Ud*_base_1.Pi*_base_1.DoInner*(_base_1.Lpipe/N)*(_base_1.Inner.Properties.Average.T - _base_1.Outer.Properties.Average.T)),
				:(_base_1.Inner.HeatTransfer.Enth(1) = _base_1.InletInner.h),
				:(_base_1.Inner.HeatTransfer.Enth(_base_1.Npoints) = _base_1.OutletInner.h),
				:(_base_1.Inner.HeatTransfer.Tlocal(1) = _base_1.InletInner.T),
				:(_base_1.Inner.HeatTransfer.Tlocal(_base_1.Npoints) = _base_1.OutletInner.T),
				:(_base_1.Inner.PressureDrop.Plocal(1) = _base_1.InletInner.P),
				:(_base_1.Inner.PressureDrop.Plocal(_base_1.Npoints) = _base_1.OutletInner.P),
				:(_base_1.Outer.HeatTransfer.Enth(_base_1.Npoints) = _base_1.InletOuter.h),
				:(_base_1.Outer.HeatTransfer.Enth(1) = _base_1.OutletOuter.h),
				:(_base_1.Outer.HeatTransfer.Tlocal(_base_1.Npoints) = _base_1.InletOuter.T),
				:(_base_1.Outer.HeatTransfer.Tlocal(1) = _base_1.OutletOuter.T),
				:(_base_1.Outer.PressureDrop.Plocal(_base_1.Npoints) = _base_1.InletOuter.P),
				:(_base_1.Outer.PressureDrop.Plocal(1) = _base_1.OutletOuter.P),
				:(_base_1.Outer.PressureDrop.Pd_fric([1:N]) = (2*_base_1.Outer.PressureDrop.fi([1:N])*_base_1.Lincr(1+_base_1.Npoints-[1:N])*_base_1.Outer.Properties.Average.rho([1:N])*_base_1.Outer.HeatTransfer.Vmean([1:N])^2)/(_base_1.Outer.PressureDrop.Dh*_base_1.Outer.HeatTransfer.Phi([1:N]))),
				:(_base_1.Outer.PressureDrop.Pd_fric(_base_1.Npoints) = 0*"kPa"),
				:(_base_1.Inner.PressureDrop.Pd_fric(2:_base_1.Npoints) = (2*_base_1.Inner.PressureDrop.fi*_base_1.Lincr(2:_base_1.Npoints)*_base_1.Inner.Properties.Average.rho*_base_1.Inner.HeatTransfer.Vmean^2)/(_base_1.DiInner*_base_1.Inner.HeatTransfer.Phi)),
				:(_base_1.Inner.PressureDrop.Pd_fric(1) = 0*"kPa"),
			],
			[
				"Energy Balance Outer Stream","Energy Balance Inner Stream","Incremental Duty","Energy Balance Hot Stream","Energy Balance Cold Stream","Incremental Duty","Enthalpy of Inner Side - Inlet Boundary","Enthalpy of inner Side - Outlet Boundary","Temperature of Inner Side - Inlet Boundary","Temperature of Inner Side - Outlet Boundary","Pressure of Inner Side - Inlet Boundary","Pressure of Inner Side - Outlet Boundary","Enthalpy of Outer Side - Inlet Boundary","Enthalpy of Outer Side - Outlet Boundary","Temperature of Outer Side - Inlet Boundary","Temperature of Outer Side - Outlet Boundary","Pressure of Outer Side - Inlet Boundary","Pressure of Outer Side - Outlet Boundary","Outer Pipe Pressure Drop for friction","Outer Pipe Pressure Drop for friction","Inner Pipe Pressure Drop for friction","Inner Pipe Pressure Drop for friction",
			],
		)
	end
	_base_1::HairpinIncr_basic
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	attributes::Dict{Symbol,Any}
end
export UpperPipe_basic
function setEquationFlow(in::UpperPipe_basic)
	let switch=HotSide
		if HotSide==InletInner.T > InletOuter.T
			set(switch,"inner")
		end
		if HotSide==InletInner.T < InletOuter.T
			set(switch,"outer")
		end
		if switch=="outer"
			addEquation(1)
			addEquation(2)
			addEquation(3)
		elseif switch=="inner"
			addEquation(4)
			addEquation(5)
			addEquation(6)
		end
	end
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	addEquation(22)
end
function atributes(in::UpperPipe_basic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Incremental Hairpin Heat Exchanger. "
	fields[:Info]="Incremental approach for Hairpin heat exchanger. "
	drive!(fields,_)
	return fields
end
UpperPipe_basic(_::Dict{Symbol,Any})=begin
	newModel=UpperPipe_basic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Documenter
import Pkg
using SummationByPartsOperators

# Define module-wide setups such that the respective modules are available in doctests
DocMeta.setdocmeta!(SummationByPartsOperators,
  :DocTestSetup, :(using SummationByPartsOperators); recursive=true)

# Copy some files from the top level directory to the docs and modify them
# as necessary
open(joinpath(@__DIR__, "src", "license.md"), "w") do io
  # Point to source license file
  println(io, """
  ```@meta
  EditURL = "https://github.com/ranocha/SummationByPartsOperators.jl/blob/main/LICENSE.md"
  ```
  """)
  # Write the modified contents
  println(io, "# License")
  println(io, "")
  for line in eachline(joinpath(dirname(@__DIR__), "LICENSE.md"))
    line = replace(line, "[LICENSE.md](LICENSE.md)" => "[License](@ref)")
    println(io, "> ", line)
  end
end

open(joinpath(@__DIR__, "src", "contributing.md"), "w") do io
  # Point to source license file
  println(io, """
  ```@meta
  EditURL = "https://github.com/ranocha/SummationByPartsOperators.jl/blob/main/CONTRIBUTING.md"
  ```
  """)
  # Write the modified contents
  println(io, "# Contributing")
  println(io, "")
  for line in eachline(joinpath(dirname(@__DIR__), "CONTRIBUTING.md"))
    line = replace(line, "[LICENSE.md](LICENSE.md)" => "[License](@ref)")
    println(io, "> ", line)
  end
end

# Make documentation
makedocs(
  modules = [SummationByPartsOperators],
  sitename="SummationByPartsOperators.jl",
  format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    canonical = "https://ranocha.github.io/SummationByPartsOperators.jl/stable"
  ),
  # Explicitly specify documentation structure
  pages = [
    "Home" => "index.md",
    "Introduction" => "introduction.md",
    "Tutorials" => [
      "tutorials/constant_linear_advection.md",
      "tutorials/advection_diffusion.md",
      "tutorials/variable_linear_advection.md",
      "tutorials/wave_equation.md",
    ],
    "Applications & references" => "applications.md",
    "Benchmarks" => "benchmarks.md",
    "API reference" => "api_reference.md",
    "Contributing" => "contributing.md",
    "License" => "license.md"
  ],
  strict = true # to make the GitHub action fail when doctests fail
)

deploydocs(
  repo = "github.com/ranocha/SummationByPartsOperators.jl",
  devbranch = "main",
  push_preview = true
)


text:
# # NNDescent on MNIST digits
# In this example, we build an approximate KNN graph over the MNIST digits and
# then use it as an index to find approximate nearest neighbors of new points.

using NearestNeighborDescent
using Distances
using BenchmarkTools
using MLDatasets
mnist_x = collect(eachcol(MNIST.convert2features(MNIST.traintensor(Float64, 1:10000))))
mnist_y = collect(eachcol(MNIST.convert2features(MNIST.testtensor(Float64, 1:1000))))

## Construct the KNN Graph
graph = nndescent(mnist_x, 20, Euclidean())

## Find neighbors of new points
inds, dists = search(graph, mnist_x, mnist_y, 20, Euclidean())


text:
# import Base.IO
using Iterators

const END_OF_DATA_SECTION = -1
const JULIA_EXCEPTION_THROWN = -2
const TIMING_DATA = -3
const END_OF_STREAM = -4
const NULL = -5

readint(io::IO) = ntoh(read(io, Int32))

"""
Read data object from a ioet. Returns code and byte array:
 * if code is negative, it's considered as a special command code
 * if code is positive, it's considered as array length
"""
function readobj(io::IO)
    len = readint(io)
    bytes = len > 0 ? readbytes(io, len) : []
    return len, bytes
end

writeint(io::IO, x::Int) = write(io, hton(Int32(x)))

"""Write length and byte array of that length to a ioet"""
function writeobj(io::IO, obj::Vector{UInt8})
    writeint(io, length(obj))
    write(io, obj)
end

writeobj(io::IO, obj) = writeobj(io, convert(Vector{UInt8}, obj))

function load_stream{T}(::Type{T}, io::IO)
    function it()
        code, _next = readobj(io)
        while code != END_OF_DATA_SECTION
            data = from_bytes(T, _next)
            produce(data)
            code, _next = readobj(io)
        end
    end
    return Task(it)
end

function dump_stream(io::IO, it)
    for v in it
        writeobj(io, to_bytes(v))
    end
end

function launch_worker()
    port = parse(Int, readline(STDIN))
    sock = connect("127.0.0.1", port)
    try
        split = readint(sock)
        # info("Julia: starting partition id: $split")
        T = deserialize(sock)
        cmd = readobj(sock)[2]
        func = deserialize(IOBuffer(cmd))
        it = load_stream(T, sock)
        dump_stream(sock, func(split, it))
        writeint(sock, END_OF_DATA_SECTION)
        writeint(sock, END_OF_STREAM)
        # info("Julia: exiting")
    catch e
        # TODO: handle the case when JVM closes connection
        Base.show_backtrace(STDERR, catch_backtrace())
        writeint(sock, JULIA_EXCEPTION_THROWN)
        rethrow()
    end
end


text:
using Test
using AgentsPlots
using Agents


text:
# version = v"1.14.1119"
version = v"1.14.1227"
deps_dir = dirname(@__FILE__)
par_dir  = abspath(joinpath(deps_dir, ".."))
sdk_dir  = joinpath(deps_dir, "sdk_v$version")
lib_target_dir = joinpath(par_dir, "lib")
lib_dir  = joinpath(sdk_dir, "lib")
lib_x64_dir = joinpath(lib_dir, "x64")
lib_x86_dir = joinpath(lib_dir, "x86")
lib_mac_dir = joinpath(lib_dir, "mac")

if Sys.islinux() || Sys.isapple()
    sdk_file = joinpath(deps_dir, "ASI_linux_mac_SDK_V$version.tar.bz2")
    if !isfile(sdk_file)
        println("Downloading ASI_linux_mac_SDK_V$version.tar.bz2 from astronomy-imaging-camera.com (ZWO)")
        download("https://astronomy-imaging-camera.com/software/ASI_linux_mac_SDK_V$version.tar.bz2", sdk_file)
    end
elseif Sys.iswindows()
    sdk_file = joinpath(deps_dir, "ASI_Windows_SDK_V$version.zip")
    if !isfile(sdk_file)
        println("Downloading ASI_Windows_SDK_V$version.zip from astronomy-imaging-camera.com (ZWO)")
        download("https://astronomy-imaging-camera.com/software/ASI_Windows_SDK_V$version.zip", sdk_file)
    end
end

if !isdir(sdk_dir) mkdir(sdk_dir) end
if !isdir(lib_target_dir) mkdir(lib_target_dir) end

#extract
if isfile(sdk_file)
    if Sys.isunix()
        unpack_cmd = `tar xjf $sdk_file --directory=$sdk_dir` end
    if Sys.iswindows()
        if isdefined(Base, :LIBEXECDIR)
          const exe7z = joinpath(Sys.BINDIR, Base.LIBEXECDIR, "7z.exe")
        else
          const exe7z = joinpath(Sys.BINDIR, "7z.exe")
        end
        unpack_cmd = `$exe7z x -o$sdk_dir -y $sdk_file`
    end
    run(unpack_cmd)
end

# copy extracted library file to lib subfolder
if Sys.islinux()
    mv(joinpath(lib_dir, "asi.rules"), joinpath(deps_dir, "..", "asi.rules"), force=true)

    if Sys.ARCH==:x86_64
        if Sys.WORD_SIZE==64
            source_file = joinpath(lib_x64_dir, "libASICamera2.so.$version")
            target_file = joinpath(lib_target_dir, "libASICamera2.so")
            mv(source_file, target_file, force=true)
        else
            source_file = joinpath(lib_x86_dir, "libASICamera2.so.$version")
            target_file = joinpath(lib_target_dir, "libASICamera2.so")
            mv(source_file, target_file, force=true)
        end

    # RaspberryPi support
    elseif Sys.ARCH==:aarch64
        source_file = joinpath(lib_dir, "armv8", "libASICamera2.so.$version")
        target_file = joinpath(lib_target_dir, "libASICamera2.so")
        mv(source_file, target_file, force=true)

    elseif Sys.ARCH==:arm || (Sys.ARCH==:aarch && Sys.WORD_SIZE==32)
        source_file = joinpath(lib_dir, "armv7", "libASICamera2.so.$version")
        target_file = joinpath(lib_target_dir, "libASICamera2.so")
        mv(source_file, target_file, force=true)
    end

elseif Sys.isapple()
    source_file = joinpath(lib_mac_dir, "libASICamera2.dylib.$version")
    target_file = joinpath(lib_target_dir, "libASICamera2.dylib")
    mv(source_file, target_file, force=true)

elseif Sys.iswindows()
    if isa(1, Int64)
        source_file = joinpath(sdk_dir, "ASI SDK", "lib", "x64", "ASICamera2.dll")
        target_file = joinpath(lib_target_dir, "ASICamera2.dll")
        mv(source_file, target_file, force=true)
    else
        source_file = joinpath(sdk_dir, "ASI SDK", "lib", "x86", "ASICamera2.dll")
        target_file = joinpath(lib_target_dir, "ASICamera2.dll")
        mv(source_file, target_file, force=true)
    end
end

# cleanup
rm(sdk_file)
rm(sdk_dir, recursive=true)

if Sys.isunix()
    rules_file = joinpath(par_dir, "asi.rules")
    println("\nPlease install the udev rules for the camera device, so that you can
    access it without root privileges. To install the rules, run
    'sudo install $rules_file /lib/udev/rules.d'\n
    or\n
    'sudo install $rules_file /etc/udev/rules.d'\n
    and reboot or relog and reconnect the camera. Then run\n
    'cat /sys/module/usbcore/parameters/usbfs_memory_mb'\n
    and make sure the result is 200.\n")
end


text:
using HDF5, JLD
using Distributions

#constructing a random bipartite graph (U = N, V= N) with degree varying from
#~0.15 N to ~0.65 N (effectively we disregard overlapping edges). Higher degree will
#result is a greater overlap and the submodular function will be further from being
#monotone.

N=100;
U_incid = Array{Array}(N);
for i=1:N
	U_incid[i] = []
end
for i = 1:N
		degree_of_i = 0.15*N + Int(round(rand()*N/2));
		for j = 1:degree_of_i
			random_neighbor = Int(round(rand()*(N-1)+1));
			push!(U_incid[i], random_neighbor);
		end
end

#assign random weights to vertices on both sides of the bipartition
m = 100*rand(N,1); #for vertices in U
weight = 100*rand(N, 1); #for vertices in V

### here is our very terrible submodular function
function jegelka_example(S)
	val = sum(m[S]);
	neighborhood = [];
	for i in S
			neighborhood = vcat(neighborhood, U_incid[i]) #compute the neighborhood of the selection of vertices from U
	end
	neigborhood = unique(neighborhood);
	val = val + 10*sqrt(sum(weight[neighborhood])); # one can vary the constant in front of the weight of the neighborhood as well.
	return val;
end

### now construct an optimization problem over the associated submodular polytope
a = rand(N)
### ??? I have no idea what the right syntax here is
g = DiffableFunction(x -> 1/2*sum((x-a).^2))
constraint(S) = AssocPoly(jegelka_example(S))
problem = CombiProblem(g, constraint)
### end ???

x = frank_wolfe_away(problem, verbose=true, maxiters=10)


text:
"""The http status module inspired from DRF status - https://github.com/encode/django-rest-framework/blob/master/rest_framework/status.py
"""
module HTTPStatus

include("codes.jl")

export is_information, is_success, is_redirect, is_client_error, is_server_error, HTTP_100_CONTINUE,
HTTP_101_SWITCHING_PROTOCOLS,
HTTP_200_OK,
HTTP_201_CREATED,
HTTP_202_ACCEPTED,
HTTP_203_NON_AUTHORITATIVE_INFORMATION,
HTTP_204_NO_CONTENT,
HTTP_205_RESET_CONTENT,
HTTP_206_PARTIAL_CONTENT,
HTTP_207_MULTI_STATUS,
HTTP_208_ALREADY_REPORTED,
HTTP_226_IM_USED,
HTTP_300_MULTIPLE_CHOICES,
HTTP_301_MOVED_PERMANENTLY,
HTTP_302_FOUND,
HTTP_303_SEE_OTHER,
HTTP_304_NOT_MODIFIED,
HTTP_305_USE_PROXY,
HTTP_306_RESERVED,
HTTP_307_TEMPORARY_REDIRECT,
HTTP_308_PERMANENT_REDIRECT,
HTTP_400_BAD_REQUEST,
HTTP_401_UNAUTHORIZED,
HTTP_402_PAYMENT_REQUIRED,
HTTP_403_FORBIDDEN,
HTTP_404_NOT_FOUND,
HTTP_405_METHOD_NOT_ALLOWED,
HTTP_406_NOT_ACCEPTABLE,
HTTP_407_PROXY_AUTHENTICATION_REQUIRED,
HTTP_408_REQUEST_TIMEOUT,
HTTP_409_CONFLICT,
HTTP_410_GONE,
HTTP_411_LENGTH_REQUIRED,
HTTP_412_PRECONDITION_FAILED,
HTTP_413_REQUEST_ENTITY_TOO_LARGE,
HTTP_414_REQUEST_URI_TOO_LONG,
HTTP_415_UNSUPPORTED_MEDIA_TYPE,
HTTP_416_REQUESTED_RANGE_NOT_SATISFIABLE,
HTTP_417_EXPECTATION_FAILED,
HTTP_418_IM_A_TEAPOT,
HTTP_422_UNPROCESSABLE_ENTITY,
HTTP_423_LOCKED,
HTTP_424_FAILED_DEPENDENCY,
HTTP_426_UPGRADE_REQUIRED,
HTTP_428_PRECONDITION_REQUIRED,
HTTP_429_TOO_MANY_REQUESTS,
HTTP_431_REQUEST_HEADER_FIELDS_TOO_LARGE,
HTTP_451_UNAVAILABLE_FOR_LEGAL_REASONS,
HTTP_500_INTERNAL_SERVER_ERROR,
HTTP_501_NOT_IMPLEMENTED,
HTTP_502_BAD_GATEWAY,
HTTP_503_SERVICE_UNAVAILABLE,
HTTP_504_GATEWAY_TIMEOUT,
HTTP_505_HTTP_VERSION_NOT_SUPPORTED,
HTTP_506_VARIANT_ALSO_NEGOTIATES,
HTTP_507_INSUFFICIENT_STORAGE,
HTTP_508_LOOP_DETECTED,
HTTP_509_BANDWIDTH_LIMIT_EXCEEDED,
HTTP_510_NOT_EXTENDED,
HTTP_511_NETWORK_AUTHENTICATION_REQUIRED

function is_information(code::Int)
    return 100 <= code <= 199
end

function is_success(code::Int)
    return 200 <= code <= 299
end

function is_redirect(code::Int)
    return 300 <= code <= 399
end

function is_client_error(code::Int)
    return 400 <= code <= 499
end

function is_server_error(code::Int)
    return 500 <= code <= 599
end

end # module

max_stars_repo_name: kracekumar/HTTPStatus

max_issues_repo_name: kracekumar/HTTPStatus

max_forks_repo_name: kracekumar/HTTPStatus


text:
using Test
using ITensors
using ITensors.ITensorNetworkMaps
using KrylovKit
using LinearAlgebra

include("utils.jl")

@testset "ITensorNetworkMaps.jl" begin
  N = 3 # Number of sites in the unit cell
  e⃗ = [n => mod1(n + 1, N) for n in 1:N]
  χ⃗ = Dict()
  χ⃗[1 => 2] = 3
  χ⃗[2 => 3] = 4
  χ⃗[3 => 1] = 5
  d = 2
  ψ = infmps(N; χ⃗=χ⃗, d=d)
  randn!.(ψ.data)
  T = transfer_matrix(ψ)
  v = randomITensor(input_inds(T))
  Tv_expected = ITensors.contract([v, reverse(ψ)..., prime(linkinds, dag(reverse(ψ)))...])
  T′v_expected = ITensors.contract([v, ψ..., prime(linkinds, dag(ψ))...])
  @test T(v) ≈ Tv_expected
  @test (2T)(v) ≈ 2Tv_expected
  @test (2T + I)(v) ≈ 2Tv_expected + v
  @test (2T + 3I)(v) ≈ 2Tv_expected + 3v
  @test (T + T)(v) ≈ 2Tv_expected
  @test T'(v) ≈ T′v_expected
  @test (T + T')(v) ≈ Tv_expected + T′v_expected
  @test (T + T)'(v) ≈ 2T′v_expected
  @test all([T T; T T]([v; v]) .≈ [2Tv_expected; 2Tv_expected])

  T⃗ = transfer_matrices(ψ)

  @test (T⃗[1] * T⃗[2] * T⃗[3])(v) ≈ Tv_expected
  @test (T⃗[1] * T⃗[2] * T⃗[3] + 3I)(v) ≈ Tv_expected + 3v
  @test (T⃗[1] * T⃗[2] * T⃗[3])'(v) ≈ T′v_expected
  @test (I * T⃗[1] * T⃗[2] * T⃗[3])(v) ≈ Tv_expected
  @test (2T⃗[1] * T⃗[2] * T⃗[3])(v) ≈ 2Tv_expected
  @test (T⃗[1] * T⃗[2] * T⃗[3] + 2I)(v) ≈ Tv_expected + 2v
  @test ((T⃗[1] * T⃗[2] + T⃗[1] * T⃗[2]) * T⃗[3])(v) ≈ 2Tv_expected
  @test ((T⃗[1] * T⃗[2] + T⃗[1] * T⃗[2]) * T⃗[3])'(v) ≈ 2T′v_expected
  @test ((T⃗[1] * T⃗[2] + T⃗[1] * T⃗[2]) * T⃗[3] + 3I)(v) ≈ 2Tv_expected + 3v

  dk, vk = eigsolve(T, v)
  for n in 1:length(dk)
    @test norm((T - dk[n]I)(vk[n])) ≈ 0 atol = 1e-12
  end
end


text:
_bdd_nextid = 1
_bdd_mutex = Threads.SpinLock()
_bdd_int_typemax = typemax(Int)

"""A Binary Decision Diagram.

- `index`: the vertex variable (-1 if terminal vertex
- `low`: low child vertex of BDD (undef if terminal vertex)
- `high`: high child vertex of BDD (undef if terminal vertex)
- `value`: terminal boolean value
- `id`: unique identifier
"""
mutable struct Bdd
  "Root vertex variable index (-1 if terminal vertex)."
  index::Int
  "Low child vertex of BDD (undef if terminal vertex)."
  low::Bdd
  "High child vertex of BDD (undef if terminal vertex)."
  high::Bdd
  "Terminal boolean value."
  value::Bool
  "Unique identifier."
  id::Int
  "Constructs a terminal."
  function Bdd(v::Bool)::Bdd
    α = new()
    α.index, α.value = -1, v
    lock(_bdd_mutex)
    global _bdd_nextid
    α.id = _bdd_nextid
    _bdd_nextid = _bdd_nextid + 1 % _bdd_int_typemax
    unlock(_bdd_mutex)
    return α
  end
  "Constructs a variable."
  function Bdd(i::Integer, low::Bdd, high::Bdd)::Bdd
    α = new()
    α.index, α.low, α.high = i, low, high
    lock(_bdd_mutex)
    global _bdd_nextid
    α.id = _bdd_nextid
    _bdd_nextid = _bdd_nextid + 1 % _bdd_int_typemax
    unlock(_bdd_mutex)
    return α
  end
  "Forcefuly construct a node with given parameters. Should not be used."
  function Bdd(i::Integer, id::Integer, low::Bdd, high::Bdd)::Bdd
    α = new()
    α.index, α.low, α.high, α.id = i, low, high, id
    return α
  end
end
export Bdd
const BinaryDecisionDiagram = Bdd
export BinaryDecisionDiagram

const ⊤ = Bdd(true)
const ⊥ = Bdd(false)
export ⊤, ⊥

"""
    shallowhash(α::Bdd, h::UInt = UInt64(0))

Returns a shallow hash for the given node (not BDD as a whole).
"""
@inline shallowhash(α::Bdd, h::UInt = UInt64(0)) = hash((α.id, α.value, α.index), h)
export shallowhash

"""
    hash(α::Bdd, h::UInt)::UInt

Returns a unique hash for the whole BDD.
"""
function Base.hash(α::Bdd, h::UInt)::UInt
  H = Tuple{Bool, Int}[]
  foreach(x -> push!(H, (is_term(x) ? x.value : false, x.index)), α)
  return hash(H, h)
end

"""
    is_term(α::Bdd)::Bool

Returns whether this Bdd node is terminal.
"""
@inline is_term(α::Bdd)::Bool = (α.index < 0) || (!isdefined(α, :low) && !isdefined(α, :high))
export is_term

"""
    is_⊤(α::Bdd)::Bool

Returns whether the given Bdd node represents a ⊤.
"""
@inline is_⊤(α::Bdd)::Bool = is_term(α) && α.value
export is_⊤

"""
    is_⊥(α::Bdd)::Bool

Returns whether the given Bdd node represents a ⊥.
"""
@inline is_⊥(α::Bdd)::Bool = is_term(α) && !α.value
export is_⊥

"""
    is_var(α::Bdd)::Bool

Returns whether the given Bdd node represents a variable.
"""
@inline is_var(α::Bdd)::Bool = (isdefined(α, :low) && is_⊥(α.low)) && (isdefined(α, :high) && is_⊤(α.high))

"""
    is_lit(α::Bdd)::Bool

Returns whether the given Bdd node represents a literal.
"""
@inline is_lit(α::Bdd)::Bool = isdefined(α, :low) && isdefined(α, :high) && is_term(α.low) && is_term(α.high)

"""
    is_atom(α::Bdd)::Bool

Returns whether the given Bdd node is an atomic formula (i.e. a variable, ⊥, ⊤, or literal).
"""
@inline is_atom(α::Bdd)::Bool = is_term(α) || is_lit(α)
export is_var, is_lit, is_atom

"""
    (¬)(α::Bdd)::Bdd
    (¬)(x::Integer)::Bdd
    (¬)(x::Bool)::Bool

Negates this boolean function.
"""
@inline (¬)(α::Bdd)::Bdd = is_⊤(α) ? ⊥ : is_⊥(α) ? ⊤ : Bdd(α.index, ¬α.low, ¬α.high)
@inline (¬)(x::Integer)::Bdd = x > 0 ? Bdd(x, ⊤, ⊥) : x < 0 ? Bdd(-x, ⊥, ⊤) : ⊥
@inline (¬)(x::Bool)::Bool = !x
export ¬

"""
    (∧)(α::Bdd, β::Bdd)::Bdd
    (∧)(x::Integer, β::Bdd)::Bdd
    (∧)(α::Bdd, x::Integer)::Bdd
    (∧)(x::Integer, y::Integer)::Bdd
    (∧)(x::Bool, y::Bool)::Bool

Returns a conjunction over the given boolean functions.
"""
@inline (∧)(α::Bdd, β::Bdd)::Bdd = apply(α, β, &)
@inline (∧)(x::Integer, β::Bdd)::Bdd = apply(bdd_var(x), β, &)
@inline (∧)(α::Bdd, x::Integer)::Bdd = apply(α, bdd_var(x), &)
@inline (∧)(x::Integer, y::Integer)::Bdd = apply(bdd_var(x), bdd_var(y), &)
@inline (∧)(x::Bool, y::Bool)::Bool = x & y
export ∧

"""
    and(α::Bdd, β::Bdd)::Bdd
    and(x::Integer, β::Bdd)::Bdd
    and(α::Bdd, x::Integer)::Bdd
    and(x::Integer, y::Integer)::Bdd
    and(Φ::Vararg{Bdd})::Bdd
    and(Φ::AbstractVector{Bdd})::Bdd
    and(Φ::Vararg{T})::Bdd where T <: Integer
    and(Φ::AbstractVector{T})::Bdd where T <: Integer
    and(Φ::Vararg{Union{T, Bdd}})::Bdd where T <: Integer

Returns a conjunction over the given boolean functions.
"""
@inline and(α::Bdd, β::Bdd)::Bdd = α ∧ β
@inline and(x::Integer, β::Bdd)::Bdd = x ∧ β
@inline and(α::Bdd, x::Integer)::Bdd = α ∧ x
@inline and(x::Integer, y::Integer)::Bdd = x ∧ y
function and(Φ::Vararg{Bdd})::Bdd
  α = first(Φ)
  for i ∈ 2:length(Φ) α = α ∧ Φ[i] end
  return α
end
@inline and(Φ::AbstractVector{Bdd})::Bdd = and(Φ...)
@inline function and(Φ::Vararg{T})::Bdd where T <: Integer and(bdd_var.(Φ)...) end
@inline function and(Φ::AbstractVector{T})::Bdd where T <: Integer and(bdd_var.(Φ)...) end
function and(Φ::Vararg{Union{T, Bdd}})::Bdd where T <: Integer
  f = first(Φ)
  α = f isa Integer ? bdd_var(f) : f
  for i ∈ 2:length(Φ) α = α ∧ Φ[i] end
  return α
end
export and

"""
    (∨)(α::Bdd, β::Bdd)::Bdd
    (∨)(x::Integer, β::Bdd)::Bdd
    (∨)(α::Bdd, x::Integer)::Bdd
    (∨)(x::Integer, y::Integer)::Bdd
    (∨)(x::Bool, y::Bool)::Bool

Returns a disjunction over the given boolean functions.
"""
@inline (∨)(α::Bdd, β::Bdd)::Bdd = apply(α, β, |)
@inline (∨)(x::Integer, β::Bdd)::Bdd = apply(bdd_var(x), β, |)
@inline (∨)(α::Bdd, x::Integer)::Bdd = apply(α, bdd_var(x), |)
@inline (∨)(x::Integer, y::Integer)::Bdd = apply(bdd_var(x), bdd_var(y), |)
@inline (∨)(x::Bool, y::Bool)::Bool = x | y
export ∨

"""
    or(α::Bdd, β::Bdd)::Bdd
    or(x::Integer, β::Bdd)::Bdd
    or(α::Bdd, x::Integer)::Bdd
    or(x::Integer, y::Integer)::Bdd
    or(Φ::Vararg{Bdd})::Bdd
    or(Φ::AbstractVector{Bdd})::Bdd
    or(Φ::Vararg{T})::Bdd where T <: Integer
    or(Φ::AbstractVector{T})::Bdd where T <: Integer
    or(Φ::Vararg{Union{T, Bdd}})::Bdd where T <: Integer

Returns a disjunction over the given boolean functions.
"""
@inline or(α::Bdd, β::Bdd)::Bdd = α ∨ β
@inline or(x::Integer, β::Bdd)::Bdd = x ∨ β
@inline or(α::Bdd, x::Integer)::Bdd = α ∨ x
@inline or(x::Integer, y::Integer)::Bdd = x ∨ y
@inline function or(Φ::Vararg{Bdd})::Bdd
  α = first(Φ)
  for i ∈ 2:length(Φ) α = α ∨ Φ[i] end
  return α
end
@inline or(Φ::AbstractVector{Bdd})::Bdd = or(Φ...)
@inline function or(Φ::Vararg{T})::Bdd where T <: Integer or(bdd_var.(Φ)...) end
@inline function or(Φ::AbstractVector{T})::Bdd where T <: Integer or(bdd_var.(Φ)...) end
function or(Φ::Vararg{Union{T, Bdd}})::Bdd where T <: Integer
  f = first(Φ)
  α = f isa Integer ? bdd_var(f) : f
  for i ∈ 2:length(Φ) α = α ∨ Φ[i] end
  return α
end
export or

"""
Returns a xor of the given boolean functions.
"""
@inline Base.:⊻(α::Bdd, β::Bdd)::Bdd = apply(α, β, ⊻)
@inline Base.:⊻(x::Integer, β::Bdd)::Bdd = apply(bdd_var(x), β, ⊻)
@inline Base.:⊻(α::Bdd, x::Integer)::Bdd = apply(α, bdd_var(x), ⊻)

@inline (→)(α::Bdd, β::Bdd)::Bdd = (¬α) ∨ β
@inline (→)(x::Integer, β::Bdd)::Bdd = (¬x) ∨ β
@inline (→)(α::Bdd, x::Integer)::Bdd = (¬α) ∨ x
@inline (→)(x::Integer, y::Integer)::Bdd = (¬x) ∨ y
@inline (→)(x::Bool, y::Bool)::Bool = (!x) | y
export →

"Returns whether the two given boolean functions are equivalent."
@inline Base.:(==)(α::Bdd, β::Bdd)::Bool = is_⊤(apply(α, β, ==))
@inline Base.:(==)(x::Integer, β::Bdd)::Bool = is_var(β) && β.index == x && ((x > 0 && β.low == ⊥ && β.high == ⊤) || (x < 0 && β.low == ⊤ && β.high == ⊥))
@inline Base.:(==)(α::Bdd, y::Integer)::Bool = y == α
@inline Base.isequal(α::Bdd, β::Bdd)::Bool = α == β
@inline Base.isequal(x::Integer, β::Bdd)::Bool = x == β
@inline Base.isequal(α::Bdd, y::Integer)::Bool = y == α

"Returns whether the two given boolean functions are not equivalent."
@inline Base.:(!=)(α::Bdd, β::Bdd)::Bool = !(α == β)
@inline Base.:(!=)(x::Integer, β::Bdd)::Bool = !(x == β)
@inline Base.:(!=)(α::Bdd, y::Integer)::Bool = !(y == α)

"Returns a new terminal node of given boolean value."
@inline terminal(v::Bool)::Bdd = Bdd(v)
export terminal

"Returns a Bdd representing a single variable. If negative, negate variable."
@inline bdd_var(i::Integer)::Bdd = i > 0 ? Bdd(i, ⊥, ⊤) : Bdd(-i, ⊤, ⊥)
export bdd_var

"Returns 0 if x is not a literal; else returns the literal's sign."
@inline Base.sign(x::Bdd)::Int = !is_lit(x) ? 0 : x.low == ⊥ ? 1 : -1
@inline Base.signbit(x::Bdd)::Bool = sign(x) == -1
"Returns 0 if x is not a literal; else returns an integer representation of x."
@inline to_int(x::Bdd)::Int = !is_lit(x) ? 0 : x.low == ⊥ ? x.index : -x.index
export to_int

"Return string representation of Bdd α."
function Base.string(α::Bdd; max_depth::Integer = 20)::String
  s = ""
  S = Tuple{Bdd, Int, Char}[(α, 0, '\0')]
  d = 0
  while !isempty(S)
    v, indent, c = pop!(S)
    for i ∈ 1:indent s *= "|  " end
    s *= c == '\0' ? '@' : c
    if is_term(v)
      s *= " (value=$(v.value), id=$(v.id))\n"
    else
      s *= " (index=$(v.index), id=$(v.id))\n"
      push!(S, (v.high, indent + 1, '+'))
      push!(S, (v.low, indent + 1, '-'))
    end
    d += 1
    if d > max_depth
      s *= "...\nOmitting some nodes. BDD is too large."
      break
    end
  end
  return s
end
Base.show(io::Core.IO, α::Bdd; kwargs...) = print(io, string(α; kwargs...))
Base.print(io::Core.IO, α::Bdd; kwargs...) = print(io, string(α; kwargs...))

let V::Set{UInt64}, Q::Vector{Bdd}
  function Base.iterate(α::Bdd, state=1)::Union{Nothing, Tuple{Bdd, Integer}}
    if state == 1
      Q = Bdd[α]
      V = Set{UInt64}(objectid(α))
    end
    if isempty(Q) return nothing end
    v = pop!(Q)
    if !is_term(v)
      l, h = v.low, v.high
      p, q = objectid(l), objectid(h)
      if q ∉ V push!(Q, h); push!(V, q) end
      if p ∉ V push!(Q, l); push!(V, p) end
    end
    return v, state+1
  end
end

function Base.foreach(f::Function, α::Bdd)
  V = Set{UInt64}(shallowhash(α))
  Q = Bdd[α]
  while !isempty(Q)
    v = pop!(Q)
    if !is_term(v)
      l, h = v.low, v.high
      p, q = objectid(l), objectid(h)
      if q ∉ V push!(Q, h); push!(V, q) end
      if p ∉ V push!(Q, l); push!(V, p) end
    end
    f(v)
  end
  nothing
end

function Base.collect(α::Bdd)::Vector{Bdd}
  V = Set{UInt64}(objectid(α))
  Q = Bdd[α]
  C = Vector{Bdd}()
  while !isempty(Q)
    v = pop!(Q)
    if !is_term(v)
      l, h = v.low, v.high
      p, q = objectid(l), objectid(h)
      if q ∉ V push!(Q, h); push!(V, q) end
      if p ∉ V push!(Q, l); push!(V, p) end
    end
    push!(C, v)
  end
  return C
end

"Returns a `Vector{Bdd}` containing all nodes in `α` in post-order."
@inline function postorder(α::Bdd)::Vector{Bdd}
  O = Vector{Bdd}()
  postorder_step(α, Set{UInt64}(), O)
  return O
end
export postorder

function postorder_step(α::Bdd, V::Set{UInt64}, O::Vector{Bdd})
  push!(V, objectid(α))
  if !is_term(α)
    l, h = α.low, α.high
    p, q = objectid(l), objectid(h)
    if p ∉ V
      push!(V, p)
      postorder_step(l, V, O)
    end
    if q ∉ V
      push!(V, q)
      postorder_step(h, V, O)
    end
  end
  push!(O, α)
  nothing
end

"""Reduce a Bdd rooted at α inplace, removing duplicate nodes and redundant sub-trees. Returns
canonical representation of α."""
function reduce!(α::Bdd)::Bdd
  if is_term(α) return α end

  V = Dict{Int, Vector{Bdd}}()
  foreach(function(v::Bdd)
            i = v.index
            haskey(V, i) ? push!(V[i], v) : V[i] = Bdd[v]
          end, α)
  nid = 0
  G = Dict{Int, Bdd}()
  I = sort!(collect(keys(V)), rev=true); pop!(I); pushfirst!(I, -1)
  for i ∈ I
    Q = Vector{Tuple{Tuple{Int, Int}, Bdd}}()
    for v ∈ V[i]
      if is_term(v) push!(Q, ((Int(v.value), -1), v))
      elseif v.low.id == v.high.id v.id = v.low.id
      else push!(Q, ((v.low.id, v.high.id), v)) end
    end
    sort!(Q, by=first)
    local oldk::Tuple{Int, Int} = (-2, -2)
    for (k, v) ∈ Q
      if k == oldk v.id = nid
      else
        nid += 1
        v.id = nid
        G[nid] = v
        if !is_term(v)
          v.low = G[v.low.id]
          v.high = G[v.high.id]
        end
        oldk = k
      end
    end
  end
  return G[α.id]
end
export reduce!

function Base.copy(α::Bdd)::Bdd
  if is_term(α) return is_⊤(α) ? Bdd(true) : Bdd(false) end
  return Bdd(α.index, α.low, α.high)
end

function Base.deepcopy(α::Bdd)::Bdd
  if is_term(α) return is_⊤(α) ? Bdd(true) : Bdd(false) end
  return Bdd(α.index, copy(α.low), copy(α.high))
end

"Returns a Bdd canonical representation of α ⊕ β, where ⊕ is some binary operator."
@inline apply(α::Bdd, β::Bdd, ⊕)::Bdd = reduce!(apply_step(α, β, ⊕, Dict{Tuple{Int, Int}, Bdd}()))
export apply

"""Recursively computes α ⊕ β. If the result was already computed as an intermediate result, return
the cached result in T."""
function apply_step(α::Bdd, β::Bdd, ⊕, T::Dict{Tuple{Int, Int}, Bdd})::Bdd
  local k::Tuple{Int, Int} = (α.id, β.id)
  if haskey(T, k) return T[k] end

  local r::Bdd
  if is_term(α) && is_term(β) r = Bdd(α.value ⊕ β.value)
  else
    local i::Int = typemax(Int)
    local j::Int = i

    if !is_term(α) i = α.index end
    if !is_term(β) j = β.index end
    m = min(i, j)

    local l1::Bdd, h1::Bdd
    if α.index == m l1, h1 = α.low, α.high
    else l1 = h1 = α end

    local l2::Bdd, h2::Bdd
    if β.index == m l2, h2 = β.low, β.high
    else l2 = h2 = β end

    l = apply_step(l1, l2, ⊕, T)
    h = apply_step(h1, h2, ⊕, T)
    r = Bdd(m, l, h)
  end

  T[k] = r
  return r
end

@inline function vec_to_dict(X::AbstractArray{T})::Dict{Integer, Bool} where T <: Integer
  D = Dict{Integer, Bool}()
  for x ∈ X if x < 0 D[-x] = false else D[x] = true end end
  return D
end

"Returns a new reduced Bdd restricted to instantiation X."
@inline function restrict(α::Bdd, X::Dict{T, Bool})::Bdd where T <: Integer reduce!(restrict_step(α, X, Dict{Int, Bdd}())) end
@inline function restrict(α::Bdd, X::AbstractArray{T})::Bdd where T <: Integer reduce!(restrict_step(α, vec_to_dict(X), Dict{Int, Bdd}())) end
@inline restrict(α::Bdd, X::Union{AbstractArray{Bool}, BitVector})::Bdd = reduce!(restrict_step(α, X, Dict{Int, Bdd}()))
@inline restrict(α::Bdd, x::Integer)::Bdd = reduce!(restrict_step(α, abs(x), x > 0, Dict{Int, Bdd}()))
export restrict

"Returns a new reduced Bdd restricted to instantiation X."
@inline function Base.:|(α::Bdd, X::Dict{T, Bool})::Bdd where T <: Integer restrict(α, X) end
@inline function Base.:|(α::Bdd, X::AbstractArray{T})::Bdd where T <: Integer restrict(α, X) end
@inline Base.:|(α::Bdd, x::Integer)::Bdd = restrict(α, x)
@inline Base.:|(α::Bdd, X::Union{AbstractArray{Bool}, BitVector})::Bdd = restrict(α, X)

@inline function (α::Bdd)(X::Dict{T, Bool})::Bool where T <: Integer is_⊤(restrict(α, X)) end
@inline function (α::Bdd)(X::AbstractArray{T})::Bool where T <: Integer is_⊤(α|X) end
@inline (α::Bdd)(X::Union{AbstractArray{Bool}, BitVector})::Bool = is_⊤(α|X)
@inline (α::Bdd)(x::Integer)::Bool = is_⊤(α|x)

"Returns a new Bdd restricted to instantiation X."
function restrict_step(α::Bdd, X::Dict{T, Bool}, V::Dict{Int, Bdd})::Bdd where T <: Integer
  if haskey(V, α.id) return V[α.id]
  elseif is_term(α) return copy(α) end
  x = α.index
  if !haskey(X, x)
    l = restrict_step(α.low, X, V)
    h = restrict_step(α.high, X, V)
    β = Bdd(x, l, h)
    V[α.id] = β
    return β
  end
  if X[x] return restrict_step(α.high, X, V) end
  return restrict_step(α.low, X, V)
end

"Returns a new Bdd restricted to instantiation X."
function restrict_step(α::Bdd, X::Union{BitVector, AbstractArray{Bool}}, V::Dict{Int, Bdd})::Bdd
  if haskey(V, α.id) return V[α.id]
  elseif is_term(α) return copy(α) end
  x = α.index
  if x > length(X)
    l = restrict_step(α.low, X, V)
    h = restrict_step(α.high, X, V)
    β = Bdd(x, l, h)
    V[α.id] = β
    return β
  end
  if X[x] return restrict_step(α.high, X, V) end
  return restrict_step(α.low, X, V)
end

function restrict_step(α::Bdd, v::Integer, y::Bool, V::Dict{Int, Bdd})::Bdd
  if haskey(V, α.id) return V[α.id]
  elseif is_term(α) return copy(α) end
  x = α.index
  if x != v
    l = restrict_step(α.low, v, y, V)
    h = restrict_step(α.high, v, y, V)
    β = Bdd(x, l, h)
    V[α.id] = β
    return β
  end
  if y return restrict_step(α.high, v, y, V) end
  return restrict_step(α.low, v, y, V)
end

struct Permutations
  V::Vector{Int}
  m::Int
end
@inline Base.length(P::Permutations)::Int = P.m

"Compute all possible valuations of scope V."
@inline valuations(V::Union{Set{T}, AbstractVector{T}, UnitRange{T}}) where T <: Integer = Permutations(collect(V), 2^length(V))
export valuations
function Base.iterate(P::Permutations, state=0)::Union{Nothing, Tuple{Dict{Int, Bool}, Int}}
  s = state + 1
  if state == 0 return Dict{Int, Bool}(broadcast(x -> abs(x) => false, P.V)), s end
  if state >= P.m return nothing end
  return Dict{Int, Bool}((i -> P.V[i] => (state >> (i-1)) & 1 == 1).(1:length(P.V))), s
end

"Computes all possible valuations of scope V and returns as a BitMatrix. Up to 64 variables."
function all_valuations(V::Union{Set{T}, AbstractVector{T}, UnitRange{T}})::BitMatrix where T <: Integer
  m = length(V)
  n = 2^m
  M = BitMatrix(undef, (n, m))
  for i ∈ 1:n
    p = i-1
    for j ∈ 1:m
      M[i, j] = p & 1
      p >>= 1
    end
  end
  return M
end
export all_valuations

struct ConjoinedPermutations
  V::Vector{Int}
  m::Int
end
@inline Base.length(P::ConjoinedPermutations)::Int = P.m

"Computes all possible valuations of scope V as conjunctions."
@inline conjunctions(V::Union{Set{T}, AbstractVector{T}, UnitRange{T}}) where T <: Integer = ConjoinedPermutations(sort!(collect(V)), 2^length(V))
export conjunctions
function Base.iterate(P::ConjoinedPermutations, state=0)::Union{Nothing, Tuple{Bdd, Int}}
  s = state + 1
  if state >= P.m return nothing end
  local V::Vector{Int}
  if state == 0 V = broadcast(-, P.V)
  else V = (i -> (state >> (i-1)) & 1 == 1 ? P.V[i] : -P.V[i]).(1:length(P.V)) end
  local α::Bdd
  f = true
  for v ∈ Iterators.reverse(V)
    if f
      α = bdd_var(v)
      f = false
    else
      α = v > 0 ? Bdd(v, ⊥, α) : Bdd(-v, α, ⊥)
    end
  end
  return α, s
end

struct ConvalPermutations
  V::Vector{Int}
  m::Int
end
@inline Base.length(P::ConvalPermutations)::Int = P.m

"Computes all possible valuations of scope V as both conjunctions and instantiation values."
@inline convals(V::Union{Set{T}, AbstractVector{T}, UnitRange{T}}) where T <: Integer = ConvalPermutations(sort!(collect(V)), 2^length(V))
export convals
function Base.iterate(P::ConvalPermutations, state=0)::Union{Nothing, Tuple{Tuple{Bdd, Dict{Int, Bool}}, Int}}
  s = state + 1
  if state >= P.m return nothing end
  local V::Vector{Int}
  if state == 0 V = broadcast(-, P.V)
  else V = (i -> (state >> (i-1)) & 1 == 1 ? P.V[i] : -P.V[i]).(1:length(P.V)) end
  local α::Bdd
  f = true
  for v ∈ Iterators.reverse(V)
    if f
      α = bdd_var(v)
      f = false
    else
      α = v > 0 ? Bdd(v, ⊥, α) : Bdd(-v, α, ⊥)
    end
  end
  return (α, Dict{Int, Bool}((i -> i > 0 ? i => true : -i => false).(V))), s
end

"Performs Shannon's Decomposition on the Bdd α, given a variable to isolate."
function shannon(α::Bdd, v::Integer)::Tuple{Bdd, Bdd, Bdd, Bdd}
  return (bdd_var(v), α|Dict{Integer, Bool}(v=>true), bdd_var(-v), α|Dict{Integer, Bool}(v=>false))
end

"Performs Shannon's Decomposition on the Bdd α, given a set of variables to isolate."
function shannon(α::Bdd, V::Union{Set{T}, AbstractVector{T}})::Vector{Tuple{Bdd, Bdd}} where T <: Integer
  Δ = Vector{Tuple{Bdd, Bdd}}()
  for (β, X) ∈ convals(V) push!(Δ, (β, α|X)) end
  return Δ
end

"""Performs Shannon's Decomposition on the Bdd α, given a set of variables to isolate. Any
decomposition that results in a ⊥ is discarded."""
function shannon!(α::Bdd, V::Union{Set{T}, AbstractVector{T}})::Vector{Tuple{Bdd, Bdd}} where T <: Integer
  Δ = Vector{Tuple{Bdd, Bdd}}()
  for (β, X) ∈ convals(V)
    ϕ = α|X
    if is_⊥(ϕ) continue end
    push!(Δ, (β, ϕ))
  end
  return Δ
end

export shannon, shannon!

"Eliminate a variable through disjunction. Equivalent to the expression (ϕ|x ∨ ϕ|¬x)."
@inline eliminate(α::Bdd, v::Integer)::Bdd = reduce!(eliminate_step(α, v))
function eliminate_step(α::Bdd, v::Integer)::Bdd
  if is_term(α) return copy(α) end
  if α.index == v return α.low ∨ α.high end
  # If idempotent (which is the case), then recursion suffices.
  l = eliminate_step(α.low, v)
  h = eliminate_step(α.high, v)
  return Bdd(α.index, l, h)
end
@inline eliminate(α::Bdd, V::Union{Set, BitSet, AbstractVector{T}}) where T <: Integer = isempty(V) ? α : reduce!(eliminate_step(α, V))
function eliminate_step(α::Bdd, V::Union{Set, BitSet, AbstractVector{T}})::Bdd where T <: Integer
  if is_term(α) return copy(α) end
  if α.index ∈ V return eliminate_step(α.low ∨ α.high, V) end
  l = eliminate_step(α.low, V)
  h = eliminate_step(α.high, V)
  return Bdd(α.index, l, h)
end
export eliminate

"""
    forget(α::Bdd, x::Integer)::Bdd = eliminate(α, x)
    forget(α::Bdd, X::Union{Set, BitSet, AbstractVector{T}}) where T <: Integer

Returns the resulting BDD after applying the `forget` operation. Equivalent to \$\\phi|_x \\vee \\phi|_{\\neg x}\$.
"""
@inline forget(α::Bdd, x::Integer)::Bdd = eliminate(α, x)
@inline forget(α::Bdd, X::Union{Set, BitSet, AbstractVector{T}}) where T <: Integer = eliminate(α, X)
export forget

"Marginalize a variable through some binary operation."
@inline marginalize(α::Bdd, v::Integer, ⊕)::Bdd = is_term(α) ? Bdd(α.value ⊕ α.value) : reduce!(marginalize_step(α, v, ⊕))
export marginalize
function marginalize_step(α::Bdd, v::Integer, ⊕)::Bdd
  if α.index == v return apply(α.low, α.high, ⊕) end
  if is_term(α) return Bdd(α.value ⊕ α.value) end
  l, h = α.low, α.high
  return Bdd(α.index, marginalize_step(l, v, ⊕), marginalize_step(h, v, ⊕))
end

"Returns all variables in this formula as a Vector{Int}."
@inline scope(α::Bdd)::Vector{Int} = collect(scopeset(α))
export scope

"Returns all variables in this formula as a Set{Int}."
function scopeset(α::Bdd)::Set{Int}
  if is_term(α) return Set{Int}() end
  if is_var(α) return Set{Int}(α.index) end
  Sc = Set{Int}()
  foreach(function(v::Bdd)
            if !is_term(v) push!(Sc, v.index) end
          end, α)
  return Sc
end
export scopeset

"Returns whether the formula (i.e. BDD) mentions a variable."
function mentions(α::Bdd, x::Integer)::Bool
  V = Set{UInt64}(objectid(α))
  Q = Bdd[α]
  while !isempty(Q)
    v = pop!(Q)
    if !is_term(v)
      if v.index == x return true end
      l, h = v.low, v.high
      p, q = objectid(l), objectid(h)
      if q ∉ V push!(Q, h); push!(V, q) end
      if p ∉ V push!(Q, l); push!(V, p) end
    end
  end
  return false
end
function mentions(α::Bdd, X::AbstractVector{T})::Bool where T <: Integer
  M = Set{Int}(X)
  V = Set{UInt64}(objectid(α))
  Q = Bdd[α]
  while !isempty(Q)
    v = pop!(Q)
    if !is_term(v)
      if v.index ∈ M return true end
      l, h = v.low, v.high
      p, q = objectid(l), objectid(h)
      if q ∉ V push!(Q, h); push!(V, q) end
      if p ∉ V push!(Q, l); push!(V, p) end
    end
  end
  return false
end
export mentions
@inline Base.:∈(x::Integer, α::Bdd)::Bool = mentions(α, x)
@inline Base.:∉(x::Integer, α::Bdd)::Bool = !mentions(α, x)
@inline Base.:∈(X::AbstractVector{T}, α::Bdd) where T <: Integer = mentions(α, X)
@inline Base.:∉(X::AbstractVector{T}, α::Bdd) where T <: Integer = !mentions(α, X)

"""Returns an approximation (does not account for some repeated nodes) of how many times each
variable is mentioned in α."""
function culledfreqs(α::Bdd)::Dict{Int, Int}
  F = Dict{Int, Int}()
  !is_term(α) && (F[α.index] = 1)
  V = Set{UInt64}(objectid(α))
  Q = Bdd[α]
  while !isempty(Q)
    v = pop!(Q)
    if !is_term(v)
      l, h = v.low, v.high
      p, q = objectid(l), objectid(h)
      !is_term(l) && (haskey(F, l.index) ? F[l.index] += 1 : F[l.index] = 1)
      !is_term(h) && (haskey(F, h.index) ? F[h.index] += 1 : F[h.index] = 1)
      if q ∉ V push!(Q, h); push!(V, q) end
      if p ∉ V push!(Q, l); push!(V, p) end
    end
  end
  return F
end
export culledfreqs

"Assumes ϕ is a full conjunction of literals. Returns ϕ as a zero-one vector and its scope."
function lit_vec(α::Bdd)::Tuple{BitVector, Vector{Int}}
  X = BitVector()
  S = Vector{Int}()
  Q = Bdd[α]
  # BDDs guarantee variables are top-down increasing.
  while !isempty(Q)
    v = pop!(Q)
    if !is_term(v)
      push!(S, v.index)
      push!(X, is_⊥(v.low) ? true : false)
      if !is_term(v.low) push!(Q, v.low) end
      if !is_term(v.high) push!(Q, v.high) end
    end
  end
  return X, S
end
export lit_vec

"Returns α as an integer literal. Assumes α is a leaf node."
@inline to_lit(α::Bdd, ::Type{T} = Int32) where T <: Integer = convert(T, is_⊥(α.low) ? α.index : -α.index)
export to_lit

"""Translates a cardinality constraint in normal pseudo-boolean constraint form into a BDD.

Since cardinality constraints correspond to having coefficients set to one, we ignore the C's.

Argument L corresponds to the vector of literals to be chosen from; b is how many literals in L are
selected.

See Eén and Sörensson 2006."""
@inline function from_npbc(L::Vector{T}, b::Integer)::Bdd where T <: Integer
  return reduce!(from_npbc_step(L, b, length(L), 0, Dict{Tuple{Integer, Integer}, Bdd}(),
                                Bdd(false), Bdd(true)))
end
function from_npbc_step(L::Vector{T}, b::Integer, n::Integer, s::Integer, M::Dict{Tuple{Integer, Integer}, Bdd},
                        reuse_⊥::Bdd, reuse_⊤::Bdd)::Bdd where T <: Integer
  if s >= b return reuse_⊤
  elseif s + n < b return reuse_⊥ end
  k = (n, s)
  if !haskey(M, k)
    v = L[n]
    positive = v > 0
    if positive
      h_s, l_s = s+1, s
    else
      h_s, l_s = s, s+1
      v = -v
    end
    m = n-1
    h = from_npbc_step(L, b, m, h_s, M, reuse_⊥, reuse_⊤)
    l = from_npbc_step(L, b, m, l_s, M, reuse_⊥, reuse_⊤)
    M[k] = Bdd(v, l, h)
  end
  return M[k]
end

"Constructs a BDD mapping to true if at least n literals in L are in the input; otherwise false."
@inline atleast!(n::Integer, L::Vector{T}) where T <: Integer = from_npbc(sort!(L, by=x->abs(x), rev=true), n)
"Constructs a BDD mapping to true if at least n literals in L are in the input; otherwise false."
@inline atleast(n::Integer, L::Vector{T}) where T <: Integer = atleast!(n, copy(L))
"Constructs a BDD mapping to true if at most n literals in L are in the input; otherwise false."
@inline atmost!(n::Integer, L::Vector{T}) where T <: Integer = (m = length(L); L .= -L; from_npbc(sort!(L, by=x->abs(x), rev=true), m-n))
"Constructs a BDD mapping to true if at most n literals in L are in the input; otherwise false."
@inline atmost(n::Integer, L::Vector{T}) where T <: Integer = atmost!(n, copy(L))
"Constructs a BDD mapping to true if exactly n literals in L are in the input; otherwise false."
@inline exactly!(n::Integer, L::Vector{T}) where T <: Integer = (α = from_npbc(sort!(L, by=x->abs(x), rev=true), n); β = from_npbc(L .= -L, length(L)-n); α ∧ β)
"Constructs a BDD mapping to true if exactly n literals in L are in the input; otherwise false."
@inline exactly(n::Integer, L::Vector{T}) where T <: Integer = exactly!(n, copy(L))
export atleast, atmost, exactly, atleast!, atmost!, exactly!

"Returns the number of nodes in the BDD graph."
@inline Base.size(α::Bdd)::Int = (n = 0; foreach(x -> n += 1, α); n)

"Returns whether a variable x appears as a positive literal in α, given that α is a conjunction of literals."
function lit_val(α::Bdd, x::Integer)::Bool
  Q = Bdd[α]
  u = abs(x)
  # BDDs guarantee variables are top-down increasing.
  while !isempty(Q)
    v = pop!(Q)
    if !is_term(v)
      if u == v.index return (x > 0) == is_⊥(v.low) end
      if !is_term(v.low) push!(Q, v.low) end
      if !is_term(v.high) push!(Q, v.high) end
    end
  end
  return false
end
export lit_val

"Pretty print a conjunction of literals BDD."
function print_conjunction(α::Bdd; out::Bool = true)::Union{String, Nothing}
  s = ""
  Q = Bdd[α]
  while !isempty(Q)
    v = pop!(Q)
    if !is_term(v)
      s *= is_⊥(v.low) ? string(v.index) : "¬$(v.index)"
      if !is_term(v.low) push!(Q, v.low) end
      if !is_term(v.high) push!(Q, v.high) end
      if !isempty(Q) s *= " ∧ " end
    end
  end
  if out
    println(s)
    return
  end
  return s
end
export print_conjunction

"""Pretty print BDD as a normal form (CNF or DNF).

Caution: may exponentially explode.

Alternatively, pretty prints using the given glyphs (default `∧`, `∨` and `¬`).

```@example
ϕ = (1 ∧ ¬2) ∨ (2 ∧ 3)
print_nf(α; out = false)
```

```@example
ϕ = (1 ∧ ¬2) ∨ (2 ∧ 3)
print_nf(α; out = false, which = "dnf", glyphs = ['+', '*', '-'])
```
"""
function print_nf(α::Bdd; out::Bool = true, which = "cnf", glyphs = nothing)::Union{String, Nothing}
  @assert which == "cnf" || which == "dnf" "BDD.print_nf: Pretty prints only as 'cnf' or 'dnf'."
  if isnothing(glyphs)
    and_glyph, or_glyph, neg_glyph = '∧', '∨', '¬'
  else
    and_glyph, or_glyph, neg_glyph = glyphs[1], glyphs[2], glyphs[3]
  end
  if which == "cnf"
    outer_glyph, inner_glyph, neg_glyph = and_glyph, or_glyph, neg_glyph
  else
    outer_glyph, inner_glyph, neg_glyph = or_glyph, and_glyph, neg_glyph
  end
  s = ""
  C = normal_form(α, which == "dnf")
  m = length(C)
  for (j, c) ∈ enumerate(C)
    n = length(c)
    s *= '('
    for (i, x) ∈ enumerate(c)
      if i == n
        s *= x > 0 ? string(x) : "$neg_glyph$(-x)"
      else
        s *= x > 0 ? "$x $inner_glyph " : "$neg_glyph$(-x) $inner_glyph "
      end
    end
    s *= j == m ? ')' : ") $outer_glyph "
  end
  if out
    println(s)
    return
  end
  return s
end
export print_nf

"Computes a mapping of the parents of each node."
function map_parents(α::Bdd)::Dict{Bdd, Vector{Tuple{Bdd, Bool}}}
  Pa = Dict{Bdd, Vector{Tuple{Bdd, Bool}}}(α => Vector{Tuple{Bdd, Bool}}())
  V = BitSet()
  Q = Bdd[α]
  while !isempty(Q)
    v = popfirst!(Q)
    if !is_term(v)
      l, h = v.low, v.high
      p, q = l.id, h.id
      if q ∉ V Pa[h] = Vector{Tuple{Bdd, Bool}}(); push!(Q, h); push!(V, q) end
      if p ∉ V Pa[l] = Vector{Tuple{Bdd, Bool}}(); push!(Q, l); push!(V, p) end
      push!(Pa[h], (v, true))
      push!(Pa[l], (v, false))
    end
  end
  return Pa
end

"""Runs a BFS on the mapping of parents, starting from either a ⊤ (true) or ⊥ (false) in order to
find the corresponding CNF or DNF encoding."""
function normal_form(α::Bdd, which::Bool)::Vector{Vector{Int}}
  # Which: DNF = true; CNF = false.
  clauses = Vector{Vector{Int}}()
  normal_form_step(which ? ⊤ : ⊥, map_parents(α), clauses, Vector{Int}(), which)
  return clauses
end

function normal_form_step(α::Bdd, Pa::Dict{Bdd, Vector{Tuple{Bdd, Bool}}},
                              all::Vector{Vector{Int}}, path::Vector{Int}, which::Bool)
  pa = Pa[α]
  for (p, e) ∈ pa
    root = isempty(Pa[p])
    x = which == e ? p.index : -p.index
    P = push!(copy(path), x)
    if root push!(all, P)
    else normal_form_step(p, Pa, all, P, which) end
  end
  nothing
end

"Save as CNF. Use the `save` function instead."
function save_cnf(α::Bdd, filename::String; kwargs...)
  S = scopeset(α)
  open(filename, "w"; kwargs...) do out
    write(out, "c CNF file format. See https://people.sc.fsu.edu/~jburkardt/data/cnf/cnf.html.\n")
    cnf = normal_form(α, false)
    write(out, "p cnf $(length(S)) $(length(cnf))\n")
    for C ∈ cnf
      clause = ""
      for x ∈ C
        clause *= "$x "
      end
      clause *= "0\n"
      write(out, clause)
    end
  end
  nothing
end

"Save as DNF. Use the `save` function instead."
function save_dnf(α::Bdd, filename::String; kwargs...)
  S = scopeset(α)
  open(filename, "w"; kwargs...) do out
    write(out, "c DNF file format. See http://gauss.ececs.uc.edu/sbsat_user_manual/node58.html.\n")
    dnf = normal_form(α, true)
    write(out, "p dnf $(length(S)) $(length(dnf))\n")
    for C ∈ dnf
      clause = ""
      for x ∈ C
        clause *= "$x "
      end
      clause *= "0\n"
      write(out, clause)
    end
  end
  nothing
end

"Save as BDD. Use the `save_bdd` function instead."
function save_bdd_(α::Bdd, filename::String; kwargs...)
  open(filename, "w"; kwargs...) do out
    write(out, "c BDD file format. See https://www.ime.usp.br/~renatolg/bdd/bdd.html.\n")
    foreach(function(ϕ::Bdd)
              if is_term(ϕ)
                write(out, "1 $(ϕ.id) $(ϕ.value ? 1 : 0)\n")
              else
                write(out, "0 $(ϕ.id) $(ϕ.index) $(ϕ.low.id) $(ϕ.high.id)\n")
              end
            end, postorder(α))
  end
  nothing
end

"""Saves a BDD as a file.

Supported file formats:
 - CNF (`.cnf`);
 - DNF (`.dnf`);
 - BDD (`.bdd`).

To save as any of these file formats, simply set the filename with the desired extension.

Keyword arguments are passed down to the `open` function.
"""
function save_bdd(α::Bdd, filename::String; kwargs...)
  @assert length(filename) > 4 "BDD.save: Filename must contain name and valid extension!"
  funcs = Dict{String, Function}(".cnf" => save_cnf, ".dnf" => save_dnf, ".bdd" => save_bdd_)
  ext = filename[end-3:end]
  @assert haskey(funcs, ext) "BDD.save: Not a valid extension!"
  funcs[ext](α, filename; kwargs...)
  nothing
end
export save_bdd

"Loads a CNF as a BDD. Use `load` instead."
function load_cnf(::Type{Bdd}, filename::String; kwargs...)::Bdd
  ϕ = ⊤
  open(filename, "r"; kwargs...) do input
    for line ∈ eachline(input)
      l = lstrip(line)
      if l[1] == 'c' || l[1] == 'p' continue end
      C = reduce(∨, map(x -> parse(Int, x), split(l[begin:end-1])))
      ϕ = ϕ ∧ C
    end
  end
  return ϕ
end

"Loads a CNF as a BDD. Use `load` instead."
function load_dnf(::Type{Bdd}, filename::String; kwargs...)::Bdd
  ϕ = ⊥
  open(filename, "r"; kwargs...) do input
    for line ∈ eachline(input)
      l = lstrip(line)
      if l[1] == 'c' || l[1] == 'p' continue end
      C = reduce(∧, map(x -> parse(Int, x), split(l[begin:end-1])))
      ϕ = ϕ ∨ C
    end
  end
  return ϕ
end

"Loads a BDD from a file. Use `load` instead."
function load_bdd(::Type{Bdd}, filename::String; kwargs...)::Bdd
  N = Dict{Int, Bdd}()
  r = -1
  open(filename, "r"; kwargs...) do input
    for line ∈ eachline(input)
      l = lstrip(line)
      if l[1] == 'c' continue end
      P = map(x -> parse(Int, x), split(l))
      if first(P) == 1
        id, value = P[2], P[3]
        if !haskey(N, id) N[id] = Bdd(value == 1) end
      else
        id, index, l_id, h_id = P[2], P[3], P[4], P[5]
        if !haskey(N, id) N[id] = Bdd(index, id, N[l_id], N[h_id]) end
        r = id
      end
    end
  end
  return N[r]
end

"""Loads a BDD from given file.

Supported file formats:
 - CNF (`.cnf`);
 - DNF (`.dnf`);
 - BDD (`.bdd`).

To load as any of these file formats, simply set the filename with the desired extension.

Keyword arguments are passed down to the `open` function.
"""
function load(::Type{Bdd}, filename::String; kwargs...)::Bdd
  @assert length(filename) > 4 "BDD.save: Filename must contain name and valid extension!"
  funcs = Dict{String, Function}(".cnf" => load_cnf, ".dnf" => load_dnf, ".bdd" => load_bdd)
  ext = filename[end-3:end]
  @assert haskey(funcs, ext) "BDD.save: Not a valid extension!"
  return funcs[ext](Bdd, filename; kwargs...)
end
export load

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
function checkAM(neq, nlag, nlead, h, b)

#   substitute the reduced form, b, into the structural
#   equations, h, to confirm that the reduced form solves
#   the model.

# Append negative identity matrix to b

b = hcat(b, -Matrix(I, neq, neq))

# Define indexes into the lagged part (minus)
# and the current and lead part (plus) of h

minus = 1:(neq * nlag)
plus  = (neq * nlag + 1):(neq * (nlag + 1 + nlead))

# Initialize q

q = zeros(neq * (nlead + 1), neq * (nlag + 1 + nlead))

# Stuff b into the upper left-hand block of q

(rb,cb) = size(b)
q[(1:rb), (1:cb)] = b[(1:rb), (1:cb)]

# Copy b into the (nlead) row blocks of q, shifting right
# by neq once more for each block (this produces a coefficient matrix that
# solves for x(t), x(t+1), ..., x(t+nlead)) in terms of x(t-nlag),...,x(t-1).

for i in (1:nlead)
   rows = i * neq .+ (1:neq)
   q[rows,:] = shiftRight!( q[(rows .- neq), :], neq )
end

# Premultiply the left block of q by the negative inverse
# of the right block of q (use equation solver to avoid
# actually computing the inverse).

q[:,minus] =  -q[:,plus] \ q[:,minus]

# Define the solution error

error = h[:,minus] +  h[:,plus] * q[:,minus]

# Take maximum absolute value of the vec of error (error(:))

err = maximum( abs.( error ) );

# space
print("Maximum absolute error: ")
println(err)

# Resize q, removing the additional neq columns and rows

rows = 1:(neq * nlead);
cols = 1:(neq * (nlag + nlead))
q = q[rows,cols];

    return q, err

end


text:
function frankwolfe(network::AbstractNetwork,
              trips::AbstractMatrix{T},
              costfn::Function;
              basedon=:link,
              errtol=1e-4) where {T<:Real}
    # initialize
    flows = allornothing(network, trips, costfn; basedon=basedon)

    # start iteration
    err = 1.
    while err > errtol
        ## find target solution
        linkcosts, = costfn(flows)
        shortflows = allornothing(network, trips, linkcosts; basedon=basedon)

        ## find stepsize
        Δflows = shortflows - flows
        μ = bisection((μ) -> sum(costfn(flows + μ * Δflows).costs .* Δflows))

        ## calculate new link travel times
        flows += μ * Δflows

        ## calculate error
        err = sum(flows .* linkcosts)/sum(shortflows .* linkcosts) - 1.
    end
    return (flows, err)
end


text:
nt = 1000
seed = 99
m = ABCDQR_simula(nt,seed)
y = m["y"]
u = m["u"]

############################### starting point
wi = collect(range(0.1, length=2, stop=3))
Ai,Bi,Ci,Di,Qi,Ri,Si,m1i,P1i = ABCDQR_em_starting(y,u,wi,1,true)

#######################
max_iter=10
tol=1e-6
txo=false
print("Testing ABCDQR_em: ")
A,B,C,D,Q,R,m1,P1,loglik,aic = ABCDQR_em(y,u,Ai,Bi,Ci,Di,Qi,Ri,m1i,P1i,max_iter,tol,txo)
if sum(A.^2) > 1e-6
    println("OK")
end


text:
using DifferentialEquations

function parameterized_lorenz(du,u,p,t)
  x,y,z = u
  σ,ρ,β = p
  du[1] = dx = σ*(y-x)
  du[2] = dy = x*(ρ-z) - y
  du[3] = dz = x*y - β*z
end

using ParameterizedFunctions
g = @ode_def begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ ρ β

# using LinearAlgebra
# f(du,u,p,t) = mul!(du,A,u)

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(g,u0,tspan,p)

using StaticArrays, DifferentialEquations
A  = @SMatrix [ 1.0  0.0 0.0 -5.0
                4.0 -2.0 4.0 -3.0
               -4.0  0.0 0.0  1.0
                5.0 -2.0 2.0  3.0]
u0 = @SMatrix rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)
using Plots;
plot(sol)

# const κ = 1; const tau = 2 * pi
const out = zeros(2) # Define a cache variable
function delay_mathieu_model(du,u,h,p,t)
  δ, ϵ, b, κ = p
  h(out, p, t-tau)
  du[1] = - κ *u[1] - (δ+ϵ*cos(t)) *u[2] -b * out[2]
  du[2] = u[1]
end

p = [1.0,2.0,1,0.1]
tau=2*pi
T=2*pi
lags = [tau]
h(out, p, t) = (out.=1.0)
# h(out, p, t, deriv::Type{Val{i}})
tspan = (0.0,T)
u0 = [1.0,1.0]
prob = DDEProblem(delay_mathieu_model,h,tspan,p; constant_lags=lags)

alg = MethodOfSteps(Tsit5())
# alg = MethodOfSteps(BS3())
# alg = MethodOfSteps(Vern6())
sol = solve(prob,alg)#,reltol=1e-8,abstol=1e-8)

using Plots
plot(sol,vars=(0,2))

# plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
#      xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
# plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")
sol[2][3]
plot(sol,vars=(0,2))