text:
push!(LOAD_PATH, "src")

using amplitudes_compass
using DelimitedFiles
function getReIm(func)
    ev = LinRange(2mπ,2.5,500)
    cal = func.(ev.^2)
    return [ev real.(cal) imag.(cal)]
end

function writeEverything2txt(path_to_write = "/tmp")
    for (func, foutname) in zip([fρ,     ff2,    fρ3,      fσ, ff0_980, ff0_1500],
                                ["rho", "f2", "rho3", "pipiS", "f0980", "f01500"])
        writedlm(joinpath(path_to_write, "isobar_"*foutname*".txt"), getReIm(func))
    end
end
writeEverything2txt()

using Plots
pyplot()
function plotEverything()
    plot(layout=grid(1,2), size=(500,350))
    for f in [fρ, ff2, fρ3]
        mtrx = getReIm(f)
        plot!(mtrx[:,1], sqrt.(mtrx[:,2].^2+mtrx[:,3].^2), subplot=1, lab="")
    end
    for f in [fσ, ff0_980, ff0_1500]
        mtrx = getReIm(f)
        plot!(mtrx[:,1], sqrt.(mtrx[:,2].^2+mtrx[:,3].^2), subplot=2, lab="")
    end
    plot!()
end

plotEverything()


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

using Markdown
using InteractiveUtils

# ╔═╡ 660bbe15-74ae-4e61-9a7e-6b2bf1b1e1f5
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
	    Pkg.PackageSpec(name="DifferentialEquations", version="6"),
	    Pkg.PackageSpec(name="Parameters", version="0.12"),
        Pkg.PackageSpec(name="LabelledArrays", version="1"),
        Pkg.PackageSpec(name="Setfield", version="0.7"),
    ])
    using Plots, PlutoUI, LinearAlgebra, DifferentialEquations, Parameters, LabelledArrays, Setfield

	Plots.gr(fmt=:png, lw=2)

	PlutoUI.TableOfContents()
end

# ╔═╡ a65fcaa1-6d49-4121-b644-2e8945763fd2
md"""
# Figure 4.15, 4.16, and 4.17

Oscillatory network.
"""

# ╔═╡ 3fa4e3bb-2fa5-4f99-9e2d-e41936cd8e6d
"""
Model of oscillatory network from Figure 4.14. This code generates Figures
4.15, 4.16, and 4.17
"""
function model!(du, u, p, t)
    @unpack s1, s2 = u
    @unpack K0, K1, K2, N = p
    v0 = K0
    v1 = K1 * s1 * (1 + s2^N)
    v2 = K2 * s2
	du.s1 = v0 - v1
	du.s2 = v1 - v2
    return du
end

# ╔═╡ cf4047c0-bdec-11eb-0bcd-d958a6ba03d6
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ 7abc2f61-a1ed-4e61-bfa8-e557f944f615
begin
	# Convenience functions
    hill(x, k) = x / (x + k)
    hill(x, k, n) = hill(x^n, k^n)
    exprel(x) = ifelse(x≈zero(x), one(x), x / expm1(x))
end

# ╔═╡ c37327cf-a798-4820-8304-36966ac64e56
function figure0415(; param = (K0 = 8.0, K1 = 1.0, K2 = 5.0, N = 2),
	                  r = LinRange(0.0, 4.0, 20),
	                  tend = 8.0,
	                  figtitle="Fig 4.15")

	u0s = ( LVector(s1=1.5, s2=1.0), LVector(s1=0.0, s2=1.0),
        	LVector(s1=0.0, s2=3.0), LVector(s1=2.0, s2=0.0))
	sols = map(u0 -> solve(ODEProblem(model!, u0, tend, param)), u0s)

	# Fig 4.15 A
	p1 = plot(sols[1], xlabel="Time", ylabel="Concentration", title ="$figtitle (A)", xlims=(0.0, 8.0))

	# Fig 4.15 B: Vetor field
	function ∂F(x, y)
		u = LVector(s1=x, s2=y)
		dxdy = model!(similar(u), u, param, 0.0)
		return dxdy ./ (norm(dxdy)^0.5 * 20)
	end

	nullcline_s1(s2, p=param) = (p.K0 / p.K1) * hill(1, s2, p.N)
	nullcline_s2(s2, p=param) = (p.K2 * s2) / (p.K1 * (1 + s2^p.N))

	xx = [x for y in r, x in r]
	yy = [y for y in r, x in r]
	p2 = quiver(xx, yy, quiver=∂F, line=(:lightblue))

	for sol in sols
		plot!(p2, sol, vars=(1, 2), label=nothing)
	end

	rMin, rMax = r[1], r[end]

	plot!(p2, nullcline_s1, identity, rMin, rMax, label="Nullcline S1", line=(:dash, :red))
	plot!(p2, nullcline_s2, identity, rMin, rMax, label="Nullcline S2", line=(:dash, :blue))
	plot!(p2, title = "$figtitle (B)", xlabel="[S1]", ylabel="[S2]",
      xlims=(rMin, rMax), ylims=(rMin, rMax), aspect_ratio=:equal, size=(700, 700))

	return (p1, p2)
end

# ╔═╡ 497e9b54-40c1-493b-aabe-f6ad66120389
fig415a, fig415b = figure0415()

# ╔═╡ c72c49a8-c718-419b-aca6-31853282b150
fig415a

# ╔═╡ 934f6dc0-5c24-466e-b07c-879efbfd1f72
fig415b

# ╔═╡ 21027a60-2903-49d1-ba10-4c19f944fbb0
fig416a, fig416b = figure0415(param = (K0 = 8.0, K1 = 1.0, K2 = 5.0, N = 2.5), tend = 1000.0, figtitle="Fig 4.16")

# ╔═╡ 39d5557b-430b-4074-8944-94f21efcab7e
fig416a

# ╔═╡ 747493eb-9196-4c6b-b70c-42c6167c0b9f
fig416b

# ╔═╡ 6553c713-9663-4260-a416-5d0c8590ea97
let
	param = (K0 = 8.0, K1 = 1.0, K2 = 5.0, N = 2.5)
	sol417 = solve(ODEProblem(model!, LVector(s1=2.0, s2=1.5), 10.0, param))
	r = LinRange(0.0, 4.0, 20)
	xx = [x for y in r, x in r]
	yy = [y for y in r, x in r]

	"Vetor field"
	function ∂F(x, y)
		u = LVector(s1=x, s2=y)
		dxdy = model!(similar(u), u, param, 0.0)
		return dxdy ./ (norm(dxdy)^0.5 * 20)
	end

	nullcline_s1(s2, p=param) = (p.K0 / p.K1) * hill(1, s2, p.N)
	nullcline_s2(s2, p=param) = (p.K2 * s2) / (p.K1 * (1 + s2^p.N))

	xx = [x for y in r, x in r]
	yy = [y for y in r, x in r]
	p2 = quiver(xx, yy, quiver=∂F, line=(:lightblue))

	quiver(xx, yy, quiver=∂F, line=(:lightblue))
	plot!(sol417, vars=(1, 2), label=nothing, line=(:black), arrow=0.4)

	plot!(nullcline_s1, identity, 0.0, 4.0, label="Nullcline S1", line=(:dash, :red))
	plot!(nullcline_s2, identity, 0.0, 4.0, label="Nullcline S2", line=(:dash, :blue))
	plot!(title = "Fig. 4.17", xlabel="[S1]", ylabel="[S2]",
		  xlims=(1.0, 3.0), ylims=(1.0, 3.0), aspect_ratio=:equal, size=(700, 700))
end

# ╔═╡ Cell order:
# ╠═a65fcaa1-6d49-4121-b644-2e8945763fd2
# ╠═3fa4e3bb-2fa5-4f99-9e2d-e41936cd8e6d
# ╠═c37327cf-a798-4820-8304-36966ac64e56
# ╠═497e9b54-40c1-493b-aabe-f6ad66120389
# ╠═c72c49a8-c718-419b-aca6-31853282b150
# ╠═934f6dc0-5c24-466e-b07c-879efbfd1f72
# ╠═21027a60-2903-49d1-ba10-4c19f944fbb0
# ╠═39d5557b-430b-4074-8944-94f21efcab7e
# ╠═747493eb-9196-4c6b-b70c-42c6167c0b9f
# ╠═6553c713-9663-4260-a416-5d0c8590ea97
# ╠═cf4047c0-bdec-11eb-0bcd-d958a6ba03d6
# ╠═7abc2f61-a1ed-4e61-bfa8-e557f944f615
# ╠═660bbe15-74ae-4e61-9a7e-6b2bf1b1e1f5

ext: jl

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#integerchallenge.jl

include("../methods/unum_optimizer.jl")

using Unums

const limit = 100_000_000

function addtolimit(T::Type)
  sum = one(T)
  for idx = 1:limit
    sum += one(T)
  end
  #println(Unums.calculate(sum))
  sum
end

r = optimize(addtolimit, 0.1, verbose = true)
describe(r)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
repo = "bsuite"

const BSUITE_DATASETS = [
    "cartpole",
    "catch",
    "mountain_car"
]

datasets = [
    "full",
    "full_train",
    "full_valid"
]

function bsuite_init()
    for env in BSUITE_DATASETS
        for stochasticity in 0.0:0.1:0.5
            for index in 0:4
                for dataset in datasets
                    register(
                        DataDep(
                            "rl-unplugged-$repo-$env-$stochasticity-$index-$dataset",
                            """
                            Dataset: RL Unplugged bsuite
                            Credits: https://arxiv.org/abs/2006.13888
                            Url: https://github.com/deepmind/deepmind-research/tree/master/rl_unplugged
                            Authors: Caglar Gulcehre, Ziyu Wang, Alexander Novikov, Tom Le Paine,
                            Sergio Gómez Colmenarejo, Konrad Zolna, Rishabh Agarwal,
                            Josh Merel, Daniel Mankowitz, Cosmin Paduraru, Gabriel
                            Dulac-Arnold, Jerry Li, Mohammad Norouzi, Matt Hoffman,
                            Ofir Nachum, George Tucker, Nicolas Heess, Nando deFreitas
                            Year: 2020

                            Data accompanying [RL Unplugged: Benchmarks for Offline Reinforcement Learning].
                            bsuite datasets are fairly light-weight and running experiments doesn't require
                            too much compute. We recommend you to try bsuite, if you are interested in
                            small-scale and easy to run offline RL datasets generated by stochastic environments
                            where the stochasticity of the environment is easy to control.
                            """,
                            "gs://rl_unplugged/$repo/$env/0_$stochasticity/$(index)_$dataset-00000-of-00001",
                            fetch_method = fetch_gc_file
                        )
                    )
                end
            end
        end
    end
end


text:
using DiffusionModels

μ = ConstDrift(1.2)
σ = ConstSigma(1)
bounds = ConstSymBounds(2,-2)
ndt = NonDecisionTime(lower = 2, upper = 3)
dm = DiffusionModel(μ, σ, bounds, ndt, Δt = 0.01)

upper, lower = fpt(dm, tmax = 10.0)
choice, rt = rand(dm, 100)

# or using the basic type constructor:

drift = [1.2]
sig = [1.0]
bound_lo = [-1.1]
bound_hi = [1.1]
ndt = (2.0, 3.0)

upper, lower = ddm_fpt(drift, sig, bound_lo, bound_hi,
                        Δt = 0.01, tmax = 10.0)

# with non-decision time
ndt = (3.0, 6.0)
choice, rt = ddm_rand(drift, sig, bound_lo, bound_hi, ndt, Δt = 0.01, n = 1000, seed = 123)

# without non-decision time
choice, rt = ddm_rand(drift, sig, bound_lo, bound_hi, Δt = 0.01, n = 1000, seed = 123)


text:
using Documenter, BioProfiling

makedocs(sitename="BioProfiling.jl",
		 pages = [
	        "BioProfiling.jl" => "index.md",
	        "API" => "api.md",
	        "Examples" => "examples.md"
         ],
         authors = "Loan Vulliard" )

deploydocs(
    repo = "github.com/menchelab/BioProfiling.jl.git",
)


text:
include("euler/euler.jl")
using .Geometry: triangle_area
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(path::String)::Int
    triangles = map(x -> parse.(Int, split(x, ",")), split(read(path, String), "\n"))
    x = [0, 0]
    result = 0
    for △ ∈ triangles
        a, b, c = △[1:2], △[3:4], △[5:6]
        if triangle_area(a, b, x) + triangle_area(a, x, c) + triangle_area(x, b, c) == triangle_area(a, b, c) result += 1 end
    end
    return result
end

compute("problems/0102/p102_triangles.txt")

@benchmark compute("problems/0102/p102_triangles.txt")

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using ADCME

PWD = pwd()
change_directory()
http_file("https://github.com/kailaix/hypre/archive/v2.19.0.tar.gz", "v2.19.0.tar.gz")
uncompress("v2.19.0.tar.gz", "hypre-2.19.0")
change_directory("hypre-2.19.0/src/build")
ROOT = joinpath(ADCME.BINDIR, "..")
run_with_env(`$(ADCME.CMAKE) -G Ninja -DCMAKE_MAKE_PROGRAM=$(ADCME.NINJA)
            -DHYPRE_SHARED:BOOL=ON -DHYPRE_INSTALL_PREFIX:PATH=$ROOT
            -DCMAKE_C_COMPILER:FILEPATH=$(ADCME.CC) -DCMAKE_CXX_COMPILER:FILEPATH=$(ADCME.CXX) ..`)
run_with_env(`$(ADCME.CMAKE) -G Ninja -DCMAKE_MAKE_PROGRAM=$(ADCME.NINJA) -L ..`)
ADCME.make()
run_with_env(`$(ADCME.NINJA) install`)
# run_with_env(`mv $(ROOT)/lib64/libHYPRE.so $(ADCME.LIBDIR)/libHYPRE.so`)
cd(PWD)


text:
module DashBootstrapComponents

using DashBase
import Pkg
using Pkg.Artifacts

DashBase.@generate_wrapper()
@generate_components

end


text:
using Interpolations
using SignalAnalysis
using DSP: hilbert

export IsoSSP, MunkSSP, SampledSSP, ConstantDepth, SampledDepth, SampledAltitude
export ReflectionCoef, FlatSurface, RayleighReflectionCoef, SurfaceLoss
export Rock, Pebbles, SandyGravel, CoarseSand, MediumSand, FineSand, VeryFineSand
export ClayeySand, CoarseSilt, SandySilt, Silt, FineSilt, SandyClay, SiltyClay, Clay
export Vacuum, SeaState0, SeaState1, SeaState2, SeaState3, SeaState4
export SeaState5, SeaState6, SeaState7, SeaState8, SeaState9
export AcousticReceiverGrid2D, AcousticReceiverGrid3D, NarrowbandAcousticSource
export RedGaussianNoise, Pinger, SampledAcousticSource

### sound speed profiles

"""
$(TYPEDEF)
Isovelocity sound speed profile.

---

    IsoSSP(c)

Create an isovelocity sound speed profile with sound speed `c`.
"""
struct IsoSSP{T} <: SoundSpeedProfile
  c::T
