text:
function find_surface(r0, dr, psi, psi0)

    r = r0[:]
    sgn = psi(r0...) - psi0
    while true
        rnew += dr
        if sgn*(psi(rnew...) - psi0) < 0
            return 0.5*(rnew + r)
        end
        r[:] = rnew
    end
end

function rootfind_bisection(f::Function, val::Float64, bracket1::Vector{Float64}, bracket2::Vector{Float64})

    a = bracket1[:]
    b = bracket2[:]
    fa = f(a) - val
    fb = f(b) - val
    c = (a + b)/2
    # c = (a*fb - b*fa) / (fb - fa)
    fc = f(c) - val

    while abs(fc) > 1e-8
        if fa*fc < 0
            b = c
        elseif fc*fb < 0
            a = c
        else
            @error "root finding algorithm has failed."
            break
        end
        fa = f(a) - val
        fb = f(b) - val
        c = (a + b)/2
        # c = (a*fb - b*fa) / (fb - fa)
        fc = f(c) - val
    end

    return c
end

function rk4_step(r::Vector{Float64}, f::Function, ds::Float64)

    k1 = f(r)
    k2 = f(r + ds*k1/2)
    k3 = f(r + ds*k2/2)
    k4 = f(r + ds*k3)

    return r + ds*(k1 + 2*k2 + 2*k3 + k4) / 6
end

function conservative_step(r::Vector{Float64}, b::Function, psi::Function, ds::Float64)

    # compute predictor step using rk4
    rp = rk4_step(r, b, ds)

    # compute correction using root finding algorithm
    en = zeros(length(rp))
    en[1:2] = [0 -1; 1 0] * b(rp)[1:2] / norm(b(rp)[1:2])
    bracket1 = rp + ds*en
    bracket2 = rp - ds*en
    rc = rootfind_bisection(psi, psi(r), bracket1, bracket2)

    # compute with correct displacement
    return rc, norm(rc - r)
end

function trace_fieldline(stop_condition::Function, r0::Vector{Float64}, b::Function, psi::Function, ds::Float64)

    r = r0[:]
    deltaS = 0.0

    prog = ProgressUnknown("Tracing fieldline... ", spinner=true)
    points = Vector{Float64}[r0]
    while !stop_condition(r)
        rnew, dsnew = conservative_step(r, b, psi, ds)
        r[:] = rnew
        deltaS += dsnew
        append!(points, [rnew])
        next!(prog)
    end

    finish!(prog)
    return hcat(points...), deltaS
end

function trace_fieldline(N::Int64, r0::Vector{Float64}, b::Function, psi::Function, ds::Float64)

    r = r0[:]
    deltaS = 0.0

    points = zeros(length(r0), N)
    points[:,1] = r[:]
    for s=2:N
        rnew, dsnew = conservative_step(r, b, psi, ds)
        r[:] = rnew
        deltaS += dsnew
        points[:,s] = r
    end

    return points, deltaS
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BezierCurve(points)

A recursive Bézier curve with control points `points`.
See https://en.wikipedia.org/wiki/Bézier_curve.
A point on the curve `b` can be evaluated by calling
`b(t)` with `t` between `0` and `1`.
The evaluation method defaults to DeCasteljau's algorithm
for accurate evaluation. Horner's method, faster with a
large number of points but less precise, can be used via
`b(t, Horner())`.

## Examples

```julia
BezierCurve(Point2[(0.,0.),(1.,-1.)])
```
"""
struct BezierCurve{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Primitive{Dim,T}
  controls::V
end

BezierCurve(points::AbstractVector{<:Tuple}) = BezierCurve(Point.(points))
BezierCurve(points::Vararg) = BezierCurve(collect(points))

ncontrols(b::BezierCurve) = length(b.controls)

degree(b::BezierCurve) = ncontrols(b) - 1

"""
Evaluation method used to obtain a point along
a Bézier curve from a parametric expression.
"""
abstract type BezierEvalMethod end

"""
Accurate evaluation using De Casteljau's algorithm.
Recommended for a small number of control points.

See https://en.wikipedia.org/wiki/De_Casteljau%27s_algorithm.
"""
struct DeCasteljau <: BezierEvalMethod end

"""
Approximate evaluation using Horner's method.
Recommended for a large number of control points,
if you can afford a precision loss.

