text:
abstract type AbstractVolumeIntegral end

get_element_variables!(element_variables, u, mesh, equations,
                       volume_integral::AbstractVolumeIntegral, dg, cache) = nothing

"""
    VolumeIntegralWeakForm

The classical weak form volume integral type for DG methods as explained in standard
textbooks such as
- Kopriva (2009)
  Implementing Spectral Methods for Partial Differential Equations:
  Algorithms for Scientists and Engineers
  [doi: 10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
"""
struct VolumeIntegralWeakForm <: AbstractVolumeIntegral end

create_cache(mesh, equations, ::VolumeIntegralWeakForm, dg, uEltype) = NamedTuple()

" BLZ: Pseudo Strong Form for Implementation of Flux Differences"
struct VolumeIntegralPseudoStrongForm <: AbstractVolumeIntegral end

create_cache(mesh, equations, ::VolumeIntegralPseudoStrongForm, dg, uEltype) = NamedTuple()

"""
    VolumeIntegralFluxDifferencing

Volume integral type for DG methods based on SBP operators and flux differencing using
symmetric two-point volume fluxes. Based upon the theory developed by
- LeFloch, Mercier, Rohde (2002)
  Fully Discrete, Entropy Conservative Schemes of Arbitrary Order
  [doi: 10.1137/S003614290240069X](https://doi.org/10.1137/S003614290240069X)
- Fisher, Carpenter (2013)
  High-order entropy stable finite difference schemes for nonlinear
  conservation laws: Finite domains
  [doi: 10.1016/j.jcp.2013.06.014](https://doi.org/10.1016/j.jcp.2013.06.014)
- Ranocha (2017)
  Comparison of Some Entropy Conservative Numerical Fluxes for the Euler Equations
  [arXiv: 1701.02264](https://arxiv.org/abs/1701.02264)
  [doi: 10.1007/s10915-017-0618-1](https://doi.org/10.1007/s10915-017-0618-1)
- Chen, Shu (2017)
  Entropy stable high order discontinuous Galerkin methods with suitable
  quadrature rules for hyperbolic conservation laws
  [doi: 10.1016/j.jcp.2017.05.025](https://doi.org/10.1016/j.jcp.2017.05.025)
"""
struct VolumeIntegralFluxDifferencing{VolumeFlux} <: AbstractVolumeIntegral
  volume_flux::VolumeFlux
end

function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxDifferencing)
  @nospecialize integral # reduce precompilation time

  if get(io, :compact, false)
    show(io, integral)
  else
    setup = [
            "volume flux" => integral.volume_flux
            ]
    summary_box(io, "VolumeIntegralFluxDifferencing", setup)
  end
end

"""
    VolumeIntegralShockCapturingHG

Shock-capturing volume integral type for DG methods proposed by
- Hennemann, Gassner (2020)
  "A provably entropy stable subcell shock capturing approach for high order split form DG"
  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
"""
struct VolumeIntegralShockCapturingHG{VolumeFluxDG, VolumeFluxFV, Indicator} <: AbstractVolumeIntegral
  volume_flux_dg::VolumeFluxDG # symmetric, e.g. split-form or entropy-conservative
  volume_flux_fv::VolumeFluxFV # non-symmetric in general, e.g. entropy-dissipative
  indicator::Indicator
end

function VolumeIntegralShockCapturingHG(indicator; volume_flux_dg=flux_central,
                                                   volume_flux_fv=flux_lax_friedrichs)
  VolumeIntegralShockCapturingHG{typeof(volume_flux_dg), typeof(volume_flux_fv), typeof(indicator)}(
    volume_flux_dg, volume_flux_fv, indicator)
end

function Base.show(io::IO, mime::MIME"text/plain", integral::VolumeIntegralShockCapturingHG)
  @nospecialize integral # reduce precompilation time

  if get(io, :compact, false)
    show(io, integral)
  else
    summary_header(io, "VolumeIntegralShockCapturingHG")
    summary_line(io, "volume flux DG", integral.volume_flux_dg)
    summary_line(io, "volume flux FV", integral.volume_flux_fv)
    summary_line(io, "indicator", integral.indicator |> typeof |> nameof)
    show(increment_indent(io), mime, integral.indicator)
    summary_footer(io)
  end
end

"""
    VolumeIntegralPureLGLFiniteVolume

A volume integral that only uses the subcell finite volume scheme from the paper
- Hennemann, Gassner (2020)
  "A provably entropy stable subcell shock capturing approach for high order split form DG"
  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
This gives a formally O(1)-accurate finite volume scheme on an LGL-type subcell mesh (LGL = Legendre-Gauss-Lobatto).
!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
struct VolumeIntegralPureLGLFiniteVolume{VolumeFluxFV} <: AbstractVolumeIntegral
  volume_flux_fv::VolumeFluxFV # non-symmetric in general, e.g. entropy-dissipative
end
# TODO: Figure out if this can also be used for Gauss nodes, not just LGL, and adjust the name accordingly

function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralPureLGLFiniteVolume)
  @nospecialize integral # reduce precompilation time

  if get(io, :compact, false)
    show(io, integral)
  else
    setup = [
            "FV flux" => integral.volume_flux_fv
            ]
    summary_box(io, "VolumeIntegralPureLGLFiniteVolume", setup)
  end
end

function get_element_variables!(element_variables, u, mesh, equations,
                                volume_integral::VolumeIntegralShockCapturingHG, dg, cache)
  # call the indicator to get up-to-date values for IO
  volume_integral.indicator(u, equations, dg, cache)
  get_element_variables!(element_variables, volume_integral.indicator, volume_integral)
end

"""
    DG(basis, basis, surface_flux, volume_integral)

Create a discontinuous Galerkin method.
If [`basis isa LobattoLegendreBasis`](@ref LobattoLegendreBasis),
this creates a [`DGSEM`](@ref).
"""
struct DG{RealT, Basis<:AbstractBasisSBP{RealT}, Mortar, SurfaceFlux, VolumeIntegral}
  basis::Basis
  mortar::Mortar
  surface_flux::SurfaceFlux
  volume_integral::VolumeIntegral
end

function Base.show(io::IO, dg::DG)
  @nospecialize dg # reduce precompilation time

  print(io, "DG{", real(dg), "}(")
  print(io,       dg.basis)
  print(io, ", ", dg.mortar)
  print(io, ", ", dg.surface_flux)
  print(io, ", ", dg.volume_integral)
  print(io, ")")
end

function Base.show(io::IO, mime::MIME"text/plain", dg::DG)
  @nospecialize dg # reduce precompilation time

  if get(io, :compact, false)
    show(io, dg)
  else
    summary_header(io, "DG{" * string(real(dg)) * "}")
    summary_line(io, "polynomial degree", polydeg(dg))
    summary_line(io, "basis", dg.basis)
    summary_line(io, "mortar", dg.mortar)
    summary_line(io, "surface flux", dg.surface_flux)
    summary_line(io, "volume integral", dg.volume_integral |> typeof |> nameof)
    if !(dg.volume_integral isa VolumeIntegralWeakForm)
      show(increment_indent(io), mime, dg.volume_integral)
    end
    summary_footer(io)
  end
end

@inline Base.real(dg::DG{RealT}) where {RealT} = RealT

# TODO: Taal refactor, use case?
# Deprecate in favor of nnodes or order_of_accuracy?
@inline polydeg(dg::DG) = polydeg(dg.basis)

@inline ndofs(mesh::TreeMesh, dg::DG, cache) = nelements(cache.elements) * nnodes(dg)^ndims(mesh)

function get_element_variables!(element_variables, u, mesh, equations, dg::DG, cache)
  get_element_variables!(element_variables, u, mesh, equations, dg.volume_integral, dg, cache)
end

# TODO: Taal performance, 1:nnodes(dg) vs. Base.OneTo(nnodes(dg)) vs. SOneTo(nnodes(dg))
@inline eachnode(dg::DG)             = Base.OneTo(nnodes(dg))
@inline eachelement(dg::DG, cache)   = Base.OneTo(nelements(dg, cache))
@inline eachinterface(dg::DG, cache) = Base.OneTo(ninterfaces(dg, cache))
@inline eachboundary(dg::DG, cache)  = Base.OneTo(nboundaries(dg, cache))
@inline eachmortar(dg::DG, cache)    = Base.OneTo(nmortars(dg, cache))
@inline eachmpiinterface(dg::DG, cache) = Base.OneTo(nmpiinterfaces(dg, cache))

@inline nnodes(dg::DG)             = nnodes(dg.basis)
@inline nelements(dg::DG, cache)   = nelements(cache.elements)
@inline nelementsglobal(dg::DG, cache) = mpi_isparallel() ? cache.mpi_cache.n_elements_global : nelements(dg, cache)
@inline ninterfaces(dg::DG, cache) = ninterfaces(cache.interfaces)
@inline nboundaries(dg::DG, cache) = nboundaries(cache.boundaries)
@inline nmortars(dg::DG, cache)    = nmortars(cache.mortars)
@inline nmpiinterfaces(dg::DG, cache) = nmpiinterfaces(cache.mpi_interfaces)

@inline function get_node_coords(x, equations, solver::DG, indices...)
  SVector(ntuple(idx -> x[idx, indices...], ndims(equations)))
end

@inline function get_node_vars(u, equations, solver::DG, indices...)
  # There is a cut-off at `n == 10` inside of the method
  # `ntuple(f::F, n::Integer) where F` in Base at ntuple.jl:17
  # in Julia `v1.5`, leading to type instabilities if
  # more than ten variables are used. That's why we use
  # `Val(...)` below.
  SVector(ntuple(v -> u[v, indices...], Val(nvariables(equations))))
end

@inline function get_surface_node_vars(u, equations, solver::DG, indices...)
  # There is a cut-off at `n == 10` inside of the method
  # `ntuple(f::F, n::Integer) where F` in Base at ntuple.jl:17
  # in Julia `v1.5`, leading to type instabilities if
  # more than ten variables are used. That's why we use
  # `Val(...)` below.
  u_ll = SVector(ntuple(v -> u[1, v, indices...], Val(nvariables(equations))))
  u_rr = SVector(ntuple(v -> u[2, v, indices...], Val(nvariables(equations))))
  return u_ll, u_rr
end

@inline function set_node_vars!(u, u_node, equations, solver::DG, indices...)
  for v in eachvariable(equations)
    u[v, indices...] = u_node[v]
  end
  return nothing
end

@inline function add_to_node_vars!(u, u_node, equations, solver::DG, indices...)
  for v in eachvariable(equations)
    u[v, indices...] += u_node[v]
  end
  return nothing
end

# Used for analyze_solution
SolutionAnalyzer(dg::DG; kwargs...) = SolutionAnalyzer(dg.basis; kwargs...)

AdaptorAMR(mesh, dg::DG) = AdaptorL2(dg.basis)

# Include utilities
include("interpolation.jl")
include("l2projection.jl")
include("basis_lobatto_legendre.jl")

"""
    DGSEM([RealT=Float64,] polydeg::Integer,
          surface_flux=flux_central,
          volume_integral::AbstractVolumeIntegral=VolumeIntegralWeakForm(),
          mortar=MortarL2(basis))

Create a discontinuous Galerkin spectral element method (DGSEM) using a
[`LobattoLegendreBasis`](@ref) with polynomials of degree `polydeg`.
"""
const DGSEM = DG{RealT, Basis, Mortar, SurfaceFlux, VolumeIntegral} where {RealT<:Real, Basis<:LobattoLegendreBasis{RealT}, Mortar, SurfaceFlux, VolumeIntegral}

function DGSEM(basis::LobattoLegendreBasis,
               surface_flux=flux_central,
               volume_integral::AbstractVolumeIntegral=VolumeIntegralWeakForm(),
               mortar=MortarL2(basis))

  return DG{real(basis), typeof(basis), typeof(mortar), typeof(surface_flux), typeof(volume_integral)}(
    basis, mortar, surface_flux, volume_integral)
end

function DGSEM(RealT, polydeg::Integer,
               surface_flux=flux_central,
               volume_integral::AbstractVolumeIntegral=VolumeIntegralWeakForm(),
               mortar=MortarL2(LobattoLegendreBasis(RealT, polydeg)))
  basis = LobattoLegendreBasis(RealT, polydeg)

  return DGSEM(basis, surface_flux, volume_integral, mortar)
end

DGSEM(polydeg, surface_flux=flux_central, volume_integral::AbstractVolumeIntegral=VolumeIntegralWeakForm()) = DGSEM(Float64, polydeg, surface_flux, volume_integral)

# The constructor using only keyword arguments is convenient for elixirs since
# it allows to modify the polynomial degree and other parameters via
# `trixi_include`.
function DGSEM(; RealT=Float64,
                 polydeg::Integer,
                 surface_flux=flux_central,
                 volume_integral=VolumeIntegralWeakForm())
  basis = LobattoLegendreBasis(RealT, polydeg)
  return DGSEM(basis, surface_flux, volume_integral)
end

"""
    pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache)

