text:
function is_Displaz_installed()
    try
        @static if VERSION < v"0.7.0"
            return Pkg.installed("Displaz") != nothing
        else
            return haskey(Pkg.installed(), "Displaz")
        end
    catch
        return false
    end
end

@static if is_Displaz_installed()
    using Displaz

    # Plot Lines
    function Displaz.plot3d(l::Line, args...; label = "Line", kwargs...)
        vertices = [l[0], l[end]]
        Displaz.plot3d(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    function Displaz.plot3d(ls::AbstractVector{<:Line}, args...; label = "Lines [$(length(ls))]", kwargs...)
        vertices = zeros(SVector{3,Float64}, 2*length(ls))
        for i = 1:length(ls)
            vertices[2*i-1] = ls[i][0.0]
            vertices[2*i] = ls[i][end]
        end
        Displaz.plot3d(vertices, args...; label = label, markershape = '-', linebreak = 2, kwargs...)
    end

    function Displaz.plot3d!(l::Line, args...; label = "Line", kwargs...)
        vertices = [l[0], l[end]]
        Displaz.plot3d!(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    function Displaz.plot3d!(ls::AbstractVector{<:Line}, args...; label = "Lines [$(length(ls))]", kwargs...)
        vertices = zeros(SVector{3,Float64}, 2*length(ls))
        for i = 1:length(ls)
            vertices[2*i-1] = ls[i][0.0]
            vertices[2*i] = ls[i][end]
        end
        Displaz.plot3d!(vertices, args...; label = label, markershape = '-', linebreak = 2, kwargs...)
    end

    # Plot LineStrings
    function Displaz.plot3d(ls::LineString, args...; label = "LineString", kwargs...)
        Displaz.plot3d(lines(ls), args...; label = label, kwargs...)
    end

    function Displaz.plot3d(lss::AbstractVector{<:LineString{N,T}}, args...; label = "LineStrings [$(length(lss))]", kwargs...) where {N,T}
        tmp = Vector{Line{N,T}}()
        for ls in lss
            append!(tmp, lines(ls))
        end
        Displaz.plot3d(tmp, args...; label = label, kwargs...)
    end

    function Displaz.plot3d!(ls::LineString, args...; label = "LineString", kwargs...)
        Displaz.plot3d!(lines(ls), args...; label = label, kwargs...)
    end

    function Displaz.plot3d!(lss::AbstractVector{<:LineString{N,T}}, args...; label = "LineStrings [$(length(lss))]", kwargs...) where {N,T}
        tmp = Vector{Line{N,T}}()
        for ls in lss
            append!(tmp, lines(ls))
        end
        Displaz.plot3d!(tmp, args...; label = label, kwargs...)
    end

    # Plot Polygon
    function Displaz.plot3d(p::Polygon, args...; label = "Polygon", kwargs...)
        Displaz.plot3d(lines(p), args...; label = label, kwargs...)
    end

    function Displaz.plot3d!(p::Polygon, args...; label = "Polygon", kwargs...)
        Displaz.plot3d!(lines(p), args...; label = label, kwargs...)
    end

    function Displaz.plot3d(ps::AbstractVector{<:Polygon{N,T}}, args...; label = "Polygons [$(length(ps))]", kwargs...) where {N,T}
        tmp = Vector{Line{N,T}}()
        for p in ps
            append!(tmp, lines(p))
        end
        Displaz.plot3d(tmp, args...; label = label, kwargs...)
    end

    function Displaz.plot3d!(p::AbstractVector{<:Polygon}, args...; label = "Polygon", kwargs...)
        Displaz.plot3d!(lines(p), args...; label = label, kwargs...)
    end

    function Displaz.plot3d!(ps::AbstractVector{<:Polygon{N,T}}, args...; label = "Polygons [$(length(ps))]", kwargs...) where {N,T}
        tmp = Vector{Line{N,T}}()
        for p in ps
            append!(tmp, lines(p))
        end
        Displaz.plot3d!(tmp, args...; label = label, kwargs...)
    end

    # Plot Quadratics
    function Displaz.plot3d(q::Quadratic, args...; label = "Quadratic", kwargs...)
        vertices = q[linspace(zero(length(q)), length(q), 20)]
        Displaz.plot3d(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    function Displaz.plot3d(qs::AbstractVector{<:Quadratic}, args...; label = "Quadratics [$(length(qs))]", kwargs...)
        vertices = zeros(SVector{3,Float64}, 20*length(qs))
        for i = 1:length(qs)
            vertices[(1+20*(i-1)):(20*i)] = qs[i][linspace(zero(length(qs[i])), length(qs[i]), 20)]
        end
        Displaz.plot3d(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    function Displaz.plot3d!(q::Quadratic, args...; label = "Quadratic", kwargs...)
        vertices = q[linspace(zero(length(q)), length(q), 20)]
        Displaz.plot3d!(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    function Displaz.plot3d!(qs::AbstractVector{<:Quadratic}, args...; label = "Quadratics [$(length(qs))]", kwargs...)
        vertices = zeros(SVector{3,Float64}, 20*length(qs))
        for i = 1:length(qs)
            vertices[(1+20*(i-1)):(20*i)] = qs[i][linspace(zero(length(qs[i])), length(qs[i]), 20)]
        end
        Displaz.plot3d!(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    # Plot Catenaries
    function Displaz.plot3d(c::Catenary, args...; label = "Catenary", kwargs...)
        vertices = c[linspace(c.lmin, c.lmax, 20)]
        Displaz.plot3d(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    function Displaz.plot3d(cs::AbstractVector{<:Catenary}, args...; label = "Catenaries [$(length(cs))]", kwargs...)
        vertices = zeros(SVector{3,Float64}, 20*length(cs))
        for i = 1:length(cs)
            vertices[(1+20*(i-1)):(20*i)] = cs[i][linspace(cs[i].lmin, cs[i].lmax, 20)]
        end
        Displaz.plot3d(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    function Displaz.plot3d!(c::Catenary, args...; label = "Catenary", kwargs...)
        vertices = c[linspace(c.lmin, c.lmax, 20)]
        Displaz.plot3d!(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    function Displaz.plot3d!(cs::AbstractVector{<:Catenary}, args...; label = "Catenaries [$(length(cs))]", kwargs...)
        vertices = zeros(SVector{3,Float64}, 20*length(cs))
        for i = 1:length(cs)
            vertices[(1+20*(i-1)):(20*i)] = cs[i][linspace(cs[i].lmin, cs[i].lmax, 20)]
        end
        Displaz.plot3d!(vertices, args...; label = label, markershape = '-', linebreak = 20, kwargs...)
    end

    # Plot BoundingBoxes
    function Displaz.plot3d(bb::BoundingBox, args...; label = "BoundingBox", kwargs...)
        lines = wireframe(bb)
        Displaz.plot3d(lines, args...; label = label, kwargs...)
    end

    function Displaz.plot3d(bbs::AbstractVector{<:BoundingBox}, args...; label = "BoundingBoxes [$(length(bb))]", kwargs...)
        lines = Vector{Line{eltype(eltype(bbs))}}()
        for i = 1:length(bbs)
            append!(lines, wireframe(bbs[i]))
        end
        Displaz.plot3d(lines, args...; label = label, kwargs...)
    end

    function Displaz.plot3d!(bb::BoundingBox, args...; label = "BoundingBox", kwargs...)
        lines = wireframe(bb)
        Displaz.plot3d!(lines, args...; label = label, kwargs...)
    end

    function Displaz.plot3d!(bbs::AbstractVector{<:BoundingBox}, args...; label = "BoundingBoxes [$(length(bb))]", kwargs...)
        lines = Vector{Line{eltype(eltype(bbs))}}()
        for i = 1:length(bbs)
            append!(lines, wireframe(bbs[i]))
        end
        Displaz.plot3d!(lines, args...; label = label, kwargs...)
    end
end


text:
using Test
using Documenter
using ExperimentalDesign

tests = ["variance_predictions.jl"]

@testset "ExperimentalDesign" begin
    for test in tests
        include(test)
    end

    @testset "Doctests" begin
        DocMeta.setdocmeta!(ExperimentalDesign,
                            :DocTestSetup,
                            :(using ExperimentalDesign, Distributions,
                              Random, StatsModels, DataFrames;
                              Random.seed!(443591););
                            recursive = true)
        doctest(ExperimentalDesign)
    end
end

ext: jl


text:
using CSV
using Plots
using StatPlots
using DataFrames
using GLM

gr()
cd(raw"D:\workspace\github\blog-py\blog\static\essay_resources\Notes_on_Computer_Age_Statistical_Inference")

kidney = CSV.read("kidney.csv", nullable=false)
@df kidney scatter(:age, :Tot)

X = hcat(ones(nrow(kidney)), kidney[:age])
y = kidney[:Tot]
OLS = fit(LinearModel, X, y)
# GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,Base.LinAlg.Cholesky{Float64,Array{Float64,2}}}}:

# Coefficients:
#        Estimate Std.Error  t value Pr(>|t|)
# x1      2.86067  0.359561  7.95603   <1e-12
# x2   -0.0786009 0.0090557 -8.67972   <1e-14

age_samples = collect(20:10:90)
Xtest = hcat(ones(length(age_samples)), age_samples)
pred = predict(OLS, Xtest, :confint)

for i in 1:size(pred, 1)
    y_pred, y_lower, y_upper = pred[i, :]
    display(plot!([age_samples[i],age_samples[i]], [ y_lower, y_upper], linewidth = 3))
end

plot!(age_samples[[1, end]], pred[[1, end], 1], legend=:none, linewidth=3)
savefig("Figure_1_1.png")


text:
# Main data structure for quadtree-like point locator
struct PointLocator{NDIMS}
  center::SVector{NDIMS, Float64}
  length::Float64
  max_point_ids::Int
  point_ids::Vector{Int}
  children::Vector{PointLocator{NDIMS}}
end

function PointLocator{NDIMS}(center::AbstractVector, length::Float64, offset=0.0) where NDIMS
  # Use offset to avoid ambiguitites for points falling on coordinate lines
  center_ = SVector{NDIMS}(center) .+ offset
  length_ = length + 2 * offset

  # Use at most 4 points per locator node
  max_point_ids = 4
  point_ids = Vector{Int}()
  children = Vector{PointLocator{NDIMS}}()

  return PointLocator{NDIMS}(center_, length_, max_point_ids, point_ids, children)
end

# Return if locator has no children
is_leaf(pl::PointLocator) = isempty(pl.children)

# Insert point into locator
function insert!(pl::PointLocator, points::Vector{P}, point::P) where {P}
  # Check if locator is leaf
  if is_leaf(pl)
    # If locator is leaf, check if point already exists
    point_id = get_point_id(pl, points, point)

    # Return point_id if found
    if point_id > 0
      return point_id
    end

    # Otherwise check if there is enough room in current locator or if it needs to be refined
    if length(pl.point_ids) < pl.max_point_ids
      # Add point to this locator and to list of points
      point_id = length(points) + 1
      push!(points, point)
      push!(pl.point_ids, point_id)
      return point_id
    else
      # Refine locator
      refine!(pl)

      # Re-add to newly refined locator
      insert!(pl, points, point)
    end
  else
    # If locator is non-lef, find appropriate child locator and add point there
    child_locator = pl.children[get_child_id(pl, point)]
    insert!(child_locator, points, point)
  end
end
function insert!(pl::PointLocator, points::Vector{P}, coordinates) where {P}
  insert!(pl, points, P(coordinates))
end

# Refine point locator and move points to child locators
function refine!(pl::PointLocator{2})
  # Store for convenience
  dx = pl.length / 2
  x = pl.center[1]
  y = pl.center[2]

  # Add lower left child
  push!(pl.children, PointLocator{2}(SVector{2}((x - dx/2, y - dx/2)), dx))

  # Add lower right child
  push!(pl.children, PointLocator{2}(SVector{2}((x + dx/2, y - dx/2)), dx))

  # Add upper left child
  push!(pl.children, PointLocator{2}(SVector{2}((x - dx/2, y + dx/2)), dx))

  # Add upper right child
  push!(pl.children, PointLocator{2}(SVector{2}((x + dx/2, y + dx/2)), dx))
end
function refine!(pl::PointLocator{3})
  # Store for convenience
  dx = pl.length / 2
  x = pl.center[1]
  y = pl.center[2]
  z = pl.center[3]

  # Add bottom lower left child
  push!(pl.children, PointLocator{3}(SVector{3}((x - dx/2, y - dx/2, z - dx/2)), dx))

  # Add bottom lower right child
  push!(pl.children, PointLocator{3}(SVector{3}((x + dx/2, y - dx/2, z - dx/2)), dx))

  # Add bottom upper left child
  push!(pl.children, PointLocator{3}(SVector{3}((x - dx/2, y + dx/2, z - dx/2)), dx))

  # Add bottom upper right child
  push!(pl.children, PointLocator{3}(SVector{3}((x + dx/2, y + dx/2, z - dx/2)), dx))

  # Add top lower left child
  push!(pl.children, PointLocator{3}(SVector{3}((x - dx/2, y - dx/2, z + dx/2)), dx))

  # Add top lower right child
  push!(pl.children, PointLocator{3}(SVector{3}((x + dx/2, y - dx/2, z + dx/2)), dx))

  # Add top upper left child
  push!(pl.children, PointLocator{3}(SVector{3}((x - dx/2, y + dx/2, z + dx/2)), dx))

  # Add top upper right child
  push!(pl.children, PointLocator{3}(SVector{3}((x + dx/2, y + dx/2, z + dx/2)), dx))
end

# Get id of child locator for given point
function get_child_id(pl::PointLocator{2}, point)
  if point[2] < pl.center[2]
    if point[1] < pl.center[1]
      # Lower left child
      return 1
    else
      # Lower right child
      return 2
    end
  else
    if point[1] < pl.center[1]
      # Upper left child
      return 3
    else
      # Upper right child
      return 4
    end
  end
end
function get_child_id(pl::PointLocator{3}, point)
  if point[3] < pl.center[3]
    if point[2] < pl.center[2]
      if point[1] < pl.center[1]
        # Bottom lower left child
        return 1
      else
        # Bottom lower right child
        return 2
      end
    else
      if point[1] < pl.center[1]
        # Bottom upper left child
        return 3
      else
        # Bottom upper right child
        return 4
      end
    end
  else
    if point[2] < pl.center[2]
      if point[1] < pl.center[1]
        # Top lower left child
        return 5
      else
        # Top lower right child
        return 6
      end
    else
      if point[1] < pl.center[1]
        # Top upper left child
        return 7
      else
        # Top upper right child
        return 8
      end
    end
  end
end

# Get point id if point exists or zero otherwise
function get_point_id(pl::PointLocator, points::Vector{P}, point::P) where {P}
  # Iterate over point ids, extract point coordinates and compare to point
  for point_id in pl.point_ids
    if isapprox(point, points[point_id], atol=1e-13)
      return point_id
    end
  end

  # If no point was found, return zero
  return 0
end


text:
module logging
using Dates: Dates
using Logging
using InteractiveUtils

import Term: Theme, theme, textlen, square_to_round_brackets, escape_brackets, reshape_text
import ..markup: has_markup
import ..box: ROUNDED
import ..style: apply_style
import ..renderables: AbstractRenderable

export TermLogger, install_term_logger

tprint = (println ∘ apply_style)

# ---------------------------- create term logger ---------------------------- #
"""
    TermLogger

Custom logger type.
"""
struct TermLogger <: Logging.AbstractLogger
    io::IO
    theme::Theme
end
TermLogger(theme::Theme) = TermLogger(stderr, theme)

# set logger beavior
Logging.min_enabled_level(logger::TermLogger) = Logging.Info
function Logging.shouldlog(logger::TermLogger, level, _module, group, id)
    return true
end
Logging.catch_exceptions(logger::TermLogger) = true

# --------------------------- handle logger message -------------------------- #

"""
    print_closing_line(color::String, width::Int)

Print the final line of a log message with style and date info
"""
function print_closing_line(color::String, width::Int = 48)
    tprint(
        "  [$color bold dim]$(ROUNDED.bottom.left)" *
        "$(ROUNDED.row.mid)"^(width) *
        "[/$color bold dim]",
    )
    _date = Dates.format(Dates.now(), "e, dd u yyyy")
    _time = Dates.format(Dates.now(), "HH:MM:SS")
    pad = width - textlen(_date * _time) - 2
    return tprint(
        " "^pad * "[dim]$_date[/dim] [bold dim underline]$_time[/bold dim underline]"
    )
end

"""
    Logging.handle_message(logger::TermLogger,

Handle printing of log messages, with style!.

In addition to the log message and info such as file/line and time of log,
it prints kwargs styled by their type.
"""
function Logging.handle_message(
    logger::TermLogger, lvl, msg, _mod, group, id, file, line; kwargs...
)

    # get name of function where logging message was called from
    fname = ""
    for frame in stacktrace()
        if "$(frame.file)" == file && frame.line == line
            fname = ".[$(logger.theme.func) underline]$(frame.func)[/$(logger.theme.func) underline]"
        end
    end

    # prepare styles
    if lvl == Logging.Info
        color = logger.theme.info
    elseif lvl == Logging.Debug
        color = logger.theme.debug
    elseif lvl == Logging.Warn
        color = logger.theme.warn
    elseif lvl == Logging.Error
        color = logger.theme.error
    else
        color = "#90CAF9"
    end

    outline_markup = "$color dim"
    hor = "[$outline_markup]▶[/$outline_markup]"
    vert = "[$outline_markup]" * ROUNDED.mid.left * "[/$outline_markup]"

    # style message
    if msg isa AbstractString
        msg = has_markup(msg) ? msg : "[#8abeff]$msg[/#8abeff]"
        msg = reshape_text(msg, 64)
    else
        msg = string(msg)
    end

    # get the first line of information
    content = "[$color underline bold]@$(string(lvl))[/$color underline bold] [#edb46f ]($(_mod)$fname):[/#edb46f ]"

    # for multi-lines message, print each line separately.
    msg_lines::Vector{AbstractString} = split(msg, "\n")
    for n in 2:length(msg_lines)
        msg_lines[n] = "  $vert   " * " "^textlen(content) * "[#8abeff]" * msg_lines[n]
    end
    content *= "  " * msg_lines[1]
    tprint(content)
    tprint.(msg_lines[2:end])

    # if no kwargs we can just quit
    if length(kwargs) == 0
        print_closing_line(color)
        return nothing
    end

    # get padding width
    _types = string.(typeof.([v for v in values(kwargs)]))
    for (i, _type) in enumerate(values(kwargs))
        if typeof(_type) <: Function
            _types[i] = string(Function)
        end
    end

    wpad = max(textlen.(_types)...) + 2
    namepad = max(textlen.(string.([v for v in keys(kwargs)]))...)

    # print all kwargs
    tprint("  $vert")
    for ((k, v), _type) in zip(kwargs, _types)
        # get line stub
        pad = wpad - textlen(_type)
        line =
            "  $vert [$(theme.type) dim]($(_type))[/$(theme.type) dim]" *
            " "^pad *
            hor *
            "  [bold #ff9959]$k[/bold #ff9959]"

        epad = namepad - textlen(string(k))
        line *= " "^epad * " [bold red]=[/bold red] "
        lpad = textlen(line) - 4

        # get value style
        if v isa Number
            _style = logger.theme.number
        elseif v isa Symbol
            _style = logger.theme.symbol
        elseif v isa AbstractString
            _style = logger.theme.string
        elseif v isa AbstractVector
            _style = logger.theme.number
            _size = size(v)
            v = square_to_round_brackets(string(v))
            v = textlen(v) > 33 ? v[1:30] * " ...)" : v
            v *= "\n [dim]size: $(_size)[/dim]"
        elseif v isa AbstractArray || v isa AbstractMatrix
            _style = logger.theme.number
            v =
                "$(typeof(v)) [dim]<: $(supertypes(typeof(v))[end-1])[/dim]" *
                "\n [dim]size: $(size(v))[/dim]"
        elseif v isa Function
            _style = logger.theme.func
        elseif v isa AbstractRenderable
            _style = "default"
            v = "$(typeof(v)) \e[2m(size: $(v.measure))\e[0m"
        else
            _style = nothing
        end

        # print value lines
        vlines = split(string(v), "\n")

        if !isnothing(_style)
            vlines = ["[$_style]$ln[/$_style]" for ln in vlines]
        end

        if length(vlines) == 1
            tprint(line * vlines[1])
        else
            tprint(line * vlines[1])
            for ln in vlines[2:end]
                tprint("  $vert " * " "^lpad * ln)
            end
        end
    end
    return print_closing_line(color)
end

# ---------------------------- install term logger --------------------------- #

"""
    install_term_logger(theme::Theme=theme)

Install `TermLogger` as the global logging system.

`theme::Theme` can be passed to specify the theme to use for styling objects.
"""
function install_term_logger(theme::Theme = theme)
    _logger = TermLogger(theme)
    return global_logger(_logger)
end

end


text:
@testset "Genus" begin
  Qx, x = QQ["x"]
  K, a = NumberField(x^2 - 2, "a")
  OK = maximal_order(K)
  Kt, t  = K["t"]

  E1, b1 = NumberField(t^2 - a, "b1") # ramified at 2
  E2, b2 = NumberField(t^2 - 5, "b2") # unramified at 2

  p = prime_decomposition(OK, 2)[1][1]
  q = prime_decomposition(OK, 3)[1][1]

  #
  # ramified & dyadic
  #

  g = @inferred genus(HermLat, E1, p, [(0, 1, 1, 0), (2, 2, -1, 1)], type = :det)
  @test isramified(g)
  @test !issplit(g)
  @test !isinert(g)
  @test isdyadic(g)
  @test scales(g) == [0, 2]
  @test ranks(g) == [1, 2]
  @test dets(g) == [1, -1]
  @test norms(g) == [0, 1]
  @test det(g) == -1
  @test discriminant(g, 1) == 1
  @test discriminant(g, 2) == 1
  @test discriminant(g) == 1
  @test g == g
  @test g != genus(HermLat, E1, p, [(0, 1, 1, 0), (2, 2, -1, 1), (6, 1, -1, 3)])
  @test g != genus(HermLat, E1, p)
  @test g != genus(HermLat, E1, p, [(0, 1, 1, 0), (2, 3, -1, 1)], type = :det)
  @test g != genus(HermLat, E1, p, [(0, 1, 1, 0), (1, 2, -1, 1)], type = :det)
  @test genus(HermLat, E1, p) != genus(HermLat, E1, q)
  @test genus(HermLat, E1, p) != genus(HermLat, E2, p)

  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E1; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E1, K(det(gram_matrix(g))) * inv(d), p)

  # empty genus
  g = @inferred genus(HermLat, E1, p)
  @test isramified(g)
  @test !issplit(g)
  @test !isinert(g)
  @test isdyadic(g)
  @test rank(g) == 0
  @test scales(g) == []
  @test ranks(g) == []
  @test dets(g) == []
  @test norms(g) == []
  @test det(g) == 1
  @test discriminant(g) == 1
  @test g == g

  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E1; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E1, K(det(gram_matrix(g))) * inv(d), p)

  g = @inferred genus(HermLat, E1, p, [(0, 1, 1, 0), (2, 2, -1, 1)], type = :disc)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  @test g == g
  @test representative(g) in g
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E1; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E1, K(det(gram_matrix(g))) * inv(d), p)

  # negative scale
  g = @inferred genus(HermLat, E1, p, [(-2, 1, 1, -1), (2, 2, -1, 1)], type = :disc)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  @test g == g
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E1; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E1, K(det(gram_matrix(g))) * inv(d), p)

  @test_throws ArgumentError genus(HermLat, E1, p, [(0, 1, 1, 1), (2, 2, -1, 0)], type = :det)
  @test_throws ArgumentError genus(HermLat, E1, p, [(0, 1, 1, 1), (2, 2, -1, 0)], type = :disc)
  @test_throws ArgumentError genus(HermLat, E1, p, [(0, 1, 1, 0), (2, 2, -1, 1)], type = :bla)
  # scale not increasing
  @test_throws ArgumentError genus(HermLat, E1, p, [(3, 1, 1, 0), (2, 2, -1, 1)])
  # negative rank
  @test_throws ArgumentError genus(HermLat, E1, p, [(1, -1, 1, 0), (2, 2, -1, 1)])
  # bad determinant class
  @test_throws ArgumentError genus(HermLat, E1, p, [(0, 1, 2, 0), (2, 2, -1, 1)])
  # wrong norm valuation (must 1 * 2 since rank is odd)
  @test_throws ArgumentError genus(HermLat, E1, p, [(1, 1, 1, 3)], type = :det)
  # wrong prime
  @test_throws ArgumentError genus(HermLat, E1, p, [(1, 1, 1)], type = :det)
  # wrong prime
  @test_throws ArgumentError genus(HermLat, E1, p, [(1, 1)])

  #
  # unramified & inert & dyadic
  #

  g = genus(HermLat, E2, p, [(0, 1, 1), (3, 5, -1)], type = :det)

  @test g == genus(HermLat, E2, p, [(0, 1), (3, 5)])

  @test !isramified(g)
  @test isinert(g)
  @test !issplit(g)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  @test g == g
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E2; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E2, K(det(gram_matrix(g))) * inv(d), p)

  # rank zero genus
  g = @inferred genus(HermLat, E2, p)
  @test rank(g) == 0
  @test !isramified(g)
  @test isinert(g)
  @test !issplit(g)
  @test scales(g) == []
  @test ranks(g) == []
  @test dets(g) == []
  @test det(g) == 1
  @test discriminant(g) == 1
  @test g == g

  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E2; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E2, K(det(gram_matrix(g))) * inv(d), p)

  g = genus(HermLat, E2, p, [(0, 1, 1), (2, 2, +1)], type = :disc)

  @test g == genus(HermLat, E2, p, [(0, 1), (2, 2)])

  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  @test g == g
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E2; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E2, K(det(gram_matrix(g))) * inv(d), p)

  # wrong det class
  @test_throws ArgumentError genus(HermLat, E2, p, [(0, 1, -1), (2, 2, +1)], type = :det)
  # wrong det class
  @test_throws ArgumentError genus(HermLat, E2, p, [(0, 1, +1), (2, 2, -1)], type = :det)
  # wrong det class
  @test_throws ArgumentError genus(HermLat, E2, p, [(0, 1, +2), (2, 2, -1)], type = :det)
  # wrong type
  @test_throws ArgumentError genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :bla)
  # negative rank
  @test_throws ArgumentError genus(HermLat, E2, p, [(0, -1, 1), (2, 2, +1)], type = :det)
  # non-increasing scale
  @test_throws ArgumentError genus(HermLat, E2, p, [(3, 1, 1), (2, 2, +1)], type = :det)
  # wrong prime
  @test_throws ArgumentError genus(HermLat, E2, p, [(0, 1, 1, 1)], type = :det)

  #
  # ramified & non-dyadic
  #

  p = prime_decomposition(OK, 5)[1][1]
  g = genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :det)
  @test rank(g) == 3
  @test isramified(g)
  @test !issplit(g)
  @test !isinert(g)
  @test g == g

  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E2; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E2, K(det(gram_matrix(g))) * inv(d), p)

  # Test some problem with non_norm_rep
  # It is cached on the genus
  @test all(1:100) do i
    g = genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :det)
    hermitian_lattice(E2; gram = gram_matrix(g)) in g
  end

  # rank zero genus
  g = @inferred genus(HermLat, E2, p)
  @test rank(g) == 0
  @test isramified(g)
  @test !issplit(g)
  @test !isinert(g)
  @test scales(g) == []
  @test ranks(g) == []
  @test dets(g) == []
  @test det(g) == 1
  @test discriminant(g) == 1
  @test g == g

  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E2; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E2, K(det(gram_matrix(g))) * inv(d), p)

  g = genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :disc)
  @test g == g
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  @test_throws ArgumentError genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :bla)
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E2; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E2, K(det(gram_matrix(g))) * inv(d), p)

  # -1 is a local square, so the determinant equals discriminant
  @test discriminant(g, 2) == -1
  @test discriminant(g) == -1
  @test det(g) == -1

  #
  # unramified & split & non-dyadic
  #

  p = prime_decomposition(OK, 3)[1][1]
  g = genus(HermLat, E2, p, [(0, 2, 1)], type = :disc)
  @test !isramified(g)
  @test !isinert(g)
  @test issplit(g)
  @test !isdyadic(g)
  @test_throws AssertionError norms(g)
  @test discriminant(g) == 1
  @test g == g
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E2; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E2, K(det(gram_matrix(g))) * inv(d), p)

  g = genus(HermLat, E2, p, [(0, 2, 1), (1, 1, 1), (2, 3, 1)], type = :disc)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E2; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E2, K(det(gram_matrix(g))) * inv(d), p)

  # Wrong det class in last component (must be 1)
  @test_throws ArgumentError g = genus(HermLat, E2, p, [(0, 2, 1), (1, 1, 1), (2, 3, -1)], type = :disc)

  #
  # ramified & non-dyadic and -1 not a local norm
  #
  K, _ = Hecke.rationals_as_number_field()
  Kt, t = K["t"]
  E, = number_field(t^2 + 7)
  p = prime_decomposition(maximal_order(K), 7)[1][1]

  # islocal_norm(E, -1, 7) == false

  g = genus(HermLat, E, p, [(0, 2, 1)], type = :disc)
  @test g == genus(HermLat, E, p, [(0, 2, -1)])
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  @test g == g
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E, K(det(gram_matrix(g))) * inv(d), p)

  g = genus(HermLat, E, p, [(0, 2, 1), (1, 4, 1), (2, 3, -1)], type = :disc)
  # swap because of the ranks
  @test g == genus(HermLat, E, p, [(0, 2, -1), (1, 4, 1), (2, 3, 1)])
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, "text/plain", g, context = :compact => true) isa String
  @test sprint(show, g) isa String
  @test sprint(show, g, context = :compact => true) isa String
  L = @inferred representative(g)
  @test L in g
  # p => g is cached on L
  @test hermitian_lattice(E; gram = gram_matrix(g)) in g
  d = @inferred det_representative(g)
  @test islocal_norm(E, K(det(gram_matrix(g))) * inv(d), p)

  for t in [(0, 1, 1), (0, 1, -1), (0, 2, 1), (0, 2, -1), (1, 2, -1)]
    g = genus(HermLat, E, p, [t])
    @test sprint(show, "text/plain", g) isa String
    @test sprint(show, "text/plain", g, context = :compact => true) isa String
    @test sprint(show, g) isa String
    @test sprint(show, g, context = :compact => true) isa String
    @test g == g
    L = @inferred representative(g)
    @test L in g
    # p => g is cached on L
    @test hermitian_lattice(E; gram = gram_matrix(g)) in g
    d = @inferred det_representative(g)
    @test islocal_norm(E, K(det(gram_matrix(g))) * inv(d), p)
  end

  #  odd scale, so rank even, and rank/2 is odd/even => det must be -1,1
  @test_throws ArgumentError genus(HermLat, E, p, [(1, 1, 1)])
  @test_throws ArgumentError genus(HermLat, E, p, [(1, 2, 1)])
  @test_throws ArgumentError genus(HermLat, E, p, [(1, 6, 1)])
  @test_throws ArgumentError genus(HermLat, E, p, [(1, 4, -1)])

  ##############################################################################
  #
  #  Now global
  #
  ##############################################################################

  Qx, x = QQ["x"]
  K, a = NumberField(x^2 - 2, "a")
  OK = maximal_order(K)
  Kt, t  = K["t"]

  E1, b1 = NumberField(t^2 - a, "b1") # ramified at 2
  E2, b2 = NumberField(t^2 - 5, "b2") # unramified at 2

  p = prime_decomposition(OK, 2)[1][1]
  q = prime_decomposition(OK, 3)[1][1]

  g = @inferred genus(HermLat, E1, p, [(0, 1, 1, 0), (2, 2, -1, 1)])
  L = representative(g)
  G = genus(L)
  @test sprint(show, "text/plain", G) isa String
  @test sprint(show, "text/plain", G, context = :compact => true) isa String
  @test sprint(show, G) isa String
  @test sprint(show, G, context = :compact => true) isa String
  g = @inferred genus(HermLat, E1, q, [(0, 2, 1)])
  M = representative(g)
  G2 = genus(M)
  @test genus(representative(orthogonal_sum(G, G2))) == orthogonal_sum(G, G2)
  @test genus(representative(orthogonal_sum(G2, G))) == orthogonal_sum(G2, G)

  rlp = real_places(K)
  sig = Dict(rlp[1] => 2, rlp[2] => 2)
  G = genus([g], sig)
  @test G == genus([g], [(rlp[1], 2), (rlp[2], 2)])

  ##############################################################################
  #
  #  Local genera hermitian
  #
  #############################################################################

  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = PolynomialRing(K, "t")
  L, b = number_field(t^2 - a * t + 1)

  p = prime_decomposition(maximal_order(K), 2)[1][1]
  G = @inferred local_genera_hermitian(L, p, 4, 2, 4)
  @test length(G) == 15
  for i in 1:length(G)
    @test rank(G[i]) == 4
    @test (@inferred representative(G[i])) in G[i]
  end

  for i in 1:10
    g1 = rand(G)
    g2 = rand(G)
    g3 = @inferred orthogonal_sum(g1, g2)
    @test genus(representative(g3), p) == genus(orthogonal_sum(representative(g1), representative(g2))[1], p)
  end

  GG = G[1]
  u = @inferred uniformizer(GG)
  @assert parent(u) == K

  p = prime_decomposition(maximal_order(K), 3)[1][1]
  G = local_genera_hermitian(L, p, 4, 2, 4)
  for i in 1:10
    g1 = rand(G)
    g2 = rand(G)
    g3 = @inferred orthogonal_sum(g1, g2)
    @test genus(representative(g3), p) == genus(orthogonal_sum(representative(g1), representative(g2))[1], p)
  end

  p = prime_decomposition(maximal_order(K), 17)[1][1]
  G = @inferred local_genera_hermitian(L, p, 5, 5, 5)
  @test length(G) == 7
  for i in 1:length(G)
    @test rank(G[i]) == 5
    @test (@inferred representative(G[i])) in G[i]
  end

  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = K["t"]
  L, b = number_field(t^2 - a * t + 1)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  l =  Vector{Tuple{Int, Int, Int, Int}}[[(0, 3, 1, 0), (4, 1, 1, 2)],
       [(0, 3, -1, 0), (4, 1, 1, 2)],
       [(0, 3, 1, 0), (4, 1, -1, 2)],
       [(0, 3, -1, 0), (4, 1, -1, 2)],
       [(0, 2, 1, 0), (2, 2, 1, 1)],
       [(0, 2, -1, 0), (2, 2, 1, 1)],
       [(0, 2, 1, 1), (2, 2, 1, 1)],
       [(0, 2, 1, 0), (2, 2, 1, 2)],
       [(0, 2, 1, 1), (2, 2, -1, 1)],
       [(0, 2, -1, 0), (2, 2, 1, 2)],
       [(0, 2, 1, 1), (2, 2, 1, 2)],
       [(0, 1, 1, 0), (1, 2, 1, 1), (2, 1, 1, 1)],
       [(0, 1, -1, 0), (1, 2, 1, 1), (2, 1, 1, 1)],
       [(1, 4, 1, 1)],
       [(1, 4, -1, 1)]]
  Gs = Hecke.LocalGenusHerm{typeof(L), typeof(p)}[ genus(HermLat, L, p, x) for x in l ]
  myG = @inferred local_genera_hermitian(L, p, 4, 2, 4)
  @test length(Gs) == length(myG)
  @test all(x -> x in Gs, myG)
  @test all(x -> x in myG, Gs)

  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = K["t"]
  L, b = number_field(t^2 - a * t + 1)
  rlp = real_places(K)
  G = @inferred genera_hermitian(L, 3, Dict(rlp[1] => 1, rlp[2] => 1), 100 * maximal_order(L))
  for i in 1:10
    g1 = rand(G)
    g2 = rand(G)
    M, = @inferred orthogonal_sum(representative(g1), representative(g2))
    @test M in (g1 + g2)
  end