See https://en.wikipedia.org/wiki/Horner%27s_method.
"""
struct Horner <: BezierEvalMethod end

(curve::BezierCurve)(t) = curve(t, DeCasteljau())

function (curve::BezierCurve)(t, ::DeCasteljau)
  if t < 0 || t > 1
    throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
  end
  ss = segments(Chain(curve.controls))
  points = [s(t) for s in ss]
  if length(points) == 1
    points[1]
  else
    BezierCurve(points)(t)
  end
end

"""
Apply Horner's method on the monomial representation of the
Bézier curve B = ∑ᵢ aᵢtⁱ with i ∈ [0, n], n the degree of the
curve, aᵢ = binomial(n, i) * pᵢ * t̄ⁿ⁻ⁱ and t̄ = (1 - t).
Horner's rule recursively reconstructs B from a sequence bᵢ
with bₙ = aₙ and bᵢ₋₁ = aᵢ₋₁ + bᵢ * t until b₀ = B.
"""
function (curve::BezierCurve{Dim,T})(t, ::Horner) where {Dim,T}
  if t < 0 || t > 1
    throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
  end
  cs = curve.controls
  t̄ = one(T) - t
  n = degree(curve)
  pₙ = coordinates(last(cs))
  aₙ = pₙ

  # initialization with i = n + 1, so bᵢ₋₁ = bₙ = aₙ
  bᵢ₋₁ = aₙ
  cᵢ₋₁ = one(T)
  t̄ⁿ⁻ⁱ = one(T)
  for i in n:-1:1
    cᵢ₋₁ *= i / (n - i + one(T))
    pᵢ₋₁ = coordinates(cs[i])
    t̄ⁿ⁻ⁱ *= t̄
    aᵢ₋₁ = cᵢ₋₁ * pᵢ₋₁ * t̄ⁿ⁻ⁱ
    bᵢ₋₁ = aᵢ₋₁ + bᵢ₋₁ * t
  end

  b₀ = bᵢ₋₁
  Point(b₀)
end

function boundary(curve::BezierCurve)
  p = curve.controls
  a = first(p)
  b = last(p)
  a == b ? nothing : PointSet([a, b])
end

ext: jl


text:
# Functional equations of state

using JLD  # for accessing data in .jld files (HDF5-based)
using Roots: fzero  # for numerical inversion
using QuadGK: quadgk  # for numerical integration

# simple EOSes
export ConstantEOS, PolytropicEOS, BME, BME3, BME4, Vinet, TFD
# temperature-dependent EOSes
export ChoukrounGrasset, IAPWS
# EOS modifiers
export BoundedEOS, MGDPressureEOS

# Shared variables

"Default range of ρ (in kg/m^3) to consider when doing numerical inversion of ρ(P)"
const inversion_ρ_range = [1e-3, 1e8]

# Functional EOS type hierarchy

"An EOS which calculates ρ by evaluating a function ρ = ρ(P[, T])"
abstract type FunctionalEOS <: EOS end
"An EOS which calculates ρ by numerically inverting a function P = P(ρ[, T])"
abstract type InverseFunctionalEOS <: FunctionalEOS end

# FIXME: this is a workaround for julia issue #14919
"Declare that an EOS type is numerically inverted when called"
macro inverseEOS(theType)
    return quote
        function (eos::$(esc(theType)))(P)
            fzero(ρ -> pressure(eos, ρ) - P, ρmin(eos), ρmax(eos))
        end

        function (eos::$(esc(theType)))(P, T)
            fzero(ρ -> pressure(eos, ρ, T) - P, ρmin(eos), ρmax(eos))
        end
    end
end

# Specific EOSes

"The ideal gas equation of state, where R is the specific gas constant"
struct IdealGas <: FunctionalEOS
    R::Float64
end

"Equation of state for a constant-density substance"
struct ConstantEOS <: FunctionalEOS
    ρ::Float64
end

"Choukroun and Grasset's cold ice EOS"
struct ChoukrounGrasset <: FunctionalEOS
    Pref::Float64
    Tref::Float64
    V₀::Float64
    aT::Vector{Float64}
    P₀::Float64
    aP::Vector{Float64}
end
function ChoukrounGrasset(Pref, Tref, V₀, aT₁, aT₂, aT₃, aT₄, P₀, aP₁, aP₂, aP₃, aP₄)
    ChoukrounGrasset(Pref, Tref, V₀, [aT₁, aT₂, aT₃, aT₄], P₀, [aP₁, aP₂, aP₃, aP₄])
end

function get_choukroungrasset_pars(phase)
    readtable("$(config.rawdata)/ChoukrounGrasset-parameters.dat")[phase]
end

"The polytropic EOS, ρ(P) = ρ₀ + aρ^n"
struct PolytropicEOS <: FunctionalEOS
    ρ₀::Float64
    a::Float64
    n::Float64
end

"The Birch-Murnaghan EOS"
abstract type BME <: InverseFunctionalEOS end

"Third-order Birch-Murnaghan EOS"
struct BME3 <: BME
    ρ₀::Float64
    K₀::Float64
    dK₀::Float64
    ρmin::Float64
    ρmax::Float64
end
@inverseEOS BME3
BME3(ρ₀, K₀, dK₀) = BME3(ρ₀, K₀, dK₀, inversion_ρ_range...)
BME(ρ₀, K₀, dK₀) = BME3(ρ₀, K₀, dK₀)

"Fourth-order Birch-Murnaghan EOS"
struct BME4 <: BME
    ρ₀::Float64
    K₀::Float64
    dK₀::Float64
    d2K₀::Float64
    ρmin::Float64
    ρmax::Float64
end
BME4(ρ₀, K₀, dK₀, d2K₀) = BME4(ρ₀, K₀, dK₀, d2K₀, inversion_ρ_range...)
BME(ρ₀, K₀, dK₀, d2K₀) = BME4(ρ₀, K₀, dK₀, d2K₀)
@inverseEOS BME4

"The Vinet EOS"
struct Vinet <: InverseFunctionalEOS
    ρ₀::Float64
    K₀::Float64
    dK₀::Float64
    ρmin::Float64
    ρmax::Float64
end
Vinet(ρ₀, K₀, dK₀) = Vinet(ρ₀, K₀, dK₀, inversion_ρ_range...)
@inverseEOS Vinet

"The Thomas-Fermi-Dirac EOS"
struct TFD <: FunctionalEOS
    Z::Vector{Int}
    A::Vector{Float64}
    n::Vector{Float64}

    function TFD(Z::Vector{Int}, A::Vector{Float64}, n::Vector{Float64})
        @assert length(Z) == length(A) == length(n)
        new(Z, A, n)
    end
end
# default constructor
function TFD(Z::Vector{Int}, A::Vector{<:Real}, n=ones(length(Z)))
    TFD(Z, Float64[A], n)
end
# coerce to arrays if necessary
function TFD(Z::Integer, A)
    TFD(Int[Z], [A])
end

"The IAPWS EOS, functional formulation"
struct IAPWS <: InverseFunctionalEOS
    c::Vector{Float64}
    d::Vector{Float64}
    t::Vector{Float64}
    n::Vector{Float64}
    α::Vector{Float64}
    β::Vector{Float64}
    γ::Vector{Float64}
    ϵ::Vector{Float64}
    a::Vector{Float64}
    b::Vector{Float64}
    B::Vector{Float64}
    C::Vector{Float64}
    D::Vector{Float64}
    A::Vector{Float64}
    ρmin::Float64
    ρmax::Float64
end
function IAPWS(ρmin, ρmax)
    coeffs = _get_iapwscoeffs()
    IAPWS([coeffs[:, i] for i=1:14]..., ρmin, ρmax)
end
function _get_iapwscoeffs()
    table = readdlm("$(config.rawdata)/IAPWS-coeffs.dat", ',', Float64, skipstart=1)
    coeffs = table[:, 2:end]
    @assert size(coeffs) == (56, 14)

    coeffs
end
@inverseEOS IAPWS

"Wrapper for indicating that some `EOS` is bounded"
struct BoundedEOS{E<:EOS} <: EOS
    eos::E
    extent::Region
end
BoundingBox(b::BoundedEOS) = BoundingBox(b.extent)

"Wrapper for including thermal pressure with a BME or Vinet EOS"
struct MGDPressureEOS <: InverseFunctionalEOS
    eos::Union{BME, Vinet}
    T₀::Float64
    θD₀::Float64
    γ₀::Float64
    q::Int
    n::Int
end
@inverseEOS MGDPressureEOS

# Helper functions for EOS evaluation

# supporting functions for IAPWS density/temperature
"Get density from dimensionless density"
ρ_from_δ(δ) = δ * ρc
"Get dimensionless density from density"
δ_from_ρ(ρ) = ρ / ρc
"Get temperature from dimensionless temperature"
T_from_τ(τ) = Tc / τ
"Get dimensionless temperature from temperature"
τ_from_T(T) = Tc / T

# supporting functions for Helmholtz energy
# residual parts (not used in this code but kept here for testing purposes)
function ϕr1(I::IAPWS, δ, τ, i)
    let n = I.n[i], d = I.d[i], t = I.t[i]

        n * δ^d * τ^t
    end
end
function ϕr2(I::IAPWS, δ, τ, i)
    let n = I.n[i], d = I.d[i], t = I.t[i], c = I.c[i]

        n * δ^d * τ^t * exp(-δ^c)
    end
end
function ϕr3(I::IAPWS, δ, τ, i)
    let n = I.n[i], d = I.d[i], t = I.t[i], α = I.α[i],
        ϵ = I.ϵ[i], β = I.β[i], γ = I.γ[i]

        n * δ^d * τ^t * exp(-α*(δ-ϵ)^2 - β*(τ-γ)^2)
    end
end
function ϕr4(I::IAPWS, δ, τ, i)
    let n = I.n[i], Δ = Δ(I, δ, τ, i), b = I.b[i], ψ = ψ(I, δ, τ, i)

        n * Δ^b * δ * ψ
    end
end

# helper functions for residual parts
function ψ(I::IAPWS, δ, τ, i)
    let C = I.C[i], D = I.D[i]
        exp(-C*(δ-1)^2 - D*(τ-1)^2)
    end
end
function θ(I::IAPWS, δ, τ, i)
    let A = I.A[i], β = I.β[i]
        (1-τ) + A*((δ-1)^2)^(1/2β)
    end
end
function Δ(I::IAPWS, δ, τ, i)
    let θ = θ(I, δ, τ, i), B = I.B[i], a = I.a[i]

        θ^2 + B*((δ-1)^2)^a
    end
end

# derivatives of residual parts
function dϕr1(I::IAPWS, δ, τ, i)
    let n = I.n[i], d = I.d[i], t = I.t[i]

        n * d * δ^(d-1) * τ^t
    end
end
function dϕr2(I::IAPWS, δ, τ, i)
    let n = I.n[i], d = I.d[i], t = I.t[i], c = I.c[i]

        n * exp(-δ^c) * (δ^(d-1) * τ^t * (d - c*δ^c))
    end
end
function dϕr3(I::IAPWS, δ, τ, i)
    let n = I.n[i], d = I.d[i], t = I.t[i], α = I.α[i],
        ϵ = I.ϵ[i], β = I.β[i], γ = I.γ[i]

        (n * δ^d * τ^t * exp(-α*(δ-ϵ)^2 - β*(τ-γ)^2)
         * (d/δ - 2α*(δ - ϵ)))
    end
end
function dϕr4(I::IAPWS, δ, τ, i)
    let n = I.n[i], Δ = Δ(I, δ, τ, i), b = I.b[i], ψ = ψ(I, δ, τ, i),
        dψdδ = dψdδ(I, δ, τ, i), dΔbdδ = dΔbdδ(I, δ, τ, i)

        n * (Δ^b * (ψ + δ*dψdδ) + dΔbdδ * δ * ψ)
    end
end

# derivatives of helper functions
function dψdδ(I::IAPWS, δ, τ, i)
    let C = I.C[i], ψ = ψ(I, δ, τ, i)

        -2C * (δ - 1) * ψ
    end
end
function dΔbdδ(I::IAPWS, δ, τ, i)
    let b = I.b[i], Δ = Δ(I, δ, τ, i), dΔdδ = dΔdδ(I, δ, τ, i)

        b * Δ^(b - 1) * dΔdδ
    end
end
function dΔdδ(I::IAPWS, δ, τ, i)
    let A = I.A[i], θ = θ(I, δ, τ, i), β = I.β[i], a = I.a[i], B = I.B[i]

        ((δ - 1) * (A * θ * 2/β * ((δ - 1)^2)^(1/2β - 1)
                    + 2B * a * ((δ - 1)^2)^(a - 1)))
    end
end

"IAPWS Helmholtz energy"
function ϕr(I::IAPWS, δ, τ)
    # we take sums over the functions dϕ1, dϕ2...
    part1 = mapreduce(i -> ϕr1(I, δ, τ, i), (+), 0, 1:7)
    part2 = mapreduce(i -> ϕr2(I, δ, τ, i), (+), 0, 8:51)
    part3 = mapreduce(i -> ϕr3(I, δ, τ, i), (+), 0, 52:54)
    part4 = mapreduce(i -> ϕr4(I, δ, τ, i), (+), 0, 55:56)

    # and then add them all together
    (part1 + part2 + part3 + part4)
end

"Derivative of the IAPWS Helmholtz energy in the density direction"
function dϕrδ(I::IAPWS, δ, τ)
    # we take sums over the functions dϕ1, dϕ2...
    part1 = mapreduce(i -> dϕr1(I, δ, τ, i), (+), 0, 1:7)
    part2 = mapreduce(i -> dϕr2(I, δ, τ, i), (+), 0, 8:51)
    part3 = mapreduce(i -> dϕr3(I, δ, τ, i), (+), 0, 52:54)
    part4 = mapreduce(i -> dϕr4(I, δ, τ, i), (+), 0, 55:56)

    # and then add them all together
    (part1 + part2 + part3 + part4)
end

"Pressure component of the Choukroun-Gras￼set ice EOS"
function ξP(cg::ChoukrounGrasset, P)
    let a = cg.aP, P₀ = cg.P₀
        a[1] + a[2]*(P - P₀) + a[3]*(P - P₀)^a[4]
    end
end

"Temperature component of the Choukroun-Grasset ice EOS"
function ξT(cg::ChoukrounGrasset, T)
    let a = cg.aT, Tref = cg.Tref
        a[1] + a[2]*(T - Tref) + a[3]*(T - Tref)^a[4]
    end
end

"Specific volume of the Choukroun-Grasset ice EOS at a given pressure and temperature"
specificvolume(cg::ChoukrounGrasset, P, T) = cg.V₀ * ξP(cg, P) * ξT(cg, T)

"Integrand in the MGD expression"
mgd_integrand(t) = t^3 / (exp(t) - 1)

"Vibrational energy in the MGD expression [J / mol]"
mgd_Evib(T, n, θD) = 9n*R*T*(T/θD)^3 * quadgk(mgd_integrand, 0, θD/T)[1]

"Thermal pressure component of the MGD"
function thermalpressure(mgd::MGDPressureEOS, ρ, T)
    let γ₀ = mgd.γ₀, θD₀ = mgd.θD₀, ρ₀ = mgd.eos.ρ₀, q = mgd.q,
        M = h2o_molar_mass, n = mgd.n, T₀ = mgd.T₀

        γ = γ₀ * (ρ/ρ₀)^(-q)
        θD = θD₀ * exp((γ₀ - γ)/q)

        ΔP = γ * ρ / M * (mgd_Evib(T, n, θD) - mgd_Evib(T₀, n, θD))
    end
end

const _tfd_g_coeffs = [0          0          0         0         0;
                       0          0          0         0         0;
                       1.512E-2   8.955E-2   1.090E-1  5.089     -5.980;
                       2.181E-3   -4.015E-1  1.698     -9.566    9.873;
                       -3.328E-4  5.167E-1   -2.369    1.349E1   -1.427E1;
                       -1.384E-2  -6.520E-1  3.529     -2.095E1  2.264E1]

"Helper function for TFD partial sum"
function β_(n::Integer, ϵ::Vector{Float64})
    n += 1 # adjust n from 2-5 to 3-6
    let g = _tfd_g_coeffs
        return 1./((g[n, 1]
                    + g[n, 2] * (ϵ.^(1/2))
                    + g[n, 3] * ϵ
                    + g[n, 4] * (ϵ.^(3/2))
                    + g[n, 5] * (ϵ.^2)).^(n-1))
    end
end

# Evaluating the EOSes

"Get the pressure for an EOS: P = pressure(ρ[, T])"
function pressure(b::BME3, ρ)
    ρ₀ = b.ρ₀
    K₀ = b.K₀
    dK₀ = b.dK₀
    η = ρ/ρ₀
    f = (η^(2/3) - 1)/2
    F = K₀*(1 + (dK₀ - 4)f)

    P = 3f*F*(1 + 2f)^(5/2)
end

function pressure(b::BME4, ρ)
    ρ₀ = b.ρ₀
    K₀ = b.K₀
    dK₀ = b.dK₀
    d2K₀ = b.d2K₀
    η = ρ/ρ₀
    f = (η^(2/3) - 1)/2
    F = K₀*(1 + (dK₀ - 4)f + (K₀*d2K₀ + (dK₀ - 4)*(dK₀ - 3) + 35/9)f^2)

    P = 3f*F*(1 + 2f)^(5/2)
end
pressure(b::BME, ρ, T) = pressure(b, ρ)

function pressure(v::Vinet, ρ)
    let ρ₀ = v.ρ₀, K₀ = v.K₀, dK₀ = v.dK₀, η = ρ/ρ₀
        return (3K₀*η^(2/3) * (1 - η^(-1/3))
                * exp(3/2*(dK₀ - 1)*(1 - η^(-1/3))))
    end
end
pressure(v::Vinet, ρ, T) = pressure(v, ρ)

"Get the base pressure (no thermal component) for an EOS with MGD thermal expansion"
basepressure(mgd::MGDPressureEOS, ρ) = pressure(mgd.eos, ρ)
pressure(mgd::MGDPressureEOS, ρ, T) = basepressure(mgd, ρ) + thermalpressure(mgd, ρ, T)

function pressure(I::IAPWS, ρ, T)
    δ = δ_from_ρ(ρ)
    τ = τ_from_T(T)

    P = ρ * R_h2o * T * (1 + δ*dϕrδ(I, δ, τ))
end

pressure(beos::BoundedEOS, ρ) = pressure(beos.eos, ρ)
pressure(beos::BoundedEOS, ρ, T) = pressure(beos.eos, ρ, T)

function (eos::TFD)(P)
    let Z = eos.Z, A = eos.A, n = eos.n
        # P is in Pa but we want it in dyn/cm^2: 1 Pa = 10 dyn/cm^2
        P *= 10

        # pre-calculations
        ζ   = @. (P / 9.524E13)^(1/5) * Z^(-2/3)
        ϵ   = @. (3 /(32π^2 * Z^2))^(1/3)
        ϕ   = @. (3^(1/3))/20 + ϵ/(4 * (3^(1/3)))
        α   = @. 1/(1.941E-2 - ϵ^(1/2)*6.277E-2 + ϵ*1.076)
        x₀0 = @. (8.884E-3 + (ϵ^(1/2))*4.998E-1
               + ϵ*5.2604E-1)^(-1)
        β₀ = @. x₀0*ϕ - 1
        β₁ = @. β₀*α + ((1 + β₀)/ϕ)
        β₂ = β_(2, ϵ)
        β₃ = β_(3, ϵ)
        β₄ = β_(4, ϵ)
        β₅ = β_(5, ϵ)

        βζ = @. β₀ + β₁*ζ + β₂*ζ^2 + β₃*ζ^3 + β₄*ζ^4 + β₅*ζ^5

        x₀ = @. 1/(ζ + ϕ) * (1 + exp(-α*ζ)*βζ)

        num = sum(n.*A)
        denom = sum(@. n*x₀^3 / Z)
        ρ = num/denom * 3.866

        # rho is in g/cm3 but we want it in kg/m3: 1 g/cm3 = 1000 kg/m3
        return 1000ρ
    end
end

(eos::TFD)(P, T) = eos(P)

function (cg::ChoukrounGrasset)(P, T)
    P = P/1e6  # Pa -> MPa
    density = 1./specificvolume(cg, P, T)
end

(eos::ConstantEOS)(P, T) = eos.ρ
(eos::IdealGas)(P, T) = P / (eos.R * T)

"Minimum density of an EOS"
ρmin(eos::EOS) = eos.ρmin
ρmin(mgd::MGDPressureEOS) = ρmin(mgd.eos)
"Maximum density of an EOS"
ρmax(eos::EOS) = eos.ρmax
ρmax(mgd::MGDPressureEOS) = ρmax(mgd.eos)

function (eos::PolytropicEOS)(P)
    let a = eos.a, ρ₀ = eos.ρ₀, n = eos.n
        return ρ₀ + a*P^n
    end
end

(eos::PolytropicEOS)(P, T) = eos(P)
(beos::BoundedEOS)(P) = beos.eos(P)
(beos::BoundedEOS)(P, T) = beos.eos(P, T)

# Defining the regions within which each function holds

# Bounded EOSes are defined by a bounding box or polygon region
Base.in(x, y, beos::BoundedEOS) = in(x, y, beos.extent)

# Other EOSes have no limit as to where they may be evaluated
Base.in(x, y, eos::EOS) = true

# Temperature dependence information

istempdependent(::BME) = false
istempdependent(::Vinet) = false
istempdependent(::TFD) = false
istempdependent(::PolytropicEOS) = false
istempdependent(::IAPWS) = true
istempdependent(::ConstantEOS) = false
istempdependent(::ChoukrounGrasset) = true
istempdependent(beos::BoundedEOS) = istempdependent(beos.eos)
istempdependent(::MGDPressureEOS) = true

# Save EOS data to file

"Load EOS data from tabular format and save to JLD files as functional representations"
function save_functional_eoses!()
    jldopen("$(config.datadir)/eos-functional.jld", "w") do file
        # bring in phase boundaries
        phaseregions = load_phase_boundaries()["regions"]

        # Choukroun and Grasset's equations of state
        ckg = let
            cgtable = readtable("$(config.rawdata)/ChoukrounGrasset-parameters.dat")
            regionI = phaseregions["I"]
            regionIII = phaseregions["III"]
            regionV = phaseregions["V"]
            regionVI = phaseregions["VI"]

            Dict(
                "I" => BoundedEOS(ChoukrounGrasset(cgtable[:ice_I]...), regionI),
                "III" => BoundedEOS(ChoukrounGrasset(cgtable[:ice_III]...), regionIII),
                "V" => BoundedEOS(ChoukrounGrasset(cgtable[:ice_V]...), regionV),
                "VI" => BoundedEOS(ChoukrounGrasset(cgtable[:ice_VI]...), regionVI))
        end
        write(file, "choukroungrasset", ckg)

        # Polytropic EOS
        pt = Dict(
            "mgsio3" => PolytropicEOS(4100., 0.00161, 0.541),
            "fe" => PolytropicEOS(8300., 0.00349, 0.528),
            "h2o" => PolytropicEOS(1460., 0.00311, 0.513),
            "graphite" => PolytropicEOS(2250., 0.00350, 0.514),
            "sic" => PolytropicEOS(3220., 0.00172, 0.537))
        write(file, "polytropic", pt)

        # Seager's EOS
        sg = Dict(
            "fe_eps" =>   Vinet(8.30e3, 156.2e9, 6.08),
            "h2o" =>        BME(1.46e3, 23.7e9,  4.15),
            "mgsio3_pv" =>  BME4(4.10e3, 247.0e9, 3.97, -0.016e-9, 3e3, 6e4),
            "fe_tfd" =>     TFD(26, 55.845),
            "h2o_tfd" =>    TFD([1, 8], [1.00794, 15.9994], [2., 1.]),
            "mgsio3_tfd" => TFD([12, 14, 8], [24.305, 28.0855, 15.9994],
                                [1., 1., 3.]))
        write(file, "seager", sg)

        # Miscellaneous
        misc = let
            # IAPWS extensions
            extent1 = BoundingBox(1e9, 1e11, 273.16, Tc)
            ρrange1 = (1e-6, 1e6)
            extent2 = BoundingBox(0.05e6, 1e9, 1273., 25000.)
            ρrange2 = (1e-6, 1e3)
            extent3 = BoundingBox(1e9, 1e11, Tc, 25000)
            ρrange3 = (1e-6, 1e5)
            iapws_highpressure = BoundedEOS(IAPWS(ρrange1...), extent1)
            iapws_hightemp = Bange1 = (1e-6, 1e6)
            extent2 = BoundingBox(0.05e6, 1e9, 1273., 25000.)
            ρrange2 = (1e-6, 1e3)
            extent3 = BoundingBox(1e9, 3e10, Tc, 25000)
            ρrange3 = (1e-6, 1e5)
            extent4 = BoundingBox(1e10, 1e12, 10000., 25000.)
            ρrange4 = (1e2, 1e6)
            iapws_highpressure = BoundedEOS(IAPWS(ρrange1...), extent1)
            iapws_hightemp = BoundedEOS(IAPWS(ρrange2...), extent2)
            iapws_highprestemp = BoundedEOS(IAPWS(ρrange3...), extent3)
            iapws_pastfrench = BoundedEOS(IAPWS(ρrange4...), extent4)

            # BME with thermal pressure for ice VII
            ρ₀ = 1464.7   # kg/m3
            K₀ = 23.9e9   # Pa
            T₀ = 300      # K
            dK₀ = 4.2     # dimensionless
            γ₀ = 1.2      # dimensionless
            θD₀ = 1470    # K
            q = -2        # dimensionless
            n = 3         # dimensionless

            bme = BME(ρ₀, K₀, dK₀)
            mgd_bme = MGDPressureEOS(bme, T₀, θD₀, γ₀, q, n)
            bounded_mgd_bme = BoundedEOS(mgd_bme, phaseregions["VII"])

            # ideal gas for low-pressure region
            ideal_gas_eos = IdealGas(R_h2o)
            ideal_gas_region = BoundingBox(0., 1e7, 0., Inf)
            ideal_gas = BoundedEOS(ideal_gas_eos, ideal_gas_region)

            fallback_eos = ConstantEOS(1.)

            Dict(
                "iapws_hightemp" => iapws_hightemp,
                "iapws_highpressure" => iapws_highpressure,
                "iapws_highprestemp" => iapws_highprestemp,
                "iapws_pastfrench" => iapws_pastfrench,
                "mgd_iceVII" => bounded_mgd_bme,
                "ideal_gas" => ideal_gas,
                "fallback" => fallback_eos)
        end
        write(file, "misc", misc)
    end
end


text:
@testset "Manual overrides" begin
    @testset "arf_get_d / arf_get_si" begin
        a = Arf(2.3)

        @test Arblib.get_d(a) isa Float64
        @test Arblib.get_d(a) == 2.3
        @test Arblib.get_d(a, rnd = RoundNearest) == 2.3
        @test Arblib.get_d(a, RoundNearest) == 2.3

        @test Arblib.get_si(a) isa Int
        @test Arblib.get_si(a) == 2
        @test Arblib.get_si(a, rnd = RoundUp) == 3
        @test Arblib.get_si(a, RoundUp) == 3
        @test Arblib.get_si(a, RoundToZero) == 2
    end
end


text:
"""
    axes2size