Given blending factors `alpha` and the solver `dg`, fill
`element_ids_dg` with the IDs of elements using a pure DG scheme and
`element_ids_dgfv` with the IDs of elements using a blended DG-FV scheme.
"""
function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg::DG, cache)
  empty!(element_ids_dg)
  empty!(element_ids_dgfv)

  for element in eachelement(dg, cache)
    # Clip blending factor for values close to zero (-> pure DG)
    dg_only = isapprox(alpha[element], 0, atol=1e-12)
    if dg_only
      push!(element_ids_dg, element)
    else
      push!(element_ids_dgfv, element)
    end
  end

  return nothing
end

function volume_jacobian(element, mesh::TreeMesh, cache)
  return inv(cache.elements.inverse_jacobian[element])^ndims(mesh)
end

# indicators used for shock-capturing and AMR
include("indicators.jl")
include("indicators_1d.jl")
include("indicators_2d.jl")
include("indicators_3d.jl")

# 1D DG implementation
include("containers_1d.jl")
include("dg_1d.jl")

# 2D DG implementation
include("containers_2d.jl")
include("dg_2d.jl")
include("dg_2d_parallel.jl")

# 3D DG implementation
include("containers_3d.jl")
include("dg_3d.jl")


text:
using Documenter
using FFTLog
using FFTLogDocs
using Plots
using PlotThemes

ENV["GKSwstype"] = "100"

push!(LOAD_PATH,"../src/")

makedocs(
    modules = [FFTLog, FFTLog],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
    sidebar_sitename=true),
    sitename = "FFTLog.jl",
    authors  = "Marco Bonici",
    pages = [
        "Home" => "index.md",
        "The algorithm" => "algorithm.md"
    ]
)

deploydocs(
    repo = "github.com/marcobonici/FFTLogDocs.jl.git",
    devbranch = "main"
)

ext: jl


text:
module CheckTests
using Base.Test
using FactCheck
using ValidatedNumerics
using IntervalStatistic
using Distributions

srand(10)

in_3_sigma_interval(mu, sigma) = (i) -> (mu in i) && (i in mu + 3*@interval(-sigma, sigma))

facts("test goodness-fit of standard (0, 1) distribution with true distribution") do
    d = Normal()
    length = 100
    values = rand(d, length)
    mu, sigma = params(d)
    average = mean(values)

    context("by chi_square check with Sturges k") do
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.SturgesChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Scott k") do
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ScottChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Taylor k") do
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.TaylorChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with FreedmanDiaconis k") do
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.FreedmanDiaconisChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Doane k") do
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.DoaneChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Wichard k") do
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.WichardChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square modified for normal dist check with 10 bins") do
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ChiSquareNormalCheck(0.05, 10, mu, sigma)
        )
        println(result)
        @fact result --> true
    end

    context("by chi_square modified for normal dist check with 11 bins") do
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ChiSquareNormalCheck(0.05, 11, mu, sigma)
        )
        println(result)
        @pending @fact result --> true
    end
end

facts("test goodness-fit of standard (0, 1) distribution") do
    d = Normal()
    length = 100
    values = rand(d, length)
    mu, sigma = params(d)
    average = mean(values)

    context("by chi_square check with Sturges k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.SturgesChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Scott k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ScottChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Taylor k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.TaylorChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with FreedmanDiaconis k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.FreedmanDiaconisChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Doane k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.DoaneChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Wichard k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.WichardChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square modified for normal dist check with 10 bins") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ChiSquareNormalCheck(0.05, 10, mu, sigma)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square modified for normal dist check with 11 bins") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ChiSquareNormalCheck(0.05, 11, mu, sigma)
        )
        println(result)
        @pending @fact result --> true
    end
end

facts("test goodness-fit of standard (0, 1) distribution with samlping=500") do
    d = Normal()
    length = 500
    values = rand(d, length)
    mu, sigma = params(d)
    average = mean(values)

    context("by chi_square check with k for Large n") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.LargeNChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end
end

facts("test goodness-fit of normal mu=3, sigma=0.1 distribution") do
    d = Normal(3, 0.1)
    length = 100
    values = rand(d, length)
    mu, sigma = params(d)
    average = mean(values)

    context("by chi_square check with Dahiya k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.SturgesChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Scott k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ScottChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Taylor k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.TaylorChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with FreedmanDiaconis k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.FreedmanDiaconisChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Doane k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.DoaneChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Wichard k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.WichardChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square modified for normal dist check with 10 bins") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ChiSquareNormalCheck(0.05, 10, mu, sigma)
        )
        println(result)
        @fact result --> true
    end

    context("by chi_square modified for normal dist check with 11 bins") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ChiSquareNormalCheck(0.05, 11, mu, sigma)
        )
        println(result)
        @fact result --> true
    end
end

facts("test goodness-fit of normal mu=3, sigma=0.1 distribution with true distribution") do
    d = Normal(3, 0.1)
    length = 100
    values = rand(d, length)
    mu, sigma = params(d)
    average = mean(values)

    context("by chi_square check with Dahiya k") do
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.SturgesChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Scott k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ScottChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Taylor k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.TaylorChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with FreedmanDiaconis k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.FreedmanDiaconisChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Doane k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.DoaneChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square check with Wichard k") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.WichardChiSquareCheck(0.05, d)
        )
        println(result)
        @pending @fact result --> true
    end

    context("by chi_square modified for normal dist check with 10 bins") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ChiSquareNormalCheck(0.05, 10, mu, sigma)
        )
        println(result)
        @fact result --> true
    end

    context("by chi_square modified for normal dist check with 11 bins") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.ChiSquareNormalCheck(0.05, 11, mu, sigma)
        )
        println(result)
        @fact result --> true
    end
end

facts("test goodness-fit of normal mu=3, sigma=0.1 distribution with sampling=500") do
    d = Normal(3, 0.1)
    length = 500
    values = rand(d, length)
    mu, sigma = params(d)
    average = mean(values)

    context("by chi_square check with k for Large n") do
        mu = mean(values)
        sigma = sqrt(var(values))
        result = IntervalStatistic.isDistribution(
            values,
            IntervalStatistic.Check.LargeNChiSquareCheck(0.05, Normal(mu, sigma))
        )
        println(result)
        @pending @fact result --> true
    end
end

end


text:
"""
Four PIFA antennas fitting inside a Bott's dot
"""
module PIFAProblem

export PIFAAnt, create_grammar, get_grammar, get_fitness,
    OutOfBoundsException, nec_run
export expr_fitness

using ExprSearch, DerivationTrees
import ExprSearch: ExprProblem, get_fitness, get_grammar
using RLESUtils, Interpreter, LogSystems, MathUtils
using NECPP

#all measurements in meters

#design bounding box
const MAXLENGTH = 0.12
const MAXWIDTH = 0.11
const MAXHEIGHT = 0.02 #fixed height

const WIRE_RADIUS = 0.0002 #meters
const FREQ_MHZ = 915.0 #MHz
const LAMBDA = 299.79 / FREQ_MHZ
const SEG_RATIO = 27 / (LAMBDA / 2)
const EXTYPE = 0 #0=voltage source
const ROUND_DIGITS = 5

const HTOL = 0.001
const WXTOL = 0.04 #buffer space around wx
const TOL = 0.01 #min dist between two wires
const WIRE_LEN_TOL = 0.2 #antenna length tolerance from lambda/2
const WIRE_LEN_MIN = (1.0 - WIRE_LEN_TOL) * LAMBDA / 2
const WIRE_LEN_MAX = (1.0 + WIRE_LEN_TOL) * LAMBDA / 2

const DIVCONST = 100.0

immutable OutOfBoundsException <: Exception end

type PIFAAnt <: ExprProblem
    departure_angle::Float64 #degrees rel horizontal
    grammar::Grammar
    Z0::Complex
end

function PIFAAnt(departure_angle::Float64=10.0,
    Z0::Complex=Complex(50.0, 0.0))
    grammar = create_grammar()
    PIFAAnt(departure_angle, grammar, Z0)
end

function create_grammar()
    @grammar grammar begin
        start = ant
        ant[make_block] = assign_drv + dir_el
        assign_drv = Expr(:(=), :lxwx, drv_el)
        drv_el = Expr(:call, :drv_wire, :nec, rn, rn, rn, rn)
        dir_el = dir_case2a | dir_case2b

        dir_case2a = Expr(:call, :dir_wire_2a, :nec, :lxwx, rn, rn, rn, rn, rn)
        dir_case2b = Expr(:call, :dir_wire_2b, :nec, :lxwx, rn, rn, rn, rn, rn)
        rn[divconst] = 0:100 #fraction of some length
    end
    grammar
end

"""
n_digits after the decimal point
"""
function rand_and_round(r::Float64, xmin::Float64, xmax::Float64, n_digits::Int64)
    x = xmin + r * (xmax - xmin)
    round(x, n_digits)
end

make_block(lst::Array) = Expr(:block, lst...)
divconst(x) = x / DIVCONST

###
#functions for grammar

function drv_wire(nec::NecContext, rn_tl::Float64, rn_w::Float64, rn_height1::Float64,
    rn_height2::Float64)
    tlx = rand_and_round(rn_tl, WIRE_LEN_MIN, WIRE_LEN_MAX, ROUND_DIGITS)
    wxmin = WXTOL
    wxmax = MAXWIDTH
    wx = rand_and_round(rn_w, wxmin, wxmax, ROUND_DIGITS)
    nec.userargs[:drv_n_segments] = round_nearest_odd(wx * SEG_RATIO) #this calculation is duplicated in wire
    lx = (tlx - wx) / 2
    height1 = rand_and_round(rn_height1, HTOL, MAXHEIGHT, ROUND_DIGITS)
    height2 = rand_and_round(rn_height2, HTOL, MAXHEIGHT, ROUND_DIGITS)
    sym_bent_wire(nec, 1, 0.0, height1, height2, WIRE_RADIUS, tlx, wx, true)
    (lx, wx)
end

function dir_wire_2a(nec::NecContext, lxwx::Tuple, rn_tl::Float64,
    rn_w::Float64, rn_xpos::Float64, rn_height1::Float64, rn_height2::Float64)

    lx, wx = lxwx
    tld = rand_and_round(rn_tl, WIRE_LEN_MIN, WIRE_LEN_MAX, ROUND_DIGITS)

    wdmin = TOL
    wdmax = wx - TOL
    wd = rand_and_round(rn_w, wdmin, wdmax, ROUND_DIGITS)
    ld = (tld - wd) / 2

    xposmin = ld + TOL
    xposmax = MAXLENGTH
    xpos = rand_and_round(rn_xpos, xposmin, xposmax, ROUND_DIGITS)

    height1 = rand_and_round(rn_height1, HTOL, MAXHEIGHT, ROUND_DIGITS)
    height2 = rand_and_round(rn_height2, HTOL, MAXHEIGHT, ROUND_DIGITS)
    sym_bent_wire(nec, 4, xpos, height1, height2, WIRE_RADIUS, tld, wd, false)
    nothing
end
function dir_wire_2b(nec::NecContext, lxwx::Tuple, rn_tl::Float64,
    rn_w::Float64, rn_xpos::Float64, rn_height1::Float64, rn_height2::Float64)

    lx, wx = lxwx
    tld = rand_and_round(rn_tl, WIRE_LEN_MIN, WIRE_LEN_MAX, ROUND_DIGITS)

    wdmin = TOL
    wdmax = wx - TOL
    wd = rand_and_round(rn_w, wdmin, wdmax, ROUND_DIGITS)
    ld = (tld - wd) / 2

    xposmin = TOL
    xposmax = MAXLENGTH - ld
    xpos = rand_and_round(rn_xpos, xposmin, xposmax, ROUND_DIGITS)

    height1 = rand_and_round(rn_height1, HTOL, MAXHEIGHT, ROUND_DIGITS)
    height2 = rand_and_round(rn_height2, HTOL, MAXHEIGHT, ROUND_DIGITS)
    sym_bent_wire(nec, 4, xpos, height1, height2, WIRE_RADIUS, tld, wd, true)
    nothing
end
function sym_bent_wire(nec::NecContext, tag_id_base::Int64, xpos::Float64,
    height1::Float64, height2::Float64, wire_radius::Float64, total_len::Float64,
    width::Float64, up::Bool)

    halfwidth = width / 2
    rdel = 1.0 #no tapering
    rrad = 1.0 #no tapering
    rad = wire_radius

    tag_id = tag_id_base
    segment_count = round_nearest_odd(width * SEG_RATIO)
    xw1 = xpos
    yw1 = +halfwidth
    xw2 = xpos
    yw2 = -halfwidth
    check_boundaries([xw1, xw2], 0.0, MAXLENGTH)
    check_boundaries([yw1, yw2], -MAXWIDTH/2, MAXWIDTH/2)
    handle_nec(nec_wire(nec, tag_id, segment_count, xw1, yw1, height1,
        xw2, yw2, height1, rad, rdel, rrad)) #base

    for y in [+halfwidth, -halfwidth]
        tag_id += 1
        xw1 = xpos
        yw1 = y
        xlen = (total_len - width)/2
        xw2 = up ? xpos + xlen : xpos - xlen
        xw2 = round(xw2, ROUND_DIGITS)
        yw2 = y
        segment_count = round_nearest_odd(xlen * SEG_RATIO)
        check_boundaries([xw1, xw2], 0.0, MAXLENGTH)
        check_boundaries([yw1, yw2], -MAXWIDTH/2, MAXWIDTH/2)
        handle_nec(nec_wire(nec, tag_id, segment_count, xw1, yw1, height1,
            xw2, yw2, height2, rad, rdel, rrad))
    end
end

function check_boundaries(xs::Vector{Float64}, xmin::Float64, xmax::Float64)
    for x in xs
        if x < xmin || x > xmax
            println("$x is outisde of ($xmin, $xmax)")
            throw(OutOfBoundsException())
        end
    end
end

const SYMTABLE = SymbolTable(
    :make_block => make_block,
    :drv_wire => drv_wire,
    :dir_wire_2a => dir_wire_2a,
    :dir_wire_2b => dir_wire_2b,
    :divconst => divconst
    )

function eval_expr(problem::PIFAAnt, expr, nec)
    SYMTABLE[:nec] = nec
    interpret(SYMTABLE, expr)
end

ExprSearch.get_grammar(problem::PIFAAnt) = problem.grammar

function ExprSearch.get_fitness(problem::PIFAAnt, derivtree::DerivationTree,
    userargs::SymbolTable)

    expr = get_expr(derivtree)

    #####################################
    # debug
    #f = open("antenna.nec","w")
    #fitness = nec_run(problem, expr, f) do nec
        #gain = nec_gain(nec, 0, 0, 0)
        #imp_re = nec_impedance_real(nec, 0)
        #imp_im = nec_impedance_imag(nec, 0)
        #Z = Complex(imp_re, imp_im)
        #@show gain
        #@show Z
        #@show vswr(Z, problem.Z0)
        #-gain + vswr_metric(vswr(Z, problem.Z0))
    #end
    #close(f)
    #####################################
    #deploy
    fitness = expr_fitness(problem, expr)
    ####################################

    fitness
end

function expr_fitness(problem::PIFAAnt, expr)
    fitness = nec_run(problem, expr) do nec
        gain = nec_gain(nec, 0, 0, 0)
        imp_re = nec_impedance_real(nec, 0)
        imp_im = nec_impedance_imag(nec, 0)
        Z = Complex(imp_re, imp_im)
        -gain + vswr_metric(vswr(Z, problem.Z0))
    end
    fitness
end

#vswr_metric(x::Float64) = 10.0 / (1+e^-(2x-5))
function vswr_metric(x::Float64)
    #passes through (0,0), (1.5,0.5), (3,100)
    if 0.0 <= x <= 1.5
        return x / 3
    elseif 1.5 < x <= 3.0
        return 66.3333x - 99.0
    elseif 3.0 < x
        return 100.0
    else
        ##println("vswr cannot be negative! $x")
        #throw(DomainError())
        return 100.0
    end
end

nec_run(f::Function, problem::PIFAAnt, expr, io::IO) = nec_run(f, problem, expr, Nullable{IO}(io))
function nec_run(f::Function, problem::PIFAAnt, expr, io::Nullable{IO}=Nullable{IO}())
    if !isnull(io)
        logsys = NECPP.logsystem()
        send_to!(get(io), logsys, "nec_cards", first)
        nec = nec_create(logsys)
    else
        nec = nec_create()
    end
    eval_expr(problem, expr, nec)
    handle_nec(nec_geometry_complete(nec, 1))

    extype = EXTYPE
    i2 = 1 #source seg
    i3 = floor(Int64, nec.userargs[:drv_n_segments]/2) + 1 #middle seg of odd number of segs
    i4 = 0
    tmp1 = 1.0 #excitation volts?
    tmp2 = tmp3 = tmp4 = tmp5 = tmp6 = 0.0
    handle_nec(nec_ex_card(nec, extype, i2, i3, i4, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6)) #feedpoint

    handle_nec(nec_fr_card(nec, 0, 1, FREQ_MHZ, 0.0))
    handle_nec(nec_gn_card(nec, 1, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
    handle_nec(nec_rp_card(nec, 0, 1, 1, 1, 5, 0, 0, (90.0-problem.departure_angle), 0.0,
        0.0, 0.0, 0.0, 0.0))

    result = f(nec)
    nec_delete(nec)

    result
end

end #module


text:
function getdevice()
    ref = Ref{Cint}()
    @apicall :cuCtxGetDevice (Ptr{Cint},) ref
    Int(ref[])
end

function setdevice(dev::Int)
    #cap = capability(dev)
    #mem = round(Int, totalmem(dev) / (1024^2))
    #@info "device[$dev]: $(devicename(dev)), capability $(cap[1]).$(cap[2]), totalmem = $(mem) MB"
    getdevice() == dev && return
    isassigned(CONTEXTS,dev) || (CONTEXTS[dev+1] = CuContext(dev))
    setcontext(CONTEXTS[dev+1])
    dev
end
function setdevice(f::Function, dev::Int)
    _dev = getdevice()
    dev == _dev || setdevice(dev)
    f()
    dev == _dev || setdevice(_dev)
end

function ndevices()
    ref = Ref{Cint}()
    @apicall :cuDeviceGetCount (Ptr{Cint},) ref
    Int(ref[])
end

function devicename(dev::Int)
    buflen = 256
    buf = Vector{Cchar}(buflen)
    @apicall :cuDeviceGetName (Ptr{Cchar},Cint,Cint) buf buflen dev
    buf[end] = 0
    unsafe_string(pointer(buf))
end

function totalmem(dev::Int)
    ref = Ref{Csize_t}()
    @apicall :cuDeviceTotalMem (Ptr{Csize_t},Cint) ref dev
    Int(ref[])
end

@enum(CUdevice_attribute, MAX_THREADS_PER_BLOCK = Cint(1),
                          MAX_BLOCK_DIM_X,
                          MAX_BLOCK_DIM_Y,
                          MAX_BLOCK_DIM_Z,
                          MAX_GRID_DIM_X,
                          MAX_GRID_DIM_Y,
                          MAX_GRID_DIM_Z,
                          MAX_SHARED_MEMORY_PER_BLOCK,
                          TOTAL_CONSTANT_MEMORY,
                          WARP_SIZE,
                          MAX_PITCH,
                          MAX_REGISTERS_PER_BLOCK,
                          CLOCK_RATE,
                          TEXTURE_ALIGNMENT,
                          GPU_OVERLAP,
                          MULTIPROCESSOR_COUNT,
                          KERNEL_EXEC_TIMEOUT,
                          INTEGRATED,
                          CAN_MAP_HOST_MEMORY,
                          COMPUTE_MODE,
                          MAXIMUM_TEXTURE1D_WIDTH,
                          MAXIMUM_TEXTURE2D_WIDTH,
                          MAXIMUM_TEXTURE2D_HEIGHT,
                          MAXIMUM_TEXTURE3D_WIDTH,
                          MAXIMUM_TEXTURE3D_HEIGHT,
                          MAXIMUM_TEXTURE3D_DEPTH,
                          MAXIMUM_TEXTURE2D_LAYERED_WIDTH,
                          MAXIMUM_TEXTURE2D_LAYERED_HEIGHT,
                          MAXIMUM_TEXTURE2D_LAYERED_LAYERS,
                          SURFACE_ALIGNMENT,
                          CONCURRENT_KERNELS,
                          ECC_ENABLED,
                          PCI_BUS_ID,
                          PCI_DEVICE_ID,
                          TCC_DRIVER,
                          MEMORY_CLOCK_RATE,
                          GLOBAL_MEMORY_BUS_WIDTH,
                          L2_CACHE_SIZE,
                          MAX_THREADS_PER_MULTIPROCESSOR,
                          ASYNC_ENGINE_COUNT,
                          UNIFIED_ADDRESSING,
                          MAXIMUM_TEXTURE1D_LAYERED_WIDTH,
                          MAXIMUM_TEXTURE1D_LAYERED_LAYERS,
                          CAN_TEX2D_GATHER,
                          MAXIMUM_TEXTURE2D_GATHER_WIDTH,
                          MAXIMUM_TEXTURE2D_GATHER_HEIGHT,
                          MAXIMUM_TEXTURE3D_WIDTH_ALTERNATE,
                          MAXIMUM_TEXTURE3D_HEIGHT_ALTERNATE,
                          MAXIMUM_TEXTURE3D_DEPTH_ALTERNATE,
                          PCI_DOMAIN_ID,
                          TEXTURE_PITCH_ALIGNMENT,
                          MAXIMUM_TEXTURECUBEMAP_WIDTH,
                          MAXIMUM_TEXTURECUBEMAP_LAYERED_WIDTH,
                          MAXIMUM_TEXTURECUBEMAP_LAYERED_LAYERS,
                          MAXIMUM_SURFACE1D_WIDTH,
                          MAXIMUM_SURFACE2D_WIDTH,
                          MAXIMUM_SURFACE2D_HEIGHT,
                          MAXIMUM_SURFACE3D_WIDTH,
                          MAXIMUM_SURFACE3D_HEIGHT,
                          MAXIMUM_SURFACE3D_DEPTH,
                          MAXIMUM_SURFACE1D_LAYERED_WIDTH,
                          MAXIMUM_SURFACE1D_LAYERED_LAYERS,
                          MAXIMUM_SURFACE2D_LAYERED_WIDTH,
                          MAXIMUM_SURFACE2D_LAYERED_HEIGHT,
                          MAXIMUM_SURFACE2D_LAYERED_LAYERS,
                          MAXIMUM_SURFACECUBEMAP_WIDTH,
                          MAXIMUM_SURFACECUBEMAP_LAYERED_WIDTH,
                          MAXIMUM_SURFACECUBEMAP_LAYERED_LAYERS,
                          MAXIMUM_TEXTURE1D_LINEAR_WIDTH,
                          MAXIMUM_TEXTURE2D_LINEAR_WIDTH,
                          MAXIMUM_TEXTURE2D_LINEAR_HEIGHT,
                          MAXIMUM_TEXTURE2D_LINEAR_PITCH,
                          MAXIMUM_TEXTURE2D_MIPMAPPED_WIDTH,
                          MAXIMUM_TEXTURE2D_MIPMAPPED_HEIGHT,
                          COMPUTE_CAPABILITY_MAJOR,
                          COMPUTE_CAPABILITY_MINOR,
                          MAXIMUM_TEXTURE1D_MIPMAPPED_WIDTH,
                          STREAM_PRIORITIES_SUPPORTED,
                          GLOBAL_L1_CACHE_SUPPORTED,
                          LOCAL_L1_CACHE_SUPPORTED,
                          MAX_SHARED_MEMORY_PER_MULTIPROCESSOR,
                          MAX_REGISTERS_PER_MULTIPROCESSOR,
                          MANAGED_MEMORY,
                          MULTI_GPU_BOARD,
                          MULTI_GPU_BOARD_GROUP_ID)
@assert Cint(MULTI_GPU_BOARD_GROUP_ID) == 85

function attribute(dev::Int, code::CUdevice_attribute)
    ref = Ref{Cint}()
    @apicall :cuDeviceGetAttribute (Ptr{Cint},Cint,Cint) ref code dev
    Int(ref[])
end

function capability(dev::Int)
    major = attribute(dev, COMPUTE_CAPABILITY_MAJOR)
    minor = attribute(dev, COMPUTE_CAPABILITY_MINOR)
    major, minor
end

warpsize(dev::Int) = attribute(dev, WARP_SIZE)


text:
module OceanBoundaryLayerParameterizations

pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using OceanTurbulenceParameterEstimation
using OceanTurbulenceParameterEstimation.Models: set!
using CairoMakie, LaTeXStrings, OrderedCollections
using Oceananigans, FileIO

export
    # lesbrary_paths.jl
    TwoDaySuite, FourDaySuite, SixDaySuite, GeneralStrat,

    # one_dimensional_ensemble_model.jl
    OneDimensionalEnsembleModel,

    # EKI_hyperparameter_search.jl
    plot_stds_within_bounds, plot_prior_variance, plot_num_ensemble_members,
    plot_observation_noise_level, plot_prior_variance_and_obs_noise_level,

    # visualize_profile_predictions.jl
    visualize!,
    visualize_and_save!

include("utils/lesbrary_paths.jl")
include("utils/one_dimensional_ensemble_model.jl")
include("EKI_hyperparameter_search.jl")
include("visualize_profile_predictions_utils.jl")
include("utils/visualize_profile_predictions.jl")

end # module


text:
# TODO:
#   - allow Vn to be an arbitrary pair potential
#   - replace the polynomial with an arbitrary family of pair potentials
#   - (or even nbody?)

module EnvIPs

using StaticArrays
using JuLIP:              AbstractCalculator
using JuLIP.Potentials:   Shift,
                          @analytic,
                          @pot
using NBodyIPs:           NBodyDescriptor
using NBodyIPs.PolyBasis: nbpolys

import Base:              Dict,
                          ==,
                          convert,
                          hash
import NBodyIPs:          NBodyIP,
                          bodyorder,
                          fast,
                          combinebasis,
                          _decode_dict,
                          descriptor,
                          degree,
                          basisname,
                          BASIS

export envpolys, info

abstract type AbstractEnvIP{N} <: AbstractCalculator end

struct EnvIP{N, P, TVR, TVN} <: AbstractEnvIP{N}
   t::Int
   Vr::TVR     # N-body potential
   Vn::TVN     # neighbour counter
   str_Vn::String  # string describing the neighbour counter
   valN::Val{N}
   valP::Val{P}
end

@pot EnvIP

==(V1::EnvIP, V2::EnvIP) = ( (V1.t == V2.t) &&
                             (V1.Vr == V2.Vr) &&
                             (V1.str_Vn == V2.str_Vn) &&
                             (V1.valN == V2.valN) )

Dict(V::EnvIP) = Dict( "__id__" => "EnvIP",
                       "t" => V.t,
                       "Vr" => Dict(V.Vr),
                       "str_Vn" => V.str_Vn,
                       "cutoff_Vn" => cutoff(V.Vn)  )

EnvIP(D::Dict) = EnvIP( D["t"],
                        _decode_dict(D["Vr"]),
                        analyse_Vn(D["str_Vn"], D["cutoff_Vn"]),
                        D["str_Vn"] )

convert(::Val{:EnvIP}, D::Dict) = EnvIP(D)

EnvIP(t::Int, Vr, Vn, str_Vn::String) =
      EnvIP(t, Vr, Vn, str_Vn, Val(bodyorder(Vr)), Val(t))

function analyse_Vn(str_Vn, cutoff_Vn)
   Vn1 = eval(Meta.parse("@analytic r -> " * str_Vn))
   Vn = Shift(Val(1), Vn1, cutoff_Vn,
              Base.invokelatest(Vn1.f,   cutoff_Vn),
              Base.invokelatest(Vn1.f_d, cutoff_Vn),
              0.0)
   return Vn
end

function EnvIP(t, Vr, str_Vn::String, cutoff_Vn::AbstractFloat)
   Vn = analyse_Vn(str_Vn, cutoff_Vn)
   return EnvIP(t, Vr, Vn, str_Vn)
end

Vn(V::EnvIP) = V.Vn
Vr(V::EnvIP) = V.Vr

bodyorder(V::AbstractEnvIP{N}) where {N} = Int(N)

descriptor(V::EnvIP) = descriptor(V.Vr)

hash(::BASIS, V::EnvIP) = hash((EnvIP,
                                hash(BASIS(), V.Vr),
                                hash(V.str_Vn),
                                hash(V.t)))

function degree(V::EnvIP)
   if length(V.Vr) == 1
      return ( degree(V.Vr), V.t )
   end
   error("`degree` is only defined for `EnvIP` basis functions, length == 1")
end

basisname(::EnvIP) = "EnvIP"

# ----------------- generate basis / IP / convert ----------------

function envpolys(N::Integer, D::NBodyDescriptor, deg_poly,
                  Vn_descr, deg_n::Integer; kwargs...)
   B_poly = nbpolys(N, D, deg_poly...; kwargs...)
   B = EnvIP[]
   str_Vn = Vn_descr[1]
   Vn = analyse_Vn(Vn_descr...)
   for deg = 0:deg_n
      append!(B, [EnvIP(deg, Vr, Vn, str_Vn) for Vr in B_poly])
   end
   return [b for b in B]
end

function combinebasis(basis::AbstractVector{TV}, coeffs) where {TV <: EnvIP}
   # @assert isleaftype(TV)
   @assert all( b.t == basis[1].t for b in basis )
   # combine the Vr components of the basis functions
   # (we get to do this because all t (=P) are the same
   Vr = combinebasis( [b.Vr for b in basis], coeffs )
   return EnvIP(basis[1].t, Vr, basis[1].Vn, basis[1].str_Vn)
end

function info(B::Vector{T}; indent = 2) where T <: EnvIP
   ind = repeat(" ", indent)
   println(ind * "EnvIP with P = $(B[1].t)")
   println(ind * "           Vn : $(B[1].str_Vn)")
   println(ind * "           Vr : ...")
   info([ b.Vr for b in B ], indent = indent+5)
end

fast(V::EnvIP) = EnvIP(V.t, fast(V.Vr), V.Vn, V.str_Vn)

# ========================================================
#    Re-Implement Fast Version of EnvIPs
# ========================================================

struct EnvPoly{N, TVR, TVN} <: AbstractEnvIP{N}
   Vr::Vector{TVR}    # N-body potentials multiplied by n^j, j = 0, 1, ...
   Vn::TVN            # neighbour counter
   str_Vn::String     # string describing the neighbour counter
   valN::Val{N}
end

@pot EnvPoly

envdegree(V::EnvPoly) = length(V.P) - 1

==(V1::EnvPoly, V2::EnvPoly) = ( (V1.Vr == V2.Vr) &&
                                 (V1.str_Vn == V2.str_Vn) &&
                                 (V1.valN == V2.valN) )

# Dict(V::EnvPoly) = Dict( "__id__" => "EnvPoly",
#                          "Vr" => Dict.(V.Vr),
#                          "str_Vn" => V.str_Vn,
#                          "cutoff_Vn" => cutoff(V.Vn)  )
#
# EnvPoly(D::Dict) = EnvPoly( _decode_dict.(D["Vr"]),
#                             analyse_Vn(D["str_Vn"], D["cutoff_Vn"]),
#                             D["str_Vn"] )

convert(::Val{:EnvPoly}, D::Dict) = EnvPoly(D)

EnvPoly(t, Vr, str_Vn::String, cutoff_Vn::AbstractFloat) =
   EnvPoly(t, Vr, analyse_Vn(str_Vn, cutoff_Vn), str_Vn)

descriptor(V::EnvPoly) = descriptor(V.Vr[1])

hash(::BASIS, V::EnvPoly) = hash(( hash(EnvPoly),
                                   hash.(Ref(BASIS()), V.Vr),
                                   hash(V.str_Vn) ))

Vn(V::EnvPoly) = V.Vn
Vr(V::EnvPoly) = V.Vr

function EnvPoly(Vr::AbstractVector, Vn, str_Vn::String)
   # for now we put the burden on the user to create a Vr vector where
   # the types match But we can still collect them here for convenience
   Vr_coll = [ v for v in Vr ]
   if !isconcretetype(typeof(Vr_coll))
      error("""For now, `EnvPoly` can only be constructed from an array of
               potentials that all share the same concrete type.""")
   end
   return EnvPoly(Vr_coll, Vn, str_Vn, Val(bodyorder(Vr[1])))
end

# function degree(V::EnvIP)
#    if length(V.Vr) == 1
#       return ( degree(V.Vr), V.t )
#    end
#    error("`degree` is only defined for `EnvIP` basis functions, length == 1")
# end

basisname(::EnvPoly) = "EnvPoly"

# all evaluation and assembly in here:
include("eval_environ.jl")

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# ---
# title: 1368. Minimum Cost to Make at Least One Valid Path in a Grid
# id: problem1368
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/minimum-cost-to-make-at-least-one-valid-path-in-a-grid/description/>
# hidden: true
# ---
#
# Given a _m_ x _n_ `grid`. Each cell of the `grid` has a sign pointing to the
# next cell you should visit if you are currently in this cell. The sign of
# `grid[i][j]` can be:
#
#   * **1** which means go to the cell to the right. (i.e go from `grid[i][j]` to `grid[i][j + 1]`)
#   * **2** which means go to the cell to the left. (i.e go from `grid[i][j]` to `grid[i][j - 1]`)
#   * **3** which means go to the lower cell. (i.e go from `grid[i][j]` to `grid[i + 1][j]`)
#   * **4** which means go to the upper cell. (i.e go from `grid[i][j]` to `grid[i - 1][j]`)
#
# Notice that there could be some **invalid signs** on the cells of the `grid`
# which points outside the `grid`.
#
# You will initially start at the upper left cell `(0,0)`. A valid path in the
# grid is a path which starts from the upper left cell `(0,0)` and ends at the
# bottom-right cell `(m - 1, n - 1)` following the signs on the grid. The valid
# path **doesn 't have to be the shortest**.
#
# You can modify the sign on a cell with `cost = 1`. You can modify the sign on
# a cell **one time only**.
#
# Return _the minimum cost_ to make the grid have at least one valid path.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/13/grid1.png)
#
#
#
#     Input: grid = [[1,1,1,1],[2,2,2,2],[1,1,1,1],[2,2,2,2]]
#     Output: 3
#     Explanation: You will start at point (0, 0).
#     The path to (3, 3) is as follows. (0, 0) --> (0, 1) --> (0, 2) --> (0, 3) change the arrow to down with cost = 1 --> (1, 3) --> (1, 2) --> (1, 1) --> (1, 0) change the arrow to down with cost = 1 --> (2, 0) --> (2, 1) --> (2, 2) --> (2, 3) change the arrow to down with cost = 1 --> (3, 3)
#     The total cost = 3.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/13/grid2.png)
#
#
#
#     Input: grid = [[1,1,3],[3,2,2],[1,1,4]]
#     Output: 0
#     Explanation: You can follow the path from (0, 0) to (2, 2).
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/13/grid3.png)
#
#
#
#     Input: grid = [[1,2],[4,3]]
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: grid = [[2,2,2],[2,2,2]]
#     Output: 3
#
#
# **Example 5:**
#
#
#
#     Input: grid = [[4]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `m == grid.length`
#   * `n == grid[i].length`
#   * `1 <= m, n <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module PostgreSQL

export extract

include("sql_exports.jl")
include("sql_imports.jl") # DBMS Structured SubQuery _to_sql _placeholder _placeholders
using .Octo.Queryable: extract #
using .Octo: @sql_keywords, @sql_functions

const DatabaseID = DBMS.PostgreSQL

"""
    to_sql(query::Structured)::String
