text:
# data loader
function load_data(path)
    file = "legendre_coefficients.json"
    f = JSON.parsefile(path * file)

    λ = get(f, "wavelengths", -1)
    λ = float.(λ)
    spec = get(f, "spectra", -1)

    num_spec = length(spec)
    coeffs_1 = get(spec[1], "legendre_coeffs", -1)
    num_pos = length(coeffs_1)
    num_coeffs = length(coeffs_1[1])

    log10_τ = zeros(num_spec)
    temp = zeros(num_pos, num_spec)
    pos_wafer = zeros(num_pos, num_spec)
    center_real = zeros(num_spec) # center of symmetry
    # dimensions are: coefficients x positions on stripe x spectrogram
    coeffs = zeros(num_coeffs, num_pos, num_spec)
    for i = 1:num_spec
        # log10 dwell time
        log10_τ[i] = get(spec[i], "logtau", -1)

        # temperature profile
        temp[:,i] = float.(get(spec[i], "temp_profile_si", -1))
        # temp[:,i] = float.(get(spec[i], "temp_profile_linear", -1))

        # positions on the wafer
        pos_wafer[:,i] = float.(get(spec[i], "pos_wafer", -1))

        center_real[i] = float(get(spec[i], "center_real", -1))

        # optical coefficients
        coeffs_i = get(spec[i], "legendre_coeffs", -1)
        for j = 1:num_pos
            coeffs[:, j, i] = float.(coeffs_i[j])
        end
    end
    T = temp
    return log10_τ, T, center_real, pos_wafer, coeffs
end

############################ data processing ###################################
# center and normalize the data
function center(x::AbstractVector)
    m, s = mean(x), std(x)
    y = @. (x - m) / s
    return y
end
# returns a linear map from [a, b] to [0, 1]
unit_transform(a, b) = x->(x-a)/(b-a)
inv_unit_transform(a, b) = x -> (b-a)*x + a # inverse unit transform

function normalize_data(data::AbstractArray{<:Real}, domain::NTuple{2, <:Real})
    unit_transform(domain...).(data)
end

######################## dimensionality reduction ##############################
# 2. calculate xrd coefficients based on PCA
using ApproxFun
# calculates Legendre coefficients
function legendre_coefficients(X::AbstractArray, m::Int)
    n = size(X, 1)
    leastsquares! = legendre_coefficients(n, m)
    coefficients = zeros(m, size(X, 2), size(X, 3))
    @threads for i in 1:size(X, 3)
        println("coefficients for stripe $i")
        Xi = @view X[:, :, i]
        Ci = @view coefficients[:, :, i]
        leastsquares!(Ci, Xi)
    end
    coefficients
end
function legendre_coefficients(n::Int, m::Int)
    x = range(-1, 1, length = n)
    A = zeros(n, m)  # Vandermonde matrix
    for k in 1:m
        A[:, k] = Fun(Legendre(), vcat(zeros(k-1), 1)).(x)
    end
    A = qr(A)
    leastsquares(C, X::AbstractMatrix) = ldiv!(C, A, X)
    leastsquares(X::AbstractMatrix) = A \ X
end

############ based on PCA
# n is number of coefficients to calculate
function svd_coefficients(X::AbstractArray, n::Int)
    C = zeros(n, size(X, 2), size(X, 3))
    @threads for i in 1:size(X, 3)
        println("coefficients for stripe $i")
        @views svd_coefficients!(C[:, :, i], X[:, :, i])
    end
    return C