end

soundspeed(ssp::IsoSSP, x, y, z) = ssp.c

"""
$(TYPEDEF)
Munk sound speed profile.
"""
struct MunkSSP <: SoundSpeedProfile end

function soundspeed(::MunkSSP, x, y, z)
  ϵ = 0.00737
  z̃ = 2.0 * (-z - 1300.0) / 1300.0
  1500.0 * (1.0 + ϵ * (z̃ - 1 + exp(-z̃)))
end

"""
$(TYPEDEF)
Sound speed profile based on measurements at discrete depths.
"""
struct SampledSSP{T1,T2,T3} <: SoundSpeedProfile
  z::Vector{T1}
  c::Vector{T2}
  interp::Symbol
  f::T3
  function SampledSSP(depth, c, interp)
    if interp === :smooth
      depth isa AbstractRange || throw(ArgumentError("depth must be sampled uniformly and specified as an AbstractRange"))
      f = CubicSplineInterpolation(depth, c; extrapolation_bc=Line())
    elseif interp === :linear
      f = LinearInterpolation(depth, c; extrapolation_bc=Line())
    else
      throw(ArgumentError("Unknown interpolation"))
    end
    new{eltype(depth),eltype(c),typeof(f)}(-depth, c, interp, f)
  end
end

"""
    SampledSSP(depth, c)
    SampledSSP(depth, c, interp)

Create a sound speed profile based on measurements at discrete depths.
`interp` may be either `:linear` or `:smooth`, and defaults to `:linear`
if unspecified.
"""
SampledSSP(depth, c) = SampledSSP(depth, c, :linear)

soundspeed(ssp::SampledSSP, x, y, z) = ssp.f(-z)

struct SampledSSPRangeDependent{T1,T2,T3,T4} <: SoundSpeedProfile
  z::Vector{T1}
  x::Vector{T2}
  c::Vector{T3}
  interp::Symbol
  f::T4
  function SampledSSPRangeDependent(depth, x, c, interp)
    if interp === :smooth
      depth isa AbstractRange || throw(ArgumentError("depth must be sampled uniformly and specified as an AbstractRange"))
      x isa AbstractRange || throw(ArgumentError("x range must be sampled uniformly and specified as an AbstractRange"))
      f = CubicSplineInterpolation((depth,x), c; extrapolation_bc=Line())
    elseif interp === :linear
      f = LinearInterpolation((depth,x), c; extrapolation_bc=Line())
    else
      throw(ArgumentError("Unknown interpolation"))
    end
    new{eltype(depth),eltype(x),eltype(c),typeof(f)}(-depth,x, c, interp, f)
  end
end

"""
  SampledSSPRangeDependent(depth, x, c)
  SampledSSPRangeDependent(depth, x, c, interp)

Create a sound speed profile based on measurements at discrete depths.
`interp` may be either `:linear` or `:smooth`, and defaults to `:linear`
if unspecified.
"""
SampledSSPRangeDependent(depth, x, c) = SampledSSPRangeDependent(depth, x, c, :linear)

soundspeed(ssp::SampledSSPRangeDependent, x, y, z) = ssp.f((-z,x))

function Base.show(io::IO, ssp::SampledSSPRangeDependent{T1,T2,T3,T4}) where {T1,T2,T3,T4}
  print(io, "SampledSSPRangeDependent{", T1, ",", T2, ",", T3, ",", ssp.interp, "}(", length(ssp.z), " points)")
end

### bathymetry models

"""
$(TYPEDEF)
Bathymetry with constant depth.

---

    ConstantDepth(depth)

Create a constant depth bathymetry.
"""
struct ConstantDepth{T} <: Bathymetry
  depth::T
end

depth(bathy::ConstantDepth, x, y) = bathy.depth
maxdepth(bathy::ConstantDepth) = bathy.depth

"""
$(TYPEDEF)
Bathymetry based on depth samples.
"""
struct SampledDepth{T1,T2,T3} <: Bathymetry
  x::Vector{T1}
  depth::Vector{T2}
  interp::Symbol
  f::T3
  function SampledDepth(x, depth, interp)
    if interp === :smooth
      x isa AbstractRange || throw(ArgumentError("x must be sampled uniformly and specified as an AbstractRange"))
      f = CubicSplineInterpolation(x, depth; extrapolation_bc=Line())
    elseif interp === :linear
      f = LinearInterpolation(x, depth; extrapolation_bc=Line())
    else
      throw(ArgumentError("Unknown interpolation"))
    end
    new{eltype(x),eltype(depth),typeof(f)}(x, depth, interp, f)
  end
end

"""
    SampledDepth(x, depth)
    SampledDepth(x, depth, interp)

Create a bathymetry given discrete depth measurements at locations given in `x`.
`interp` may be either `:linear` or `:smooth`, and defaults to `:linear`
if unspecified.
"""
SampledDepth(x, depth) = SampledDepth(x, depth, :linear)

depth(bathy::SampledDepth, x, y) = bathy.f(x)
maxdepth(bathy::SampledDepth) = maximum(bathy.depth)

function Base.show(io::IO, b::SampledDepth{T1,T2,T3}) where {T1,T2,T3}
  print(io, "SampledDepth{", T1, ",", T2, ",", b.interp, "}(", length(b.x), " points)")
end

### altimetry models

"""
$(TYPEDEF)
Altimetry for a flat surface with constant altitude of zero.
"""
struct FlatSurface <: Altimetry end

altitude(::FlatSurface, x, y) = 0.0

"""
$(TYPEDEF)
Altimetry based on altitude samples.
"""
struct SampledAltitude{T1,T2,T3} <: Altimetry
  x::Vector{T1}
  altitude::Vector{T2}
  interp::Symbol
  f::T3
  function SampledAltitude(x, altitude, interp)
    if interp === :smooth
      x isa AbstractRange || throw(ArgumentError("x must be sampled uniformly and specified as an AbstractRange"))
      f = CubicSplineInterpolation(x, altitude; extrapolation_bc=Line())
    elseif interp === :linear
      f = LinearInterpolation(x, altitude; extrapolation_bc=Line())
    else
      throw(ArgumentError("Unknown interpolation"))
    end
    new{eltype(x),eltype(altitude),typeof(f)}(x, altitude, interp, f)
  end
end

"""
    SampledAltitude(x, altitude)
    SampledAltitude(x, altitude, interp)

Create an altimetry given discrete altitude measurements at locations given in `x`.
`interp` may be either `:linear` or `:smooth`, and defaults to `:linear`
if unspecified.
"""
SampledAltitude(x, altitude) = SampledAltitude(x, altitude, :linear)

altitude(a::SampledAltitude, x, y) = a.f(x)

function Base.show(io::IO, a::SampledAltitude{T1,T2,T3}) where {T1,T2,T3}
  print(io, "SampledAltitude{", T1, ",", T2, ",", a.interp, "}(", length(a.x), " points)")
end

### reflection models

"""
$(TYPEDEF)
Reflection model for a surface with a constant reflection coefficient.

---

    ReflectionCoef(coef)

Create a reflection model for a surface with a constant reflection coefficient `coef`.
"""
struct ReflectionCoef{T<:Number} <: ReflectionModel
  coef::T
end

reflectioncoef(rm::ReflectionCoef, f, θ) = rm.coef

"""
$(TYPEDEF)
Reflection model for a surface with a Rayleigh reflection coefficient.

---

    RayleighReflectionCoef(ρᵣ, cᵣ, δ)
    RayleighReflectionCoef(ρᵣ, cᵣ)

Create a reflection model for a surface with a Rayleigh reflection coefficient
with relative density `ρᵣ`, relative sound speed `cᵣ`, and dimensionless
attentuation `δ`. If attentuation `δ` is unspecified, it is assumed to be zero.

See `reflectioncoef()` for more details.
"""
struct RayleighReflectionCoef{T1,T2,T3} <: ReflectionModel
  ρᵣ::T1
  cᵣ::T2
  δ::T3
end

RayleighReflectionCoef(ρᵣ, cᵣ) = RayleighReflectionCoef(ρᵣ, cᵣ, 0.0)

# from APL-UW TR 9407 (1994), IV-6 Table 2
const Rock = RayleighReflectionCoef(2.5, 2.5, 0.01374)
const Pebbles = RayleighReflectionCoef(2.5, 1.8, 0.01374)
const SandyGravel = RayleighReflectionCoef(2.492, 1.3370, 0.01705)
const CoarseSand = RayleighReflectionCoef(2.231, 1.2503, 0.01638)
const MediumSand = RayleighReflectionCoef(1.845, 1.1782, 0.01624)
const FineSand = RayleighReflectionCoef(1.451, 1.1073, 0.01602)
const VeryFineSand = RayleighReflectionCoef(1.268, 1.0568, 0.01875)
const ClayeySand = RayleighReflectionCoef(1.224, 1.0364, 0.02019)
const CoarseSilt = RayleighReflectionCoef(1.195, 1.0179, 0.02158)
const SandySilt = RayleighReflectionCoef(1.169, 0.9999, 0.01261)
const Silt = RayleighReflectionCoef(1.149, 0.9873, 0.00386)
const FineSilt = RayleighReflectionCoef(1.148, 0.9861, 0.00306)
const SandyClay = RayleighReflectionCoef(1.147, 0.9849, 0.00242)
const SiltyClay = RayleighReflectionCoef(1.146, 0.9824, 0.00163)
const Clay = RayleighReflectionCoef(1.145, 0.98, 0.00148)

reflectioncoef(rm::RayleighReflectionCoef, f, θ) = reflectioncoef(θ, rm.ρᵣ, rm.cᵣ, rm.δ)

"""
$(TYPEDEF)
Reflection model for a water surface affected by wind.

---

    SurfaceLoss(windspeed)

Create a reflection model for a surface affected by wind. `windspeed` is
given in m/s.
"""
struct SurfaceLoss{T} <: ReflectionModel
  windspeed::T
end

reflectioncoef(rm::SurfaceLoss, f, θ) = complex(-surfaceloss(rm.windspeed, f, θ), 0.0)

const Vacuum = ReflectionCoef(complex(-1.0, 0.0))

# WMO sea states
# from APL-UW TR 9407 (1994), II-4 Table 2 (median windspeed)
const SeaState0 = SurfaceLoss(0.8)
const SeaState1 = SurfaceLoss(2.6)
const SeaState2 = SurfaceLoss(4.4)
const SeaState3 = SurfaceLoss(6.9)
const SeaState4 = SurfaceLoss(9.8)
const SeaState5 = SurfaceLoss(12.6)
const SeaState6 = SurfaceLoss(19.3)
const SeaState7 = SurfaceLoss(26.5)
const SeaState8 = SurfaceLoss(30.6)
const SeaState9 = SurfaceLoss(32.9)

### basic environmental model

Base.@kwdef struct BasicUnderwaterEnvironment{T1<:Altimetry, T2<:Bathymetry, T3<:SoundSpeedProfile, T4<:Number, T5<:ReflectionModel, T6<:ReflectionModel, T7} <: UnderwaterEnvironment
  altimetry::T1 = FlatSurface()
  bathymetry::T2 = ConstantDepth(20.0)
  ssp::T3 = IsoSSP(soundspeed())
  salinity::T4 = 35.0
  seasurface::T5 = SeaState1
  seabed::T6 = SandySilt
  noise::T7 = RedGaussianNoise(db2amp(120.0))
end

"""
    UnderwaterEnvironment(; altimetry, bathymetry, ssp, salinity, seasurface, seabed, noise)

Create an underwater environment.
"""
UnderwaterEnvironment(; kwargs...) = BasicUnderwaterEnvironment(; kwargs...)

altimetry(env::BasicUnderwaterEnvironment) = env.altimetry
bathymetry(env::BasicUnderwaterEnvironment) = env.bathymetry
ssp(env::BasicUnderwaterEnvironment) = env.ssp
salinity(env::BasicUnderwaterEnvironment) = env.salinity
seasurface(env::BasicUnderwaterEnvironment) = env.seasurface
seabed(env::BasicUnderwaterEnvironment) = env.seabed
noise(env::BasicUnderwaterEnvironment) = env.noise

### noise models

"""
$(TYPEDEF)
Ambient noise model with Gaussian noise with a `1/f²` power spectral density.

---

    RedGaussianNoise(σ)

Create an ambient noise model with variance `σ²` and `1/f²` power spectral density.
"""
struct RedGaussianNoise{T} <: NoiseModel
  σ::T
end

function record(noisemodel::RedGaussianNoise, duration, fs; start=0.0)
  analytic(signal(rand(RedGaussian(; n=round(Int, duration*fs), σ=noisemodel.σ)), fs))
end

### basic source & recevier models

"""
$(TYPEDEF)
Narrowband acoustic source.
"""
struct NarrowbandAcousticSource{T1,T2,T3,T4} <: AcousticSource
  pos::NTuple{3,T1}
  f::T2
  A::T3
  ϕ::T4
end

"""
$(SIGNATURES)
Create a narrowband acoustic source with frequency `f` Hz at location (`x`, `y`, `z`).
The `sourcelevel` is in µPa @ 1m. A phase `ϕ` may be optionlly specified.
"""
NarrowbandAcousticSource(x, y, z, f; sourcelevel=db2amp(180.0), ϕ=0.0) = NarrowbandAcousticSource(promote(x, y, z), f, sourcelevel, ϕ)

"""
$(SIGNATURES)
Create a narrowband acoustic source with frequency `f` Hz at location (`x`, z`).
The `sourcelevel` is in µPa @ 1m. A phase `ϕ` may be optionlly specified.
"""
NarrowbandAcousticSource(x, z, f; sourcelevel=db2amp(180.0), ϕ=0.0) = NarrowbandAcousticSource(promote(x, 0, z), f, sourcelevel, ϕ)

"""
$(SIGNATURES)
Create a narrowband acoustic source with frequency `f` Hz at location (`x`, `y`, `z`).
The `sourcelevel` is in µPa @ 1m. A phase `ϕ` may be optionlly specified.
"""
AcousticSource(x, y, z, f; sourcelevel=db2amp(180.0), ϕ=0.0) = NarrowbandAcousticSource(promote(x, y, z), f, sourcelevel, ϕ)

"""
$(SIGNATURES)
Create a narrowband acoustic source with frequency `f` Hz at location (`x`, `y`, `z`).
The `sourcelevel` is in µPa @ 1m. A phase `ϕ` may be optionlly specified.
"""
AcousticSource(x, z, f; sourcelevel=db2amp(180.0), ϕ=0.0) = NarrowbandAcousticSource(promote(x, 0, z), f, sourcelevel, ϕ)

location(tx::NarrowbandAcousticSource) = tx.pos
nominalfrequency(tx::NarrowbandAcousticSource) = tx.f
phasor(tx::NarrowbandAcousticSource) = tx.A * cis(tx.ϕ)
record(tx::NarrowbandAcousticSource, duration, fs; start=0.0) = signal(tx.A .* cis.(2π .* tx.f .* (start:1/fs:start+duration-1/fs) .+ tx.ϕ), fs)