"""
to_sql(query::Structured)::String = _to_sql(DatabaseID(), query)

"""
    to_sql(subquery::SubQuery)::String
"""
to_sql(subquery::SubQuery)::String = _to_sql(DatabaseID(), subquery)

placeholder(nth::Int) = _placeholder(DatabaseID(), nth)
placeholders(dims::Int) = _placeholders(DatabaseID(), dims)

_placeholder(db::DatabaseID, nth::Int) = PlaceHolder("\$$nth")
_placeholders(db::DatabaseID, dims::Int) = Enclosed([PlaceHolder("\$$x") for x in 1:dims])

export         AUTOCOMMIT, BIGSERIAL, COPY, CURRENT_DATE, DESCRIBE, EXPLAIN, FALSE, LATERAL, SEQUENCE, SERIAL, TRUE, WINDOW
@sql_keywords  AUTOCOMMIT  BIGSERIAL  COPY  CURRENT_DATE  DESCRIBE  EXPLAIN  FALSE  LATERAL  SEQUENCE  SERIAL  TRUE  WINDOW

export         COALESCE, NOW
@sql_functions COALESCE  NOW

using .Octo.AdapterBase: FromItem, SqlPart, _sqlrepr
import .Octo.AdapterBase: sqlrepr
function sqlrepr(db::DatabaseID, clause::FromItem)::SqlPart
    _sqlrepr(db, clause; with_as=false)
end

end # Octo.Adapters.PostgreSQL


text:
"""
    integer_contractor(x::Interval)