"""
axes2size(x::AbstractArray) = axes2size(axes(x))
axes2size(x::Tuple) = length.(x)

"""
    axes2length
"""
axes2length(x::AbstractArray) = axes2length(axes(x))
axes2length(x::Tuple{Any,Vararg{Any}}) = length(first(x)) * axes2length(tail(x))
axes2length(x::Tuple{Any}) = length(first(x))

_prod(x::Tuple{Integer,Vararg}) = first(x) * _prod(tail(x))
_prod(x::Tuple{Colon,Vararg}) = _prod(tail(x))
_prod(x::Tuple{Integer}) = first(x)
_prod(x::Tuple{Colon}) = 1


text:
# TODO: Add a more thorough explanation of the dependency confusion
# vulnerability and how this guideline mitigates it.

const guideline_dependency_confusion = Guideline(;
    info="No UUID conflict with other registries.",
    docs=nothing,
    check=data ->
        has_no_dependency_confusion(data.pkg, data.registry_head, data.public_registries),
)

# TODO: Needs a strategy to handle connection failures for the public
# registries. Preferably they should also be cloned only once and then
# just updated to mitigate the effect of them being temporarily
# offline. This could be implemented with the help of the Scratch
# package, but requires Julia >= 1.5.
function has_no_dependency_confusion(pkg, registry_head, public_registries)
    uuid, package_repo = parse_registry_pkg_info(registry_head, pkg)
    for repo in public_registries
        try
            registry = clone_repo(repo)
            registry_toml = TOML.parsefile(joinpath(registry, "Registry.toml"))
            packages = registry_toml["packages"]
            if haskey(packages, uuid)
                message = string(
                    "UUID $uuid conflicts with the package ",
                    packages[uuid]["name"],
                    " in registry ",
                    registry_toml["name"],
                    " at $repo. ",
                    "This could be a dependency confusion attack.",
                )
                # Conflict detected. This is benign if the package name
                # *and* the package URL matches.
                if packages[uuid]["name"] != pkg
                    return false, message
                end
                package_path = packages[uuid]["path"]
                other_package_repo = TOML.parsefile(
                    joinpath(registry, package_path, "Package.toml")
                )["repo"]
                if package_repo != other_package_repo
                    return false, message
                end
            end
        catch
            message = string(
                "Failed to clone public registry $(repo) for a check against dependency confusion.\n",
                "This is an internal issue with the AutoMerge process and has nothing to do with "."the package being registered but requires manual intervention before AutoMerge ",
                "can be resumed.",
            )
            return false, message
        end
    end

    return true, ""
end


text:
using Symbolics, SparseArrays, Test
@variables a b c1 c2 c3 d e g

# Multiple argument matrix
h = [a + b + c1 + c2,
     c3 + d + e + g,
     0] # uses the same number of arguments as our application
h_julia(a, b, c, d, e, g) = [a[1] + b[1] + c[1] + c[2],
                             c[3] + d[1] + e[1] + g[1],
                             0]
function h_julia!(out, a, b, c, d, e, g)
    out .= [a[1] + b[1] + c[1] + c[2], c[3] + d[1] + e[1] + g[1], 0]
end

h_str = Symbolics.build_function(h, [a], [b], [c1, c2, c3], [d], [e], [g])
h_oop = eval(h_str[1])
h_str_par = Symbolics.build_function(h, [a], [b], [c1, c2, c3], [d], [e], [g], parallel=Symbolics.MultithreadedForm())
h_oop_par = eval(h_str_par[1])
h_par_rgf = Symbolics.build_function(h, [a], [b], [c1, c2, c3], [d], [e], [g], parallel=Symbolics.MultithreadedForm(), expression=false)
h_ip! = eval(h_str[2])
h_ip_skip! = eval(Symbolics.build_function(h, [a], [b], [c1, c2, c3], [d], [e], [g], skipzeros=true, fillzeros=false)[2])
h_ip_skip_par! = eval(Symbolics.build_function(h, [a], [b], [c1, c2, c3], [d], [e], [g], skipzeros=true, parallel=Symbolics.MultithreadedForm(), fillzeros=false)[2])
inputs = ([1], [2], [3, 4, 5], [6], [7], [8])

@test h_oop(inputs...) == h_julia(inputs...)
@test h_oop_par(inputs...) == h_julia(inputs...)
@test h_par_rgf[1](inputs...) == h_julia(inputs...)
out_1 = similar(h, Int)
out_2 = similar(out_1)
h_ip!(out_1, inputs...)
h_julia!(out_2, inputs...)
@test out_1 == out_2
out_1 = similar(h, Int)
h_par_rgf[2](out_1, inputs...)
@test out_1 == out_2
fill!(out_1, 10)
h_ip_skip!(out_1, inputs...)
@test out_1[3] == 10
out_1[3] = 0
@test out_1 == out_2

fill!(out_1, 10)
h_ip_skip_par!(out_1, inputs...)
@test out_1[3] == 10
out_1[3] = 0
@test out_1 == out_2

# Multiple input matrix, some unused arguments
h_skip = [a + b + c1; c2 + c3 + g] # skip d, e
h_julia_skip(a, b, c, d, e, g) = [a[1] + b[1] + c[1]; c[2] + c[3] + g[1]]
function h_julia_skip!(out, a, b, c, d, e, g)
    out .= [a[1] + b[1] + c[1]; c[2] + c[3] + g[1]]
end

h_str_skip = Symbolics.build_function(h_skip, [a], [b], [c1, c2, c3], [], [], [g], checkbounds=true)
h_oop_skip = eval(h_str_skip[1])
h_ip!_skip = eval(h_str_skip[2])
inputs_skip = ([1], [2], [3, 4, 5], [], [], [8])

@test h_oop_skip(inputs_skip...) == h_julia_skip(inputs_skip...)
out_1_skip = Array{Int64}(undef, 2)
out_2_skip = similar(out_1_skip)
h_ip!_skip(out_1_skip, inputs_skip...)
h_julia_skip!(out_2_skip, inputs_skip...)
@test out_1_skip == out_2_skip

# Same as above, except test ability to call with non-matrix arguments (i.e., for `nt`)
inputs_skip_2 = ([1], [2], [3, 4, 5], [], (a = 1, b = 2), [8])
@test h_oop_skip(inputs_skip_2...) == h_julia_skip(inputs_skip_2...)
out_1_skip_2 = Array{Int64}(undef, 2)
out_2_skip_2 = similar(out_1_skip_2)
h_ip!_skip(out_1_skip_2, inputs_skip_2...)
h_julia_skip!(out_2_skip_2, inputs_skip_2...)
@test out_1_skip_2 == out_2_skip_2

# Multiple input scalar
h_scalar = a + b + c1 + c2 + c3 + d + e + g
h_julia_scalar(a, b, c, d, e, g) = a[1] + b[1] + c[1] + c[2] + c[3] + d[1] + e[1] + g[1]
h_str_scalar = Symbolics.build_function(h_scalar, [a], [b], [c1, c2, c3], [d], [e], [g])
h_oop_scalar = eval(h_str_scalar)
@test h_oop_scalar(inputs...) == h_julia_scalar(inputs...)

@variables z[1:100]
@test isequal(simplify(Symbolics.unflatten_long_ops(sum(z))),
              simplify(sum(z)))

@test isequal(simplify(Symbolics.unflatten_long_ops(prod(z))),
              simplify(prod(z)))

@variables t x(t) y(t) k
f = eval(build_function((x+y)/k, [x,y,k]))
@test f([1,1,2]) == 1

f = eval(build_function([(x+y)/k], [x,y,k])[1])
@test f([1,1,2]) == [1]

f = eval(build_function([(x+y)/k], [x,y,k])[2])
z = [0.0]
f(z, [1,1,2])
@test z == [1]

f = eval(build_function(sparse([1],[1], [(x+y)/k], 10,10), [x,y,k])[1])

@test size(f([1.,1.,2])) == (10,10)
@test f([1.,1.,2])[1,1] == 1.0
@test sum(f([1.,1.,2])) == 1.0


text:
function VariationalEquation(arg0::ExpandableODE, arg1::ODEJacobiansProvider)
    return VariationalEquation((ExpandableODE, ODEJacobiansProvider), arg0, arg1)
end

function VariationalEquation(arg0::ExpandableODE, arg1::OrdinaryDifferentialEquation, arg2::Vector{jdouble}, arg3::ParametersController, arg4::Vector{ParameterConfiguration})
    return VariationalEquation((ExpandableODE, OrdinaryDifferentialEquation, Vector{jdouble}, ParametersController, Vector{ParameterConfiguration}), arg0, arg1, arg2, arg3, arg4)
end

function extract_main_set_jacobian(obj::VariationalEquation, arg0::ODEState)
    return jcall(obj, "extractMainSetJacobian", Vector{Vector{jdouble}}, (ODEState,), arg0)
end

function extract_parameter_jacobian(obj::VariationalEquation, arg0::ODEState, arg1::JString)
    return jcall(obj, "extractParameterJacobian", Vector{jdouble}, (ODEState, JString), arg0, arg1)
end

function set_initial_main_state_jacobian(obj::VariationalEquation, arg0::Vector{Vector{jdouble}})
    return jcall(obj, "setInitialMainStateJacobian", void, (Vector{Vector{jdouble}},), arg0)
end

function set_initial_parameter_jacobian(obj::VariationalEquation, arg0::JString, arg1::Vector{jdouble})
    return jcall(obj, "setInitialParameterJacobian", void, (JString, Vector{jdouble}), arg0, arg1)
end

function set_up_initial_state(obj::VariationalEquation, arg0::ODEState)
    return jcall(obj, "setUpInitialState", ODEState, (ODEState,), arg0)
end


text:
module Interferometry

using Conv
import JuMIT.Acquisition
import JuMIT.Data

"""
enhance diffractions in the `TD`