end


text:
module TestFuncTests

using DiffBase
using Base.Test

n = rand()
x, y = rand(5, 5), rand(26)
A, B = rand(5, 5), rand(5, 5)

# f returns Number

for f in DiffBase.NUMBER_TO_NUMBER_FUNCS
    @test isa(f(n), Number)
end

for f in DiffBase.VECTOR_TO_NUMBER_FUNCS
    @test isa(f(y), Number)
end

for f in DiffBase.MATRIX_TO_NUMBER_FUNCS
    @test isa(f(x), Number)
end

for f in DiffBase.TERNARY_MATRIX_TO_NUMBER_FUNCS
    @test isa(f(A, B, x), Number)
end

# f returns Array

for f in DiffBase.NUMBER_TO_ARRAY_FUNCS
    @test isa(f(n), Array)
end

for f in DiffBase.ARRAY_TO_ARRAY_FUNCS
    @test isa(f(A), Array)
    @test isa(f(y), Array)
end

for f in DiffBase.MATRIX_TO_MATRIX_FUNCS
    @test isa(f(A), Array)
end

for f in DiffBase.BINARY_MATRIX_TO_MATRIX_FUNCS
    @test isa(f(A, B), Array)
end

# f! returns Void

for f! in DiffBase.INPLACE_ARRAY_TO_ARRAY_FUNCS
    @test isa(f!(y, x), Void)
