text:
int_rules_9_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 9.2*Derivative*integration*rules =#
    @apply_utils Antiderivative(((Derivative(~n))(~f))(~x), ~x) => ((Derivative(~n - 1))(~f))(~x) <-- FreeQ([~f, ~n], ~x)
    @apply_utils Antiderivative((~(c') * (~F) ^ (~(a') + ~(b') * ~x)) ^ ~(p') * ((Derivative(~n))(~f))(~x), ~x) => (~c * (~F) ^ (~a + ~b * ~x)) ^ ~p * ((Derivative(~n - 1))(~f))(~x) - ~b * ~p * log(~F) * Antiderivative((~c * (~F) ^ (~a + ~b * ~x)) ^ ~p * ((Derivative(~n - 1))(~f))(~x), ~x) <-- FreeQ([~a, ~b, ~c, ~f, ~F, ~p], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(c') * (~F) ^ (~(a') + ~(b') * ~x)) ^ ~(p') * ((Derivative(~n))(~f))(~x), ~x) => ((~c * (~F) ^ (~a + ~b * ~x)) ^ ~p * ((Derivative(~n))(~f))(~x)) / (~b * ~p * log(~F)) - (1 / (~b * ~p * log(~F))) * Antiderivative((~c * (~F) ^ (~a + ~b * ~x)) ^ ~p * ((Derivative(~n + 1))(~f))(~x), ~x) <-- FreeQ([~a, ~b, ~c, ~f, ~F, ~p], ~x) && ILtQ(~n, 0)
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) * ((Derivative(~n))(~f))(~x), ~x) => sin(~a + ~b * ~x) * ((Derivative(~n - 1))(~f))(~x) - ~b * Antiderivative(cos(~a + ~b * ~x) * ((Derivative(~n - 1))(~f))(~x), ~x) <-- FreeQ([~a, ~b, ~f], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) * ((Derivative(~n))(~f))(~x), ~x) => cos(~a + ~b * ~x) * ((Derivative(~n - 1))(~f))(~x) + ~b * Antiderivative(sin(~a + ~b * ~x) * ((Derivative(~n - 1))(~f))(~x), ~x) <-- FreeQ([~a, ~b, ~f], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) * ((Derivative(~n))(~f))(~x), ~x) => (-(cos(~a + ~b * ~x)) * ((Derivative(~n))(~f))(~x)) / ~b + (1 / ~b) * Antiderivative(cos(~a + ~b * ~x) * ((Derivative(~n + 1))(~f))(~x), ~x) <-- FreeQ([~a, ~b, ~f], ~x) && ILtQ(~n, 0)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) * ((Derivative(~n))(~f))(~x), ~x) => (sin(~a + ~b * ~x) * ((Derivative(~n))(~f))(~x)) / ~b - (1 / ~b) * Antiderivative(sin(~a + ~b * ~x) * ((Derivative(~n + 1))(~f))(~x), ~x) <-- FreeQ([~a, ~b, ~f], ~x) && ILtQ(~n, 0)
    @apply_utils Antiderivative(~u * ((Derivative(~n))(~f))(~x), ~x) => Subst(Antiderivative(SimplifyIntegrand(SubstFor(((Derivative(~n - 1))(~f))(~x), ~u, ~x), ~x), ~x), ~x, ((Derivative(~n - 1))(~f))(~x)) <-- FreeQ([~f, ~n], ~x) && FunctionOfQ(((Derivative(~n - 1))(~f))(~x), ~u, ~x)
    @apply_utils Antiderivative(~u * (~(a') * ((Derivative(1))(~f))(~x) * (~g)(~x) + ~(a') * (~f)(~x) * ((Derivative(1))(~g))(~x)), ~x) => ~a * Subst(Antiderivative(SimplifyIntegrand(SubstFor((~f)(~x) * (~g)(~x), ~u, ~x), ~x), ~x), ~x, (~f)(~x) * (~g)(~x)) <-- FreeQ([~a, ~f, ~g], ~x) && FunctionOfQ((~f)(~x) * (~g)(~x), ~u, ~x)
    @apply_utils Antiderivative(~u * (~(a') * ((Derivative(~m))(~f))(~x) * (~g)(~x) + ~(a') * ((Derivative(~m1))(~f))(~x) * ((Derivative(1))(~g))(~x)), ~x) => ~a * Subst(Antiderivative(SimplifyIntegrand(SubstFor(((Derivative(~m - 1))(~f))(~x) * (~g)(~x), ~u, ~x), ~x), ~x), ~x, ((Derivative(~m - 1))(~f))(~x) * (~g)(~x)) <-- FreeQ([~a, ~f, ~g, ~m], ~x) && (EqQ(~m1, ~m - 1) && FunctionOfQ(((Derivative(~m - 1))(~f))(~x) * (~g)(~x), ~u, ~x))
    @apply_utils Antiderivative(~u * (~(a') * ((Derivative(~m))(~f))(~x) * ((Derivative(~n1))(~g))(~x) + ~(a') * ((Derivative(~m1))(~f))(~x) * ((Derivative(~n))(~g))(~x)), ~x) => ~a * Subst(Antiderivative(SimplifyIntegrand(SubstFor(((Derivative(~m - 1))(~f))(~x) * ((Derivative(~n - 1))(~g))(~x), ~u, ~x), ~x), ~x), ~x, ((Derivative(~m - 1))(~f))(~x) * ((Derivative(~n - 1))(~g))(~x)) <-- FreeQ([~a, ~f, ~g, ~m, ~n], ~x) && (EqQ(~m1, ~m - 1) && (EqQ(~n1, ~n - 1) && FunctionOfQ(((Derivative(~m - 1))(~f))(~x) * ((Derivative(~n - 1))(~g))(~x), ~u, ~x)))
    @apply_utils Antiderivative(~u * (~f)(~x) ^ ~(p') * (~(a') * ((Derivative(1))(~f))(~x) * (~g)(~x) + ~(b') * (~f)(~x) * ((Derivative(1))(~g))(~x)), ~x) => ~b * Subst(Antiderivative(SimplifyIntegrand(SubstFor((~f)(~x) ^ (~p + 1) * (~g)(~x), ~u, ~x), ~x), ~x), ~x, (~f)(~x) ^ (~p + 1) * (~g)(~x)) <-- FreeQ([~a, ~b, ~f, ~g, ~p], ~x) && (EqQ(~a, ~b * (~p + 1)) && FunctionOfQ((~f)(~x) ^ (~p + 1) * (~g)(~x), ~u, ~x))
    @apply_utils Antiderivative(~u * ((Derivative(~m1))(~f))(~x) ^ ~(p') * (~(a') * ((Derivative(~m))(~f))(~x) * (~g)(~x) + ~(b') * ((Derivative(~m1))(~f))(~x) * ((Derivative(1))(~g))(~x)), ~x) => ~b * Subst(Antiderivative(SimplifyIntegrand(SubstFor(((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * (~g)(~x), ~u, ~x), ~x), ~x), ~x, ((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * (~g)(~x)) <-- FreeQ([~a, ~b, ~f, ~g, ~m, ~p], ~x) && (EqQ(~m1, ~m - 1) && (EqQ(~a, ~b * (~p + 1)) && FunctionOfQ(((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * (~g)(~x), ~u, ~x)))
    @apply_utils Antiderivative(~u * (~g)(~x) ^ ~(q') * (~(a') * ((Derivative(~m))(~f))(~x) * (~g)(~x) + ~(b') * ((Derivative(~m1))(~f))(~x) * ((Derivative(1))(~g))(~x)), ~x) => ~a * Subst(Antiderivative(SimplifyIntegrand(SubstFor(((Derivative(~m - 1))(~f))(~x) * (~g)(~x) ^ (~q + 1), ~u, ~x), ~x), ~x), ~x, ((Derivative(~m - 1))(~f))(~x) * (~g)(~x) ^ (~q + 1)) <-- FreeQ([~a, ~b, ~f, ~g, ~m, ~q], ~x) && (EqQ(~m1, ~m - 1) && (EqQ(~a * (~q + 1), ~b) && FunctionOfQ(((Derivative(~m - 1))(~f))(~x) * (~g)(~x) ^ (~q + 1), ~u, ~x)))
    @apply_utils Antiderivative(~u * ((Derivative(~m1))(~f))(~x) ^ ~(p') * (~(a') * ((Derivative(~m))(~f))(~x) * ((Derivative(~n1))(~g))(~x) + ~(b') * ((Derivative(~m1))(~f))(~x) * ((Derivative(~n))(~g))(~x)), ~x) => ~b * Subst(Antiderivative(SimplifyIntegrand(SubstFor(((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * ((Derivative(~n - 1))(~g))(~x), ~u, ~x), ~x), ~x), ~x, ((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * ((Derivative(~n - 1))(~g))(~x)) <-- FreeQ([~a, ~b, ~f, ~g, ~m, ~n, ~p], ~x) && (EqQ(~m1, ~m - 1) && (EqQ(~n1, ~n - 1) && (EqQ(~a, ~b * (~p + 1)) && FunctionOfQ(((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * ((Derivative(~n - 1))(~g))(~x), ~u, ~x))))
    @apply_utils Antiderivative(~u * (~f)(~x) ^ ~(p') * (~g)(~x) ^ ~(q') * (~(a') * ((Derivative(1))(~f))(~x) * (~g)(~x) + ~(b') * (~f)(~x) * ((Derivative(1))(~g))(~x)), ~x) => (~a / (~p + 1)) * Subst(Antiderivative(SimplifyIntegrand(SubstFor((~f)(~x) ^ (~p + 1) * (~g)(~x) ^ (~q + 1), ~u, ~x), ~x), ~x), ~x, (~f)(~x) ^ (~p + 1) * (~g)(~x) ^ (~q + 1)) <-- FreeQ([~a, ~b, ~f, ~g, ~p, ~q], ~x) && (EqQ(~a * (~q + 1), ~b * (~p + 1)) && FunctionOfQ((~f)(~x) ^ (~p + 1) * (~g)(~x) ^ (~q + 1), ~u, ~x))
    @apply_utils Antiderivative(~u * ((Derivative(~m1))(~f))(~x) ^ ~(p') * (~g)(~x) ^ ~(q') * (~(a') * ((Derivative(~m))(~f))(~x) * (~g)(~x) + ~(b') * ((Derivative(~m1))(~f))(~x) * ((Derivative(1))(~g))(~x)), ~x) => (~a / (~p + 1)) * Subst(Antiderivative(SimplifyIntegrand(SubstFor(((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * (~g)(~x) ^ (~q + 1), ~u, ~x), ~x), ~x), ~x, ((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * (~g)(~x) ^ (~q + 1)) <-- FreeQ([~a, ~b, ~f, ~g, ~m, ~p, ~q], ~x) && (EqQ(~m1, ~m - 1) && (EqQ(~a * (~q + 1), ~b * (~p + 1)) && FunctionOfQ(((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * (~g)(~x) ^ (~q + 1), ~u, ~x)))
    @apply_utils Antiderivative(~u * ((Derivative(~m1))(~f))(~x) ^ ~(p') * ((Derivative(~n1))(~g))(~x) ^ ~(q') * (~(a') * ((Derivative(~m))(~f))(~x) * ((Derivative(~n1))(~g))(~x) + ~(b') * ((Derivative(~m1))(~f))(~x) * ((Derivative(~n))(~g))(~x)), ~x) => (~a / (~p + 1)) * Subst(Antiderivative(SimplifyIntegrand(SubstFor(((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * ((Derivative(~n - 1))(~g))(~x) ^ (~q + 1), ~u, ~x), ~x), ~x), ~x, ((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * ((Derivative(~n - 1))(~g))(~x) ^ (~q + 1)) <-- FreeQ([~a, ~b, ~f, ~g, ~m, ~n, ~p, ~q], ~x) && (EqQ(~m1, ~m - 1) && (EqQ(~n1, ~n - 1) && (EqQ(~a * (~q + 1), ~b * (~p + 1)) && FunctionOfQ(((Derivative(~m - 1))(~f))(~x) ^ (~p + 1) * ((Derivative(~n - 1))(~g))(~x) ^ (~q + 1), ~u, ~x))))
    @apply_utils Antiderivative(((~f)')(~x) * (~g)(~x) + (~f)(~x) * ((~g)')(~x), ~x) => (~f)(~x) * (~g)(~x) <-- FreeQ([~f, ~g], ~x)
    @apply_utils Antiderivative((((~f)')(~x) * (~g)(~x) - (~f)(~x) * ((~g)')(~x)) / (~g)(~x) ^ 2, ~x) => (~f)(~x) / (~g)(~x) <-- FreeQ([~f, ~g], ~x)
    @apply_utils Antiderivative((((~f)')(~x) * (~g)(~x) - (~f)(~x) * ((~g)')(~x)) / ((~f)(~x) * (~g)(~x)), ~x) => log((~f)(~x) / (~g)(~x)) <-- FreeQ([~f, ~g], ~x)
end

max_stars_repo_path: src/rules/9 Miscellaneous/9.2 Derivative integration rules.jl


text:
mutable struct CheapTrickOption
    q1::Cdouble
    f0_floor::Cdouble
    fft_size::Cint  # (1, 2, ..., 12)

    # allow uninitialized creation
    CheapTrickOption() = new()
end

function cheaptrick(x::Vector{Float64}, fs::Int, temporal_positions::Vector{Float64}, f0::Vector{Float64})
    cheaptrick_option = Ref(CheapTrickOption())
    @ccall libworld.InitializeCheapTrickOption(
        fs::Cint,
        cheaptrick_option::Ptr{CheapTrickOption})::Cvoid

    fft_size = @ccall libworld.GetFFTSizeForCheapTrick(
        fs::Cint,
        cheaptrick_option::Ptr{CheapTrickOption})::Cint

    f0_floor = @ccall libworld.GetF0FloorForCheapTrick(fs::Cint, fft_size::Cint)::Cdouble

    nfreqs = fft_size ÷ 2 + 1
    nbins = length(f0)

    spectrogram = Matrix{Float64}(undef, nfreqs, nbins) # transposed for C

    spectrogram_pointer_vector = [Ref(spectrogram, offset)
        for offset in 1:nfreqs:length(spectrogram)]

    @ccall libworld.CheapTrick(
        x::Ptr{Cdouble},
        length(x)::Cint,
        fs::Cint,
        temporal_positions::Ptr{Cdouble},
        f0::Ptr{Cdouble},
        length(f0)::Cint,
        cheaptrick_option::Ptr{CheapTrickOption},
        spectrogram_pointer_vector::Ptr{Ptr{Cdouble}}
    )::Cvoid

    spectrogram
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 3e0c149a-3350-11eb-0b17-7bea5d1a28d4
begin
	DAY = 1
	md"# 2020 Day $DAY"
end

# ╔═╡ 18080fe2-339b-11eb-34f0-f393638506fc
md"## Part 1"

# ╔═╡ 4e084a4e-339b-11eb-3ffc-41118215f226
md"*<Insert Description>*"

# ╔═╡ 24e1613c-339b-11eb-2976-59884053d1ee
input = read("../data/day_$DAY.txt", String)

# ╔═╡ df1d1be2-36c6-11eb-07c6-493a612dd182
test_input = """
some texts
"""

# ╔═╡ 5f797e22-339d-11eb-1ff3-dda90f946838
md"## Part 2"

# ╔═╡ e3177eae-339e-11eb-3f9f-1f7230976832
md"*<Insert Description>*"

# ╔═╡ Cell order:
# ╠═3e0c149a-3350-11eb-0b17-7bea5d1a28d4
# ╟─18080fe2-339b-11eb-34f0-f393638506fc
# ╠═4e084a4e-339b-11eb-3ffc-41118215f226
# ╠═24e1613c-339b-11eb-2976-59884053d1ee
# ╠═df1d1be2-36c6-11eb-07c6-493a612dd182
# ╟─5f797e22-339d-11eb-1ff3-dda90f946838
# ╟─e3177eae-339e-11eb-3f9f-1f7230976832

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
#=##############################################################################################
Filename: solve_rotor.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: define an `Action` object to solve a CCBlade rotor
=###############################################################################################

"""
    solve_rotor(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

Action function.

# Arguments:

* `aircraft::Aircraft` : `Aircraft` system object
* `parameters<:Parameters` `Parameters` struct
* `freestream::Freestream` : `Freestream` object
* `environment::Environment` `Environment` object
* `step_range::AbstractArray` : array of times for which the simulation is run
* `stepi::Int` : index of the current step
* `step_symbol::String` : defines the step, e.g. `alpha` or `time`

NOTE: `omegas` is commanded while `Js` is calculated

`parameters <: Parameters` requires the following elements:

* `omegas::Vector{Float64}` : a vector of rotational speeds in rad/s at the current step (commanded)
* `Js::Array{Float64,2}` : each i,jth element is the advance ratio of the ith rotor at the jth step (calculated)
* `Ts::Array{Float64,2}` : each i,jth element is the thrust of the ith rotor at the jth step
* `Qs::Array{Float64,2}` : each i,jth element is the torque of the ith rotor at the jth step
* `Ps::Array{Float64,2}`: [i,j]th element is the power of the ith rotor at the jth step
* `us::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the axial induced velocity at ith step of the jth rotor at the kth radial section
* `vs::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the swirl induced velocity at ith step of the jth rotor at the kth radial section

"""
function solve_rotor(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

    # extract parameters
    omegas = parameters.omegas[stepi]
    Js = parameters.Js[:,stepi]
    Ts = parameters.Ts[:,stepi]
    Qs = parameters.Qs[:,stepi]
    Ps = parameters.Ps[:,stepi]
    us = parameters.us[stepi]
    vs = parameters.vs[stepi]

    solve_rotors!(Js, Ts, Qs, Ps, us, vs, aircraft.rotor_system, omegas, freestream, environment)

    return false
end

"""
    solve_rotor(aircraft, step_range)

Method returns initialized elements required for the `parameters <: Parameters` struct during simulation.

# Arguments:

* `aircraft::Aircraft`: system to be simulated
* `step_range::AbstractArray`: defines each step of the simulation

# Returns:

* `omegas::Array{Float64,2}`: [i,j]th elemnt is the rotational speeds in rad/s of the ith rotor at the jth step
* `Js::Array{Float64,2}`: [i,j]th element is the advance ratio of the ith rotor at the jth step
* `Ts::Array{Float64,2}`: [i,j]th element is the thrust of the ith rotor at the jth step
* `Qs::Array{Float64,2}`: [i,j]th element is the torque of the ith rotor at the jth step
* `Ps::Array{Float64,2}`: [i,j]th element is the power of the ith rotor at the jth step
* `us::Vector{Vector{Vector{Float64}}}`: each [i][j][k]th element is the axial induced velocity at ith step of the jth rotor at the kth radial section
* `vs::Vector{Vector{Vector{Float64}}}`: each [i][j][k]th element is the swirl induced velocity at ith step of the jth rotor at the kth radial section

"""
function solve_rotor(aircraft, step_range)

    nrotors = length(aircraft.rotor_system.index) # number of rotors
    omegas = ones(Float64, nrotors, length(step_range))
    Js = zeros(nrotors, length(step_range))
    Ts = zeros(nrotors, length(step_range))
    Qs = zeros(nrotors, length(step_range))
    Ps = zeros(nrotors, length(step_range))
    us = [[zeros(length(aircraft.rotor_system.rlists[i])) for i in aircraft.rotor_system.index] for _ in 1:length(step_range)]
    vs = deepcopy(us)

    return omegas, Js, Ts, Qs, Ps, us, vs
end

# # prepare solutionkeys and solutioninits for this particular system and simulation
# solutionkeys = [
#     "thrust",
#     "torque",
#     "efficiency",
#     "u",
#     "v"
# ]

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using Plots

"""
`plot_phichain(outputfile::String, burnin::Int64, thin::Int64)`

Plots the phi values at each iteration resulting from pmdi output
## Input
- `outputFile::String` a string referring to the location on disk of pmdi output
- `burnin::Int64` an integer of the number of initial iterations to discard as burn-in
- `thin::Int64` an integer for the rate at which to thin output, `thin = 2` discards
every second iteration

## Output
Outputs a line plot of phi values resulting from pmdi output
"""

function get_nclust(outputFile::String, burnin::Int64 = 0, thin::Int64 = 1)
    outputNames = split(readline(outputFile), ',')
    K = sum(map(x -> occursin(r"MassParameter", x), outputNames))
    output = readdlm(outputFile, ',', header = false, skipstart = burnin + 1)
    hyperCols = Int64(K * (K - 1) / 2 + K + 1 + (K == 1))
    output = Matrix{Int64}(output[1:thin:end, (hyperCols + 1):end])
    n_obs = Int64((size(output, 2)) / K)
    dataNames = unique(map(x -> split(x, "_")[1], outputNames[(hyperCols + 1):end]))

    plot_matrix = Matrix{Int64}(undef, size(output, 1), K)

    for k in 1:K
        start_col = (1:n_obs:size(output, 2))[k]
        end_col = (n_obs:n_obs:size(output, 2))[k]
        for i in 1:size(output, 1)
            plot_matrix[i, k] = length(unique(output[i, start_col:end_col]))
        end
    end
    return plot_matrix, dataNames, K
end

function plot_nclust_hist(outputFile::String, burnin::Int64 = 0, thin::Int64 = 1)
    plot_matrix, dataNames, K = get_nclust(outputFile, burnin, thin)

    clust_range = (minimum(plot_matrix), maximum(plot_matrix))
    Plots.histogram(plot_matrix,
                    bins = clust_range[1]:clust_range[2],
                    layout = K,
                    c = :viridis,
                    legend = false,
                    title = [dataNames[i] for j = 1:1, i = 1:K],
                    titlefont = Plots.font(family = "serif", pointsize = 12))
end

function plot_nclust_chain(outputFile::String, burnin::Int64 = 0, thin::Int64 = 1)
    plot_matrix, dataNames, K = get_nclust(outputFile, burnin, thin)
    clust_range = (minimum(plot_matrix), maximum(plot_matrix))
    Plots.plot(plot_matrix,
                    legend = false,
                    layout = K,
                    c = :viridis,
                    title = [dataNames[i] for j = 1:1, i = 1:K],
                    left_margin = 20px,
                    xticks = false)
end


text:
@auto_hash_equals mutable struct Transition
    id::Int
end

mutable struct Place
    id::Int
    marking::Int
end

mutable struct PetriNet
    places::Array{Place}
    transitions::Array{Transition}
    arcs::Array{Union{Tuple{Place, Transition}, Tuple{Transition, Place}}}
end

function fire!(transition::Transition, model::PetriNet)
    if is_enabled(transition, model)
        update_markings!(transition, model)
    else
        print("Transition is not enabled\nDid not fire")
    end
    return model
end

function is_enabled(transition::Transition, model::PetriNet)
    predecessors = get_predecessors(transition, model)
    return sum([place.marking > 0 for place in predecessors]) == length(predecessors)
end

function get_predecessors(transition::Transition, model::PetriNet)
    return [arc[1] for arc in model.arcs if arc[2] == transition]
end

function get_successors(transition::Transition, model::PetriNet)
    return [arc[2] for arc in model.arcs if arc[1] == transition]
end

function update_markings!(transition::Transition, model::PetriNet)
    for predecessor in get_predecessors(transition, model)
        predecessor.marking -= 1
    end
    for successor in get_successors(transition, model)
        successor.marking += 1
    end
    return model
end

function add_place!(model::PetriNet, marking::Int)
    new_place_id = length(model.places) + 1
    push!(model.places, Place(new_place_id, marking))
    return model
end

function add_transition!(model::PetriNet)
    new_transition_id = length(model.transitions) + 1
    push!(model.transitions, Transition(new_transition_id))
    return model
end

function add_arc!(model::PetriNet, from::Place, to::Transition)
    push!(model.arcs, (from, to))
    return model
end

function add_arc!(model::PetriNet, from::Transition, to::Place)
    push!(model.arcs, (from, to))
    return model
end

mutable struct WorkFlowNet <: PetriNet
    source::Int
    sink::Int
    # places::Array{Place}
    # transitions::Array{Transition}
    # arcs::Array{Union{Tuple{Place, Transition}, Tuple{Transition, Place}}}
end

# light weight version test
# places = [Place(1, 2), Place(2, 0)]
# transitions = [Transition(1)]
# arcs = [(places[1], transitions[1]), (transitions[1], places[2])]
# model = PetriNet(places, transitions, arcs)
# fire(model.transitions[1], model)
# model
# add_place!(model, 3)
# add_transition!(model)
# add_arc!(model, model.places[1], model.transitions[2])
# add_arc!(model, model.transitions[2], model.places[2])
# print(model)

# ----- MORE IDEAS ----- #
# reachability graph
# boundedness
# safeness
# deadlocks
# liveness


text:
## @lc code=start
using LeetCode
using DataStructures

function minimum_operations(nums::Vector{Int})::Int
    function get_max_and_premax_data(dict::DefaultDict)::Tuple{Int, Int, Int}
        if length(dict) > 1
            (a, b), (c,) = sort!([(freq, num) for (num, freq) in dict], by=i -> i[1], rev=true)
        else
            (a, b), = sort!([(freq, num) for (num, freq) in dict], by=i -> i[1], rev=true)
            c = 0
        end
        return a, b, c
    end
    n = length(nums)
    n == 1 && return 0
    odd, even = DefaultDict{Int,Int}(0), DefaultDict{Int,Int}(0)
    for (i, num) in enumerate(nums)
        if i % 2 == 1
            odd[num] += 1
        else
            even[num] += 1
        end
    end
    odd_max_freq, odd_max_num, odd_premax_freq = get_max_and_premax_data(odd)
    even_max_freq, even_max_num, even_premax_freq = get_max_and_premax_data(even)
    if even_max_num != odd_max_num
        return n - even_max_freq - odd_max_freq
    else
        return n - max(even_max_freq + odd_premax_freq, odd_max_freq + even_premax_freq)
    end
end
## @lc code=end

## @lc test=start
@testset "2170.minimum-operations-to-make-the-array-alternating.jl" begin
    @test minimum_operations([3,1,3,2,4,3]) == 3
    @test minimum_operations([2,2,2,2,2]) == 2
    @test minimum_operations([2,2,2,2,1,1,1,1]) == 4
    @test minimum_operations([3,1,3,2,4,3]) == 3
    @test minimum_operations([1,2,2,2,2]) == 2
end
## @lc test=end

## @lc info=start
# link: [solution to 2170](https://leetcode-cn.com/problems/minimum-operations-to-make-the-array-alternating/solution/6005-zui-da-pin-ci-he-ci-da-pin-ci-julia-6oz2/)
## @lc info=end


text:
"""
    lag(v::AbstractArray, n = 1; default = missing)

Return a `ShiftedArray` object, with underlying data `v`. The second argument gives the amount
to shift in each dimension. If it is an integer, it is assumed to refer to the first dimension.
`default` specifies a default value when you are out of bounds.

# Examples

```jldoctest lag
julia> v = [1, 3, 5, 4];

julia> lag(v)
4-element ShiftedArray{Int64,Missing,1,Array{Int64,1}}:
  missing
 1
 3
 5

julia> w = 1:2:9
1:2:9

julia> s = lag(w, 2)
5-element ShiftedArray{Int64,Missing,1,StepRange{Int64,Int64}}:
  missing
  missing
 1
 3
 5

julia> copy(s)
5-element Array{Union{Missing, Int64},1}:
  missing
  missing
 1
 3
 5

julia> v = reshape(1:16, 4, 4);

julia> s = lag(v, (0, 2))
4×4 ShiftedArray{Int64,Missing,2,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}}:
 missing  missing  1  5
 missing  missing  2  6
 missing  missing  3  7
 missing  missing  4  8
```
"""
lag(v::AbstractArray, n = 1; default = missing) = ShiftedArray(v, n; default = default)

"""
    lead(v::AbstractArray, n = 1; default = missing)

Return a `ShiftedArray` object, with underlying data `v`. The second argument gives the amount
to shift negatively in each dimension. If it is an integer, it is assumed to refer
to the first dimension. `default` specifies a default value when you are out of bounds.

# Examples

```jldoctest lead
julia> v = [1, 3, 5, 4];

julia> lead(v)
4-element ShiftedArray{Int64,Missing,1,Array{Int64,1}}:
 3
 5
 4
  missing

julia> w = 1:2:9
1:2:9

julia> s = lead(w, 2)
5-element ShiftedArray{Int64,Missing,1,StepRange{Int64,Int64}}:
 5
 7
 9
  missing
  missing

julia> copy(s)
5-element Array{Union{Missing, Int64},1}:
 5
 7
 9
  missing
  missing

julia> v = reshape(1:16, 4, 4);

julia> s = lead(v, (0, 2))
4×4 ShiftedArray{Int64,Missing,2,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}}:
  9  13  missing  missing
 10  14  missing  missing
 11  15  missing  missing
 12  16  missing  missing
```
"""
lead(v::AbstractArray, n = 1; default = missing) = ShiftedArray(v, map(-, n); default = default)

"""
    lag(v::AbstractVector, times::AbstractVector, period = oneunit(eltype(times)); default = missing) -> Vector

Shifts with respect to a times given in the vector `times`.  The third variable `period` gives the period by which to shift.
`default` specifies a default value when the shifted time is not in `times`.
Elements in `times` must all be distinct.

# Examples

```jldoctest lead
julia> v = [1, 3, 5, 4];
julia> times = [1990, 1992, 1993];
julia> lag(v, times)
3-element Array{Union{Missing, Int64},1}:
  missing
  missing
 3
julia> using Dates
julia> times = [Date(1990, 1, 1), Date(1990, 1, 3), Date(1990, 1, 4)]
julia> lag(v, times, Day(1))
3-element Array{Union{Missing, Int64},1}:
 missing
 missing
3
"""
function lag(v::AbstractVector, times::AbstractVector, period = oneunit(eltype(times)); default = missing)
   	# Code follows the function indexin
    inds = keys(times)
    timesdict = Dict{eltype(times),eltype(inds)}()
    for (val, ind) in zip(times, inds)
         out = get!(timesdict, val, ind)
         out != ind && error("Times must be distinct")
    end
    return Union{eltype(v), typeof(default)}[
         (i = get(timesdict, x - period, nothing); i !== nothing ? v[i] : default) for x in times
    ]
end

"""
    lead(v::AbstractVector, times::AbstractVector, period = oneunit(eltype(times)); default = missing) -> Vector

Shifts with respect to a vector of times `times`. The third variable `period` gives the period by which to shift.
`default` specifies a default value when the shifted time is not in `times`.
Elements in `times` must all be distinct.

# Examples

```jldoctest lead
julia> v = [1, 3, 5, 4];
julia> times = [1990, 1992, 1993];
julia> lead(v, times, 1)
3-element Array{Union{Missing, Int64},1}:
  missing
  5
  missing
"""
function lead(v::AbstractVector, times::AbstractVector, period = oneunit(eltype(times)); default = missing)
    lag(v, times, -period; default = default)
end


text:
let
	conc_CH4_preind 	= 710.0		# ppb
	conc_N2O_preind		= 273.0		# ppb

	cacheterm = 0.47 * log(1 + 2.01e-5*(conc_CH4_preind*conc_N2O_preind)^0.75 +
							+ 5.31e-15*conc_CH4_preind*(conc_CH4_preind*conc_N2O_preind)^1.52)

	global function radiativeforcing!(Concentration, RadiativeForcing)
		# Saturation of CH4 absorption bands minus overlap with N2O
		# IPCC TAR WG1 table 6.2 (p358), or Tanaka eq 2.2.4 & 2.2.5 (sign error in 2.2.4)
		# Can skip the entire 5e-15*() terms, error ~0.1% for 2400 ppb
		RadiativeForcing[:CH4] =
			0.036 * (sqrt(Concentration[:CH4]) - sqrt(conc_CH4_preind)) +
			- 0.47 * log(1 + 2.01e-5*(Concentration[:CH4]*conc_N2O_preind)^0.75 +
							+ 5.31e-15*Concentration[:CH4]*(Concentration[:CH4]*conc_N2O_preind)^1.52) +
			+ cacheterm

		# ditto for N2O
		RadiativeForcing[:N2O] =
			0.12 * (sqrt(Concentration[:N2O]) - sqrt(conc_N2O_preind)) +
			- 0.47 * log(1 + 2.01e-5*(conc_CH4_preind*Concentration[:N2O])^0.75 +
							+ 5.31e-15*conc_CH4_preind*(conc_CH4_preind*Concentration[:N2O])^1.52) +
			+ cacheterm

		# Stratospheric H2O = 15% of radiative forcing from CH4, see IPCC AR4 WG1 sections 2.3.7 and 2.10.3.1
		# Large uncertainty here - it was judged to be 2-5% in TAR (see WG1 section 6.6.4 and Tanaka 2.2.29).
		RadiativeForcing[:H2O] = 0.15 * RadiativeForcing[:CH4]

		# Tropospheric O3 from CH4, Tanaka eq 2.2.21, based on IPCC TAR WG1 tables 4.11 (footnote p269) & 4.9 (p261)
		# Forcing due to ozone precursors (NOx, CO, VOC) is aggregated with RF_other.
		RadiativeForcing[:O3] = 0.042 * 5.0 * log(Concentration[:CH4]/conc_CH4_preind)

		# IPCC TAR WG1 table 6.2 (p358), or Tanaka eq 2.1.73,    5.35 ≈ 3.71/log(2)
		RadiativeForcing[:CO2] = 5.35 * log(Concentration[:CO2]/conc_CO2_preind)
	end

	global function radiativeforcing!(s::ClimateState)
		radiativeforcing!(s.Concentration, s.RadiativeForcing)
	end

	global function radiativeforcing!(t, s::ClimateState, p::ClimateParams, rcp)
		radiativeforcing!(s)
		s.TotalRadiativeForcing = sum(s.RadiativeForcing[g] for g in GAS) +
			interpolate(t, otherforcing[rcp]) + p.aerosolforcingfactor * interpolate(t, aerosolforcing[rcp])
	end

	global function init_concentration_and_forcing!(s::ClimateState, p::ClimateParams)
		@unpack Concentration = s
		RadiativeForcing = Dict{Symbol,Float64}(g => 0.0 for g in GAS)
		TotalRadiativeForcing = 0.0

		Concentration[:CH4] = conc_CH4_preind
		Concentration[:N2O] = conc_N2O_preind

		@pack! s = Concentration, RadiativeForcing, TotalRadiativeForcing
		radiativeforcing!(float(YEARS[1]), s, p, "RCP3PD")
	end
end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
@with_kw struct precipParam
  # --------------------
  # Physical Parameters
  # --------------------
  μ0::Float64  = 1.0  # viscosity
  μ1::Float64  = 1.0  # brinkman viscosity
  α0::Float64  = 0.0  # initial value of thetam/thetaf
  G::Float64   = 1e0  # pressure inflow
  ξϵ::Float64  = 1e0  # regularization of friction coefficient

  PeA::Float64 = 30    # Peclet number for chemical A
  PeB::Float64 = 30    # Peclet number for chemical B
  PeC::Float64 = 1     # Peclet number for chemical C

  ρf::Float64 = 1.0     # mass density of fluid
  ρs::Float64 = 4.1     # mass density of solid

  # computed parameters
  α::Float64 = 1
  h::Float64 = 1

  # ---------------------
  # Chemistry Parameters
  # ---------------------
  a::Int = 1
  b::Int = 1
  c::Int = 1
  Ma::Float64 = 58.693
  Mb::Float64 = 17.008
  Mc::Float64 = 0.0

  # ----------------------------
  # Time-integration parameters
  # ----------------------------
  Δt::Float64        = 5e-2
  Tfinal::Float64    = 80.0
  integrator::Symbol = :cranknicholson

  # ---------------
  # vtk Parameters
  # ---------------
  ψAname::String = "aqueous_A"
  ψBname::String = "aqueous_B"
  ψCname::String = "aqueous_C"
  θsname::String = "solid_vf"
  θfname::String = "fluid_vf"
  printSkipIndex::Int = 10
  printBool::Bool     = true
  SaveFolder::String  = string(homedir(),"/data/6-7-2020-y-dom")
  SaveName::String    = "y-dom"

  # ----------------
  # Mesh Parameters
  # ----------------
  MeshFolder::String="/home/peastham/Dropbox/projects/project-precipitation---numerics/PrecipSim.jl/meshes/y-dom"
  MeshFile::String = "y-dom6.msh"

  # -------------
  # run ID (1-7)
  # -------------
  runID::Int = 1
end


text:
myseek(x)=pipeline(`echo seek $x absolute`,`socat - /tmp/mpvsocket`)

set_gamma(x)=pipeline(`echo set gamma $x`,`socat - /tmp/mpvsocket`)

set_brightness(x)=pipeline(`echo set brightness $x`,`socat - /tmp/mpvsocket`)

mpv_open(x) = readandwrite(`mpv --hr-seek=always --input-ipc-server=/tmp/mpvsocket --quiet --osdlevel=0 $x`)

const pause_cmd = pipeline(`echo 'set pause yes'`,`socat - /tmp/mpvsocket`)


text:
using LinearMixingModels
using Documenter

DocMeta.setdocmeta!(
    LinearMixingModels,
    :DocTestSetup,
    :(using LinearMixingModels, KernelFunctions, FillArrays, LinearAlgebra, AbstractGPs);
    recursive=true,
)

makedocs(;
    modules=[LinearMixingModels],
    authors="Invenia Technical Computing Corporation",
    repo="https://github.com/invenia/LinearMixingModels.jl/blob/{commit}{path}#{line}",
    sitename="LinearMixingModels.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://invenia.github.io/LinearMixingModels.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
    checkdocs=:exports,
    strict=true,
)

deploydocs(;
    repo="github.com/invenia/LinearMixingModels.jl",
)


text:
# Metatheory.options[:printiter] = true
using Metatheory
using Metatheory.EGraphs
using Metatheory.Classic
using Test

or_alg = @theory begin
    ((p ∨ q) ∨ r)       ==  (p ∨ (q ∨ r))
    (p ∨ q)             ==  (q ∨ p)
    (p ∨ p)             =>  p
    (p ∨ true)          =>  true
    (p ∨ false)         =>  p
end

and_alg = @theory begin
    ((p ∧ q) ∧ r)       ==  (p ∧ (q ∧ r))
    (p ∧ q)             ==  (q ∧ p)
    (p ∧ p)             =>  p
    (p ∧ true)          =>  p
    (p ∧ false)         =>  false
end

comb = @theory begin
    # DeMorgan
    ¬(p ∨ q)            ==  (¬p ∧ ¬q)
    ¬(p ∧ q)            ==  (¬p ∨ ¬q)
    # distrib
    (p ∧ (q ∨ r))       ==  ((p ∧ q) ∨ (p ∧ r))
    (p ∨ (q ∧ r))       ==  ((p ∨ q) ∧ (p ∨ r))
    # absorb
    (p ∧ (p ∨ q))       =>  p
    (p ∨ (p ∧ q))       =>  p
    # complement
    (p ∧ (¬p ∨ q))      =>  p ∧ q
    (p ∨ (¬p ∧ q))      =>  p ∨ q
end

negt = @theory begin
    (p ∧ ¬p)            =>  false
    (p ∨ ¬(p))          =>  true
    ¬(¬p)               ==  p
end

impl = @theory begin
    (p == ¬p)           =>  false
    (p == p)            =>  true
    (p == q)            =>  (¬p ∨ q) ∧ (¬q ∨ p)
    (p => q)            =>  (¬p ∨ q)
end

fold = @theory begin
    (p::Bool == q::Bool)   |> (p == q)
    (p::Bool ∨ q::Bool)    |> (p || q)
    (p::Bool => q::Bool)   |> ((p || q) == q)
    (p::Bool ∧ q::Bool)    |> (p && q)
    ¬(p::Bool)             |> (!p)
end

t = or_alg ∪ and_alg ∪ comb ∪ negt ∪ impl ∪ fold


text:
const CACHE_MODE = Symbol(uppercase(@load_preference("cache_mode", "DEFAULT")))
const CACHE_DIR = @load_preference(
    "cache_dir", joinpath(DEPOT_PATH[1], "datadeps", "JuliaConSchedule")
)
const TIMEOUT = parse(Float64, @load_preference("timeout", "5.0"))
const TERMINAL_LINKS = parse(Bool, @load_preference("terminal_links", "false"))
const JULIACON_TIMEZONE = tz"UTC"
const LOCAL_TIMEZONE = begin
    tzstr = @load_preference("local_timezone", "")
    if !istimezone(tzstr)
        @has_preference("local_timezone") &&
            @warn "Preference \"local_timezone\" isn't a valid timezone. Check out TimeZones.timezone_names() for possible options. Falling back to localzone()."
        localzone()
    else
        TimeZone(tzstr)
    end
end

const PRETALX_JSON_URL = "https://pretalx.com/juliacon2021/schedule/export/schedule.json"
const DATA_ARCHIVE_JSON_URL = "https://raw.githubusercontent.com/JuliaCon/JuliaConDataArchive/master/juliacon2021_schedule/schedule.json"
const jcon = Ref{DataFrame}()

function set_cachemode(mode::Symbol)
    @assert mode in (:DEFAULT, :NEVER, :ALWAYS)
    @set_preferences!("cache_mode" => uppercase(string(mode)))
    @info("New cache mode set; restart your Julia session for this change to take effect!")
end
function set_cachedir(path::String)
    @set_preferences!("cache_dir" => string(path))
    @info("New cache dir set; restart your Julia session for this change to take effect!")
end
function set_timeout(val::Union{Float64,Int})
    @assert val > 0
    @set_preferences!("timeout" => string(val))
    @info("New timeout set; restart your Julia session for this change to take effect!")
end
function set_terminallinks(on::Bool)
    @set_preferences!("terminal_links" => string(on))
    @info(
        "New terminal links preference set; restart your Julia session for this change to take effect!"
    )
end
function set_local_timezone(tz::AbstractString)
    if !istimezone(tz)
        throw(ArgumentError("\"$tz\" isn't a valid timezone. Check out TimeZones.timezone_names() for possible options."))
    else
        @set_preferences!("local_timezone" => tz)
        @info(
            "New local timezone set; restart your Julia session for this change to take effect!"
        )
        return nothing
    end
end
function reset_local_timezone()
    tz = string(localzone())
    tz = tz == "Etc/UTC" ? "UTC" : tz # Make CI runners happy....
    set_local_timezone(tz)
end

default_json_url() = DATA_ARCHIVE_JSON_URL
default_now() = TimeZones.now(LOCAL_TIMEZONE)

"""
    debugmode(on::Bool=true)

Simulates that we are live / in the middle of JuliaCon.
"""
function debugmode(on::Bool=true)
    if on
        # @eval JuliaCon default_now() = ZonedDateTime(Dates.DateTime("2021-07-20T18:00:00.000"), tz"MET") # JuliaCon2021 workshops
        # @eval JuliaCon default_now() = ZonedDateTime(Dates.DateTime("2021-07-30T22:00:00.000"), tz"MET") # JuliaCon2021 posters
        @eval JuliaCon function default_now()
            return ZonedDateTime(Dates.DateTime("2021-07-30T21:30:00.000"), tz"MET")
            # return ZonedDateTime(Dates.DateTime("2021-07-24T18:15:00.000"), tz"America/Los_Angeles")
            # return ZonedDateTime(Dates.DateTime("2021-07-24T20:00:00.000"), tz"MET")
        end # JuliaCon2021
    else
        @eval JuliaCon default_now() = TimeZones.now(LOCAL_TIMEZONE)
    end
    return nothing
end

"""
    set_json_source(src::Symbol)

Anticipated input: `:pretalx`, `:github` (JuliaConDataArchive)
"""
function set_json_source(src::Symbol)
    if src == :pretalx
        @eval JuliaCon default_json_url() = PRETALX_JSON_URL
    else
        @eval JuliaCon default_json_url() = DATA_ARCHIVE_JSON_URL
    end
    return nothing
end


text:
# Book Crossing Dataset (http://www.informatik.uni-freiburg.de/~cziegler/BX)
#
# Collected by Cai-Nicolas Ziegler in a 4-week crawl (August / September 2004)
# from the Book-Crossing community with kind permission from Ron Hornbaker, CTO
# of Humankind Systems. Contains 278,858 users (anonymized but with demographic
# information) providing 1,149,780 ratings (explicit / implicit) about 271,379
# books.
#
# Freely available for research use when acknowledged with the following
# reference (further details on the dataset are given in this publication):
# Improving Recommendation Lists Through Topic Diversification,
# Cai-Nicolas Ziegler, Sean M. McNee, Joseph A. Konstan, Georg Lausen;
# Proceedings of the 14th International World Wide Web Conference (WWW '05),
# May 10-14, 2005, Chiba, Japan.
#
# Ratings are 0-10, but 0 means that the user didn't rate the book. There are
# also lots of users who have very few ratings and lots of books that have
# only been rated a few times, so these rating need a lot of cleaning before
# they're usuable. With the default parameters, we exclude any user who has
# rated less than 30 books and any book that's been rated less than 30 times,
# as well as any user who hasn't used at least 6 different ratings out of the
# range 0-10. Finally, we re-scale the ratings to the range 1-6 by mapping 0-5
# to 1 and scaling all of the other ratings down by 4.

function load_book_crossing_dataset(;
                                    user_ratings_threshold=30,
                                    book_ratings_threshold=30,
                                    user_distinct_ratings_threshold=6,
                                    book_distinct_ratings_threshold=1)
    zipfile_name = "BX-CSV-Dump.zip"
    temp_dir = tempdir()
    download_path = joinpath(temp_dir, zipfile_name)
    if !isfile(download_path)
        println("Downloading book ratings data...")
        download("http://www.informatik.uni-freiburg.de/~cziegler/BX/$(zipfile_name)", download_path)
    else
        println("Reusing existing downloaded files...")
    end
    ratings_file_name = joinpath(temp_dir, "BX-Book-Ratings.csv")
    if !isfile(ratings_file_name)
        println("Extracting ratings data...")
        run(`unzip $download_path BX-Book-Ratings.csv -d $temp_dir`)
    end
    books_file_name = joinpath(temp_dir, "BX-Books.csv")
    if !isfile(books_file_name)
        println("Extracting books data...")
        run(`unzip $download_path BX-Books.csv -d $temp_dir`)
    end

    books_count = countlines(books_file_name) - 1
    progress = Progress(books_count, 1, "Loading ISBN-to-book mapping ")
    book_file = open(books_file_name, "r")
    readline(book_file) # skip CSV header
    isbn_to_title = Dict{AbstractString,AbstractString}()
    for line in eachline(book_file)
        delim1 = search(line, "\";\"")
        delim2 = search(line, "\";\"", last(delim1))
        isbn, title = line[2:first(delim1)-1], line[last(delim1)+1:first(delim2)-1]
        isbn_to_title[isbn] = title
        next!(progress)
    end
    close(book_file)

    ratings_count = countlines(ratings_file_name) - 1
    progress = Progress(ratings_count, 1, "Collecting counts of books and users for filtering ")
    user_ratings_count = Dict{AbstractString,Int32}()
    book_ratings_count = Dict{AbstractString,Int32}()
    user_distinct_ratings_set = Dict{AbstractString,Set{AbstractString}}()
    book_distinct_ratings_set = Dict{AbstractString,Set{AbstractString}}()
    ratings_file = open(ratings_file_name, "r")
    readline(ratings_file) # skip CSV header
    for line in eachline(ratings_file)
        try ascii(line) catch; continue end # Some of these lines (~17 out of a million) are non-ASCII garbage.
        delim1 = search(line, "\";\"")
        delim2 = search(line, "\";\"", last(delim1))
        user, isbn, rating = line[2:first(delim1)-1], line[last(delim1)+1:first(delim2)-1], line[last(delim2)+1:end-3]
        user_ratings_count[user] = get(user_ratings_count, user, 0) + 1
        book_ratings_count[isbn] = get(book_ratings_count, isbn, 0) + 1
        if !haskey(user_distinct_ratings_set, user)
            user_distinct_ratings_set[user] = Set{AbstractString}()
        end
        push!(user_distinct_ratings_set[user], rating)
        if !haskey(book_distinct_ratings_set, isbn)
            book_distinct_ratings_set[isbn] = Set{AbstractString}()
        end
        push!(book_distinct_ratings_set[isbn], rating)
        next!(progress)
    end
    close(ratings_file)

    # Filter ratings while translating them into Rating structures.
    progress = Progress(ratings_count, 1, "Loading ratings ")
    ratings = (Rating)[]
    item_to_index = Dict{AbstractString, Int32}()
    user_to_index = Dict{AbstractString, Int32}()
    ratings_file = open(ratings_file_name, "r")
    readline(ratings_file) # skip CSV header
    for line in eachline(ratings_file)
        try ascii(line) catch; continue end # Filter non-ASCII garbage like we did in the loop above.
        delim1 = search(line, "\";\"")
        delim2 = search(line, "\";\"", last(delim1))
        user, isbn, rating = line[2:first(delim1)-1], line[last(delim1)+1:first(delim2)-1], line[last(delim2)+1:end-3]
        if haskey(isbn_to_title, isbn) && # some isbns aren't in BX-Books.csv
           user_ratings_count[user] >= user_ratings_threshold &&
           book_ratings_count[isbn] >= book_ratings_threshold &&
           length(user_distinct_ratings_set[user]) >= user_distinct_ratings_threshold &&
           length(book_distinct_ratings_set[isbn]) >= book_distinct_ratings_threshold
            item = isbn_to_title[isbn]
            rating = max(float(rating) - 4, 1.0)
            user_index = get(user_to_index, user, 0)
            if user_index == 0
                user_index = user_to_index[user] = length(user_to_index) + 1
            end
            item_index = get(item_to_index, item, 0)
            if item_index == 0
                item_index = item_to_index[item] = length(item_to_index) + 1
            end
            push!(ratings, Rating(user_index, item_index, rating))
        end
        next!(progress)
    end
    close(ratings_file)
    training_set, test_set = split_ratings(ratings, 0.10)
    RatingSet(training_set, test_set, user_to_index, item_to_index)
end


text:
export BrickletTemperatureTemperatureCallbackThreshold
struct BrickletTemperatureTemperatureCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletTemperatureIdentity
struct BrickletTemperatureIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletTemperature
"""
Measures ambient temperature with 0.5°C accuracy
"""
mutable struct BrickletTemperature <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletTemperature(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_temperature")
        deviceInternal = package.BrickletTemperature(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_temperature
"""
    $(SIGNATURES)

Returns the temperature of the sensor.

If you want to get the temperature periodically, it is recommended
to use the :cb:`Temperature` callback and set the period with
:func:`Set Temperature Callback Period`.
"""
function get_temperature(device::BrickletTemperature)
    return device.deviceInternal.get_temperature()
end

export set_temperature_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Temperature` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Temperature` callback is only triggered if the temperature has changed
since the last triggering.
"""
function set_temperature_callback_period(device::BrickletTemperature, period)
    device.deviceInternal.set_temperature_callback_period(period)
end

export get_temperature_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Temperature Callback Period`.
"""
function get_temperature_callback_period(device::BrickletTemperature)
    return device.deviceInternal.get_temperature_callback_period()
end

export set_temperature_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Temperature Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the temperature is *outside* the min and max values"
 "'i'",    "Callback is triggered when the temperature is *inside* the min and max values"
 "'<'",    "Callback is triggered when the temperature is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the temperature is greater than the min value (max is ignored)"
"""
function set_temperature_callback_threshold(device::BrickletTemperature, option, min, max)
    device.deviceInternal.set_temperature_callback_threshold(option, min, max)
end

export get_temperature_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Temperature Callback Threshold`.
"""
function get_temperature_callback_threshold(device::BrickletTemperature)
    return device.deviceInternal.get_temperature_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Temperature Reached`

is triggered, if the threshold

* :func:`Set Temperature Callback Threshold`

keeps being reached.
"""
function set_debounce_period(device::BrickletTemperature, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletTemperature)
    return device.deviceInternal.get_debounce_period()
end

export set_i2c_mode
"""
    $(SIGNATURES)

Sets the I2C mode. Possible modes are:

* 0: Fast (400kHz)
* 1: Slow (100kHz)

If you have problems with obvious outliers in the
Temperature Bricklet measurements, they may be caused by EMI issues.
In this case it may be helpful to lower the I2C speed.

It is however not recommended to lower the I2C speed in applications where
a high throughput needs to be achieved.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_i2c_mode(device::BrickletTemperature, mode)
    device.deviceInternal.set_i2c_mode(mode)
end

export get_i2c_mode
"""
    $(SIGNATURES)

Returns the I2C mode as set by :func:`Set I2C Mode`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_i2c_mode(device::BrickletTemperature)
    return device.deviceInternal.get_i2c_mode()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletTemperature)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletTemperature, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
## Overridden behaviours
"""
    WrapperLogType

Abstract wrapper type to override the default backend to which a type would be
logged in TensorBoard.

Every concrete implementation should implement the method `content` for the
concrete type, returning the content.
"""
abstract type WrapperLogType end

"""
    content(x::WrapperLogType)

Returns the wrapped data inside of `x`.
"""
content(x::WrapperLogType) = @error "You should define `content($(typeof(x)))`"

# When we hit a wrapped type, we dispatch based on it's type.
preprocess(name, val::WrapperLogType, data) = push!(data, name=>val)

# When logging to text or console, unwrap the LogType to fix formatting.
Base.show(io::IO, mime::AbstractString, x::WrapperLogType) =
    Base.show(io, mime, content(x))

########## For things going to LogAudio ##############################
"""
    TBAudios(data, samplerate)

Forces elements of Array `data` to be serialized as Audio to TensorBoard.
"""
struct TBAudios <:WrapperLogType
    data::AbstractArray
    samplerate::Real
end
content(x::TBAudios) = x.data
function preprocess(name, val::TBAudios, data)
    for (n, audio) in enumerate(val.data)
        preprocess(name*"/$n", TBAudio(audio, val.samplerate), data)
    end
end

"""
    TBAudio(data, samplerate)

Forces `data` to be serialized as Audio to TensorBoard.
"""
struct TBAudio <: WrapperLogType
    data::AbstractArray
    samplerate::Real
end
content(x::TBAudio) = x.data
preprocess(name, val::TBAudio, data) = push!(data, name=>val)
summary_impl(name, val::TBAudio) = audio_summary(name, val.data, val.samplerate)

########## For things going to LogImage ########################
"""
    TBImages(data, format)

Forces elements of Array `data` to be serialized as an Image to TensorBoard.
"""
struct TBImages <:WrapperLogType
    data::AbstractArray
    format::ImageFormat
end
content(x::TBImages) = x.data
function preprocess(name, val::TBImages, data)
    for (n, img) in enumerate(val.data)
        preprocess(name*"/$n", TBImage(img, val.format), data)
    end
end

"""
    TBImage(data, format)

Forces `data` to be serialized as an Image to TensorBoard.
"""
struct TBImage <: WrapperLogType
    data::AbstractArray
    format::ImageFormat
end
content(x::TBImage) = x.data

function Base.convert(T::Type{PNG}, tb_image::TBImage)
    imgArray = tb_image.data
    format   = tb_image.format

    #unpack RGB, RGBA value to channels using channelview
    imgArray = channelview(imgArray)
    #if data contains integer numbers, scale them to 0-255
    if isa(first(imgArray), Integer)
        imgArray = (imgArray./255)
    end
    #convert all values to `Float64` for uniformity
    imgArray = Float64.(imgArray)
    #scale all values to 0-1
    imgArray = (imgArray./(max(maximum(imgArray), 1)))
    #convert any format to CHW
    imgArray =
        format == L   ? reshape(imgArray, (1, 1, size(imgArray, 1))) :
        format == CL  ? reshape(imgArray, (size(imgArray, 1), 1, size(imgArray, 2))) :
        format == LC  ? reshape(transpose(imgArray), (size(imgArray, 2), 1, size(imgArray, 1))) :
        format == HW  ? reshape(imgArray, (1, size(imgArray, 1), size(imgArray, 2))) :
        format == WH  ? reshape(transpose(imgArray), (1, size(imgArray, 2), size(imgArray, 1))) :
        format == HWC ? permutedims(imgArray, (3, 1, 2)) :
        format == WHC ? permutedims(imgArray, (3, 2, 1)) :
        format == CHW ? imgArray :
        format == CWH ? permutedims(imgArray, (1, 3, 2)) :
        #== else ==# throw("Invalid format")
    channels, height, width = size(imgArray)
    color =
        channels == 1 ? Gray :
        channels == 2 ? GrayA :
        channels == 3 ? RGB :
        channels == 4 ? RGBA :
        #== else ==# throw("Too many channels")
    #if it is a single channel Array, convert it to HW
    if color == Gray
        imgArray = imgArray[1, :, :]
    end
    #convert Array to PNG pass it to image_summary
    img = colorview(color, imgArray)
    return convert(PNG, img)
end

function preprocess(name, val::TBImage, data)
    imgArray = val.data
    format = val.format
    imgArray = channelview(imgArray)
    dims = ndims(imgArray)
    @assert dims == expected_ndims(format)
    obsdim = obs_dim(format)
    if iszero(obsdim)
        preprocess(name, convert(PNG, val), data)
    else
        format = strip_obs[format]
        index = collect("[:"* ",:"^(dims-1) *"]")
        index[2*obsdim] = 'g'
        index = join(index)
        global gimgArray = imgArray
        nth_img = Meta.parse("gimgArray$(index)")
        for n in 1:size(imgArray, obsdim)
            global g = n
            preprocess(name*"/$n", convert(PNG, TBImage(eval(nth_img), format)), data)
        end
    end
    return data
end

########## For things going to LogText ##############################
"""
    TBText(data)

Forces `data` to be serialized as text to TensorBoard.
"""
struct TBText <: WrapperLogType data; end
content(x::TBText) = x.data
summary_impl(name, val::TBText) = text_summary(name, val.data)

"""
    RealArrayWrapperLogType{T,N} <: WrapperLogType

A wrapper type that accepts only arrays and converts complex-valued arrays
into two real value ones
"""
abstract type RealArrayWrapperLogType{T,N} <: WrapperLogType end

########## For things going to LogHistograms ########################
"""
    TBHistogram(data)

Forces `data` to be serialized as an histogram to TensorBoard.
"""
struct TBHistogram{T,N} <: RealArrayWrapperLogType{T,N}
    data::AbstractArray{T,N};
end
content(x::TBHistogram) = x.data
preprocess(name, val::TBHistogram{T,N}, data) where {T<:Complex,N} =
    push!(data, name*"/re"=>TBHistogram(real.(content(val))), name*"/im"=>TBHistogram(imag.(content(val))))
function summary_impl(name, val::TBHistogram)
    hvals = fit(Histogram, val.data)
    return histogram_summary(name, collect(hvals.edges[1]), hvals.weights)
end

"""
    TBVector(data)

Forces `data` to be serialized as a vector in the histogram backend of
TensorBoard.
"""
struct TBVector{T,N} <: RealArrayWrapperLogType{T,N} data::AbstractArray{T,N}; end
content(x::TBVector) = x.data
preprocess(name, val::TBVector{T,N}, data) where {T<:Complex,N} =
    push!(data, name*"/re"=>TBVector(real.(content(val))), name*"/im"=>TBVector(imag.(content(val))))
summary_impl(name, val::TBVector) = histogram_summary(name, collect(0:length(val.data)),
                                                            val.data)


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# urbanplan.jl
#
# An "urban planning" problem.
# Based on
#  http://puzzlor.editme.com/urbanplanning
#############################################################################

using JuMP

function SolveUrban()

    m = Model()

    # x is indexed by row and column
    @variable(m, 0 <= x[1:5,1:5] <= 1, Int)

    # y is indexed by R or C, and the points
    # JuMP allows indexing on arbitrary sets
    rowcol = ["R","C"]
    points = [+5,+4,+3,-3,-4,-5]
    @variable(m, 0 <= y[rowcol,points,i=1:5] <= 1, Int)

    # Objective - combine the positive and negative parts
    @objective(m, Max, sum{
      3*(y["R", 3,i] + y["C", 3,i])
    + 1*(y["R", 4,i] + y["C", 4,i])
    + 1*(y["R", 5,i] + y["C", 5,i])
    - 3*(y["R",-3,i] + y["C",-3,i])
    - 1*(y["R",-4,i] + y["C",-4,i])
    - 1*(y["R",-5,i] + y["C",-5,i]), i=1:5})

    # Constrain the number of residential lots
    @constraint(m, sum{x[i,j], i=1:5, j=1:5} == 12)

    # Add the constraints that link the auxiliary y variables
    # to the x variables
    # Rows
    for i = 1:5
        @constraint(m, y["R", 5,i] <=   1/5*sum{x[i,j], j=1:5}) # sum = 5
        @constraint(m, y["R", 4,i] <=   1/4*sum{x[i,j], j=1:5}) # sum = 4
        @constraint(m, y["R", 3,i] <=   1/3*sum{x[i,j], j=1:5}) # sum = 3
        @constraint(m, y["R",-3,i] >= 1-1/3*sum{x[i,j], j=1:5}) # sum = 2
        @constraint(m, y["R",-4,i] >= 1-1/2*sum{x[i,j], j=1:5}) # sum = 1
        @constraint(m, y["R",-5,i] >= 1-1/1*sum{x[i,j], j=1:5}) # sum = 0
    end
    # Columns
    for j = 1:5
        @constraint(m, y["C", 5,j] <=   1/5*sum{x[i,j], i=1:5}) # sum = 5
        @constraint(m, y["C", 4,j] <=   1/4*sum{x[i,j], i=1:5}) # sum = 4
        @constraint(m, y["C", 3,j] <=   1/3*sum{x[i,j], i=1:5}) # sum = 3
        @constraint(m, y["C",-3,j] >= 1-1/3*sum{x[i,j], i=1:5}) # sum = 2
        @constraint(m, y["C",-4,j] >= 1-1/2*sum{x[i,j], i=1:5}) # sum = 1
        @constraint(m, y["C",-5,j] >= 1-1/1*sum{x[i,j], i=1:5}) # sum = 0
    end

    # Solve it with the default solver (CBC)
    status = solve(m)
    if status == :Infeasible
        error("Solver couldn't find solution!")
    end

    # Print results
    println("Best objective: $(round(getobjectivevalue(m)))")
    println(getvalue(x))
end

SolveUrban()


text:
module Word2Vec

# using Base.Collections      # for priority queue
using Base.Cartesian        # for @nexprs
using Distances
# using NumericExtensions
using DataStructures
using StatsFuns

export LinearClassifier, train_one, WordEmbedding, train, accuracy
export save, restore
export find_nearest_words

include("utils.jl")
include("softmax_classifier.jl")
include("tree.jl")
include("word_stream.jl")
include("train.jl")
include("query.jl")

end # module


text:
dtr = 0.4;                                                                          # replanning time interval [s]
tcalc = 0.2;                                                                        # pre-allocated control computation time [s] (< dtr)

model_dir = normpath(joinpath(@__DIR__, "../../CrowdNav/crowd_nav/data/output_om_sarl_radius_0.4")) # directory of the trained policy
env_config = "env.config"                                                           # environment config file name
policy_config = "policy.config"                                                     # policy config file name
policy_name = "sarl"                                                                # policy name


text:
function interpolate_at_reference_points(
    nodalvalues,
    dofspernode,
    basis,
    refpoints,
    refcellids,
    levelsetsign,
    mesh,
)
    dim, numpts = size(refpoints)
    @assert length(refcellids) == numpts

    interpolatedvals = zeros(dofspernode, numpts)

    for idx = 1:numpts
        cellid = refcellids[idx]
        nodeids = nodal_connectivity(mesh, levelsetsign, cellid)
        celldofs = element_dofs(nodeids, dofspernode)
        cellvals = nodalvalues[celldofs]

        vals = basis(refpoints[:, idx])
        NI = interpolation_matrix(vals, dofspernode)

        interpolatedvals[:, idx] = NI * cellvals
    end

    return interpolatedvals
end

function displacement_at_reference_points(
    nodaldisplacement,
    basis,
    refpoints,
    refcellids,
    levelsetsign,
    mesh,
)

    dim = dimension(mesh)
    return interpolate_at_reference_points(
        nodaldisplacement,
        dim,
        basis,
        refpoints,
        refcellids,
        levelsetsign,
        mesh,
    )
end

function plane_strain_stress(
    celldisp,
    basis,
    stiffness,
    referencepoint,
    jac,
    vectosymmconverter,
)

    dim = length(vectosymmconverter)

    grad = transform_gradient(gradient(basis, referencepoint), jac)
    NK = sum([make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim])

    symmdispgrad = NK * celldisp

    inplanestress = stiffness * symmdispgrad

    return inplanestress
end

function stress_at_reference_points(
    nodaldisplacement,
    basis,
    stiffness,
    referencepoints,
    referencecellids,
    levelsetsign,
    mesh,
)

    dim = dimension(mesh)

    nphase, dim2, numpts = size(referencepoints)
    @assert dim2 == dim
    @assert size(referencecellids) == (nphase, numpts)
    @assert nphase == 2

    vectosymmconverter = vector_to_symmetric_matrix_converter()
    jac = jacobian(mesh)

    stress = zeros(3, numpts)
    row = cell_sign_to_row(levelsetsign)

    for i = 1:numpts
        cellid = referencecellids[row, i]
        nodeids = nodal_connectivity(mesh, levelsetsign, cellid)
        celldofs = element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        referencepoint = referencepoints[row, :, i]

        stress[:, i] .= plane_strain_stress(
            celldisp,
            basis,
            stiffness[levelsetsign],
            referencepoint,
            jac,
            vectosymmconverter,
        )
    end

    return stress
end

function symmetric_displacement_gradient(
    celldisp,
    basis,
    point,
    jac,
    vectosymmconverter,
)
    dim = length(vectosymmconverter)
    grad = transform_gradient(gradient(basis, point), jac)
    NK = sum([make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim])
    symmdispgrad = NK * celldisp
    return symmdispgrad
end

function plane_strain_at_reference_points(
    nodaldisplacement,
    basis,
    referencepoints,
    referencecellids,
    levelsetsign,
    mesh,
)

    @assert levelsetsign == +1 || levelsetsign == -1
    dim = dimension(mesh)
    dim2, numpts = size(referencepoints)
    @assert dim == dim2
    @assert length(referencecellids) == numpts

    strain = zeros(3, numpts)
    vectosymmconverter = vector_to_symmetric_matrix_converter()
    jac = jacobian(mesh)

    for i = 1:numpts
        cellid = referencecellids[i]
        nodeids = nodal_connectivity(mesh, levelsetsign, cellid)
        celldofs = element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        point = referencepoints[:, i]

        strain[:, i] .= symmetric_displacement_gradient(
            celldisp,
            basis,
            point,
            jac,
            vectosymmconverter,
        )
    end
    return strain
end

function parent_strain(
    nodaldisplacement,
    basis,
    referencepoints,
    referencecellids,
    mesh,
)

    strain = plane_strain_at_reference_points(
        nodaldisplacement,
        basis,
        referencepoints,
        referencecellids,
        -1,
        mesh,
    )
    ndofs, numpts = size(strain)
    return vcat(strain, zeros(numpts)')
end

function product_elastic_strain(
    nodaldisplacement,
    basis,
    theta0,
    referencepoints,
    referencecellids,
    mesh,
)

    strain = plane_strain_at_reference_points(
        nodaldisplacement,
        basis,
        referencepoints,
        referencecellids,
        +1,
        mesh,
    )
    ndofs,numpts = size(strain)
    planestrain = vcat(strain,zeros(numpts)')
    planestrain[[1,2,4],:] .-= theta0/3

    return planestrain
end

function parent_stress(symmdispgrad, stiffness)
    lambda, mu = lame_coefficients(stiffness, -1)

    inplanestress = stiffness[-1] * symmdispgrad[1:3, :]
    s33 = lambda * (symmdispgrad[1, :] + symmdispgrad[2, :])

    return vcat(inplanestress, s33')
end

function product_stress(elasticstrain, stiffness, theta0)
    lambda, mu = lame_coefficients(stiffness, +1)

    inplanestress = stiffness[+1] * elasticstrain[1:3, :]
    inplanecorrection = lambda * elasticstrain[4, :]

    inplanestress[1, :] .+= inplanecorrection
    inplanestress[2, :] .+= inplanecorrection

    s33 =
        (lambda + 2mu) * elasticstrain[4, :] +
        lambda * (elasticstrain[1, :] + elasticstrain[2, :])

    return vcat(inplanestress,s33')
end

function strain_energy(stress,strain)
    product = stress .* strain
    return 0.5*vec(sum(product,dims=1))
end

function traction_force(stressvector, normal)
    return [
        stressvector[1] * normal[1] + stressvector[3] * normal[2],
        stressvector[3] * normal[1] + stressvector[2] * normal[2],
    ]
end

function traction_force_at_points(stresses, normals)
    npts = size(stresses)[2]
    @assert size(normals)[2] == npts

    tractionforce = zeros(2, npts)
    for i = 1:npts
        tractionforce[:, i] .= traction_force(stresses[:, i], normals[:, i])
    end
    return tractionforce
end

function traction_component(traction, normal)
    component = traction .* normal
    return vec(sum(component, dims = 1))
end

function pressure_at_points(stress)
    return -1.0 / 3.0 * (stress[1, :] + stress[2, :] + stress[4, :])
end

function collect_cell_quadratures(cellquads, mesh, cellsign, cellids)
    totalnumqps =
        sum([length(cellquads[cellsign, cellid]) for cellid in cellids])
    dim = dimension(mesh)

    referencepoints = zeros(dim, totalnumqps)
    spatialpoints = zeros(dim, totalnumqps)
    referencecellids = zeros(Int, totalnumqps)

    start = 1
    for cellid in cellids
        quad = cellquads[cellsign, cellid]
        numqps = length(quad)

        stop = start + numqps - 1
        cellmap = cell_map(mesh, cellsign, cellid)

        referencepoints[:, start:stop] = points(quad)
        spatialpoints[:, start:stop] = cellmap(points(quad))
        referencecellids[start:stop] = repeat([cellid], numqps)

        start = stop + 1
    end
    return referencepoints, spatialpoints, referencecellids
end

function collect_cell_quadratures(cellquads, mesh, cellsign)
    ncells = number_of_cells(mesh)
    cellsigns = [cell_sign(mesh, cellid) for cellid = 1:ncells]
    cellids = findall(x -> x == cellsign || x == 0, cellsigns)
    collect_cell_quadratures(cellquads, mesh, cellsign, cellids)
end

function stress_inner_product(stress)
    return (stress[1, :] .* stress[1, :]) +
           (stress[2, :] .* stress[2, :]) +
           (stress[4, :] .* stress[4, :]) +
           2.0 * (stress[3, :] .* stress[3, :])
end

function dilatation(strain)
    return strain[1, :] + strain[2, :] + strain[4, :]
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: account
using AWS.Compat
using AWS.UUIDs

"""
    delete_alternate_contact(alternate_contact_type)
    delete_alternate_contact(alternate_contact_type, params::Dict{String,<:Any})

Deletes the specified alternate contact from an Amazon Web Services account. For complete
details about how to use the alternate contact operations, see Access or updating the
alternate contacts.

# Arguments
- `alternate_contact_type`: Specifies which of the alternate contacts to delete.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AccountId"`: Specifies the 12 digit account ID number of the Amazon Web Services
  account that you want to access or modify with this operation. If you do not specify this
  parameter, it defaults to the Amazon Web Services account of the identity used to call the
  operation. To use this parameter, the caller must be an identity in the organization's
  management account or a delegated administrator account, and the specified account ID must
  be a member account in the same organization. The organization must have all features
  enabled, and the organization must have trusted access enabled for the Account Management
  service, and optionally a delegated admin account assigned.  The management account can't
  specify its own AccountId; it must call the operation in standalone context by not
  including the AccountId parameter.  To call this operation on an account that is not a
  member of an organization, then don't specify this parameter, and call the operation using
  an identity belonging to the account whose contacts you wish to retrieve or modify.
"""
function delete_alternate_contact(
    AlternateContactType; aws_config::AbstractAWSConfig=global_aws_config()
)
    return account(
        "POST",
        "/deleteAlternateContact",
        Dict{String,Any}("AlternateContactType" => AlternateContactType);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_alternate_contact(
    AlternateContactType,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return account(
        "POST",
        "/deleteAlternateContact",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("AlternateContactType" => AlternateContactType),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_alternate_contact(alternate_contact_type)
    get_alternate_contact(alternate_contact_type, params::Dict{String,<:Any})

Retrieves the specified alternate contact attached to an Amazon Web Services account. For
complete details about how to use the alternate contact operations, see Access or updating
the alternate contacts.

# Arguments
- `alternate_contact_type`: Specifies which alternate contact you want to retrieve.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AccountId"`: Specifies the 12 digit account ID number of the Amazon Web Services
  account that you want to access or modify with this operation. If you do not specify this
  parameter, it defaults to the Amazon Web Services account of the identity used to call the
  operation. To use this parameter, the caller must be an identity in the organization's
  management account or a delegated administrator account, and the specified account ID must
  be a member account in the same organization. The organization must have all features
  enabled, and the organization must have trusted access enabled for the Account Management
  service, and optionally a delegated admin account assigned.  The management account can't
  specify its own AccountId; it must call the operation in standalone context by not
  including the AccountId parameter.  To call this operation on an account that is not a
  member of an organization, then don't specify this parameter, and call the operation using
  an identity belonging to the account whose contacts you wish to retrieve or modify.
"""
function get_alternate_contact(
    AlternateContactType; aws_config::AbstractAWSConfig=global_aws_config()
)
    return account(
        "POST",
        "/getAlternateContact",
        Dict{String,Any}("AlternateContactType" => AlternateContactType);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_alternate_contact(
    AlternateContactType,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return account(
        "POST",
        "/getAlternateContact",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("AlternateContactType" => AlternateContactType),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_alternate_contact(alternate_contact_type, email_address, name, phone_number, title)
    put_alternate_contact(alternate_contact_type, email_address, name, phone_number, title, params::Dict{String,<:Any})

Modifies the specified alternate contact attached to an Amazon Web Services account. For
complete details about how to use the alternate contact operations, see Access or updating
the alternate contacts.

# Arguments
- `alternate_contact_type`: Specifies which alternate contact you want to create or update.
- `email_address`: Specifies an email address for the alternate contact.
- `name`: Specifies a name for the alternate contact.
- `phone_number`: Specifies a phone number for the alternate contact.
- `title`: Specifies a title for the alternate contact.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AccountId"`: Specifies the 12 digit account ID number of the Amazon Web Services
  account that you want to access or modify with this operation. If you do not specify this
  parameter, it defaults to the Amazon Web Services account of the identity used to call the
  operation. To use this parameter, the caller must be an identity in the organization's
  management account or a delegated administrator account, and the specified account ID must
  be a member account in the same organization. The organization must have all features
  enabled, and the organization must have trusted access enabled for the Account Management
  service, and optionally a delegated admin account assigned.  The management account can't
  specify its own AccountId; it must call the operation in standalone context by not
  including the AccountId parameter.  To call this operation on an account that is not a
  member of an organization, then don't specify this parameter, and call the operation using
  an identity belonging to the account whose contacts you wish to retrieve or modify.
"""
function put_alternate_contact(
    AlternateContactType,
    EmailAddress,
    Name,
    PhoneNumber,
    Title;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return account(
        "POST",
        "/putAlternateContact",
        Dict{String,Any}(
            "AlternateContactType" => AlternateContactType,
            "EmailAddress" => EmailAddress,
            "Name" => Name,
            "PhoneNumber" => PhoneNumber,
            "Title" => Title,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_alternate_contact(
    AlternateContactType,
    EmailAddress,
    Name,
    PhoneNumber,
    Title,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return account(
        "POST",
        "/putAlternateContact",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "AlternateContactType" => AlternateContactType,
                    "EmailAddress" => EmailAddress,
                    "Name" => Name,
                    "PhoneNumber" => PhoneNumber,
                    "Title" => Title,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
export global_norm,
    clip_by_global_norm!,
    find_all_max,
    discount_rewards,
    discount_rewards!,
    discount_rewards_reduced,
    generalized_advantage_estimation,
    generalized_advantage_estimation!,
    flatten_batch,
    orthogonal

using FillArrays: Trues

#####
# Zygote
#####

global_norm(gs, ps) = sqrt(sum(mapreduce(x -> x^2, +, gs[p]) for p in ps))

function clip_by_global_norm!(gs, ps, clip_norm::Float32)
    gn = global_norm(gs, ps)
    if clip_norm <= gn
        for p in ps
            gs[p] .*= clip_norm / max(clip_norm, gn)
        end
    end
    gn
end

#####
# Flux
#####

# https://github.com/FluxML/Flux.jl/pull/1171/
# https://www.tensorflow.org/api_docs/python/tf/keras/initializers/Orthogonal
function orthogonal_matrix(rng::AbstractRNG, nrow, ncol)
    shape = reverse(minmax(nrow, ncol))
    a = randn(rng, Float32, shape)
    q, r = qr(a)
    q = Matrix(q) * diagm(sign.(diag(r)))
    nrow < ncol ? permutedims(q) : q
end

function orthogonal(rng::AbstractRNG, d1, rest_dims...)
    m = orthogonal_matrix(rng, d1, *(rest_dims...))
    reshape(m, d1, rest_dims...)
end

orthogonal(dims...) = orthogonal(Random.GLOBAL_RNG, dims...)
orthogonal(rng::AbstractRNG) = (dims...) -> orthogonal(rng, dims...)

#####
# MLUtils
#####

"""
    flatten_batch(x::AbstractArray)

Merge the last two dimension.

# Example

```julia-repl
julia> x = reshape(1:12, 2, 2, 3)
2×2×3 reshape(::UnitRange{Int64}, 2, 2, 3) with eltype Int64:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

[:, :, 3] =
  9  11
 10  12

julia> flatten_batch(x)
2×6 reshape(::UnitRange{Int64}, 2, 6) with eltype Int64:
 1  3  5  7   9  11
 2  4  6  8  10  12
```
"""
flatten_batch(x::AbstractArray) = reshape(x, size(x)[1:end-2]..., :)

#####
# RLUtils
#####

function find_all_max(x)
    v = maximum(x)
    v, findall(==(v), x)
end

find_all_max(x, mask::Trues) = find_all_max(x)

function find_all_max(x, mask::AbstractVector{Bool})
    v = maximum(view(x, mask))
    v, [k for (m, k) in zip(mask, keys(x)) if m && x[k] == v]
end

# !!! watch https://github.com/JuliaLang/julia/pull/35316#issuecomment-622629895
# Base.findmax(f, domain) = mapfoldl(x -> (f(x), x), _rf_findmax, domain)
# _rf_findmax((fm, m), (fx, x)) = isless(fm, fx) ? (fx, x) : (fm, m)

# !!! type piracy
Base.findmax(A::AbstractVector{T}, mask::AbstractVector{Bool}) where {T} =
    findmax(ifelse.(mask, A, typemin(T)))

Base.findmax(A::AbstractVector, mask::Trues) = findmax(A)

const VectorOrMatrix = Union{AbstractMatrix,AbstractVector}

"""
    discount_rewards(rewards::VectorOrMatrix, γ::Number;kwargs...)

Calculate the gain started from the current step with discount rate of `γ`.
`rewards` can be a matrix.

# Keyword arguments

- `dims=:`, if `rewards` is a `Matrix`, then `dims` can only be `1` or `2`.
- `terminal=nothing`, specify if each reward follows by a terminal. `nothing` means the game is not terminated yet. If `terminal` is provided, then the size must be the same with `rewards`.
- `init=nothing`, `init` can be used to provide the the reward estimation of the last state.

# Example
"""
function discount_rewards(rewards::VectorOrMatrix, γ::T; kwargs...) where {T<:Number}
    res = similar(rewards, promote_type(eltype(rewards), T))
    discount_rewards!(res, rewards, γ; kwargs...)
    res
end

discount_rewards!(new_rewards, rewards, γ; terminal = nothing, init = nothing, dims = :) =
    _discount_rewards!(new_rewards, rewards, γ, terminal, init, dims)

function _discount_rewards!(
    new_rewards::AbstractMatrix,
    rewards::AbstractMatrix,
    γ,
    terminal::Nothing,
    init::Nothing,
    dims::Int,
)
    dims = ndims(rewards) - dims + 1
    for (r′, r) in zip(eachslice(new_rewards, dims = dims), eachslice(rewards, dims = dims))
        _discount_rewards!(r′, r, γ, nothing, nothing)
    end
end

function _discount_rewards!(
    new_rewards::AbstractMatrix,
    rewards::AbstractMatrix,
    γ,
    terminal::Nothing,
    init,
    dims::Int,
)
    dims = ndims(rewards) - dims + 1
    for (i, (r′, r)) in
        enumerate(zip(eachslice(new_rewards, dims = dims), eachslice(rewards, dims = dims)))
        _discount_rewards!(r′, r, γ, nothing, init[i])
    end
end

function _discount_rewards!(
    new_rewards::AbstractMatrix,
    rewards::AbstractMatrix,
    γ,
    terminal,
    init::Nothing,
    dims::Int,
)
    dims = ndims(rewards) - dims + 1
    for (r′, r, t) in zip(
        eachslice(new_rewards, dims = dims),
        eachslice(rewards, dims = dims),
        eachslice(terminal, dims = dims),
    )
        _discount_rewards!(r′, r, γ, t, nothing)
    end
end

function _discount_rewards!(
    new_rewards::AbstractMatrix,
    rewards::AbstractMatrix,
    γ,
    terminal,
    init,
    dims::Int,
)
    dims = ndims(rewards) - dims + 1
    for (i, (r′, r, t)) in enumerate(
        zip(
            eachslice(new_rewards, dims = dims),
            eachslice(rewards, dims = dims),
            eachslice(terminal, dims = dims),
        ),
    )
        _discount_rewards!(r′, r, γ, t, init[i])
    end
end

_discount_rewards!(
    new_rewards::AbstractVector,
    rewards::AbstractVector,
    γ,
    terminal,
    init,
    dims::Colon,
) = _discount_rewards!(new_rewards, rewards, γ, terminal, init)

"assuming rewards and new_rewards are Vector"
_discount_rewards!(new_rewards, rewards, γ, terminal, init::Nothing) =
    _discount_rewards!(new_rewards, rewards, γ, terminal, zero(eltype(new_rewards)))

function _discount_rewards!(new_rewards, rewards, γ, terminal, init)
    gain = init
    for i in length(rewards):-1:1
        is_continue = isnothing(terminal) ? true : (!terminal[i])
        gain = rewards[i] + γ * gain * is_continue
        new_rewards[i] = gain
    end
    new_rewards
end

discount_rewards_reduced(rewards::AbstractVector, γ; terminal = nothing, init = nothing) =
    _discount_rewards_reduced(rewards, γ, terminal, init)

function discount_rewards_reduced(
    rewards::AbstractMatrix,
    γ::T;
    terminal = nothing,
    init = nothing,
    dims,
) where {T<:Number}
    dims = ndims(rewards) - dims + 1
    res = Array{promote_type(eltype(rewards), T)}(undef, size(rewards, dims))
    _discount_rewards_reduced!(res, rewards, γ, terminal, init, dims)
    res
end

_discount_rewards_reduced(rewards, γ, terminal, init::Nothing) =
    _discount_rewards_reduced(rewards, γ, terminal, zero(eltype(rewards)))

function _discount_rewards_reduced(rewards, γ, terminal, init)
    gain = init
    for i in length(rewards):-1:1
        is_continue = isnothing(terminal) ? true : (!terminal[i])
        gain = rewards[i] + γ * gain * is_continue
    end
    gain
end

discount_rewards_reduced!(
    reduced_rewards::AbstractVector,
    rewards::AbstractMatrix,
    γ;
    terminal = nothing,
    init = nothing,
    dims,
) = _discount_rewards_reduced!(reduced_rewards, rewards, γ, terminal, init, dims)

function _discount_rewards_reduced!(
    reduced_rewards,
    rewards,
    γ,
    terminal::Nothing,
    init::Nothing,
    dims::Int,
)
    for (i, r) in enumerate(eachslice(rewards, dims = dims))
        reduced_rewards[i] = _discount_rewards_reduced(r, γ, nothing, nothing)
    end
end

function _discount_rewards_reduced!(
    reduced_rewards,
    rewards,
    γ,
    terminal::Nothing,
    init,
    dims::Int,
)
    for (i, r) in enumerate(eachslice(rewards, dims = dims))
        reduced_rewards[i] = _discount_rewards_reduced(r, γ, nothing, init[i])
    end
end

function _discount_rewards_reduced!(
    reduced_rewards,
    rewards,
    γ,
    terminal,
    init::Nothing,
    dims::Int,
)
    for (i, (r, t)) in
        enumerate(zip(eachslice(rewards, dims = dims), eachslice(terminal, dims = dims)))
        reduced_rewards[i] = _discount_rewards_reduced(r, γ, t, nothing)
    end
end

function _discount_rewards_reduced!(reduced_rewards, rewards, γ, terminal, init, dims::Int)
    for (i, (r, t)) in
        enumerate(zip(eachslice(rewards, dims = dims), eachslice(terminal, dims = dims)))
        reduced_rewards[i] = _discount_rewards_reduced(r, γ, t, init[i])
    end
end

"""
    generalized_advantage_estimation(rewards::VectorOrMatrix, values::VectorOrMatrix, γ::Number, λ::Number;kwargs...)

Calculate the generalized advantage estimate started from the current step with discount rate of `γ` and a lambda for GAE-Lambda of 'λ'.
`rewards` and 'values' can be a matrix.

# Keyword arguments

- `dims=:`, if `rewards` is a `Matrix`, then `dims` can only be `1` or `2`.
- `terminal=nothing`, specify if each reward follows by a terminal. `nothing` means the game is not terminated yet. If `terminal` is provided, then the size must be the same with `rewards`.

# Example
"""
function generalized_advantage_estimation(
    rewards::VectorOrMatrix,
    values::VectorOrMatrix,
    γ::T,
    λ::T;
    kwargs...,
) where {T<:Number}
    res = similar(rewards, promote_type(eltype(rewards), T))
    generalized_advantage_estimation!(res, rewards, values, γ, λ; kwargs...)
    res
end

generalized_advantage_estimation!(
    advantages,
    rewards,
    values,
    γ,
    λ;
    terminal = nothing,
    dims = :,
) = _generalized_advantage_estimation!(advantages, rewards, values, γ, λ, terminal, dims)

function _generalized_advantage_estimation!(
    advantages::AbstractMatrix,
    rewards::AbstractMatrix,
    values::AbstractMatrix,
    γ,
    λ,
    terminal::Nothing,
    dims::Int,
)
    dims = ndims(rewards) - dims + 1
    for (r′, r, v) in zip(
        eachslice(advantages, dims = dims),
        eachslice(rewards, dims = dims),
        eachslice(values, dims = dims),
    )
        _generalized_advantage_estimation!(r′, r, v, γ, λ, nothing)
    end
end

function _generalized_advantage_estimation!(
    advantages::AbstractMatrix,
    rewards::AbstractMatrix,
    values::AbstractMatrix,
    γ,
    λ,
    terminal,
    dims::Int,
)
    dims = ndims(rewards) - dims + 1
    for (r′, r, v, t) in zip(
        eachslice(advantages, dims = dims),
        eachslice(rewards, dims = dims),
        eachslice(values, dims = dims),
        eachslice(terminal, dims = dims),
    )
        _generalized_advantage_estimation!(r′, r, v, γ, λ, t)
    end
end

_generalized_advantage_estimation!(
    advantages::AbstractVector,
    rewards::AbstractVector,
    values::AbstractVector,
    γ,
    λ,
    terminal,
    dims::Colon,
) = _generalized_advantage_estimation!(advantages, rewards, values, γ, λ, terminal)

"assuming rewards and advantages are Vector"
function _generalized_advantage_estimation!(advantages, rewards, values, γ, λ, terminal)
    gae = 0
    for i in length(rewards):-1:1
        is_continue = isnothing(terminal) ? true : (!terminal[i])
        delta = rewards[i] + γ * values[i+1] * is_continue - values[i]
        gae = delta + γ * λ * is_continue * gae
        advantages[i] = gae
    end
    advantages
end


text:
struct GreedyHillClimbing <: GraphSearchStrategy
    cache::ScoreComponentCache
    max_n_parents::Int
    prior::DirichletPrior

    function GreedyHillClimbing(
        cache::ScoreComponentCache;
        max_n_parents::Int=3,
        prior::DirichletPrior=UniformPrior(1.0),
        )

        new(cache, max_n_parents, prior)
    end
end
function Distributions.fit(::Type{DiscreteCPD},
    data::DataFrame,
    target::NodeName,
    prior::DirichletPrior;
    ncategories::Int = infer_number_of_instantiations(data[target]),
    )

    prior_counts = get(prior, ncategories)
    for v in data[target]
        prior_counts[v] += 1.0
    end

    d = Categorical{Float64}(prior_counts ./ sum(prior_counts))
    CategoricalCPD(target, NodeName[], Int[], Categorical{Float64}[d])
end
function Distributions.fit(::Type{DiscreteCPD},
    data::DataFrame,
    target::NodeName,
    parents::NodeNames,
    prior::DirichletPrior;
    parental_ncategories::Vector{Int} = map!(p->infer_number_of_instantiations(data[p]), Array{Int}(length(parents)), parents),
    target_ncategories::Int = infer_number_of_instantiations(data[target]),
    )

    # with parents

    if isempty(parents)
        return fit(DiscreteCPD, data, target, prior, ncategories=target_ncategories)
    end

    nparents = length(parents)
    dims = [1:parental_ncategories[i] for i in 1:nparents]
    distributions = Array{Categorical{Float64}}(prod(parental_ncategories))
    for (q, parent_instantiation) in enumerate(product(dims...))

        prior_counts = get(prior, target_ncategories)
        for i in 1 : nrow(data)
            if all(j->data[i,parents[j]]==parent_instantiation[j], 1:nparents) # parental instantiation matches
                prior_counts[data[i, target]] += 1.0
            end
        end
        distributions[q] = Categorical{Float64}(prior_counts ./ sum(prior_counts))
    end

    CategoricalCPD(target, parents, parental_ncategories, distributions)
end
function Distributions.fit(::Type{DiscreteBayesNet}, data::DataFrame, params::GreedyHillClimbing;
    ncategories::Vector{Int} = map!(i->infer_number_of_instantiations(data[i]), Array{Int}(ncol(data)), 1:ncol(data)),
    )

    n = ncol(data)
    parent_list = map!(i->Int[], Array{Vector{Int}}(n), 1:n)
    datamat = convert(Matrix{Int}, data)'
    score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)

    while true
        best_diff = 0.0
        best_parent_list = parent_list
        for i in 1:n

            # 1) add an edge (j->i)
            if length(parent_list[i]) < params.max_n_parents
                for j in deleteat!(collect(1:n), parent_list[i])
                    if adding_edge_preserves_acyclicity(parent_list, j, i)
                        new_parents = sort!(push!(copy(parent_list[i]), j))
                        new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
                        if new_component_score - score_components[i] > best_diff
                            best_diff = new_component_score - score_components[i]
                            best_parent_list = deepcopy(parent_list)
                            best_parent_list[i] = new_parents
                        end
                    end
                end
            end

            # 2) remove an edge
            for (idx, j) in enumerate(parent_list[i])

                new_parents = deleteat!(copy(parent_list[i]), idx)
                new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
                if new_component_score - score_components[i] > best_diff
                    best_diff = new_component_score - score_components[i]
                    best_parent_list = deepcopy(parent_list)
                    best_parent_list[i] = new_parents
                end

                # 3) flip an edge
                new_parent_list = deepcopy(parent_list) # TODO: make this more efficient
                deleteat!(new_parent_list[i], idx)

                if adding_edge_preserves_acyclicity(new_parent_list, i, j)
                    sort!(push!(new_parent_list[j], i))
                    new_diff = bayesian_score_component(i, new_parent_list[i], ncategories, datamat, params.prior, params.cache) - score_components[i]
                    new_diff += bayesian_score_component(j, new_parent_list[j], ncategories, datamat, params.prior, params.cache) - score_components[j]
                    if new_diff > best_diff
                        best_diff = new_diff
                        best_parent_list = new_parent_list
                    end
                end
            end
        end

        if best_diff > 0.0
            parent_list = best_parent_list
            score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)
        else
            break
        end
    end

    # construct the BayesNet
    cpds = Array{DiscreteCPD}(n)
    varnames = names(data)
    for i in 1:n
        name = varnames[i]
        parents = varnames[parent_list[i]]
        cpds[i] = fit(DiscreteCPD, data, name, parents, params.prior, parental_ncategories=ncategories[parent_list[i]], target_ncategories=ncategories[i])
    end
    BayesNet(cpds)
end


text:
using CompScienceMeshes
using SauterSchwabQuadrature

pI = point(1,5,3)
pII = point(2,5,3)
pIII = point(7,1,0)
pIV = point(5,1,-3)

Sourcechart = simplex(pII,pI,pIII)
Testchart = simplex(pII,pI,pIV)

Accuracy = 12

function integrand(x,y)
			return(((x-pI)'*(y-pII))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

result = sauterschwabintegral(Sourcechart, Testchart, integrand, Accuracy, Accuracy)
println(result)

#=For those who want to test the sauterschwab_nonparameterized() function,
may uncomment the following five lines=#

#sourcechart = simplex(pI,pIII,pII)
#testchart = simplex(pI,pIV,pII)
#ce = CommonEdge(Accuracy)
#result2 = sauterschwab_nonparameterized(sourcechart, testchart, integrand, ce)
#println(result2)

#=The first argument and the third argument of both simplex() functions are
equal, hence the required conditions are fulfilled. The user may also compare
the two results, and see that both are equal=#


text:
using Parameters

Dim = 3
include("utility/constant.jl")
import QuantumStatistics:Grid
# const QS=QuantumStatistics

_rs = 1.0
_kf = (Dim == 3) ? ((9π / 4.0)^(1.0 / 3) / _rs) : (sqrt(2) / _rs)
_beta, _mass2, _lambda2 = 40.0, 0.0, 1.0
const Mom = MVector{Dim,Float}

@with_kw struct Para
    dim::Int = Dim
    spin::Int = 2
    order::Int = 1
    diagrams = [Diag.Π, ]
    β::Float = _beta / _kf^2
    Rs::Float = _rs
    m²::Float = _mass2
    λ²::Float = _lambda2
    Kf::Float = _kf
    Ef::Float = _kf^2
    μ::Float = _kf^2
    boldG::Bool = false
end

@with_kw struct MCPara
    block::Int = 101
    reweight::Array{Float,1} = [1.0, 0.1, 30.0, 1.0, 0.2, 0.1, 0.01, 0.01]
    printTime::Int = 5
    saveTime::Int = 10
    reweightTime::Int = 10
    messageTime::Int = 10
end

const para = Para()
const mcPara = MCPara()
const grids = (tau = Grid.tau(para.β, 3 * para.Kf^2, 128),
    fermiK = Grid.fermiK(para.Kf, 3.0 * para.Kf, 0.3 * para.Kf, 64),
    boseK = Grid.boseK(para.Kf, 3.0 * para.Kf, 0.3 * para.Kf, 64))

# const groups=[(∅, 0), ]
# for order in 1:para.order
#     for diag in para.diagrams
#         push!(groups, (diag, order))
#     end
# end

@flush printstyled(
    "Rs: $(para.Rs), kF: $(para.Kf), EF: $(para.Ef), β: $(para.β), T/T_F: $(round(1.0 / para.β / para.Ef, digits = 4)), TotalBlock: $(mcPara.block)\n",
    color = :green,
)

max_stars_repo_name: kunyuan/DiagramCalculator

max_issues_repo_name: kunyuan/DiagramCalculator

max_forks_repo_name: kunyuan/DiagramCalculator


text:
module IntegerSmithNormalForm

export SNF!, SNF, SNFWithoutTransform

"""
    (S,B,T) = negateRow(A,i)

  Negate the i-th row of the integer matrix A such that B = SAT holds true
"""
function negateRow(A::AbstractArray{I,2},i) where I
  S = eye(I,size(A,1))
  T = eye(I,size(A,2))
  S[i,i] = -1
  B = S*A*T
  (S,B,T)
end

"""
    (S,B,T) = negateCol(A,i)

  Negate the i-th column of the integer matrix A such that B = SAT holds true
"""
function negateCol(A::AbstractArray{I,2},i) where I
  (S,B,T) = negateRow(A',i)
  (T',B',S')
end

"""
    (S,B,T) = swapRows(A,i,j)

  Negate the i-th and j-th row of the integer matrix A such that B = SAT holds true
"""
function swapRows(A::AbstractArray{I,2},i,j) where I
  S = eye(I,size(A,1))
  T = eye(I,size(A,2))

  S[i,i] = 0
  S[j,j] = 0
  S[i,j] = 1
  S[j,i] = 1

  B = S*A*T
  (S,B,T)
end

"""
    (S,B,T) = swapCols(A,i,j)

  Negate the i-th and j-th column of the integer matrix A such that B = SAT holds true
"""
function swapCols(A::AbstractArray{I,2},i,j) where I
  (S,B,T) = swapRows(A',i,j)
  (T',B',S')
end

"""
    (S,B,T) = addRow(A,i,j,a)

  Add the j-th row a times to the i-th row of the integer matrix A such that B = SAT holds true
"""
function addRow(A::AbstractArray{I,2},i,j,a::I) where I
  S = eye(I,size(A,1))
  T = eye(I,size(A,2))

  S[i,j] = a

  B = S*A*T
  (S,B,T)
end

"""
    (S,B,T) = addCol(A,i,j,a)

  Add the j-th column a times to the i-th column of the integer matrix A such that B = SAT holds true
"""
function addCol(A::AbstractArray{I,2},i,j,a::I) where I
  (S,B,T) = addRow(A',i,j,a)
  (T',B',S')
end

"""
    (S,B,T) = SNF!(A)

Computes the Smith normal form according to
https://www.charite.de/sysbio/people/hoppe/Diplomarbeit_Hoppe.pdf
Algorithmus 1

The matrices fulfill:
* B is a diagonal matrix
* B[i,i] >= 0 for all i
* B[i,i] divides all B[j,j] for all j > i
"""
function SNF!(A::AbstractArray{I,2}) where {I<:Integer}
  d1 = size(A,1)
  d2 = size(A,2)

  S = eye(I,d1)
  T = eye(I,d2)

  for subMatrixInd in 1:(min(d1,d2))
    ASub = A[subMatrixInd:end,subMatrixInd:end]
    d1Sub = size(ASub,1)
    d2Sub = size(ASub,2)

    SSub = eye(I,d1Sub)
    TSub = eye(I,d2Sub)

    nextPivot = findnext(ASub,1)

    if nextPivot != 0

      pivot = ind2sub(size(ASub),nextPivot)

      #Step 2: Put pivot to (1,1)-position
      (S_,ASub,T_) = swapRows(ASub,1,pivot[1])
      SSub = S_ * SSub
      TSub = TSub * T_
      (S_,ASub,T_) = swapCols(ASub,1,pivot[2])
      SSub = S_ * SSub
      TSub = TSub * T_

      divideMatrix = false
      while !divideMatrix
        #check if ASub[1,1] divides all elements in the
        #column and row
        divideCol = false
        divideRow = false
        while !(divideRow & divideCol)
          divideCol=true
          #check if ASub[1,1] divides all elements in its column
          if d1Sub > 1
            for i in 2:d1Sub
              #does not divide
              if mod.(ASub[i,1],ASub[1,1]) != 0
                divideCol=false
                q = div(ASub[i,1],ASub[1,1])
                #Add row 1 q-times to row i
                (S_,ASub,T_) = addRow(ASub,i,1,-q)
                SSub = S_ * SSub
                TSub = TSub * T_
                #Make (1,i) the new pivot on position (1,1)
                (S_,ASub,T_) = swapRows(ASub,1,i)
                SSub = S_ * SSub
                TSub = TSub * T_
                break
              end
            end
            divideCol || break
          end

          divideRow=true
          #check if ASub[1,1] divides all elements in its row
          if d2Sub > 1
            for i in 1:d2Sub
              #does not divide
              if mod.(ASub[1,i],ASub[1,1]) != 0
                divideRow=false
                q = div(ASub[1,i],ASub[1,1])
                #Add col 1 q-times to col i
                (S_,ASub,T_) = addCol(ASub,i,1,-q)
                SSub = S_ * SSub
                TSub = TSub * T_
                #Make (1,i) the new pivot on position (1,1)
                (S_,ASub,T_) = swapCols(ASub,1,i)
                SSub = S_ * SSub
                TSub = TSub * T_
                break
              end
            end
          end
        end

        #eliminate all entries in the pivot row and column
        if d1Sub > 1
          for i = 2:d1Sub
            q = div(ASub[i,1],ASub[1,1])
            (S_,ASub,T_) = addRow(ASub,i,1,-q)
            SSub = S_ * SSub
            TSub = TSub * T_
          end
        end

        if d2Sub > 1
          for i = 2:d2Sub
            q = div(ASub[1,i],ASub[1,1])
            (S_,ASub,T_) = addCol(ASub,i,1,-q)
            SSub = S_ * SSub
            TSub = TSub * T_
          end
        end

        #check if the pivot element divides all elements of the matrix
        divideMatrix = mod.(ASub,ASub[1,1]) == zeros(size(ASub))

        if !divideMatrix
          #Add row of non-dividing index to pivot row
          nonDividing = ind2sub(size(ASub),findnext(mod.(ASub,ASub[1,1]),1))
          (S_,ASub,T_) = addRow(ASub,1,nonDividing[1],1)
          SSub = S_ * SSub
          TSub = TSub * T_
        end
      end

      S_ = eye(I,d1)
      T_ = eye(I,d2)
      S_[subMatrixInd:end,subMatrixInd:end] = SSub
      T_[subMatrixInd:end,subMatrixInd:end] = TSub
      S = S_ * S
      T = T * T_
      A[subMatrixInd:end,subMatrixInd:end] = ASub
    end
  end

  #make diagonal element positive
  for i in 1:min(d1,d2)
    if A[i,i] < 0
      (S_,A,T_) = negateRow(A,i)
      S = S_ * S
      T = T * T_
    end
  end

  (S,A,T)
end

"""
    (S,B,T) = SNF(A)

Computes the Smith normal form according to
https://www.charite.de/sysbio/people/hoppe/Diplomarbeit_Hoppe.pdf
Algorithmus 1

The matrices fulfill:
* B = SAT
* B is a diagonal matrix
* B[i,i] >= 0 for all i
* B[i,i] divides all B[j,j] for all j > i
"""
function SNF(A)
  B = copy(A)
  SNF!(B)
end

"""
    B = SNFWithoutTransform(A)

Computes the Smith normal form according to
https://www.charite.de/sysbio/people/hoppe/Diplomarbeit_Hoppe.pdf
Algorithmus 1

The matrix fulfills:
* There are matrices S,T with |det(S)| = |det(T)| = 1 such that B = SAT
* B is a diagonal matrix
* B[i,i] >= 0 for all i
* B[i,i] divides all B[j,j] for all j > i
"""
function SNFWithoutTransform(A)
  (S,B,T) = SNF(A)
  B
end

end # module