"""
$(TYPEDEF)
Narrowband pulsed acoustic source.
"""
Base.@kwdef struct Pinger{T1,T2,T3,T4,T5,T6,T7,T8} <: AcousticSource
  pos::NTuple{3,T1}
  frequency::T2
  sourcelevel::T3 = db2amp(180.0)
  phase::T4 = 0.0
  duration::T5 = 0.02
  start::T6 = 0.0
  interval::T7 = 1.0
  window::T8 = nothing
end

"""
    Pinger(x, y, z, f; sourcelevel, phase, duration, start, interval, window)

Create a pulsed narrowband acoustic source with frequency `f` Hz at location (`x`, `y`, `z`).
Additional parameters that may be specified:
- `sourcelevel` in µPa @ 1m (default 180 dB)
- `phase` of the narrowband signal (default 0)
- `duration` of the pulse in seconds (default 20 ms)
- `start` time of one of the pulses in seconds (default 0)
- pulse repetition `interval` in seconds (default 1 second)
- `window` type (from `DSP.jl`) (default `nothing`)
"""
Pinger(x, y, z, f; kwargs...) = Pinger(; pos=promote(x, y, z), frequency=f, kwargs...)

"""
    Pinger(x, z, f; sourcelevel, phase, duration, start, interval, window)

Create a pulsed narrowband acoustic source with frequency `f` Hz at location (`x`, `z`).
Additional parameters that may be specified:
- `sourcelevel` in µPa @ 1m (default 180 dB)
- `phase` of the narrowband signal (default 0)
- `duration` of the pulse in seconds (default 20 ms)
- `start` time of one of the pulses in seconds (default 0)
- pulse repetition `interval` in seconds (default 1 second)
- `window` type (from `DSP.jl`) (default `nothing`)
"""
Pinger(x, z, f; kwargs...) = Pinger(; pos=promote(x, 0, z), frequency=f, kwargs...)

location(tx::Pinger) = tx.pos
nominalfrequency(tx::Pinger) = tx.frequency
phasor(tx::Pinger) = tx.sourcelevel * cis(tx.phase)

function record(pinger::Pinger, duration, fs; start=0.0)
  nsamples = round(Int, duration * fs)
  ping = cw(pinger.frequency, pinger.duration, fs; phase=pinger.phase, window=pinger.window)
  x = zeros(eltype(ping), nsamples)
  k1 = ceil(Int, (start - pinger.start - pinger.duration) / pinger.interval)
  k2 = floor(Int, (start - pinger.start + duration) / pinger.interval)
  n = length(ping)
  for k ∈ k1:k2
    ndx = round(Int, (pinger.start + k * pinger.interval - start) * fs) + 1
    if ndx > 0 && ndx + n ≤ nsamples
      x[ndx:ndx+n-1] .= ping
    elseif ndx > 0
      x[ndx:end] .= ping[1:nsamples-ndx+1]
    elseif ndx + n ≤ nsamples
      x[1:n+ndx-1] .= ping[2-ndx:end]
    else
      x .= ping[2-ndx:1-ndx+nsamples]
    end
  end
  signal(pinger.sourcelevel * x, fs)
end

"""
$(TYPEDEF)
Acoustic source transmitting a sampled signal.
"""
struct SampledAcousticSource{T1,T2,T3<:AbstractArray} <: AcousticSource
  pos::NTuple{3,T1}
  frequency::T2
  signal::T3
end

"""
    SampledAcousticSource(x, y, z, sig; fs, frequency)

Create a sampled acoustic source transmitting signal `sig` at location (`x`, `y`, `z`). The
samples are assumed to be in µPa @ 1m from the source.

Additional parameters that may be specified:
- `fs` is the sampling rate of the signal `sig`
- nominal `frequency` in Hz (`nothing` to auto-estimate)
"""
function SampledAcousticSource(x, y, z, sig::AbstractArray; fs=framerate(sig), frequency=nothing)
  s = signal(sig, fs)
  SampledAcousticSource(promote(x, y, z), frequency === nothing ? meanfrequency(s) : frequency, s)
end

"""
    SampledAcousticSource(x, z, sig; fs, frequency)

Create a sampled acoustic source transmitting signal `sig` at location (`x`, `z`). The
samples are assumed to be in µPa @ 1m from the source.

Additional parameters that may be specified:
- `fs` is the sampling rate of signal `sig`
- nominal `frequency` in Hz (`nothing` to auto-estimate)
"""
function SampledAcousticSource(x, z, sig::AbstractArray; fs=framerate(sig), frequency=nothing)
  s = signal(sig, fs)
  SampledAcousticSource(promote(x, 0, z), frequency === nothing ? meanfrequency(s) : frequency, s)
end

location(tx::SampledAcousticSource) = tx.pos
nominalfrequency(tx::SampledAcousticSource) = tx.frequency
phasor(tx::SampledAcousticSource) = throw(ArgumentError("phasor is undefined for a SampledAcousticSource"))

function record(tx::SampledAcousticSource, duration, fs=framerate(tx.signal); start=0.0)
  fs == framerate(tx.signal) || throw(ArgumentError("Resampling of SampledAcousticSource is not supported"))
  t1 = toframe(start, tx.signal)
  t2 = toframe(start+duration, tx.signal) - 1
  y = padded(tx.signal, (max(t1-1, 0), max(t2-t1+1-length(tx.signal), 0)); delay=-t1+1)
  @view y[1:t2-t1+1]
end

"""
$(TYPEDEF)
Omnidirectional acoustic receiver.
"""
struct BasicAcousticReceiver{T} <: AcousticReceiver
  pos::NTuple{3,T}
end

"""
$(SIGNATURES)
Create an omnidirectional acoustic receiver at location (`x`, `y`, `z`).
"""
BasicAcousticReceiver(x, y, z) = BasicAcousticReceiver(promote(x, y, z))

"""
$(SIGNATURES)
Create an omnidirectional acoustic receiver at location (`x`, `z`).
"""
BasicAcousticReceiver(x, z) = BasicAcousticReceiver(promote(x, 0, z))

"""
$(SIGNATURES)
Create an omnidirectional acoustic receiver at location (`x`, `y`, `z`).
"""
AcousticReceiver(x, y, z) = BasicAcousticReceiver(promote(x, y, z))

"""
$(SIGNATURES)
Create an omnidirectional acoustic receiver at location (`x`, `z`).
"""
AcousticReceiver(x, z) = BasicAcousticReceiver(promote(x, 0, z))

location(rx::BasicAcousticReceiver) = rx.pos

### receiver grids for transmission loss computation

"""
$(TYPEDEF)
A 2D Cartesian grid of omnidirectional acoustic receivers.
"""
struct AcousticReceiverGrid2D{T} <: AbstractArray{BasicAcousticReceiver{T},2}
  xrange::StepRangeLen{T,T,T}
  zrange::StepRangeLen{T,T,T}
end

"""
$(SIGNATURES)
Create a 2D Cartesian grid of omnidirectional acoustic receivers with `nx` × `nz`
receviers starting (`xmin`, `zmin`) with step sizes `xstep` and `zstep`.
"""
function AcousticReceiverGrid2D(xmin, xstep, nx, zmin, zstep, nz)
  xmin, xstep, zmin, zstep = promote(xmin, xstep, zmin, zstep)
  AcousticReceiverGrid2D(StepRangeLen(xmin, xstep, nx), StepRangeLen(zmin, zstep, nz))
end

Base.size(g::AcousticReceiverGrid2D) = (g.xrange.len, g.zrange.len)
Base.getindex(g::AcousticReceiverGrid2D, I::Vararg{Int,2}) = AcousticReceiver(g.xrange[I[1]], g.zrange[I[2]])
Base.setindex!(g::AcousticReceiverGrid2D, v, I::Vararg{Int,2}) = throw(ArgumentError("AcousticReceiverGrid2D is readonly"))

"""
$(TYPEDEF)
A 3D Cartesian grid of omnidirectional acoustic receivers.
"""
struct AcousticReceiverGrid3D{T} <: AbstractArray{BasicAcousticReceiver{T},3}
  xrange::StepRangeLen{T,T,T}
  yrange::StepRangeLen{T,T,T}
  zrange::StepRangeLen{T,T,T}
end

"""
$(SIGNATURES)
Create a 3D Cartesian grid of omnidirectional acoustic receivers with `nx` × `ny` × `nz`
receviers starting (`xmin`, `ymin`, `zmin`) with step sizes `xstep`, `ystep`, and `zstep`.
"""
function AcousticReceiverGrid3D(xmin, xstep, nx, ymin, ystep, ny, zmin, zstep, nz)
  xmin, xstep, ymin, ystep, zmin, zstep = promote(xmin, xstep, ymin, ystep, zmin, zstep)
  AcousticReceiverGrid3D(StepRangeLen(xmin, xstep, nx), StepRangeLen(ymin, ystep, ny), StepRangeLen(zmin, zstep, nz))
end

Base.size(g::AcousticReceiverGrid3D) = (g.xrange.len, g.yrange.len, g.zrange.len)
Base.getindex(g::AcousticReceiverGrid3D, I::Vararg{Int,3}) = AcousticReceiver(g.xrange[I[1]], g.yrange[I[2]], g.zrange[I[3]])
Base.setindex!(g::AcousticReceiverGrid3D, v, I::Vararg{Int,3}) = throw(ArgumentError("AcousticReceiverGrid3D is readonly"))


text:
using SchroedingersSmoke, CLArrays
using Colors, GPUArrays

vol_size = (4,2,2)# box size
dims = (64, 32, 32) .* 2 # volume resolution
hbar = 0.1f0      # Planck constant
dt = 1f0/48f0     # time step

jet_velocity = (1f0, 0f0, 0f0)
nozzle_cen = Float32.((2-1.7, 1-0.034, 1+0.066))
nozzle_len = 0.5f0
nozzle_rad = 0.5f0
n_particles = 1000   # number of particles

ArrayType = CLArray

isf2 = ISF{ArrayType, UInt32, Float32}(vol_size, dims, hbar, dt);

# initialize psi
psi = ArrayType([(one(Complex64), one(Complex64) * 0.01f0) for i=1:dims[1], j=1:dims[2], k=1:dims[3]]);

psi .= normalize_psi.(psi);

kvec = jet_velocity ./ hbar;
omega = sum(jet_velocity.^2f0) / (2f0*hbar);
isjetarr = isjet.(isf2.positions, (nozzle_cen,), (nozzle_len,), (nozzle_rad,))
# constrain velocity
for iter = 1:10
    restrict_velocity!(isf2, psi, kvec, isjetarr, 0f0)
    pressure_project!(isf2, psi)
end

particles = ArrayType(map(x-> (0f0, 0f0, 0f0), 1:(10^6) * 3))

add_particles!(particles, 1:n_particles, nozzle_cen, nozzle_rad)

using GLVisualize;
w = glscreen(color = RGBA(0f0, 0f0, 0f0, 0f0));
@async renderloop(w)
particle_vis = visualize(
    (GeometryTypes.Circle(GeometryTypes.Point2f0(0), 0.002f0), GeometryTypes.Point3f0.(Array(particles))),
    boundingbox = nothing, # don't waste time on bb computation
    color = fill(RGBA{Float32}(0, 0, 0, 0.09), length(particles)),
    billboard = true
).children[]
_view(particle_vis, camera = :perspective)
particle_vis[:color][1:n_particles] = map(1:n_particles) do i
    xx = (i / n_particles) * 2pi
    RGBA{Float32}((sin(xx) + 1) / 2, (cos(xx) + 1.0) / 2.0, 0.0, 0.1)
end
function simloop(
        N, isf, psi, kvec, omega, n_particles, isjetarr,
        nozzle_rad, nozzle_cen, particle, particle_vis
    )
    dt = isf.dt; d = isf.d
    for iter = 1:N
        isopen(w) || break
        t = iter * dt
        # incompressible Schroedinger flow
        schroedinger_flow!(isf, psi)
        psi .= normalize_psi.(psi)
        pressure_project!(isf, psi)

        start = mod((iter - 1) * n_particles + 1, length(particle))
        stop = start + n_particles - 1
        add_particles!(particle, start:stop, nozzle_cen, nozzle_rad)
        particle_vis[:color][start:stop] = map(1:n_particles) do i
            xx = (i / n_particles) * 2pi
            RGBA{Float32}((sin(xx) + 1) / 2, (cos(xx) + 1.0) / 2.0, iter / N, 0.1)
        end
        # constrain velocity
        restrict_velocity!(isf, psi, kvec, isjetarr, omega*t)
        pressure_project!(isf, psi)
        velocity_one_form!(isf, psi, isf.hbar)
        # inplace StaggeredSharp
        dinv = inv.(isf.d)
        broadcast!((x, y)-> x .* y, isf.velocity, isf.velocity, (dinv,))
        staggered_advect!(particle, isf)

        #TODO only update active particles
        GLAbstraction.set_arg!(particle_vis, :position, reinterpret(GeometryTypes.Point3f0, Array(particle)))

        yield()
    end
end

@time simloop(
    300, isf2, psi, kvec, omega, n_particles, isjetarr,
    nozzle_rad, nozzle_cen, particles, particle_vis
)


text:
if visualtests
  @testset "Plot recipes" begin
    Random.seed!(2018)
    data₁ = rand(LogNormal(2, 2), 5000)
    data₂ = rand(Pareto(.9, 10), 5000)

    Random.seed!(2018)
    xs = rand(LogNormal(0, 1), 5000)

    @testset "Return levels" begin
      @test_reference "data/ReturnLevels.png" returnplot(xs,label="log-normal")
    end

    @testset "Mean excess" begin
      plt1 = excessplot(data₁, label="Lognormal")
      plt2 = excessplot(data₂, label="Pareto")
      plt = plot(plt1, plt2, size=(600,800), layout=(2,1))
      @test_reference "data/MeanExcess.png" plt
    end

    @testset "Pareto quantile" begin
      plt1 = paretoplot(data₁, label="Lognormal")
      plt2 = paretoplot(data₂, label="Pareto")
      plt = plot(plt1, plt2, size=(600,800), layout=(2,1))
      @test_reference "data/ParetoQuantile.png" plt
    end
  end
end


text:
"""
    SurfaceFluxes

  Surface flux functions, e.g., for buoyancy flux,
  friction velocity, and exchange coefficients.

## Sub-modules
  - module Byun1990
  - module Nishizawa2018

## Interface
  - [`compute_buoyancy_flux`](@ref) computes the buoyancy flux
  - In addition, each sub-module has the following functions:
    - [`monin_obukhov_len`](@ref) computes the Monin-Obukhov length
    - [`compute_friction_velocity`](@ref) computes the friction velocity
    - [`compute_exchange_coefficients`](@ref) computes the exchange coefficients

## References

@article{nishizawa2018surface,
  title={A Surface Flux Scheme Based on the Monin-Obukhov Similarity for Finite Volume Models},
  author={Nishizawa, S and Kitamura, Y},
  journal={Journal of Advances in Modeling Earth Systems},
  volume={10},
  number={12},
  pages={3159--3175},
  year={2018},
  publisher={Wiley Online Library}
}

@article{byun1990analytical,
  title={On the analytical solutions of flux-profile relationships for the atmospheric surface layer},
  author={Byun, Daewon W},
  journal={Journal of Applied Meteorology},
  volume={29},
  number={7},
  pages={652--657},
  year={1990}
}

"""
module SurfaceFluxes

