text:
#
# Functions to find the appropriate ARMA model from data.
#

using FFTW
using Statistics
include("RandomMatrix.jl")
include("optimize_exponentials.jl")

function padded_length(N::Integer)
    powerof2 = round(Int, 2^ceil(log2(N)))
    if N == powerof2
        return N
    elseif N >= 0.75 * powerof2
        return powerof2
    elseif N >= 0.625 * powerof2
        return div(3*powerof2, 4)
    else
        return div(5*powerof2, 8)
    end
end

"""
    estimate_covariance(timeseries::Vector, nsamp::Int, chunklength::Int)

Take a noise sequence (as a 1D vector or as columns of a 2D matrix) and
estimate the sample-sample covariance vector from these.

`nsamp` and `chunklength` are both optional arguments.

`nsamp` is the number of samples to compute from the covariance. If not a positive
number, then compute all possible values.

`chunklength` is the amount of data to be considered at once. If omitted, then this
will be either the length of the timeseries or 15*nsamp, whichever is shorter.

The idea is that a full computation is very inefficient. Making chunklength
much smaller than the timeseries length allows computation more efficient, but
at the cost of losing a small amount of correlation between values in successive
chunks. The recommended value for long data sets is to let chunklength be at
least 10-20x nsamp.
"""
estimate_covariance(timeseries::AbstractVector) = estimate_covariance(timeseries, length(timeseries))

function estimate_covariance(timeseries::AbstractVector, nsamp::Int)
    N = length(timeseries)
    ideal_chunksize = 15*nsamp
    nchunks = div(N-1+ideal_chunksize, ideal_chunksize)
    chunklength = div(N, nchunks)
    estimate_covariance(timeseries, nsamp, chunklength)
end

function estimate_covariance(timeseries::Vector, nsamp::Int, chunklength::Int)
    N = length(timeseries)
    if nsamp > N
        error("Cannot compute $(nsamp) covariance values from a length-$N data.")
    end
    paddedsize = padded_length(chunklength+nsamp)
    padded_data = fill(0.0, paddedsize)
    result = fill(0.0, nsamp)

    i=0
    chunks_consumed = 0
    while i+chunklength <= N
        datamean = mean(timeseries[i+1:i+chunklength])
        padded_data[1:chunklength] = timeseries[i+1:i+chunklength] .- datamean
        chunks_consumed += 1

        power = abs2.(rfft(padded_data))
        acsum = irfft(power, length(padded_data))
        result += acsum[1:nsamp]
        i += chunklength
    end
    result ./ (chunks_consumed * collect(chunklength:-1:chunklength+1-nsamp))
end

"""
    main_exponentials(data, nexp; minexp=nothing)

Find the `nexp` "main exponentials" in the time stream `data`.

Specifically, follow the prescription of P. De Groen & B. De Moor (1987).
"The fit of a sum of exponentials to noisy data." J. Computational & Applied
Mathematics, vol. 20, pages 175–187.

We build a Hankel matrix H whose columns are contiguous segments of the time stream.
Perform a singular value decomposition to find the rank-`nexp` decomposition of H,
from which we can construct the square matrix A (of size `nexp`) which is known
to be similar to the "system matrix" (i.e., a time-step-advance matrix). The
latter has eigenvalues equal to the exponentials comprising the time stream, so
A (by similarity) has the same eigenvalues.

That's the best explanation I have. It might be more honest to say that we
find the exponentials by magic.

Note that the current implementation does not construct the exact SVD of H,
but rather uses a randomized matrix technique to compute an *approximate* SVD
very efficiently, containing only a specified number of leading singular vectors.

If `minexp==nothing` (the default), then only the length-`nexp` model is found.
If `minexp ≤ nexp`, then this function returns an array of arrays, one for each
value `p` with `minexp ≤ p ≤ nexp`, so that the length of the first is `minexp`
and the length of the last is `nexp`.

TODO: we need a better heuristic for constructing H out of a sufficient number
of segments from the data, without allowing it to get insanely large when `data`
is very, very long.
"""
function main_exponentials(data::Vector, nexp::Int; minexp=nothing)
    N = length(data)
    if 2nexp > N
        error("Cannot compute $(nexp) exponentials from data with fewer than twice as many data values.")
    end
    if (minexp != nothing) && (minexp > nexp)
        error("The value of `minexp` must be nothing or a number no larger than `nexp`.")
    end

    ncol = min(40 + 5nexp, div(N, 2))
    H = fill(0.0, N+1-ncol, ncol)
    for c=1:ncol
        H[:,c] = data[c:c+N-ncol]
    end
    U,s,V = find_svd_randomly(H[1:end-1,:], nexp)
    W = Diagonal(s .^ (-0.5))
    A = W*U'*H[2:end,:]*V*W

    # By default, return the exponential set of size nexp only.
    if minexp == nothing
        return sortbases!(eigvals(A))
    end

    # The minexp has been set, so return all sizes minexp...nexp, inclusive.
    exps = []
    for p = minexp:nexp
        A = W[1:p, 1:p]*U[:,1:p]'*H[2:end,:]*V[:,1:p]*W[1:p,1:p]
        push!(exps, sortbases!(eigvals(A)))
    end
    exps
end

"""
    sortbases!(b)

Re-order the elements of `b` in place, such that for every pair b[1],b[2] and so
on, either both are real, or they are complex conjugates of one another.
Elements of `b` with negative imaginary parts are simply assumed to be
conjugates of those with positive imaginary parts. The complex pairs will come
first, followed by the real elements.

Returns the sorted `b`.
"""
function sortbases!(b::Vector)
    real_ones = real(b[imag(b) .== 0])
    sort!(real_ones, rev=true)
    if length(real_ones) < length(b)
        imag_ones = b[imag(b) .> 0]
        if length(imag_ones) != sum(imag(b) .< 0)
            error("sortbases!(b) requires that b contain equal numbers of elements with imag(x) > 0 and < 0.")
        end
        for i=1:length(imag_ones)
            b[2i-1] = imag_ones[i]
            b[2i] = conj(imag_ones[i])
        end
        b[2*length(imag_ones)+1:end] = real_ones[:]
    else
        b[:] = real_ones[:]
    end
    b
end

"""
    B2C(B)

Converts an array `B` of exponential bases to the coefficients of the set of
monic quadratic polynomials.

It is assumed that B are arranged such that B[1] and B[2] are either both real
or are complex conjugates of one another, and similarly for elements (3,4), (5,6),
and so on. Use `sortbases!(B)` to achieve this.

If `B` violates these assumptions, then the computed `C` will make no sense.

Returns `C` such that B[1:2] are the roots of x^2+C[1]*X+C[2], and similarly for
[3:4] and all other pairs. If length(B) is odd, then B[end] is the root of x+C[end],
thus `-1 = B[end]*C[end]`.
"""
function B2C(B::AbstractVector{T}) where {T<:Number}
    C = fill(0.0, length(B))
    n = length(B)
    for i=1:2:n-1
        C[i] = -real(B[i]+B[i+1])
        C[i+1] = real(B[i]*B[i+1])
    end
    if n%2 > 0
        C[end] = -1/real(B[end])
    end
    C
end

"""
    C2B(C)

Converts a real array of polynomial coefficients `C` to their roots. See `B2C(B)`
for more.

Returns `B`, the possibly complex roots.
"""
function C2B(C::AbstractVector{T}) where {T<:Real}
    B = fill(zero(Complex{eltype(C)}), size(C))
    n = length(C)
    iscomplex = false
    for i = 1:2:n-1
        x = -0.5C[i]
        disc = x*x-C[i+1]
        if disc >= 0
            y::ComplexF64 = sqrt(disc)
        else
            y = 1im*sqrt(-disc)
            iscomplex = true
        end
        B[i] = x+y
        B[i+1] = x-y
    end
    if n%2 > 0
        B[end] = -1/C[end]
    end
    if !iscomplex
        return real(B)
    end
    B
end

"""
    findA(t, r, B, [w=weights])

Compute the amplitudes `A` for a linear model of the form

r_t ≈ f(t) = Sum_{i=1...p} A[i] B[i]^t

by minimizing

Sum_t [w_t (r_t - f(t))^2].

This is a linear problem, given `B`. If `w` is not given, then equal weights on
all data are assumed.

Returns the array of amplitudes `A`.
"""
function findA(t::AbstractVector, r::Vector, B::Vector{T}; w=nothing) where {T<:Number}
    @assert length(t) == length(r)
    if w==nothing
        w = ones(r)
    end
    wr = w.*r
    p = length(B)
    M = fill(zero(T), p, p)
    D = fill(zero(T), p)
    for i=1:p
        for j=1:i
            M[j,i] = M[i,j] = sum(w .* (B[i]*B[j]).^t)
        end
        D[i] = sum(wr .* B[i].^t)
    end

    # A = M\D is theoretically correct, but don't do this: can be ill-conditioned in certain cases
    A = pinv(M) * D
end

"""
    exponential_model(t, A, B)`

Computes and returns the sum-of-exponentials model with amplitudes `A` and
exponential bases `B` at time steps `t`.
"""
function exponential_model(t::AbstractVector, A::Vector, B::Vector)
    r = fill(0.0, length(t))
    for i=1:length(A)
        r .+= real(A[i]* B[i].^t)
    end
    r
end

"""
    squashgrowingexp!(bases::Vector, epsfrom1)

Replace any growing exponentials, that is, where  `abs(bases[j]) >= 1`, with
random decaying exponentials. These will be between `1-epsfrom1` and 1 in absolute
value (suggest `epsfrom1 <= 0.01`).
"""
function squashgrowingexp!(bases::Vector, epsfrom1::Float64)
    lastscale = 1 - 0.5*epsfrom1
    for (i,b) = enumerate(bases)
        if abs(b) >= 1 - 0.01*epsfrom1
            if i%2==1 || imag(b) == 0.0
                lastscale = (1-rand(1)[1]*epsfrom1) / abs(b)
            end
            bases[i] *= lastscale
        end
    end
end

"""
    fit_exponentials(data; pmin=0, pmax=6, w=nothing, deltar=nothing, good_enough=0.0)

Fit the time series `data` (regularly sampled) as a sum of possibly complex
exponentials numbering at least `pmin` and at most `pmax`. Uses
`main_exponentials` to determine a starting guess for the exponential bases.
Then a nonlinear fit is attempted, using `NLopt` with the `:LD_MMA` method.

Returns `A,B`, the amplitudes and exponential bases of the best-fit model. The
length of either gives the order of the best-fit model.

# Optional Arguments:
The penalty function is a weighted sum-of-square-errors. If the optional
argument `w` is given, it must be a vector of the same length as `data`, and
then the penalty is the sum over all samples of `w[i]*(model[i]-data[i])^2`.

If `w` is absent, then the weights are taken to be all equal to `deltar^(-2)`.

If `deltar` is also not given, then it is estimated from the standard deviation
of the successive-differences from the last 1/4 of the samples in `data`, on the
assumption that these are noise-dominated. If they aren't, then the caller must
offer a better figure for `w` or `deltar`.

If `pmin` < `pmax`, then models of more than one order will be tried, and the one
with the lowest penalty function will be returned.

If `good_enough` is assigned a positive value, and if any model order `p<pmax`
produces a penalty function less than this value, then this function returns
early with the lowest-order model that meets the criterion.

See also: [`fitARMA`](@ref)
"""
function fit_exponentials(data::Vector; pmin::Int=0, pmax::Int=6,
    w=nothing, deltar=nothing, good_enough::Float64=0.0)

    guess_exponentials = main_exponentials(data, pmax, minexp=pmin)

    N = length(data)
    if w == nothing
        # If deltar isn't given, use the heuristic that the std dev of the
        # diff of the last values in data estimates sqrt(2) times deltar.
        if deltar == nothing
            n = div(N, 4)
            deltar = std(diff(data[end-n:end]), mean=0) / sqrt(2)
        end

        w = ones(Float64, N) / float(deltar^2)
    end

    best_cost = +Inf
    best_fit = nothing
    for (guessB, p) = zip(guess_exponentials, pmin:pmax)
        squashgrowingexp!(guessB, 0.25/N)
        guessC = B2C(guessB)
        cost, fitC = optimize_exponentials(data, w, guessC)
        if (cost < best_cost)
            best_cost = cost
            best_fit = fitC
            # Stop early if cost is small enough.
            if cost <= good_enough
                break
            end
        end
    end

    bestB = C2B(best_fit)
    squashgrowingexp!(bestB, 0.25/N)

    t = 0:(length(data)-1)
    bestA = findA(t, data, bestB, w=w)
    bestA, bestB
end

"""
    fitARMA(covariance, p, q; kwargs)

Fit an `ARMAModel` of order `(p,q)` or lower to the `covariance` data. The model
order with the lowest cost function will be returned.

# Optional Arguments
- `pmin`: the minimum allowed ARMA order.
- `w`: the vector of weights to apply to the data (default: nothing, meaning equal weights)
- `deltar`: rms uncertainty on the covariance values.
- `good_enough`: if given, then the lowest-order model with a cost function at least
   this small will be used, rather than trying all orders up to `p`.

See also: [`fit_exponentials`](@ref)
"""
function fitARMA(covariance::Vector, p::Int, q::Int;
    w=nothing, deltar=nothing, good_enough::Float64=0.0, pmin::Int=0)
    nspecial = max(1+q-p, 0)
    amplitudes, bases = fit_exponentials(covariance[1+nspecial:end], pmin=pmin, pmax=p,
        w=w, deltar=deltar, good_enough=good_enough)
    if nspecial > 0
        amplitudes ./= (bases .^ nspecial)
    end
    ARMAModel(bases, amplitudes, covariance[1:nspecial])