end
function svd_coefficients!(C::AbstractMatrix, X::AbstractMatrix)
    U, S, V = svd(X)
    n = size(C, 1)
    Un = @view U[:, 1:n]
    mul!(C, Un', X) # projection of X onto reduced orthonormal basis
end
# svd coefficients but based on relative tolerance, instead of number of singular values
# returns vector of vectors, i.e. coefficients for for each position
function svd_coefficients(Y::AbstractMatrix, rtol::Real = 1e-2)
    U, S, V = svd(Y)
    tol = maximum(norm, eachcol(Y)) * rtol
    i = findlast(≥(tol), S)
    Vt = @view V'[1:i, :]
    @. Vt = S[1:i]
    y = [v for v in eachrow(Vt)] # each element is vector of singular coefficients by position
end
#############################
# generalized coefficient of determination, where errors are weighted by their standard deviation
# McFadden (in the context of logistic regression), (log?) likelihood ratio index
# x is data, y is model prediction, σ is expected noise for each data point
function likelihood_ratio(x::AbstractVector, y::AbstractVector, σ::AbstractVector = ones(length(x)))
    δ = (x .- mean(x)) ./ σ
    tss = sum(abs2, δ)
    δ = (x .- y) ./ σ
    ess = sum(abs2, δ)
    r2 = 1 - ess / tss
end

# interpolating xrd patterns on common q-grid
function interpolate_xrd_stripe(x::AbstractVector, Y::AbstractMatrix, q::AbstractVector)
    dx = diff(x)
    all(x->isapprox(x, dx[1], atol = 1e-5), dx) || throw("q not uniformly spaced")
    x = range(extrema(x)..., length = length(x))
    Yq = similar(Y)
    for i in 1:size(Y, 2)
        y = @views Y[:, i]
        f = CubicSplineInterpolation(x, y)
        @. Yq[:, i] = f(q)
    end
    return Yq
end


text:
using ITensorsGrad
using Documenter

DocMeta.setdocmeta!(ITensorsGrad, :DocTestSetup, :(using ITensorsGrad); recursive=true)

makedocs(;
  modules=[ITensorsGrad],
  authors="Matthew Fishman <mfishman@flatironinstitute.org>",
  repo="https://github.com/mtfishman/ITensorsGrad.jl/blob/{commit}{path}#{line}",
  sitename="ITensorsGrad.jl",
  format=Documenter.HTML(;
    prettyurls=get(ENV, "CI", "false") == "true",
    canonical="https://mtfishman.github.io/ITensorsGrad.jl",
    assets=String[],
  ),
  pages=["Home" => "index.md"],
)

deploydocs(;
  repo="github.com/mtfishman/ITensorsGrad.jl.git", devbranch="main", push_preview=true
)

ext: jl

max_stars_repo_name: mtfishman/ITensorsGrad

max_issues_repo_name: mtfishman/ITensorsGrad

max_forks_repo_name: mtfishman/ITensorsGrad


text:
# """
#     toDataFrame(x::AbstractVector{<:SLArray})

# Convert an Vector{SLVector} to DataFrame.

# Such objects are returned when broadcasting on a function
# that retuns an SLVector.
# """
# function toDataFrame(x::AbstractVector{<:SLArray})
#     names = collect(keys(first(x)))
#     df = DataFrame(Tables.table(VectorOfArray(x)'))
#     rename!(df, names)
# end

# function vec_tuple_to_DataFrame(data; names = string.(keys(first(data))))
#     DataFrame(collect(map(idx -> getindex.(data, idx), eachindex(collect(first(data))))), collect(names))
# end

"""
    frac_hour(float::Number)
    frac_hour(period::Type{<:Period}, float::Number)

Create a period in given type (defaults to `Nanosecond`) from
fractional hours.

```jldoctest; output = false
using Dates
frac_hour(1+1/60) == Hour(1) + Minute(1)
# output
true
```
"""
function frac_hour(period::Type{<:Period}, float::Number)
    #adapted from https://stackoverflow.com/a/51448499
    full_hour, Δ = divrem(float, 1)
    partial = period(round(Dates.value(period(Hour(1))) * Δ))
    Hour(full_hour) + partial
end
frac_hour(float::Number) = frac_hour(Nanosecond, float)

"""
    moving_average(vs,n; nmin = n ÷ 2)

Compute the moving average over a vector allowing for missings.

# Arguments
- vs: numeric vector to average over
- n: window size: number of items to average
- nmin: minimum number of non-missing records

Values outside the edges are assumed missing.

If the number of non-missing records within a window is smaller than nmin
then the averaged value is assumed missing. This avoids average at edges of
periods with many missings to be very sensitive to the edge values.
"""
function moving_average(vs,n; nmin = n ÷ 2)
    kernel = -trunc(Int,(n-1)/2):ceil(Int,(n-1)/2)
    #[mean(skipmissing(@view vs[i:(i+n-1)])) for i in 1:(length(vs)-(n-1))]
    vsp = PaddedView(missing, vs, (-n:length(vs)+n,))
    fagg = function(x)
        sum(.!ismissing.(x)) < nmin && return(missing)
        mean(skipmissing(x))
    end
    #ans = [fagg(@view vsp[i.+kernel]) for i in 100:101]
    ans = [fagg(@view vsp[i.+kernel]) for i in 1:length(vs)]
end

"""
    get_nonoverlapping_periods(dfo::DataFrame)

Get non-overlapping periods.

# Arguments
- `dfo`: DataFrame with columns `p_start` and `p_end` sorted by increasing `p_start`

# Value
- DataFrame with columns `p_start` and `p_end` with `p_start[i] > p_end[i-1]`
"""
function get_nonoverlapping_periods(dfo::DataFrame)
    nrow(dfo) == 0 && return(dfo)
    dfno = DataFrame()
    p_start, p_end = dfo.p_start[1], dfo.p_end[1]
    for i in 2:(nrow(dfo)-1)
        if dfo.p_start[i] <= p_end
            # if current row is in previous intervals, just extend the previous interval
            p_end = max(p_end, dfo.p_end[i])
        else
            # push the previous interval and take this row as interval
            push!(dfno, (p_start = p_start, p_end = p_end))
            p_start, p_end = dfo.p_start[i], dfo.p_end[i]
        end
    end
    if dfo.p_start[end] <= p_end
        # last row start within previous interval: push extended interval
        push!(dfno, (p_start = p_start, p_end = max(p_end, dfo.p_end[end])))
    else
        # last row not in previous interval: push both interval and last row
        push!(dfno, (p_start = p_start, p_end = p_end))
        push!(dfno, (p_start = dfo.p_start[end], p_end = dfo.p_end[end]))
    end
    dfno
end

"""
    set_datetime_ydh!(df, timezone = tz"UTC+0";
        year = df.year, doy = df.doy, hour = df.hour)

Update column DateTime column datetime according to year, dayOfYear, and fractional hour.

# Value
`df` with updated or added column `:datetime` moved to the first position.
"""
function set_datetime_ydh!(df, timezone = tz"UTC+0";
    year = df.year, doy = df.doy, hour = df.hour)
    df[!,:datetime] = @. ZonedDateTime(
        DateTime(year) + Day(doy-1) + frac_hour(Second, hour), timezone)
    select!(df, :datetime, :)
end


text:
using RecipesBase

function planar_contour(p::Polyhedron)
    if fulldim(p) != 2
        if fulldim(p) == 3
            error("Plotting 3-dimensional polyhedron with Plots is not supported, use Makie or MeshCat.")
        else
            error("Plotting $(fulldim(p))-dimensional polyhedron with Plots is not supported.")
        end
    end
    removevredundancy!(p)
    if hasallrays(p)
        error("Rays not supported yet in the 2D plotting recipe.")
    end
    ps = collect(points(p))
    if isempty(ps)
        error("Plotting empty polyhedron is not supported.")
    end
    sort!(ps, by = x -> x[1])
    counterclockwise(p1, p2) = dot(cross([p1; 0], [p2; 0]), [0, 0, 1])
    top = getsemihull(ps,  1, counterclockwise)
    bot = getsemihull(ps, -1, counterclockwise)
    if !isempty(top) && !isempty(bot)
        @assert top[end] == bot[1]
        pop!(top)
    end
    if !isempty(bot) && !isempty(top)
        @assert bot[end] == top[1]
        pop!(bot)
    end
    hull = [top; bot]
    push!(hull, hull[1])  # ensure shape is closed
    return [p[1] for p in hull], [p[2] for p in hull]
end

@recipe function f(p::Polyhedron)
    seriestype --> :shape
    legend --> false
    planar_contour(p)
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
using Plots
using Statistics

"""
    PlotSpectroTemporal(data, sample_rate)

This function plots the time, spectrogram, and periodogram of a signal.

"""
function PlotSpectroTemporal(x::AbstractArray, sample_rate::Number;
                                figure_size::Tuple=(750, 400),
                                window = hamming,
                                amplitude_limits = nothing,
                                power_limits = nothing,
                                time_limits = nothing,
                                frequency_limits = [0, 1500],
                                correlation_annotate = true, kwargs...)

    # Generate time vector
    t = 1:size(x, 1)
    t = t ./ sample_rate
    if time_limits == nothing
        time_limits = [0, maximum(t)]
    end

    # Calculate signal transforms
    spec = spectrogram(x[:, 1], 1024, 256, fs = sample_rate, window = window)
    peri1 = welch_pgram(x[:, 1], 2048, 512, fs = sample_rate)
    peri1_power = vector_pow2db(power(peri1))
    if size(x, 2)>1; peri2 = welch_pgram(x[:, 2], 2048, 512, fs = sample_rate); end
    tfft = fft(x[:, 1])

    # Extract required stats
    if isnothing(amplitude_limits)
        amplitude_limits = maximum(abs.(x))
        amplitude_limits = (-amplitude_limits, amplitude_limits)
    end
    if isnothing(power_limits)
        idxs = (freq(peri1) .> frequency_limits[1]) .& (freq(peri1) .< frequency_limits[2])
        power_limits = (minimum(peri1_power[idxs]), maximum(peri1_power[idxs]))
    end

    # Create plots
    spec_plot = heatmap(time(spec), freq(spec), power(spec), colorbar = false, xlab = "Time (s)", ylab = "Frequency (Hz)", ylims = frequency_limits, xlims = time_limits)

    peri_plot = plot(peri1_power, freq(peri1), yticks = [], xlab = "Power (dB)", lab = "", ylims = frequency_limits, xlims = power_limits)
    if size(x, 2)>1; peri_plot = plot!(vector_pow2db(power(peri2)), freq(peri2), yticks = [], xlab = "Power (dB)", lab = "", ylims = frequency_limits); end

    time_plot = plot(t, x, xticks = [], leg = false, ylab = "Amplitude", lab = "", ylims = amplitude_limits, xlims = time_limits)

    corr_plot = histogram(x, orientation = :h,  ticks = [], leg = false, framestyle = :none, link = :none, bins=LinRange(amplitude_limits[1], amplitude_limits[2], 15), ylims = amplitude_limits)
    if ((size(x, 2)>1) & correlation_annotate)

        maximum_hist_val = maximum(corr_plot.series_list[2].plotattributes[:y][.~isnan.(corr_plot.series_list[2].plotattributes[:y])])

        d = annotate!(0, 0.9 * maximum(amplitude_limits), text(string("Corr = ", round(cor(x)[2, 1], digits=3)),:left,8))
        d = annotate!(0, -0.9 * maximum(amplitude_limits), text(string("Std = ", round(std(x), digits=3)),:left,8))
    end

    # Return all plots in layout
    l = Plots.@layout [c{0.6w, 0.3h} d  ; a{0.8w} b]
    return plot(time_plot, corr_plot, spec_plot, peri_plot, layout = l, size = figure_size, link = :none; kwargs...)

end

function vector_pow2db(a::Vector)
    for index = 1:size(a, 1)
        a[index] = pow2db(a[index])
    end
    return a
end

function plot_cross_correlation(x::T; lags::AbstractQuantity=0u"s") where {T<:SampledSignals.SampleBuf}

    lags_seconds = lags |> u"s" |> ustrip
    lags_samples = Int(lags_seconds * x.samplerate)
    plot_cross_correlation(x.data, lags_samples, x.samplerate)
end

function plot_cross_correlation(x::T; lags::AbstractQuantity=0u"s") where {T<:DummySampleSink}

    lags_seconds = lags |> u"s" |> ustrip
    lags_samples = Int(lags_seconds * x.samplerate)
    plot_cross_correlation(x.buf, lags_samples, x.samplerate)
end

"""
    plot_cross_correlation(x::SampleBuf, lags::Unitful.Time)

Plot the cross correlation of a two channel audio signal.

Inputs
------
* `x` data in the form of SampledSignals.SampleBuf. Must be two channels of audio.
* `lags` time range of lags to be used when plotting the cross correlation function.
  If lags=0, then the entire function will be used, effecively same as lags=Inf.

Example
-------
```julia
correlation = 0.6
source = CorrelatedNoiseSource(Float64, 48u"kHz", 2, 0.1, correlation)
a = read(source, 3u"s")
plot_cross_correlation(a, lags=4u"ms")
```
"""
function plot_cross_correlation(x::Array{T, 2}, lags::Int, samplerate::Number) where {T<:Number}

    if lags == 0
        lags = size(x, 1) - 1
    end
    lags = round.(Int, -lags:1:lags)

    lag_times = lags ./ samplerate
    lag_times = lag_times .* 1.0u"s"

    plot(lag_times, crosscor(x[:, 1], x[:, 2], lags),
         label="", ylab="Cross Correlation", xlab="Lag", ylims=(-1, 1))

end


text:
using BenchmarkTools
using Trixi

const SUITE = BenchmarkGroup()

for elixir in [joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_extended.jl"),
               joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_amr_nonperiodic.jl"),
               joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_ec.jl"),
               joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_vortex_mortar.jl"),
               joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_vortex_mortar_shockcapturing.jl"),
               joinpath(examples_dir(), "tree_2d_dgsem", "elixir_mhd_ec.jl"),
               joinpath(examples_dir(), "structured_2d_dgsem", "elixir_advection_extended.jl"),
               joinpath(examples_dir(), "structured_2d_dgsem", "elixir_advection_nonperiodic.jl"),
               joinpath(examples_dir(), "structured_2d_dgsem", "elixir_euler_ec.jl"),
               joinpath(examples_dir(), "structured_2d_dgsem", "elixir_euler_source_terms_nonperiodic.jl"),
               joinpath(examples_dir(), "structured_2d_dgsem", "elixir_mhd_ec.jl"),
               joinpath(examples_dir(), "unstructured_2d_dgsem", "elixir_euler_wall_bc.jl"), # this is the only elixir working for polydeg=3
               joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_extended.jl"),
               joinpath(@__DIR__, "elixir_2d_euler_vortex_tree.jl"),
               joinpath(@__DIR__, "elixir_2d_euler_vortex_structured.jl"),
               joinpath(@__DIR__, "elixir_2d_euler_vortex_unstructured.jl"),
               joinpath(@__DIR__, "elixir_2d_euler_vortex_p4est.jl"),
               joinpath(examples_dir(), "tree_3d_dgsem", "elixir_advection_extended.jl"),
               joinpath(examples_dir(), "tree_3d_dgsem", "elixir_euler_ec.jl"),
               joinpath(examples_dir(), "tree_3d_dgsem", "elixir_euler_mortar.jl"),
               joinpath(examples_dir(), "tree_3d_dgsem", "elixir_euler_shockcapturing.jl"),
               joinpath(examples_dir(), "tree_3d_dgsem", "elixir_mhd_ec.jl"),
               joinpath(examples_dir(), "structured_3d_dgsem", "elixir_advection_nonperiodic.jl"),
               joinpath(examples_dir(), "structured_3d_dgsem", "elixir_euler_ec.jl"),
               joinpath(examples_dir(), "structured_3d_dgsem", "elixir_euler_source_terms_nonperiodic.jl"),
               joinpath(examples_dir(), "structured_3d_dgsem", "elixir_mhd_ec.jl"),
               joinpath(examples_dir(), "p4est_3d_dgsem", "elixir_advection_basic.jl"),]
  benchname = basename(dirname(elixir)) * "/" * basename(elixir)
  SUITE[benchname] = BenchmarkGroup()
  for polydeg in [3, 7]
    trixi_include(elixir, tspan=(0.0, 1.0e-10); polydeg)
    SUITE[benchname]["p$(polydeg)_rhs!"] = @benchmarkable Trixi.rhs!(
      $(similar(sol.u[end])), $(copy(sol.u[end])), $(semi), $(first(tspan)))
    SUITE[benchname]["p$(polydeg)_analysis"] = @benchmarkable ($analysis_callback)($sol)
  end
end

let
  SUITE["latency"] = BenchmarkGroup()
  SUITE["latency"]["default_example"] = @benchmarkable run(
    `$(Base.julia_cmd()) -e 'using Trixi; trixi_include(default_example())'`) seconds=60
  for polydeg in [3, 7]
    command = "using Trixi; trixi_include(joinpath(examples_dir(), \"tree_2d_dgsem\", \"elixir_advection_extended.jl\"), tspan=(0.0, 1.0e-10), polydeg=$(polydeg), save_restart=TrivialCallback(), save_solution=TrivialCallback())"
    SUITE["latency"]["polydeg_$polydeg"] = @benchmarkable run($`$(Base.julia_cmd()) -e $command`) seconds=60
  end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the Euler angles.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                                  Operations
################################################################################

"""
    *(Θ₂::EulerAngles, Θ₁::EulerAngles)

Compute the composed rotation of `Θ₁ -> Θ₂`.

The rotation will be represented by Euler angles (see [`EulerAngles`](@ref))
with the same rotation sequence as `Θ₂`.

# Examples

```jldoctest
julia> ea1 = EulerAngles(deg2rad(35), 0, 0, :XYZ)
EulerAngles{Float64}:
  R(X) :  0.610865 rad  ( 35.0°)
  R(Y) :  0.0      rad  ( 0.0°)
  R(Z) :  0.0      rad  ( 0.0°)

julia> ea2 = EulerAngles(0, 0, deg2rad(25), :ZYX)
EulerAngles{Float64}:
  R(Z) :  0.0      rad  ( 0.0°)
  R(Y) :  0.0      rad  ( 0.0°)
  R(X) :  0.436332 rad  ( 25.0°)

julia> ea2 * ea1
EulerAngles{Float64}:
  R(Z) :  0.0    rad  ( 0.0°)
  R(Y) : -0.0    rad  (-0.0°)
  R(X) :  1.0472 rad  ( 60.0°)
```
"""
@inline function *(Θ₂::EulerAngles, Θ₁::EulerAngles)
    # Convert to quaternions, compute the composition, and convert back to Euler
    # angles.
    q₁ = angle_to_quat(Θ₁)
    q₂ = angle_to_quat(Θ₂)

    return quat_to_angle(q₁ * q₂, Θ₂.rot_seq)
end

"""
    inv(Θ::EulerAngles)

Return the Euler angles that represent the inverse rotation of `Θ`.

The rotation sequence of the result will be the inverse of the input. Hence, if
the input rotation sequence is, for example, `:XYZ`, then the result will be
represented using `:ZYX`.

# Examples

```julia-repl
julia> ea = EulerAngles(π / 3, π / 6,  2 / 3 * π, :ZYX)
EulerAngles{Float64}:
  R(Z) :  1.0472   rad  ( 60.0°)
  R(Y) :  0.523599 rad  ( 30.0°)
  R(X) :  2.0944   rad  ( 120.0°)

julia> inv(ea)
EulerAngles{Float64}:
  R(X) : -2.0944   rad  (-120.0°)
  R(Y) : -0.523599 rad  (-30.0°)
  R(Z) : -1.0472   rad  (-60.0°)
```
"""
function inv(Θ::EulerAngles)
    # Check what will be the inverse rotation.
    if Θ.rot_seq == :XYZ
        inv_rot_seq = :ZYX
    elseif Θ.rot_seq == :XZY
        inv_rot_seq = :YZX
    elseif Θ.rot_seq == :YXZ
        inv_rot_seq = :ZXY
    elseif Θ.rot_seq == :YZX
        inv_rot_seq = :XZY
    elseif Θ.rot_seq == :ZXY
        inv_rot_seq = :YXZ
    elseif Θ.rot_seq == :ZYX
        inv_rot_seq = :XYZ
    else
        inv_rot_seq = Θ.rot_seq
    end

    # Return the Euler angle that represented the inverse rotation.
    return EulerAngles(-Θ.a3, -Θ.a2, -Θ.a1, inv_rot_seq)
end

################################################################################
#                                      IO
################################################################################

function show(io::IO, Θ::EulerAngles{T}) where T
    # Get if `io` request a compact printing, defaulting to true.
    compact_printing = get(io, :compact, true)

    # Convert the values using `print` and compact printing.
    θ₁_str = sprint(print, Θ.a1; context = :compact => compact_printing)
    θ₂_str = sprint(print, Θ.a2; context = :compact => compact_printing)
    θ₃_str = sprint(print, Θ.a3; context = :compact => compact_printing)
    rot_seq = String(Θ.rot_seq)

    print(io, "EulerAngles{$T}:")
    print(io, " R($rot_seq)  " * θ₁_str * "  " * θ₂_str * "  " * θ₃_str * " rad")

    return nothing
end

function show(io::IO, mime::MIME"text/plain", Θ::EulerAngles{T}) where T
    # Check if the user wants colors.
    color = get(io, :color, false)::Bool

    # Check if the user wants compact printing, defaulting to `true`.
    compact_printing = get(io, :compact, true)::Bool

    # Assemble the context.
    context = IOContext(io, :color => color, :compact => compact_printing)

    d = color ? string(_d) : ""
    g = color ? string(_g) : ""
    y = color ? string(_y) : ""
    u = color ? string(_u) : ""

    # Convert the values using `print` and compact printing.
    θ₁_alg = !signbit(Θ.a1) ? " " : ""
    θ₁_str = θ₁_alg * sprint(print, Θ.a1; context = context)
    lθ₁ = length(θ₁_str)

    θ₂_alg = !signbit(Θ.a2) ? " " : ""
    θ₂_str = θ₂_alg * sprint(print, Θ.a2; context = context)
    lθ₂ = length(θ₂_str)

    θ₃_alg = !signbit(Θ.a3) ? " " : ""
    θ₃_str = θ₃_alg * sprint(print, Θ.a3; context = context)
    lθ₃ = length(θ₃_str)

    alignment_pad = max(lθ₁, lθ₂, lθ₃)

    θ₁_str *= " "^(alignment_pad - lθ₁) * " rad  (" *
        θ₁_alg * sprint(print, rad2deg(Θ.a1); context = context) * "°)"

    θ₂_str *= " "^(alignment_pad - lθ₂) * " rad  (" *
        θ₂_alg * sprint(print, rad2deg(Θ.a2); context = context) * "°)"

    θ₃_str *= " "^(alignment_pad - lθ₃) * " rad  (" *
        θ₃_alg * sprint(print, rad2deg(Θ.a3); context = context) * "°)"

    rot_seq = String(Θ.rot_seq)

    println(io, "EulerAngles{$T}:")
    println(io, "$g  R($(rot_seq[1])) : $d" * θ₁_str)
    println(io, "$y  R($(rot_seq[2])) : $d" * θ₂_str)
    print(io,   "$u  R($(rot_seq[3])) : $d" * θ₃_str)

    return nothing
end


text:
function LColorbar(parent::Scene, plot::AbstractPlot; kwargs...)

    LColorbar(parent;
        colormap = plot.colormap,
        limits = plot.colorrange,
        kwargs...
    )

end

function LColorbar(parent::Scene; bbox = nothing, kwargs...)
    attrs = merge!(Attributes(kwargs), default_attributes(LColorbar, parent).attributes)

    default_attrs = default_attributes(LColorbar, parent).attributes
    theme_attrs = subtheme(parent, :LColorbar)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)

    @extract attrs (
        label, labelcolor, labelsize, labelvisible, labelpadding, ticklabelsize,
        ticklabelspace, labelfont, ticklabelfont, ticklabelcolor,
        ticklabelsvisible, ticks, tickformat, ticksize, ticksvisible, ticklabelpad, tickalign,
        tickwidth, tickcolor, spinewidth, topspinevisible,
        rightspinevisible, leftspinevisible, bottomspinevisible, topspinecolor,
        leftspinecolor, rightspinecolor, bottomspinecolor, colormap, limits,
        halign, valign, vertical, flipaxisposition, ticklabelalign, flip_vertical_label,
        nsteps)

    decorations = Dict{Symbol, Any}()

    protrusions = Node(GridLayoutBase.RectSides{Float32}(0, 0, 0, 0))
    layoutobservables = LayoutObservables(LColorbar, attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        halign, valign, attrs.alignmode; suggestedbbox = bbox, protrusions = protrusions)

    framebox = layoutobservables.computedbbox

    colorlinepoints = lift(framebox, nsteps) do fb, nsteps
        fbw = fb.widths[1]
        fbh = fb.widths[2]

        if vertical[]
            [Point2f0(0.5f0 * fbw, y * fbh) for y in LinRange(0f0, 1f0, nsteps)]
        else
            [Point2f0(x * fbw, 0.5f0 * fbh) for x in LinRange(0f0, 1f0, nsteps)]
        end
    end

    linewidth = lift(framebox) do fb
        fbw = fb.widths[1]
        fbh = fb.widths[2]

        if vertical[]
            fbw
        else
            fbh
        end
    end

    xrange = lift(framebox) do fb
        range(left(fb), right(fb), length = 2)
    end
    yrange = lift(framebox) do fb
        range(bottom(fb), top(fb), length = 2)
    end

    colorcells = lift(vertical, nsteps) do v, nsteps
        if v
            reshape(collect(1:nsteps), 1, :)
        else
            reshape(collect(1:nsteps), :, 1)
        end
    end

    hm = heatmap!(parent, xrange, yrange, colorcells, colormap = colormap, raw = true)[end]
    decorations[:heatmap] = hm

    ab, al, ar, at = axislines!(
        parent, framebox, spinewidth, topspinevisible, rightspinevisible,
        leftspinevisible, bottomspinevisible, topspinecolor, leftspinecolor,
        rightspinecolor, bottomspinecolor)
    decorations[:topspine] = at
    decorations[:leftspine] = al
    decorations[:rightspine] = ar
    decorations[:bottomspine] = ab

    axispoints = lift(framebox, vertical, flipaxisposition) do scenearea,
            vertical, flipaxisposition

        if vertical
            if flipaxisposition
                (bottomright(scenearea), topright(scenearea))
            else
                (bottomleft(scenearea), topleft(scenearea))
            end
        else
            if flipaxisposition
                (topleft(scenearea), topright(scenearea))
            else
                (bottomleft(scenearea), bottomright(scenearea))
            end
        end

    end

    axis = LineAxis(parent, endpoints = axispoints, flipped = flipaxisposition,
        limits = limits, ticklabelalign = ticklabelalign, label = label,
        labelpadding = labelpadding, labelvisible = labelvisible, labelsize = labelsize,
        labelfont = labelfont, ticklabelfont = ticklabelfont, ticks = ticks, tickformat = tickformat,
        ticklabelsize = ticklabelsize, ticklabelsvisible = ticklabelsvisible, ticksize = ticksize,
        ticksvisible = ticksvisible, ticklabelpad = ticklabelpad, tickalign = tickalign,
        tickwidth = tickwidth, tickcolor = tickcolor, spinewidth = spinewidth,
        ticklabelspace = ticklabelspace, ticklabelcolor = ticklabelcolor,
        spinecolor = :transparent, spinevisible = :false, flip_vertical_label = flip_vertical_label)
    decorations[:axis] = axis

    onany(axis.protrusion, vertical, flipaxisposition) do axprotrusion,
            vertical, flipaxisposition

        left, right, top, bottom = 0f0, 0f0, 0f0, 0f0

        if vertical
            if flipaxisposition
                right += axprotrusion
            else
                left += axprotrusion
            end
        else
            if flipaxisposition
                top += axprotrusion
            else
                bottom += axprotrusion
            end
        end

        protrusions[] = GridLayoutBase.RectSides{Float32}(left, right, bottom, top)
    end

    # trigger protrusions with one of the attributes
    vertical[] = vertical[]

    # trigger bbox
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    LColorbar(parent, layoutobservables, attrs, decorations)
end