using RootSolvers
using ..MoistThermodynamics
using CLIMAParameters
using CLIMAParameters.Planet: molmass_ratio, grav

"""
    compute_buoyancy_flux(param_set, shf, lhf, T_b, q, α_0)

Computes buoyancy flux given
 - `shf` sensible heat flux
 - `lhf` latent heat flux
 - `T_b` surface boundary temperature
 - `q` phase partition (see [`PhasePartition`](@ref))
 - `α_0` specific volume
"""
function compute_buoyancy_flux(
    param_set::AbstractParameterSet,
    shf,
    lhf,
    T_b,
    q,
    α_0,
)
    FT = typeof(shf)
    _molmass_ratio::FT = molmass_ratio(param_set)
    _grav::FT = grav(param_set)
    cp_ = cp_m(q, param_set)
    lv = latent_heat_vapor(T_b, param_set)
    temp1 = (_molmass_ratio - 1)
    temp2 = (shf + temp1 * cp_ * T_b * lhf / lv)
    return (_grav * α_0 / cp_ / T_b * temp2)
end

module Byun1990

using RootSolvers
using ...MoistThermodynamics
using CLIMAParameters
using CLIMAParameters.SubgridScale: von_karman_const

""" Computes ψ_m for stable case. See Eq. 12 """
ψ_m_stable(ζ, ζ_0, β_m) = -β_m * (ζ - ζ_0)

""" Computes ψ_h for stable case. See Eq. 13 """
ψ_h_stable(ζ, ζ_0, β_h) = -β_h * (ζ - ζ_0)

""" Computes ψ_m for unstable case. See Eq. 14 """
function ψ_m_unstable(ζ, ζ_0, γ_m)
    x(ζ) = sqrt(sqrt(1 - γ_m * ζ))
    temp(ζ, ζ_0) = log1p((ζ - ζ_0) / (1 + ζ_0))  # log((1 + ζ)/(1 + ζ_0))
    ψ_m = (
        2 * temp(x(ζ), x(ζ_0)) + temp(x(ζ)^2, x(ζ_0)^2) - 2 * atan(x(ζ)) +
        2 * atan(x(ζ_0))
    )
    return ψ_m
end

""" Computes ψ_h for unstable case. See Eq. 15 """
function ψ_h_unstable(ζ, ζ_0, γ_h)
    y(ζ) = sqrt(1 - γ_h * ζ)
    ψ_h = 2 * log1p((y(ζ) - y(ζ_0)) / (1 + y(ζ_0))) # log((1 + y(ζ))/(1 + y(ζ_0)))
    return ψ_h
end

"""
    monin_obukhov_len(param_set, u, flux)

Computes the Monin-Obukhov length (Eq. 3)
"""
function monin_obukhov_len(param_set::AbstractParameterSet, u, flux)
    FT = typeof(u)
    _von_karman_const::FT = von_karman_const(param_set)
    return -u^3 / (flux * _von_karman_const)
end

"""
    compute_friction_velocity(param_set, flux, z_0, z_1, β_m, γ_m, tol_abs, iter_max)

Computes roots of friction velocity equation (Eq. 10)

`u_ave = u_* ( ln(z/z_0) - ψ_m(z/L, z_0/L) ) /κ`
"""
function compute_friction_velocity(
    param_set::AbstractParameterSet,
    u_ave::FT,
    flux::FT,
    z_0::FT,
    z_1::FT,
    β_m::FT,
    γ_m::FT,
    tol_abs::FT,
    iter_max::IT,
) where {FT, IT}

    _von_karman_const = FT(von_karman_const(param_set))

    ustar_0 = u_ave * _von_karman_const / log(z_1 / z_0)
    ustar = ustar_0
    let u_ave = u_ave,
        flux = flux,
        z_0 = z_0,
        z_1 = z_1,
        β_m = β_m,
        γ_m = γ_m,
        param_set = param_set

        # use neutral condition as first guess
        stable = z_1 / monin_obukhov_len(param_set, ustar_0, flux) >= 0
        function compute_ψ_m(u)
            L_MO = monin_obukhov_len(param_set, u, flux)
            ζ = z_1 / L_MO
            ζ_0 = z_0 / L_MO
            return stable ? ψ_m_stable(ζ, ζ_0, β_m) : ψ_m_unstable(ζ, ζ_0, γ_m)
        end
        function compute_u_ave_over_ustar(u)
            _von_karman_const = FT(von_karman_const(param_set))
            return (log(z_1 / z_0) - compute_ψ_m(u)) / _von_karman_const # Eq. 10
        end
        compute_ustar(u) = u_ave / compute_u_ave_over_ustar(u)

        if (abs(flux) > 0)
            ustar_1 = compute_ustar(ustar_0)
            sol = RootSolvers.find_zero(
                u -> u_ave - u * compute_u_ave_over_ustar(u),
                ustar_0,
                ustar_1,
                SecantMethod(),
                CompactSolution(),
                tol_abs,
                iter_max,
            )
            ustar = sol.root
        end

    end

    return ustar
end

"""
    compute_exchange_coefficients(param_set, Ri, z_b, z_0, γ_m, γ_h, β_m, β_h, Pr_0)

Computes exchange transfer coefficients:
 - `C_D`  momentum exchange coefficient      (Eq. 36)
 - `C_H`  thermodynamic exchange coefficient (Eq. 37)
 - `L_mo` Monin-Obukhov length               (re-arranged Eq. 3)

TODO: `Pr_0` should come from CLIMAParameters
"""
function compute_exchange_coefficients(
    param_set::AbstractParameterSet,
    Ri::FT,
    z_b::FT,
    z_0::FT,
    γ_m::FT,
    γ_h::FT,
    β_m::FT,
    β_h::FT,
    Pr_0::FT,
) where {FT}
    logz = log(z_b / z_0)
    zfactor = z_b / (z_b - z_0) * logz
    s_b = Ri / Pr_0
    _von_karman_const::FT = von_karman_const(param_set)
    if Ri > 0
        temp = ((1 - 2 * β_h * Ri) - sqrt(1 + 4 * (β_h - β_m) * s_b))
        ζ = zfactor / (2 * β_h * (β_m * Ri - 1)) * temp                    # Eq. 19
        L_mo = z_b / ζ                                             # LHS of Eq. 3
        ζ_0 = z_0 / L_mo
        ψ_m = ψ_m_stable(ζ, ζ_0, β_m)
        ψ_h = ψ_h_stable(ζ, ζ_0, β_h)
    else
        Q_b = FT(1 / 9) * (1 / (γ_m^2) + 3 * γ_h / γ_m * s_b^2)           # Eq. 31
        P_b = FT(1 / 54) * (-2 / (γ_m^3) + 9 / γ_m * (-γ_h / γ_m + 3) * s_b^2) # Eq. 32
        crit = Q_b^3 - P_b^2
        if crit < 0
            T_b = cbrt(sqrt(-crit) + abs(P_b))                     # Eq. 34
            ζ = zfactor * (1 / (3 * γ_m) - (T_b + Q_b / T_b))              # Eq. 29
        else
            θ_b = acos(P_b / sqrt(Q_b^3))                            # Eq. 33
            ζ = zfactor * (-2 * sqrt(Q_b) * cos(θ_b / 3) + 1 / (3 * γ_m))  # Eq. 28
        end
        L_mo = z_b / ζ                                             # LHS of Eq. 3
        ζ_0 = z_0 / L_mo
        ψ_m = ψ_m_unstable(ζ, ζ_0, γ_m)
        ψ_h = ψ_h_unstable(ζ, ζ_0, γ_h)
    end
    cu = _von_karman_const / (logz - ψ_m)            # Eq. 10, solved for u^*
    cth = _von_karman_const / (logz - ψ_h) / Pr_0    # Eq. 11, solved for h^*
    C_D = cu^2                                       # Eq. 36
    C_H = cu * cth                                   # Eq. 37
    return C_D, C_H, L_mo
end

end # Byun1990 module

module Nishizawa2018
using RootSolvers
using ...MoistThermodynamics
using CLIMAParameters
using CLIMAParameters.Planet: grav
using CLIMAParameters.SubgridScale: von_karman_const

""" Computes `R_z0` expression, defined after Eq. 15 """
compute_R_z0(z_0, Δz) = 1 - z_0 / Δz

""" Computes `f_m` in Eq. A7 """
compute_f_m(ζ) = sqrt(sqrt(1 - 15 * ζ))

""" Computes `f_h` in Eq. A8 """
compute_f_h(ζ) = sqrt(1 - 9 * ζ)

""" Computes `ψ_m` in Eq. A3 """
function compute_ψ_m(ζ, L, a)
    FT = typeof(ζ)
    f_m = compute_f_m(ζ)
    return L >= 0 ? -a * ζ :
           log((1 + f_m)^2 * (1 + f_m^2) / 8) - 2 * atan(f_m) + FT(π / 2)
end

""" Computes `ψ_h` in Eq. A4 """
compute_ψ_h(ζ, L, a, Pr) =
    L >= 0 ? -a * ζ / Pr : 2 * log((1 + compute_f_h(ζ)) / 2)

""" Computes `Ψ_m` in Eq. A5 """
function compute_Ψ_m(ζ, L, a, tol)
    FT = typeof(ζ)
    if ζ < tol
        return ζ >= 0 ? -a * ζ / 2 : -FT(15) * ζ / FT(8) # Computes Ψ_m in Eq. A13
    else
        f_m = compute_f_m(ζ)
        # Note that "1-f^3" in is a typo, it is supposed to be "1-f_m^3".
        # This was confirmed by communication with the author.
        return L >= 0 ? -a * ζ / 2 :
               log((1 + f_m)^2 * (1 + f_m^2) / 8) - 2 * atan(f_m) + FT(π / 2) -
               1 + (1 - f_m^3) / (12 * ζ)
    end
end

""" Computes `Ψ_h` in Eq. A6 """
function compute_Ψ_h(ζ, L, a, Pr, tol)
    FT = typeof(ζ)
    if ζ < tol
        return ζ >= 0 ? -a * ζ / (2 * Pr) : -9 * ζ / 4 # Computes Ψ_h in Eq. A14
    else
        f_h = compute_f_h(ζ)
        return L >= 0 ? -a * ζ / (2 * Pr) :
               2 * log((1 + f_h) / 2) + 2 * (1 - f_h) / (9 * ζ)
    end
end

"""
    monin_obukhov_len(param_set, u, θ, flux)

Computes Monin-Obukhov length. Eq. 3
"""
function monin_obukhov_len(param_set::AbstractParameterSet, u, θ, flux)
    FT = typeof(u)
    _von_karman_const::FT = von_karman_const(param_set)
    _grav::FT = grav(param_set)
    return -u^3 * θ / (_von_karman_const * _grav * flux)
end

"""
    compute_friction_velocity(param_set, u_ave, θ, flux, Δz, z_0, a, Ψ_m_tol, tol_abs, iter_max)

Computes friction velocity, in Eq. 12 in, by solving the
non-linear equation:

`u_ave = ustar/κ * ( ln(Δz/z_0) - Ψ_m(Δz/L) + z_0/Δz * Ψ_m(z_0/L) + R_z0 [ψ_m(z_0/L) - 1] )`

where `L` is a non-linear function of `ustar` (see [`monin_obukhov_len`](@ref)).
"""
function compute_friction_velocity(
    param_set::AbstractParameterSet,
    u_ave::FT,
    θ::FT,
    flux::FT,
    Δz::FT,
    z_0::FT,
    a::FT,
    Ψ_m_tol::FT,
    tol_abs::FT,
    iter_max::IT,
) where {FT, IT}
    _von_karman_const::FT = von_karman_const(param_set)
    ustar_0 = u_ave * _von_karman_const / log(Δz / z_0)
    ustar = ustar_0
    let u_ave = u_ave,
        _von_karman_const = _von_karman_const,
        θ = θ,
        flux = flux,
        Δz = Δz,
        z_0 = z_0,
        a = a,
        Ψ_m_tol = Ψ_m_tol,
        tol_abs = tol_abs,
        iter_max = iter_max
        # Note the lowercase psi (ψ) and uppercase psi (Ψ):
        Ψ_m_closure(ζ, L) = compute_Ψ_m(ζ, L, a, Ψ_m_tol)
        ψ_m_closure(ζ, L) = compute_ψ_m(ζ, L, a)
        function compute_u_ave_over_ustar(u)
            L = monin_obukhov_len(param_set, u, θ, flux)
            R_z0 = compute_R_z0(z_0, Δz)
            temp1 = log(Δz / z_0)
            temp2 = -Ψ_m_closure(Δz / L, L)
            temp3 = z_0 / Δz * Ψ_m_closure(z_0 / L, L)
            temp4 = R_z0 * (ψ_m_closure(z_0 / L, L) - 1)
            return (temp1 + temp2 + temp3 + temp4) / _von_karman_const
        end
        compute_ustar(u) = u_ave / compute_u_ave_over_ustar(u)
        ustar_1 = compute_ustar(ustar_0)
        sol = RootSolvers.find_zero(
            u -> u_ave - u * compute_u_ave_over_ustar(u),
            ustar_0,
            ustar_1,
            SecantMethod(),
            CompactSolution(),
            tol_abs,
            iter_max,
        )
        ustar = sol.root
    end
    return ustar
end

"""
    compute_exchange_coefficients(param_set, z, F_m, F_h, a, u_star, θ, flux, Pr)

Computes exchange transfer coefficients:

  - `K_D`  momentum exchange coefficient
  - `K_H`  thermodynamic exchange coefficient
  - `L_mo` Monin-Obukhov length

TODO: `Pr` should come from CLIMAParameters
"""
function compute_exchange_coefficients(
    param_set::AbstractParameterSet,
    z::FT,
    F_m::FT,
    F_h::FT,
    a::FT,
    u_star::FT,
    θ::FT,
    flux::FT,
    Pr::FT,
) where {FT}

    _von_karman_const::FT = von_karman_const(param_set)
    L_mo = monin_obukhov_len(param_set, u_star, θ, flux)
    ψ_m = compute_ψ_m(z / L_mo, L_mo, a)
    ψ_h = compute_ψ_h(z / L_mo, L_mo, a, Pr)

    K_m = -F_m * _von_karman_const * z / (u_star * ψ_m) # Eq. 19
    K_h = -F_h * _von_karman_const * z / (Pr * θ * ψ_h) # Eq. 20

    return K_m, K_h, L_mo
end

end # Nishizawa2018 module

end # SurfaceFluxes module

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libwatcher"
version = v"2.0.5"

# Collection of sources required to build libwatcher
sources = [
    GitSource("https://github.com/JuliaPluto/watcher.git", "18d62d71091a0d9c1fb5c28375bcd7e417d39d66")
]