Return the integers enclosed in the interval `x`.
"""

function integer_contractor(x::Interval)
    a = floor(x.lo)+1
    b = floor(x.hi)

    a > b && return emptyinterval(x)

    return Interval(a, b)
end

### Transformations on IntervalBoxes

"""Reflect in mirror at position x_mirror"""
function reflect_x(x_mirror)
    X -> begin
        x, y = X
        x = 2*x_mirror - x

        return IntervalBox(x, y)
    end
end

"""
    translate(α)

Returns a function that shifts (translates) a 2D `IntervalBox` in `x` (the first coordinate).
"""
function translate(α)
    X -> begin
        x, y = X
        return IntervalBox(x - α, y)
    end
end

odd(X::IntervalBox) = ( (x,y) = X; IntervalBox(-x, -y) )

## Transformations on Contractors

"""
Translation of a Contractor `C` by `α`.
Uses `inv(op) ∘ C ∘ op`
"""
translate(C, α) = translate(-α) ∘ C ∘ translate(α)

"""
Symmetric part of a Contractor, via an involution `op`
(i.e. such that `inv(op) == op`).
"""
symmetrise(C, op) = op ∘ C ∘ op

"Periodize the contractor C"
function periodise(C, period)

    X -> begin
        x, y = X

        x2 = -∞..∞
        x2, y = C(IntervalBox(x2, y))

        isempty(IntervalBox(x2, y)) && return(IntervalBox(∅, ∅))

        # periods where the periodization of x intersects with x2:
        periods = integer_contractor((x - x2) / period)

        isempty(periods) && return(IntervalBox(∅, ∅))

        x3 = x2 + periods*period
        x = x ∩ x3

        return IntervalBox(x, y)
    end

end


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

@doc raw"""RuntimeClassStrategyOptions define the strategy that will dictate the allowable RuntimeClasses for a pod.

    IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions(;
        allowedRuntimeClassNames=nothing,
        defaultRuntimeClassName=nothing,
    )

    - allowedRuntimeClassNames::Vector{String} : allowedRuntimeClassNames is a whitelist of RuntimeClass names that may be specified on a pod. A value of \"*\" means that any RuntimeClass name is allowed, and must be the only item in the list. An empty list requires the RuntimeClassName field to be unset.
    - defaultRuntimeClassName::String : defaultRuntimeClassName is the default RuntimeClassName to set on the pod. The default MUST be allowed by the allowedRuntimeClassNames list. A value of nil does not mutate the Pod.