function tight_ticklabel_spacing!(lc::LColorbar)
    tight_ticklabel_spacing!(lc.decorations[:axis])
end


text:
isorthogonal(::AbstractQ) = true
isorthogonal(q) = q'q ≈ I

"""
    QLHessenberg(factors, q)

represents a Hessenberg QL factorization where factors contains L in its
lower triangular components and q is a vector of 2x2 orthogonal transformations
whose product gives Q.
"""
struct QLHessenberg{T,S<:AbstractMatrix{T},QT<:AbstractVector{<:AbstractMatrix{T}}} <: Factorization{T}
    factors::S
    q::QT

    function QLHessenberg{T,S,QT}(factors, q) where {T,S<:AbstractMatrix{T},QT<:AbstractVector{<:AbstractMatrix{T}}}
        require_one_based_indexing(factors)
        new{T,S,QT}(factors, q)
    end
end

QLHessenberg(factors::AbstractMatrix{T}, q::AbstractVector{<:AbstractMatrix{T}}) where {T} = QLHessenberg{T,typeof(factors),typeof(q)}(factors, q)
QLHessenberg{T}(factors::AbstractMatrix, q::AbstractVector) where {T} =
    QLHessenberg(convert(AbstractMatrix{T}, factors), convert.(AbstractMatrix{T}, q))

# iteration for destructuring into components
Base.iterate(S::QLHessenberg) = (S.Q, Val(:L))
Base.iterate(S::QLHessenberg, ::Val{:L}) = (S.L, Val(:done))
Base.iterate(S::QLHessenberg, ::Val{:done}) = nothing

QLHessenberg{T}(A::QLHessenberg) where {T} = QLHessenberg(convert(AbstractMatrix{T}, A.factors), convert.(AbstractMatrix{T}, A.q))
Factorization{T}(A::QLHessenberg{T}) where {T} = A
Factorization{T}(A::QLHessenberg) where {T} = QLHessenberg{T}(A)
AbstractMatrix(F::QLHessenberg) = F.Q * F.L
AbstractArray(F::QLHessenberg) = AbstractMatrix(F)
Matrix(F::QLHessenberg) = Array(AbstractArray(F))
Array(F::QLHessenberg) = Matrix(F)

function show(io::IO, mime::MIME{Symbol("text/plain")}, F::QLHessenberg)
    summary(io, F); println(io)
    println(io, "Q factor:")
    show(io, mime, F.Q)
    println(io, "\nL factor:")
    show(io, mime, F.L)
end

@inline function getL(F::QLHessenberg, _)
    m, n = size(F)
    tril!(getfield(F, :factors)[end-min(m,n)+1:end, 1:n], max(n-m,0))
end
@inline getQ(F::QLHessenberg, _) = LowerHessenbergQ(F.q)

getL(F::QLHessenberg) = getL(F, size(F.factors))
getQ(F::QLHessenberg) = getQ(F, size(F.factors))

function getproperty(F::QLHessenberg, d::Symbol)
    if d == :L
        return getL(F)
    elseif d == :Q
        return getQ(F)
    else
        getfield(F, d)
    end
end

Base.propertynames(F::QLHessenberg, private::Bool=false) =
    (:L, :Q, (private ? fieldnames(typeof(F)) : ())...)

abstract type AbstractHessenbergQ{T} <: LayoutQ{T} end

for Typ in (:UpperHessenbergQ, :LowerHessenbergQ)
    @eval begin
        struct $Typ{T, QT<:AbstractVector{<:AbstractMatrix{T}}} <: AbstractHessenbergQ{T}
            q::QT
            function $Typ{T,QT}(q::QT) where {T,QT<:AbstractVector{<:AbstractMatrix{T}}}
                all(isorthogonal.(q)) || throw(ArgumentError("input must be orthogonal"))
                all(size.(q) .== Ref((2,2))) ||  throw(ArgumentError("input must be 2x2"))
                new{T,QT}(q)
            end
        end

        $Typ(q::AbstractVector{<:AbstractMatrix{T}}) where T =
            $Typ{T,typeof(q)}(q)
    end
end

size(Q::AbstractHessenbergQ, k::Integer) = length(Q.q)+1
axes(Q::AbstractHessenbergQ, k::Integer) = oneto(length(Q.q)+1)
size(F::QLHessenberg, dim::Integer) = size(getfield(F, :factors), dim)
size(F::QLHessenberg) = size(getfield(F, :factors))

bandwidths(Q::UpperHessenbergQ) = (1,size(Q,2)-1)
bandwidths(Q::LowerHessenbergQ) = (size(Q,1)-1,1)

adjoint(Q::UpperHessenbergQ) = LowerHessenbergQ(adjoint.(Q.q))
adjoint(Q::LowerHessenbergQ) = UpperHessenbergQ(adjoint.(Q.q))

check_mul_axes(A::AbstractHessenbergQ, B, C...) =
    axes(A,2) == axes(B,1) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and first axis of B, $(axes(B,1)) must match"))

struct HessenbergQLayout{UPLO} <: AbstractQLayout end

MemoryLayout(::Type{<:UpperHessenbergQ}) = HessenbergQLayout{'U'}()
MemoryLayout(::Type{<:LowerHessenbergQ}) = HessenbergQLayout{'L'}()

function materialize!(L::MatLmulVec{<:HessenbergQLayout{'L'}})
    Q, x = L.A,L.B
    T = eltype(Q)
    t = Array{T}(undef, 2)
    nz = nzzeros(x,1)
    for n = 1:length(Q.q)
        v = view(x, n:n+1)
        mul!(t, Q.q[n], v)
        v .= t
        n > nz && norm(t) ≤ 10floatmin(real(T)) && return x
    end
    x
end

function materialize!(L::MatLmulVec{<:HessenbergQLayout{'U'}})
    Q, x = L.A,L.B
    T = eltype(Q)
    t = Array{T}(undef, 2)
    for n = min(length(Q.q),nzzeros(x,1)):-1:1
        v = view(x, n:n+1)
        mul!(t, Q.q[n], v)
        v .= t
    end
    x
end

function materialize!(L::MatLmulMat{<:HessenbergQLayout})
    Q,X = L.A,L.B
    for j in axes(X,2)
        ArrayLayouts.lmul!(Q, view(X,:,j))
    end
    X
end

###
# getindex
####

function getindex(Q::UpperHessenbergQ, i::Int, j::Int)
    y = zeros(eltype(Q), size(Q, 2))
    y[j] = 1
    lmul!(Q, y)[i]
end