end

for f! in DiffBase.INPLACE_NUMBER_TO_ARRAY_FUNCS
    @test isa(f!(y, n), Void)
end

end # module


text:
#=
   From BLOG example/galaxy.blog

   numCluster: 10.0817
   Distributions of variable origCluster1 (num:0)
   2.00000 =>    1149  (0.114900)
   1.00000 =>    1133  (0.113300)
   4.00000 =>    1108  (0.110800)
   3.00000 =>    1090  (0.109000)
   5.00000 =>     990  (0.099000)
   6.00000 =>     969  (0.096900)
   7.00000 =>     845  (0.084500)
   8.00000 =>     771  (0.077100)
   9.00000 =>     586  (0.058600)
   10.00000 =>     436  (0.043600)
   11.00000 =>     339  (0.033900)
   12.00000 =>     242  (0.024200)
   13.00000 =>     152  (0.015200)
   14.00000 =>      87  (0.008700)
   15.00000 =>      52  (0.005200)
   16.00000 =>      31  (0.003100)
   17.00000 =>      11  (0.001100)
   18.00000 =>       5  (0.000500)
   19.00000 =>       3  (0.000300)
   20.00000 =>       1  (0.000100)

   g1: 27623.544339264394
   g6: 27443.067174186628
   g10: 27617.655492791007
   g15: 27490.9282783939
   g16: 27577.219002335823
   g32: 27520.492984190976
   g51: 27606.068946048606
   g71: 27560.49783845427
   g81: 27531.756338358115

   Cf ~/webppl/galaxy.wppl
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function galaxy(data)
    n = length(data)

    numCluster ~ Truncated(Poisson(10),1,Inf)
    clusVelocity ~ filldist(Uniform(5000.0, 50000.0),n)
    origCluster ~ filldist(DiscreteUniform(1,numCluster),n)

    velocity = tzeros(n)
    for i in 1:n
        velocity[i] ~ Normal(clusVelocity[origCluster[i]], sqrt(10000))
    end

    for i in 1:n
       data[i] ~ Dirac(velocity[i])
    end

    origCluster1 ~ Dirac(origCluster[1])
    g1 ~ Dirac(clusVelocity[origCluster[1]])
    g6 ~ Dirac(clusVelocity[origCluster[6]])
    g10 ~ Dirac(clusVelocity[origCluster[10]])
    g15 ~ Dirac(clusVelocity[origCluster[15]])
    g16 ~ Dirac(clusVelocity[origCluster[16]])
    g32 ~ Dirac(clusVelocity[origCluster[32]])
    g51 ~ Dirac(clusVelocity[origCluster[51]])
    g71 ~ Dirac(clusVelocity[origCluster[71]])
    g81 ~ Dirac(clusVelocity[origCluster[81]])