end

# Absent further information, we find that q=p is a good choice most of the time,
# so let that be a default. Also, if we need a maximum order, we can supply
# the rather arbitrary guess of 6.
fitARMA(covariance::Vector, p::Int; kwargs...) = fitARMA(covariance, p, p; kwargs...)
fitARMA(covariance::Vector; kwargs...) = fitARMA(covariance, 6; kwargs...)


text:
# Basic canopy architecture parameters, 10/10/00 S.Kim
# modified to represent heterogeneous canopies
# Uniform continuouse canopy: Width2 = 0
# Hedgerow canopy : Width1 = row width, Width2 = interrow, Height1 = hedge height, Height2 = 0
# Intercropping canopy: Height1, Width1, LA1 for Crop1, and so on
# Rose bent canopy: Height1=Upright canopy, Height2 = bent portion height, 10/16/02 S.Kim

# absorptance not explicitly considered here because it's a leaf characteristic not canopy
# Scattering is considered as canopy characteristics and reflectance is computed based on canopy scattering
# 08-20-12, SK

@enum LeafAngle begin
    spherical = 1
    horizontal = 2
    vertical = 3
    diaheliotropic = 4
    empirical = 5
    ellipsoidal = 6
end

#HACK not used
@enum Cover begin
    glass = 1
    acrylic = 2
    polyethyl = 3
    doublepoly = 4
    whitewashed = 5
    no_cover = 6
end

@enum WaveBand begin
    photosynthetically_active_radiation = 1
    near_infrared = 2
    longwave = 3
end