# Keyword Arguments

`λdom::Float64=0.0` : distance between receivers must be greater than twice central wavelength, 2*λdom (Shapiro 2005)
`tlag::Float64=data.tgrid[end]-data.tgrid[1]` : maximum lag time in the output traces

"""
function TD_virtual_diff(
			 data::Data.TD;
			 λdom::Float64=0.0,
			 tlag::Float64=data.tgrid[end]-data.tgrid[1]
			)

	nr = maximum(data.acqgeom.nr);	nss = data.acqgeom.nss;	nt = length(data.tgrid)
	fields = data.fields;
	# normalize the records in data
	datan = Data.TD_normalize(data,:recrms)

	# get unique receiver postions; a virtual source at each
	urpos = Acquisition.Geom_get([data.acqgeom],:urpos)
	nur = Acquisition.Geom_get([data.acqgeom],:nur)

	# central wavelength (use zero for testing)
	println(string("dominant wavelength in the data:\t",λdom))

	rx = Array{Vector{Float64}}(nur); rz = Array{Vector{Float64}}(nur);
	datmat = zeros(nur, length(fields), 2*nt-1, nur);
	for ifield =1:length(fields)
		# loop over virtual sources
		for irs = 1:nur

			irvec = [];
			# loop over second receiver
			for ir = 1:nur
				rpos = [urpos[1][irs], urpos[2][irs]];
				rpos0 = [urpos[1][ir], urpos[2][ir]];
				δrpos = sqrt((rpos[1] - rpos0[1])^2 + (rpos[2] - rpos0[2])^2)

				# the distance between receivers must be greater than 2λdom
				# here λdom is the central wavelength (Shapiro 2005)
				if(δrpos > 2.0*λdom)

					# find sources that shoot at these two receivers
					sson = Acquisition.Geom_find(data.acqgeom; rpos=rpos, rpos0=rpos0)
					nsson = count(x->x!=[0],sson);
					if(nsson!=0)
						push!(irvec, ir)
					end
					# stacking over these sources
					for isson=1:length(sson)
						if(sson[isson] != [0])
							datmat[ir, ifield,:,irs] +=
							-1.0 * xcorr(
							  datan.d[isson,ifield][:, sson[isson][2]],
							  datan.d[isson,ifield][:, sson[isson][1]])
						end
					end
					# normalize depending on the stack
					nsson != 0 ? datmat[ ir,  ifield, :, irs] /= nsson : nothing
				end
			end
			if(irvec != [])
				rx[irs] = [urpos[2][i] for i in irvec]
				rz[irs] = [urpos[1][i] for i in irvec]
			end
		end
	end
	# virtual source positions
	sx = [[urpos[2][irs]] for irs in 1:nur];
	sz = [[urpos[1][irs]] for irs in 1:nur];

	# bool array acoording to undef
	ars = [isassigned(rx,irs) ? true : false for irs=1:nur];

	# select only positions that are active
	rx = rx[ars];	rz = rz[ars];
	sx = sx[ars];	sz = sz[ars];

	# number of virtual sources
	nvs = length(rx) == length(sx) ? length(sx) : error("some error")

	# geom
	vacqgeom = Acquisition.Geom(sx, sz, rx, rz, nvs, fill(1,nvs),
			     [length(rx[ir]) for ir=1:length(rx)])

	# tgrid after correlation
	dt = data.tgrid[end]-data.tgrid[1]
	tgridxcorr=range(-dt, stop=+dt, step=step(data.tgrid))

	tlag > 0.0 ? tgridcut=range(-tlag,stop=+tlag,step=step(data.tgrid)) : error("tlag < 0")

	return Data.TD_resamp(
		       Data.TD_urpos(datmat,data.fields,tgridxcorr,vacqgeom,nur,urpos),
		       tgridcut)

end

"""
Correlating noise records with reference records
* `irref` : reference receiver
"""
function TD_noise_corr(data::Data.TD; tlagfrac=0.5, irref=[1],)

	error("fix this")
	# tgrid after correlation
	tt = data.tgrid[end]-data.tgrid[1]
	#tgridc = lag(tlagfrac*tt, step(data.tgrid))

	# allocate TD
	dataout=Data.TD_zeros(data.fields, tgridc[1], data.acqgeom)

	TD_noise_corr!(dataout, data, irref)
end

function TD_noise_corr!(dataout, data, irref)
	nr = data.acqgeom.nr;	nss = data.acqgeom.nss;
	fields=data.fields
	ntlag=Int((length(dataout.tgrid)-1)/2)

	nto = length(dataout.tgrid);
	isodd(nto) || error("odd samples in output data required")
	for ifield = 1:length(fields), iss = 1:nss
		# need param for each supersource, because of different receivers
		paxcorr=Conv.P_xcorr(length(data.tgrid), nr[iss],
		   cglags=[ntlag, ntlag], norm_flag=false,cg_indices=irref)
		dd=data.d[iss, ifield]
		ddo=[dataout.d[iss, ifield]]

		Conv.mod!(paxcorr, cg=ddo, g=dd)

	end

	return dataout

end

end # module


text:
export isempty

function isempty(str::AbstractString, ignore_whitespace::Bool=false)::Bool
    return (ignore_whitespace == true ? length(strip(str)) : length(str)) === 0
end


text:
function initialize!(contact::ContactConstraint{T,N,Nc,Cs}) where {T,N,Nc,Cs}
    initialize_positive_orthant!(contact.impulses[1], contact.impulses_dual[1])
    initialize_positive_orthant!(contact.impulses[2], contact.impulses_dual[2])
    return nothing
end

function initialize!(contact::ContactConstraint{T,N,Nc,Cs}) where {T,N,Nc,Cs<:NonlinearContact{T,N}}
	γort, sort = initialize_positive_orthant!(contact.impulses[1][1:1], contact.impulses_dual[1][1:1])
	γsoc, ssoc = initialize_second_order_cone!(contact.impulses[1][2:4], contact.impulses_dual[1][2:4])
	contact.impulses[1] = [γort; γsoc]
	contact.impulses_dual[1] = [sort; ssoc]
	γort, sort = initialize_positive_orthant!(contact.impulses[2][1:1], contact.impulses_dual[2][1:1])
	γsoc, ssoc = initialize_second_order_cone!(contact.impulses[2][2:4], contact.impulses_dual[2][2:4])
	contact.impulses[2] = [γort; γsoc]
	contact.impulses_dual[2] = [sort; ssoc]
    return nothing
end

function initialize_positive_orthant!(γ::AbstractVector{T}, s::AbstractVector{T};
    ϵ::T = 1e-20) where T

    δs = max(-1.5 * minimum(s), 0)
    δγ = max(-1.5 * minimum(γ), 0)

    sh = s .+ δs
    γh = γ .+ δγ
    δhs = 0.5 * transpose(sh) * γh / (sum(γh) + ϵ)
    δhγ = 0.5 * transpose(sh) * γh / (sum(sh) + ϵ)
    s0 = sh .+ δhs
    γ0 = γh .+ δhγ
	return γ0, s0
end

function initialize_second_order_cone!(γ::AbstractVector{T}, s::AbstractVector{T};
    ϵ::T = 1e-20) where T

    e = [1.0; zeros(length(γ) - 1)] # identity element
    δs = max(-1.5 * (s[1] - norm(s[2:end])), 0)
    δγ = max(-1.5 * (γ[1] - norm(γ[2:end])), 0)

    sh = s + δs * e
    γh = γ + δγ * e
    δhs = 0.5 * transpose(sh) * γh / ((γh[1] + norm(γh[2,end])) + ϵ)
    δhγ = 0.5 * transpose(sh) * γh / ((sh[1] + norm(sh[2,end])) + ϵ)

    s0 = sh + δhs * e
    γ0 = γh + δhγ * e
	return γ0, s0
end


text:
# standard scales
const CHROMATIC      = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
const NATURAL_MINOR  = [0, 2, 3, 5, 7, 8, 10]
const HARMONIC_MINOR = [0, 2, 3, 5, 7, 8, 11]
const MELODIC_MINOR  = [0, 2, 3, 5, 7, 9, 11]
const MAJOR          = [0, 2, 4, 5, 7, 9, 11]
const MINOR          = NATURAL_MINOR

# modes
const IONIAN     = [0, 2, 4, 5, 7, 9, 11] # (I) same as major
const DORIAN     = [0, 2, 3, 5, 7, 9, 10] # (II)
const PHRYGIAN   = [0, 1, 3, 5, 7, 8, 10] # (III)
const LYDIAN     = [0, 2, 4, 6, 7, 9, 9]  # (IV)
const MIXOLYDIAN = [0, 2, 4, 5, 7, 9, 10] # (V)
const AEOLIAN    = [0, 2, 3, 5, 7, 8, 10] # (VI) same as natural minor
const LOCRIAN    = [0, 1, 3, 5, 6, 8, 10] # (VII)

# dictionaries containing all elements
const SCALES = Dict("Chromatic" => CHROMATIC,
                    "Major" => MAJOR,
                    "Minor" => NATURAL_MINOR,
                    "Harmonic Minor" => HARMONIC_MINOR)
const MODES = Dict("Ionian" => IONIAN,
                   "Dorian" => DORIAN,
                   "Phrygian" => PHRYGIAN,
                   "Lydian" => LYDIAN,
                   "Mixolydian" => MIXOLYDIAN,
                   "Aeolian" => AEOLIAN,
                   "Locrian" => LOCRIAN)


text:
using MPI

MPI.Init()

include(joinpath("src", "MPISum.jl"))
include(joinpath("src", "MatrixMaps.jl"))

n = 4
v = ones(n)
A = MPISum(MatrixMap(randn(n, n)))
Av = A(v)

@show MPI.Comm_rank(A.comm), Av

MPI.Finalize()

max_stars_repo_name: mtfishman/ITensorParallel


text:
import FLOWFarm; const ff = FLOWFarm

# set initial turbine x and y locations
diam = 80.0
data = readdlm("inputfiles/layout_38turb_round.txt",  ' ', skipstart=1)
turbine_x = data[:, 1].*diam
nturbines = length(turbine_x)
turbine_y = data[:, 2].*diam
turbine_z = zeros(nturbines)

turbine_x = turbine_x .- turbine_x[1]
turbine_y = turbine_y .- turbine_y[1]

# calculate the number of turbines
nturbines = length(turbine_x)

# set turbine base heights
turbine_z = zeros(nturbines) .+ 0.0

# set turbine yaw values
turbine_yaw = zeros(nturbines)

# set turbine design parameters
rotor_diameter = zeros(nturbines) .+ 80.0 # m
hub_height = zeros(nturbines) .+ 70.0   # m
cut_in_speed = zeros(nturbines) .+4.  # m/s
cut_out_speed = zeros(nturbines) .+25.  # m/s
rated_speed = zeros(nturbines) .+16.  # m/s
rated_power = zeros(nturbines) .+2.0E6  # W
generator_efficiency = zeros(nturbines) .+0.944

# rotor swept area sample points (normalized by rotor radius)
rotor_points_y = [0.0]
rotor_points_z = [0.0]

# set flow parameters
data = readdlm("inputfiles/windrose_nantucket_12dir.txt",  ' ', skipstart=1)
winddirections = data[:, 1].*pi/180.0
windspeeds = data[:,2]
windprobabilities = data[:, 3]
nstates = length(windspeeds)

air_density = 1.1716  # kg/m^3
ambient_ti = 0.077
shearexponent = 0.15
ambient_tis = zeros(nstates) .+ ambient_ti
measurementheight = zeros(nstates) .+ hub_height[1]

# load power curve
powerdata = readdlm("inputfiles/niayifar_vestas_v80_power_curve_observed.txt",  ',', skipstart=1)
velpoints = powerdata[:,1]
powerpoints = powerdata[:,2]*1E6

# initialize power model
power_model = ff.PowerModelPowerPoints(velpoints, powerpoints)
power_models = Vector{typeof(power_model)}(undef, nturbines)
for i = 1:nturbines
    power_models[i] = power_model
end

# load thrust curve
ctdata = readdlm("inputfiles/predicted_ct_vestas_v80_niayifar2016.txt",  ',', skipstart=1)
velpoints = ctdata[:,1]
ctpoints = ctdata[:,2]

# initialize thurst model
ct_model = ff.ThrustModelCtPoints(velpoints, ctpoints)
ct_models = Vector{typeof(ct_model)}(undef, nturbines)
for i = 1:nturbines
    ct_models[i] = ct_model
end

# initialize wind shear model
wind_shear_model = ff.PowerLawWindShear(shearexponent)

# get sorted indecies
sorted_turbine_index = sortperm(turbine_x)

# initialize the wind resource definition
windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheight, air_density, ambient_tis, wind_shear_model)

# set up wake and related models
wakedeficitmodel = ff.GaussYawVariableSpread()

wakedeflectionmodel = ff.GaussYawVariableSpreadDeflection()
wakecombinationmodel = ff.LinearLocalVelocitySuperposition()
localtimodel = ff.LocalTIModelNoLocalTI()

# initialize model set
model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--folder_name"
            help = "The name of the folder containing the images and other data"
            arg_type = String
            default = "labrun5"
        "--previous"
            help = "Compute within a previous results folder."
            arg_type = String
            default = ""
        "--jldfile"
            help = "Which JLD file to load from the previous results folder"
            arg_type = String
            default = ""
        "--failsafe"
            help = "Single process and recursive Bayes tree solving, slower by but more robust with better error stack traces."
            action = :store_true
        "--cxoffset"
            help = "Center x offset"
            arg_type = Float64
            default = 0.0
        "--cyoffset"
            help = "Center y offset"
            arg_type = Float64
            default = 0.0
        "--focalscale"
            help = "Focal distance scaling"
            arg_type = Float64
            default = 1.0
        "--tagscale"
            help = "Tag sized scaling"
            arg_type = Float64
            default = 1.0
        "--batch_size"
            help = "Batch solve every n poses"
            arg_type = Int64
            default = 20
        "--show"
            help = "Show pdfs during computation using `evince`"
            action = :store_true
    end

    return parse_args(s)
end

parsed_args = parse_commandline()

# global folderName, cx, cy, fx, fy, batchSize
#
# for (arg, val) in parsed_args
#     if arg == "folder_name"
#         folderName = val
#     end
#     if arg == "cx"
#         cx = val
#     end
#     if arg == "cy"
#         cy = val
#     end
#     if arg == "fx"
#         fx = val
#     end
#     if arg == "fy"
#         fy = val
#     end
#     if arg == "batch_size"
#         batchSize = val
#     end
# end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 1000,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, 2, false),),
    ((4, 1, false),),
    ((4, 3, false),),
    ((4, 2, true), nothing, (default_tol_relax = 100,)),
    ]
insts["fast"] = [
    ((1000, 1, false),),
    ((200, 1, true), nothing, relaxed_tols),
    ((1000, 2, false),),
    ((200, 2, true), nothing, relaxed_tols),
    ((1000, 3, false), nothing, relaxed_tols),
    ]
insts["various"] = [
    ((1000, 1, false),),
    ((7000, 1, false),),
    ((500, 1, true), nothing, relaxed_tols),
    ((1000, 1, true), nothing, relaxed_tols),
    ((1000, 2, false),),
    ((500, 2, true), nothing, relaxed_tols),
    ((2000, 3, false), nothing, relaxed_tols),
    ((100, 3, true), nothing, relaxed_tols),
    ((400, 3, true), nothing, relaxed_tols),
    ((4000, 3, false), nothing, relaxed_tols),
    ]
return (NonparametricDistrJuMP, insts)


text:
using Printf
using BenchmarkTools

println("-------------------------------")
println(@sprintf "Compute the first four Munchausen numbers")
println("-------------------------------")

"""
    raisedto(x)