end

data = [9172,9350,9483,9558,9775,10227,10406,16084,16170,18419,18552,18600,18927,19052,
        19070,19330,19343,19349,19440,19473,19529,19541,19547,19663,19846,19856,19863,
        19914,19918,19973,19989,20166,20175,20179,20196,20215,20221,20415,20629,20795,
        20821,20846,20875,20986,21137,21492,21701,21814,21921,21960,22185,22209,22242,
        22249,22314,22374,22495,22746,22747,22888,22914,23206,23241,23263,23484,23538,
        23542,23666,23706,23711,24129,24285,24289,24366,24717,24990,25633,26960,26995,
        32065,32789,34279];
println("len(data): ", length(data))
model = galaxy(data)

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

println("numCluster: ", mean_val(chns,:numCluster))
show_var_dist_pct(chns,:numCluster)
show_var_dist_pct(chns,:origCluster1)

println("g1: ", mean_val(chns,:g1))
println("g6: ", mean_val(chns,:g6))
println("g10: ", mean_val(chns,:g10))
println("g15: ", mean_val(chns,:g15))
println("g16: ", mean_val(chns,:g16))
println("g32: ", mean_val(chns,:g32))
println("g51: ", mean_val(chns,:g51))
println("g71: ", mean_val(chns,:g71))
println("g81: ", mean_val(chns,:g81))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text: @deprecate SVGP SparseVariationalApproximation


text:
@testset "Inital params" begin
    @testset "stationary initial params" begin
        gas_1_1 = Model(1, 1, Normal, 0.0)
        gas_1_1.ω = [1; 1]
        gas_1_1.B[1] = [0.8 0;0 0.5]
        initial_params_tilde = stationary_initial_params_tilde(gas_1_1)
        initial_params = stationary_initial_params(gas_1_1)
        @test initial_params_tilde[1, 1] ≈ 5 atol = 1e-2 rtol = 1e-2
        @test initial_params_tilde[1, 2] ≈ 2 atol = 1e-2 rtol = 1e-2
        @test initial_params[1, 1] ≈ 5 atol = 1e-2 rtol = 1e-2
        @test initial_params[1, 2] ≈ exp(2) atol = 1e-2 rtol = 1e-2

        gas_1_2 = Model(1, 2, Normal, 0.0)
        gas_1_2.ω = [1; 1]
        gas_1_2.B[1] = [0.4 0;0 0.25]
        gas_1_2.B[2] = [0.4 0;0 0.25]
        initial_params_tilde = stationary_initial_params_tilde(gas_1_2)
        initial_params = stationary_initial_params(gas_1_2)
        @test initial_params_tilde[1, 1] ≈ 5 atol = 1e-2 rtol = 1e-2
        @test initial_params_tilde[1, 2] ≈ 2 atol = 1e-2 rtol = 1e-2
        @test initial_params_tilde[2, 1] ≈ 5 atol = 1e-2 rtol = 1e-2
        @test initial_params_tilde[2, 2] ≈ 2 atol = 1e-2 rtol = 1e-2
        @test initial_params[1, 1] ≈ 5 atol = 1e-2 rtol = 1e-2
        @test initial_params[1, 2] ≈ exp(2) atol = 1e-2 rtol = 1e-2
        @test initial_params[2, 1] ≈ 5 atol = 1e-2 rtol = 1e-2
        @test initial_params[2, 2] ≈ exp(2) atol = 1e-2 rtol = 1e-2
    end
    @testset "Dynamic initial params" begin
        test_dynamic(0.2, 12)
        test_dynamic(0.1, 3)
    end
end


text:
abstract type AbstractAttenOp end
abstract type AbstractAttenScoreOp end
abstract type AbstractMixingOp end

struct DotProductScore <: AbstractAttenScoreOp end
(::DotProductScore)(args...) = dot_product_score(args...)

struct ScaledDotProductScore <: AbstractAttenScoreOp end
(::ScaledDotProductScore)(args...) = scaled_dot_product_score(args...)

struct NormalizedScoreOp{S, N} <: AbstractAttenScoreOp
    score::S
    norm::N
end
(op::NormalizedScoreOp)(args...) = normalized_score(op.norm, op.score, args...)

struct MaskedScoreOp{M, S} <: AbstractAttenScoreOp
    maskop::M
    score::S
end
MaskedScoreOp(score) = MaskedScoreOp(nothing, score)
(op::MaskedScoreOp)(mask, args...) = masked_score(op.score, op.maskop, mask, args...)

struct WeightedSumMixing <: AbstractMixingOp end
(op::WeightedSumMixing)(args...) = weighted_sum_mixing(args...)

struct NaiveQKVAttenOp <: AbstractAttenOp end
(::NaiveQKVAttenOp)(args...) = naive_qkv_attention(args...)

struct MultiheadQKVAttenOp <: AbstractAttenOp
    head::Int
end
(op::MultiheadQKVAttenOp)(args...) = multihead_qkv_attention(op.head, args...)


text:
using Documenter, NaturalSelection

makedocs(
    modules = [NaturalSelection],
    format = :html,
    sitename = "NaturalSelection.jl",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "dNdS" => "man/dNdS.md",
            "Tajima's D" => "man/tajimad.md"
        ],
        "Contributing" => "contributing.md"
    ],
    authors = "Ben J. Ward, The BioJulia organisation, and contributors."
)