@system Radiation begin
    sun ~ ::Sun(override)
    LAI: leaf_area_index ~ track(override)

    leaf_angle => ellipsoidal ~ preserve::LeafAngle(parameter)

    # ratio of horizontal to vertical axis of an ellipsoid
    LAF: leaf_angle_factor => begin
        #1
        # leaf angle factor for corn leaves, Campbell and Norman (1998)
        #1.37
        # leaf angle factor for Radish canopy, from Rizzalli et al. (2002),  X factor in Campbell and Norman (1998)
        0.7
    end ~ preserve(parameter)

    wave_band => photosynthetically_active_radiation ~ preserve::WaveBand(parameter)

    # scattering coefficient (reflectance + transmittance)
    s: scattering => 0.15 ~ preserve(parameter)

    # clumping index
    clumping => 1.0 ~ preserve(parameter)

    #FIXME reflectance?
    #r_h => 0.05 ~ preserve(parameter)

    # Forward from Sun

    #TODO better name to make it drive?
    current_zenith_angle(sun.zenith_angle) ~ track(u"°")
    elevation_angle(sun.αs) ~ track(u"°")
    I0_dr(sun.PARdir): directional_photosynthetic_radiation ~ track(u"μmol/m^2/s" #= Quanta =#)
    I0_df(sun.PARdif): diffusive_photosynthetic_radiation ~ track(u"μmol/m^2/s" #= Quanta =#)

    # Leaf angle stuff?

    #TODO better name?
    leaf_angle_coeff(a=leaf_angle, leaf_angle_factor; zenith_angle(u"°")) => begin
        elevation_angle = 90u"°" - zenith_angle
        #FIXME need to prevent zero like sin_beta / cot_beta?
        α = elevation_angle
        t = zenith_angle
        # leaf angle distribution parameter
        x = leaf_angle_factor
        if a == spherical
            # When Lt accounts for total path length, division by sin(elev) isn't necessary
            1 / (2sin(α))
        elseif a == horizontal
            1.
        elseif a == vertical
            1 / (tan(α) * π/2)
        elseif a == empirical
            0.667
        elseif a == diaheliotropic
            1 / sin(α)
        elseif a == ellipsoidal
            sqrt(x^2 + tan(t)^2) / (x + 1.774 * (x+1.182)^-0.733)
        else
            1.
        end
    end ~ call

    #TODO make it @property if arg is not needed
    # Kb: Campbell, p 253, Ratio of projected area to hemi-surface area for an ellipsoid
    #TODO rename to extinction_coeff?
    # extinction coefficient assuming spherical leaf dist
    Kb_at(leaf_angle_coeff, clumping; zenith_angle(u"°")): projection_ratio_at => begin
        leaf_angle_coeff(zenith_angle) * clumping
    end ~ call

    Kb(Kb_at, current_zenith_angle): projection_ratio => begin
        Kb_at(current_zenith_angle)
    end ~ track

    # diffused light ratio to ambient, integrated over all incident angles from 0 to 90
    Kd_F(leaf_angle_coeff, LAI; a): diffused_fraction_for_Kd => begin
        c = leaf_angle_coeff(a)
        x = exp(-c * LAI)
        # Why multiplied by 2?
        2x * sin(a) * cos(a)
    end ~ integrate(from=0, to=π/2)

    # Kd: K for diffuse light, the same literature as above
    Kd(F=Kd_F, LAI, clumping): diffusion_ratio => begin
        K = -log(F) / LAI
        K * clumping
    end ~ track

    ###############################
    # de Pury and Farquhar (1997) #
    ###############################

    # Kb1: Kb prime in de Pury and Farquhar(1997)
    #TODO better name
    Kb1(Kb, s): projection_ratio_prime => (Kb * sqrt(1 - s)) ~ track

    # Kd1: Kd prime in de Pury and Farquhar(1997)
    #TODO better name
    Kd1(Kd, s): diffusion_ratio_prime => (Kd * sqrt(1 - s)) ~ track

    ################
    # Reflectivity #
    ################

    #TODO better name?
    reflectivity(rho_h, Kb, Kd) => begin
        rho_h * (2Kb / (Kb + Kd))
    end ~ track

    # canopy reflection coefficients for beam horizontal leaves, beam uniform leaves, and diffuse radiations

    # rho_h: canopy reflectance of beam irradiance on horizontal leaves, de Pury and Farquhar (1997)
    # also see Campbell and Norman (1998) p 255 for further info on potential problems
    rho_h(s): canopy_reflectivity_horizontal_leaf => begin
        (1 - sqrt(1 - s)) / (1 + sqrt(1 - s))
    end ~ track

    #TODO make consistent interface with siblings
    # rho_cb: canopy reflectance of beam irradiance for uniform leaf angle distribution, de Pury and Farquhar (1997)
    rho_cb_at(rho_h, Kb_at; zenith_angle): canopy_reflectivity_uniform_leaf_at => begin
        Kb = Kb_at(zenith_angle)
        1 - exp(-2rho_h * Kb / (1 + Kb))
    end ~ call

    rho_cb(rho_cb_at, current_zenith_angle): canopy_reflectivity_uniform_leaf => begin
        rho_cb_at(current_zenith_angle)
    end ~ track

    rho_cd_F(rho_cb_at; a): diffused_fraction_for_rho_cd => begin
        x = rho_cb_at(a)
        # Why multiplied by 2?
        2x * sin(a) * cos(a)
    end ~ integrate(from=0, to=π/2)

    # rho_cd: canopy reflectance of diffuse irradiance, de Pury and Farquhar (1997) Table A2
    rho_cd(I0_df, rho_cd_F): canopy_reflectivity_diffusion => begin
        # Probably the eqn A21 in de Pury is missing the integration terms of the angles??
        iszero(I0_df) ? 0. : rho_cd_F
    end ~ track

    # rho_soil: soil reflectivity for PAR band
    rho_soil: soil_reflectivity => 0.10 ~ preserve(parameter)

    #######################
    # I_l?: dePury (1997) #
    #######################

    # I_lb: dePury (1997) eqn A3
    I_lb(I0_dr, rho_cb, Kb1; L): irradiance_lb => begin
        I0_dr * (1 - rho_cb) * Kb1 * exp(-Kb1 * L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # I_ld: dePury (1997) eqn A5
    I_ld(I0_df, rho_cb, Kd1; L): irradiance_ld => begin
        I0_df * (1 - rho_cb) * Kd1 * exp(-Kd1 * L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # I_l: dePury (1997) eqn A5
    I_l(; L): irradiance_l => (I_lb(L) + I_ld(L)) ~ call(u"μmol/m^2/s" #= Quanta =#)

    # I_lbSun: dePury (1997) eqn A5
    I_lbSun(I0_dr, s, Kb, I_lSh; L): irradiance_l_sunlit => begin
        I_lb_sunlit = I0_dr * (1 - s) * Kb
        #TODO: check name I_lbSun vs. I_l_sunlit?
        I_l_sunlit = I_lb_sunlit + I_lSh(L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # I_lSh: dePury (1997) eqn A5
    I_lSh(I_ld, I_lbs; L): irradiance_l_shaded => begin
        I_ld(L) + I_lbs(L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # I_lbs: dePury (1997) eqn A5
    #FIXME: check name I_lbs vs. I_lbSun
    I_lbs(I0_dr, rho_cb, s, Kb1, Kb; L): irradiance_lbs => begin
        I0_dr * ((1 - rho_cb) * Kb1 * exp(-Kb1 * L) - (1 - s) * Kb * exp(-Kb * L))
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # I0tot: total irradiance at the top of the canopy,
    # passed over from either observed PAR or TSolar or TIrradiance
    I0_tot(I0_dr, I0_df): irradiance_I0_tot => (I0_dr + I0_df) ~ track(u"μmol/m^2/s" #= Quanta =#)

    ########
    # I_c? #
    ########

    # I_tot, I_sun, I_shade: absorbed irradiance integrated over LAI per ground area
    #FIXME: not used, but seems producing very low values, need to check equations

    # I_c: Total irradiance absorbed by the canopy, de Pury and Farquhar (1997)
    I_c(rho_cb, I0_dr, I0_df, Kb1, Kd1, LAI): canopy_irradiance => begin
        #I_c = I_cSun + I_cSh
        I(I0, K) = (1 - rho_cb) * I0 * (1 - exp(-K * LAI))
        I_tot = I(I0_dr, Kb1) + I(I0_df, Kd1)
    end ~ track(u"μmol/m^2/s" #= Quanta =#)

    # I_cSun: The irradiance absorbed by the sunlit fraction, de Pury and Farquhar (1997)
    # should this be the same os Qsl? 03/02/08 SK
    I_cSun(s, rho_cb, rho_cd, I0_dr, I0_df, Kb, Kb1, Kd1, LAI): canopy_sunlit_irradiance => begin
        I_c_sunlit = begin
            I0_dr * (1 - s) * (1 - exp(-Kb * LAI)) +
            I0_df * (1 - rho_cd) * (1 - exp(-(Kd1 + Kb) * LAI)) * Kd1 / (Kd1 + Kb) +
            I0_dr * ((1 - rho_cb) * (1 - exp(-(Kb1 + Kb) * LAI)) * Kb1 / (Kb1 + Kb) - (1 - s) * (1 - exp(-2Kb * LAI)) / 2)
        end
    end ~ track(u"μmol/m^2/s" #= Quanta =#)

    # I_cSh: The irradiance absorbed by the shaded fraction, de Pury and Farquhar (1997)
    I_cSh(I_c, I_cSun): canopy_shaded_irradiance => (I_c - I_cSun) ~ track(u"μmol/m^2/s" #= Quanta =#)

    ######
    # Q? #
    ######

    # sunlit_photon_flux_density(_sunlit_Q) ~ track
    # shaded_photon_flux_density(_shaded_Q) ~ track

    # Qtot: total irradiance (dir + dif) at depth L, simple empirical approach
    Q_tot(I0_tot, s, Kb, Kd; L): irradiance_Q_tot => begin
        I0_tot * exp(-sqrt(1 - s) * ((Kb + Kd) / 2) * L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # Qbt: total beam radiation at depth L
    Q_bt(I0_dr, s, Kb; L): irradiance_Q_bt => begin
        I0_dr * exp(-sqrt(1 - s) * Kb * L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # net diffuse flux at depth of L within canopy
    Q_d(I0_dr, s, Kd; L): irradiance_Q_d => begin
        I0_df * exp(-sqrt(1 - s) * Kd * L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # weighted average absorbed diffuse flux over depth of L within canopy
    # accounting for exponential decay, Campbell p261
    Q_dm(LAI, I0_df, s, Kd): irradiance_Q_dm => begin
        # Integral Qd / Integral L
        Q = I0_df * (1 - exp(-sqrt(1 - s) * Kd * LAI)) / (sqrt(1 - s) * Kd * LAI)
        isnan(Q) ? zero(Q) : Q
    end ~ track(u"μmol/m^2/s" #= Quanta =#)

    # unintercepted beam (direct beam) flux at depth of L within canopy
    Q_b(I0_dr, Kb; L): irradiance_Q_b => begin
        I0_dr * exp(-Kb * L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # mean flux density on sunlit leaves
    Q_sun(I0_dr, Kb, Q_sh): irradiance_Q_sunlit => begin
        I0_dr * Kb + Q_sh
    end ~ track(u"μmol/m^2/s" #= Quanta =#)

    # flux density on sunlit leaves at depth L
    Q_sun_at(I0_dr, Kb; L): irradiance_Q_sunlit_at => begin
        I0_dr * Kb + Q_sh_at(L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # mean flux density on shaded leaves over LAI
    Q_sh(Q_dm, Q_scm): irradiance_Q_shaded => begin
        # It does not include soil reflection
        Q_dm + Q_scm
    end ~ track(u"μmol/m^2/s" #= Quanta =#)

    # diffuse flux density on shaded leaves at depth L
    Q_sh_at(Q_d, Q_sc; L): irradiance_Q_shaded_at => begin
        # It does not include soil reflection
        Q_d(L) + Q_sc(L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # weighted average of Soil reflectance over canopy accounting for exponential decay
    Q_soilm(LAI, rho_soil, s, Kd, Q_soil): irradiance_Q_soilm => begin
        # Integral Qd / Integral L
        Q = Q_soil * rho_soil * (1 - exp(-sqrt(1 - s) * Kd * LAI)) / (sqrt(1 - s) * Kd * LAI)
        isnan(Q) ? zero(Q) : Q
    end ~ track(u"μmol/m^2/s" #= Quanta =#)

    # weighted average scattered radiation within canopy
    Q_scm(LAI, I0_dr, s, Kb): irradiance_Q_scm => begin
        # total beam including scattered absorbed by canopy
        #FIXME should the last part be multiplied by LAI like others?
        #TODO simplify by using existing variables (i.e. Q_bt, Q_b)
        total_beam = I0_dr * (1 - exp(-sqrt(1 - s) * Kb * LAI)) / (sqrt(1 - s) * Kb)
        # non scattered beam absorbed by canopy
        nonscattered_beam = I0_dr * (1 - exp(-Kb * LAI)) / Kb
        Q = (total_beam - nonscattered_beam) / LAI
        # Campbell and Norman (1998) p 261, Average between top (where scattering is 0) and bottom.
        #(self.Q_bt(LAI) - Q_b(LAI)) / 2
        isnan(Q) ? zero(Q) : Q
    end ~ track(u"μmol/m^2/s" #= Quanta =#)

    # scattered radiation at depth L in the canopy
    Q_sc(Q_bt, Q_b; L): irradiance_Q_sc => begin
        # total beam - nonscattered beam at depth L
        Q_bt(L) - Q_b(L)
    end ~ call(u"μmol/m^2/s" #= Quanta =#)

    # total PFD at the soil surface under the canopy
    Q_soil(LAI, Q_tot): irradiance_Q_soil => Q_tot(LAI) ~ track(u"μmol/m^2/s" #= Quanta =#)

    ###################
    # Leaf Area Index #
    ###################

    sunrisen(elevation_angle, minimum_elevation_angle=5u"°") => begin
        elevation_angle > minimum_elevation_angle
    end ~ flag

    # sunlit LAI assuming closed canopy; thus not accurate for row or isolated canopy
    LAI_sunlit(sunrisen, Kb, LAI): sunlit_leaf_area_index => begin
        sunrisen ? (1 - exp(-Kb * LAI)) / Kb : 0.
    end ~ track

    # shaded LAI assuming closed canopy
    LAI_shaded(LAI, LAI_sunlit): shaded_leaf_area_index => begin
        LAI - LAI_sunlit
    end ~ track

    # sunlit fraction of current layer
    sunlit_fraction(sunrisen, Kb; L) => begin
        sunrisen ? exp(-Kb * L) : 0.
    end ~ call

    shaded_fraction(sunlit_fraction; L) => begin
        1 - sunlit_fraction(L)
    end ~ call
end

ext: jl


text:
using Wild
using Test

@testset "Wild.jl" begin
    # Write your tests here.
end


text:
# We will always reprenset a point as an `SVector`. Some aliases are defined
# below for convenience, but otherwise there is not much to defined regarding
# the representation of points.

"""
    const Point1D
    const Point1D(x1)

A point in 1D space, stored in a StaticArray.
Point1D = SVector{1, Float64}.
"""
const Point1D = SVector{1, Float64}

"""
    const Point2D
    const Point2D(x1, x2)
    const Point2D(x::NTuple{2, Float64})

A point in 2D space, stored in a StaticArray.
Point2D = SVector{2, Float64}.
"""
const Point2D = SVector{2, Float64}

"""
    const Point3D
    const Point3D(x1, x2, x3)
    const Point3D(x::NTuple{3, Float64})

A point in 3D space, stored in a StaticArray.
Point3D = SVector{3, Float64}.
"""
const Point3D = SVector{3, Float64}

"""
    const ComplexPoint3D
    const ComplexPoint3D(x1, x2, x3)
    const ComplexPoint3D(x::NTuple{3, ComplexF64})

A complex 3D point, stored in a StaticArray.
ComplexPoint3D = SVector{3, ComplexF64}.
"""
const ComplexPoint3D = SVector{3, ComplexF64}

"""
    const ComplexPoint2D
    const ComplexPoint2D(x1, x2)
    const ComplexPoint2D(x::NTuple{2, ComplexF64})

A complex 2D point, stored in a StaticArray.
ComplexPoint2D = SVector{2, ComplexF64}.
"""
const ComplexPoint2D = SVector{2, ComplexF64}

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
using BrainForth: @bf, @run, stack
using Base.Test

@testset "BrainForth" begin

@testset "Basics" begin
  @test stack(@run([8, 5, dup]))  == [8, 5, 5]
  @test stack(@run([8, 5, swap])) == [5, 8]
  @test stack(@run([8, 5, drop])) == [8]

  @test stack(@run([0, !]))     == [1]
  @test stack(@run([8, !]))     == [0]
  @test stack(@run([8, 5, +]))  == [13]
  @test stack(@run([8, 5, -]))  == [3]
  @test stack(@run([8, 5, *]))  == [40]

  @test stack(@run([8, 5, [1, +], dip]))         == [9, 5]
  @test stack(@run([8, 5, nip]))                 == [5]
  @test stack(@run([8, 5, over]))                == [8, 5, 8]
  @test stack(@run([8, 5, [sq], keep]))          == [8, 25, 5]
  @test stack(@run([1, 2, 3, rotl]))             == [2, 3, 1]
  @test stack(@run([8, 5, [1, +], [1, -], bi]))  == [8, 6, 4]
  @test stack(@run([8, 5, [1, +], [1, -], bi_])) == [9, 4]
  @test stack(@run([8, 5, [1, +], bia]))         == [9, 6]

  @test stack(@run([5, [[1, +], call, [1, +], call], call])) == [7]

  @test stack(@run([10, 1, [1, +], [1, -], iff])) == [11]
  @test stack(@run([10, 0, [1, +], [1, -], iff])) == [ 9]
  @test stack(@run([10, 1, [[[1, +], call], [[1, -], call], iff], call])) == [11]
end

@testset "Recursion" begin
  @bf factorial = [dup, 1, ==, [dup, 1, -, factorial, *], unless]
  @test stack(@run([5, factorial])) == [120]

  @bf fib = [dup, [1, ==], [0, ==], bi, or,
                  [[1, -, fib], [2, -, fib], bi, +], unless]
  @test stack(@run([7, fib])) == [13]
end

@testset "Arrays" begin
  @test stack(@run [5, iota, sum]) == [15]
  @test stack(@run [5, iota, prod]) == [120]
  @test stack(@run [5, iota, [sq], map, sum]) == [55]
  @test stack(@run [5, iota, reverse]) == [1, 2, 3, 4, 5, 5]
  @test @run(["foo", "bar", swapvv]) == @run(["bar", "foo"])
end

end

max_stars_repo_name: JuliaTagBot/BrainForth

max_issues_repo_name: JuliaTagBot/BrainForth

max_forks_repo_name: JuliaTagBot/BrainForth


text:
using JuLIP:               decode_dict
using SHIPs.JacobiPolys:   Jacobi

import Base:   Dict, convert, ==
import JuLIP:  cutoff

export PolyTransform, rbasis, PolyCutoff1s, PolyCutoff2s,
       PolyTransformCut

abstract type DistanceTransform end

abstract type DistanceTransformCut <: DistanceTransform end

poly_trans(p, r0, r) = @fastmath(((1+r0)/(1+r))^p)

poly_trans_d(p, r0, r) = @fastmath((-p/(1+r0)) * ((1+r0)/(1+r))^(p+1))

poly_trans_inv(p, r0, x) = ( (1+r0)/(x^(1/p)) - 1 )

"""
Implements the distance transform
```
r -> ( (1+r0)/(1+r))^p
```

Constructor:
```
PolyTransform(p, r0)
```
"""
struct PolyTransform{TP, T} <: DistanceTransform
   p::TP
   r0::T
end

Dict(T::PolyTransform) =
   Dict("__id__" => "SHIPs_PolyTransform", "p" => T.p, "r0" => T.r0)

PolyTransform(D::Dict) = PolyTransform(D["p"], D["r0"])

convert(::Val{:SHIPs_PolyTransform}, D::Dict) = PolyTransform(D)

transform(t::PolyTransform, r::Number) = poly_trans(t.p, t.r0, r)

transform_d(t::PolyTransform, r::Number) = poly_trans_d(t.p, t.r0, r)

inv_transform(t::PolyTransform, x::Number) = poly_trans_inv(t.p, t.r0, x)

"""
Implements the distance transform
```
r -> ( (1+r0)/(1+r))^p - c0 - c1 (r - rcut)
```

Constructor:
```
PolyTransformCut(p, r0)
```
"""
struct PolyTransformCut{TP, T} <: DistanceTransformCut
   p::TP
   r0::T
   c0::T
   c1::T
   rcut::T
end

Dict(T::PolyTransformCut) =
   Dict("__id__" => "SHIPs_PolyTransform",
        "p" => T.p, "r0" => T.r0, "rcut" => rcut)

PolyTransformCut(D::Dict) = PolyTransformCut(D["p"], D["r0"], D["rcut"])

convert(::Val{:SHIPs_PolyTransformCut}, D::Dict) = PolyTransformCut(D)

transform(t::PolyTransformCut, r::Number) =
       (poly_trans(t.p, t.r0, r) + t.c0 + t.c1 * (r - t.rcut)) * (r < t.rcut)

transform_d(t::PolyTransformCut, r::Number) =
       (poly_trans_d(t.p, t.r0, r) + t.c1) * (r < t.rcut)

function PolyTransformCut(p, r0, rcut)
   c0 = - poly_trans(p, r0, rcut)
   c1 = - poly_trans_d(p, r0, rcut)
   return PolyTransformCut(p, r0, c0, c1, rcut)
end

cutoff(trans::PolyTransformCut) = trans.rcut

abstract type PolyCutoff end

"""
Implements the one-sided cutoff
```
r -> (x - xu)^p
```
Constructor:
```
PolyCutoff1s(p, rcut)
```
"""
struct PolyCutoff1s{P} <: PolyCutoff
   valP::Val{P}
   rl::Float64
   ru::Float64
   PolyCutoff1s(valP::Val{P}, rl::Real, ru::Real) where {P} = (
         ((P isa Integer) && (P > 0))  ? new{P}(valP, Float64(rl), Float64(ru))
                                       : error("P must be a positive integer") )
end

PolyCutoff1s(p::Integer, ru) = PolyCutoff1s(Val(Int(p)), 0.0, ru)
PolyCutoff1s(p::Integer, rl, ru) = PolyCutoff1s(Val(Int(p)), rl, ru)

Dict(C::PolyCutoff1s{P}) where {P} =
   Dict("__id__" => "SHIPs_PolyCutoff1s",
        "P" => P, "rl" => C.rl, "ru" => C.ru)
PolyCutoff1s(D::Dict) = PolyCutoff1s(D["P"], D["rl"],  D["ru"])
convert(::Val{:SHIPs_PolyCutoff1s}, D::Dict) = PolyCutoff1s(D)

# what happened to @pure ??? => not exported anymore
fcut(C::PolyCutoff1s{P}, r::T, x::T) where {P, T} =
      r < C.ru ? @fastmath( (1 - x)^P ) : zero(T)
fcut_d(C::PolyCutoff1s{P}, r::T, x::T) where {P, T} =
      r < C.ru ? @fastmath( - P * (1 - x)^(P-1) ) : zero(T)

"""
Implements the two-sided cutoff
```
r -> (x - xu)^p (x-xl)^p
```
Constructor:
```
PolyCutoff2s(p, rl, ru)
```
where `rl` is the inner cutoff and `ru` the outer cutoff.
"""
struct PolyCutoff2s{P} <: PolyCutoff
   valP::Val{P}
   rl::Float64
   ru::Float64
   PolyCutoff2s(valP::Val{P}, rl::Real, ru::Real) where {P} = (
         ((P isa Integer) && (P > 0))  ? new{P}(valP, Float64(rl), Float64(ru))
                                       : error("P must be a positive integer") )
end

PolyCutoff2s(p::Integer, rl, ru) = PolyCutoff2s(Val(Int(p)), rl, ru)

Dict(C::PolyCutoff2s{P}) where {P} =
   Dict("__id__" => "SHIPs_PolyCutoff2s", "P" => P,
        "rl" => C.rl, "ru" => C.ru)
PolyCutoff2s(D::Dict) = PolyCutoff2s(D["P"], D["rl"], D["ru"])
convert(::Val{:SHIPs_PolyCutoff2s}, D::Dict) = PolyCutoff2s(D)

fcut(C::PolyCutoff2s{P}, r::T, x) where {P, T} =
      C.rl < r < C.ru ? @fastmath( (1 - x^2)^P ) : zero(T)
fcut_d(C::PolyCutoff2s{P}, r::T, x) where {P, T} =
      C.rl < r < C.ru ? @fastmath( -2*P * x * (1 - x^2)^(P-1) ) : zero(T)

struct OneCutoff
   rcut::Float64
end
fcut(C::OneCutoff, r, x) = r < rcut ? one(r) : zero(r)
fcut_d(C::OneCutoff, r, x) = zero(r)

# Transformed Jacobi Polynomials
# ------------------------------
# these define the radial components of the polynomials

struct TransformedJacobi{T, TT, TM}
   J::Jacobi{T}
   trans::TT      # coordinate transform
   mult::TM       # a multiplier function (cutoff)
   rl::T          # lower bound r
   ru::T          # upper bound r
   tl::T          #  bound t(ru)
   tu::T          #  bound t(rl)
end

==(J1::TransformedJacobi, J2::TransformedJacobi) = (
   (J1.J == J2.J) &&
   (J1.trans == J2.trans) &&
   (J1.mult == J2.mult) &&
   (J1.rl == J2.rl) &&
   (J1.ru == J2.ru) )

TransformedJacobi(J, trans, mult, rl, ru) =
   TransformedJacobi(J, trans, mult, rl, ru,
                     transform(trans, rl), transform(trans, ru) )

Dict(J::TransformedJacobi) = Dict(
      "__id__" => "SHIPs_TransformedJacobi",
      "a" => J.J.α,
      "b" => J.J.β,
      "deg" => length(J.J) - 1,
      "rl" => J.rl,
      "ru" => J.ru,
      "trans" => Dict(J.trans),
      "cutoff" => Dict(J.mult),
      "skip0" => J.J.skip0
   )

TransformedJacobi(D::Dict) =
   TransformedJacobi(
      Jacobi(D["a"], D["b"], D["deg"],
             skip0 = haskey(D, "skip0")  ? D["skip0"]  : false),
      decode_dict(D["trans"]),
      decode_dict(D["cutoff"]),
      D["rl"],
      D["ru"]
   )

convert(::Val{:SHIPs_TransformedJacobi}, D::Dict) = TransformedJacobi(D)

Base.length(J::TransformedJacobi) = length(J.J)

cutoff(J::TransformedJacobi) = J.ru
transform(J::TransformedJacobi, r) = transform(J.trans, r)
transform_d(J::TransformedJacobi, r) = transform_d(J.trans, r)
fcut(J::TransformedJacobi, r, x) = fcut(J.mult, r, x)
fcut_d(J::TransformedJacobi, r, x) = fcut_d(J.mult, r, x)

SHIPs.alloc_B( J::TransformedJacobi{T}, args...) where {T} =
      Vector{T}(undef, length(J))

SHIPs.alloc_dB(J::TransformedJacobi{T}, args...) where {T} =
      Vector{T}(undef, length(J))

function eval_basis!(P, tmp, J::TransformedJacobi, r)
   N = length(J)-1
   @assert length(P) >= N+1
   # transform coordinates
   t = transform(J.trans, r)
   x = -1 + 2 * (t - J.tl) / (J.tu-J.tl)
   # evaluate the cutoff multiplier
   # the (J.tu-J.tl) / 2 factor makes the basis orthonormal
   # (just for the kick of it...)
   fc = fcut(J, r, x) * sqrt(abs(2 / (J.tu-J.tl)))
   if fc == 0
      fill!(P, 0.0)
   else
      # evaluate the actual Jacobi polynomials
      eval_basis!(P, nothing, J.J, x)
      for n = 1:N+1
         @inbounds P[n] *= fc
      end
   end
   return P
end

function eval_basis_d!(P, dP, tmp, J::TransformedJacobi, r)
   N = length(J)-1
   @assert length(P) >= N+1
   # transform coordinates
   t = transform(J.trans, r)
   x = -1 + 2 * (t - J.tl) / (J.tu-J.tl)
   dx = (2/(J.tu-J.tl)) * transform_d(J.trans, r)
   # evaluate the cutoff multiplier and chain rule
   fc = fcut(J, r, x) * sqrt(abs(2 / (J.tu-J.tl)))
   if fc ==  0
      fill!(P, 0.0)
      fill!(dP, 0.0)
   else
      fc_d = fcut_d(J, r, x) * sqrt(abs(2 / (J.tu-J.tl)))
      # evaluate the actual Jacobi polynomials + derivatives w.r.t. x
      eval_basis_d!(P, dP, nothing, J.J, x)
      for n = 1:N+1
         @inbounds p = P[n]
         @inbounds dp = dP[n]
         @inbounds P[n] = p * fc
         @inbounds dP[n] = (dp * fc + p * fc_d) * dx
      end
   end
   return dP
end

TransformedJacobi(maxdeg::Integer,
                  trans::DistanceTransform,
                  cut::PolyCutoff1s{P}) where {P} =
      TransformedJacobi( Jacobi(2*P,   0, maxdeg), trans, cut, cut.rl, cut.ru)

TransformedJacobi(maxdeg::Integer,
                  trans::DistanceTransform,
                  cut::PolyCutoff2s{P}) where {P} =
      TransformedJacobi( Jacobi(2*P, 2*P, maxdeg), trans, cut, cut.rl, cut.ru)

TransformedJacobi(maxdeg::Integer,
                  trans::DistanceTransformCut,
                  rl = 0.0) =
      TransformedJacobi( Jacobi(0, 0, maxdeg, skip0=true), trans,
                         OneCutoff(cutoff(trans)),
                         rl, cutoff(trans))

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
using StanSample, Tables, Test

# Testing

mat = [1 4.0 "7"; 2 5.0 "8"; 3 6.0 "9"]
stantbl = Tables.table(mat)

# first, create a MatrixTable from our matrix input
mattbl = Tables.table(mat)
# test that the MatrixTable `istable`
@test Tables.istable(typeof(mattbl))
# test that it defines row access
@test Tables.rowaccess(typeof(mattbl))
#@test Tables.rows(mattbl) === mattbl
# test that it defines column access
@test Tables.columnaccess(typeof(mattbl))
@test Tables.columns(mattbl) === mattbl
# test that we can access the first "column" of our matrix table by column name
@test mattbl.Column1 == [1,2,3]
# test our `Tables.AbstractColumns` interface methods
@test Tables.getcolumn(mattbl, :Column1) == [1,2,3]
@test Tables.getcolumn(mattbl, 1) == [1,2,3]
@test Tables.columnnames(mattbl) == [:Column1, :Column2, :Column3]
# now let's iterate our MatrixTable to get our first MatrixRow
matrow = first(mattbl)
#@test eltype(mattbl) == typeof(matrow)
# now we can test our `Tables.AbstractRow` interface methods on our MatrixRow
#@test matrow.Column1 == 1
#@test Tables.getcolumn(matrow, :Column1) == 1
#@test Tables.getcolumn(matrow, 1) == 1
#@test propertynames(mattbl) == propertynames(matrow) == [:Column1, :Column2, :Column3]

rt = [(a=1, b=4.0, c="7"), (a=2, b=5.0, c="8"), (a=3, b=6.0, c="9")]
ct = (a=[1,2,3], b=[4.0, 5.0, 6.0])

# let's turn our row table into a plain Julia Matrix object
mat = Tables.matrix(rt)
# test that our matrix came out like we expected
@test mat[:, 1] == [1, 2, 3]
@test size(mat) == (3, 3)
@test eltype(mat) == Any
# so we successfully consumed a row-oriented table,
# now let's try with a column-oriented table
mat2 = Tables.matrix(ct)
@test eltype(mat2) == Float64

# now let's take our matrix input, and make a column table out of it
tbl = Tables.table(mat) |> Tables.columntable
@test keys(tbl) == (:Column1, :Column2, :Column3)
@test tbl.Column1 == [1, 2, 3]
# and same for a row table
tbl2 = Tables.table(mat2) |> Tables.rowtable
@test length(tbl2) == 3
@test map(x->x.Column1, tbl2) == [1.0, 2.0, 3.0]

@test Tables.istable(tbl2) == true

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
# ---
# title: Faceting
# cover: assets/faceting.png
# description: Generating a grid of plots.
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

# ## Facet grid

df = let
    N = 100
    x0 = rand(1:10, N)
    i = rand(["α", "β"], N)
    j = rand(["a", "b", "c"], N)

    x = map(zip(x0, j)) do (xx, jj)
        shift = jj == "a" ? -2.9 : jj == "c" ? 2.9 : 0.0
        xx + shift
    end

    y = map(zip(x0, i)) do (xx, ii)
        shift = ii == "α" ? -3.9 : 3.9
        xx + 2 + shift + rand()
    end

    (; x, y, i, j)
end

plt = data(df) * mapping(:x, :y, row=:i, col=:j)

draw(plt)

# ## Facet grid with minimal axes linking needed to remove ticks

draw(plt, facet = (; linkxaxes = :minimal, linkyaxes = :minimal))

# ## Facet grid with unlinked x-axes

draw(plt, facet = (; linkxaxes = :none))

# ## Facet wrap

df = (x=rand(100), y=rand(100), l=rand(["a", "b", "c", "d", "e"], 100))
plt = data(df) * mapping(:x, :y, layout=:l)
draw(plt)

# ## Facet wrap with unlinked axes

draw(plt, facet = (; linkxaxes = :none, linkyaxes = :none))

# ## Adding traces to only some subplots

df1 = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
df2 = (x=[0, 1], y=[0.5, 0.5], i=fill("a", 2), j=fill("e", 2))
layers = data(df1) * visual(Scatter) + data(df2) * visual(Lines)
fg = draw(layers * mapping(:x, :y, col=:i, row=:j))

# save cover image #src
mkpath("assets") #src
save("assets/faceting.png", fg) #src


text:
#####
##### generate the output in the README
#####

using Fretboards, Crayons, ANSIColoredPrinters
using Fretboards.Intervals

fb = Fretboard(STANDARD_TUNING, 12);
let root = pitch"E"
    for (mark, semitones, color) in [("1", 0, :red),
                                     ("4", P4, :yellow),
                                     ("5", P5, :green)]
        annotate_matches!(x -> x.pitch ≂ root + semitones,
                          Annotation(mark, Crayon(foreground = color)), fb)
    end
end
buf = IOBuffer()
Fretboards.print_ascii(buf, fb, Fretboards.UNICODE_CHARS)
show(stdout, MIME"text/html"(), HTMLPrinter(buf))
end


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

"""
    UniformPartitioner(k, [shuffle])

A method for partitioning spatial data uniformly into `k` subsets
of approximately equal size. Optionally `shuffle` the data (default
to true).
"""
struct UniformPartitioner <: AbstractPartitioner
  k::Int
  shuffle::Bool
end

UniformPartitioner(k::Int) = UniformPartitioner(k, true)

function partition(object::AbstractSpatialObject{T,N},
                   partitioner::UniformPartitioner) where {N,T}
  npts = npoints(object)
  nset = partitioner.k

  @assert nset ≤ npts "number of subsets must be smaller than number of points"

  inds = partitioner.shuffle ? shuffle(1:npts) : collect(1:npts)
  subsets = collect(Iterators.partition(inds, npts ÷ nset))

  SpatialPartition(object, subsets)
end


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

"""
    SimplificationMethod

A method for simplifying geometries.
"""
abstract type SimplificationMethod end

"""
    simplify(geometry, method)

Simplify `geometry` with given `method`.
"""
function simplify end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("simplification/douglaspeucker.jl")


text:
# # Sum-of-Squares matrices

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Sum-of-Squares Matrices.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Sum-of-Squares Matrices.ipynb)
# **Adapted from**: Examples 3.77 of [BPT12]
#
# [BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
# *Semidefinite Optimization and Convex Algebraic Geometry*.
# Society for Industrial and Applied Mathematics, **2012**.

# ### Introduction

# Consider the symmetric polynomial matrix
# $$P(x) =
# \begin{bmatrix}
#     x^2 - 2x + 2 & x\\
#     x            & x^2
# \end{bmatrix}.$$
# We could like to know whether $P(x)$ is positive semidefinite for all $x \in \mathbb{R}$.

using Test #src
using DynamicPolynomials
@polyvar x
P = [x^2 - 2x + 2 x
            x     x^2]

# A sufficient condition for a symmetric polynomial matrix $P(x)$ to be positive semidefinite for all $x$ is to the existence of a matrix $M(x)$ such that $P(x) = M^\top(x) M(x)$. If such matrix $M$ exists, we say that the matrix is an \emph{sos matrix} (see [Definition 3.76, BPT13]).
# While determining whether $P(x)$ is positive semidefinite for all $x$, is NP-hard (checking nonnegativity of a polynomial is reduced to this problem for $1 \times 1$ matrices), checking whether $P(x)$ is an sos matrix is an sos program.

using SumOfSquares

# We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices.

import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)

model = SOSModel(solver)
mat_cref = @constraint(model, P in PSDCone())
optimize!(model)
@test termination_status(model) == MOI.OPTIMAL #src
termination_status(model) #!src

# While the reformulation of sos matrix to sos polynomial is rather simple, as explained in the "Sum-of-Squares reformulation" section below, there is a technical subtelty about the Newton polytope that if not handled correctly may result in an SDP of large size with bad numerical behavior. For this reason, it is recommended to model sos *matrix* constraints as such as will be shown in this notebook and not do the formulation manually unless there is a specific reason to do so.
#
# As we can verify as follows, only 3 monomials are used using the sos *matrix* constraint.

@test length(certificate_monomials(mat_cref)) == 3 #src
certificate_monomials(mat_cref) #!jl

# ### Sum-of-Squares reformulation
#
# One way to obtain the reduction to an sos program is to create an intermediate vector of variable $y$ and check whether the polynomial $p(x, y) = y^\top P(x) y$ is sos.
# However, special care is required when approximating the Newton polytope of $p(x, y)$.
# Indeed, for instance if the entries of $P(x)$ are quadratic forms then the Newton polytope of $p(x, y)$ is the cartesian product between the Newton polytope of $y^\top y$ and the Newton polytope of $x^\top x$.
# In other words, $p(x, y)$ belongs to a family of quartic forms called biquadratic forms.
# This fact is important when generating the semidefinite program so that only bilinear monomials are used.
# So if the cheap outer approximation is used (instead of the exact polyhedral computation) for the newton polytope then it is important to use a multipartite computation approximation of the newton polytope.
# The multipartie exact approach may perform worse compared to the unipartite exact in certain cases though.
# Consider for instance the polynomial matrix $\mathrm{Diag}(x_1^1, x_2^2)$ for which $p(x, y) = x_1^2y_1^2 + x_2^2y_2^2$.
# For this polynomial, only the monomials $x_1y_1$ and $x_2y_2$ are needed in the SDP reformulation while the multipartite approach,
# as it will compute the Newton polytope as a cartesian product, will not see the dependence between $x$ and $y$ in the presence of monomials and will also select the monomials $x_1y_2$ and $x_2y_1$.

@polyvar y[1:2]
p = vec(y)' * P * vec(y)

# We can see above that `p` is biquadratic polynomial in the variables `x` and `y`.
# Computing the Newton polytope with the cheap outer approximation
# without exploiting this multipartite structure gives the following 6 monomials.

X = monomials(p)
@test Certificate.monomials_half_newton_polytope(X, tuple(), apply_post_filter = false) == [x * y[1], x * y[2], y[1] * y[2], x, y[1], y[2]] #src
Certificate.monomials_half_newton_polytope(X, tuple(), apply_post_filter = false) #!jl

# Exploiting the multipartite structure gives 4 monomials.

@test Certificate.monomials_half_newton_polytope(X, ([x], y), apply_post_filter = false) == [x * y[1], x * y[2], y[1], y[2]] #src
Certificate.monomials_half_newton_polytope(X, ([x], y), apply_post_filter = false) #!jl

# In the example above, there were only 3 monomials, where does the difference come from ?
# Using the monomial basis, the only product of two monomials that is equal to
# `y[2]^2` is `y[2] * y[2]`. As `y[2]^2` is not a monomial of `p`, we can conclude
# that the diagonal entry with row and column corresponding to `y[2]` will be zero
# hence the whole column and row will be zero as well.
# Therefore, we can remove this monomial.

@test Certificate.monomials_half_newton_polytope(X, ([x], y)) == [x * y[1], x * y[2], y[1]] #src
Certificate.monomials_half_newton_polytope(X, ([x], y)) #!jl

# The same reasoning can be used for monomials `y[1]y[2]` and `x` therefore whether
# we exploit the multipartite structure or not, we get only 3 monomials thanks
# to this post filter.

@test Certificate.monomials_half_newton_polytope(X, tuple()) == [x * y[1], x * y[2], y[1]] #src
Certificate.monomials_half_newton_polytope(X, tuple()) #!jl

max_stars_repo_path: examples/Sum-of-Squares Matrices.jl

max_issues_repo_path: examples/Sum-of-Squares Matrices.jl

max_forks_repo_path: examples/Sum-of-Squares Matrices.jl


text:
using Zeros
if VERSION < v"0.7-"
    using Base.Test
else
    using Test
end

# Real
Z = Zero()
@test Zero(0) === Z
@test Z === Zero()
@test Z == Z
@test Z == 0
@test sizeof(Z) == 0
@test float(Z) === 0.0
@test -Z === Z
@test +Z === Z
@test 2*Z === Z
@test 2.0*Z === Z
@test Z*3 === Z
@test Z/2 === Z
@test Z-Z === Z
@test Z+Z === Z
@test Z*Z === Z
@test 1-Z === 1
@test 1.0-Z === 1.0
@test Z-1 == -1
@test Z+1 == 1
@test 2+Z == 2
@test (Z < Z) == false
@test (Z > Z) == false
@test Z <= Z
@test Z >= Z
@test Z < 3
@test Z > -2.0
@test ldexp(Z, 3) === Z
@test copysign(Z, 3) === Z
@test copysign(Z, -1) === Z
@test flipsign(Z, -1) === Z
@test sign(Z) === Z
@test round(Z) === Z
@test floor(Z) === Z
@test ceil(Z) === Z
@test trunc(Z) === Z
@test significand(Z) === Z
@test !isodd(Z)
@test iseven(Z)
@test string(Z) == "0̸"
@test fma(Z,1,Z) === Z
@test muladd(Z,1,Z) === Z
@test fma(Z,1,3) === 3
@test muladd(Z,1,3) === 3
@test fma(Z,Z,Z) === Z
@test muladd(Z,Z,Z) === Z
@test mod(Z, 3) === Z
@test rem(Z, 3) === Z
@test modf(Z) === (Z, Z)

#Complex
@test Z*im === Z
@test im*Z === Z
@test real(Z) === Z
@test imag(Z) === Z
@test Z+(2+3im) === 2+3im
@test (2+3im)+Z === 2+3im
@test Z-(2+3im) === -2-3im
@test (2+3im)-Z === 2+3im
@test Z*(2+3im) === Z
@test (2.0+3.0im)*Z === Z
@test Z/(2+3im) === Z

@test Complex(Z,Z) == Z
@test Complex(1,Z) == 1
@test Complex(1.0,Z) == 1.0
@test Complex(true,Z) == true

# testzero()
@test testzero(3) === 3
@test testzero(3+3im) === 3+3im
@test testzero(0) === Z
@test testzero(0+0im) === Z

# Array functions
for cf in [(T)->T, (T)->Complex{T}]
    for T in [UInt8, UInt16, UInt32, UInt64, UInt128, Int8, Int16, Int32, Int64, Int128, BigInt, Float16, Float32, Float64, BigFloat]
        A = ones(cf(T),10)
        @test zero!(A) === A
        @test all(A .== zero(cf(T)))
    end
end

# Test error handling
@test_throws InexactError Zero(1)
@test_throws InexactError convert(Zero, 0.1)
@test_throws DivideError 1.0/Z
@test_throws DivideError (1.0+2.0im)/Z
@test_throws DivideError Z/Z

# Test `MyComplex` example type
include("mycomplex_example.jl")
@test MyImaginary(2)*MyImaginary(3) === MyReal(-6)


text:
using MatrixColorings
using Ragged
using Base.Test

# g00: nothing
adj = RaggedArray(Vector{Int}[])
g00 = MatrixColorings.Graph(1:0, adj)
# g0: 1
adj = RaggedArray(Vector{Int}[[]])
g0 = MatrixColorings.Graph(1:1, adj)

# g01: 1 -- 2
adj = RaggedArray(Vector{Int}[[2],[1]])
g01 = MatrixColorings.Graph(1:2, adj)

# g1:
# 1 -- 2 -- 3 -- 5
#    \ | /
#      4
adj = RaggedArray(Vector{Int}[[2,4], [1,3,4], [2,4,5], [1,2,3], [3]])
g1 = MatrixColorings.Graph(1:5, adj)

# g2:
# 1 -- 2    3 -- 5
#    \ |
#      4
adj = RaggedArray(Vector{Int}[[2,4], [1,4], [5], [1,2], [3]])
g2 = MatrixColorings.Graph(1:5, adj)

# gn
n = 1000
adj = Vector{Int}[]
for i=1:n
    len = rand(1:n-1)
    push!(adj, randperm(n)[1:len])
end
adj = RaggedArray(adj)
gn = MatrixColorings.Graph(1:n, adj)

# test iteration over nearest neighbors
g,i = 1,1
for g in [g00, g0, g01, g1, g2]#, gn] #uncommenting leads to an error!?
    for i =1:length(g)
        @test Set(g.adjlist[:,i]) == Set(collect(MatrixColorings.Neigh1(g, i)))
    end
end

@test collect(MatrixColorings.Neigh2(g0, 1))==Any[]
@test collect(MatrixColorings.Neigh2(g01, 1))==Any[]

# # test iteration over 2nd neighbors
for g in [g1, g2]
    for v=1:length(g)
        nset = Int[]
        for n in MatrixColorings.Neigh1(g, v)
            append!(nset, g.adjlist[:,n])
        end
        nset = Set(nset)
        setdiff!(nset, v)
        @test Set{Int}(collect(MatrixColorings.Neigh2(g, v))) == nset
    end
end

###
# Jacobians
###

function check_color_jac(color, A)
    # Checks whether a coloring col is valid for A along dimension dim
    # Takes a bit of time for large A.
    dim = color.dim
    color = color.c
    dim2 = dim==1 ? 2 : 1
    out = zeros(eltype(A), maximum(color), size(A,dim2))
    if dim2==1
        for (j,c) in enumerate(color)
            for jj=1:size(A,1)
                out[c,jj] +=  A[jj,j]
            end
        end
    else
        for (j,c) in enumerate(color)
            for jj=1:size(A,2)
                out[c,jj] +=  A[j,jj]
            end
        end
    end
    @test maximum(out)<=1
    out
end
if VERSION < v"0.4.0-dev"
    # TODO: remove these two for 0.4
    rowvals(S::SparseMatrixCSC) = S.rowval
    nzrange(S::SparseMatrixCSC, col::Integer) = S.colptr[col]:(S.colptr[col+1]-1)
end
function check_color_jac(color, A::SparseMatrixCSC)
    # Checks whether a coloring col is valid for A along dimension dim
    # Takes a bit of time for large A.
    dim = color.dim
    color = color.c
    if dim==1
        A = A.'
    end
    out = zeros(eltype(A), maximum(color), size(A,1))
    for (j,c) in enumerate(color)
        for jj in nzrange(A, j)
            out[c,rowvals(A)[jj]] += nonzeros(A)[jj]
        end
    end
    @test maximum(out)<=1
    out
end

function check_jac(A, cols, adj)
    bi = MatrixColorings.matrix2bipartie_graph(A)
    @test bi.v1[end] == size(A,1)
    @test length(bi.v1)==size(A,1)
    @test length(bi.v2)==size(A,2)
    @test bi.adjlist==adj
    # color graph
    col,i = 1,1
    for dim=1:2
        col = MatrixColorings.partial_dist2_coloring(bi,dim)
        @test col.dim==dim
        @test length(col.c)==size(A,dim)
        @test col.c==cols[dim]
        check_color_jac(col, A)
        S=color_jac(A,dim)
        s=sum(S,dim)
        @test maximum(s)==minimum(s)==1
        B = dim==1 ? S*A : A*S
    end
end
println("warming up @time:")
@time 1
function check_jac(A)
    # tests when colors and adjlist are unknown. Assumes all non-zeros of A are ==1.
    # Also does timeings

    # warm-ups for timings below
    AA = A[1:4,1:5]
    bi = MatrixColorings.matrix2bipartie_graph(AA)
    col = MatrixColorings.partial_dist2_coloring(bi,2)
    S=color_jac(AA,2)
    jacout = deepcopy(AA)
    B = AA*S
    recover_jac!(jacout, B, S)

    # do tests
    println("@time for constructing a bipartite graph from $(n)x$n+1 sparse matrix:")
    @time bi = MatrixColorings.matrix2bipartie_graph(A)
    @test bi.v1[end] == size(A,1)
    @test length(bi.v1)==size(A,1)
    @test length(bi.v2)==size(A,2)
    # color graph
    out = 1
    j,jj,c,dim = 1,1,1,1
    for dim=1:2
        dim2 = dim==1 ? 2 : 1
        println("@time for coloring that graph $(["row","columns"][dim])-wise:")
        gc()
        @time col = MatrixColorings.partial_dist2_coloring(bi,dim)
        @test length(col.c)==size(A,dim)
        # sum all columns of same color, no entry should be >1
        check_color_jac(col, A)
        println("@time for creating a seedmatrix $(["row","columns"][dim])-wise:")
        @time S=color_jac(A,dim)
        s=sum(S,dim)
        @test maximum(s)==minimum(s)==1
        B = dim==1 ? S*A : A*S
        if dim==2
            jacout = deepcopy(A)
            fill!(nonzeros(jacout), 0)
            println("@time for recovery:")
            @time recover_jac!(jacout, B, S)
            @test isequal_sparse(jacout,A)
        end
    end
end

# generate graph:
A = [0.0374053  0.0       0.1
     0.0        0.140752  0.380924]
cols = Vector{Int}[[1,2], [1,1,2]]
adj = RaggedArray(Vector{Int}[[1+2, 3+2], [2+2,3+2],
                              [1], [2], [1,2]])
check_jac(A,cols,adj)

# generate graph:
A = [0.0374053  1.0       0.1
     1.0        0.140752  0.380924]
cols = Vector{Int}[[1,2], [1,2,3]]
adj = RaggedArray(Vector{Int}[[1+2, 2+2, 3+2], [1+2,2+2,3+2],
                              [1,2], [1,2], [1,2]])
check_jac(A,cols,adj)

# Generate a big, random graph

# A = sparse([0.0374053  0.0       0.1
#             0.0        0.140752  0.380924])

# S=color_jac(A,2)
# B = A*S
# jacout = deepcopy(A)
# fill!(nonzeros(jacout), 0)
# recover_jac!(jacout, B, S)
# jacout==A
# error("dsaf")

bi = MatrixColorings.matrix2bipartie_graph(sprand(1,1+1, 0.001)) # warmup
n = 10^4
A = sprand(n,n+1, 0.001)
A.nzval[A.nzval.>0] = 1.0
check_jac(A)


text:
"""
    print the err text
"""
_printerr(io::IO, err) = print(io, _err_str(err))
_printerr(err) = _printerr(stdout, err)

function _dict(d = Dict{String, Any}(); kwargs...)
	for (k, v) in kwargs
		d[string(k)] = v
	end
	return d
end


text:
using Dates

function hfun_bar(vname)
  val = Meta.parse(vname[1])
  return round(sqrt(val), digits=2)
end

function hfun_m1fill(vname)
  var = vname[1]
  return pagevar("index", var)
end

html_img_click(src, alt) = "[![$(Franklin.htmlesc(alt))]($src)]($src)"

function lx_baz(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return uppercase(brace_content)
end

function lx_figclick(lx, _)
  # keep this first line
  brace_content = Franklin.content(lx.braces[1]) # input string
  alt, rpath = strip.(split(brace_content, ','))
  path  = Franklin.parse_rpath(rpath; canonical=false, code=true)
  return html_img_click(path, alt)
end

function hfun_show_tags()
  tags = locvar("tags")
  html_tags = map(x->"$x", tags)
  return join(html_tags, ", ")
end

function hfun_formatted_date()
  date = locvar("date")
  return Dates.format(date, "dd.mm.Y")
end


text:
struct SemCor{S}
    filepaths::Vector{S}
end

function SemCor(dirpath)
    @assert(isdir(dirpath), dirpath)

    if "brown1" ∈ readdir(dirpath) && "brown2" ∈ readdir(dirpath)
        tagfile_paths = joinpath.(dirpath, ["brown1", "brown2"], "tagfiles")

        paths = joinpath.(tagfile_paths[1], readdir(tagfile_paths[1]))
        append!(paths, joinpath.(tagfile_paths[2], readdir(tagfile_paths[2])))

    elseif "tagfiles" ∈ readdir(dirpath)
        paths = joinpath.(dirpath,"tagfiles", readdir(joinpath(dirpath, "tagfiles")))
    else
        paths = joinpath.(dirpath, readdir(dirpath))
    end
    SemCor(paths)
end

SemCor() = SemCor(datadep"SemCor 3.0")

MultiResolutionIterators.levelname_map(::Type{SemCor}) = [
    :doc=>1, :contextfile=>1, :context=>1,
    :para=>2, :paragraph=>2,
    :sent=>3, :sentence=>3,
    :word=>4, :token=>4,
    :char=>5, :character=>5
	]

function parse_sense_annotated_word(line::AbstractString)
    captures = match(r"<wf cmd=done.* pos=(.*) lemma=(.*) wnsn=(.*) lexsn=(\d.*:\d*).*>(.*).*</wf>", line).captures
	pos, lemma, wnsn, lexsn, word = captures
	if ';' in wnsn
		# Discard Extra Senses
		wnsn = split(wnsn, ';') |> first
		lexsn = split(lexsn, ';') |> first
	end

	SenseAnnotatedWord(pos, lemma, parse(Int, wnsn), lexsn, word)
end

function parse_tagged_word(line::AbstractString)
    captures = match(r"<wf cmd=.* pos=(.*).*>(.*)</wf>", line).captures
    PosTaggedWord(captures...)
end

function parse_punc(line::AbstractString)
    captures = match(r"<punc>(.*)</punc>", line).captures
    PosTaggedWord("PUNC", first(captures))
end

function parse_semcorfile(filename)
    local sent
    local para
    paras = @NestedVector(TaggedWord,3)()
    context = Document(intern(basename(filename)), paras)

    # structure
    function new_paragraph(line)
        para = @NestedVector(TaggedWord,2)()
        push!(paras, para)
    end

    function new_sentence(line)
        sent = @NestedVector(TaggedWord,1)()
        push!(para, sent)
    end

    # words
    get_tagged(line) = push!(sent, parse_tagged_word(line))
    get_sense_annotated(line) = push!(sent, parse_sense_annotated_word(line))
	get_punc(line) = push!(sent, parse_punc(line))

    # Parse
    subparsers = [
        "<wf cmd=tag"=>          get_tagged,
        "<wf cmd=ignore"=>       get_tagged,
        "<wf cmd=done"=>         line -> contains(line,"lemma=") ? get_sense_annotated(line) : get_tagged(line),
        "<punc>"=>               get_punc,
        "<context"=>             ignore,
        "</context" =>           ignore,
        "<p" =>                  new_paragraph,
        "</p" =>                 ignore,
        "<s" =>                  new_sentence,
        "</s" =>                 ignore
    ]

    apply_subparsers(filename,subparsers)

    return context
end

function load(corpus::SemCor, doc_buffersize=16)
    Channel(;ctype=Document{@NestedVector(TaggedWord, 3), String}, csize=doc_buffersize) do ch
        for fn in corpus.filepaths
            doc = parse_semcorfile(fn)
            put!(ch, doc)
        end
    end
end


text:
struct ZygoteRuleConfig{CTX<:AContext} <: RuleConfig{Union{HasReverseMode,NoForwardsMode}}
  context::CTX
end
ZygoteRuleConfig() = ZygoteRuleConfig(Context())

_is_rrule_redispatcher(m::Method) = m.sig == Tuple{typeof(rrule), RuleConfig, Vararg}

"""
  has_chain_rrule(T)

For a type-tuple `T` e.g. `Tuple{typeof(f), Int, Float64}`, checks if there is a `rrule` defined for it.
Excluding the generic fallback.
The first return value is `true` if the `rrule` exists, `false` otherwise.
If it does not, then the second argument is a list of edges to attach to the CodeInfo for a generated function,
such that if a suitable rule is defined later, the generated function will recompile.
"""
function has_chain_rrule(T)
  config_T, arg_Ts = Iterators.peel(T.parameters)
  configured_rrule_m = meta(Tuple{typeof(rrule), config_T, arg_Ts...})
  if _is_rrule_redispatcher(configured_rrule_m.method)
    # The config is not being used:
    # it is being redispatched without config, so we need the method it redispatches to
    rrule_m = meta(Tuple{typeof(rrule), arg_Ts...})
    # Thus any no_rrule that might apply must also not have a config because if there was a
    # no_rrule with a config that applied then there would also be a rrule with config that applied
    no_rrule_m = meta(Tuple{typeof(ChainRulesCore.no_rrule), arg_Ts...})
  else
    # Not being redispatched: it does have a config
    rrule_m = configured_rrule_m
    # Thus any no_rrule that might apply must also have a config because if it applied
    # it will be identical, and if it doesn't we don't care what it is.
    no_rrule_m = meta(Tuple{typeof(ChainRulesCore.no_rrule), config_T, arg_Ts...})
  end

  # To understand why we only need to check if the sigs match between no_rrule_m and rrule_m
  # in order to decide if to use, one must consider the following facts:
  # - for every method in `no_rrule` there is a identical one in `rrule` that returns nothing
  # - this includes the general fallback `rrule(::Any...)=nothing`.
  # - a configured rrule/no_rrule is always more specific than a otherwise equivalent unconfigured rrule/no_rrule
  #
  # Consider the following truth table, for what can occur:
  # rrule: fallback, no_rrule: fallback =>  matches => do not use rrule.
  # rrule: specific, no_rrule: fallback => !matches => do use rrule, as haven't opted out.
  # rrule: fallback, no_rrule: specific =>  IMPOSSIBLE, every no_rule is identical to some rrule
  # rrule: specific, no_rrule: specific =>  matches => do not use rrule as opted out
  # rrule: specific, no_rrule: general  => !matches => do use rrule as a more specific rrule takes preciedent over more general opted out
  # rrule: general , no_rrule: specific =>  IMPOSSIBLE, every no_rule us identical to some rrule so can't have a more general rrule being hit, as the specific one would hit first
  #
  # Note that the fallback cases are the same outcome as the general cases as fallback is just most general.
  # It can be seen that checking if it matches is the correct way to decide if we should use the rrule or not.

  do_not_use_rrule = matching_cr_sig(no_rrule_m, rrule_m)
  if do_not_use_rrule
    # Return instance for configured_rrule_m as that will be invalidated
    # directly if configured rule added, or indirectly if unconfigured rule added
    # Do not need an edge for `no_rrule` as no addition of methods to that can cause this
    # decision to need to be revisited (only changes to `rrule`), since we are already not
    # using the rrule, so not using more rules wouldn't change anything.
    return false, configured_rrule_m.instance
  else
    # Otherwise found a rrule, no need to add any edges for `rrule`, as it will generate
    # code with natural edges if a new method is defined there.
    # We also do not need an edge to `no_rrule`, as any time a method is added to `no_rrule`
    # a corresponding method is added to `rrule` (to return `nothing`), thus we will already
    # be revisiting this decision when a new opt-out is added.
    return true, nothing
  end
end

matching_cr_sig(t, s) = matching_cr_sig(t.method.sig, s.method.sig)
matching_cr_sig(::DataType, ::UnionAll) = false
matching_cr_sig(::UnionAll, ::DataType) = false
matching_cr_sig(t::Type, s::Type) = type_tuple_tail(t) == type_tuple_tail(s)

type_tuple_tail(d::DataType) = Tuple{d.parameters[2:end]...}
function type_tuple_tail(d::UnionAll)
    body = Base.unwrap_unionall(d)
    body_tt = type_tuple_tail(body)
    return Base.rewrap_unionall(body_tt, d)
end

"""
    is_kwfunc(sigt...)

Determines if `sigt` is the type signature of a kwfunction.
Each element of `sigt` should be a type.
Either the first 3 types are a kwfunc type, a NamedTuple and the matching base function type,
or the first argument is the base function type and it is not a kwfunction.
the remaining types in `sigt` are the types of the argument.

"""
is_kwfunc(::Vararg{Any}) = false
is_kwfunc(k, ::Type{<:NamedTuple}, f, args...) = k===Core.kwftype(f)

"""
    wrap_chainrules_output(x)

Convert `x` from the differentials types ChainRules uses to the format Zygote uses internally.
"""
@inline wrap_chainrules_output(x) = x
@inline wrap_chainrules_output(x::AbstractThunk) = wrap_chainrules_output(unthunk(x))  # For now we are just not going to deal with thunks
@inline wrap_chainrules_output(x::Tuple) = map(wrap_chainrules_output, x)
# Zygote convention: even if many AbstractZero partials (i.e. multi-input function), make just 1 nothing.
@inline wrap_chainrules_output(x::Tuple{Vararg{ChainRules.AbstractZero}}) = nothing
@inline wrap_chainrules_output(x::ChainRules.AbstractZero) = nothing
for T_outer in (:Tuple, :NamedTuple)
  # we create separate methods rather than using a `Union` + an `if` so that we avoid a
  # branch that changes output type, because nested AD on that kinda thing makes Zygote less
  # than happy.
  @eval @inline function wrap_chainrules_output(x::ChainRules.Tangent{P, T}) where {P, T<:$T_outer}
    xp = map(wrap_chainrules_output, canonicalize(x))
    ChainRulesCore.backing(xp)  # this is accessing ChainRulesCore internals, but it is prob safe enough, and it is fastest
  end
end

"""
    wrap_chainrules_input(x)

Convert `x` from the format Zygote uses internally to differentials types ChainRules uses.
"""
@inline wrap_chainrules_input(x) = x
@inline wrap_chainrules_input(::Nothing) = ChainRules.ZeroTangent()
@inline function wrap_chainrules_input(xs::Union{Tuple, NamedTuple})
  xp = map(wrap_chainrules_input, xs)
  ChainRules.Tangent{Any, typeof(xp)}(xp)
end

"""
  ZBack{F}(back) <: Function

Wrapper for a ChainRules pullback `back`, that causes it to follow Zygote conventions.
(A functor here is used rather than a closure to avoid boxing issues);
"""
struct ZBack{F} <: Function
  back::F
end
@inline (s::ZBack)(dy) = wrap_chainrules_output(s.back(wrap_chainrules_input(dy)))
# `nothing->nothing` can be deleted after https://github.com/FluxML/Zygote.jl/issues/603
# though it might be worth keeping as a performance optimization (benchmarking pending)
@inline (s::ZBack)(::Nothing) = nothing

"""
    chain_rrule(config, f, args...)

Returns a the (primal) value of `f(args...)` and a pullback, by invoking `ChainRulesCore.rrule(f, args...)`.
The pullback is appropriately wrapped up to follow Zygote conventions.
"""
@inline function chain_rrule(config, f, args...)
  y, back = rrule(config, f, args...)
  return y, ZBack(back)
end

"""
  chain_rrule_kw(config, kwf, kwargs, f, args...)

As per [`chain_rrule`](@ref) but with support for kwargs.
`kwf` should be the kwfunc matching to `f`, and `kwargs` are a `NamedTuple` of keyword arguments.
"""
@inline function chain_rrule_kw(config, kwf, kwargs, f, args...)
  y, back = rrule(config, f, args...; kwargs...)
  function kw_zpullback(dy)
    dxs = ZBack(back)(dy)
    if dxs === nothing  # if dxs is nothing, then all partiaols are nothing
      # Zygote convention is a single nothing no mather how partials, if all are nothing
      return nothing
    else
      return (nothing, nothing, dxs...)  # first two nothings are for kwfunc and kwargs
    end
  end
  return y, kw_zpullback
end

function ChainRulesCore.rrule_via_ad(config::ZygoteRuleConfig, f, args...)
    y, pb = _pullback(config.context, f, args...)
    ad_pullback(Δ) = zygote2differential(pb(wrap_chainrules_output(Δ)), (f, args...))
    return y, ad_pullback
end

"""
    zygote2differential(x)

Convert input `x` from the Zygote format to the ChainRules differential types.
"""
zygote2differential(x, primal) = z2d(x, primal)
zygote2differential(::Nothing, ::Any) = NoTangent()
zygote2differential(t::Tuple, primal::Tuple) = map(z2d, t, primal)
zygote2differential(t::Tuple, primal) = (@warn "primal should be a tuple, not $primal"; return t)
z2d(x, ::Any) = x
z2d(::Nothing, ::Any) = NoTangent()
z2d(a::AbstractArray{<:Number}, primal::AbstractArray{T}) where T = a
z2d(a::AbstractArray, primal::AbstractArray{T}) where T = z2d.(a, primal)
# Note: this should never be hit if we are converting things right, but it seems to be
# happening in the wild for sufficiently weird functions/types.
# This fixes most (all?) cases, but it would be good to find what we miss.
z2d(x::Union{AbstractZero, Tangent}, ::Any) = return x
function z2d(t::Tuple, primal::Tuple)
  tp::Tuple = map(z2d, t, primal)
  primal_type = typeof(primal)
  return canonicalize(Tangent{primal_type, typeof(tp)}(tp))
end

function z2d(t::NamedTuple, primal)
  primal_type = typeof(primal)
  fnames = fieldnames(primal_type)
  complete_t = NamedTuple{fnames}(fn in keys(t) ? t[fn] : nothing for fn in fnames)
  primals = NamedTuple{fnames}(getfield(primal, fn) for fn in fnames)
  tp::NamedTuple = map(z2d, complete_t, primals)
  return canonicalize(Tangent{primal_type, typeof(tp)}(tp))
end


text:
# deprecations
@deprecate iparameters(args...) getiparams(args...)
@deprecate setiparameters!(args...) setiparams!(args...)
@deprecate niparameters(args...) niparams(args...)
@deprecate parameter_type(args...) parameters_eltype(args...)


text:
# De Marchi, Themelis, "Proximal gradient algorithms under local Lipschitz
# gradient continuity: a convergence and robustness analysis of PANOC",
# arXiv:2112.13000 (2021).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

"""
    PANOCplusIteration(; <keyword-arguments>)

Instantiate the PANOCplus algorithm (see [1]) for solving optimization problems
of the form

    minimize f(Ax) + g(x),

where `f` is locally smooth and `A` is a linear mapping (for example, a matrix).

# Arguments
- `x0`: initial point.
- `f=Zero()`: smooth objective term.
- `A=I`: linear operator (e.g. a matrix).
- `g=Zero()`: proximable objective term.
- `Lf=nothing`: Lipschitz constant of the gradient of x ↦ f(Ax).
- `gamma=nothing`: stepsize to use, defaults to `1/Lf` if not set (but `Lf` is).
- `adaptive=false`: forces the method stepsize to be adaptively adjusted.
- `minimum_gamma=1e-7`: lower bound to `gamma` in case `adaptive == true`.
- `max_backtracks=20`: maximum number of line-search backtracks.
- `directions=LBFGS(5)`: strategy to use to compute line-search directions.

# References
1. De Marchi, Themelis, "Proximal gradient algorithms under local Lipschitz gradient continuity: a convergence and robustness analysis of PANOC", arXiv:2112.13000 (2021).
"""
Base.@kwdef struct PANOCplusIteration{R,Tx,Tf,TA,Tg,TLf,Tgamma,D}
    f::Tf = Zero()
    A::TA = I
    g::Tg = Zero()
    x0::Tx
    alpha::R = real(eltype(x0))(0.95)
    beta::R = real(eltype(x0))(0.5)
    Lf::TLf = nothing
    gamma::Tgamma = Lf === nothing ? nothing : (alpha / Lf)
    adaptive::Bool = gamma === nothing
    minimum_gamma::R = real(eltype(x0))(1e-7)
    max_backtracks::Int = 20
    directions::D = LBFGS(5)
end

Base.IteratorSize(::Type{<:PANOCplusIteration}) = Base.IsInfinite()

Base.@kwdef mutable struct PANOCplusState{R,Tx,TAx,TH}
    x::Tx             # iterate
    Ax::TAx           # A times x
    f_Ax::R           # value of smooth term
    grad_f_Ax::TAx    # gradient of f at Ax
    At_grad_f_Ax::Tx  # gradient of smooth term
    gamma::R          # stepsize parameter of forward and backward steps
    y::Tx             # forward point
    z::Tx             # forward-backward point
    g_z::R            # value of nonsmooth term (at z)
    res::Tx           # fixed-point residual at iterate (= x - z)
    H::TH             # variable metric
    tau::R = zero(gamma)
    x_prev::Tx = similar(x)
    res_prev::Tx = similar(x)
    d::Tx = similar(x)
    Az::TAx = similar(Ax)
    grad_f_Az::TAx = similar(Ax)
    At_grad_f_Az::Tx = similar(x)
end

f_model(iter::PANOCplusIteration, state::PANOCplusState) = f_model(state.f_Ax, state.At_grad_f_Ax, state.res, iter.alpha / state.gamma)

function Base.iterate(iter::PANOCplusIteration{R}) where {R}
    x = copy(iter.x0)
    Ax = iter.A * x
    grad_f_Ax, f_Ax = gradient(iter.f, Ax)
    gamma = iter.gamma === nothing ? iter.alpha / lower_bound_smoothness_constant(iter.f, iter.A, x, grad_f_Ax) : iter.gamma
    At_grad_f_Ax = iter.A' * grad_f_Ax
    y = x - gamma .* At_grad_f_Ax
    z, g_z = prox(iter.g, y, gamma)
    state = PANOCplusState(
        x=x, Ax=Ax, f_Ax=f_Ax, grad_f_Ax=grad_f_Ax, At_grad_f_Ax=At_grad_f_Ax,
        gamma=gamma, y=y, z=z, g_z=g_z, res=x-z, H=initialize(iter.directions, x),
    )
    if (iter.gamma === nothing || iter.adaptive == true)
        state.gamma, state.g_z, _, _ = backtrack_stepsize!(
            state.gamma, iter.f, iter.A, iter.g,
            state.x, state.f_Ax, state.At_grad_f_Ax, state.y, state.z, state.g_z, state.res,
            state.Az, state.grad_f_Az,
            alpha = iter.alpha, minimum_gamma = iter.minimum_gamma,
        )
    else
        mul!(state.Az, iter.A, state.z)
        gradient!(state.grad_f_Az, iter.f, state.Az)
    end
    mul!(state.At_grad_f_Az, adjoint(iter.A), state.grad_f_Az)
    return state, state
end

set_next_direction!(::QuasiNewtonStyle, ::PANOCplusIteration, state::PANOCplusState) = mul!(state.d, state.H, -state.res)
set_next_direction!(::NoAccelerationStyle, ::PANOCplusIteration, state::PANOCplusState) = state.d .= .-state.res
set_next_direction!(iter::PANOCplusIteration, state::PANOCplusState) = set_next_direction!(acceleration_style(typeof(iter.directions)), iter, state)

update_direction_state!(::QuasiNewtonStyle, ::PANOCplusIteration, state::PANOCplusState) = update!(state.H, state.x - state.x_prev, state.res - state.res_prev)
update_direction_state!(::NoAccelerationStyle, ::PANOCplusIteration, state::PANOCplusState) = return
update_direction_state!(iter::PANOCplusIteration, state::PANOCplusState) = update_direction_state!(acceleration_style(typeof(iter.directions)), iter, state)

reset_direction_state!(::QuasiNewtonStyle, ::PANOCplusIteration, state::PANOCplusState) = reset!(state.H)
reset_direction_state!(::NoAccelerationStyle, ::PANOCplusIteration, state::PANOCplusState) = return
reset_direction_state!(iter::PANOCplusIteration, state::PANOCplusState) = reset_direction_state!(acceleration_style(typeof(iter.directions)), iter, state)

function Base.iterate(iter::PANOCplusIteration{R}, state::PANOCplusState) where R
    f_Az, a, b, c = R(Inf), R(Inf), R(Inf), R(Inf)

    # store iterate and residual for metric update later on
    state.x_prev .= state.x
    state.res_prev .= state.res

    # compute FBE
    FBE_x = f_model(iter, state) + state.g_z

    sigma = iter.beta * (0.5 / state.gamma) * (1 - iter.alpha)
    tol = 10 * eps(R) * (1 + abs(FBE_x))
    threshold = FBE_x - sigma * norm(state.res)^2 + tol

    tau_backtracks = 0
    can_update_direction = true

    while true

        if can_update_direction
            # compute direction
            set_next_direction!(iter, state)
            # backtrack tau 1 → 0
            state.tau = R(1)
            state.x .= state.x_prev .+ state.d
        else
            state.x .= (1 - state.tau) * (state.x_prev .- state.res_prev) + state.tau * (state.x_prev .+ state.d)
        end

        mul!(state.Ax, iter.A, state.x)
        state.f_Ax = gradient!(state.grad_f_Ax, iter.f, state.Ax)
        mul!(state.At_grad_f_Ax, adjoint(iter.A), state.grad_f_Ax)

        state.y .= state.x .- state.gamma .* state.At_grad_f_Ax
        state.g_z = prox!(state.z, iter.g, state.y, state.gamma)
        state.res .= state.x .- state.z

        f_Az_upp = f_model(iter, state)

        if (iter.gamma === nothing || iter.adaptive == true)
            mul!(state.Az, iter.A, state.z)
            f_Az = gradient!(state.grad_f_Az, iter.f, state.Az)
            tol = 10 * eps(R) * (1 + abs(f_Az))
            if f_Az > f_Az_upp + tol && state.gamma >= iter.minimum_gamma
                state.gamma *= 0.5
                if state.gamma < iter.minimum_gamma
                    @warn "stepsize `gamma` became too small ($(state.gamma))"
                end
                can_update_direction = true
                reset_direction_state!(iter, state)
                continue
            end
        else
            mul!(state.Az, iter.A, state.z)
            f_Az = gradient!(state.grad_f_Az, iter.f, state.Az)
        end
        mul!(state.At_grad_f_Az, adjoint(iter.A), state.grad_f_Az)

        FBE_x_new = f_Az_upp + state.g_z
        if FBE_x_new <= threshold
            # update metric
            update_direction_state!(iter, state)
            return state, state
        end
        state.tau *= 0.5
        if tau_backtracks > iter.max_backtracks
            @warn "stepsize `tau` became too small ($(state.tau))"
            return nothing
        end
        tau_backtracks += 1
        can_update_direction = false

    end

end

# Solver

struct PANOCplus{R, K}
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int
    kwargs::K
end

function (solver::PANOCplus)(x0; kwargs...)
    stop(state::PANOCplusState) = norm((state.res / state.gamma) - state.At_grad_f_Ax + state.At_grad_f_Az, Inf) <= solver.tol
    disp((it, state)) = @printf(
        "%5d | %.3e | %.3e | %.3e\n",
        it,
        state.gamma,
        norm(state.res, Inf) / state.gamma,
        state.tau,
    )
    iter = PANOCplusIteration(; x0=x0, solver.kwargs..., kwargs...)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.z, num_iters
end

PANOCplus(; maxit=1_000, tol=1e-8, verbose=false, freq=10, kwargs...) =
    PANOCplus(maxit, tol, verbose, freq, kwargs)


text:
export add_walsh, add_residue

function add_walsh(mgr::Ptr{Manager}, f_pp::Ptr{Ptr{Node}}, g_pp::Ptr{Ptr{Node}}, n::Int)
    result = Cudd_addWalsh(mgr, f_pp, g_pp, n)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_residue(mgr::Ptr{Manager}, n::Int, m::Int, options::Int, top::Int)
    result = Cudd_addResidue(mgr, n, m, options, top)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end


text:
export getPaddedMesh, jInvObj

"""
function getPaddedMesh(M::RegularMesh)

return regular mesh padded by h/2 on all sided and with
one cell added, so that the cell-centered points on M
are nodal points on the padded mesh.
"""
function getPaddedMesh(M::RegularMesh)
	domainPad = copy(M.domain);
	domainPad[1:2:end] -= M.h/2
	domainPad[2:2:end] += M.h/2
	return getRegularMesh(domainPad,M.n.+1)
end

"""
function jInvObj(mc,pMis,pInv)

evaluates the jInv objective function defined by pMis and pInv
"""
function jInvObj(mc::BitArray{1},pMis,pInv)
   return jInvObj(Array{Float64}(mc),pMis,pInv);
end
function jInvObj(mc::Array{Float64},pMis,pInv)

    sig,dsig              = pInv.modelfun(mc)
	to = @elapsed begin
		Dc,F,dF,d2F,pMis,tMis = computeMisfit(sig,pMis,true)
	end
	timeF = to
	dF                    = dsig'*dF

	# Compute regularizer
	to = @elapsed begin
	R,dR,d2R = computeRegularizer(pInv.regularizer,mc,pInv.mref,pInv.MInv,pInv.alpha)
	end
	timeR = to

	# Objective function
	Jc = F  + R
	gc = dF + dR

    return Jc,gc,Dc,F,R,timeF,timeR
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
using DiffEqSensitivity
using DifferentialEquations
using ForwardDiff
using DiffResults
using Zygote
using Plots
using FiniteDifferences

# define function for RHS of lorenz system
function lorenz!(du, u, p, t)
    x,y,z = u
    σ,ρ,β = p

    du[1] = σ*(y - x)
    du[2] = x*(ρ - z) - y
    du[3] = x*y - β*z
end

# set up initial conditions for truth
u0 = [1.0; 0.0; 0.0]
p = [10.0, 28.0, (8/3)]
# tspan = (0.0, 100.0)
tspan = (0.0, 0.5)

# set up the ODE problem
prob = ODEProblem(lorenz!, u0, tspan, p)

# solve and plot
sol = solve(prob, Tsit5(), reltol=1e-6,abstol=1e-6)

# integrate model from k to k+1
function model_forward(uₖ)
    _prob = remake(prob, u0=uₖ)
    sol = solve(_prob, Tsit5(), reltol=1e-6,abstol=1e-6)
    return sol[:, end]
end

# # test that it works
J_forwarddiff = ForwardDiff.jacobian(model_forward, u0)

u_test = [1.0, 2.0, 2.0]

res = DiffResults.JacobianResult(u_test)
ForwardDiff.jacobian!(res, model_forward, u_test)
DiffResults.value(res)
DiffResults.jacobian(res)

# ----- Now try with zygote and compare
J_reversediff = Zygote.jacobian(model_forward, u0)[1]

res = Zygote.withjacobian(model_forward, u_test)
res.val
res.grad[1]

# ---- Let's verify with ForwardDiff
J_finitediff = FiniteDifferences.jacobian(central_fdm(5,1), model_forward, u0)[1]

J_forwarddiff ≈ J_reversediff


text:
@testset "Interaction" begin

    @testset "PowerLaw" begin
        p1 = PowerLaw(1)
        p2 = PowerLaw(2)

        @test p2 == PowerLaw(2.0) == PowerLaw(2//1)

        geom = Box([2,2])
        points = sample_blockaded(geom, 2)
        distances = distance_matrix(geom, points)

        expected = (1 ./ distances)
        expected[1,1] = expected[2,2] = 0
        @test interaction_matrix(p1, distances) ≈ expected
        @test interaction_matrix(p1, geom, points) ≈ interaction_matrix(p1, distances)
        @test isapprox(interaction_matrix(p2, geom, points), interaction_matrix(p1, distances) .^ 2; atol=1e-10) # need a bit lower precision
    end

end


text:
######################################################################
# Mesh
#######################################################################

abstract type AbstractMesh{D} end

struct Mesh{D,T,V<:AbstractArray{SVector{D,T}}} <: AbstractMesh{D}   # D is dimension of parameter space
    vertices::V                         # Iterable vertex container with SVector{D,T} eltype
    adjmat::SparseMatrixCSC{Bool,Int}   # Undirected graph: both dest > src and dest < src
end

# const Mesh{D,T} = Mesh{D,T,Vector{SVector{D,T}},Vector{Tuple{Int,Vararg{Int,D}}}}

# Mesh{D,T}() where {D,T} = Mesh(SVector{D,T}[], sparse(Int[], Int[], Bool[]), NTuple{D+1,Int}[])

function Base.show(io::IO, mesh::Mesh{D}) where {D}
    i = get(io, :indent, "")
    print(io,
"$(i)Mesh{$D}: mesh of a $D-dimensional manifold
$i  Vertices   : $(nvertices(mesh))
$i  Edges      : $(nedges(mesh))")
end

nvertices(m::Mesh) = length(m.vertices)

nedges(m::Mesh) = div(nnz(m.adjmat), 2)

nsimplices(m::Mesh) = length(simplices(m))

vertices(m::Mesh) = m.vertices

edges(m::Mesh, src) = nzrange(m.adjmat, src)

edgedest(m::Mesh, edge) = rowvals(m.adjmat)[edge]

edgevertices(m::Mesh) =
    ((vsrc, m.vertices[edgedest(m, edge)]) for (i, vsrc) in enumerate(m.vertices) for edge in edges(m, i))

function minmax_edge_length(m::Mesh{D,T}) where {D,T<:Real}
    minlen2 = typemax(T)
    maxlen2 = zero(T)
    verts = vertices(m)
    for src in eachindex(verts), edge in edges(m, src)
        dest = edgedest(m, edge)
        dest > src || continue # Need only directed graph
        vec = verts[dest] - verts[src]
        norm2 = vec' * vec
        norm2 < minlen2 && (minlen2 = norm2)
        norm2 > maxlen2 && (maxlen2 = norm2)
    end
    return sqrt(minlen2), sqrt(maxlen2)
end

######################################################################
# Compute N-simplices (N = number of vertices)
######################################################################
function simplices(mesh::Mesh{D}, ::Val{N} = Val(D+1)) where {D,N}
    N > 0 || throw(ArgumentError("Need a positive number of vertices for simplices"))
    N == 1 && return Tuple.(1:nvertices(mesh))
    simps = NTuple{N,Int}[]
    buffer = (NTuple{N,Int}[], NTuple{N,Int}[], Int[])
    for src in eachindex(vertices(mesh))
        append!(simps, _simplices(buffer, mesh, src))
    end
    N > 2 && alignnormals!(simps, vertices(mesh))
    return simps
end

# Add (greater) neighbors to last vertex of partials that are also neighbors of scr, till N
function _simplices(buffer::Tuple{P,P,V}, mesh, src) where {N,P<:AbstractArray{<:NTuple{N}},V}
    partials, partials´, srcneighs = buffer
    resize!(srcneighs, 0)
    resize!(partials, 0)
    for edge in edges(mesh, src)
        srcneigh = edgedest(mesh, edge)
        srcneigh > src || continue # Directed graph, to avoid simplex duplicates
        push!(srcneighs, srcneigh)
        push!(partials, padright((src, srcneigh), 0, Val(N)))
    end
    for pass in 3:N
        resize!(partials´, 0)
        for partial in partials
            nextsrc = partial[pass - 1]
            for edge in edges(mesh, nextsrc)
                dest = edgedest(mesh, edge)
                dest > nextsrc || continue # If not directed, no need to check
                dest in srcneighs && push!(partials´, modifyat(partial, pass, dest))
            end
        end
        partials, partials´ = partials´, partials
    end
    return partials
end

modifyat(s::NTuple{N,T}, ind, el) where {N,T} = ntuple(i -> i === ind ? el : s[i], Val(N))

function alignnormals!(simplices, vertices)
    for (i, s) in enumerate(simplices)
        volume = elementvolume(vertices, s)
        volume < 0 && (simplices[i] = switchlast(s))
    end
    return simplices
end

# Project N-1 edges onto (N-1)-dimensional vectors to have a deterministic volume
elementvolume(verts, s::NTuple{N,Int}) where {N} =
    elementvolume(hcat(ntuple(i -> padright(SVector(verts[s[i+1]] - verts[s[1]]), Val(N-1)), Val(N-1))...))
elementvolume(mat::SMatrix{N,N}) where {N} = det(mat)

switchlast(s::NTuple{N,T}) where {N,T} = ntuple(i -> i < N - 1 ? s[i] : s[2N - i - 1] , Val(N))

######################################################################
# Special meshes
######################################################################
"""
    marchingmesh(npoints::Integer...; axes = 1.0 * I, shift = missing)

Creates a L-dimensional marching-tetrahedra `Mesh`. The mesh is confined to the box defined
by the rows of `axes`, shifted by `shift` (an `NTuple` or `Number`) if not `missing`, and
contains `npoints[i]` along each axis `i`.

    marchingmesh(ranges::AbstractRange...; axes = 1.0 * I, shift = missing)

The same as above, but allows to specify the points in axis `i` by a range `ranges[i]`, such
as e.g. `-1.0:0.1:1.0`.

Note that the size of `axes` should match the number `L` of elements in `npoints` or
`ranges`. The `eltype` of points is given by that of `ranges` or `axes`.

    marchingmesh(h::Hamiltonian{<:Lattice}; npoints = 13, shift = missing)

Equivalent to `marchingmesh(ntuple(_ -> range(-π, π; length = npoints), Val(L))...)` where
`L` is the dimension of the Hamiltonian's lattice.

# External links

- Marching tetrahedra (https://en.wikipedia.org/wiki/Marching_tetrahedra) in Wikipedia
"""
marchingmesh(ranges::Vararg{AbstractRange,L}; axes = 1.0 * I, shift = missing) where {L} =
    _marchingmesh(shiftranges(ranges, shift), SMatrix{L,L}(axes))

function marchingmesh(npoints::Vararg{Integer,L}; axes = 1.0 * I, shift = missing) where {L}
    ranges = meshranges(npoints, Val(L))
    ranges´ = shiftranges(ranges, shift)
    return _marchingmesh(ranges´, SMatrix{L,L}(axes))
end

function marchingmesh(h::Hamiltonian{<:Lattice,L}; npoints = 13, shift = missing) where {L}
    ranges = meshranges(npoints, Val(L))
    ranges´ = shiftranges(ranges, shift)
    return _marchingmesh(ranges´, SMatrix{L,L}(I))
end

meshranges(n::Int, ::Val{L}) where {L} = meshranges(filltuple(n, Val(L)), Val(L))
meshranges(ns::NTuple{L,Int}, ::Val{L}) where {L} = (n -> range(-π, π; length = n)).(ns)

shiftranges(rs, shift::Missing) = rs
shiftranges(rs::NTuple{L}, shift::Number) where {L} = shiftranges(rs, filltuple(shift, Val(L)))
shiftranges(rs::NTuple{L}, shifts::NTuple{L}) where {L} = ((r, s) -> r .+ s).(rs, shifts)

function _marchingmesh(ranges::NTuple{D,AbstractRange}, axes::SMatrix{D,D}) where {D,T<:AbstractFloat}
    npoints = length.(ranges)
    projection = axes' # ./ (SVector(npoints) - 1)' # Projects binary vector to m box with npoints
    cs = CartesianIndices(ntuple(n -> 1:npoints[n], Val(D)))
    ls = LinearIndices(cs)
    csinner = CartesianIndices(ntuple(n -> 1:npoints[n]-1, Val(D)))

    # edge vectors for marching tetrahedra in D-dimensions (skip zero vector [first])
    uedges = [c for c in CartesianIndices(ntuple(_ -> 0:1, Val(D)))][2:end]
    # tetrahedra built from the D unit-length uvecs added in any permutation
    perms = permutations(
            ntuple(i -> CartesianIndex(ntuple(j -> i == j ? 1 : 0, Val(D))), Val(D)))
    utets = [cumsum(pushfirst!(perm, zero(CartesianIndex{D}))) for perm in perms]

    # We don't use generators because their non-inferreble eltype causes problems later
    verts = [projection * SVector(getindex.(ranges, Tuple(c))) for c in cs]

    s = SparseMatrixBuilder{Bool}(length(cs), length(cs))
    for c in cs
        for u in uedges
            dest = c + u    # dest > src
            dest in cs && pushtocolumn!(s, ls[dest], true)
            dest = c - u    # dest < src
            dest in cs && pushtocolumn!(s, ls[dest], true)
        end
        finalizecolumn!(s)
    end
    adjmat = sparse(s)

    return Mesh(verts, adjmat)
end