# Bash recipe for building across all platforms
script = raw"""
cd watcher*/
bash build.sh
install -m 755 build/libwatcher.${dlext} ${libdir}/libwatcher.${dlext}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("aarch64", "macos"),
    Platform("x86_64", "windows"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libwatcher", :libwatcher),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("LibUV_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
#
# This example assumes below packages are available in your environment
#

using StatsModelComparisons, StanSample
using StatsFuns, CSV, Random

ProjDir = @__DIR__

df = CSV.read(joinpath(ProjDir, "..", "..", "data", "cars.csv"), DataFrame)

cars_stan = "
data {
    int N;
    vector[N] speed;
    vector[N] dist;
}
parameters {
    real a;
    real b;
    real sigma;
}
transformed parameters{
    vector[N] mu;
    mu = a + b * speed;
}
model {
    a ~ normal(0, 100);
    b ~ normal(0, 10);
    sigma ~ exponential(1);
    dist ~ normal(mu, sigma)    ;
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(dist[i] | mu[i], sigma);
}
"

Random.seed!(1)
cars_stan_model = SampleModel("cars.model", cars_stan)
data = (N = size(df, 1), speed = df.Speed, dist = df.Dist)
rc = stan_sample(cars_stan_model; data)
println()

if success(rc)
    #stan_summary(cars_stan_model, true)
    nt_cars = read_samples(cars_stan_model);
end

log_lik = nt_cars.log_lik'
n_sam, n_obs = size(log_lik)
lppd = reshape(logsumexp(log_lik .- log(n_sam); dims=1), n_obs)

pwaic = [var(log_lik[:, i]) for i in 1:n_obs]
@show -2(sum(lppd) - sum(pwaic))
println()

@show waic(log_lik)
println()

loo, loos, pk = psisloo(log_lik)
@show -2loo


text:
list_test_str = split(
"..#.#..#####.#.#.#.###.##.....###.##.#..###.####..#####..#....#..#..##..###..######.###...####..#..#####..##..#.#####...##.#.#..#.##..#.#......#.###.######.###.####...#.##.##..#..#..#####.....#.#....###..#.##......#.....#..#..#..##..#...##.######.####.####.#.#...#.......#..#.#.#...####.##.#......#..#...##.#.##..#...##.#.##..###.#......#.#.......#.#.#.####.###.##...#.....####.#..#..#.##.#....##..#.####....##...##..#...#......#.#.......#.......##..####..#...#.#.#...##..#.#..###..#####........#..####......#..#

#..#.
#....
##..#
..#..
..###",
'\n'
)

list_str = let input = @__DIR__() * "/day20_input.txt"
    isfile(input) ? readlines(input) : list_test_str
end

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

function enhance!(img, aux, algorithm, num, k)
    n, m = size(img)
    aux .= ifelse(img[1] == '.', algorithm[1], algorithm[end])
    for j in 2:m-1
        for i in 2:m-1
            r = mapreduce(
                u -> last(u) * 2^(9-first(u)),
                +,
                enumerate(reshape([ c == '#' for c in permutedims(@view(img[i-1:i+1, j-1:j+1]))], 1, :))
            )
            aux[i, j] = algorithm[r+1]
        end
    end
    img .= aux
end

function part1(list, num)
    algorithm = list[1]
    n = length(list) + 2 * num
    m = length(list[end]) + 2 * (num + 1)
    img = fill('.', n, m)
    aux = copy(img)
    img[num+2:end-num-1, num+2:end-num-1] .= permutedims(reduce(hcat, collect.(list[3:end])))
    for k in 1:num
        enhance!(img, aux, algorithm, num, k)
    end
    return img, count(img .== '#')
end

testit("Test Part 1", part1(list_test_str, 2)[2], ==(35))
testit("Challenge Part 1", part1(list_str, 2)[2], ==(5884))
# 5944 too high

# Part 2

function part2(list)
end

testit("Test Part 2", part1(list_test_str, 50)[2], ==(3351))
testit("Challenge Part 2", part1(list_str, 50)[2], ==(19043))


text:
import DandelionWebSockets: on_text, on_binary,
    state_connecting, state_open, state_closing, state_closed,
    write, handle, FrameFromServer

import Base.==

#
# A lot of tests use WebSocket frames, naturally, so these are common frames that all tests can use.
#

nomask = Array{UInt8,1}()
testmask = b"\x37\xfa\x21\x3d"

# A single frame with payload "Hello"
test_frame1 = Frame(true,  OPCODE_TEXT,         false, 5, 0, nomask, b"Hello")

# Two fragments of a text message with payload "Hello"
test_frame2 = Frame(false, OPCODE_TEXT,         false, 3, 0, nomask, b"Hel")
test_frame3 = Frame(true,  OPCODE_CONTINUATION, false, 2, 0, nomask, b"lo")

# A single text frame, masked, with body "Hello"
test_frame4 = Frame(true,  OPCODE_TEXT, true, 5, 0, testmask, b"\x7f\x9f\x4d\x51\x58")
test_bin_frame4 = Frame(true,  OPCODE_BINARY, true, 5, 0, testmask, b"\x7f\x9f\x4d\x51\x58")

mask2 = b"\x17\x42\x03\x7f"

# Two masked fragments, one initial and one final. They are masked by two different masks.
test_frame5 = Frame(false, OPCODE_TEXT, true, 3, 0, testmask, b"\x7f\x9f\x4d")
test_frame6 = Frame(true, OPCODE_CONTINUATION, true, 2, 0,  mask2, b"\x7b\x2d")

# Two binary fragments, one initial and one final.
frame_bin_start = Frame(false, OPCODE_BINARY,       false, 3, 0, nomask, b"Hel")
frame_bin_final = Frame(true,  OPCODE_CONTINUATION, false, 2, 0, nomask, b"lo")
frame_bin_1     = Frame(true,  OPCODE_BINARY,       false, 5, 0, nomask, b"Hello")

server_close_frame = Frame(true, OPCODE_CLOSE, false, 0, 0, nomask, b"")
client_close_reply = Frame(true, OPCODE_CLOSE, true, 0, 0, testmask, b"")
server_ping_frame = Frame(true, OPCODE_PING, false, 0, 0, nomask, b"")
client_pong_frame = Frame(true, OPCODE_PONG, true, 0, 0, testmask, b"")
server_pong_frame = Frame(true, OPCODE_PONG, false, 0, 0, nomask, b"")
client_ping_frame = Frame(true, OPCODE_PING, true, 0, 0, testmask, b"")
server_ping_frame_w_pay = Frame(true, OPCODE_PING, false, 5, 0, nomask, b"Hello")
client_pong_frame_w_pay = Frame(true, OPCODE_PONG, true, 5, 0, testmask, b"\x7f\x9f\x4d\x51\x58")

zero256 = Array{UInt8}([UInt8(0) for x in range(1, length=256)])
zero65k = Array{UInt8}([UInt8(0) for x in range(1, length=65536 + 1024)])
zero256_masked = Array{UInt8}([zero256[i] ⊻ testmask[(i-1)%4 + 1] for i in 1:length(zero256)])
zero65k_masked = Array{UInt8}([zero65k[i] ⊻ testmask[(i-1)%4 + 1] for i in 1:length(zero65k)])

# Binary message, payload is 256 bytes, single masked
test_bin_frame_256 = Frame(true, OPCODE_BINARY, true, 126, 256, testmask, zero256_masked)

# Binary message, payload is 64KiB, masked
test_bin_frame_65k = Frame(true, OPCODE_BINARY, true, 127, 65536 + 1024, testmask, zero65k_masked)

#
# To accurately test a fake TCPSocket I need a blocking stream.
# The implementation below is meant to be simple, not performant or good.
#
struct BlockingStream <: IO
    buf::IOBuffer
end

function blocking_read(s::BlockingStream)
    x = nothing
    while true
        try
            x = read(s.buf, UInt8)
            return x
        catch ex
            if !isa(ex, EOFError)
                rethrow(ex)
            end
        end
        sleep(0.05)
    end
    x
end

function Base.read(s::BlockingStream, ::Type{UInt8})
    blocking_read(s)
end

function Base.read(s::BlockingStream, ::Type{Array{UInt8, 1}}, n::Int)
    buf = Array{UInt8, 1}(n)
    for i in 1:n
        buf[i] = blocking_read(s)
    end
    buf
end

function Base.read(s::BlockingStream, ::Type{UInt64})
    buf = read(s, Array{UInt8, 1}, 8)
    x::UInt64 =
        buf[1] << 56 | buf[2] << 48 << buf[3] << 40 | buf[4] << 32 |
        buf[5] << 24 | buf[6] << 16 << buf[7] << 8  | buf[8]
    x
end

function Base.read(s::BlockingStream, ::Type{UInt16})
    buf = read(s, Array{UInt8, 1}, 2)
    x::UInt16 = buf[1] << 8 | buf[2]
    x
end

#
# MockClientProtocol mocks ClientProtocol, and should have used the @mock macro, except that there are
# issues with FactCheck and doing asserts in other tasks. This custom mock ensures that all asserts
# are done afterwards, in the same task that created it.
#

==(a::FrameFromServer, b::FrameFromServer) = a.frame == b.frame

const MockLogicCall = Tuple{Symbol, Vector{Any}}

mutable struct MockClientProtocol <: AbstractClientProtocol
    actuals::Vector{MockLogicCall}
    expected::Vector{MockLogicCall}

    MockClientProtocol(expected::Vector{MockLogicCall}) = new([], expected)
end

function call(m::MockClientProtocol, s::Symbol, args...)
    push!(m.actuals, (s, collect(args)))
end

function check(m::MockClientProtocol)
    #@fact m.actuals --> m.expected
end

handle(m::MockClientProtocol, args...) = call(m, :handle, args...)

#
# A fake stream for checking that we read and write the right frames.
#

struct FakeFrameStream <: IO
    reading::Vector{Frame}
    writing::Vector{Frame}
    close_on_empty::Bool
    stop_chan::Channel{Symbol}

    FakeFrameStream(reading::Vector{Frame}, writing::Vector{Frame}, close_on_empty::Bool) =
        new(reading, writing, close_on_empty, Channel{Symbol}(32))
end

function Base.read(s::FakeFrameStream, ::Type{Frame})
    if isempty(s.reading)
        if s.close_on_empty
            throw(EOFError())
        else
            take!(s.stop_chan)
            throw(EOFError())
        end
    end
    sleep(0.2)
    shift!(s.reading)
end

function Base.write(s::FakeFrameStream, frame::Frame)
    push!(s.writing, frame)
    if frame.opcode == DandelionWebSockets.OPCODE_CLOSE
        put!(s.stop_chan, :stop)
    end
end

Base.close(s::FakeFrameStream) = nothing


text:
# Utilities for benchmarking NN Descent

recall(nn, true_nn) = sum(_recall(nn[:,i], true_nn[:,i]) for i in 1:size(nn,2))/size(nn,2)
_recall(π, πₜ) = length(intersect(π, πₜ))/length(πₜ)


text:
using GalerkinSparseGrids
using Test
using HCubature
using StaticArrays
using SparseArrays

#--------------------------------------
# Testing differentiation
#--------------------------------------

@testset "differentiation.jl" begin
    @info "Testing differentiation 1-D DG basis... "
    @testset "1D DG Differentiation" begin
        k=3
        for l in 2:5
            frefVD = V2Dref(Val(1), k, l, Val(:sparse))
            frefDV = D2Vref(Val(1), k, l, Val(:sparse))
            D_op = D_matrix(1, k, l, frefVD, frefDV; scheme="full")
            vcoeffs = vcoeffs_DG(1, k, l, x->cos(2*pi*x[1]); scheme="full")
            dvcoeffs = D_op * vcoeffs
            dict= V2D(1, k, l, vcoeffs; scheme="full")
            ddict = V2D(1, k, l, dvcoeffs; scheme="full")
            err = x->(reconstruct_DG(ddict,[x[1]])+2*pi*sin(2*pi*x[1]))^2
            @test hquadrature(err, 0, 1, atol=1.0e-10, maxevals=500)[1] < 1/(1<<(k+l-2))
        end
    end

    @info "Testing differentiation 2-D full DG basis... "
    @testset "2D Full DG Differentiation" begin
        D = 2
        k = 3
        for l in 2:5
            D_op = D_matrix(D, 1, k, l; scheme="full")
            vcoeffs = vcoeffs_DG(D, k, l, x->cos(2*pi*x[1])*cos(2*pi*x[2]); scheme="full")
            dvcoeffs = D_op * vcoeffs
            dict= V2D(D, k, l, vcoeffs; scheme="full")
            ddict = V2D(D, k, l, dvcoeffs; scheme="full")
            err = x->(reconstruct_DG(ddict,[x[1],x[2]])+2*pi*sin(2*pi*x[1])*cos(2*pi*x[2]))^2
            @test hcubature(err, [0,0], [1,1], atol=1.0e-10, maxevals=500)[1] < 1/(1<<(k+l-2))
        end
    end
end


text:
"""
Analytic ITER-like Solov'ev equilibrium with X-point.
"""
module GuidingCenter4dSolovevIterXpoint

    using ElectromagneticFields.Solovev

    export initial_conditions_barely_passing, initial_conditions_barely_trapped,
           initial_conditions_deeply_passing, initial_conditions_deeply_trapped,
           initial_conditions_trapped

    # export toroidal_momentum

    Solovev.@code_iter_xpoint() # inject magnetic field code

    include("coordinate_transformations.jl")
    include("gc_common.jl")
    include("gc_equations.jl")

    function solovev_xpoint_iter_initial_conditions(t₀, q₀, μ)
        local ω₀ = ωabs(t₀, q₀)
        local params = (μ=μ, R₀=R₀, ω₀=ω₀)
        (transform_q_to_q̃(t₀, q₀, params), params)
    end

    const x₀ = from_cartesian(0, [2.5, 0., 0.])
    const μ₀ = 1E-2

    initial_conditions_barely_passing() = solovev_xpoint_iter_initial_conditions(0, [x₀..., 3.425E-1], μ₀)
    initial_conditions_barely_trapped() = solovev_xpoint_iter_initial_conditions(0, [x₀..., 3.375E-1], μ₀)
    initial_conditions_deeply_passing() = solovev_xpoint_iter_initial_conditions(0, [x₀..., 5.0E-1  ], μ₀)
    initial_conditions_deeply_trapped() = solovev_xpoint_iter_initial_conditions(0, [x₀..., 1.0E-1  ], μ₀)

    initial_conditions_trapped() = solovev_xpoint_iter_initial_conditions(0, [from_cartesian(0, [7.0, 0., 0.])..., -2E-3], 1.88E-7)

    include("../guiding_center_4d/guiding_center_4d_diagnostics.jl")

end


text:
mutable struct EdaxPlayer <: Player
    proc::Base.Process
    # TODO options

    function EdaxPlayer()
        c = Cmd(`./mEdax -n 1 -ponder off`, dir="/Users/sglyon/Downloads/edax/4.3/bin")
        proc = open(c, "w+")
        if !process_running(proc)
            error("error starting edax process")
        end
        read(proc, 1)
        read(proc.out, nb_available(proc.out))
        out = new(proc)
        finalizer(out) do x
            close(x.proc)
        end
        out
    end
end

function edax_board(s::State)
    b = s.board
    junk = fill('-', 8, 8)
    for (r, c) in bits_to_tuples(b.p1_placed)
        junk[c, r] = 'X'
    end
    for (r, c) in bits_to_tuples(b.p2_placed)
        junk[c, r] = 'O'
    end
    player = s.player == 1 ? " X" : " O"
    string(join(junk, ""), player)
end

function Base.write(e::EdaxPlayer, s::String)
    write(e.proc.in, s)
    out = read(e.proc.out, 1)
    append!(out, read(e.proc.out, nb_available(e.proc.out)))
    String(out)
end

set_board!(e::EdaxPlayer, s::State) = write(e, "setboard $(edax_board(s))\n")

function _get_edax_board(e::EdaxPlayer)
    write(e.proc.in, "\n")
    msg = readuntil(e.proc.out, "\n>")
    parse_edax_board(msg)
end

function parse_edax_board(msg::String)
    lines = split(msg, "\n")

    board_start = findfirst(line -> startswith(line, "  A"), lines)
    board = lines[board_start:board_start+9]

    markers_in_row = 3:2:17

    # extract tile positions
    p1_placed = UInt64(0)
    p2_placed = UInt64(0)
    for row in 1:8
        data = board[row+1]
        markers = data[markers_in_row]
        for col in 1:8
            if markers[col] == '*'
                p1_placed += encode(row, col)
            end
            if markers[col] == 'O'
                p2_placed += encode(row, col)
            end
        end
    end

    # extract turn
    turn = board[6][24] == 'B' ? 1 : 2

    # assume round is
    round = score(p1_placed) + score(p2_placed) + 1

    State(Board(p1_placed, p2_placed), turn, round)
end

function select_action(e::EdaxPlayer, s::State)
    if s.round < 4
        return rand(legal_actions(s))
    end
    set_board!(e, s)
    msg = write(e, "go\n")
    my_match = match(r"Edax plays (\w)(\d)", msg)
    row = parse(Int, my_match[2])
    col = Int(my_match[1][1]) - 64
    (row, col)
end


text:
function get_module_2_excersise_3_author_name()
    return "YOUR NAME HERE"
end

function wave_function(x, quantum_number, box_length)
    return 0 #implement the wave function here use the built in sin(x) function
end

# Numerical integration procedure for one dimensional particle in a box
# PIB = Particle In a Box
function one_dimensional_PIB_numerical_integration(quantum_number,
    box_length,
    number_of_rectangles,
    integration_start,
    integration_end)

    area_under_curve = 0
    # implement the function here with a loop
    return area_under_curve
end


text:
module DiffPrivacy
using Random, Distributions

function get_one_parallel_database(A, remove_index)

    # Returns an array A, without the element at position remove_index, of size (1, length(A) - 1)

    return vcat(A[1:remove_index-1], A[remove_index+1:length(A)])'
end

function get_multiple_parallel_databases(A)

    # Returns Array of length(A) with all the parallel databases generated and stored in parallel_databases

    parallel_databases = []

    for i in 1:length(A)
        x = get_one_parallel_database(A, i)
        push!(parallel_databases, x)
    end

    return parallel_databases
end

function sum_sensitivity(A)

    db = A
    pdbs = get_multiple_parallel_databases(A)

    full_db_result = sum(db)

    max_distance = 0
    for pdb in pdbs
        pdb_result = sum(pdb)

        db_distance = abs(pdb_result - full_db_result)

        if(db_distance > max_distance)
            max_distance = db_distance
        end
    end

    return max_distance
end

function mean_sensitivity(A)

    db = A
    pdbs = get_multiple_parallel_databases(A)

    full_db_result = sum(db) / length(db)

    max_distance = 0
    for pdb in pdbs
        pdb_result = sum(pdb) / length(pdb)

        db_distance = abs(pdb_result - full_db_result)

        if(db_distance > max_distance)
            max_distance = db_distance
        end
    end

    return max_distance
end

function sum_sensitive_query(A, threshold)
    if sum(A) > threshold
        return 1
    else
        return 0
    end
end

function sum_sensitivity_threshold(A, threshold)

    db = A
    pdbs = get_multiple_parallel_databases(A)

    full_db_result = sum_sensitive_query(db, threshold)

    max_distance = 0
    for pdb in pdbs
        pdb_result = sum_sensitive_query(pdb, threshold)

        db_distance = abs(pdb_result - full_db_result)

        if(db_distance > max_distance)
            max_distance = db_distance
        end
    end

    return max_distance
end

function add_gaussian_noise(A, mu, sigma)
    d = rand(Gaussian(mu, sigma), size(A))
    A = A + d
    return A
end

function add_laplacian_noise(A, loc, scale)
    d = rand(Laplacian(loc, scale), size(A))
    A = A + d
    push!(A, d)
    return A
end

function add_normal_noise(A, mu, sigma)
    d = rand(Normal(mu, sigma), size(A))
    A = A + d
    return A
end

export DiffPrivacy

end # module


text:
using PolynomialFactors
using Test

include("test-utils.jl")
include("test-factor.jl")


text:
using Base.Test
using Compat
using Images, FixedPointNumbers
import VideoIO

testdir = joinpath(Pkg.dir("VideoIO"), "test")
videodir = joinpath(Pkg.dir("VideoIO"), "videos")

VideoIO.TestVideos.available()
VideoIO.TestVideos.download_all()

swapext(f, new_ext) = "$(splitext(f)[1])$new_ext"

println(STDERR, "Testing file reading...")

function notblank(img)
    all(green(img) .== 0x00uf8) || all(blue(img) .== 0x00uf8) || all(red(img) .== 0x00uf8) || maximum(reinterpret(Ufixed8, img)) < 0xcfuf8
end

for name in VideoIO.TestVideos.names()
    @osx_only startswith(name, "crescent") && continue
    println(STDERR, "   Testing $name...")

    first_frame_file = joinpath(testdir, swapext(name, ".png"))
    first_frame = imread(first_frame_file) # comment line when creating png files

    f = VideoIO.testvideo(name)
    v = VideoIO.openvideo(f)

    img = read(v, Image)

    # Find the first non-trivial image
    while notblank(img)
        read!(v, img)
    end

    #imwrite(img, first_frame_file)        # uncomment line when creating png files

    @test img == first_frame               # comment line when creating png files

    while !eof(v)
        read!(v, img)
    end

    # read first frames again, and compare
    seekstart(v)

    read!(v, img)

    while notblank(img)
        read!(v, img)
    end

    @test img == first_frame

    close(v)
end

println(STDERR, "Testing IO reading...")
for name in VideoIO.TestVideos.names()
    @osx_only startswith(name, "crescent") && continue
    # TODO: fix me?
    (startswith(name, "ladybird") || startswith(name, "NPS")) && continue

    println(STDERR, "   Testing $name...")
    first_frame_file = joinpath(testdir, swapext(name, ".png"))
    first_frame = imread(first_frame_file) # comment line when creating png files

    filename = joinpath(videodir, name)
    v = VideoIO.openvideo(open(filename))

    img = read(v, Image)

    # Find the first non-trivial image
    while notblank(img)
        read!(v, img)
    end

    #imwrite(img, first_frame_file)        # uncomment line when creating png files

    @test img == first_frame               # comment line when creating png files

    while !eof(v)
        read!(v, img)
    end
end

VideoIO.testvideo("ladybird") # coverage testing
@test_throws ErrorException VideoIO.testvideo("rickroll")
@test_throws ErrorException VideoIO.testvideo("")

#VideoIO.TestVideos.remove_all()


text:
export cross_section, radial_sample, punch, bulk, multi_region_sampling

function cross_section(tumor; x = nothing, y = nothing, z = nothing, width = 2.0, reduce_dim = false)
    for (d, coord) in enumerate((x, y, z))
        if !isnothing(coord)
            rows = -width .< getindex.(tumor.position, d) .- coord .< width
            if reduce_dim
                plane = tumor[rows, :]
                plane.position = plane.position .|> p -> reduce_dimension(p, d)
                return plane
            else
                return tumor[rows, :]
            end
        end
    end
end

# this function returns a copy without entries at dimension d
function reduce_dimension(position, d)
    dim = position |> length
    select = mod.((d-1:2:d+1) .- 1, dim) .+ 1 |> unique
    return getindex(position, select)
end

radial_sample(tumor; r = 2.0, width = 2.0) = tumor[map(p -> abs(norm(p) - r) < width, tumor.position), :] |> deepcopy

punch(tumor; pos, r = 10.0) = tumor[map(p -> norm(p .- pos) < r, tumor.position), :] |> deepcopy

bulk(tumor; pos, box) = tumor[map(p -> all(abs.(p .- pos) .< box ./ 2), tumor.position), :] |> deepcopy

struct Square
    anchorpoint
    width
end

Base.in(point, c::Circle) = norm(point .- c.center) < c.r

function triangular_lattice(sq::Square; n = 0, a = sqrt(π / n) * radius)
    v, w = [1.0, 0.0], [0.5, -cosd(30)]

    steps_x = ceil(Int, sq.width / a)
    steps_y = ceil(Int, sq.width / (a * cosd(30)))

    lattice = Matrix{Vector{Float64}}(undef, steps_x, steps_y)
    lattice[:, 1] .= [sq.anchorpoint .+ a * i * v for i = 1:steps_x]
    for j = 2:steps_y
        @. lattice[:, j] = [p + a * (w - iseven(j) * v) for p in lattice[:, j-1]]
    end
    return lattice
end

function multi_region_sampling(tumor; n = 0, a = 0.0, cells_per_sample = 0, sample_r = a / 2)

    cm = mean(tumor.position)
    r = mean(norm.(p .- cm for p in tumor.position))
    density = size(punch(tumor; pos = cm, r = r), 1) / (π * r^2)
    R = sqrt(size(tumor, 1) / (π * density))

    sample_r = iszero(cells_per_sample) ? a / 2 : sqrt(cells_per_sample / (π * density))
    if !iszero(n)
        a = sqrt(π / n) * R
        sample_r = iszero(sample_r) ? a / 2 : sample_r
    elseif iszero(a)
        a = 2 * sample_r
    end
    iszero(sample_r) && error("specify parameters")
    sample_r > a / 2 && error("samples too large for given spacing")

    square = Square(cm .- R .* (1, -1), 2 * R)
    lattice = [triangular_lattice(square; n = n, a = a)...]
    filter!(p -> p in Circle(Point{2}(cm), R - a / 2), lattice)

    samples = lattice .|> p -> punch(tumor; pos = p, r = sample_r)

    return (lattice = lattice, samples = samples, sample_r = sample_r)
end

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
module Hop
using StaticArrays, LinearAlgebra

export TightBindingModel, sethopping!,
    gethamiltonian, geteig, geteigvals,
    KPath, getband

struct TightBindingModel
    "number of orbits"
    norbits::Int
    "lattice vector"
    lat::Matrix{Float64}
    "reciprocal lattice vector"
    rlat::Matrix{Float64}
    "atom positions"
    positions::Matrix{Float64}
    "hoppings"
    hoppings::Dict{SVector{3,Int64},Matrix{ComplexF64}}
end

"""
```julia
TightBindingModel(lat::Matrix{Float64}, positions::Matrix{Float64}; spinful::Bool=false)
```
A tight binding model.