getindex(Q::LowerHessenbergQ, i::Int, j::Int) = (Q')[j,i]'

###
# QLPackedQ, QRPackedQ <-> Lower/UpperHessenbergQ
###

UpperHessenbergQ(Q::LinearAlgebra.QRPackedQ) = UpperHessenbergQ(QRPackedQ(Q))

function UpperHessenbergQ(Q::QRPackedQ{T}) where T
    @assert bandwidth(Q.factors,1) == 1
    q = Vector{Matrix{T}}()
    for j = 1:length(Q.τ)-2
        push!(q, QRPackedQ(Q.factors[j:j+1,j:j+1], [Q.τ[j],zero(T)]))
    end
    push!(q, QRPackedQ(Q.factors[end-1:end,end-1:end], Q.τ[end-1:end]))
    UpperHessenbergQ(q)
end

function LowerHessenbergQ(Q::QLPackedQ{T}) where T
    @assert bandwidth(Q.factors,2) == 1
    q = Vector{Matrix{T}}()
    push!(q, QLPackedQ(Q.factors[1:2,1:2], Q.τ[1:2]))
    for j = 2:length(Q.τ)-1
        push!(q, QLPackedQ(Q.factors[j:j+1,j:j+1], [zero(T),Q.τ[j+1]]))
    end
    LowerHessenbergQ(q)
end


text:
using DataStructures

head = HalfPath()
tail = HalfPath()
for x in 3:-1:1; enqueue!(head, x); end
for x in 4:6; enqueue!(tail, x); end
path = Path(head, tail)
@test collect(path) == [x for x in 1:6]


text:
using Plots
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

x̄, ū = unpack(z̄_nom, prob_nominal)
x, u = unpack(z, prob_nominal)

t_nominal = range(0, stop = sum([ū[t][end] for t = 1:T-1]), length = T)
plot(t_nominal[1:end-1], hcat(ū...)[1:end-1, :]', linetype = :steppost)
@show sum([ū[t][end] for t = 1:T-1]) # works when 2.72

# COM traj
xx_nom = [x̄[t][1] for t = 1:T]
zz_nom = [x̄[t][2] for t = 1:T]
x̄_dpo, ū_dpo = unpack(z[prob_dpo.prob.idx.nom], prob_dpo.prob.prob.nom)
t_dpo = range(0, stop = sum([ū_dpo[t][end] for t = 1:T-1]), length = T)

xx = [x̄_dpo[t][1] for t = 1:T]
zz = [x̄_dpo[t][2] for t = 1:T]

plot(xx_nom, zz_nom, color = :cyan)
plot!(xx, zz, color = :orange)

plot(t_nominal, hcat(x̄...)', color = :cyan)
plot!(t_dpo, hcat(x̄_dpo...)', color = :orange)

# orientation tracking
plt_x = plot(t_nominal, hcat(x̄...)[3,:],legend=:topright,color=:red,
    label="",width=2.0,xlabel="time (s)",
    title="Rocket",ylabel="state")
plt_x = plot!(t_sim_nom,hcat(z_lqr...)[3,:],color=:black,label="",
    width=1.0)

plt_x = plot(t_dpo,hcat(x̄_dpo...)[3,:],legend=:topright,color=:red,
    label="",width=2.0,xlabel="time (s)",
    title="Rocket",ylabel="state")
plt_x = plot!(t_sim_dpo,hcat(z_dpo...)[3,:],color=:black,label="",
    width=1.0)

visualize!(vis, model_sl, z_lqr, Δt = dt_sim_nom)
visualize!(vis, model_sl, z_dpo, Δt = dt_sim_dpo)

using PGFPlots
const PGF = PGFPlots

p_traj_nom = PGF.Plots.Linear(hcat(x̄...)[1,:],hcat(x̄...)[2,:],
	mark="none",style="color=cyan, line width = 2pt",legendentry="TO")
p_traj_sample = PGF.Plots.Linear(hcat(x...)[1,:],hcat(x...)[2,:],
	mark="none",style="color=orange, line width = 2pt",legendentry="DPO")

a = Axis([p_traj_nom;p_traj_sample],
    axisEqualImage=false,
    hideAxis=false,
	ylabel="z",
	xlabel="y",
	legendStyle="{at={(0.01,0.99)},anchor=north west}")

# Save to tikz format
PGF.save(joinpath(pwd(),
	"examples/direct_policy_optimization/figures/rocket_traj.tikz"), a, include_preamble=false)

# orientation tracking
p_nom_orientation = PGF.Plots.Linear(t_nom,hcat(x̄...)[3,:],
	mark="none",style="color=cyan, line width = 2pt",legendentry="TO")
p_nom_sim_orientation = PGF.Plots.Linear(t_sim_nom,hcat(z_lqr...)[3,:],
	mark="none",style="color=black, line width = 1pt")

p_sample_orientation = PGF.Plots.Linear(t_dpo,hcat(x...)[3,:],
	mark="none",style="color=orange, line width = 2pt",legendentry="DPO")
p_sample_sim_orientation = PGF.Plots.Linear(t_sim_dpo,hcat(z_dpo...)[3,:],
	mark="none",style="color=black, line width = 1pt")

a = Axis([p_nom_orientation;p_nom_sim_orientation],
    axisEqualImage=false,
    hideAxis=false,
	ylabel="orientation",
	xlabel="time",
	legendStyle="{at={(0.5,0.99)},anchor=north}")

# Save to tikz format
dir = joinpath(pwd(),
	"examples/direct_policy_optimization/figures/rocket_orientation_lqr.tikz")
PGF.save(dir, a, include_preamble=false)

a = Axis([p_sample_orientation;p_sample_sim_orientation],
    axisEqualImage=false,
    hideAxis=false,
	ylabel="orientation",
	xlabel="time",
	legendStyle="{at={(0.5,0.99)},anchor=north}")

# Save to tikz format
dir = joinpath(pwd(),
	"examples/direct_policy_optimization/figures/rocket_orientation_dpo.tikz")
PGF.save(dir, a, include_preamble=false)

# # eigen value analysis
# C = Diagonal(ones(model_sl.nx))[1:model_nom.nx,:]
# # nominal
# A_nom, B_nom = nominal_jacobians(model_nom,X_nom,U_nom,H_nom)
# A_nom_cl = [(A_nom[t] - B_nom[t]*K[t]) for t = 1:T-1]
# sv_nom = [norm.(eigen(A_nom_cl[t]).values) for t = 1:T-1]
# plt_nom = plot(hcat(sv_nom...)',xlabel="time step t",ylabel="eigen value norm",
# 	title="TVLQR nominal model",linetype=:steppost,
# 	ylims=(-3,3),labels="")
#
# # slosh nominal
# X_nom_slosh = [[copy(X_nom[t]);0.0;0.0] for t = 1:T]
# A_nom_slosh, B_nom_slosh = nominal_jacobians(model_sl,X_nom_slosh,U_nom,H_nom)
# A_nom_slosh_cl = [(A_nom_slosh[t] - B_nom_slosh[t]*K[t]*C) for t = 1:T-1]
# sv_nom_slosh = [norm.(eigen(A_nom_slosh_cl[t]).values) for t = 1:T-1]
# plt_nom_slosh = plot(hcat(sv_nom_slosh...)',xlabel="time step t",ylabel="eigen value norm",
# 	title="TVLQR slosh model",linetype=:steppost,
# 	ylims=(-3,3),labels="")
#
# # slosh
# A_dpo, B_dpo = nominal_jacobians(model_nom,X_nom_sample,U_nom_sample,H_nom_sample)
# A_dpo_cl = [(A_dpo[t] - B_dpo[t]*Θ_mat[t]) for t = 1:T-1]
# sv_dpo = [norm.(eigen(A_dpo_cl[t]).values) for t = 1:T-1]
# plt_dpo_nom = plot(hcat(sv_dpo...)',xlabel="time step t",ylabel="singular value",
# 	title="DPO nominal model",linetype=:steppost,
# 	ylims=(-1,3),labels="")
#
# X_dpo_slosh = [[copy(X_nom_sample[t]);0.0;0.0] for t = 1:T]
# A_dpo_slosh, B_dpo_slosh = nominal_jacobians(model_sl,X_dpo_slosh,U_nom_sample,H_nom_sample)
# A_dpo_slosh_cl = [(A_dpo_slosh[t] - B_dpo_slosh[t]*Θ_mat[t]*C) for t = 1:T-1]
# sv_dpo_slosh = [norm.(eigen(A_dpo_slosh_cl[t]).values) for t = 1:T-1]
# plt_dpo_slosh = plot(hcat(sv_dpo_slosh...)',xlabel="time step t",ylabel="eigen value norm",
# 	title="DPO slosh model",linetype=:steppost,
# 	ylims=(-1.0,3.0),labels="")
#
# plot(plt_nom,plt_dpo_nom,layout=(2,1))
#
# plot(plt_nom_slosh,plt_dpo_slosh,layout=(2,1))
#
#
# plot(t_nom,hcat(X_nom...)[1:6,:]',color=:cyan,label=["y" "z" "ϕ"])
# plot!(t_nom_sample, hcat(X_nom_sample...)[1:6,:]',color=:orange,label=["y" "z" "ϕ"])

# Animation plot
vis = Visualizer()
# render(vis)
open(vis)

pts_nom = collect(eachcol(hcat([[-p[1]; -1.0; p[2]] for p in x]...)))
material_nom = LineBasicMaterial(color = colorant"orange", linewidth = 10.0)
setobject!(vis["dpo_traj"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["dpo_traj"],true)

pts_nom = collect(eachcol(hcat([[-p[1]; 0.01; p[2]] for p in z_dpo]...)))
material_nom = LineBasicMaterial(color = colorant"grey", linewidth = 2.5)
setobject!(vis["dpo_sim"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["dpo_sim"],true)
# visualize_rocket!(vis, model_slosh, pad_trajectory(z_dpo, zeros(model_slosh.n), 100) , x1, xT,
# 	Δt = dt_sim_dpo, T_off = length(z_dpo) + 100)

pts_nom = collect(eachcol(hcat([[-p[1]; -1.1; p[2]] for p in x̄]...)))
material_nom = LineBasicMaterial(color = colorant"cyan", linewidth = 10.0)
setobject!(vis["to_traj"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["to_traj"],true)

pts_nom = collect(eachcol(hcat([[-p[1]; 0.01; p[2]] for p in z_lqr]...)))
material_nom = LineBasicMaterial(color = colorant"grey", linewidth = 2.5)
setobject!(vis["to_sim"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["to_sim"],true)
# visualize_rocket!(vis, model_slosh, pad_trajectory(z_lqr, zeros(model_slosh.n), 100) , x1, xT,
# 	Δt = dt_sim_nom, T_off = length(z_lqr) + 100)

x_ghost = [x[1], x[7], x[15], x[25], x[T]]
visualize_rocket_ghost!(vis, model_slosh, x_ghost)


text:
using GooglePlaces: callapi
using JSON3, Tables, DataFrames

function nearbysearch(
    ;
    key::String,
    location::Tuple,
    radius::Int64,
    rankby::String,
    language::String="en-GB",
    keyword::String="",
    name::String="",
    oftype::String="",
    pricerange::Tuple=(),
    pages::Int64=3
)
    @assert radius <= 50000 && radius >= 0 "accepted range for radius 0 => 50000"
    @assert rankby == "distance" || rankyby == "prominence" "accepted values for rankyby distance, promonence"
    @assert pages >= 1 && pages <= 3 "accepted values for pages 1 => 3"
    """
    if rankby == "distance"
        @assert !isempty(oftype) "value for oftype required when rankby==distance"
    end
    """

    urlbase = "https://maps.googleapis.com/maps/api/place/nearbysearch/json"
    urllocation = "?location=$(join(repr.(location), ","))"
    urlradius = "&radius=$radius"
    urlrankby = "&rankBy=$rankby"
    urllanguage = "&language=$language"
    urlkeyword = "&keyword=$keyword"
    urlname = "&name=$name"
    urltype = "&type=$oftype"
    urlkey = "&key=$key"

    url = (
        urlbase
        * urllocation
        * urlradius
        * urlrankby
        * urlkey
    )

    if !isempty(keyword)
        url *= urlkeyword
    end

    if !isempty(oftype)
        url *= urltype
    end

    if !isempty(name)
        url *= urlname
    end

    json_flattened = transformresults(url, pages)

    table = Tables.dictrowtable(json_flattened) # treat a json object of arrays or array of objects as a "table"
    df = DataFrame(table) # turn table into DataFrame

    return df
end

function transformresults(url::String, pages::Int64)
    pagecount = 0
    urlnextpagetoken = ""
    d = []

    while true
        response = String(callapi(url * urlnextpagetoken).body)
        json = JSON3.read(response) # read JSON string
        json_flattened = flatten_json.(json.results)
        append!(d, json_flattened)
        pagecount += 1

        "next_page_token" in keys(json) && pagecount < pages || break
        urlnextpagetoken = "&pagetoken=$(json["next_page_token"])"
    end
    return d
end


text:
@everywhere const roots_sin = Complex64[-2*pi,-pi,0.0,pi,2*pi]
@everywhere const colors_sin = ([1.0,1.0,0.0],[0.0,1.0,1.0],[1.0,0.0,1.0],[0.0,0.0,1.0],[1.0,0.0,0.0])

@everywhere function fsin(z0, roots, colors, maxiter, eps)
	n = 1
	z = z0 - sin(z0)/cos(z0)
	while n <= maxiter
		for j in 1:length(roots)
			if abs2(z - roots[j]) <= eps*eps
				return (1.0 - (n-1)/(maxiter-1))*colors[j]
			end
		end
		z -= sin(z)/cos(z)
		n += 1
	end
	return [0.0,0.0,0.0]
end

max_stars_repo_name: kalvotom/julia-examples

max_issues_repo_name: kalvotom/julia-examples

max_forks_repo_name: kalvotom/julia-examples


text:
##################################################################
"""
$(SIGNATURES)

Find the circumcenter of a triangle.

Derived from C source of Jonathan R Shewchuk <jrs@cs.cmu.edu>

Modified to return absolute coordinates.
"""
function tricircumcenter!(circumcenter,a,b,c)

    # Use coordinates relative to point `a' of the triangle.
    xba = b[1] - a[1]
    yba = b[2] - a[2]
    xca = c[1] - a[1]
    yca = c[2] - a[2]

    # Squares of lengths of the edges incident to `a'.
    balength = xba * xba + yba * yba
    calength = xca * xca + yca * yca

    # Calculate the denominator of the formulae.
    # if EXACT
    #    Use orient2d() from http://www.cs.cmu.edu/~quake/robust.html
    #    to ensure a correctly signed (and reasonably accurate) result,
    #    avoiding any possibility of division by zero.
    #  denominator = 0.5 / orient2d((double*) b, (double*) c, (double*) a)

    # Take your chances with floating-point roundoff
    denominator = 0.5 / (xba * yca - yba * xca)

    # Calculate offset (from `a') of circumcenter.
    xcirca = (yca * balength - yba * calength) * denominator
    ycirca = (xba * calength - xca * balength) * denominator

# The result is returned both in terms of x-y coordinates and xi-eta
# coordinates, relative to the triangle's point `a' (that is, `a' is
# the origin of both coordinate systems).	 Hence, the x-y coordinates
# returned are NOT absolute; one must add the coordinates of `a' to
# find the absolute coordinates of the circumcircle.  However, this means
# that the result is frequently more accurate than would be possible if
# absolute coordinates were returned, due to limited floating-point
# precision.  In general, the circumradius can be computed much more
# accurately.

    circumcenter[1] = xcirca+a[1]
    circumcenter[2] = ycirca+a[2]

    return circumcenter
end

"""
$(TYPEDEF)

Centers of voronoi cell facets (currently 1D, 2D).
"""
abstract type VoronoiFaceCenters <: AbstractGridFloatArray2D  end

function instantiate(grid::ExtendableGrid{Tc,Ti}, ::Type{VoronoiFaceCenters}) where {Tc,Ti}
    coord=grid[Coordinates]
    en=grid[EdgeNodes]
    nedges=size(en,2)
    dim=size(coord,1)
    xsigma=zeros(Tc,dim,nedges)

    if dim==1
        for iedge=1:nedges
            en1=en[1,iedge]
            en2=en[2,iedge]
            xsigma[1,iedge]=0.5*(coord[1,en1]+coord[1,en2])
        end
    elseif dim==2
        cn=grid[CellNodes]
        ec=grid[EdgeCells]
        cc1=zeros(2)
	cc2=zeros(2)
        for iedge=1:nedges
            icell1=ec[1,iedge]
            icell2=ec[2,iedge]

            # calculate voronoi face ends: either  cellcenter -- cellcenter or cellcenter --edgecenter
            @views tricircumcenter!(cc1,
                                    coord[:,cn[1,icell1]],
                                    coord[:,cn[2,icell1]],
                                    coord[:,cn[3,icell1]])
            if icell2>0
                @views tricircumcenter!(cc2,
                                        coord[:,cn[1,icell2]],
                                        coord[:,cn[2,icell2]],
                                        coord[:,cn[3,icell2]])

                @views xsigma[:,iedge].=0.5*(cc1+cc2)
            else
                inode1=en[1,iedge]
                inode2=en[2,iedge]
                @views cc2.=0.5*(coord[:,inode1]+coord[:,inode2])
                @views xsigma[:,iedge].=cc2
            end
            # We observe artifacts at domain corners if we use
            # this expression which would be consistent to the formal approach.
            # @views xsigma[:,iedge].=0.5*(cc1+cc2)
        end
    else
        error("3D Voronoi Face Centers not implemented yet")
    end
    xsigma
end


text:
# """
# Source Term Function
# ```julia
# source_term = SourceTerm(λ, μ)
# source_term(x, p)
# ```
# where x is a discritzed price point, source_params is a vector of parameters for
# the source term and p is the last "observed" mid price point
#
# The function returns a single scalar output
#
# """

mutable struct SourceTerm
    λ::Float64
    μ::Float64
end

function (st::SourceTerm)(x::Float64,p::Float64)
    return st.λ*(p-x)*exp(-st.μ*(p-x)^2)
end

# function (st::SourceTerm)(x::Float64,p::Float64)
#     return st.λ*tanh(st.μ*(p-x))
# end


text:
const BASEDIR="../input"
const DOCDIR="../doc"

const rules_url="https://rulebasedintegration.org/NotebookRuleFiles/Rubi-4.16.1.0.zip"
const doc_url="https://rulebasedintegration.org/PdfRuleFiles/Rubi-4.16.1.0.zip"

if Sys.isunix() && Sys.KERNEL != :FreeBSD
    function unpack_cmd(file, directory, extension, secondary_extension)
        if extension == ".zip"
            return (`unzip -x $file -d $directory`)
        end
        throw(ArguementError("I don't know how to unpack $file"))
    end
end

if Sys.KERNEL == :FreeBSD
    # The `tar` on FreeBSD can auto-detect the archive format via libarchive.
    # The supported formats can be found in libarchive-formats(5).
    # For NetBSD and OpenBSD, libarchive is not available.
    # For macOS, it is. But the previous unpack function works fine already.
    function unpack_cmd(file, dir, ext, secondary_ext)
        tar_args = ["--no-same-owner", "--no-same-permissions"]
        return pipeline(
            `/bin/mkdir -p $dir`,
            `/usr/bin/tar -xf $file -C $dir $tar_args`)
    end
end

if Sys.iswindows()
    const exe7z = joinpath(Sys.BINDIR, "7z.exe")

    function unpack_cmd(file,directory,extension,secondary_extension)
        if (extension == ".zip" || extension== ".gz" || extension == ".7z" || extension == ".tar" ||
                (extension == ".exe" && secondary_extension == ".7z"))
            return (`$exe7z x $file -y -o$directory`)
        end
        throw(ArguementError("I don't know how to unpack $file"))
    end
end

"""
    unpack(f; keep_originals=false)
Extracts the content of an archive in the current directory;
deleting the original archive, unless the `keep_originals` flag is set.
"""
function unpack(f; keep_originals=false)
    run(unpack_cmd(f, pwd(), last(splitext(f)), last(splitext(first(splitext(f))))))
    !keep_originals && rm(f)
end

function fix_names(root, names)
    for name in names
        old_name = joinpath(root, name)
        new_name = joinpath(root, replace(name, ' ' => '_'))
        if old_name != new_name
            mv(old_name, new_name)
        end
    end
end

function fix_names(dir)
    for (root, dirs, files) in walkdir(dir, topdown=false)
        fix_names(root, [dirs; files])
    end
end

function download_rules()
    rm(BASEDIR, recursive=true, force=true)
    mkdir(BASEDIR)
    cd(BASEDIR)
    filename = splitpath(rules_url)[end]
    download(rules_url, filename)
    unpack(filename)
    # cd("..")
end

function download_doc()
    rm(DOCDIR, recursive=true, force=true)
    mkdir(DOCDIR)
    cd(DOCDIR)
    filename = splitpath(doc_url)[end]
    download(doc_url, filename)
    unpack(filename)
    # cd("..")
end

download_rules()
download_doc()

fix_names(BASEDIR)
fix_names(DOCDIR)


text:
"""
    isna(x)

Indicate whether `x` is `missing` or `NaN`.
"""
isna(x::Any) = isnan(x) ? true : false
isna(::Missing) = true

# """
#     skipnas(itr)

# Return an iterator over the elements in `itr` skipping all non-applicable values (`missing` and `NaN`). The returned object can be indexed using indices of `itr`.
# """
# skipnas(itr) = SkipNAs(itr)

# struct SkipNAs{T}
#     x::T
# end
# IteratorSize(::Type{<:SkipNAs}) = SizeUnknown()
# IteratorEltype(::Type{SkipNAs{T}}) where {T} = IteratorEltype(T)
# eltype(::Type{SkipNAs{T}}) where {T} = nonmissingtype(eltype(T))

# function iterate(itr::SkipNAs, state...)
#     item, state = iterate(itr.x, state...)
#     while isna(item)
#         item, state = iterate(itr.x, state)
#     end
#     item, state
# end

# IndexStyle(::Type{<:SkipNAs{T}}) where {T} = IndexStyle(T)
# eachindex(itr::SkipNAs) =
#     Iterators.filter(i -> !isna(@inbounds(itr.x[i])), eachindex(itr.x))
# keys(itr::SkipNAs) =
#     Iterators.filter(i -> !isna(@inbounds(itr.x[i])), keys(itr.x))
# Base.@propagate_inbounds function getindex(itr::SkipNAs, I...)
#     v = itr.x[I...]
#     isna(v) && throw(MissingException("the value at index $I is missing or NaN"))
#     v
# end

"""
    lchoose(n, j)

    Return the log of the j-combination of a set that has n elements.
"""
function lchoose(n::Int, j::Int)
    return loggamma(n+1) - loggamma(j+1) - loggamma(n-j+1)
end


text:
zero(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat{T}(zero(T), zero(T))
one(::Type{DoubleFloat{T}}) where {T<:AbstractFloat} = DoubleFloat{T}(one(T), zero(T))

nan(::Type{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat{T}(T(NaN), T(NaN))
inf(::Type{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat{T}(T(Inf), T(NaN))
posinf(::Type{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat{T}(T(Inf), T(NaN))
neginf(::Type{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat{T}(T(-Inf), T(NaN))

typemax(::Type{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat{T}(typemax(T))
typemin(::Type{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat{T}(typemin(T))

floatmin(::Type{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat{T}(floatmin(T))

floatmax(::Type{DoubleFloat{Float64}}) = Double64(ldexp(Float64(0.9999999999999999), 1024), ldexp(Float64(0.9999999999999999), 970))
floatmax(::Type{DoubleFloat{Float32}}) = Double32(ldexp(Float32(0.99999994f0), 128), ldexp(Float32(0.99999994f0), 103))
floatmax(::Type{DoubleFloat{Float16}}) = Double16(ldexp(Float16(0.9995),16), ldexp(Float16(0.9996),4))

maxintfloat(::Type{DoubleFloat{Float64}}) = DoubleFloat{Float64}(ldexp(0.5, 107), ldexp(0.5, 54))
maxintfloat(::Type{DoubleFloat{Float32}}) = DoubleFloat{Float32}(ldexp(0.5f0, 49), ldexp(0.5f0, 25))
maxintfloat(::Type{DoubleFloat{Float16}}) = DoubleFloat{Float16}(65519.0)

isnan(x::DoubleFloat{T}) where {T<:IEEEFloat} = isnan(HI(x))
isinf(x::DoubleFloat{T}) where {T<:IEEEFloat} = isinf(HI(x))
isposinf(x::DoubleFloat{T}) where {T<:IEEEFloat} = isinf(HI(x)) && !signbit(HI(x))
isneginf(x::DoubleFloat{T}) where {T<:IEEEFloat} = isinf(HI(x)) && signbit(HI(x))


text:
using Documenter
using ContextVariablesX

makedocs(
    sitename = "ContextVariablesX",
    format = Documenter.HTML(),
    modules = [ContextVariablesX],
)

deploydocs(; repo = "github.com/tkf/ContextVariablesX.jl", push_preview = true)


text:
function larger(x, y)
    if (x>y)
        return x
    end
    return y
end
println(larger(7,8))

ismatch(r"^\([0-9]{3}\)[0-9]{3}-[0-9]{4}$", "(781)244-1212")

ismatch(r"^\([0-9]{3}\)[0-9]{3}-[0-9]{4}$", "-781-244-1212")

using FactCheck

f(x) = x^3
facts("cubes") do
    @fact f(2) --> 8
    @fact f(2) --> 7
end

using Base.Test
F(x) = x^3
@test F(2) == 8
@test F(2) == 7

max_stars_repo_path: Section 4/Control+Flow.jl

max_issues_repo_path: Section 4/Control+Flow.jl

max_forks_repo_path: Section 4/Control+Flow.jl


text:
export snr, smooth, smooth!, abs_max, abs_max!, standardize, standardize!
export mad, std_threshold

"""
    snr(A)

Signal to noise ratio of cross-correlations in matrix `A`.

Follows method of Clarke et. al, 2011. Measures SNR at each point.
"""
function snr(A::AbstractArray)
    Nrows,Ncols = size(A)
    A_mean = mean(A,dims=2)

    # calculate noise and envelope functions
    sigma = mean(A.^2,dims=2) .- A_mean.^2
    sigma .= sqrt.(sigma./ (Ncols-1))
    s = abs.(A_mean .+ im .* hilbert(A_mean))
    return s ./ sigma
end
snr(C::CorrData) = snr(C.corr)

"""
    smooth(A, half_win)

Smooth array `A` with half-window `half_win` (defaults to 3).
"""
function smooth!(A::AbstractArray, half_win::Int=3, dims::Int=1)
    T = eltype(A)
    window_len = 2 * half_win + 1
    csumsize = tuple(collect(size(A)) .+ [i==1 ? 2 * (window_len - 1) + 1 : 0 for i in 1:ndims(A)]...)
    csum = similar(A,T,csumsize)
    csum[1:window_len,:] .= zero(T)
    csum[end - window_len + 1:end,:] .= zero(T)
    csum[window_len+1:end-window_len + 1,:] .= A
    csum .= cumsum(csum,dims=dims)
    weight = similar(A,T,size(A,1))
    weight[1:half_win] = T.(window_len ÷ 2 + 1:window_len - 1)
    weight[half_win + 1: end - half_win] .= T(window_len)
    weight[end-half_win:end] = T.(window_len:-1:window_len ÷ 2 + 1)
    A[:,:] .= (csum[window_len+half_win+1:end-half_win,:] .- csum[half_win+1:end-window_len-half_win,:]) ./ weight
    return nothing
end
smooth(A::AbstractArray,half_win::Int=3, dims::Int=1) =
      (U = deepcopy(A);smooth!(U,half_win,dims);return U)

"""
    abs_max!(A)

Returns array `A` divided by its absolute maximum value.
"""
function abs_max!(A::AbstractArray)
    maxabs = maximum(abs.(A),dims=1)
    if any(maxabs .== 0)
        throw(DomainError("All zero column leads to NaN"))
    end
    A ./= maxabs
    return nothing
end
abs_max(A::AbstractArray) = (U = deepcopy(A);abs_max!(U);return U)
abs_max!(C::CorrData) = abs_max!(C.corr)
abs_max(C::CorrData) = (U = deepcopy(C);abs_max!(U);return U)

"""
    standardize!(A)

Demean and standardize array `A` to unit std.
"""
function standardize!(A::AbstractArray)
    A .-= mean(A,dims=1)
    A ./= std(A,dims=1)
    return nothing
end
standardize(A::AbstractArray) = (U = deepcopy(A);standardize!(U);return U)
standardize!(C::CorrData) = standardize!(C.corr)
standardize(C::CorrData) = (U = deepcopy(C);standardize!(U);return U)

"""
    mad(A)

Median Absolute Deviation of array `A`.
MAD = median(|Xi- median(A)|)
"""
function mad(A::AbstractArray)
    return median(abs.(A .- median(A,dims=1)),dims=1)
end
mad(C::CorrData) = mad(C.corr)

"""

  std_threshold(A,thresh)

Returns indices of cols of `A` where max(abs.A[:,col]))/std(A[:,col]) < `max_std`.
"""
function std_threshold(A::AbstractArray, max_std::Real)
    stds = std(A,dims=1)[1,:]
    maxs = maximum(abs.(A),dims=1)[1,:]
    threshs = maxs ./ stds
    ind = []
    for ii = 1:length(threshs)
        if threshs[ii] < max_std
            append!(ind,ii)
        end
    end
    return ind
end


text:
# abstract type declarations

abstract type AbstractEnsemble end

abstract type AbstractBins end

abstract type AbstractWEsampler end


text:
# This file is a part of ArraysOfArrays.jl, licensed under the MIT License (MIT).

Base.sum(X::AbstractVectorOfSimilarArrays{T,M}, w::StatsBase.AbstractWeights) where {T,M} = sum(flatview(X), w, M + 1)

Statistics.mean(X::AbstractVectorOfSimilarArrays{T,M}, w::StatsBase.AbstractWeights) where {T,M} =
    vec(mean(flatview(X), w, dims = M + 1))

Statistics.var(X::AbstractVectorOfSimilarArrays{T,M}, w::StatsBase.AbstractWeights; corrected::Bool = true) where {T,M} =
    vec(var(flatview(X), w, M + 1; corrected = corrected))

Statistics.cov(X::AbstractVectorOfSimilarVectors, w::StatsBase.AbstractWeights; corrected::Bool = true) =
    cov(flatview(X), w, 2; corrected = corrected)

Statistics.cor(X::AbstractVectorOfSimilarVectors, w::StatsBase.AbstractWeights) =
    cor(flatview(X), w, 2)


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

# Run the distributed test outside of the main driver since it needs its own
# set of dedicated workers.

inline_flag = Base.JLOptions().can_inline == 1 ? `` : `--inline=no`
cov_flag = ``
if Base.JLOptions().code_coverage == 1
    cov_flag = `--code-coverage=user`
elseif Base.JLOptions().code_coverage == 2
    cov_flag = `--code-coverage=all`
end

cmd = `$(Base.julia_cmd()) $inline_flag $cov_flag --check-bounds=yes --startup-file=no --depwarn=error distributed_exec.jl`

if !success(pipeline(cmd; stdout=STDOUT, stderr=STDERR)) && ccall(:jl_running_on_valgrind,Cint,()) == 0
    error("Distributed test failed, cmd : $cmd")
end


text:
include("tools.jl")
include("loadFiles.jl")

df_classes."Matched"   .= ""
df_classes."Unmatched" .= ""
# Testing class matches
#counter = 1
# Iterate over all rows in faculty sheet:
for row in eachrow(df_faculty)
    println("#######################################################################################")
    println(row["Your Name"])
    # Look for possible courses:
    for course in keys(FacTA_class)
        #filteredClass = dropmissing(df_faculty,FacTA_class[course])
        if ismissing(row[FacTA_class[course]]) == false
            println(row[FacTA_class[course]])
            ind_matched = findall(df_classes."Class Name".==row[FacTA_class[course]])[1]
            println(ind_matched)
            df_classes[ind_matched,3] = "X"
            println("Selectable candidates (chosen by faculty): ", row[FacTA_student[course]])
            students = split(row[FacTA_student[course]],",")
            # @show students
            println("Student Choices (Name, PriorityRank)")
            for a in keys(StuClass)
                matchedStudent = filter(df -> df[StuClass[a]]  == row[FacTA_class[course]], dropmissing(df_students, StuClass[a]))
                for i=1:size(matchedStudent,1)
                    sName = matchedStudent[!,:"Your Name"][i]
                    matches, _ = matchNames(sName,students)
                    @show sName
                    if matches>1
                        df_classes[ind_matched,4] = df_classes[ind_matched,4] *  sName * ", "
                        println(sName, " ", matchedStudent[!,StuInterest[a]][i], ", NameMatches: ", matches)
                    else
                        df_classes[ind_matched,5] = df_classes[ind_matched,5]  * sName *  ", "
                    end
                end
            end
            println("____________________________________________________________________________________________")
            #counter +=1
        end
    end
end
outFile = "/Users/cfranken/GDrive/work/Caltech/OptionRepWork/TA2021/testMatch.xlsx"
XLSX.writetable(outFile, collect(DataFrames.eachcol(df_classes)), DataFrames.names(df_classes))

max_stars_repo_name: cfranken/TA-matching

max_issues_repo_name: cfranken/TA-matching

max_forks_repo_name: cfranken/TA-matching


text:
using StaticArrays, LightGraphs, MetaGraphs,
LazySets, Random, POMDPs, Combinatorics,
Parameters, POMDPModelTools, Distributions, StaticArrays
include("ndgrid.jl")

import POMDPs: actions, states, discount, isterminal, generate_s, generate_o, reward
import POMDPModelTools: obs_weight

const Vec2 = SVector{2, Int64}
const AGENT_ACTIONS = [Vec2(i, j) for i in -1:1 for j in -1:1]
const SENSOR_MODEL = Dict(:o1s1 => BoolDistribution(0.9),
                          :o1s0 => BoolDistribution(0.1))
# Just an alias
const Region = Vector{Vec2}

############################################
#########  Environment -related  ###########
############################################
struct Environment
    size::NTuple{2, Int}
    obstacles::Vector{<:Hyperrectangle}
    visibility::Vector{Vector{Tuple}}
    G::MetaGraph
    pos_to_ind::Dict{Vec2, Int64}
    R::Int64

    function Environment(size::NTuple{2, Int}, obstacles::Vector{<:Hyperrectangle}, R::Int64)
        G = create_graph(size...)
        remove_obstacles(G, obstacles)
        visibility = limited_visible(G, obstacles, R, inf_visible(G, obstacles))
        node_mapping = Dict(get_pos(G, i) => i for i in 1:nv(G))
        return new(size, obstacles, visibility, G, node_mapping, R)
    end
end

Base.size(env::Environment) = env.size
graph(env::Environment) = env.G
obstacles(env::Environment) = env.obstacles
visibility(env::Environment) = env.visibility
visibility(env::Environment, ind::Int64) = env.visibility[ind]
visibility(env::Environment, pos::Vec2)  = env.visibility[pos_to_ind(env, pos)]

pos_to_ind(env::Environment, pos::Vec2) = env.pos_to_ind[pos]
get_pos(G::MetaGraph, i::Int64) = Vec2(get_pos_tup(G, i))

in_obstacle(env::Environment, x) = in_obstacle(env::Environment, Vector{Float64}(x))
in_obstacle(env::Environment, x::Vector{Float64}) = any(x ∈ o for o in obstacles(env))

Base.rand(G::AbstractGraph, rng::AbstractRNG = Random.GLOBAL_RNG) = rand(1:nv(G))
function rand_pos(G::AbstractGraph, rng::AbstractRNG = Random.GLOBAL_RNG)
    i = rand(G, rng)
    Vec2(get_pos_tup(G, i))
end
rand_pos(G::AbstractGraph) = rand_pos(G, Random.GLOBAL_RNG)

function rand_pos(G::AbstractGraph, n::Integer, rng::AbstractRNG = Random.GLOBAL_RNG; replace = false)
    points = Set{Vec2}()
    while length(points) < n
        push!(points, rand_pos(G, rng))
    end
    return collect(points)
end

####################################################
#################### S, A, O #######################
####################################################
struct Agent
    pos::Vec2
    partition::Region
end

partition(ag::Agent) = ag.partition # not sure how much this will be used.
pos(ag::Agent)       = ag.pos
Agent(ag::Agent, pos::Vec2) = Agent(pos, ag.partition)
Base.:+(ag::Agent, a::Vec2) = Agent(ag, ag.pos+a)
isvalid_position(s::Agent, pos::Vec2) = pos ∈ partition(s) && !in_obstacle(env(p), pos)
isvalid_position(pos::Vec2) = !in_obstacle(env(p), pos)

struct SNState
    patrollers::Vector{Agent}
    adversaries::Vector{Vec2}
end

struct Alarm
    guess::Vec2
end

SNAction = Union{Vector{Vec2}, Alarm}
SNObservation = Vector{Pair{Vec2, Bool}}

agents(s::SNState)      = (s.patrollers, s.adversaries)
patrollers(s::SNState)  = s.patrollers
adversaries(s::SNState) = s.adversaries
# occupied_probability(o::SNObservation, pos::Vec2) = o.occupied[pos]

####################################################
################ Adversary related #################
####################################################

# need to define these before the pomdp type, even
# though their implementations only happen later on
abstract type AdversaryModel end

struct RandomActions <: AdversaryModel end
struct MaximizeDistance <: AdversaryModel end
struct MaximizeSafety <: AdversaryModel end
# struct Omniscient{APB<:AbstractParticleBelief} <: AdversaryModel
#     belief::APB
# end
# There was another one wasn't there?

####################################################
############# POMDP type + functions ###############
####################################################
"""
S :: vector of agents and pos of Adversary/s
A :: a rearrangement of the agents or an Alarm

O :: A boolean for each node being observed

R ::
    seeing Adversary (closer is better),
    minimize uncertainty,
    maximize overall overage,
    maximize agent utilization
"""
@with_kw mutable struct SkyNet <: POMDP{SNState, SNAction, SNObservation}
    size::NTuple{2, Int64}         = (20, 20)
    R::Int64                       = 3
    n_agents::Int8                 = 5
    discount_factor::Float64       = 0.9
    env::Environment               =    begin
                                            obs = [Hyperrectangle(low=[3, 3],  high=[6, 8]),
                                                   Hyperrectangle(low=[4, 9],  high=[8, 10]),
                                                   Hyperrectangle(low=[15,15], high=[16,19]),
                                                   Hyperrectangle(low=[16,3],  high=[18,10]),
                                                   Hyperrectangle(low=[3,16],  high=[10,18])]
                                            Environment(size, obs, R)
                                        end
    partition_points::Vector{Vec2} = [CartesianIndex((2, 3)),
                                       CartesianIndex((9, 2)),
                                       CartesianIndex((5,13)),
                                       CartesianIndex((14,11)),
                                       CartesianIndex((13,7))]
                                            # rand_pos(graph(env), n_agents) # maybe don't need
    partitions::Vector{Region}     =    begin
                                            voronoi_regions = voronoi_cells(graph(env), obstacles(env), partition_points)
                                            Region.(voronoi_regions)
                                        end
    adversarial_model              = RandomActions()
    # max_adversaries::Int8               = 1      # maybe don't need
    # n_adversaries::Int8                 = 1      # maybe don't need
    # intrusion_prob::Float64             = 1.0    # maybe don't need
    # sensor_distr::Truncated{Normal{Float64},Continuous} = Truncated(Normal(0, -R, R))
end

environment(p::SkyNet)       = p.env
env(p::SkyNet)               = p.env
graph(p::SkyNet)             = p.env.G
obstacles(p::SkyNet)         = p.env.obstacles
adversarial_model(p::SkyNet) = p.adversarial_model
discount(p::SkyNet)          = p.discount_factor

"""
    isterminal checks if the nodes are outside the environment,
    which is used as a flag for being terminal. Should only happen
    if they were intentionally set that way at the end of generate_s,
    if they chose to sound the alarm.
"""
isterminal(p::SkyNet, s::SNState) = any(pat->pos(pat) > Vec2(size(env(p))), patrollers(s))

####################################################
############### More Adversary stuff ###############
####################################################

heuristic_update(p, s, rng) = heuristic_update(adversarial_model(p), p, s, rng)
function heuristic_update(::MaximizeDistance, p::SkyNet, s::SNState, rng)
end
function heuristic_update(::RandomActions, p::SkyNet, s::SNState, rng)
    G = graph(p)
    newadv = similar(adversaries(s))
    for (i, adv) in enumerate(adversaries(s))
        adv2 = adv + rand(rng, AGENT_ACTIONS)
        if all(adv2 .> 0) && all(adv2 .<= Vec2(size(G))) && !in_obstacle(env(p), adv2)
            newadv[i] = adv2
        else
            newadv[i] = adv
        end
    end
    return newadv
end

function actions(p::SkyNet, s::SNState)
    agentwise_action_space = actions.([p], patrollers(s))
    moves = tuple.(vec.(ndgrid(agentwise_action_space...))...)
    alarms = alarm_actions(p, s)
    A = SNAction[collect.(moves); alarms]
    println("A size: $(length(A))")
    println("$(length(moves)) moves")
    println("$(length(alarms)) alarms")
    A
end

function actions(p::SkyNet, s::Agent)
    x0 = pos(s)
    A = filter(dir->isvalid_position(s, x0+dir), AGENT_ACTIONS)
end

function alarm_actions(p::SkyNet, s::SNState)
    env = environment(p)
    # can only ring an alarm of a cell you can see
    agent_inds    = pos_to_ind.([env], pos.(patrollers(s))) # hacky broadcast
    visible_cells = vcat(visibility.([env], vcat(agent_inds...))...)
    Alarm.(Set(visible_cells))
end

function generate_s(p::SkyNet, s::SNState, a::SNAction, rng::AbstractRNG)

    # @show a

    if isterminal(p, s)
        error("State is terminal in generate_s. Why is this happening? \ns = $s, \n a = $a")
        # return s
    end
    if a isa Alarm
        out_of_bounds_pats = patrollers(s) .+ [Vec2(size(env(p)))]
        return SNState(out_of_bounds_pats, adversaries(s))
    end

    pats, advs = agents(s)
    updated_pats = pats .+ a
    # update adversary's position somehow:
    new_adv = heuristic_update(p, s, rng)

    return SNState(updated_pats, new_adv)
end

function obs_weight(p::SkyNet, s::SNState, a::SNAction, sp::SNState, o::SNObservation)

    advs = adversaries(s)
    prob = 1.0
    # if a isa Alarm
        # return 1.0 * (a.guess ∈ advs)
        # return 0.5
    # end

    for oᵢ in o
        pos, val = oᵢ
        isin = pos ∈ advs ? :o1s1 : :o1s0
        p_accurate = pdf(SENSOR_MODEL[isin], val)

        prob *= p_accurate
    end
    # @show prob
    return prob
end

function generate_o(p::SkyNet, s::SNState, a::SNAction, sp::SNState, rng::AbstractRNG)
    # println(pos.(patrollers))
    # println("adv: ", adversaries)
    # println(a)
    o = Vector{Pair{Vec2, Bool}}()
    if a isa Alarm
        # println(a)
        return o
        # sp = s
    end
    patrollers, adversaries = agents(sp)
    for pat in patrollers
        visibles = visibility(env(p), pos(pat))
        for v in Vec2.(visibles)
            positive_rate = v ∈ adversaries ? :o1s1 : :o1s0
            detected = rand(rng, SENSOR_MODEL[positive_rate])
            push!(o, v=>detected)
        end
    end
    return o
end

function reward(p::SkyNet, s::SNState, a::SNAction, sp::SNState)
    # subtract some points for each adversary in the scene.
    # @show typeof(a)
    r = -20.0 * length(adversaries(s))
    if a isa Alarm
        # Alarms are very costly.
        r -= 10000.0
        dists = norm.([a.guess] .- adversaries(s))
        # @show dists
        if any(dists .== 0)
            r += 600.0
        elseif any(dists .<= 1.415) # √2 is diagonal distance
            r += 300.0
        end
    else
        travelled = norm.(a)  # gives, 0, 1, or √2 depeding on direction of travel
        r -= 10*sum(travelled)
        # @show r
        # old_distance_to_adversary = norm.(pos.(patrollers(s)) .- adversaries(s))
        # new_distance_to_adversary = norm.(pos.(patrollers(sp)) .- adversaries(sp))
        # r += 200*maximum(new_distance_to_adversary .- old_distance_to_adversary)
    end
    return r
end


text:
developing = false; # are you developing the package or just using?
if developing
	# ] add Revise # if you don't have it
    path = pwd()
    if path[end-6:end] != "MatLang"
        error("cd(to the MatLang path)")
    end
    push!(LOAD_PATH, ".")
    using Revise
    using MatLang
else
    using MatLang
end
################################################################
################################################################
nFloat64 = doubleM(12.0)
t1 = typeof(nFloat64) # Float64

nFloat32 = singleM(12.0)
t2 = typeof(nFloat32) # Float32

nInt8 = int8M(-10)
t3 = typeof(nInt8) # Int8

nInt16 = int16M(-10)
t4 = typeof(nInt16) # Int16

nInt32 = int32M(-10)
t5 = typeof(nInt32) # Int32

nInt64 = int64M(-10)
t6 = typeof(nInt64) # Int64

nUInt8 = uint8M(10)
t7 = typeof(nUInt8) # UInt8

nUInt16 = uint16M(10)
t8 = typeof(nUInt16) # UInt16

nUInt32 = uint32M(10)
t9 = typeof(nUInt32) # UInt32

nUInt64 = uint64M(10)
t10 = typeof(nUInt64) # UInt64
################################################################


text:
@testset "switch combinator" begin

    # ------------ Trace ------------ #

    @gen function swtrg()
        z ~ normal(3.0, 5.0)
        return z
    end

    @testset "switch trace" begin
        tr = simulate(swtrg, ())
        swtr = Gen.SwitchTrace(swtrg, tr, get_retval(tr), (1, ), get_score(tr), 0.0)
        @test swtr[:z] == tr[:z]
        @test project(swtr, AllSelection()) == project(swtr.branch, AllSelection())
        @test project(swtr, EmptySelection()) == swtr.noise
    end

    # ------------ Bare combinator ------------ #

    # Model chunk.
    @gen (grad) function bang0((grad)(x::Float64), (grad)(y::Float64))
        std::Float64 = 3.0
        z = @trace(normal(x + y, std), :z)
        return z
    end

    @gen (grad) function fuzz0((grad)(x::Float64), (grad)(y::Float64))
        std::Float64 = 3.0
        z = @trace(normal(x + 2 * y, std), :z)
        return z
    end
    sc = Switch(bang0, fuzz0)
    # ----.

    @testset "simulate" begin
        tr = simulate(sc, (1, 5.0, 3.0))
        @test isapprox(get_score(tr), logpdf(normal, tr[:z], 5.0 + 3.0, 3.0))
        tr = simulate(sc, (2, 5.0, 3.0))
        @test isapprox(get_score(tr), logpdf(normal, tr[:z], 5.0 + 2 * 3.0, 3.0))
    end

    @testset "generate" begin
        chm = choicemap()
        chm[:z] = 5.0
        tr, w = generate(sc, (2, 5.0, 3.0), chm)
        assignment = get_choices(tr)
        @test assignment[:z] == 5.0
        @test isapprox(w, logpdf(normal, 5.0, 5.0 + 2 * 3.0, 3.0))
    end

    @testset "assess" begin
        chm = choicemap()
        chm[:z] = 5.0
        w, ret = assess(sc, (2, 5.0, 3.0), chm)
        @test isapprox(w, logpdf(normal, 5.0, 5.0 + 2 * 3.0, 3.0))
    end

    @testset "propose" begin
        chm, w = propose(sc, (2, 5.0, 3.0))
        @test isapprox(w, logpdf(normal, chm[:z], 5.0 + 2 * 3.0, 3.0))
    end

    @testset "update" begin
        tr = simulate(sc, (1, 5.0, 3.0))
        old_sc = get_score(tr)
        chm = choicemap((:x => :z, 5.0))
        new_tr, w, rd, discard = update(tr, (2, 5.0, 3.0),
                                        (UnknownChange(), NoChange(), NoChange()),
                                        chm)
        @test isapprox(old_sc, get_score(new_tr) - w)
        chm = choicemap((:x => :z, 10.0))
        new_tr, w, rd, discard = update(tr, (1, 5.0, 3.0),
                                        (UnknownChange(), NoChange(), NoChange()),
                                        chm)
        @test isapprox(old_sc, get_score(new_tr) - w)
    end

    @testset "regenerate" begin
        tr = simulate(sc, (2, 5.0, 3.0))
        old_sc = get_score(tr)
        sel = select(:z)
        new_tr, w, rd = regenerate(tr, (2, 5.0, 3.0),
                                   (UnknownChange(), NoChange(), NoChange()),
                                   sel)
        @test isapprox(old_sc, get_score(new_tr) - w)
        new_tr, w, rd = regenerate(tr, (1, 5.0, 3.0),
                                   (UnknownChange(), NoChange(), NoChange()),
                                   sel)
        @test isapprox(old_sc, get_score(new_tr) - w)
        args = map(x -> 2 * x, get_args(new_tr))
        argdiffs = map(_ -> UnknownChange(), args)
        old_sc = get_score(new_tr)
        new_tr, w, rd = regenerate(new_tr, args, argdiffs, sel)
        @test isapprox(old_sc, get_score(new_tr) - w)
    end

    @testset "choice gradients" begin
        for z in [1.0, 3.0, 5.0, 10.0]
            chm = choicemap((:z, z))
            tr, _ = generate(sc, (1, 5.0, 3.0), chm)
            sel = select(:z)
            input_grads, choices, gradients = choice_gradients(tr, sel)
            expected_choice_grad = logpdf_grad(normal, z, 5.0 + 3.0, 3.0)
            @test isapprox(gradients[:z], expected_choice_grad[1])
            tr, _ = generate(sc, (2, 5.0, 3.0), chm)
            input_grads, choices, gradients = choice_gradients(tr, sel)
            expected_choice_grad = logpdf_grad(normal, z, 5.0 + 2 * 3.0, 3.0)
            @test isapprox(gradients[:z], expected_choice_grad[1])
        end
    end

    # ------------ Hierarchy ------------ #

    # Model chunk.
    @gen (grad) function bang1((grad)(x::Float64), (grad)(y::Float64))
        @param(std::Float64)
        z = @trace(normal(x + y, std), :z)
        return z
    end
    init_param!(bang1, :std, 3.0)
    @gen (grad) function fuzz1((grad)(x::Float64), (grad)(y::Float64))
        @param(std::Float64)
        z = @trace(normal(x + 2 * y, std), :z)
        return z
    end
    init_param!(fuzz1, :std, 3.0)
    sc = Switch(bang1, fuzz1)
    @gen (grad) function bam(s::Int)
        x ~ sc(s, 5.0, 3.0)
        return x
    end
    # ----.

    @testset "simulate" begin
        tr = simulate(bam, (2, ))
        @test isapprox(get_score(tr), logpdf(normal, tr[:x => :z], 5.0 + 2 * 3.0, 3.0))
    end

    @testset "generate" begin
        chm = choicemap()
        chm[:x => :z] = 5.0
        tr, w = generate(bam, (2, ), chm)
        assignment = get_choices(tr)
        @test assignment[:x => :z] == 5.0
        @test isapprox(w, logpdf(normal, 5.0, 5.0 + 2 * 3.0, 3.0))
    end

    @testset "assess" begin
        chm = choicemap()
        chm[:x => :z] = 5.0
        w, ret = assess(bam, (2, ), chm)
        @test isapprox(w, logpdf(normal, 5.0, 5.0 + 2 * 3.0, 3.0))
    end

    @testset "propose" begin
        chm, w = propose(bam, (2, ))
        @test isapprox(w, logpdf(normal, chm[:x => :z], 5.0 + 2 * 3.0, 3.0))
    end

    @testset "update" begin
        tr = simulate(bam, (2, ))
        old_sc = get_score(tr)
        chm = choicemap((:x => :z, 5.0))
        new_tr, w, rd, discard = update(tr, (2, ), (UnknownChange(), ), chm)
        @test isapprox(old_sc, get_score(new_tr) - w)
        chm = choicemap((:x => :z, 10.0))
        new_tr, w, rd, discard = update(tr, (1, ), (UnknownChange(), ), chm)
        @test isapprox(old_sc, get_score(new_tr) - w)
    end

    @testset "regenerate" begin
        tr = simulate(bam, (2, ))
        old_sc = get_score(tr)
        new_tr, w = regenerate(tr, (1, ), (UnknownChange(), ), select())
        @test isapprox(old_sc, get_score(new_tr) - w)
        new_tr, w = regenerate(tr, (2, ), (UnknownChange(), ), select())
        @test isapprox(old_sc, get_score(new_tr) - w)
        new_tr, w = regenerate(tr, (1, ), (UnknownChange(), ), select())
        @test isapprox(old_sc, get_score(new_tr) - w)
    end

    @testset "choice gradients" begin
        for z in [1.0, 3.0, 5.0, 10.0]
            chm = choicemap((:x => :z, z))
            tr, _ = generate(bam, (1, ), chm)
            sel = select(:x => :z)
            input_grads, choices, gradients = choice_gradients(tr, sel)
            expected_choice_grad = logpdf_grad(normal, z, 5.0 + 3.0, 3.0)
            @test isapprox(gradients[:x => :z], expected_choice_grad[1])
            chm = choicemap((:x => :z, z))
            tr, _ = generate(bam, (2, ), chm)
            sel = select(:x => :z)
            input_grads, choices, gradients = choice_gradients(tr, sel)
            expected_choice_grad = logpdf_grad(normal, z, 5.0 + 2 * 3.0, 3.0)
            @test isapprox(gradients[:x => :z], expected_choice_grad[1])
        end
    end

    @testset "accumulate parameter gradients" begin
        for z in [1.0, 3.0, 5.0, 10.0]
            chm = choicemap((:z, z))
            tr, _ = generate(bam, (1, ), chm)
            zero_param_grad!(bang1, :std)
            input_grads = accumulate_param_gradients!(tr, 1.0)
            expected_std_grad = logpdf_grad(normal, tr[:x => :z], 5.0 + 3.0, 3.0)[3]
            @test isapprox(get_param_grad(bang1, :std), expected_std_grad)
            tr, _ = generate(bam, (2, ), chm)
            zero_param_grad!(fuzz1, :std)
            input_grads = accumulate_param_gradients!(tr, 1.0)
            expected_std_grad = logpdf_grad(normal, tr[:x => :z], 5.0 + 2 * 3.0, 3.0)[3]
            @test isapprox(get_param_grad(fuzz1, :std), expected_std_grad)
        end
    end

    # ------------ (More complex) hierarchy ------------ #

    # Model chunk.
    @gen (grad) function bang2((grad)(x::Float64), (grad)(y::Float64))
        std::Float64 = 3.0
        z = @trace(normal(x + y, std), :z)
        return z
    end
    @gen (grad) function fuzz2((grad)(x::Float64), (grad)(y::Float64))
        std::Float64 = 3.0
        z = @trace(normal(x + 2 * y, std), :z)
        q = @trace(bang2(z, y), :q)
        return z
    end
    sc2 = Switch(bang2, fuzz2)
    @gen (grad) function bam2(s::Int)
        x ~ sc2(s, 5.0, 3.0)
        return x
    end
    # ----.

    @testset "simulate" begin
        tr = simulate(bam2, (1, ))
        @test isapprox(get_score(tr), logpdf(normal, tr[:x => :z], 5.0 + 3.0, 3.0))
        tr = simulate(bam2, (2, ))
        @test isapprox(get_score(tr), logpdf(normal, tr[:x => :z], 5.0 + 2 * 3.0, 3.0) + logpdf(normal, tr[:x => :q => :z], tr[:x => :z] + 3.0, 3.0))
    end

    @testset "generate" begin
        chm = choicemap()
        chm[:x => :z] = 5.0
        tr, w = generate(bam2, (1, ), chm)
        assignment = get_choices(tr)
        @test assignment[:x => :z] == 5.0
        @test isapprox(w, logpdf(normal, 5.0, 5.0 + 3.0, 3.0))
        tr, w = generate(bam2, (2, ), chm)
        assignment = get_choices(tr)
        @test assignment[:x => :z] == 5.0
        @test isapprox(w, logpdf(normal, tr[:x => :z], 5.0 + 2 * 3.0, 3.0))
    end

    @testset "assess" begin
        chm = choicemap()
        chm[:x => :z] = 5.0
        w, ret = assess(bam2, (1, ), chm)
        @test isapprox(w, logpdf(normal, 5.0, 5.0 + 3.0, 3.0))
        chm[:x => :q => :z] = 5.0
        w, ret = assess(bam2, (2, ), chm)
        @test isapprox(w, logpdf(normal, 5.0, 5.0 + 2 * 3.0, 3.0) + logpdf(normal, 5.0, 5.0 + 3.0, 3.0))
    end

    @testset "propose" begin
        chm, w = propose(bam2, (1, ))
        @test isapprox(w, logpdf(normal, chm[:x => :z], 5.0 + 3.0, 3.0))
        chm, w = propose(bam2, (2, ))
        @test isapprox(w, logpdf(normal, chm[:x => :z], 5.0 + 2 * 3.0, 3.0) + logpdf(normal, chm[:x => :q => :z], chm[:x => :z] + 3.0, 3.0))
    end

    @testset "update" begin
        tr = simulate(bam2, (2, ))
        old_sc = get_score(tr)
        chm = choicemap((:x => :z, 5.0))
        new_tr, w, rd, discard = update(tr, (2, ), (UnknownChange(), ), chm)
        @test isapprox(old_sc, get_score(new_tr) - w)
        chm = choicemap((:x => :z, 10.0))
        new_tr, w, rd, discard = update(tr, (1, ), (UnknownChange(), ), chm)
        @test isapprox(old_sc, get_score(new_tr) - w)
    end

    @testset "regenerate" begin
        tr = simulate(bam2, (2, ))
        old_sc = get_score(tr)
        new_tr, w = regenerate(tr, (1, ), (UnknownChange(), ), select())
        @test isapprox(old_sc, get_score(new_tr) - w)
        new_tr, w = regenerate(tr, (2, ), (UnknownChange(), ), select())
        @test isapprox(old_sc, get_score(new_tr) - w)
        new_tr, w = regenerate(tr, (1, ), (UnknownChange(), ), select())
        @test isapprox(old_sc, get_score(new_tr) - w)
        new_tr, w = regenerate(tr, (1, ), (UnknownChange(), ), select(:x => :z))
        @test isapprox(old_sc, get_score(new_tr) - w)
    end

    @testset "choice gradients" begin
        for z in [1.0, 3.0, 5.0, 10.0]
            chm = choicemap((:x => :z, z))
            tr, _ = generate(bam2, (1, ), chm)
            sel = select(:x => :z)
            input_grads, choices, gradients = choice_gradients(tr, sel)
            expected_choice_grad = logpdf_grad(normal, z, 5.0 + 3.0, 3.0)
            @test isapprox(gradients[:x => :z], expected_choice_grad[1])
        end
    end

    # ------------ (More complex) hierarchy to test discard ------------ #

    # Model chunk.
    @gen (grad) function bang3((grad)(x::Float64), (grad)(y::Float64))
        std::Float64 = 3.0
        z = @trace(normal(x + y, std), :z)
        q = @trace(bang2(z, y), :q)
        return z
    end
    @gen (grad) function fuzz3((grad)(x::Float64), (grad)(y::Float64))
        std::Float64 = 3.0
        z = @trace(normal(x + 2 * y, std), :z)
        m = @trace(normal(x + 3 * y, std), :m)
        q = @trace(bang3(z, y), :q)
        return z
    end
    sc3 = Switch(bang3, fuzz3)
    @gen (grad) function bam3(s::Int)
        x ~ sc3(s, 5.0, 3.0)
        return x
    end
    # ----.

    @testset "update" begin
        tr = simulate(bam3, (2, ))
        old_sc = get_score(tr)
        chm = choicemap((:x => :z, 5.0))
        future_discarded = tr[:x => :z]
        new_tr, w, rd, discard = update(tr, (2, ), (UnknownChange(), ), chm)
        @test discard[:x => :z] == future_discarded
        @test isapprox(old_sc, get_score(new_tr) - w)
        chm = choicemap((:x => :z, 10.0))
        future_discarded = tr[:x => :q => :q => :z]
        new_tr, w, rd, discard = update(tr, (1, ), (UnknownChange(), ), chm)
        @test discard[:x => :q => :q => :z] == future_discarded
        @test isapprox(old_sc, get_score(new_tr) - w)
    end
end


text:
using Base64

## Text Utilities

function find_yaml_header(lines)
    yaml = []
    yaml_start, yaml_stop = 1, length(lines)

    # read out YAML block.
    if strip(lines[1]) == "---"
        for i = 2:length(lines)
            if strip(lines[i]) == "---"
                yaml = lines[1:i]
                yaml_stop = i
                break
            end
        end
    end
    return yaml, yaml_start, yaml_stop
end

function remove_yaml(file, key=nothing)
    lines = readlines(file, keep=true)
    yaml, yaml_start, yaml_stop = find_yaml_header(lines)
    if !isempty(yaml)
        open(file, "w+") do f
            # write back YAML lines.
            if key !== nothing
                for line in yaml
                    startswith(line, key * ":") || write(f, line)
                end
            end
            # write back non-YAML lines
            for line in lines[yaml_stop + 1 : end]
                write(f, line)
            end
        end
    end

    return yaml
end

function tidy_api(file)
	lines = readlines(file, keep=true)

    # Find the ID sections.
    for i = 1:(length(lines)-1)
		first = lines[i]
		second = lines[i+1]
        if startswith(first, "<a id=") && startswith(second, "**")
            first = replace(first, "\n" => "")
			second = replace(second, "\n" => "")
			# final_line = replace(first, ">" => ">$second", count = 1)

			lines[i] = "### $first $second"
			lines[i+1] = ""
        end
    end

    # Write lines back.
    open(file, "w+") do f
        for line in lines
            write(f, line)
        end
    end
end

## Utilities for Tutorial docs preprocess

"""
    fix_header_1

This function is used to add a first-level header into the
markdown file.

There are no H1 titles for documents in tutorials, so we need
this.
"""
function fix_header_1(file)
    lines = readlines(file, keep=true)
    yaml, yaml_start, yaml_stop = find_yaml_header(lines)
    isempty(yaml) && return

    for line in lines[yaml_stop + 1 : end]
        strip(line) == "" && continue
        if startswith(line, "# ") # has a h1 title
            return
        else # no h1 title
            break
        end
    end

    title = "Untitled"
    for line in yaml
        if startswith(line, "title:")
            title = line[7:end] |> strip
            break
        end
    end

    open(file, "w+") do f
        for line in yaml
            write(f, line)
        end
        write(f, "\n# $title\n")
        # write back non-YAML lines
        for line in lines[yaml_stop + 1 : end]
            write(f, line)
        end
    end
end

function fix_image_path(file)
    # ![](/tutorials/figures/image.png)
    # to
    # ![](../figures/image.png)
    lines = readlines(file, keep=true)
    open(file, "w+") do f
        for line in lines
            line = replace(line, "](/tutorials/" => "](../")
            write(f, line)
        end
    end
end

## Building Utilities

function preprocess_markdown(folder)
    yaml_dict = Dict()

    try
        for (root, dirs, files) in walkdir(folder)
            for file in files
				if endswith(file, ".md")
	                full_path = joinpath(root, file)
	                yaml_dict[full_path] = remove_yaml(full_path)
				end
            end
        end
    catch e
        # println("Markdown copy error: $e")
        rethrow(e)
    end

    return yaml_dict
end

function postprocess_markdown(folder, yaml_dict; original = "")
    try
        for (root, dirs, files) in walkdir(folder)
            for file in files
                full_path = joinpath(root, file)
                original_path = full_path

                if length(original) > 0
                    # original_path = abspath(original, file)
                    original_path = replace(full_path, folder => original)
                    # original_path = replace(
                    #     full_path,
                    #     joinpath("docs", "site") => joinpath("docs", "src")
                    # )
                end

                if haskey(yaml_dict, original_path)
                    # println("Original: $original_path => Full path: $full_path")

                    txt = open(f -> read(f, String), full_path)
                    open(full_path, "w+") do f
                        # Add in the yaml block.
                        for line in yaml_dict[original_path]
                            write(f, line)
                        end

                        # This is needed to replace the hyperlinks Documenter.jl generates
                        # for the API pages.
						txt = replace(txt, "api.md" => "{{site.baseurl}}/docs/library/")
						txt = replace(txt, "bijectors.md" => "{{site.baseurl}}/docs/library/bijectors/")
						txt = replace(txt, "advancedhmc.md" => "{{site.baseurl}}/docs/library/advancedhmc/")

                        # Add the rest of the text.
						if original_path == full_path
							write(f, txt)
						else
							write(f, replace(txt, "![](figures" => "![](/{{site.baseurl}}/tutorials/figures"))
						end
                    end
                elseif endswith(file, ".md")
                    println("Original: $original_path")
                    println("Full:     $full_path \n")
                end

				# Make specific api items headers.
				if file == "api.md" && original_path != full_path
					tidy_api(full_path)
				end
            end
        end
    catch e
        # println("Markdown copy error: $e")
        rethrow(e)
    end

    return yaml_dict
end

function with_clean_docs(func, source, target)
    src_temp = mktempdir()
    cp(source, src_temp, force=true)
    yaml_dict = preprocess_markdown(source)

    try
        func(source, target)
    catch e
        rethrow(e)
    finally
        # Put back the original files in the event of an error.
        cp(src_temp, source_path, force=true)
        rm(src_temp, recursive=true)
    end
    postprocess_markdown(build_path, yaml_dict, original=source_path)
end

function copy_tutorial(tutorial_path)
    isdir(tutorial_path) || mkpath(tutorial_path)
    # Clone TuringTurorials
    tmp_path = tempname()
    mkdir(tmp_path)
    clone("https://github.com/TuringLang/TuringTutorials", tmp_path)

    # Move to markdown folder.
    md_path = joinpath(tmp_path, "markdown")

    # Copy the .md versions of all examples.
    try
        @debug(md_path)
        for file in readdir(md_path)
            full_path = joinpath(md_path, file)
            target_path = joinpath(tutorial_path, file)
            println("Copying $full_path to $target_path")
            cp(full_path, target_path, force=true)
            if endswith(target_path, ".md")
                # remove_yaml(target_path, "permalink")
                fix_header_1(target_path)
                fix_image_path(target_path)
            end
        end
        index = joinpath(@__DIR__, "src/tutorials/index.md")
        cp(index, tutorial_path * "/index.md", force=true)
    catch e
        rethrow(e)
    finally
        rm(tmp_path, recursive=true)
    end
end

function with_baseurl(func, baseurl)
    jekyll_config = joinpath(@__DIR__, "_config.yml")
    lines = readlines(jekyll_config, keep=true)
    open(jekyll_config, "w+") do f
        for line in lines
            if startswith(line, "baseurl:")
                write(f, "baseurl: \"$(baseurl)\"\n")
                continue
            end
            write(f, line)
        end
    end
    ret = func()
    open(jekyll_config, "w+") do f
        for line in lines
            write(f, line)
        end
    end
    return ret
end


text:
using MPI

using Printf

function do_reduce()
    comm = MPI.COMM_WORLD

    MPI.Barrier(comm)

    root = 0
    r = [1,2,3]
    op = (x, y) -> broadcast(*, x, 2y)
    sr = MPI.Reduce(r, op, root, comm)
    if MPI.Comm_rank(comm) == root
        @printf("sum of ranks: %s\n", sr)
    end
end

function main()
    MPI.Init()
    do_reduce()
    MPI.Finalize()
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# Windows COM Structures constants and types

struct GUID
    d1::Cuint
    d2::Cushort
    d3::Cushort
    d4::NTuple{8, Cuchar}
end

const IID_ITypeInfo  = GUID(0x00020401, 0x0000, 0x0000, (0xc0,0x00, 0x00,0x00,0x00,0x00,0x00,0x46))
const IID_IDispatch  = GUID(0x00020400, 0x0000, 0x0000, (0xc0,0x00,0x00,0x00,0x00,0x00,0x00,0x46))
const IID_NULL       = GUID(0, 0, 0, (0, 0, 0, 0, 0, 0, 0, 0))

const CLSCTX_ALL = 0x1 | 0x2 | 0x4 | 0x10
const CLSCTX_SERVER = 0x1 | 0x4 | 0x10
const LCID = ccall((:GetUserDefaultLCID, "kernel32.dll"), Cint, ())

const PARAMFLAG_OUT  = UInt16(0x2)
const VT_USERDEFINED = UInt16(29)
const VT_VARIANT     = UInt16(12)
const VT_ARRAY       = UInt16(0x2000)
const VT_BYREF       = UInt16(0x4000)
const VT_TYPEMASK    = UInt16(0xfff)

abstract type IUnknown end
abstract type IDispatch end
abstract type IDispaxxx end
abstract type ITypeInfo end

struct Variant
    vt::Cshort
    r1::Cshort
    r2::Cshort
    r3::Cshort
    data::UInt128
end

struct SafeArrayBound
    cElements::Culong
    lLbound::Clong
end

struct SafeArray
    cDims::Cushort
    fFeatures::Cushort
    cbElements::Culong
    cLocks::Culong
    pvData::Ptr{Cvoid}
    bounds::SafeArrayBound
end

struct DISPPARAMS
    rgvarg::Ptr{Variant}
    rgdispidNamedArgs::Ptr{Int32}
    cArgs::Cuint
    cNamedArgs::Cuint
end

struct TYPEDESC
    union::NTuple{14, Cuchar}
    vt::UInt16
end

struct IDLDESC
    dwReserved::Ptr{Cvoid}
    wIDLFlags::UInt16
end

struct ELEMDESC
    tdesc::TYPEDESC
    paramdesc::UInt128 # 10 padded to 16
end

struct TYPEATTR
    guid::GUID
    lcid::UInt32
    dwReserved::UInt32
    memidConstructor::UInt32
    memidDestructor::UInt32
    lpstrSchema::Ptr{Cvoid}
    cbSizeInstance::UInt32
    typekind::UInt32
    cFuncs::UInt32
    cVars::UInt32
    cImplTypes::UInt32
    cbSizeVft::UInt32
    cbAlignment::UInt32
    wTypeFlags::UInt32
    wMajorVerNum::UInt32
    wMinorVerNum::UInt32
    tdescAlias::TYPEDESC
    idldescType::IDLDESC
end

struct EXCEPINFO
    wCode::Cushort
    wReserved::Cushort
    bstrSource::Ptr{UInt16}
    bstrDescription::Ptr{UInt16}
    bstrHelpFile::Ptr{UInt16}
    dwHelpContext::Cushort
    pvReserved::Ptr{Cvoid}
    pfnDeferredFillIn::Ptr{Cvoid}
    scode::Clong
end

struct FUNCDESC
    memid::Clong
    lprgscode::Ptr{Cvoid}
    lprgelemdescParam::Ptr{ELEMDESC}
    funckind::UInt32
    invkind::UInt32
    callconv::UInt32
    cParams::UInt16
    cParamsOpt::UInt16
    oVft::UInt16
    cScodes::UInt16
    elemdescFunc::ELEMDESC
    wFuncFlags::UInt16
end

mutable struct FuncData
    memid::Clong
    invoke_mask::UInt32
    param_flags::Array{UInt16}
end

mutable struct TypeData
    info::Ptr{ITypeInfo}
    func::Dict{Symbol, FuncData}
end

struct ComObject
    ptr::Ptr{IDispatch}
    dat::TypeData
end

# Global Dictionaries

type_cache = Dict{UInt64, TypeData}()

const typeindex = Dict(
    Ptr{Cvoid}      => [26 1 24],
    Ptr{IDispatch}  => [9],
    Ptr{IUnknown}   => [13],
    Ptr{Int32}      => [37],
    Ptr{UInt32}     => [38],
    Ptr{UInt16}     => [8 30 31],
    Variant         => [12],
    Array           => [27],
    Int8            => [16],
    Int16           => [2],
    Int32           => [3 22],
    Int64           => [20],
    UInt8           => [17],
    UInt16          => [18 29],
    UInt32          => [19 23 25],
    UInt64          => [21],
    Float32         => [4],
    Float64         => [5 7],
    GUID            => [72],
    Bool            => [11]
)

function reversedict(dict::Dict)
    rev = Dict()
    for (k, v) in dict
        for n in v
            rev[n] = k
        end
    end
    return rev
end

const indextype = reversedict(typeindex)

function gettype(idx::Cshort)
    if (idx & VT_ARRAY) > 0
        return Array
    else
        return get(indextype, idx, Nothing)
    end
end

function getindex(::Type{T}) where {T}
    if T == Int64
        # Int64 not suppported in CATIA.
        return typeindex[Int32][1]
    else
        return typeindex[T][1]
    end
end

max_stars_repo_name: zaika-generative/ShapeFactory

max_issues_repo_name: zaika-generative/ShapeFactory

max_forks_repo_name: zaika-generative/ShapeFactory