"""
mutable struct IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions <: SwaggerModel
    allowedRuntimeClassNames::Any # spec type: Union{ Nothing, Vector{String} } # spec name: allowedRuntimeClassNames
    defaultRuntimeClassName::Any # spec type: Union{ Nothing, String } # spec name: defaultRuntimeClassName

    function IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions(;allowedRuntimeClassNames=nothing, defaultRuntimeClassName=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions, Symbol("allowedRuntimeClassNames"), allowedRuntimeClassNames)
        setfield!(o, Symbol("allowedRuntimeClassNames"), allowedRuntimeClassNames)
        validate_property(IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions, Symbol("defaultRuntimeClassName"), defaultRuntimeClassName)
        setfield!(o, Symbol("defaultRuntimeClassName"), defaultRuntimeClassName)
        o
    end
end # type IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions

const _property_map_IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions = Dict{Symbol,Symbol}(Symbol("allowedRuntimeClassNames")=>Symbol("allowedRuntimeClassNames"), Symbol("defaultRuntimeClassName")=>Symbol("defaultRuntimeClassName"))
const _property_types_IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions = Dict{Symbol,String}(Symbol("allowedRuntimeClassNames")=>"Vector{String}", Symbol("defaultRuntimeClassName")=>"String")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions)
    (getproperty(o, Symbol("allowedRuntimeClassNames")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1RuntimeClassStrategyOptions }, name::Symbol, val)
end


text:
LARGE_ANCESTOR = [
    NOP_1(),
    NOP_1(),
    NOP_1(),
    NOP_1(),
    ZERO(),
    OR_1(),
    SH_L(),
    SH_L(),
    MOV_CD(),
    ADR_BACK_A(),
    NOP_0(),
    NOP_0(),
    NOP_0(),
    NOP_0(),
    SUB_AC(),
    MOV_AB(),
    ADR_FORW_A(),
    NOP_0(),
    NOP_0(),
    NOP_0(),
    NOP_1(),
    INC_A(),
    SUB_AB(),
    NOP_1(),
    NOP_1(),
    NOP_0(),
    NOP_1(),
    MALLOC(),
    CALL(),
    NOP_0(),
    NOP_0(),
    NOP_1(),
    NOP_1(),
    DIVIDE(),
    JMP(),
    NOP_0(),
    NOP_0(),
    NOP_1(),
    NOP_0(),
    IF_0(),
    NOP_1(),
    NOP_1(),
    NOP_0(),
    NOP_0(),
    PUSH_A(),
    PUSH_B(),
    PUSH_C(),
    NOP_1(),
    NOP_0(),
    NOP_1(),
    NOP_0(),
    MOVI_BA(),
    DEC_C(),
    IF_0(),
    JMP(),
    NOP_0(),
    NOP_1(),
    NOP_0(),
    NOP_0(),
    INC_A(),
    INC_B(),
    JMP(),
    NOP_0(),
    NOP_1(),
    NOP_0(),
    NOP_1(),
    IF_0(),
    NOP_1(),
    NOP_0(),
    NOP_1(),
    NOP_1(),
    POP_C(),
    POP_B(),
    POP_A(),
    RET(),
    NOP_1(),
    NOP_1(),
    NOP_1(),
    NOP_0(),
    IF_0()
]

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using MolecularGraph: lookahead, forward!, backtrack!

@testset "smarts.base" begin

    state1 = SmilesParser("C1SC(C=O)CCC1", true)
    state2 = SmartsParser(raw"*OC$([Cv4H2+0])", false)
    @test read(state1) == 'C'
    @test lookahead(state1, 2) == 'S'
    @test read(state2) == '*'
    forward!(state2)
    @test read(state2) == 'O'
    forward!(state2, 2)
    @test read(state2) == '$'
    backtrack!(state2)
    @test lookahead(state2, 1) == '$'
    @test !state2.done
    forward!(state2, 13)
    @test state2.done

end # smarts.base


text:
module GeoCost

using FillArrays
using StaticArrays
using Distances
using OffsetArrays
using DataStructures
using Statistics

const sqrt2 = sqrt(2.0)
const neib_8 = @SMatrix[1. 1 1; 1 0 1; 1 1 1]
const distance_8 = @SMatrix[sqrt2 1 sqrt2; 1 Inf 1; sqrt2 1 sqrt2]
const distance_4 = @SMatrix[Inf 1 Inf; 1 Inf 1; Inf 1 Inf]
const Δ = CartesianIndex(1, 1)

# Write your package code here.
"""Total friction distance spread from `points`."""
function spread(points::Matrix{<:Real}, initial::Matrix{<:AbstractFloat}, friction::Matrix{<:AbstractFloat}; res=1, limit=Inf)

    ofriction = OffsetMatrix(fill(Inf, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    ofriction[begin + 1:end - 1,begin + 1:end - 1] .= friction

    result = OffsetMatrix(fill(limit, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    r = @view result[1:end - 1,1:end - 1]
    locations = points .> 0
    r[locations] .= initial[locations]

    # Construct stack for locations
    mask = OffsetMatrix(trues(size(points) .+ 2), UnitRange.(0, size(points) .+ 1))
    mask[begin + 1:end - 1,begin + 1:end - 1] .= false

    II = CartesianIndices(size(points))
    stack = Deque{CartesianIndex}()
    for I in II[locations]
        push!(stack, I)
    end

    # Step 1: Set the distance of the starting node to 0 and the distances of all other nodes to the highest value possible.
    sdata = zeros(MMatrix{3,3})
    mcell = MMatrix{3,3}(false, false, false, false, false, false, false, false, false)

    # fcdata = zeros(MMatrix{3,3})
    # rcdata = zeros(MMatrix{3,3})

    # Step 3: For each of the active node’s adjacent neighbors, set its distance to whichever is
    # less: its current distance value or the sum of the distance of the active node plus the
    # weight of the arc from the active node to that neighbor.
    while !isempty(stack)
        spread!(stack, mask, result, ofriction, sdata, mcell, res)
    end
    r
end

function spread2(points::Matrix{<:Real}, initial::Matrix{<:AbstractFloat}, friction::Matrix{<:AbstractFloat}; res=1, limit=Inf)

    ofriction = OffsetMatrix(fill(Inf, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    ofriction[begin + 1:end - 1,begin + 1:end - 1] .= friction

    result = OffsetMatrix(fill(limit, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    r = @view result[1:end - 1,1:end - 1]
    locations = points .> 0
    r[locations] .= initial[locations]

    mask = OffsetMatrix(trues(size(points) .+ 2), UnitRange.(0, size(points) .+ 1))
    mask[begin + 1:end - 1,begin + 1:end - 1] .= false

    minval, minidx = [0.], [CartesianIndex(1, 1)]
    x = @MMatrix zeros(3, 3)

    II = CartesianIndices(size(points))
    for I ∈ II
        patch = I - Δ:I + Δ

        rdata = view(result, patch)
        fdata = view(ofriction, patch)

        x .= (fdata .+ fdata[2,2]) .* res ./ 2 .* distance_8 .+ rdata
        findmin!(minval, minidx, x)
        rdata[2,2] = min(rdata[2,2], minval[1])
    end
    for I ∈ reverse(II)
        patch = I - Δ:I + Δ

        rdata = view(result, patch)
        fdata = view(ofriction, patch)
        x .= (fdata .+ fdata[2,2]) .* res ./ 2 .* distance_8 .+ rdata
        findmin!(minval, minidx, x)
        rdata[2,2] = min(rdata[2,2], minval[1])
    end
    for I ∈ II
        patch = I - Δ:I + Δ

        rdata = view(result, patch)
        fdata = view(ofriction, patch)

        x .= (fdata .+ fdata[2,2]) .* res ./ 2 .* distance_8 .+ rdata
        findmin!(minval, minidx, x)
        rdata[2,2] = min(rdata[2,2], minval[1])
    end
    r
end

function spread3(points::Matrix{<:Real}, initial::Matrix{<:AbstractFloat}, friction::Matrix{<:AbstractFloat}; res=1, limit=Inf)

    ofriction = OffsetMatrix(fill(Inf, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    ofriction[begin + 1:end - 1,begin + 1:end - 1] .= friction

    result = OffsetMatrix(fill(limit, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    r = @view result[1:end - 1,1:end - 1]
    locations = points .> 0
    r .= initial

    mask = OffsetMatrix(trues(size(points) .+ 2), UnitRange.(0, size(points) .+ 1))
    mask[begin + 1:end - 1,begin + 1:end - 1] .= false

    minval, minidx = [0.], [CartesianIndex(1, 1)]
    x = @MMatrix zeros(3, 3)

    II = CartesianIndices(size(points))
    for I ∈ II
        patch = I - Δ:I + Δ

        rdata = view(result, patch)
        fdata = view(ofriction, patch)

        x .= (fdata .+ fdata[2,2]) .* res ./ 2 .* distance_8 .+ rdata
        findmin!(minval, minidx, x)
        rdata[2,2] = min(rdata[2,2], minval[1])
    end
    for I ∈ reverse(II)
        patch = I - Δ:I + Δ

        rdata = view(result, patch)
        fdata = view(ofriction, patch)
        x .= (fdata .+ fdata[2,2]) .* res ./ 2 .* distance_8 .+ rdata
        findmin!(minval, minidx, x)
        rdata[2,2] = min(rdata[2,2], minval[1])
    end
    for I ∈ II
        patch = I - Δ:I + Δ

        rdata = view(result, patch)
        fdata = view(ofriction, patch)

        x .= (fdata .+ fdata[2,2]) .* res ./ 2 .* distance_8 .+ rdata
        findmin!(minval, minidx, x)
        rdata[2,2] = min(rdata[2,2], minval[1])
    end
    r
end

function spread!(stack, mask, result, ofriction, sdata, mcell, res)
    I = popfirst!(stack)
    mask[I] = true
    patch = I - Δ:I + Δ

    rdata = view(result, patch)
    fdata = view(ofriction, patch)

    # fcdata .= fdata[2,2]
    # rcdata .= rdata[2,2]

    # New distance is cell_distance + average friction values
    for i ∈ eachindex(sdata)
        sdata[i] = muladd(fdata[i] + fdata[2,2], res / 2 * distance_8[i], rdata[2,2])
        mcell[i] = sdata[i] < rdata[i]  # cells where new distance is lower
    end
    rdata[mcell] .= sdata[mcell]
    result[patch] .= rdata

    # Add new cells to stack
    for I in patch[mcell]
        mask[I] || push!(stack, I)
    end
end

"""Optimized (and more accurate) function based on the same friction everywhere."""
function spread(points::Matrix{<:AbstractFloat}, initial::AbstractFloat, friction::AbstractFloat; distance=Euclidean(), res=1.0)
    locations = points .> 0
    I = CartesianIndices(size(points))

    result = fill(Inf, size(points))
    for location ∈ I[locations]
        for cell ∈ I
            result[cell] = min(evaluate(distance, location.I, cell.I) * res * friction + initial, result[cell])
        end
    end
    # result .+ initial
    m = .~isfinite.(points)
    result[m] = points[m]
    return result
end

"""Optimized (and more accurate) function based on the same friction everywhere."""
function spread(points::Matrix{<:AbstractFloat}, initial::Matrix{<:AbstractFloat}, friction::Real; distance=Euclidean(), res=1.0)
    locations = points .> 0
    I = CartesianIndices(size(points))

    result = fill(Inf, size(points))
    for location ∈ I[locations]
        for cell ∈ I
            result[cell] = min(evaluate(distance, location.I, cell.I) * res * friction + initial[location], result[cell])
        end
    end
    m = .~isfinite.(points)
    result[m] .= points[m]
    return result
end

"""
    roughness(dem::Matrix{<:AbstractFloat})