deploydocs(
    repo = "github.com/BioJulia/NaturalSelection.jl.git",
    julia = "0.6",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
# Note: poor performance

# this may not be working properly

cd(@__DIR__)
using DifferentialEvolutionMCMC, Random, StatsBase, Distributions

N = 30
# Parameters for each cluster, we assume that each cluster is Gaussian distributed in the example.
μs = [-3.5, 0.0]
# Construct the data points.
data = mapreduce(c -> rand(MvNormal([μs[c], μs[c]], 1.), N), hcat, 1:2)

priors = (
    idx = (Categorical([.5,.5]),N),
    μs = (Normal(0, 1),2),
)

bounds = ((1,2),(0.0,Inf))

function loglike(data, idx, μs)
    LL = 0.0
    for (i,id) in enumerate(idx)
        LL += logpdf(MvNormal([μs[id], μs[id]], 1.0), data[:,i])
    end
    return LL
end

model = DEModel(priors=priors, model=loglike, data=data)

de = DE(bounds=bounds, burnin=1000, priors=priors, Np=15)
n_iter = 2000
@elapsed chains = sample(model, de, n_iter, progress=true)


text:
module IrisApp

using Matte

const title = "Matte Example App :: Iris Dataset"

function ui()
    sidebar_layout(
        side_panel(
            h1("What species would you like?"),
            br(),
            radio("species", ["setosa","versicolor","virginica"], ["Setosa","Versicolor","Virginica"])
        ),
        main_panel(
            h1("Iris data: ", text_output("selected_species")),
            br(),
            datatable_output("iris_df")
        )
    )
end

module Server

using RDatasets
using DataFrames
using InvertedIndices
iris = dataset("datasets", "iris")

function iris_df(species)
    if species != nothing
        iris[iris.Species .== species, Not(:Species)]
    else
        nothing
    end
end

function selected_species(species)
    if species == nothing
        "No species selected"
    else
        species
    end
end

end

end


text:
"""
    systematic_sample([rng::AbstractRNG,] N, d=Normal(0,1); permute=true)

returns a `Vector` of length `N` sampled systematically from the distribution `d`. If `permute=false`, this vector will be sorted.
"""
function systematic_sample(rng::AbstractRNG, N, d=Normal(0,1); permute=true)
    T = eltype(d)
    e = T(0.5/N) # rand()/N
    y = e:1/N:1
    o = quantile.(d,y)
    permute && permute!(o, randperm(rng, N))
    return eltype(o) == T ? o : T.(o)
end

function systematic_sample(N, d=Normal(0,1); kwargs...)
    return systematic_sample(Random.GLOBAL_RNG, N, d; kwargs...)
end

"""
    ess(p::AbstractParticles{T,N})

Calculates the effective sample size. This is useful if particles come from MCMC sampling and are correlated in time. The ESS is a number between [0,N].

Initial source: https://github.com/tpapp/MCMCDiagnostics.jl
"""
function ess(p::AbstractParticles)
    ac    = autocor(p.particles,1:min(250, length(p)÷2))
    N     = length(ac)
    τ_inv = 1 + 2ac[1]
    K     = 2
    while K < N - 2
        Δ = ac[K] + ac[K+1]
        Δ < 0 && break
        τ_inv += 2Δ
        K += 2
    end
    min(1 / τ_inv, one(τ_inv))*length(p)
end


text:
function besselroots(nu::Float64, n::Int64)
#BESSELROOTS    The first N roots of the function J_v(x)

# DEVELOPERS NOTES:
#   V = 0 --> Full double precision for N <= 20 (Wolfram Alpha), and very
#     accurate approximations for N > 20 (McMahon's expansion)
#   -1 <= V <= 5 : V ~= 0 --> 12 decimal figures for the 6 first zeros
#     (Piessens's Chebyshev series approximations), and very accurate
#     approximations for the others (McMahon's expansion)
#   V > 5 --> moderately accurate for the 6 first zeros and good
#     approximations for the others (McMahon's expansion)

# This code was originally written by L. L. Peixoto in MATLAB.
# Later modified by A. Townsend to work in Julia

if ( n == 0 )
    x = []
elseif( n > 0 && nu == 0 )
    x = McMahon( nu, n )
    correctFirstFew = Tabulate( )
    x[1:min(n,20)] = correctFirstFew[1:min(n,20)]
    x
elseif ( n>0 && nu >= -1 && nu <= 5 )
    x = McMahon( nu, n )
    correctFirstFew = Piessens( nu )
    x[1:min(n,6)] = correctFirstFew[1:min(n,6)]
    x
else
    error("Input N must be a positive integer")
end

end

function McMahon( nu::Float64, n::Int64 )
# McMahon's expansion. This expansion gives very accurate approximation
# for the sth zero (s >= 7) in the whole region V >=- 1, and moderate
# approximation in other cases.
b = Array(Float64,n)
mu = 4*nu^2
a1 = 1. / 8.
a3 = (7.*mu-31.) / 384.
a5 = 4.*(3779.+mu*(-982.+83.*mu)) / 61440. # Evaluate via Horner's method.
a7 = 6.*(-6277237.+mu*(1585743.+mu*(-153855.+6949.*mu))) / 20643840.;
a9 = 144.*(2092163573.+mu*(-512062548.+mu*(48010494.+mu*(-2479316.+70197.*mu)))) / 11890851840.
a11 = 720.*(-8249725736393.+mu*(1982611456181.+mu*(-179289628602.+mu*(8903961290. +
      mu*(-287149133.+5592657.*mu))))) / 10463949619200.
a13 = 576.*(423748443625564327. + mu*(-100847472093088506.+mu*(8929489333108377. +
    mu*(-426353946885548.+mu*(13172003634537.+mu*(-291245357370. + mu*4148944183.)))))) / 13059009124761600.
for k=1:n b[k] = .25*(2*nu+4*k-1)*pi end # beta
# Evaluate using Horner's scheme:
x = b - (mu-1)*( ((((((a13./b.^2 + a11)./b.^2 + a9)./b.^2 + a7)./b.^2 + a5)./b.^2 + a3)./b.^2 + a1)./b)
end

function Tabulate( )
# First 20 roots of J0(x) are precomputed (using Wolfram Alpha):
y = Array(Float64,20)
y = [   2.4048255576957728
        5.5200781102863106
        8.6537279129110122
        11.791534439014281
        14.930917708487785
        18.071063967910922
        21.211636629879258
        24.352471530749302
        27.493479132040254
        30.634606468431975
        33.775820213573568
        36.917098353664044
        40.058425764628239
        43.199791713176730
        46.341188371661814
        49.482609897397817
        52.624051841114996
        55.765510755019979
        58.906983926080942
        62.048469190227170  ]
end

function Piessens( nu::Float64 )
# Piessens's Chebyshev series approximations (1984). Calculates the 6 first
# zeros to at least 12 decimal figures in region -1 <= V <= 5:
    y = Array(Float64,6);

    C = [
       2.883975316228  8.263194332307 11.493871452173 14.689036505931 17.866882871378 21.034784308088
       0.767665211539  4.209200330779  4.317988625384  4.387437455306  4.435717974422  4.471319438161
      -0.086538804759 -0.164644722483 -0.130667664397 -0.109469595763 -0.094492317231 -0.083234240394
       0.020433979038  0.039764618826  0.023009510531  0.015359574754  0.011070071951  0.008388073020
      -0.006103761347 -0.011799527177 -0.004987164201 -0.002655024938 -0.001598668225 -0.001042443435
       0.002046841322  0.003893555229  0.001204453026  0.000511852711  0.000257620149  0.000144611721
      -0.000734476579 -0.001369989689 -0.000310786051 -0.000105522473 -0.000044416219 -0.000021469973
       0.000275336751  0.000503054700  0.000083834770  0.000022761626  0.000008016197  0.000003337753
      -0.000106375704 -0.000190381770 -0.000023343325 -0.000005071979 -0.000001495224 -0.000000536428
       0.000042003336  0.000073681222  0.000006655551  0.000001158094  0.000000285903  0.000000088402
      -0.000016858623 -0.000029010830 -0.000001932603 -0.000000269480 -0.000000055734 -0.000000014856
       0.000006852440  0.000011579131  0.000000569367  0.000000063657  0.000000011033  0.000000002536
      -0.000002813300 -0.000004672877 -0.000000169722 -0.000000015222 -0.000000002212 -0.000000000438
       0.000001164419  0.000001903082  0.000000051084  0.000000003677  0.000000000448  0.000000000077
      -0.000000485189 -0.000000781030 -0.000000015501 -0.000000000896 -0.000000000092 -0.000000000014
       0.000000203309  0.000000322648  0.000000004736  0.000000000220  0.000000000019  0.000000000002
      -0.000000085602 -0.000000134047 -0.000000001456 -0.000000000054 -0.000000000004               0.
       0.000000036192  0.000000055969  0.000000000450  0.000000000013               0.               0.
      -0.000000015357 -0.000000023472 -0.000000000140 -0.000000000003               0.               0.
       0.000000006537  0.000000009882  0.000000000043  0.000000000001               0.               0.
      -0.000000002791 -0.000000004175 -0.000000000014               0.               0.               0.
       0.000000001194  0.000000001770  0.000000000004               0.               0.               0.
      -0.000000000512 -0.000000000752              0.               0.              0.               0.
       0.000000000220  0.000000000321               0.               0.               0.               0.
      -0.000000000095 -0.000000000137               0.               0.               0.               0.
       0.000000000041  0.000000000059               0.               0.               0.               0.
      -0.000000000018 -0.000000000025               0.               0.               0.               0.
       0.000000000008  0.000000000011               0.               0.               0.               0.
      -0.000000000003 -0.000000000005               0.               0.               0.               0.
       0.000000000001  0.000000000002               0.               0.               0.               0.]

    T = Array(Float64,size(C,1))
    pt = (nu-2)/3
    T[1], T[2] = 1., pt
    for k = 2:size(C,1)-1
        T[k+1] = 2*pt*T[k] - T[k-1]
    end
    y[1:6] = T'*C;
    y[1] *= sqrt(nu+1);                  # Scale the first root.
    y
end


text:
immutable BlendMode
    color_src_factor::Cint
    color_dst_factor::Cint
    color_equation::Cint
    alpha_src_factor::Cint
    alpha_dst_factor::Cint
    alpha_equation::Cint
end

const blend_alpha = BlendMode(6, 7, 0, 1, 7, 0)
const blend_add = BlendMode(6, 1, 0, 1, 1, 0)
const blend_multiply = BlendMode(4, 0, 0, 4, 0, 0)
const blend_none = BlendMode(1, 0, 0, 1, 0, 0)


text:
# ##
# Surface Quasi-Geostrophy:
# This code solves an exmaple problem from Held et al JFM 1995
# Imposed condition is the 'surface buoyancy' given by a smooth elliptic vortex.
# This equations are Eq.(2) of the paper.

using FourierFlows, Plots, Statistics, Printf, Random

using FFTW: irfft, rfft
using Statistics: mean
using Random: seed!

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum

using CUDA
#    Reexport

#nothingfunction(args...) = nothing

# parameters
#
# Let's define a struct which contain the
# main parameters of the problem
# struct Params{T} <: AbstractParams
#      ν :: T          # hyperviscosity coefficient
#     nν :: Int       # Order of the hyperviscosity
# end
# nothing #hide

struct Params{T} <: AbstractParams
    ν :: T         # Hyperviscosity coefficient
   nν :: Int       # Order of the hyperviscous operator
end
nothing #hide

# variables
# conatins all the problem variables in physical
# and Fourier (transformed) space :  Aphys, Atrans
struct Vars{Aphys, Atrans} <: AbstractVars
    b  :: Aphys
    u  :: Aphys
    v  :: Aphys
    bh :: Atrans
    uh :: Atrans
    vh :: Atrans
end
nothing #hide

# constructor creates the size of the variables
# based on the 'grids'
function Vars(::Dev, grid) where Dev
    T = eltype(grid)
    @devzeros Dev T (grid.nx, grid.ny) b u v
    @devzeros Dev Complex{T} (grid.nkr, grid.nl) bh uh vh
    return Vars(b, u, v, bh, uh, vh)
end
nothing #hide

# constructing the equation:
# here the linear part is described ny hyperviscosity part,
# and the nonlinear part is computed by the 'calcN!'
# the equations are in Fourier space
function Equation(paramas::Params, grid::AbstractGrid)
    #T = eltype(grid)
    L = @. -paramas.ν*grid.Krsq^paramas.nν
    CUDA.@allowscalar L[1,1] = 0
    return FourierFlows.Equation(L, calcN!, grid)
end
nothing #hide

# constructing the nonlinear part of the problem
# ofcourse in Fourier space.
function calcN!( N, sol, t, clock, vars, params, grid )
    @. vars.bh = sol
    @. vars.uh =   im * grid.l  * sqrt( grid.invKrsq ) * sol
    @. vars.vh = - im * grid.kr * sqrt( grid.invKrsq ) * sol

    # doing inverse-fft (I suppose!!)
    ldiv!( vars.u, grid.rfftplan, vars.uh )
    ldiv!( vars.v, grid.rfftplan, vars.vh )
    ldiv!( vars.b, grid.rfftplan, vars.bh )

    ub, ubh = vars.u, vars.uh
    vb, vbh = vars.v, vars.vh

    @. ub *= vars.b
    @. vb *= vars.b

    # compute the advection term: u*bx + v*by
    mul!( ubh, grid.rfftplan, ub )
    mul!( vbh, grid.rfftplan, vb )

    @. N = - im * grid.kr * ubh - im * grid.l * vbh
    return nothing
end
nothing #hide

# All are ready to go, my master!
# Hold on to the moment my apprentice, let me create some helpers to do thing smoothly.

# function which updates the variables in current time!
function updatevars!(prob)
    vars, grid, sol = prob.vars, prob.grid, prob.sol

    @. vars.bh = sol
    @. vars.uh =   im * grid.l  * sqrt( grid.invKrsq ) * sol
    @. vars.vh = - im * grid.kr * sqrt( grid.invKrsq ) * sol

    # doing inverse-fft
    ldiv!( vars.u, grid.rfftplan, deepcopy(vars.uh) )
    ldiv!( vars.v, grid.rfftplan, deepcopy(vars.vh) )
    ldiv!( vars.b, grid.rfftplan, deepcopy(vars.bh) )

    return nothing
end

# set buoyancy as the solution and upadate all other variables.
function set_b!(prob, b)
    mul!( prob.sol, prob.grid.rfftplan, b )
    CUDA.@allowscalar prob.sol[1,1] = 0  # domain average gives zero

    updatevars!(prob)

    return nothing
end

# calculating surface-averaged kinetic energy.
# In SQG formulation, this is exactly half of the buoyancy variance.
@inline function kinetic_energy(prob)

    vars, grid, sol = prob.vars, prob.grid, prob.sol

    @. vars.uh =   im * grid.l  * sqrt( grid.invKrsq ) * sol
    @. vars.vh = - im * grid.kr * sqrt( grid.invKrsq ) * sol

    ke_h = vars.uh

    @. ke_h = 0.5*( abs2(vars.uh) + abs2(vars.vh) )

    return 1/(grid.Lx * grid.Ly) * parsevalsum( ke_h, grid )
end

# ## Choosing a device: CPU or GPU
dev = CPU()    # Device (CPU/GPU)
nothing # hide

# numerical parameters
     nx = 512                           # grid resolution
     ny = 512
stepper = "FilteredRK4"                 # time-stepping scheme
     dt = 0.01                          # timesep (in seconds)
     Tf = 80.                           # total time of the simulation
 nsteps = Int( Tf/dt )                  # total iteration of the simulation
  nplot = 10 #round( Int, nsteps/100 )      # step to save data
nothing # hide

# physical parameters
     Lx = 2π                # domain size
     Ly = 2π
     nν = 4
      ν = 1.e-18
nothing # hide

# creating the grids struct
grid = TwoDGrid( dev, nx, Lx, ny, Ly; T=Float64)

# creating parameters struct
params = Params(ν, nν)

vars = Vars(dev, grid)
equation = Equation(params, grid)
prob = FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)