Raise x to the power of itself.
"""
function raisedto(x)
    if x == 0
       return 0
    else
       return x^x
    end
end

function find_munchausen_numbers()

   power_of_digits = [raisedto(i) for i in 0:9]

   num = 0
   i = 0
   while num < 4
         thisnumber = i
         # find the sum of the digits raised to themselves
         sumofpowers = 0
         while thisnumber > 0
               digit = thisnumber % 10
               sumofpowers += power_of_digits[digit+1]
               thisnumber ÷= 10
         end
         if i == sumofpowers
            num += 1
            println("Munchausen number: ", num, " ", i)
         end
         i += 1
   end
end

@btime find_munchausen_numbers()

println(" ")


text:
## types ##

const (<:) = issubtype

super(T::DataType) = T.super

## generic comparison ##

==(x,y) = x === y

isequal(x, y) = x == y
isequal(x::FloatingPoint, y::FloatingPoint) = (isnan(x) & isnan(y)) | (signbit(x) == signbit(y)) & (x == y)
isequal(x::Real,          y::FloatingPoint) = (isnan(x) & isnan(y)) | (signbit(x) == signbit(y)) & (x == y)
isequal(x::FloatingPoint, y::Real         ) = (isnan(x) & isnan(y)) | (signbit(x) == signbit(y)) & (x == y)

isless(x::FloatingPoint, y::FloatingPoint) = (!isnan(x) & isnan(y)) | (signbit(x) & !signbit(y)) | (x < y)
isless(x::Real,          y::FloatingPoint) = (!isnan(x) & isnan(y)) | (signbit(x) & !signbit(y)) | (x < y)
isless(x::FloatingPoint, y::Real         ) = (!isnan(x) & isnan(y)) | (signbit(x) & !signbit(y)) | (x < y)

# avoid ambiguity with isequal(::Tuple, ::Tuple)
==(T::(Type...), S::(Type...)) = typeseq(T, S)
==(T::Type, S::Type) = typeseq(T, S)

## comparison fallbacks ##

!=(x,y) = !(x==y)
const ≠ = !=
const ≡ = is
!==(x,y) = !is(x,y)
const ≢ = !==

<(x,y) = isless(x,y)
>(x,y) = y < x
<=(x,y) = !(y < x)
const ≤ = <=
>=(x,y) = (y <= x)
const ≥ = >=
.>(x,y) = y .< x
.>=(x,y) = y .<= x
const .≥ = .>=

# this definition allows Number types to implement < instead of isless,
# which is more idiomatic:
isless(x::Real, y::Real) = x<y
lexcmp(x::Real, y::Real) = isless(x,y) ? -1 : ifelse(isless(y,x), 1, 0)

ifelse(c::Bool, x, y) = Intrinsics.select_value(c, x, y)

cmp(x,y) = isless(x,y) ? -1 : ifelse(isless(y,x), 1, 0)
lexcmp(x,y) = cmp(x,y)
lexless(x,y) = lexcmp(x,y)<0

# cmp returns -1, 0, +1 indicating ordering
cmp(x::Real, y::Real) = int(sign(x-y))

max(x,y) = ifelse(y < x, x, y)
min(x,y) = ifelse(x < y, x, y)

scalarmax(x,y) = max(x,y)
scalarmax(x::AbstractArray, y::AbstractArray) = error("ordering is not well-defined for arrays")
scalarmax(x               , y::AbstractArray) = error("ordering is not well-defined for arrays")
scalarmax(x::AbstractArray, y               ) = error("ordering is not well-defined for arrays")

scalarmin(x,y) = min(x,y)
scalarmin(x::AbstractArray, y::AbstractArray) = error("ordering is not well-defined for arrays")
scalarmin(x               , y::AbstractArray) = error("ordering is not well-defined for arrays")
scalarmin(x::AbstractArray, y               ) = error("ordering is not well-defined for arrays")

## definitions providing basic traits of arithmetic operators ##

+(x::Number) = x
*(x::Number) = x
(&)(x::Integer) = x
(|)(x::Integer) = x
($)(x::Integer) = x

for op = (:+, :*, :&, :|, :$, :min, :max, :kron)
    @eval begin
        # note: these definitions must not cause a dispatch loop when +(a,b) is
        # not defined, and must only try to call 2-argument definitions, so
        # that defining +(a,b) is sufficient for full functionality.
        ($op)(a, b, c)        = ($op)(($op)(a,b),c)
        ($op)(a, b, c, xs...) = ($op)(($op)(($op)(a,b),c), xs...)
        # a further concern is that it's easy for a type like (Int,Int...)
        # to match many definitions, so we need to keep the number of
        # definitions down to avoid losing type information.
    end
end

\(x::Number,y::Number) = y/x

# .<op> defaults to <op>
./(x::Number,y::Number) = x/y
.\(x::Number,y::Number) = y./x
.*(x::Number,y::Number) = x*y
.^(x::Number,y::Number) = x^y
.+(x,y) = x+y
.-(x,y) = x-y

.==(x::Number,y::Number) = x == y
.!=(x::Number,y::Number) = x != y
.< (x::Real,y::Real) = x < y
.<=(x::Real,y::Real) = x <= y
const .≤ = .<=
const .≠ = .!=

# core << >> and >>> takes Int32 as second arg
<<(x,y::Int32)    = no_op_err("<<", typeof(x))
>>(x,y::Int32)    = no_op_err(">>", typeof(x))
>>>(x,y::Int32)   = no_op_err(">>>", typeof(x))
<<(x,y::Integer)  = x << convert(Int32,y)
>>(x,y::Integer)  = x >> convert(Int32,y)
>>>(x,y::Integer) = x >>> convert(Int32,y)

# fallback div and fld implementations
# NOTE: C89 fmod() and x87 FPREM implicitly provide truncating float division,
# so it is used here as the basis of float div().
div{T<:Real}(x::T, y::T) = convert(T,round((x-rem(x,y))/y))
fld{T<:Real}(x::T, y::T) = convert(T,round((x-mod(x,y))/y))
#rem{T<:Real}(x::T, y::T) = convert(T,x-y*trunc(x/y))
#mod{T<:Real}(x::T, y::T) = convert(T,x-y*floor(x/y))

# operator alias
const % = rem
.%(x::Real, y::Real) = x%y
const ÷ = div

# mod returns in [0,y) whereas mod1 returns in (0,y]
mod1{T<:Real}(x::T, y::T) = y-mod(y-x,y)
rem1{T<:Real}(x::T, y::T) = rem(x-1,y)+1
fld1{T<:Real}(x::T, y::T) = fld(x-1,y)+1

# transpose
transpose(x) = x
ctranspose(x) = conj(transpose(x))
conj(x) = x

# transposed multiply
Ac_mul_B (a,b) = ctranspose(a)*b
A_mul_Bc (a,b) = a*ctranspose(b)
Ac_mul_Bc(a,b) = ctranspose(a)*ctranspose(b)
At_mul_B (a,b) = transpose(a)*b
A_mul_Bt (a,b) = a*transpose(b)
At_mul_Bt(a,b) = transpose(a)*transpose(b)

# transposed divide
Ac_rdiv_B (a,b) = ctranspose(a)/b
A_rdiv_Bc (a,b) = a/ctranspose(b)
Ac_rdiv_Bc(a,b) = ctranspose(a)/ctranspose(b)
At_rdiv_B (a,b) = transpose(a)/b
A_rdiv_Bt (a,b) = a/transpose(b)
At_rdiv_Bt(a,b) = transpose(a)/transpose(b)

Ac_ldiv_B (a,b) = ctranspose(a)\b
A_ldiv_Bc (a,b) = a\ctranspose(b)
Ac_ldiv_Bc(a,b) = ctranspose(a)\ctranspose(b)
At_ldiv_B (a,b) = transpose(a)\b
A_ldiv_Bt (a,b) = a\transpose(b)
At_ldiv_Bt(a,b) = transpose(a)\transpose(b)

oftype{T}(::Type{T},c) = convert(T,c)
oftype{T}(x::T,c) = convert(T,c)

widen{T<:Number}(x::T) = convert(widen(T), x)

sizeof(T::Type) = error(string("size of type ",T," unknown"))
sizeof(T::DataType) = if isleaftype(T) T.size else error("type does not have a native size") end
sizeof(::Type{Symbol}) = error("type does not have a native size")
sizeof{T<:Array}(::Type{T}) = error("type $(T) does not have a native size")
sizeof(x) = sizeof(typeof(x))

# copying immutable things
copy(x::Union(Symbol,Number,String,Function,Tuple,LambdaStaticData,
              TopNode,QuoteNode,DataType,UnionType)) = x

# function pipelining
|>(x, f::Callable) = f(x)

# array shape rules

function promote_shape(a::(Int,), b::(Int,))
    if a[1] != b[1]
        error("dimensions must match")
    end
    return a
end

function promote_shape(a::(Int,Int), b::(Int,))
    if a[1] != b[1] || a[2] != 1
        error("dimensions must match")
    end
    return a
end

promote_shape(a::(Int,), b::(Int,Int)) = promote_shape(b, a)

function promote_shape(a::(Int, Int), b::(Int, Int))
    if a[1] != b[1] || a[2] != b[2]
        error("dimensions must match")
    end
    return a
end

function promote_shape(a::Dims, b::Dims)
    if length(a) < length(b)
        return promote_shape(b, a)
    end
    for i=1:length(b)
        if a[i] != b[i]
            error("dimensions must match")
        end
    end
    for i=length(b)+1:length(a)
        if a[i] != 1
            error("dimensions must match")
        end
    end
    return a
end

# shape of array to create for getindex() with indexes I
# drop dimensions indexed with trailing scalars
index_shape(I::Real...) = ()
index_shape(i, I...) = tuple(length(i), index_shape(I...)...)

function throw_setindex_mismatch(X, I)
    if length(I) == 1
        e = DimensionMismatch("tried to assign $(length(X)) elements to $(length(I[1])) destinations")
    else
        e = DimensionMismatch("tried to assign $(dims2string(size(X))) array to $(dims2string(map(length,I))) destination")
    end
    throw(e)
end

# check for valid sizes in A[I...] = X where X <: AbstractArray
# we want to allow dimensions that are equal up to permutation, but only
# for permutations that leave array elements in the same linear order.
# those are the permutations that preserve the order of the non-singleton
# dimensions.
function setindex_shape_check(X::AbstractArray, I...)
    li = ndims(X)
    lj = length(I)
    i = j = 1
    while true
        ii = size(X,i)
        jj = length(I[j])::Int
        if i == li || j == lj
            while i < li
                i += 1
                ii *= size(X,i)
            end
            while j < lj
                j += 1
                jj *= length(I[j])::Int
            end
            if ii != jj
                throw_setindex_mismatch(X, I)
            end
            return
        end
        if ii == jj
            i += 1
            j += 1
        elseif ii == 1
            i += 1
        elseif jj == 1
            j += 1
        else
            throw_setindex_mismatch(X, I)
        end
    end
end

setindex_shape_check(X::AbstractArray) =
    (length(X)==1 || throw_setindex_mismatch(X,()))

setindex_shape_check(X::AbstractArray, i) =
    (length(X)==length(i) || throw_setindex_mismatch(X, (i,)))

setindex_shape_check{T}(X::AbstractArray{T,1}, i) =
    (length(X)==length(i) || throw_setindex_mismatch(X, (i,)))

setindex_shape_check{T}(X::AbstractArray{T,1}, i, j) =
    (length(X)==length(i)*length(j) || throw_setindex_mismatch(X, (i,j)))

function setindex_shape_check{T}(X::AbstractArray{T,2}, i, j)
    li, lj = length(i), length(j)
    if length(X) != li*lj
        throw_setindex_mismatch(X, (i,j))
    end
    sx1 = size(X,1)
    if !(li == 1 || li == sx1 || sx1 == 1)
        throw_setindex_mismatch(X, (i,j))
    end
end

# convert to integer index
to_index(i::Int) = i
to_index(i::Real) = convert(Int,i)::Int
to_index(r::UnitRange{Int}) = r
to_index(r::Range{Int}) = r
to_index(I::UnitRange{Bool}) = find(I)
to_index(I::Range{Bool}) = find(I)
to_index{T<:Real}(r::UnitRange{T}) = to_index(first(r)):to_index(last(r))
to_index{T<:Real}(r::StepRange{T}) = to_index(first(r)):to_index(step(r)):to_index(last(r))
to_index(I::AbstractArray{Bool}) = find(I)
to_index(A::AbstractArray{Int}) = A
to_index{T<:Real}(A::AbstractArray{T}) = [to_index(x) for x in A]
to_index(i1, i2)         = to_index(i1), to_index(i2)
to_index(i1, i2, i3)     = to_index(i1), to_index(i2), to_index(i3)
to_index(i1, i2, i3, i4) = to_index(i1), to_index(i2), to_index(i3), to_index(i4)
to_index(I...) = to_index(I)
to_index(I::(Any,))            = (to_index(I[1]), )
to_index(I::(Any,Any,))        = (to_index(I[1]), to_index(I[2]))
to_index(I::(Any,Any,Any))     = (to_index(I[1]), to_index(I[2]), to_index(I[3]))
to_index(I::(Any,Any,Any,Any)) = (to_index(I[1]), to_index(I[2]), to_index(I[3]), to_index(I[4]))
to_index(I::Tuple) = map(to_index, I)
to_index(i) = error("invalid index: $i")

# Addition/subtraction of ranges
for f in (:+, :-)
    @eval begin
        function $f(r1::OrdinalRange, r2::OrdinalRange)
            r1l = length(r1)
            r1l == length(r2) || error("argument dimensions must match")
            range($f(r1.start,r2.start), $f(step(r1),step(r2)), r1l)
        end

        function $f{T<:FloatingPoint}(r1::FloatRange{T}, r2::FloatRange{T})
            len = r1.len
            len == r2.len || error("argument dimensions must match")
            divisor1, divisor2 = r1.divisor, r2.divisor
            if divisor1 == divisor2
                FloatRange{T}($f(r1.start,r2.start), $f(r1.step,r2.step),
                              len, divisor1)
            else
                d1 = int(divisor1)
                d2 = int(divisor2)
                d = lcm(d1,d2)
                s1 = div(d,d1)
                s2 = div(d,d2)
                FloatRange{T}($f(r1.start*s1, r2.start*s2),
                              $f(r1.step*s1, r2.step*s2),  len, d)
            end
        end

        $f(r1::FloatRange, r2::FloatRange) = $f(promote(r1,r2)...)
        $f(r1::FloatRange, r2::OrdinalRange) = $f(promote(r1,r2)...)
        $f(r1::OrdinalRange, r2::FloatRange) = $f(promote(r1,r2)...)
    end
end

# vectorization

macro vectorize_1arg(S,f)
    S = esc(S); f = esc(f); T = esc(:T)
    quote
        ($f){$T<:$S}(x::AbstractArray{$T,1}) = [ ($f)(x[i]) for i=1:length(x) ]
        ($f){$T<:$S}(x::AbstractArray{$T,2}) =
            [ ($f)(x[i,j]) for i=1:size(x,1), j=1:size(x,2) ]
        ($f){$T<:$S}(x::AbstractArray{$T}) =
            reshape([ ($f)(x[i]) for i=1:length(x) ], size(x))
    end
end

macro vectorize_2arg(S,f)
    S = esc(S); f = esc(f); T1 = esc(:T1); T2 = esc(:T2)
    quote
        ($f){$T1<:$S, $T2<:$S}(x::($T1), y::AbstractArray{$T2}) =
            reshape([ ($f)(x, y[i]) for i=1:length(y) ], size(y))
        ($f){$T1<:$S, $T2<:$S}(x::AbstractArray{$T1}, y::($T2)) =
            reshape([ ($f)(x[i], y) for i=1:length(x) ], size(x))

        function ($f){$T1<:$S, $T2<:$S}(x::AbstractArray{$T1}, y::AbstractArray{$T2})
            shp = promote_shape(size(x),size(y))
            reshape([ ($f)(x[i], y[i]) for i=1:length(x) ], shp)
        end
    end
end

# vectorized ifelse

function ifelse(c::AbstractArray{Bool}, x, y)
    reshape([ifelse(ci, x, y) for ci in c], size(c))
end

function ifelse(c::AbstractArray{Bool}, x::AbstractArray, y::AbstractArray)
    shp = promote_shape(size(c), promote_shape(size(x), size(y)))
    reshape([ifelse(c[i], x[i], y[i]) for i = 1 : length(c)], shp)
end

function ifelse(c::AbstractArray{Bool}, x::AbstractArray, y)
    shp = promote_shape(size(c), size(c))
    reshape([ifelse(c[i], x[i], y) for i = 1 : length(c)], shp)
end

function ifelse(c::AbstractArray{Bool}, x, y::AbstractArray)
    shp = promote_shape(size(c), size(y))
    reshape([ifelse(c[i], x, y[i]) for i = 1 : length(c)], shp)
end

# some operators not defined yet
global //, .>>, .<<, >:, <|, |>, hcat, hvcat, ⋅, ×, ∈, ∉, ∋, ∌, ⊆, ⊈, ⊊, ∩, ∪, √, ∛

module Operators

export
    !,
    !=,
    !==,
    $,
    %,
    .%,
    &,
    *,
    +,
    -,
    .!=,
    .+,
    .-,
    .*,
    ./,
    .<,
    .<=,
    .==,
    .>,
    .>=,
    .\,
    .^,
    /,
    //,
    <,
    <:,
    >:,
    <<,
    <=,
    ==,
    >,
    >=,
    ≥,
    ≤,
    ≠,
    .≥,
    .≤,
    .≠,
    >>,
    .>>,
    .<<,
    >>>,
    \,
    ^,
    |,
    |>,
    <|,
    ~,
    ÷,
    ⋅,
    ×,
    ∈,
    ∉,
    ∋,
    ∌,
    ⊆,
    ⊈,
    ⊊,
    ∩,
    ∪,
    √,
    ∛,
    colon,
    hcat,
    vcat,
    hvcat,
    getindex,
    setindex!,
    transpose,
    ctranspose

import Base: !, !=, $, %, .%, &, *, +, -, .!=, .+, .-, .*, ./, .<, .<=, .==, .>,
    .>=, .\, .^, /, //, <, <:, <<, <=, ==, >, >=, >>, .>>, .<<, >>>,
    <|, |>, \, ^, |, ~, !==, >:, colon, hcat, vcat, hvcat, getindex, setindex!,
    transpose, ctranspose,
    ≥, ≤, ≠, .≥, .≤, .≠, ÷, ⋅, ×, ∈, ∉, ∋, ∌, ⊆, ⊈, ⊊, ∩, ∪, √, ∛

end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
#=
The area of a circle is defined as πr^2. Estimate π to 3 decimal places using a Monte Carlo method.

Hint: The basic equation of a circle is x2 + y2 = r2
=#
function estimate_pi(;number_of_decimal::Int64 = 3)
    in = 0
    out = 0
    pi_estimated = 0.0
    n_points = 0
    while round(pi_estimated; digits = number_of_decimal) != round(pi; digits = number_of_decimal)
        x, y = rand(2)
        if (x^2 + y^2) <= 1
            in = in + 1
        end
        n_points = n_points + 1
        pi_estimated = in/n_points * 4
    end
    return pi_estimated
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using Distributed
using Test

@assert nprocs() == 1
addprocs(3, exeflags = "--project=$(Base.active_project())")

@everywhere begin
    using Distributed
    using SharedArrays
    using ClimaComms
    using ClimaCommsSA
    include(joinpath(@__DIR__, "..", "..", "..", "test", "stencil.jl"))
end

@testset "Stencil" begin
    @everywhere stencil_test(ClimaCommsSA.SACommsContext)
end