Generally, every tight binding model is three dimensional.
Therefore, every lattice vector should have three components.
Lower dimensional models should be simulated by vacuum layer.

# Constructor Arguments
- `lat`: lattice vector.
  Lattice vectors should be provided in columns.
- `positions`: atom positions in reduced coordinate.
  Atom positions should be provided in columns.
- `spinful`: false for spinless systems and true for spinful systems. If `spinful`
  is true, `norbits` will be twice the number of `size(positions, 1)`. Orbits are
  ordered as (|1↑⟩, |1↓⟩, |2↑⟩, |2↓⟩, ...).

# Fields
- `norbits::Int`: number of orbits.
- `lat::Matrix{Float64}`: lattice vectors stored in columns.
- `rlat::Matrix{Float64}`: reciprocal lattice vectors stored in columns.
- `positions::Matrix{Float64}`: position of orbits in reduced coordinate stored in columns.
- `hoppings::Dict{SVector{3,Int64},Matrix{ComplexF64}}`: hoppings stored as
   R->⟨0m|H|Rn⟩.
"""
function TightBindingModel(lat::Matrix{Float64}, positions::Matrix{Float64}; spinful::Bool=false)
    @assert size(lat) == (3, 3) "Size of lat is not correct."
    @assert size(positions, 1) == 3 "Size of positions is not correct."
    rlat = 2*π*inv(lat)'
    if !spinful
        return TightBindingModel(size(positions, 2), lat, rlat, positions, Dict())
    else
        spinpositions = zeros(3, size(positions, 2)*2)
        for i in 1:size(positions, 2)
            spinpositions[:, 2*i-1] = positions[:, i]
            spinpositions[:, 2*i] = positions[:, i]
        end
        return TightBindingModel(size(positions, 2)*2, lat, rlat, spinpositions, Dict())
    end
end

function Base.show(io::IO, t::TightBindingModel)
    print(io, "TightBindingModel: $(t.norbits) orbitals.")
end

"""
```julia
sethopping!(t::TightBindingModel, m::Int64, n::Int64, R, hopping;
    mode::Symbol=:a)
```

Set ⟨0m|H|Rn⟩ to `hopping`. `hopping::Number` for spinless models and
`hopping::Matrix{<:Number}` for spinful models. For spinful models,
`size(hopping)` should be (2, 2) and the basis for `hopping` is (|↑⟩, |↓⟩).
`mode` has two possible values: :a for add mode and :s for set or reset mode.
"""
function sethopping!(t::TightBindingModel, m::Int64, n::Int64, R,
    hopping::Number; mode::Symbol=:a)
    @assert (m in 1:t.norbits) && (n in 1:t.norbits) "No such orbit."
    @assert length(R) == 3 "Size of R is not correct."
    if !(R in keys(t.hoppings))
        @assert !(-R in keys(t.hoppings))
        t.hoppings[R] = zeros(ComplexF64, (t.norbits, t.norbits))
        t.hoppings[-R] = zeros(ComplexF64, (t.norbits, t.norbits))
    end
    if mode == :s
        t.hoppings[R][m, n] = hopping
        t.hoppings[-R][n, m] = conj(hopping)
    elseif mode == :a
        t.hoppings[R][m, n] += hopping
        if (m, n, R) != (n, m, -R) # not onsite energy
            t.hoppings[-R][n, m] += conj(hopping)
        end
    else
        error("Unknown mode.")
    end
    return nothing
end

# For spinful model
function sethopping!(t::TightBindingModel, m::Int64, n::Int64, R,
    hopping::Matrix{<:Number}; mode::Symbol=:a)
    @assert iseven(t.norbits) "Not a spinful model."
    @assert (m in 1:(t.norbits÷2)) && (n in 1:(t.norbits÷2)) "No such orbit."
    @assert size(hopping) == (2, 2) "Size of hopping is not correct."
    sethopping!(t, 2m-1, 2n-1, R, hopping[1, 1], mode=mode)
    sethopping!(t, 2m, 2n-1, R, hopping[2, 1], mode=mode)
    if (m, n, R) != (n, m, -R) # not onsite energy
        sethopping!(t, 2m-1, 2n, R, hopping[1, 2], mode=mode)
    end
    sethopping!(t, 2m, 2n, R, hopping[2, 2], mode=mode)
    return nothing
end

"""
```julia
gethamiltonian(t::TightBindingModel, k::Vector{<:Real}) --> Matrix{ComplexF64}
```