# Let's create some shortcuts.
sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = prob.grid.x, prob.grid.y
nothing # hide

# Here comes the initial ccondition, and set to fly.
X, Y = gridpoints(grid)
b0 = @. exp( -(X^2 + 4*Y^2) )

set_b!(prob, b0)
nothing #hide

#### plot the initial condition
heatmap(x, y, prob.vars.b',
    aspectratio = 1,
             c  = :deep,
           clim = [0, 1],
           xlim = ( -grid.Lx/2, grid.Lx/2 ),
           ylim = ( -grid.Ly/2, grid.Ly/2 ),
         xticks = -3:3,
         yticks = -3:3,
         xlabel = "x",
         ylabel = "y",
          title = "Initial distributon of buoyancy",
     framestyle = :box)

# creaing diganosis
KE = Diagnostic(kinetic_energy, prob; nsteps=nsteps)
diags = KE
nothing #

### ouput
base_filename = string( "surfaceqg_Held_", nx )
datapath = "./"
plotpath = "./"

dataname = joinpath(datapath, base_filename)
plotname = joinpath(datapath, base_filename)
nothing #hide

if !isdir(plotpath); mkdir(plotpath); end
if !isdir(datapath); mkdir(datapath); end
nothing #hide

# get the output
get_sol(prob) = sol  # extract the fft-solution
get_u(prob) = irfft( im * grid.l  .* sqrt( grid.invKrsq ) .* sol, grid.nx )
out = Output(prob, dataname, (:sol, get_sol), (:u, get_u))
nothing # hide

### Visulation
function plot_output(prob)
    br = prob.vars.b
    ur = prob.vars.u
    vr = prob.vars.v

    pb = heatmap(x, y, br',
         aspectratio = 1,
                   c = :Reds,
                clim = (0, 1),
               xlims = (-grid.Lx/2, grid.Lx/2),
               ylims = (-grid.Ly/2, grid.Ly/2),
              xticks = -3:3,
              yticks = -3:3,
              xlabel = "x",
              ylabel = "y",
               title = "buoyancy",
          framestyle = :box)

    # pKE = plot(1,
    #            label = "kinetic energy",
    #        linewidth = 2,
    #           legend = :bottomright,
    #            alpha = 0.7,
    #            xlims = (0, Tf),
    #            ylims = (0, 1e-2),
    #           xlabel = "t")

    # pb² = plot(1,
    #            label = "buoyancy variance ∫bₛ²dxdy/L²",
    #        linecolor = :red,
    #           legend = :bottomright,
    #        linewidth = 2,
    #            alpha = 0.7,
    #            xlims = (0, tf),
    #            ylims = (0, 2e-2),
    #           xlabel = "t")

    # layout = @layout [a{0.5w} Plots.grid(2, 1)]
    # layout = @layout grid(2, 1)
    p = plot(pb, size = (500, 500), dpi=350)

    return p
end
nothing # hide

# Time-stepping forward in time

startwalltime = time()

p = plot_output(prob)

anim = @animate for j=0:Int(nsteps/nplot)

    # CFL condition
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    #if j%(500/nplot)==0
    if j%nplot==0
        log1 = @sprintf("step: %04d, t: %.1f, cfl: %.3f, walltime: %.2f min",
            clock.step, clock.t, cfl, (time()-startwalltime)/60)

        log2 = @sprintf("kinetic energy: %.2e", KE.data[KE.i])

        println(log1)
        println(log2)
    end

    p[1][1][:z] = vars.b
    p[1][:title] = "buoyancy, t="*@sprintf("%.2f", clock.t)
    #push!(p[2][1], KE.t[KE.i], KE.data[KE.i])
    #push!(p[3][1], B.t[B.i], B.data[B.i])

    stepforward!(prob, diags, nplot)
    updatevars!(prob)

end

mp4(anim, "SQG_Held1995.mp4", fps=8)


text:
# IIF testCompareVariableFactors repeated in IIF

using DistributedFactorGraphs
using IncrementalInference
using Test

@testset "testing compare functions for variables and factors..." begin

    fg = initfg()

    v1 = addVariable!(fg, :x0, ContinuousScalar)
    f1 = addFactor!(fg, [:x0;], Prior(Normal()))

    @test compareVariable(v1,v1)
    @test compareFactor(f1,f1)

    v2 = addVariable!(fg, :x1, ContinuousScalar)
    f2 = addFactor!(fg, [:x0;:x1], LinearRelative(Normal(2.0, 0.1)))

    fg2 = deepcopy(fg)

    @test !compareVariable(v1,v2)

    # not testing different factors in this way
    # @test !compareFactor(f1,f2)

    @test compareAllVariables(fg, fg)

    @test compareAllVariables(fg, fg2)

    @test compareSimilarVariables(fg, fg)
    @test compareSimilarVariables(fg, fg2)

    @test compareSimilarFactors(fg, fg)
    @test compareSimilarFactors(fg, fg2)

    @test compareFactorGraphs(fg, fg)
    @test compareFactorGraphs(fg, fg2)

    tree = solveTree!(fg)

    x1a = getVariable(fg, :x0)
    x1b = getVariable(fg2, :x0)

    @test !compareVariable(x1a, x1b, skipsamples=false)

    @test !compareSimilarVariables(fg, fg2, skipsamples=false)
    @test !compareSimilarFactors(fg, fg2, skipsamples=false)

    @test compareFactorGraphs(fg, fg)
    @test !compareFactorGraphs(fg, fg2, skipsamples=false)

    ensureAllInitialized!(fg2)

    @test compareSimilarVariables(fg, fg2, skipsamples=true, skip=Symbol[:initialized;:inferdim;:ppeDict;:solvedCount])
    # fg2 has been solved, so it should fail on the estimate dictionary
    @test !compareSimilarVariables(fg, fg2, skipsamples=true, skip=Symbol[:initialized;:inferdim])

    tree = buildTreeReset!(fg2)

    @test compareSimilarFactors(fg, fg2, skipsamples=true, skipcompute=true)

    @test !compareSimilarFactors(fg, fg2, skipsamples=true, skipcompute=false)

    @test compareFactorGraphs(fg, fg2, skipsamples=true, skipcompute=true, skip=[:initialized;:inferdim;:ppeDict; :solvedCount])

end

@testset "test subgraph functions..." begin

    fg = initfg()

    addVariable!(fg, :x0, ContinuousScalar)
    addFactor!(fg, [:x0;], Prior(Normal()))

    addVariable!(fg, :x1, ContinuousScalar)
    addFactor!(fg, [:x0;:x1], LinearRelative(Normal(2.0, 0.1)))

    addVariable!(fg, :x2, ContinuousScalar)
    addFactor!(fg, [:x1;:x2], LinearRelative(Normal(4.0, 0.1)))

    addVariable!(fg, :l1, ContinuousScalar)
    addFactor!(fg, [:x1;:l1], LinearRelative(Rayleigh()))

    sfg = buildSubgraph(LightDFG, fg, [:x0;:x1])

    @warn "FIXME This is NOT supposed to pass"
    @test_skip compareFactorGraphs(fg, sfg, skip=[:labelDict;:addHistory;:logpath])
    # drawGraph(sfg)

end


text:
"""
A package for handling observables in a Markov Chain Monte Carlo simulation.

See http://github.com/crstnbr/MonteCarloObservable.jl for more information.
"""
module MonteCarloObservable

    using Statistics
    using JLD, EllipsisNotation, Lazy, Reexport
    import HDF5

    @reexport using BinningAnalysis

    abstract type AbstractObservable end

    include("helpers.jl")
    include("shared.jl")
    include("observable.jl")
    include("lightobservable.jl")

    # general
    export Observable, DiskObservable, LightObservable
    export @obs, @diskobs

    # statistics
    export tau, iswithinerrorbars
    export std_error
    export jackknife
    export mean, var, std

    # interface
    export push!, append!, reset!
    export timeseries, ts
    export rename!, name
    export inmemory, isinmemory, length, eltype, getindex, view, isempty, ndims, size, iterate

    # io
    export saveobs, loadobs, listobs, rmobs
    export export_result, export_error
    export loadobs_frommemory
    export timeseries_frommemory, timeseries_frommemory_flat, mean_frommemory, error_frommemory
    export timeseries_flat, ts_flat
    export getfrom
    export flush
    export ObservableResult, load_result

    # deprecations
    include("deprecated.jl")