Roughness is the largest inter-cell difference of a central pixel and its surrounding cell, as defined in Wilson et al (2007, Marine Geodesy 30:3-35).
"""
function roughness(dem::Matrix{<:AbstractFloat})

    ex_dem = OffsetMatrix(fill(Inf, size(dem) .+ 2), UnitRange.(0, size(dem) .+ 1))
    # Update center
    ex_dem[begin + 1:end - 1,begin + 1:end - 1] .= dem
    # Set edges to mirror center
    ex_dem[begin, begin + 1:end - 1] .= dem[begin, :]
    ex_dem[end, begin + 1:end - 1] .= dem[end, :]
    ex_dem[begin + 1:end - 1, begin] .= dem[:, begin]
    ex_dem[begin + 1:end - 1, end] .= dem[:, end]
    # Set corners to mirror corners of center
    ex_dem[begin, begin] = dem[begin, begin]
    ex_dem[begin, end] = dem[begin, end]
    ex_dem[end, begin] = dem[end, begin]
    ex_dem[end, end] = dem[end, end]

    roughness = similar(dem)

    x = @MMatrix zeros(3, 3)

    @inbounds for I ∈ CartesianIndices(size(roughness))
        patch = I - Δ:I + Δ
        rdata = view(ex_dem, patch)

        x .= rdata .- rdata[2,2]
        roughness[I] = maximum(abs.(x))
    end
    roughness
end

"""
    TPI(dem::Matrix{<:AbstractFloat})

TPI stands for Topographic Position Index, which is defined as the difference between a central pixel and the mean of its surrounding cells (see Wilson et al 2007, Marine Geodesy 30:3-35).
"""
function TPI(dem::Matrix{<:AbstractFloat})

    ex_dem = OffsetMatrix(fill(Inf, size(dem) .+ 2), UnitRange.(0, size(dem) .+ 1))
    # Update center
    ex_dem[begin + 1:end - 1,begin + 1:end - 1] .= dem
    # Set edges to mirror center
    ex_dem[begin, begin + 1:end - 1] .= dem[begin, :]
    ex_dem[end, begin + 1:end - 1] .= dem[end, :]
    ex_dem[begin + 1:end - 1, begin] .= dem[:, begin]
    ex_dem[begin + 1:end - 1, end] .= dem[:, end]
    # Set corners to mirror corners of center
    ex_dem[begin, begin] = dem[begin, begin]
    ex_dem[begin, end] = dem[begin, end]
    ex_dem[end, begin] = dem[end, begin]
    ex_dem[end, end] = dem[end, end]

    tpi = similar(dem)

    x = @MMatrix zeros(3, 3)

    @inbounds for I ∈ CartesianIndices(size(tpi))
        patch = I - Δ:I + Δ
        rdata = view(ex_dem, patch)

        x .= rdata .* neib_8
        tpi[I] = rdata[2,2] - mean(x)
    end
    tpi
end

"""
    TRI(dem::Matrix{<:AbstractFloat})