Calculate Hamiltonian of a TightBindingModel t for a specific k point. k should
be provided in reduced coordinate.
"""
function gethamiltonian(t::TightBindingModel, k::Vector{<:Real})
    @assert size(k) == (3,) "Size of k is not correct."
    h = zeros(ComplexF64, (t.norbits, t.norbits))
    for (R, hopping) = t.hoppings
        h += exp(2π*im*(k⋅R))*hopping
    end
    @assert ishermitian(h)
    return h
end

"""
```julia
geteig(t::TightBindingModel, k::Vector{<:Real})
```

Calculate eigenvalues and eigenvectors of t. k should be provided in reduced coordinate.

# Return
`(egvals::Vector{Float64}, egvecs::Matrix{ComplexF64})`,
Eigenvectors are stored in columns and eigenvalues are sorted from small to large.
"""
function geteig(t::TightBindingModel, k::Vector{<:Real})
    @assert size(k) == (3,) "Size of k is not correct."
    hamiltonian = gethamiltonian(t, k)
    egvals, egvecs = eigen(Hermitian(hamiltonian))
    perm = sortperm(egvals)
    return (egvals[perm], egvecs[:, perm])
end

"""
```julia
geteigvals(t::TightBindingModel, k::Vector{<:Real})
  --> Vector{Float64}
```

Calculate eigenvalues and eigenvectors of t. k should be provided in reduced coordinate.

Eigenvalues are sorted from small to large.
"""
function geteigvals(t::TightBindingModel, k::Vector{<:Real})
    @assert size(k) == (3,) "Size of k is not correct."
    hamiltonian = gethamiltonian(t, k)
    return sort(eigvals(Hermitian(hamiltonian)))
end

struct KPath
    rlat::Matrix{Float64}
    nodes::Matrix{Float64}
    distances::Vector{Float64}
    kpts::Matrix{Float64}
    ndiv::Int64
    nkpts::Int64
end

"""
```julia
KPath(rlat::Matrix{<:Real}, nodes::Matrix{<:Real}, ndiv::Int64)
```

A path in reciprocal space.

# Constructor Arguments
 - `rlat`: reciprocal lattice vector stored in columns.
 - `nodes`: start and end points of the path stored in columns in reduced coordinates.
 - `ndiv`: number of divisions on each section of the path.

# Fields
 - `rlat::Matrix{Float64}`: see explanations above.
 - `nodes::Matrix{Float64}`: see explanations above.
 - `distances::Vector{Float64}`: distance along the path for each k point.
 - `kpts::Matrix{Float64}`: k points stored in columns  in reduced coordinates.
 - `ndiv::Int64`: see explanations above.
 - `nkpts::Int64`: number of k points on the path.
"""
function KPath(rlat::Matrix{<:Real}, nodes::Matrix{<:Real}, ndiv::Int64)
    @assert iseven(size(nodes, 2))
    npaths = size(nodes, 2)÷2
    nkpts = ndiv*npaths
    distances = zeros(nkpts)
    kpts = zeros(3, nkpts)
    for ipath = 1:npaths
        dk = (nodes[:, 2*ipath]-nodes[:, 2*ipath-1])/(ndiv-1)
        dkn = norm(rlat*dk) # real norm of dk
        if ipath == 1
            distance0 = 0
        else
            distance0 = distances[(ipath-1)*ndiv]
        end
        for ikpt = 1:ndiv
            distances[(ipath-1)*ndiv+ikpt] = dkn*(ikpt-1) + distance0
            kpts[:, (ipath-1)*ndiv+ikpt] = nodes[:, 2*ipath-1]+dk*(ikpt-1)
        end
    end
    return KPath(rlat, nodes, distances, kpts, ndiv, nkpts)
end

function Base.show(io::IO, kp::KPath)
    print(io, "KPath: $(kp.nkpts÷kp.ndiv) line(s).")
end

"""
```julia
getband(t::TightBindingModel, kp::KPath)
    --> Matrix{Float64}
```

Calculate bands of `t` on `kp`.

Return energies at each k point stored in columns.
"""
function getband(t::TightBindingModel, kp::KPath)
    bands = zeros(t.norbits, kp.nkpts)
    for ik=1:kp.nkpts
        bands[:, ik] = geteigvals(t, kp.kpts[:, ik])
    end
    return bands
end

include("floquet.jl")

end


text:
@testset "vector projections" begin
    n̂  = SVector{3,Float64}(0, 0, 1)
    v1 = SVector{3,Float64}(1, 0, 0)
    v2 = SVector{3,Float64}(0, 1, 1)

    @test vec_proj(v1, n̂) == 0.0 * n̂
    @test vec_proj(n̂, n̂) == 1.0 * n̂
    @test vec_proj(v2, n̂) == 1.0 * n̂

    @test vec_sub_vec_proj(v1, n̂) == v1
    @test vec_sub_vec_proj(n̂, n̂) == 0.0 * n̂
    @test vec_sub_vec_proj(v2, n̂) == SVector{3,Float64}(0, 1, 0)
end

@testset "a_dot_one_pad_b" begin
	a = SVector{4,Float64}(1.0, 2.0, 3.0, 4.0)
	b1 = SVector{3,Float64}(2.0, 0.0, 0.0)
	b2 = SVector{3,Float64}(1.0, 2.0, 0.0)
	b3 = SVector{3,Float64}(1.0, 2.0, 3.0)

	@test a_dot_one_pad_b(a, b1) == 6.0
	@test a_dot_one_pad_b(a, b2) == 9.0
	@test a_dot_one_pad_b(a, b3) == 18.0
end


text:
@testset "centering" begin

    p1 = Point([4.831961989402771, 45.75764678012361])
    l1 = LineString([[
        4.86020565032959,
        45.76884015325622
      ],
      [
        4.85994815826416,
        45.749558161214516
      ]])

    poly = Polygon([[[
          4.8250579833984375,
          45.79398056386735
        ],
        [
          4.882392883300781,
          45.79254427435898
        ],
        [
          4.910373687744141,
          45.76081677972451
        ],
        [
          4.894924163818359,
          45.7271539426975
        ],
        [
          4.824199676513671,
          45.71337148333104
        ],
        [
          4.773387908935547,
          45.74021417890731
        ],
        [
          4.778022766113281,
          45.778418789239055
        ],
        [
          4.8250579833984375,
          45.79398056386735
        ]]])

    @test centroid(p1).coordinates == [4.831961989402771, 45.75764678012361]
    @test centroid(l1).coordinates == [4.860076904296875,45.75919915723537]
    @test centroid(poly).coordinates == [4.839177131652832,45.76256007199914]

    poly2 = Polygon([[[0, 0], [0, 0], [0, 0], [0, 0]]])

    @test meancenter(poly).coordinates ≈ [4.839177,45.76256]
    @test meancenter(poly2).coordinates == [0, 0]
    @test meancenter(l1).coordinates ≈ [4.860077,45.759199]
    @test meancenter(p1).coordinates == p1.coordinates

    @test masscenter(poly2).coordinates == [0, 0]
    @test masscenter(poly).coordinates == [4.840728965137111,45.75581209996416]

    coll = FeatureCollection([Feature(Point([0, 0])),
      Feature(Point([1, 0])), Feature(Point([0, 1])),
      Feature(Point([5, 8]))])

    coll2 = FeatureCollection([Feature(Point([0, 0])),Feature(Point([9, 9])), Feature(Point([9.25, 9.25])),
      Feature(Point([9.5, 9.5])), Feature(Point([9.75, 9.75])), Feature(Point([10, 10]))])

    @test mediancenter(coll).coordinates ≈ [0.383876,0.616989] atol=0.0001
    @test mediancenter(coll2).coordinates ≈ [9.254246,9.254246] atol=0.0001

    fc = FeatureCollection([Feature(Point([4.833351373672485,45.760809294695534])),
      Feature(Point([4.8331475257873535,45.760296567821456])), Feature(Point([4.833984374999999,45.76073818687033])),
      Feature(Point([4.834005832672119,45.76022171678877]))])

    @test centroid(fc).coordinates ≈ [4.8336222767829895,45.76051644154402] atol=0.001
    @test center(fc).coordinates ≈ centroid(fc).coordinates atol=0.01
end


text:
using PSFModels

# TODO - figure out why using `ones(1, 1)` fails only in test
@testset "image injection" begin
    kernel = [0 0 0; 0 1 0; 0 0 0]
    frame = zeros(3, 3)
    @test inject(frame, kernel; x=1, y=3)[1, 3] ≈ 1
    @test inject(frame, kernel; x=2, y=2, amp=2)[2, 2] ≈ 2
    @test inject(frame, kernel; x=10, y=10) ≈ zeros(3, 3) atol = sqrt(eps(eltype(frame)))

    cube = zeros(3, 3, 10)
    angles = 90 .* ones(10)
    @test inject(cube, kernel; x=2, y=2)[2, 2, :] ≈ ones(10)
    @test inject(cube, kernel, angles; x=2, y=2)[2, 2, :] ≈ ones(10)

    @testset "synthetic psf - $mod" for mod in [gaussian, moffat, airydisk]
        frame = zeros(3, 3)
        @test inject(frame, mod; x=1, y=3, fwhm=3)[1, 3] ≈ 1
        @test inject(frame, mod; x=2, y=2, fwhm=3, amp=2)[2, 2] ≈ 2

        cube = zeros(3, 3, 10)
        @test inject(cube, mod; x=2, y=2, fwhm=3)[2, 2, :] ≈ ones(10)
        @test inject(cube, mod, angles; x=2, y=2, fwhm=3)[2, 2, :] ≈ ones(10)
    end
end

@testset "inject - AnnulusView" begin
    cube = AnnulusView(zeros(21, 21, 10))
    angles = 90 .* ones(10)

    kernel = [0 0 0; 0 1 0; 0 0 0]
    c = inject(cube, kernel, angles; x=11, y=11)
    @test c isa AnnulusView
    @test all(c[11, 11, :] .≈ 1)

    @testset "synthetic psf - $mod" for mod in [gaussian, moffat, airydisk]
        c = inject(cube, mod, angles; x=11, y=11, fwhm=3)
        @test c isa AnnulusView
        @test all(c[11, 11, :] .≈ 1)
    end
end

@testset "inject - MultiAnnulusView" begin
    cube = MultiAnnulusView(zeros(21, 21, 10), 5)
    angles = 90 .* ones(10)

    kernel = [0 0 0; 0 1 0; 0 0 0]
    c = inject(cube, kernel, angles; x=11, y=11)
    @test c isa MultiAnnulusView
    @test all(c[11, 11, :] .≈ 1)

    @testset "synthetic psf - $mod" for mod in [gaussian, moffat, airydisk]
        c = inject(cube, mod, angles; x=11, y=11, fwhm=3)
        @test c isa MultiAnnulusView
        @test all(c[11, 11, :] .≈ 1)
    end
end


text:
"""
    setline(linewidth)

Set the line width and multiply it with the current multiplier
which is i.e. set by [`appear`](@ref) and [`disappear`](@ref).

Normal behavior without any animation is the same as `Luxor.setline`.

# Example
```
setline(10)
line(O, Point(10, 10))
```

# Arguments:
- `linewidth`: the line width in pixel
"""
function setline(linewidth)
    action = CURRENT_ACTION[1]
    cs  = action.current_setting
    cs.line_width  = linewidth
    current_line_width = cs.line_width * cs.mul_line_width
    Luxor.setline(current_line_width)
end

"""
    setopacity(linewidth)

Set the opacity and multiply it with the current multiplier
which is i.e. set by [`appear`](@ref) and [`disappear`](@ref).

Normal behavior without any animation is the same as `Luxor.setopacity`.

# Example
```
setopacity(0.5)
circle(O, 20, :fill)
```

# Arguments:
- `opacity`: the opacity between 0.0 and 1.0
"""
function setopacity(opacity)
    action = CURRENT_ACTION[1]
    cs  = action.current_setting
    cs.opacity  = opacity
    current_opacity = cs.opacity * cs.mul_opacity
    Luxor.setopacity(current_opacity)
end


text:
###########################################################
# Sylvester based preconditioner for the Waveguide eigenvalue problem
# Implementing the Preconditioner in the optimized way described in Ringh et al.

"""
    WEP_preconditioner
 A wrapper struct to be able to use the preconditioner. Update in GMRES for Julia v0.6 remove possibility to use a function directly.
"""
    struct WEP_preconditioner
        scratch_pad_for_FFT::Matrix{ComplexF64}
        scratch_pad_for_transpose::Matrix{ComplexF64}
        scratch_pad_for_Z::Matrix{ComplexF64}
        nep::WEP_FD
        M
        σ
        function WEP_preconditioner(nep::NEP, M, σ)
            nx = nep.nx
            nz = nep.nz
            A = zeros(ComplexF64, 2*(nx+1), nz)
            B = zeros(ComplexF64, nx, nz)
            C = zeros(ComplexF64, nz, nx)
            return new(A, B, C, nep, M, σ)
        end
    end

    function ldiv!(A::WEP_preconditioner, B)
        C = reshape(B, A.nep.nz, A.nep.nx)
        B[:] = vec( solve_smw( A.nep, A.M, C, A.σ, A.scratch_pad_for_FFT, A.scratch_pad_for_transpose, A.scratch_pad_for_Z) )
    end

# Ringh - Algorithm 2, step 10
"""
    wep_generate_preconditioner(nep::WEP_FD, N::Integer, σ)
 Given a nep of type WEP_FD, the number of domains in z-direction N, and a fixed shift σ,\\
 this computes a functor (struct acting like a function) that acts as a Preconditioner for the WEP.
"""
    function wep_generate_preconditioner(nep::WEP_FD, N::Integer, σ)
        M = generate_smw_matrix(nep, N, σ)
        return WEP_preconditioner(nep, M, σ)
    end

"""
    generate_smw_matrix( nep::WEP_FD, N::Integer, σ)
 Given a nep of type WEP_FD, this computes the Sylvester-SMW matrix with N domains in z-direction\\
 and for fixed shift σ.