end


text:
abstract type AbstractMedium{T<:Real} end

"""
    refractive_index(medium :: AbstractMedium{T}, λ :: T) where T

Return the refractive index of medium at the given wavelength λ in μm.
"""
function refractive_index end

"""
    ConstantIndex{T}

Definition of a medium with no dispersion.

# Fields
- `n :: T`: The refractive index of the medium
"""
struct ConstantIndex{T} <: AbstractMedium{T}
    n::T
end

function refractive_index(medium::ConstantIndex{T}, λ::T) where {T}
    medium.n
end

"""
    Sellmeier_1{T}

Definition of a material's dispersion using the Sellmeier equation. This definition is
equivalent to the Zemax Sellmeier_1 dispersion equation.

``n^2(λ) = 1 + \\frac{k_1 λ^2}{λ^2 - l_1} + \\frac{k_2 λ^2}{λ^2 - l_1} + \\frac{k_3 λ^2}{λ^2 - l_3}``
"""
struct Sellmeier_1{T} <: AbstractMedium{T}
    k1::T
    l1::T
    k2::T
    l2::T
    k3::T
    l3::T
end

function refractive_index(medium::Sellmeier_1{T}, λ::T) where {T}
    λ2 = λ^2
    a = medium.k1 * λ2 / (λ2 - medium.l1)
    b = medium.k2 * λ2 / (λ2 - medium.l2)
    c = medium.k3 * λ2 / (λ2 - medium.l3)
    sqrt(a + b + c + 1)
end

#
# Definition of a few materials
#

"""
    Silica

Dispersion of fused silica, defined as Sellmeier equation coefficients.
The coefficients match Zemax.
"""
const Silica = Sellmeier_1(
    6.96166300e-1,
    4.67914800e-3,
    4.07942600e-1,
    1.35120600e-2,
    8.974794e-1,
    9.7934e+1,
)

"""
    Air

The refractive index of air is defined to be 1.0 at all wavelengths. This follows the usual
convention to define refractive indices relative to air.
"""
const Air = ConstantIndex(1.0)


text:
pages = Any[
         "DifferentialEquations.jl: Scientific Machine Learning (SciML) Enabled Simulation and Estimation" => "index.md",
         "Tutorials" => Any[
           "tutorials/ode_example.md",
           "tutorials/faster_ode_example.md",
           "tutorials/advanced_ode_example.md",
           "tutorials/sde_example.md",
           "tutorials/rode_example.md",
           "tutorials/dde_example.md",
           "tutorials/dae_example.md",
           "tutorials/discrete_stochastic_example.md",
           "tutorials/jump_diffusion.md",
           "tutorials/bvp_example.md",
           "tutorials/additional.md"
         ],
         "Basics" => Any[
           "basics/overview.md",
           "basics/common_solver_opts.md",
           "basics/solution.md",
           "basics/plot.md",
           "basics/integrator.md",
           "basics/problem.md",
           "basics/faq.md",
           "basics/compatibility_chart.md"
         ],
         "Problem Types" => Any[
           "types/discrete_types.md",
           "types/ode_types.md",
           "types/nonautonomous_linear_ode.md",
           "types/dynamical_types.md",
           "types/split_ode_types.md",
           "types/steady_state_types.md",
           "types/bvp_types.md",
           "types/sde_types.md",
           "types/sdae_types.md",
           "types/rode_types.md",
           "types/dde_types.md",
           "types/sdde_types.md",
           "types/dae_types.md",
           "types/jump_types.md",
         ],
         "Solver Algorithms" => Any[
           "solvers/discrete_solve.md",
           "solvers/ode_solve.md",
           "solvers/nonautonomous_linear_ode.md",
           "solvers/dynamical_solve.md",
           "solvers/split_ode_solve.md",
           "solvers/steady_state_solve.md",
           "solvers/bvp_solve.md",
           "solvers/jump_solve.md",
           "solvers/sde_solve.md",
           "solvers/sdae_solve.md",
           "solvers/rode_solve.md",
           "solvers/dde_solve.md",
           "solvers/sdde_solve.md",
           "solvers/dae_solve.md",
           "solvers/benchmarks.md"
         ],
         "Additional Features" => Any[
           "features/performance_overloads.md",
           "features/diffeq_arrays.md",
           "features/diffeq_operator.md",
           "features/noise_process.md",
           "features/linear_nonlinear.md",
           "features/callback_functions.md",
           "features/callback_library.md",
           "features/ensemble.md",
           "features/io.md",
           "features/low_dep.md",
           "features/progress_bar.md"
         ],
         "Analysis Tools" => Any[
           "analysis/parameterized_functions.md",
           "analysis/parameter_estimation.md",
           "analysis/bifurcation.md",
           "analysis/sensitivity.md",
           "analysis/global_sensitivity.md",
           "analysis/uncertainty_quantification.md",
           "analysis/neural_networks.md",
           "analysis/dev_and_test.md"
         ],
         "Domain Modeling Tools" => Any[
             "models/multiscale.md",
             "models/physical.md",
             "models/financial.md",
             "models/chemical_reactions.md",
             "models/external_modeling.md"
         ],
         "Extra Details" => Any[
             "extras/timestepping.md",
         ]
         ]

max_stars_repo_name: JuliaDiffEq/DiffEqDocs

max_issues_repo_name: JuliaDiffEq/DiffEqDocs

max_forks_repo_name: JuliaDiffEq/DiffEqDocs


text:
#///////////////////////////////////////
#// File Name: optimal_control_problems.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2020/10/28
#// Description: Optimal control problem definitions.
#///////////////////////////////////////
"""
    OptimalControlProblem

Abstract base type for an Optimal Control Problem.
"""
abstract type OptimalControlProblem end

# Optimal Control Problem for Risk-Sensitive iLEQG (or Risk-Neutral iLQG)
"""
    FiniteHorizonRiskSensitiveOptimalControlProblem(f, c, h, W, N) <: OptimalControlProblem

A finite horizon, stochastic optimal control problem where the dynamics function is subject to additive Gaussian noise w ~ N(0, W).

# Arguments
- `f(x, u, f_returns_jacobian=false)` -- deterministic dynamics function
  - `x` is a state vector and `u` is a control input vector.
  - The third positional argument `f_returns_jacobian` determines whether the user computes
    and returns the Jacobians, and should default to `false`. If `true`, the return value
    must be augmented with matrices `A` and `B`, where `A = dx_next/dx` and `B = dx_next/du`.
    Otherwise the return value is the (noiseless) next state `x_next`.
- `c(k, x, u)` -- stage cost function
  - `k::Int >= 0` is a time index where `k == 0` is the initial time.
  -  We assume that `c` is non-negative.
- `h(x)` -- terminal cost function
  -  We assume that `h` is non-negative.
- `W(k)` -- covariance matrix function
  - Returns a symmetric positive semidefinite matrix that represents the covariance matrix for
    additive Gaussian noise w ~ N(0, W).
  - `k::Int >= 0` is a time index where `k == 0` is the initial time.
- `N::Int64` -- final time index
  - Note that `0` is the initial time index.

# Notes
- Functions `f`, `c`, and `h` should be written generically enough to accept
  the state `x` and the input `u` of type `Vector{<:Real}`. This is to ensure that
  ForwardDiff can compute Jacobians and Hessians for iLQG/iLEQG.

# Example
```julia
import LinearAlgebra;

function f(x, u, f_returns_jacobian=false)
    x_next = x + u; # 2D single integrator dynamics
    if f_returns_jacobian
        A = Matrix(1.0LinearAlgebra.I, 2, 2); # dx_next/dx
        B = Matrix(1.0LinearAlgebra.I, 2, 2); # dx_next/du
        return x_next, A, B;
    else
        return x_next;
    end
end

c(k, x, u) = k/2*x'*x + k/2*u'*u  # time-dependent quadratic stage cost
N = 10;
h(x) = N/2*x'*x; # quadratic terminal cost
W(k) = Matrix(0.1LinearAlgebra.I, 2, 2);

problem = FiniteHorizonRiskSensitiveOptimalControlProblem(f, c, h, W, N);
```
"""
mutable struct FiniteHorizonRiskSensitiveOptimalControlProblem <: OptimalControlProblem
    f::Function # Dynamics function x_{t+1} = f(x_t, u_t)
    c::Function # Stage cost function c(k, x, u) where the time index k == 0 is the initial timestep
    h::Function # Terminal cost function h(x)
    W::Function # Additive process noise covariance matrix W(k) where the time index k == 0 is the initial timestep
    N::Int64 # horizon length > 0
end

# Optimal Control Problem for Sampling-based MPC methods
"""
    FiniteHorizonGenerativeOptimalControlProblem(f_stochastic, c, h, N) <: OptimalControlProblem

A finite horizon, stochastic optimal control problem where the dynamics function is stochastic and generative.

# Arguments
- `f_stochastic(x, u, rng, use_true_model=false)` -- stochastic dynamics function
  - `x` is a state vector and `u` is a control input vector.
  - The third positional argument `rng` is a random seed.
  - The fourth positional argument `use_true_model` determines whether a solver has access
    to the true stochastic dynamics and defaults to `false`.
  - The return value is the (noisy) next state `x_next`.

- `c(k, x, u)` -- stage cost function
  - `k::Int >= 0` is a time index where `k == 0` is the initial time.
  -  We assume that `c` is non-negative.
- `h(x)` -- terminal cost function
  -  We assume that `h` is non-negative.
- `N::Int64` -- final time index
  - Note that `0` is the initial time index.

# Example
```julia
import Distributions;
import LinearAlgebra;

function f_stochastic(x, u, rng, use_true_model=false)
    Σ_1 = Matrix(0.5LinearAlgebra.I, 2, 2);

    if use_true_model  # accurate GMM model
        Σ_2 = Matrix(1.0LinearAlgebra.I, 2, 2)
        d = Distributions.MixtureModel([Distributions.MvNormal(zeros(2), Σ_1),
                                        Distributions.MvNormal(ones(2), Σ_2)],
                                        [0.5, 0.5]);
    else  # inaccurate Gaussian model
        d = Distributions.MvNormal(zeros(2), Σ_1);
    end

    x_next = x + u + Distributions.rand(rng, d); # 2D single integrator dynamics
    return x_next;
end

c(k, x, u) = k/2*x'*x + k/2*u'*u  # time-dependent quadratic stage cost
N = 10;
h(x) = N/2*x'*x; # quadratic terminal cost
W(k) = Matrix(0.1LinearAlgebra.I, 2, 2);

problem = FiniteHorizonGenerativeOptimalControlProblem(f_stochastic, c, h, N);
```
"""
mutable struct FiniteHorizonGenerativeOptimalControlProblem <: OptimalControlProblem
    f_stochastic::Function # Stochastic dynamics function x_{t+1} = f_stochastic(x_t, u_t, rng, use_true_model=false)
    c::Function # Stage cost function c(k, x, u) where the time index k == 0 is the initial timestep
    h::Function # Terminal cost function h(x)
    N::Int64 # horizon length > 0
end