TRI stands for Terrain Ruggedness Index, which measures the difference between a central pixel and its surrounding cells.
This algorithm uses the square root of the sum of the square of the difference between a central pixel and its surrounding cells.
This is recommended for terrestrial use cases.
"""
function TRI(dem::Matrix{<:AbstractFloat})

    ex_dem = OffsetMatrix(fill(Inf, size(dem) .+ 2), UnitRange.(0, size(dem) .+ 1))

    # Update center
    ex_dem[begin + 1:end - 1,begin + 1:end - 1] .= dem

    # Set edges to mirror center
    ex_dem[begin, begin + 1:end - 1] .= dem[begin, :]
    ex_dem[end, begin + 1:end - 1] .= dem[end, :]
    ex_dem[begin + 1:end - 1, begin] .= dem[:, begin]
    ex_dem[begin + 1:end - 1, end] .= dem[:, end]

    # Set corners to mirror corners of center
    ex_dem[begin, begin] = dem[begin, begin]
    ex_dem[begin, end] = dem[begin, end]
    ex_dem[end, begin] = dem[end, begin]
    ex_dem[end, end] = dem[end, end]

    tri = similar(dem)

    x = @MMatrix zeros(3, 3)

    @inbounds for I ∈ CartesianIndices(size(tri))
        patch = I - Δ:I + Δ
        rdata = view(ex_dem, patch)

        x .= (rdata .- rdata[2,2]).^2
        tri[I] = sqrt(sum(x))
    end
    tri
end

"""
```
B, flags = pmf(A; ωₘ, slope, dhₘ, dh₀, cellsize)
```
Applies the progressive morphological filter by [Zhang et al. (2003)] to `A`.
# Output
- `B::Array{T,2}` Maximum allowable values
- `flags::Array{Float64,2}` A sized array with window sizes if filtered, zero if not filtered.
Afterwards, one can retrieve the resulting mask for `A` by `A .<= B` or `flags .== 0.`.
# Arguments
- `A::Array{T,2}` Input Array
- `ωₘ::Float64=20.` Maximum window size [m]
- `slope::Float64=0.01` Terrain slope [m/m]
- `dhₘ::Float64=2.5` Maximum elevation threshold [m]
- `dh₀::Float64=0.2` Initial elevation threshold [m]
- `cellsize::Float64=1.` Cellsize in [m]
[Zhang et al. (2003)] Zhang, Keqi, Shu-Ching Chen, Dean Whitman, Mei-Ling Shyu, Jianhua Yan, and Chengcui Zhang. “A Progressive Morphological Filter for Removing Nonground Measurements from Airborne LIDAR Data.” IEEE Transactions on Geoscience and Remote Sensing 41, no. 4 (2003): 872–82. [https://doi.org/10.1109/TGRS.2003.810682].
"""
function pmf(A::Array{T,2};
    ωₘ::Float64=20.,
    slope::Float64=0.01,
    dhₘ::Float64=2.5,
    dh₀::Float64=0.2,
    cellsize::Float64=1.0) where T <: Real

    # Compute windowsizes and thresholds
    ωₘ = round(Int, ωₘ / cellsize)
    κ_max = floor(Int, log2(ωₘ - 1))  # determine # iterations based on exp growth
    windowsizes = Int.(exp2.(1:κ_max)) .+ 1

    # Compute tresholds
    dwindows = vcat(windowsizes[1], windowsizes)  # prepend first element so we get 0 as diff
    window_diffs = [dwindows[i] - dwindows[i - 1] for i in 2:length(dwindows)]
    height_tresholds = [min(dhₘ, slope * window_diff * cellsize + dh₀) for window_diff in window_diffs]

    # Set up arrays
    Af = copy(A)  # array to be morphed
    nan_mask = isnan.(Af)
    Af[nan_mask] .= Inf  # Replace NaN with Inf, as to always filter these

    B = copy(A)  # max_elevation raster
    out = copy(A)  # max_elevation raster

    flags = zeros(size(A))  # 0 = ground, other values indicate window size
    flags[nan_mask] .= NaN

    mask = falses(size(A))

    # Iterate over window sizes and height tresholds
    for (ωₖ, dhₜ) in zip(windowsizes, height_tresholds)
        opening!(Af, ωₖ, out)
        for I in eachindex(A)
            mask[I] = (A[I] - Af[I]) > dhₜ
        end
        for I in eachindex(flags)
            if mask[I] && flags[I] == 0
                flags[I] = ωₖ
            end
        end
        # @info "PMF with window size $ωₖ and threshold $dhₜ filters $(sum(mask)) cells."
        B .= min.(B, Af .+ dhₜ)
    end

    B, flags
end

# First discussed here https://github.com/JuliaImages/ImageFiltering.jl/issues/179
function mapwindow!(f, img, window, out)
    R = CartesianIndices(img)
    I_first, I_last = first(R), last(R)
    Δ = CartesianIndex(ntuple(x -> window ÷ 2, ndims(img)))
    @inbounds @simd for I in R
        patch = max(I_first, I - Δ):min(I_last, I + Δ)
        out[I] = f(view(img, patch))
    end
    out
end

"""Apply the opening operation to `A` with window size `ω`."""
function opening!(A::Array{T,2}, ω::Integer, out::Array{T,2}) where T <: Real
    mapwindow!(minimum, A, ω, out)  # erosion
    mapwindow!(maximum, out, ω, A)  # dilation
    A
end

export spread, spread2, spread3, roughness, TPI, TRI, pmf

end  # module


text:
"Shared docstring part for the MCMC API."
const DOC_MCMC_ARGS =
"""
# Arguments
- `rng`: the random number generator, eg `Random.GLOBAL_RNG`.
- `ℓ`: the log density, supporting the API of the `LogDensityProblems` package
- `N`: the number of samples for inference, after the warmup.
# Keyword arguments
- `initialization`: see below.
- `warmup_stages`: a sequence of warmup stages. See [`default_warmup_stages`](@ref) and
  [`fixed_stepsize_warmup_stages`](@ref); the latter requires an `ϵ` in initialization.
- `algorithm`: see [`NUTS`](@ref). It is very unlikely you need to modify
  this, except perhaps for the maximum depth.
- `reporter`: how progress is reported. By default, verbosely for interactive sessions using
  the log message mechanism (see [`LogProgressReport`](@ref), and no reporting for
  non-interactive sessions (see [`NoProgressReport`](@ref)).
# Initialization
The `initialization` keyword argument should be a `NamedTuple` which can contain the
following fields (all of them optional and provided with reasonable defaults):
$(DOC_INITIAL_WARMUP_ARGS)
"""

# """
# $(SIGNATURES)
# Perform MCMC with NUTS, keeping the warmup results. Returns a `NamedTuple` of
# - `initial_warmup_state`, which contains the initial warmup state
# - `warmup`, an iterable of `NamedTuple`s each containing fields
#     - `stage`: the relevant warmup stage
#     - `results`: results returned by that warmup stage (may be `nothing` if not applicable,
#       or a chain, with tree statistics, etc; see the documentation of stages)
#     - `warmup_state`: the warmup state *after* the corresponding stage.
# - `final_warmup_state`, which contains the final adaptation after all the warmup
# - `inference`, which has `chain` and `tree_statistics`, see [`mcmc_with_warmup`](@ref).
# !!! warning
#     This function is not (yet) exported because the the warmup interface may change with
#     minor versions without being considered breaking. Recommended for interactive use.
# $(DOC_MCMC_ARGS)
# """
# function mcmc_keep_warmup(rng::AbstractRNG, tree::Tree, ℓ, N::Integer;
#                           initialization = (),
#                           warmup_stages = default_warmup_stages(),
#                           algorithm = NUTS(),
#                           reporter = default_reporter())
#     sampling_logdensity = SamplingLogDensity(rng, ℓ, algorithm, reporter)
#     sptr, initial_warmup_state = initialize_warmup_state(rng, sptr, ℓ; initialization...)
#     warmup, warmup_state = _warmup(sampling_logdensity, warmup_stages, initial_warmup_state)
#     inference = mcmc(sampling_logdensity, N, warmup_state)
#     (initial_warmup_state = initial_warmup_state, warmup = warmup,
#      final_warmup_state = warmup_state, inference = inference)
# end

# """
# $(SIGNATURES)
# Perform MCMC with NUTS, including warmup which is not returned. Return a `NamedTuple` of
# - `chain`, a vector of positions from the posterior
# - `tree_statistics`, a vector of tree statistics
# - `κ` and `ϵ`, the adapted metric and stepsize.
# $(DOC_MCMC_ARGS)
# # Usage examples
# Using a fixed stepsize:
# ```julia
# mcmc_with_warmup(rng, ℓ, N;
#                  initialization = (ϵ = 0.1, ),
#                  warmup_stages = fixed_stepsize_warmup_stages())
# ```
# Starting from a given position `q₀` and kinetic energy scaled down (will still be adapted):
# ```julia
# mcmc_with_warmup(rng, ℓ, N;
#                  initialization = (q = q₀, κ = GaussianKineticEnergy(5, 0.1)))
# ```
# Using a dense metric:
# ```julia
# mcmc_with_warmup(rng, ℓ, N;
#                  warmup_stages = default_warmup_stages(; M = Symmetric))
# ```
# Disabling the optimization step:
# ```julia
# mcmc_with_warmup(rng, ℓ, N;
#                  warmup_stages = default_warmup_stages(; local_optimization = nothing,
#                                                          M = Symmetric))
# ```
# """
# function mcmc_with_warmup(rng, ℓ, N; initialization = (),
#                           warmup_stages = default_warmup_stages(),
#                           algorithm = NUTS(), reporter = default_reporter())
#     @unpack final_warmup_state, inference =
#         mcmc_keep_warmup(rng, ℓ, N; initialization = initialization,
#                          warmup_stages = warmup_stages, algorithm = algorithm,
#                          reporter = reporter)
#     @unpack κ, ϵ = final_warmup_state
#     (inference..., κ = κ, ϵ = ϵ)
# end

function mcmc_with_warmup!(
    rng::AbstractRNG, sptr::StackPointer, chain::AbstractMatrix{T}, tree_statistics::AbstractVector{TreeStatisticsNUTS}, ℓ::AbstractProbabilityModel{D}, N::Int = size(chain,2);
    initialization = (), warmup_stages = default_warmup_stages(), algorithm = NUTS(), reporter = default_reporter()
) where {D,T}

    # We allocate the tree here.
    sampling_logdensity = SamplingLogDensity(rng, ℓ, algorithm, reporter)
    sptr, ϵs = DynamicPtrVector{Float64}(sptr, maximum(length, warmup_stages))
    tree, initial_warmup_state = initialize_warmup_state(rng, sptr, ℓ; initialization...)
    warmup_state = _warmup!(tree, chain, tree_statistics, ϵs, sampling_logdensity, warmup_stages, initial_warmup_state)
    mcmc!(tree, chain, tree_statistics, sampling_logdensity, N, warmup_state)

end

function mcmc_with_warmup(
    ℓ::AbstractProbabilityModel{D}, N; δ::Float64 = 0.8, initialization = (),
    warmup_stages = default_warmup_stages(stepsize_adaptation=DualAveraging(δ=δ)),
    algorithm = NUTS(), reporter = default_reporter()
) where {D}
    sptr = stackpointer()
    nwarmup = maximum(length, warmup_stages)
    NS = max(N, nwarmup)
    L = VectorizationBase.align(D, Float64)
    chain = DynamicPaddedMatrix{Float64}(undef, (D, NS), L)
    tree_statistics = DynamicPaddedVector{TreeStatisticsNUTS}(undef, NS)

    chain_m = DynamicPtrMatrix{Float64}(pointer(chain), (D, NS), L)
    tree_statistic = DynamicPtrVector{TreeStatisticsNUTS}(pointer(tree_statistics), NS)
    mcmc_with_warmup!(
        first(GLOBAL_PCGs), sptr, chain_m, tree_statistic, ℓ, N;
        initialization = initialization, warmup_stages = warmup_stages, algorithm = algorithm, reporter = reporter
    )
    chain, tree_statistics
end

function threaded_mcmc(
    ℓ::AbstractProbabilityModel{D}, N; δ::Float64 = 0.8, initialization = (),
    warmup_stages = default_warmup_stages(stepsize_adaptation=DualAveraging(δ=δ)),
    algorithm = NUTS(), reporter = NoProgressReport(),
    nchains = NTHREADS[]
) where {D}
    sptr = STACK_POINTER_REF[]
    LSS = LOCAL_STACK_SIZE[]
    nwarmup = maximum(length, warmup_stages)

    NS = max(N,nwarmup)
    L = VectorizationBase.align(D, Float64)
    NSA = VectorizationBase.align(NS, Float64)
    chains = DynamicPaddedArray{Float64}(undef, (D, NS, nchains), L)
    tree_stride = NSA + VectorizationBase.CACHELINE_SIZE >> 3
    tree_statistics = DynamicPaddedMatrix{TreeStatisticsNUTS}(undef, (NS, nchains), tree_stride)

    chain_ptr = pointer(chains)
    stat_ptr = pointer(tree_statistics)

    Threads.@threads for t in 0:nchains-1
        chain = DynamicPtrMatrix{Float64}(chain_ptr + t*8NS*L, (D, NS), L)
        tree_statistic = DynamicPtrVector{TreeStatisticsNUTS}(stat_ptr + t*tree_stride*sizeof(TreeStatisticsNUTS), (NS,), NSA)
        mcmc_with_warmup!(
            GLOBAL_PCGs[t+1], sptr + t*LSS, chain, tree_statistic, ℓ, N;
            initialization = initialization, warmup_stages = warmup_stages, algorithm = algorithm, reporter = reporter
        )
    end
    chains, tree_statistics
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://twitter.com/mkashi/status/1434847564256268291

# %%
ENV["LINES"] = 1000

# %%
function msb64(x)
    ε = eps(Float64)
    q = (1/ε + 1) * x
    m = q - (1 - ε/2) * q
    Int64(m)
end

@show eps(Float64) == 1/2^52

[(x = 2^k; (k, log2(msb64(x - 2)), log2(msb64(x - 1)), log2(msb64(x)))) for k in 1:54]

# %%
using Quadmath

function msb128(x)
    ε = eps(Float128)
    q = (1/ε + 1) * x
    m = q - (1 - ε/2) * q
    Int128(m)
end

@show eps(Float128) == 1/Int128(2)^112

[(x = Int128(2)^k; (k, log2(msb128(x - 2)), log2(msb128(x - 1)), log2(msb128(x)))) for k in 1:114]

# %%
function msb(x)
    ε = eps(float(typeof(x)))
    q = (1/ε + 1) * x
    m = q - (1 - ε/2) * q
    oftype(x, m)
end

# %%
@show precision(Float64)
@show eps(Float64) == 1/2^(precision(Float64)-1)
[(x = 2^k; (k, log2(msb(x - 2)), log2(msb(x - 1)), log2(msb(x)))) for k in 1:precision(Float64)+1]

# %%
@show precision(BigFloat)
@show eps(BigFloat) == 1/big(2)^(precision(BigFloat)-1)
[(x = big(2)^k; (k, log2(msb(x - 2)), log2(msb(x - 1)), log2(msb(x)))) for k in 1:precision(BigFloat)+1]

# %%
@code_typed msb(12345)

# %%
@code_native debuginfo=:none msb(12345)

# %%

max_stars_repo_path: 0019/msb trick.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0019/msb trick.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0019/msb trick.jl

max_forks_repo_name: genkuroki/public


text:
include("../quantum_espresso.jl")
include("../builder.jl")
config = Dict{String, Any}()
include("config.jl")

function main()
    cfg = config[@__DIR__]
    kwargs_scf = cfg.kwargs_scf
    mixing_beta = 0.7

    lattice, atoms = load_system(cfg.case, cfg.n_repeat)
    fileprefix = joinpath(cfg.directory, "TF/espresso")
    run_qe(lattice, atoms, fileprefix=fileprefix, mixing_beta=mixing_beta,
           temperature=kwargs_scf.temperature, kspacing=kwargs_scf.kspacing,
           smearing="gaussian", maxiter=50)
end


text:
# Prototype ParseStream interface
#
# Here we test the ParseStream interface, by taking input code and checking
# that the correct sequence of emit() and bump() produces a valid parse tree.

using JuliaSyntax: ParseStream,
    peek, peek_token,
    bump, bump_trivia, bump_invisible,
    emit, emit_diagnostic, TRIVIA_FLAG, INFIX_FLAG

code = """
for i = 1:10
    xx[i] + 2
    # hi
    yy
end
"""

st = ParseStream(code)

# Here we manually issue parse events in the order the Julia parser would issue
# them
@testset "ParseStream" begin
    p1 = position(st)
        @test peek(st) == K"for"
        bump(st, TRIVIA_FLAG)
        p2 = position(st)
            @test peek(st) == K"Identifier"    # 'i'
            bump(st)
            @test peek(st) == K"="
            bump(st, TRIVIA_FLAG)
            p3 = position(st)
                @test peek(st) == K"Integer"   # 1
                bump(st)
                @test peek(st) == K":"
                bump(st) # :
                @test peek(st) == K"Integer"   # 10
                bump(st) # 10
            emit(st, p3, K"call", INFIX_FLAG)
        emit(st, p2, K"=")
        @test peek(st) == K"NewlineWs"
        bump(st, TRIVIA_FLAG)
        p4 = position(st)
            p5 = position(st) # [call]
                p6 = position(st) # [ref]
                    @test peek(st) == K"Identifier" # 'xx'
                    bump(st)
                    @test peek(st) == K"["
                    bump(st, TRIVIA_FLAG)
                    @test peek(st) == K"Identifier" # 'i'
                    bump(st)
                    @test peek(st) == K"]"
                    bump(st, TRIVIA_FLAG)
                emit(st, p6, K"ref")
                @test peek(st) == K"+"
                bump(st)
                @test peek(st) == K"Integer"        # 2
                bump(st)
            emit(st, p5, K"call", INFIX_FLAG)
            @test peek(st) == K"NewlineWs"
            bump(st, TRIVIA_FLAG)
            @test peek(st) == K"NewlineWs"
            bump(st, TRIVIA_FLAG)
            @test peek(st) == K"Identifier" # 'yy'
            bump(st)
        emit(st, p4, K"block")
        @test peek(st) == K"NewlineWs"
        bump(st, TRIVIA_FLAG)
        bump(st, TRIVIA_FLAG) # end
    emit(st, p1, K"for")
    @test peek(st) == K"NewlineWs"
    bump(st, TRIVIA_FLAG)
    emit(st, p1, K"toplevel")
end

@test JuliaSyntax.build_tree(GreenNode, st) isa JuliaSyntax.GreenNode

# ## Input code
#=
println("-----------------------")
print(code)
println()

# ## Output tree
show(stdout, MIME"text/plain"(), t, code, show_trivia=true)
=#


text:
using Distributed

include("utils.jl")
include("pairwise_ui.jl")
include("advanced_ui.jl")
include("output_ui.jl")

const logo = joinpath(dirname(@__FILE__), "..", "assets", "cs_logo.ico")

function log_window()
    lw = node(:pre, "", id = "log",
              attributes = Dict(:style => "height: 900px; border: 2px solid Black;
                                width: 800px; border-style: solid;
                                overflow: auto"))

    s = Scope()
    s.dom = lw
    s["log"] = Observable("")
    s["clear"] = Observable(rand())
    onjs(s["log"], JSExpr.@js function (msg)
             JSExpr.@var el = this.dom.querySelector("#log")
             el.textContent += ("\n" + msg)
         end)
    onjs(s["clear"], JSExpr.@js function (msg)
             JSExpr.@var el = this.dom.querySelector("#log")
             el.textContent = ""
         end)
    s
end

function showsome(uis, which)
    s = Scope()
    s["visible"] = which
    s.dom = node(:div, id="cont", uis...)
    onjs(s["visible"], JSExpr.@js function (visbl)
                 JSExpr.@var cont = this.dom.querySelector("cont")
                 JSExpr.@var cs = cont.children
                 for i = 1:cs.length
                     if visbl.indexOf(i) >= 0
                         cs[i].style.display = "block"
                     else
                         cs[i].style.display = "none"
                     end
                 end
             end)
    s
end

function ui_logger(logging)
    function (msg, typ)
        logging["log"][] = msg
    end
end

function generate_ui(w)

    heading = node(:div, tachyons_css, "Circuitscape 5.0") |>
                    class"f2 lh-copy tc"

    picture = HTML("<img src = $(AssetRegistry.register(logo))>")

    section1 = node(:div, tachyons_css, "Data Type and Modelling Mode") |>
                    class"f4 lh-title"

    data_type = "Raster"
    scenario = "Pairwise"
    focal = Observable("")
    source = Observable("")
    ground = Observable("")
    points_input = Observable{Any}(Scope())
    num_parallel = Observable(1)
    solver_choice = "cg+amg"

    on(points_input) do x
        on(x["focal"]) do y
            focal[] = y
        end
        on(x["source"]) do z
            source[] = z
        end
        on(x["ground"]) do v
            ground[] = v
        end
    end
    points_input[] = pairwise_input_ui()

    # First drop down
    dt = get_data_type()

    # Next drop down
    mod_mode = map(dt["value"]) do v
        @show v
        data_type = v
        scenario = "Pairwise"
        points_input[] = pairwise_input_ui()
        if v == "Network"
            get_mod_mode_network()
        else
            get_mod_mode_raster()
        end
    end

    # Get the input raster/graph
    input_section = node(:div, tachyons_css, "Input Resistance Data") |>
                    class"f4 lh-title"
    input = input_ui()

    input_graph = Observable("")
    is_res = Observable(false)
    on(input["filepath"]) do x
        println("filepath changed to $x")
        input_graph[] = x
    end
    on(input["check"]) do x
       println("check changed to $x")
       is_res[] = x
    end

    # Focal points or advanced mode
    pair = pairwise_input_ui()
    adv = advanced_input_ui()

    on(mod_mode) do x
        v = x["value"]
        on(v) do s
            scenario = s
            if s == "Advanced"
                points_input[] = adv
                focal[] = ""
            else
                points_input[] = pair
            end
        end
    end
    dt["value"][] = "Raster"

    # Output options
    write_cur_maps = Observable(false)
    write_volt_maps = Observable(false)
    out = Observable("")
    output = output_ui()

    on(output["cur"]) do x
        write_cur_maps[] = x
    end
    on(output["volt"]) do x
        write_volt_maps[] = x
    end
    on(output["out"]) do x
        out[] = x
    end
    logging = log_window()
    Circuitscape.ui_interface[] = ui_logger(logging)

    o = Observable(1)
    parallel = spinbox(value = o)
    on(o) do x
        @show x
        num_parallel[] = x
    end

    o2 = Observable("cg+amg")
    solver = dropdown(["cg+amg", "cholmod"], value = o2)
    on(o2) do x
        @show x
        solver_choice = x
    end

    compute_section = node(:div, tachyons_css, "Compute Options") |>
                    class"f4 lh-title"
    parallel_option = vbox(node(:div, "Number of parallel processes to use: ",
                      attributes = Dict(:style => "margin-top: 12px")),
                 parallel)
    solver_option = vbox(node(:div, "Solver: ",
                      attributes = Dict(:style => "margin-top: 12px")),
                 solver)

    # Run button
    run, ob = run_button()
    on(ob) do x
        @show data_type
        @show scenario
       @show input_graph[]
       @show is_res[]
       @show focal[]
       @show source[]
       @show ground[]
       @show out[]
       @show write_cur_maps[]
       @show write_volt_maps[]
       @show num_parallel[]
       cfg = Dict{String,String}()
       cfg["data_type"] = data_type
       cfg["scenario"] = scenario
       cfg["habitat_file"] = input_graph[]
       cfg["habitat_map_is_resistances"] = string(is_res[])
       cfg["point_file"] = focal[]
       cfg["source_file"] = source[]
       cfg["ground_file"] = ground[]
       cfg["output_file"] = out[]
       cfg["write_cur_maps"] = string(write_cur_maps[])
       cfg["write_volt_maps"] = string(write_volt_maps[])
       if num_parallel[] > 1
           println("attempting to start up with $(num_parallel[]) processes")
           cfg["parallelize"] = "true"
           cfg["max_parallel"] = string(num_parallel[])
       end
       cfg["solver"] = solver_choice

       logging["clear"][] = rand()
       try
           compute(cfg)
           Blink.@js_ w alert("Run completed!")
       catch e
           @show e
           Blink.js(w, Blink.JSString("""alert("Error: $(e.msg)")"""), callback = false)
       end

    end

    runtests_prompt = """
    It is highly recommended you run the tests on startup.
    This warms up the app and checks whether it works as expected:
    """
    runtests = runtests_button(w)

    left = vbox(section1,
                dt,
                mod_mode,
                input_section,
                input,
                points_input,
                compute_section,
                parallel_option,
                solver_option,
                output,
                vskip(1em),
                run)

    right = vbox(node(:div, "Logging") |> class"f4 lh-title",
                 vskip(1em),
                 logging)

    page = vbox(alignself(:center,
    				hbox(picture,
    					hskip(3em),
    					heading,
    					hskip(3em),
    					vbox(vskip(0.5em), runtests))),
                hline(style = :solid, w=5px)(style = Dict(:margin => 20px)),
            	hbox(left,
                	hskip(0.5em),
                	vline(style = :solid, w=3px)(style = Dict(:margin => 10px)),
                	hskip(0.5em),
                	right)
            	)

    page = page |> class"pa3 system-sans-serif"

    # body!(w, page)

end

function run_ui()
    w = Window(Dict(:show=>false))
    wait(w.content)
    p = generate_ui(w)
   # tools(w)
    Blink.body!(w, p)
    w1 = Window(Dict(:title => "Circuitscape",
                    :width => 1300,
                    :height => 800,
                    :icon => logo))
    p = generate_ui(w1)
    wait(w1.content)
    Blink.body!(w1, p)
    w1
end


text:
mutable struct GaussianDist
    _mu::Float64
    _var::Float64 # sigma^2
    _precision::Float64 # 1/sigma^2
    function GaussianDist(mu::Float64, var::Float64)
        new(mu, var, 1.0 / var)
    end
end

mutable struct GaussianBayesMeanDist
    _mu_gauss::GaussianDist
    _var::Float64
    _precision::Float64 # 1/sigma^2
    function GaussianBayesMeanDist(mu_gauss::GaussianDist, var::Float64)
        new(mu_gauss, var, 1.0 / var)
    end
end

mutable struct GaussianBayesVarDist
    _var_gamma::GammaDist
end

function pdf(gauss::GaussianDist, x::Float64)
    d = x - gauss._mu
    return exp(-0.5 * d^2 / gauss._var) / sqrt(2 * pi * gauss._var)
end

function pdf(gauss::GaussianDist, X::AbstractArray{Float64,1})
    d = X .- gauss._mu
    return exp.(-0.5 * d .^ 2 / gauss._var) / sqrt(2 * pi * gauss._var)
end

function fitting(gauss::GaussianDist, X::AbstractArray{Float64,1})
    gauss._mu = mean(X)
    gauss._var = var(X)
    gauss._precision = 1.0 / gauss._var
end

function fitting(gauss::GaussianBayesMeanDist, X::AbstractArray{Float64,1})
    # mu and sigma for prior
    mu_0 = gauss._mu_gauss._mu
    precision_0 = gauss._mu_gauss._precision
    precision = gauss._precision # precision of ground-truth distribution

    N = size(X)[1]
    mu_ML = mean(X)

    precision_N = precision_0 + N * precision
    mu_N = (precision_0 / precision_N) * mu_0 + (N * precision / precision_N) * mu_ML

    gauss._mu_gauss._mu = mu_N
    gauss._mu_gauss._var = 1.0 / precision_N
    gauss._mu_gauss._precision = precision_N
end

function fitting(gauss::GaussianBayesVarDist, X::AbstractArray{Float64,1})
    sigma2_ML = var(X)
    N = size(X)[1]

    if N != 1
        gauss._var_gamma._a += N / 2.0
        gauss._var_gamma._b += N / 2.0 * sigma2_ML
    end
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
#----------------------------------------------------------------------------
"""
    OlsPs(y,x,ExciseIt=false,UnExciseIt=false,SkipCovbIt=false)

Do OLS of y on x, for one dependent variable or SURE with same regressors

# Usage
`(b,res,yhat,Covb,R2a,T) = OlsPs(y,x,ExciseIt,UnExciseIt,SkipCovbIt)`

# Input
- `y::Array`:         Tx1 or Txn matrix of the dependent variables
- `x::Array`:         Txk matrix of regressors (including deterministic ones)
- `ExciseIt::Bool`:   (optional) bool, true to doing excise([y,x])
- `UnExciseIt::Bool`: (optional) bool, true to put residuals and fitted values in Txn matrix
- `SkipCovbIt::Bool`: (optional) bool, true to skip calculation of Covb (sensitive inversion)

# Output
- `b::Array`:     kxn matrix, regression coefficients
- `res::Array`:   Tx1 or Txn matrix, residuals y - yhat
- `yhat::Array`:  Tx1 or Txn matrix, fitted values x*b
- `Covb::Array`:  matrix, covariance matrix of vec(b) = [beq1;beq2;...]
- `R2a::Array`:   1xn vector, R2 values
- `Ty::Int`:      scalar, number of obs (after excise)

# Requires
- FindNoNaNPs

"""
function OlsPs(y,x,ExciseIt=false,UnExciseIt=false,SkipCovbIt=false)

  (T,n) = (size(y,1),size(y,2))

  if ExciseIt
     vvNoNaNR = FindNNPs(y,x)                #find rows with no NaNs
    (y,x)     = (y[vvNoNaNR,:],x[vvNoNaNR,:])
  end

  (Ty,Tx,k) = (size(y,1),size(x,1),size(x,2))

  if Tx != Ty
    error("y and x must have same number of observations")
  end
  if any(isnan.([y x]))
    error("NaN in x or y")
  end
  if ndims(x) == 1               #T vector to (T,1) matrix, could also do reshape(x,T,1)
    x = x[:,:]
  end

  if Ty >= k
    b = x\y
    #b = (x'x)\(x'y)
    #ndims(x) == 1 ? yhat2  = x*b[1] : yhat2  = x*b  #if? then: else. to handle x is vector
    yhat2  = x*b
    res2   = y - yhat2
    Covres = cov(res2)*(Ty-1)/Ty
    if SkipCovbIt
      Covb = NaN
    else
      Covb   = kron(Covres,inv(x'x))
    end
    R2a = 1 .- var(res2,dims=1)./var(y,dims=1)
  else
    (b,yhat2,res2,Covb,R2a) = (NaN,NaN,NaN,NaN,NaN)
  end

  if UnExciseIt && ExciseIt   #put fitted value and residual in Tyxn matrix
    (yhat,res)       = [fill(NaN,(T,n)) for i=1:2]
    yhat[vvNoNaNR,:] = yhat2
    res[vvNoNaNR,:]  = res2
  else
    yhat = yhat2
    res  = res2
  end

  if n == 1                #if only one regression
    R2a = R2a[1]
  end

  return b,res,yhat,Covb,R2a,Ty

end
#------------------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
invalid_imgs = [
    ("too_few_dimensions", MethodError, rand(127)),
    ("too_many_channels", AssertionError, rand(127, 257, 5)),
    ("too_many_dimensions", MethodError, rand(127, 257, 3, 1)),
    ("palette_too_large", ArgumentError, IndirectArray(rand(1:257, 127, 256), rand(RGB{N0f8}, 257))),
    ("palette_index_too_few_dimensions", MethodError, IndirectArray(rand(1:4, 127), rand(RGB{N0f8}, 4))),
    ("palette_index_too_many_dimensions", ArgumentError, IndirectArray(rand(1:4, 127, 256, 1), rand(RGB{N0f8}, 4))),
    ("palette_eltype_too_wide", ArgumentError, IndirectArray(rand(1:4, 127, 256), rand(RGB{N0f16}, 4))),
    ("palette_non_RGB_colorant", ArgumentError, IndirectArray(rand(1:4, 127, 256), rand(Gray{N0f8}, 4))),
]

@testset "invalid inputs" begin
for (case, exception, image) in invalid_imgs
        @testset "$(case) throws" begin
            @test_throws exception PNGFiles.save(
                joinpath(PNG_TEST_PATH, "test_img_err_$(case).png"),
                image
            )
        end
    end
end


text:
# Autogenerated wrapper script for Dierckx_jll for x86_64-apple-darwin-libgfortran5
export libddierckx

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Dierckx")
JLLWrappers.@declare_library_product(libddierckx, "@rpath/libddierckx.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libddierckx,
        "lib/libddierckx.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BitTorrent, Test

@testset "Bencode/Bdecode" begin
    # BitTorrent.bencode
    @test sprint(BitTorrent.bencode, 42) == "i42e"
    @test sprint(BitTorrent.bencode, [42, [1, 2], Dict("hello"=>42), "hello world"]) ==
          "li42eli1ei2eed5:helloi42ee11:hello worlde"
    @test sprint(BitTorrent.bencode, Dict("int"=>42, "list"=>[1, 2], "dict"=>Dict("key1"=>1), "string"=>"hello world")) ==
          "d4:dictd4:key1i1ee3:inti42e4:listli1ei2ee6:string11:hello worlde"
    @test sprint(BitTorrent.bencode, "hello world") == "11:hello world"
    @test sprint(BitTorrent.bencode, "\x8e\xc8\xfb\xd3") == "4:\x8e\xc8\xfb\xd3"
    # BitTorrent.bdecode
    @test BitTorrent.bdecode(IOBuffer("i42e")) == 42
    @test BitTorrent.bdecode(IOBuffer("i42e")) == 42
end