"""
    function generate_smw_matrix(nep::WEP_FD, N::Integer, σ)
        # OBS: n = nz, and nz = nx + 4;
        if( (nep.nz+4) != nep.nx)
            error("This implementation requires nx = nz + 4. Provided NEP has nz = ", nep.nz, " and nx = ", nep.nx)
        end
        if( !isinteger(nep.nz/N) )
            error("This implementation is uniform in the blocking and therefore requires nz/N to be an integer. Provided data is nz = ", nep.nz, " with N = ", N, " and hence nz/N = ", nep.nz/N, " which is not deemed to be numerically equal to an integer." )
        end

        nz::Integer = nep.nz
        nx::Integer = nep.nx

        dd1 = nep.d1/nep.hx^2;
        dd2 = nep.d2/nep.hx^2;

        # Sylvester solver
        scratch_pad_for_FFT::Matrix{ComplexF64} = zeros(ComplexF64, 2*(nx+1), nz)
        scratch_pad_for_transpose::Matrix{ComplexF64} = zeros(ComplexF64, nx, nz)
        scratch_pad_for_Z::Matrix{ComplexF64} = zeros(ComplexF64, nz, nx)
        Linv! = function(rhs)
            return solve_wg_sylvester_fft!(rhs, σ, nep.k_bar, nep.hx, nep.hz, scratch_pad_for_FFT, scratch_pad_for_transpose, scratch_pad_for_Z)
        end

        # OBS: MINUS sign as in Ringh - (4.10)
        Pm(v) = -P_inv_m(nep, σ, v)
        Pp(v) = -P_inv_p(nep, σ, v)

        return generate_smw_matrix(nz, N, Linv!, dd1, dd2, Pm, Pp, nep.K )
    end

"""
    solve_smw(nep::WEP_FD, M, C, σ, scratch_pad_for_FFT, scratch_pad_for_transpose, scratch_pad_for_Z)
 Given a nep of type WEP_FD, an SMW-system matrix M computed with shift σ, and a right hand side C,\\
 This computes the solution to the SMW-matrix equation.
"""
    solve_smw(nep::WEP_FD, M, C, σ) = solve_smw(nep::WEP_FD, M, C, σ, zeros(ComplexF64, 2*(size(C,2)+1), size(C,1)), zeros(ComplexF64, size(C,2), size(C,1)), zeros(ComplexF64, size(C,1), size(C,2)))

    function solve_smw(nep::WEP_FD, M, C, σ, scratch_pad_for_FFT, scratch_pad_for_transpose, scratch_pad_for_Z)

        C_copy::Matrix{ComplexF64} = copy(C) #Make sure it is complex since that is how FFT works. Also take copy since WG_FFT solver works in place.

        nz::Integer = nep.nz
        nx::Integer = nep.nx

        dd1 = nep.d1/nep.hx^2;
        dd2 = nep.d2/nep.hx^2;

        # Sylvester solver
        Linv! = function(rhs)
            return solve_wg_sylvester_fft!(rhs, σ, nep.k_bar, nep.hx, nep.hz, scratch_pad_for_FFT, scratch_pad_for_transpose, scratch_pad_for_Z)
        end

        # OBS: MINUS sign as in Ringh - (4.10)
        Pm(v) = -P_inv_m(nep, σ, v)
        Pp(v) = -P_inv_p(nep, σ, v)

        return solve_smw(M, C_copy, Linv!, dd1, dd2, Pm, Pp, nep.K)
    end

#FFT diagonalization and solution to WEP matrix equation.
#Ringh - Section 5.3
"""
    solve_wg_sylvester_fft( C, λ, k_bar, hx, hz, scratch_pad_for_FFT, scratch_pad_for_transpose )
 Solves the Sylvester equation for the WEP, with C as right hand side.
 OBS: Writes the solution to the varaible C.
 Last three arguments (scratch pad:s) are optional and there to allow reuse of memory allocation. They will be overwritten in the process!
"""
solve_wg_sylvester_fft!( C, λ, k_bar, hx, hz) = solve_wg_sylvester_fft!( C, λ, k_bar, hx, hz, zeros(ComplexF64, 2*(size(C,2)+1), size(C,1)), zeros(ComplexF64, size(C,2), size(C,1)), zeros(ComplexF64, size(C,1), size(C,2)))

function solve_wg_sylvester_fft!( C, λ, k_bar, hx, hz, scratch_pad_for_FFT, scratch_pad_for_transpose, scratch_pad_for_Z)

    nz = size(C,1)
    nx = size(C,2)

    alpha = λ^2+k_bar;

    #eigenvalues of A
    v=zeros(ComplexF64, nz);   v[1]=-2;    v[2]=1;     v[nz]=1;   v=v/(hz^2);
    w=zeros(ComplexF64, nz);   w[2]=1;     w[nz]=-1;   w=w*(λ/hz);
    D = fft(v+w) .+ alpha

    # eigenvalues of B = Dxx
    S = -((4+0.0im)/hx^2) * sin.(pi*(1:nx)/(2*(nx+1))).^2
#    S=S.'

    # change variables
    adjoint!(scratch_pad_for_transpose, C)  # scratch_pad_for_transpose = C'
    adjoint!(C, Wh(scratch_pad_for_transpose, scratch_pad_for_FFT )) # C = (Wh(C'))'
    Vh!( C )  #In effect: C = Vh( Wh(C')' )

    # solve the diagonal matrix equation
    scratch_pad_for_Z[:,:] .= 0.0im
    for k=1:nx
        scratch_pad_for_Z[:,k] += C[:,k] ./ (D .+ S[k])
    end

    # change variables
    adjoint!(scratch_pad_for_transpose, scratch_pad_for_Z)  # scratch_pad_for_transpose = Z'
    adjoint!(C, W(scratch_pad_for_transpose, scratch_pad_for_FFT ))  #C = (W(Z'))'
    V!( C ) #In effect: C = V( W(Z')' )

    return C

end

# Start: Auxiliary computations of eigenvector actions using FFT
#Connects to the solving of the Sylvester equation by FFT-diagonalization
#Ringh - Section 5.3
    function V!(X::Matrix{ComplexF64})
    # Compute the action of the eigenvectors of A = Dzz + Dz + c*I
        nx::ComplexF64 = size(X,2)
        fft!(X,1)#/sqrt(nx)
        rmul!(X, (1+0.0im)/sqrt(nx))
    end

    function Vh!(X::Matrix{ComplexF64})
    # Compute the action of the transpose of the eigenvectors of A = Dzz + Dz + c*I
        nx::ComplexF64 = size(X,2)
        ifft!(X,1)#* sqrt(nx)
        rmul!(X, (1+0.0im)*sqrt(nx))
    end

    function W(X::Matrix{ComplexF64}, scratch_pad::Matrix{ComplexF64})
    #W Compute the action of the matrix W
    #   W is the matrix of the eigenvectors of the second derivative Dxx
    #   W*X can be computed with FFTs

        nz::ComplexF64 = size(X,1)
        WX::Matrix{ComplexF64} = F(X,scratch_pad)-Fh(X,scratch_pad)

        rmul!(WX, (1.0im/2.0) * 1/sqrt((nz+1)/2.0) )
        return WX

    end

    function Wh(X::Matrix{ComplexF64}, scratch_pad::Matrix{ComplexF64})
    #Wh Compute the action of the matrix Wh
    #   Wh is the transpose of the matrix of the eigenvectors of the second derivative Dxx
    #   Wh*X can be computed with FFTs

        nz::ComplexF64 = size(X,1)
        WX::Matrix{ComplexF64} = F(X,scratch_pad)-Fh(X,scratch_pad)

        rmul!(WX, (1.0im/2.0) * 1/sqrt((nz+1)/2.0) )
        return WX

    end

    function F(v::Matrix{ComplexF64}, scratch_pad::Matrix{ComplexF64})
    #F is an auxiliary function for W and Wh

        m=size(v,2)
        n=size(v,1) + 1

        scratch_pad[:,:] .= 0.0im
        scratch_pad[2:n,:] = v
        fft!(scratch_pad,1)
        return scratch_pad[2:n,:]
    end

    function Fh( v::Matrix{ComplexF64}, scratch_pad::Matrix{ComplexF64})
    #Fh is an auxiliary function for W and Wh

        m=size(v,2)
        n=size(v,1) + 1

        scratch_pad[:,:] .= 0.0im
        scratch_pad[2:n,:] = v
        ifft!(scratch_pad,1)
        return scratch_pad[2:n,:]*2*n

    end
# End: Auxiliary computations of eigenvector actions using FFT

###########################################################
# Sylvester SMW
# Ringh - Section 4
"""
    generate_smw_matrix(n::Integer, N::Integer, Linv!::Function, dd1, dd2, Pm, Pp, K)
 Computes the SMW matrix for the Sylvester SMW on rectangular domains, with n points in z-direction, n+4 points in x-direction,\\
 N domains in z-direction, N+4 domains in x-direction, and fixed shift.\\
 Obs: Linv! works in place on the matrix rhs
"""
function generate_smw_matrix(n::Integer, N::Integer, Linv!::Function, dd1, dd2, Pm::Function, Pp::Function, K::Union{Matrix{ComplexF64}, Matrix{Float64}})

    # OBS: n = nz, and nz = nx + 4
    nz::Integer = n
    nx::Integer = n + 4

    L::Integer = n/N             # Number of points in one dimanesion of the regions
    LL::Integer = L*L            # Number of points in the "square interior regions"
    mm::Integer = (N^2 + 4*N)    # Number of elements in SMW-matrix

    # block extract index
    II = function (i::Integer) # z-direction, always equal
        return (i-1)*L+1:i*L
    end
    JJ = function(j::Integer) #x-direction, different if boundary or interior
        return ((j-3)*L+1:(j-2)*L) .+ 2
    end
    JJ_2 = function(j::Integer)
        return (j==1)*1 + (j==2)*2 + (j==N+3)*(n+3) + (j==N+4)*(n+4)
    end

    # convert a single index k to two indeces
    # reading the matrix X by rows; left to right and top to down
    k2ij = function(k::Integer)
        j::Integer = rem(k,N+4)+(rem(k,N+4)==0)*(N+4)
        i::Integer = (k-j)/(N+4)+1
        return(i,j)
    end

    # compute matrix M
    M::Matrix{ComplexF64} = zeros(ComplexF64, mm, mm)

    EEk::Matrix{ComplexF64} = zeros(ComplexF64, nz, nx)
    ek::Vector{ComplexF64} = zeros(ComplexF64, nz)

    for k=1:mm

        i,j = k2ij(k)

        # Ek tilde
        EEk[:,:] .= 0.0im
        if j==1
            EEk[II(i), JJ_2(j)] = K[II(i), JJ_2(j)]
            ek[:] .= 0.0im
            ek[II(i)] .= dd1
            EEk[:, 1] += Pm(ek)
        elseif j==2
            EEk[II(i), JJ_2(j)] = K[II(i), JJ_2(j)]
            ek[:] .= 0.0im
            ek[II(i)] .= dd2
            EEk[:, 1] += Pm(ek)
        elseif j==N+4
            EEk[II(i),JJ_2(j)] = K[II(i),JJ_2(j)]
            ek[:] .= 0.0im
            ek[II(i)] .= dd1
            EEk[:, nx] += Pp(ek)
        elseif j==N+3
            EEk[II(i), JJ_2(j)] = K[II(i), JJ_2(j)]
            ek[:] .= 0.0im
            ek[II(i)] .= dd2
            EEk[:, nx] += Pp(ek)
        else
            EEk[II(i), JJ(j)] = K[II(i), JJ(j)]
        end

        # Sylvester solve of E tilde
        Linv!(EEk) # In effect EEk = Linv!(EEk) = Fk

        # Build this matrix element

        for kk=1:mm
            i,j = k2ij(kk)
            # evaluate the linear functional
            if((j==1)||(j==2)||(j==N+3)||(j==N+4))
                M[kk, k] = sum(sum(EEk[II(i), JJ_2(j)] ))/L
            else
                M[kk, k] = sum(sum(EEk[II(i), JJ(j)] ))/LL
            end
        end
    end

    M += Matrix{ComplexF64}(I, mm, mm)

    return factorize(M)

end

"""
    solve_smw( M, C, Linv!::Function, dd1, dd2, Pm, Pp, K)
 Solves the matrix equation SMW system and computes the solution, on rectangular domains.\n
 With SMW-system matrix M, right hand side C, and matrix equation solver Linv! which was used to compute M.
 Obs: Linv! works in place on the matrix rhs
"""
function solve_smw( M, C::Matrix{ComplexF64}, Linv!::Function, dd1, dd2, Pm::Function, Pp::Function, K::Union{Matrix{ComplexF64}, Matrix{Float64}})

    mm::Integer = size(M,1)
    N::Integer = sqrt(mm+4)-2      #OBS: N^2 + 4N = length(M)

    nz::Integer = size(C,1)
    nx::Integer = size(C,2)
    n::Integer = nz
    L::Integer = nz/N
    LL::Integer = L*L

    # block extract index
    II = function (i::Integer) # z-direction, always equal
        return (i-1)*L+1:i*L
    end
    JJ = function(j::Integer) #x-direction, different if boundary or interior
        return ((j-3)*L+1:(j-2)*L) .+ 2
    end
    JJ_2 = function(j::Integer)
        return (j==1)*1 + (j==2)*2 + (j==N+3)*(n+3) + (j==N+4)*(n+4)
    end

    # convert a single index k to two indeces
    # reading the matrix X by rows; left to right and top to down
    k2ij = function(k::Integer)
        j::Integer = rem(k,N+4)+(rem(k,N+4)==0)*(N+4)
        i::Integer = (k-j)/(N+4)+1
        return(i,j)
    end

    # compute the right hand side
    Linv!(C) # Now C = LinvC

    b=zeros(ComplexF64,mm)
    for k=1:mm
        i, j = k2ij(k)
        # evaluate the linear functional
        if((j==1)||(j==2)||(j==N+3)||(j==N+4))
            b[k] = sum(sum( C[II(i), JJ_2(j)] ))/L
        else
            b[k] = sum(sum( C[II(i), JJ(j)] ))/LL
        end
    end

    # solve for the coefficients
    alpha = M\b;

    # build the solution
    Y::Matrix{ComplexF64} = zeros(ComplexF64, nz, nx);
    ek::Vector{Float64} = zeros(ComplexF64, nz);
    for k=1:mm

        i, j = k2ij(k)

        # Ek tilde implicitly created
        if j==1
            Y[II(i), 1] += alpha[k]*K[II(i), 1]
            ek[:] .= 0.0im
            ek[II(i)] .= dd1
            Y[:, 1] += alpha[k]*Pm(ek)
        elseif j==2
            Y[II(i), 2] += Y[II(i),2] + alpha[k]*K[II(i), 2]
            ek[:] .= 0.0im
            ek[II(i)] .= dd2
            Y[:, 1] += alpha[k]*Pm(ek)
        elseif j==N+4
            Y[II(i), nx] += alpha[k]*K[II(i), nx]
            ek[:] .= 0.0im
            ek[II(i)] .= dd1
            Y[:, nx] += alpha[k]*Pp(ek)
        elseif j==N+3
            Y[II(i), nx-1] += alpha[k]*K[II(i), nx-1]
            ek[:] .= 0.0im
            ek[II(i)] .= dd2
            Y[:, nx] += alpha[k]*Pp(ek)
        else
            Y[II(i),JJ(j)] += alpha[k]*K[II(i), JJ(j)]
        end

    end
    Linv!(Y) # Now Y = LinvY

    return C-Y # Which is LinvC - LinvY

end