text:
"Solution of a partitioned ordinary differential equation."
immutable SolutionPODE{T} <: Solution{T,4}
    nd::Int
    nt::Int
    n0::Int
    t::Timeseries{T}
    x::Array{T,4}
    q::AbstractArray{T,3}
    p::AbstractArray{T,3}
    ntime::Int
    nsave::Int

    function SolutionPODE(nd, n0, ntime, nsave, Δt)
        @assert T <: Real
        @assert nd > 0
        @assert n0 > 0
        @assert nsave > 0
        @assert ntime ≥ nsave
        @assert mod(ntime, nsave) == 0

        nt = div(ntime, nsave)
        t = Timeseries{T}(nt, Δt, nsave)
        x = zeros(T, 2, nd, nt+1, n0)
        q = view(x, 1, :, :, :)
        p = view(x, 2, :, :, :)
        new(nd, nt, n0, t, x, q, p, ntime, nsave)
    end
end

function SolutionPODE{T}(equation::Union{PODE{T},SODE{T}}, Δt::T, ntime::Int, nsave::Int=1)
    s = SolutionPODE{T}(equation.d, equation.n, ntime, nsave, Δt)
    set_initial_conditions!(s, equation)
    return s
end

function set_initial_conditions!(solution::SolutionPODE, equation::Union{PODE,SODE})
    for k in 1:solution.n0
        for i in 1:solution.nd
            solution[1, i, 0, k] = equation.q₀[i,k]
            solution[2, i, 0, k] = equation.p₀[i,k]
        end
    end
    solution.t[0] = equation.t₀
    compute_timeseries!(solution.t)
end

function reset!(s::SolutionPODE)
    for k in 1:solution.n0
        for i in 1:solution.nd
            solution[1, i, 0, k] = solution[1, i, end, k]
            solution[2, i, 0, k] = solution[2, i, end, k]
        end
    end
end

Base.indices(s::SolutionPODE) = (1:2, 1:s.nd, 0:s.nt, 1:s.n0)
Base.strides(s::SolutionPODE) = (1, 2, 2*s.nd, 2*s.nd*s.nt)
# Base.linearindexing{T<:SolutionPODE}(::Type{T}) = LinearFast()

@inline function Base.getindex(s::SolutionPODE, i::Int, j::Int, k::Int, m::Int)
    @boundscheck checkbounds(s.x, i, j, k+1, m)
    @inbounds r = getindex(s.x, i, j, k+1, m)
    return r
end

@inline function Base.getindex(s::SolutionPODE, j::Int, k::Int, m::Int)
    if s.n0 == 1
        @boundscheck checkbounds(s.x, j, k, m+1, 1)
        @inbounds r = getindex(s.x, j, k, m+1, 1)
    else
        @boundscheck checkbounds(s.x, :, j, k+1, m)
        @inbounds r = getindex(s.x, :, j, k+1, m)
    end
    return r
end

@inline function Base.getindex(s::SolutionPODE, k::Int, m::Int)
    if s.n0 == 1
        @boundscheck checkbounds(s.x, :, k, m+1, 1)
        @inbounds r = getindex(s.x, :, k, m+1, 1)
    else
        @boundscheck checkbounds(s.x, :, 1:s.nd, k+1, m)
        @inbounds r = getindex(s.x, :, 1:s.nd, k+1, m)
    end
    return r
end

@inline function Base.getindex(s::SolutionPODE, m::Int)
    if s.n0 == 1
        @boundscheck checkbounds(s.x, :, 1:s.nd, m, 1)
        @inbounds r = getindex(s.x, :, 1:s.nd, m, 1)
    else
        @boundscheck checkbounds(s.x, :, 1:s.nd, :, m)
        @inbounds r = getindex(s.x, :, 1:s.nd, :, m)
    end
    return r
end

@inline function Base.setindex!(s::SolutionPODE, x, i::Int, j::Int, k::Int, m::Int)
    @boundscheck checkbounds(s.x, i, j, k+1, m)
    @inbounds setindex!(s.x, x, i, j, k+1, m)
end

@inline function Base.setindex!(s::SolutionPODE, x, j::Int, k::Int, m::Int)
    if s.n0 == 1
        @assert length(x) == 1
        @boundscheck checkbounds(s.x, j, k, m+1, 1)
        @inbounds setindex!(s.x, x, j, k, m+1, 1)
    else
        @assert length(x) == 2
        @boundscheck checkbounds(s.x, :, j, k+1, m)
        @inbounds setindex!(s.x, x, :, j, k+1, m)
    end
end

@inline function Base.setindex!(s::SolutionPODE, x, k::Int, m::Int)
    if s.n0 == 1
        @assert length(x) == 2
        @boundscheck checkbounds(s.x, :, k, m+1, 1)
        @inbounds setindex!(s.x, x, :, k, m+1, 1)
    else
        @assert size(x, 1) == 2
        @assert size(x, 2) == s.nd
        @boundscheck checkbounds(s.x, :, :, k+1, m)
        @inbounds setindex!(s.x, x, :, :, k+1, m)
    end
end

@inline function Base.setindex!(s::SolutionPODE, x, m::Int)
    if s.n0 == 1
        @assert size(x, 1) == 2
        @assert size(x, 2) == s.nd
        @boundscheck checkbounds(s.x, :, :, m+1, 1)
        @inbounds setindex!(s.x, x, :, :, m+1, 1)
    else
        @assert size(x, 1) == 2
        @assert size(x, 2) == s.nd
        @assert size(x, 3) == s.nt
        @boundscheck checkbounds(s.x, :, :, :, m)
        @inbounds setindex!(s.x, x, :, :, :, m)
    end
end


text:
@testset "Deterministic" begin
    img = generate_surface(SynthSphere(), 0.5, [0,0,1], radius = 50, background = false)
    p,q = synthetic_gradient(SynthSphere(), radius = 40)
    p,q = retrieve_surface(DeterministicCustom(), img, p, q)
    @test p[76,76] ≈ -2.3853543385909196e-5
    @test p[1,1] ≈ 0.0
    @test p[151,151] ≈ 0.0
    @test q[90,90] ≈ 0.4154339757326621
    @test q[101,36] ≈ 0.14605145430763825
    @test q[22, 136] ≈ 0.0
    @test q[42,42] ≈ -0.07890922635723528
end


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

abstract type MethodTableView; end

struct MethodLookupResult
    # Really Vector{Core.MethodMatch}, but it's easier to represent this as
    # and work with Vector{Any} on the C side.
    matches::Vector{Any}
    valid_worlds::WorldRange
    ambig::Bool
end
length(result::MethodLookupResult) = length(result.matches)
function iterate(result::MethodLookupResult, args...)
    r = iterate(result.matches, args...)
    r === nothing && return nothing
    match, state = r
    return (match::MethodMatch, state)
end
getindex(result::MethodLookupResult, idx::Int) = getindex(result.matches, idx)::MethodMatch

"""
    struct InternalMethodTable <: MethodTableView

A struct representing the state of the internal method table at a
particular world age.
"""
struct InternalMethodTable <: MethodTableView
    world::UInt
end

"""
    struct OverlayMethodTable <: MethodTableView

Overlays the internal method table such that specific queries can be redirected to an
external table, e.g., to override existing method.
"""
struct OverlayMethodTable <: MethodTableView
    world::UInt
    mt::Core.MethodTable
end

"""
    struct CachedMethodTable <: MethodTableView

Overlays another method table view with an additional local fast path cache that
can respond to repeated, identical queries faster than the original method table.
"""
struct CachedMethodTable{T} <: MethodTableView
    cache::IdDict{Any, Union{Missing, MethodLookupResult}}
    table::T
end
CachedMethodTable(table::T) where T =
    CachedMethodTable{T}(IdDict{Any, Union{Missing, MethodLookupResult}}(),
        table)

"""
    findall(sig::Type, view::MethodTableView; limit=typemax(Int))

Find all methods in the given method table `view` that are applicable to the
given signature `sig`. If no applicable methods are found, an empty result is
returned. If the number of applicable methods exceeded the specified limit,
`missing` is returned.
"""
function findall(@nospecialize(sig::Type), table::InternalMethodTable; limit::Int=typemax(Int))
    _min_val = RefValue{UInt}(typemin(UInt))
    _max_val = RefValue{UInt}(typemax(UInt))
    _ambig = RefValue{Int32}(0)
    ms = _methods_by_ftype(sig, nothing, limit, table.world, false, _min_val, _max_val, _ambig)
    if ms === false
        return missing
    end
    return MethodLookupResult(ms::Vector{Any}, WorldRange(_min_val[], _max_val[]), _ambig[] != 0)
end

function findall(@nospecialize(sig::Type), table::OverlayMethodTable; limit::Int=typemax(Int))
    _min_val = RefValue{UInt}(typemin(UInt))
    _max_val = RefValue{UInt}(typemax(UInt))
    _ambig = RefValue{Int32}(0)
    ms = _methods_by_ftype(sig, table.mt, limit, table.world, false, _min_val, _max_val, _ambig)
    if ms === false
        return missing
    elseif isempty(ms)
        # fall back to the internal method table
        _min_val[] = typemin(UInt)
        _max_val[] = typemax(UInt)
        ms = _methods_by_ftype(sig, nothing, limit, table.world, false, _min_val, _max_val, _ambig)
        if ms === false
            return missing
        end
    end
    return MethodLookupResult(ms::Vector{Any}, WorldRange(_min_val[], _max_val[]), _ambig[] != 0)
end

function findall(@nospecialize(sig::Type), table::CachedMethodTable; limit::Int=typemax(Int))
    box = Core.Box(sig)
    return get!(table.cache, sig) do
        findall(box.contents, table.table; limit=limit)
    end
end

"""
    findsup(sig::Type, view::MethodTableView)::Union{Tuple{MethodMatch, WorldRange}, Nothing}

Find the (unique) method `m` such that `sig <: m.sig`, while being more
specific than any other method with the same property. In other words, find
the method which is the least upper bound (supremum) under the specificity/subtype
relation of the queried `signature`. If `sig` is concrete, this is equivalent to
asking for the method that will be called given arguments whose types match the
given signature. This query is also used to implement `invoke`.

Such a method `m` need not exist. It is possible that no method is an
upper bound of `sig`, or it is possible that among the upper bounds, there
is no least element. In both cases `nothing` is returned.
"""
function findsup(@nospecialize(sig::Type), table::InternalMethodTable)
    min_valid = RefValue{UInt}(typemin(UInt))
    max_valid = RefValue{UInt}(typemax(UInt))
    result = ccall(:jl_gf_invoke_lookup_worlds, Any, (Any, UInt, Ptr{Csize_t}, Ptr{Csize_t}),
                   sig, table.world, min_valid, max_valid)::Union{MethodMatch, Nothing}
    result === nothing && return nothing
    (result.method, WorldRange(min_valid[], max_valid[]))
end

# This query is not cached
findsup(@nospecialize(sig::Type), table::CachedMethodTable) = findsup(sig, table.table)

isoverlayed(::MethodTableView)     = error("unsatisfied MethodTableView interface")
isoverlayed(::InternalMethodTable) = false
isoverlayed(::OverlayMethodTable)  = true
isoverlayed(mt::CachedMethodTable) = isoverlayed(mt.table)

max_stars_repo_name: fkastner/julia

max_issues_repo_name: fkastner/julia

max_forks_repo_name: fkastner/julia


text:
J = 400
B = 10.0
θ = 0.0024
M = 0.0025
Gt = 2.0
L = 1000.0
n0 = 0.033
nt = 0.0133
Q = 30
ηu = 0.9
ηd = 0.95
f = 2
Hu = 3
Hd = 3
Vmin = 1.0
Vmax = 2.8
N_turbine = 5
Pmax = 5000

upstream = false
downstream = false
enumerate = false
solver = bonmin
parameters = Dict()

max_stars_repo_name: Thiago-NovaesB/alpheus

max_issues_repo_name: Thiago-NovaesB/alpheus

max_forks_repo_name: Thiago-NovaesB/alpheus