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

name = "at_spi2_core"
version = v"2.34.0"

# Collection of sources required to build at-spi2-core
sources = [
    ArchiveSource("http://ftp.gnome.org/pub/gnome/sources/at-spi2-core/$(version.major).$(version.minor)/at-spi2-core-$(version).tar.xz",
                  "d629cdbd674e539f8912028512af583990938c7b49e25184c126b00121ef11c6"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/at-spi2-core-*/
mkdir build && cd build

# Get a local gettext for msgfmt cross-building
apk add gettext

meson .. --cross-file="${MESON_TARGET_TOOLCHAIN}" \
    -Dintrospection=no \
    -Dx11=yes \
    -Dsystemd_user_dir=no
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if Sys.islinux(p) || Sys.isfreebsd(p)]

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Dbus_jll"),
    Dependency("Glib_jll", v"2.59.0"; compat="2.59.0"),
    Dependency("Xorg_libXtst_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using JHmat
using Test

@testset "JHmat.jl" begin

    @test greet()==1
end

ext: jl


text:
struct HydroDispatchRunOfRiverOutage <: PSI.AbstractHydroDispatchFormulation end
struct HydroDispatchReservoirStorageOutage <: PSI.AbstractHydroReservoirFormulation end
struct HydroDispatchPumpedStoragewReservationOutage <: PSI.AbstractHydroReservoirFormulation end

function outage_constraints!(
    optimization_container::PSI.OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::PSI.DeviceModel{T, D},
    ::Type{S},
    feedforward::Union{Nothing, PSI.AbstractAffectFeedForward},
) where {T <: PSY.HydroGen, S <: PM.AbstractPowerModel, D <: PSI.AbstractHydroFormulation}
    parameters = PSI.model_has_parameters(optimization_container)
    resolution = PSI.model_resolution(optimization_container)

    forecast_label = "outage"
    constraint_infos = Vector{DeviceOutageConstraintInfo}()
    for d in devices
        name = PSY.get_name(d)
        info = DeviceOutageConstraintInfo(
            name,
            nothing,
            1.0,
            PSI.get_time_series(optimization_container, d, forecast_label),
        )
        push!(constraint_infos, info)
    end

    if !(isempty(devices))
        if parameters
            device_outage_ub_parameter!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(OUTAGE, T),
                PSI.make_variable_name(PSI.ActivePowerVariable, T),
                PSI.UpdateRef{T}(OUTAGE, forecast_label),
            )
        else
            device_outage_ub!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(OUTAGE, T),
                PSI.make_variable_name(PSI.ActivePowerVariable, T),
            )
        end
    else
        @warn "Data doesn't contain generators of type $T, consider adjusting your formulation"
    end

    return
end

function outage_constraints!(
    optimization_container::PSI.OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::PSI.DeviceModel{T, D},
    ::Type{S},
    feedforward::Union{Nothing, PSI.AbstractAffectFeedForward},
) where {
    T <: PSY.HydroPumpedStorage,
    S <: PM.AbstractPowerModel,
    D <: PSI.AbstractHydroFormulation,
}
    parameters = PSI.model_has_parameters(optimization_container)
    resolution = PSI.model_resolution(optimization_container)

    forecast_label = "outage"
    constraint_infos = Vector{DeviceOutageConstraintInfo}()
    for d in devices
        name = PSY.get_name(d)
        info = DeviceOutageConstraintInfo(
            name,
            nothing,
            1.0,
            PSI.get_time_series(optimization_container, d, forecast_label),
        )
        push!(constraint_infos, info)
    end

    if !(isempty(devices))
        if parameters
            device_outage_ub_parameter!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(OUTAGE, T),
                PSI.make_variable_name(PSI.ActivePowerInVariable, T),
                PSI.UpdateRef{T}(OUTAGE, forecast_label),
            )
            device_outage_ub_parameter!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(OUTAGE, T),
                PSI.make_variable_name(PSI.ActivePowerOutVariable, T),
                PSI.UpdateRef{T}(OUTAGE, forecast_label),
            )
        else
            device_outage_ub!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(OUTAGE, T),
                PSI.make_variable_name(PSI.ActivePowerInVariable, T),
            )
            device_outage_ub!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(OUTAGE, T),
                PSI.make_variable_name(PSI.ActivePowerOutVariable, T),
            )
        end
    else
        @warn "Data doesn't contain generators of type $T, consider adjusting your formulation"
    end

    return
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

@doc raw"""
    IdealGlmMhdMulticomponentEquations1D

The ideal compressible multicomponent GLM-MHD equations in one space dimension.
"""
mutable struct IdealGlmMhdMulticomponentEquations1D{NVARS, NCOMP, RealT<:Real} <: AbstractIdealGlmMhdMulticomponentEquations{1, NVARS, NCOMP}
  gammas            ::SVector{NCOMP, RealT}
  gas_constants     ::SVector{NCOMP, RealT}
  cv                ::SVector{NCOMP, RealT}
  cp                ::SVector{NCOMP, RealT}

  function IdealGlmMhdMulticomponentEquations1D{NVARS, NCOMP, RealT}(gammas       ::SVector{NCOMP, RealT},
                                                                     gas_constants::SVector{NCOMP, RealT}) where {NVARS, NCOMP, RealT<:Real}

    NCOMP >= 1 || throw(DimensionMismatch("`gammas` and `gas_constants` have to be filled with at least one value"))

    cv = gas_constants ./ (gammas .- 1)
    cp = gas_constants + gas_constants ./ (gammas .- 1)

    new(gammas, gas_constants, cv, cp)
  end
end

function IdealGlmMhdMulticomponentEquations1D(; gammas, gas_constants)

  _gammas        = promote(gammas...)
  _gas_constants = promote(gas_constants...)
  RealT          = promote_type(eltype(_gammas), eltype(_gas_constants))

  NVARS = length(_gammas) + 7
  NCOMP = length(_gammas)

  __gammas        = SVector(map(RealT, _gammas))
  __gas_constants = SVector(map(RealT, _gas_constants))

  return IdealGlmMhdMulticomponentEquations1D{NVARS, NCOMP, RealT}(__gammas, __gas_constants)
end

@inline Base.real(::IdealGlmMhdMulticomponentEquations1D{NVARS, NCOMP, RealT}) where {NVARS, NCOMP, RealT} = RealT

have_nonconservative_terms(::IdealGlmMhdMulticomponentEquations1D) = Val(false)

function varnames(::typeof(cons2cons), equations::IdealGlmMhdMulticomponentEquations1D)

  cons  = ("rho_v1", "rho_v2", "rho_v3", "rho_e", "B1", "B2", "B3")
  rhos  = ntuple(n -> "rho" * string(n), Val(ncomponents(equations)))
  return (cons..., rhos...)
end

function varnames(::typeof(cons2prim), equations::IdealGlmMhdMulticomponentEquations1D)

  prim  = ("v1", "v2", "v3", "p", "B1", "B2", "B3")
  rhos  = ntuple(n -> "rho" * string(n), Val(ncomponents(equations)))
  return (prim..., rhos...)
end

"""
    initial_condition_convergence_test(x, t, equations::IdealGlmMhdMulticomponentEquations1D)

An Alfvén wave as smooth initial condition used for convergence tests.
"""
function initial_condition_convergence_test(x, t, equations::IdealGlmMhdMulticomponentEquations1D)
  # smooth Alfvén wave test from Derigs et al. FLASH (2016)
  # domain must be set to [0, 1], γ = 5/3

  rho = 1.0
  prim_rho  = SVector{ncomponents(equations), real(equations)}(2^(i-1) * (1-2)/(1-2^ncomponents(equations)) * rho for i in eachcomponent(equations))
  v1 = 0.0
  si, co = sincos(2 * pi * x[1])
  v2 = 0.1 * si
  v3 = 0.1 * co
  p = 0.1
  B1 = 1.0
  B2 = v2
  B3 = v3
  prim_other = SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3)
  return prim2cons(vcat(prim_other, prim_rho), equations)
end

"""
    initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdMulticomponentEquations1D)

Compound shock tube test case for one dimensional ideal MHD equations. It is bascially an
MHD extension of the Sod shock tube. Taken from Section V of the article
- Brio and Wu (1988)
  An Upwind Differencing Scheme for the Equations of Ideal Magnetohydrodynamics
  [DOI: 10.1016/0021-9991(88)90120-9](https://doi.org/10.1016/0021-9991(88)90120-9)
"""
function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdMulticomponentEquations1D)
  # domain must be set to [0, 1], γ = 2, final time = 0.12
  if x[1] < 0.5
    rho       = 1.0
    prim_rho  = SVector{ncomponents(equations), real(equations)}(2^(i-1) * (1-2)/(1-2^ncomponents(equations)) * rho for i in eachcomponent(equations))
  else
    rho       = 0.125
    prim_rho  = SVector{ncomponents(equations), real(equations)}(2^(i-1) * (1-2)/(1-2^ncomponents(equations)) * rho for i in eachcomponent(equations))
  end
  v1 = 0.0
  v2 = 0.0
  v3 = 0.0
  p = x[1] < 0.5 ? 1.0 : 0.1
  B1 = 0.75
  B2 = x[1] < 0.5 ? 1.0 : -1.0
  B3 = 0.0

  prim_other = SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3)
  return prim2cons(vcat(prim_other, prim_rho), equations)
end

# Calculate 1D flux in for a single point
@inline function flux(u, orientation::Integer, equations::IdealGlmMhdMulticomponentEquations1D)
  rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u

  rho = density(u, equations)

  v1 = rho_v1/rho
  v2 = rho_v2/rho
  v3 = rho_v3/rho
  kin_en = 0.5 * rho * (v1^2 + v2^2 + v3^2)
  mag_en = 0.5*(B1^2 + B2^2 + B3^2)
  gamma = totalgamma(u, equations)
  p = (gamma - 1) * (rho_e - kin_en - mag_en)

  f_rho = densities(u, v1, equations)
  f1 = rho_v1*v1 + p + mag_en - B1^2
  f2 = rho_v1*v2 - B1*B2
  f3 = rho_v1*v3 - B1*B3
  f4 = (kin_en + gamma*p/(gamma - 1) + 2*mag_en)*v1 - B1*(v1*B1 + v2*B2 + v3*B3)
  f5 = 0.0
  f6 = v1*B2 - v2*B1
  f7 = v1*B3 - v3*B1

  f_other  = SVector{7, real(equations)}(f1, f2, f3, f4, f5, f6, f7)

  return vcat(f_other, f_rho)
end

"""
    flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D)

Entropy conserving two-point flux adapted by
- Derigs et al. (2018)
  Ideal GLM-MHD: About the entropy consistent nine-wave magnetic field
  divergence diminishing ideal magnetohydrodynamics equations for multicomponent
  [DOI: 10.1016/j.jcp.2018.03.002](https://doi.org/10.1016/j.jcp.2018.03.002)
"""
function flux_derigs_etal(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdMulticomponentEquations1D)
  # Unpack left and right states to get velocities, pressure, and inverse temperature (called beta)
  rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll = u_ll
  rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr = u_rr
  @unpack gammas, gas_constants, cv = equations

  rho_ll = density(u_ll, equations)
  rho_rr = density(u_rr, equations)

  gamma_ll = totalgamma(u_ll, equations)
  gamma_rr = totalgamma(u_rr, equations)

  rhok_mean   = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i+7], u_rr[i+7]) for i in eachcomponent(equations))
  rhok_avg    = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i+7] + u_rr[i+7]) for i in eachcomponent(equations))

  v1_ll = rho_v1_ll/rho_ll
  v2_ll = rho_v2_ll/rho_ll
  v3_ll = rho_v3_ll/rho_ll
  v1_rr = rho_v1_rr/rho_rr
  v2_rr = rho_v2_rr/rho_rr
  v3_rr = rho_v3_rr/rho_rr
  vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
  vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
  mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
  mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
  # for convenience store v⋅B
  vel_dot_mag_ll = v1_ll*B1_ll + v2_ll*B2_ll + v3_ll*B3_ll
  vel_dot_mag_rr = v1_rr*B1_rr + v2_rr*B2_rr + v3_rr*B3_rr

  # Compute the necessary mean values needed for either direction
  v1_avg = 0.5*(v1_ll+v1_rr)
  v2_avg = 0.5*(v2_ll+v2_rr)
  v3_avg = 0.5*(v3_ll+v3_rr)
  v_sum  = v1_avg + v2_avg + v3_avg
  B1_avg = 0.5*(B1_ll+B1_rr)
  B2_avg = 0.5*(B2_ll+B2_rr)
  B3_avg = 0.5*(B3_ll+B3_rr)
  vel_norm_avg = 0.5*(vel_norm_ll+vel_norm_rr)
  mag_norm_avg = 0.5*(mag_norm_ll+mag_norm_rr)
  vel_dot_mag_avg = 0.5*(vel_dot_mag_ll+vel_dot_mag_rr)

  enth      = zero(v_sum)
  help1_ll  = zero(v1_ll)
  help1_rr  = zero(v1_rr)

  for i in eachcomponent(equations)
    enth      += rhok_avg[i] * gas_constants[i]
    help1_ll  += u_ll[i+7] * cv[i]
    help1_rr  += u_rr[i+7] * cv[i]
  end

  T_ll        = (rho_e_ll - 0.5*rho_ll * (vel_norm_ll) - 0.5*mag_norm_ll) / help1_ll
  T_rr        = (rho_e_rr - 0.5*rho_rr * (vel_norm_rr) - 0.5*mag_norm_rr) / help1_rr
  T           = 0.5 * (1.0/T_ll + 1.0/T_rr)
  T_log       = ln_mean(1.0/T_ll, 1.0/T_rr)

  # Calculate fluxes depending on orientation with specific direction averages
  help1       = zero(T_ll)
  help2       = zero(T_rr)

  f_rho       = SVector{ncomponents(equations), real(equations)}(rhok_mean[i]*v1_avg for i in eachcomponent(equations))
  for i in eachcomponent(equations)
    help1     += f_rho[i] * cv[i]
    help2     += f_rho[i]
  end
  f1 = help2 * v1_avg + enth/T + 0.5 * mag_norm_avg - B1_avg*B1_avg
  f2 = help2 * v2_avg - B1_avg*B2_avg
  f3 = help2 * v3_avg - B1_avg*B3_avg
  f5 = 0.0
  f6 = v1_avg*B2_avg - v2_avg*B1_avg
  f7 = v1_avg*B3_avg - v3_avg*B1_avg

  # total energy flux is complicated and involves the previous eight components
  v1_mag_avg = 0.5*(v1_ll*mag_norm_ll + v1_rr*mag_norm_rr)

  f4 = (help1/T_log) - 0.5 * (vel_norm_avg) * (help2) + f1 * v1_avg + f2 * v2_avg + f3 * v3_avg +
        f5 * B1_avg + f6 * B2_avg + f7 * B3_avg - 0.5*v1_mag_avg +
        B1_avg * vel_dot_mag_avg

  f_other  = SVector{7, real(equations)}(f1, f2, f3, f4, f5, f6, f7)

  return vcat(f_other, f_rho)
end

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdMulticomponentEquations1D)
  rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll = u_ll
  rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr = u_rr

  rho_ll   = density(u_ll, equations)
  rho_rr   = density(u_rr, equations)
  gamma_ll = totalgamma(u_ll, equations)
  gamma_rr = totalgamma(u_rr, equations)

  # Calculate velocities and fast magnetoacoustic wave speeds
  # left
  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v3_ll = rho_v3_ll / rho_ll
  v_mag_ll = sqrt(v1_ll^2 + v2_ll^2 + v3_ll^2)
  cf_ll = calc_fast_wavespeed(u_ll, orientation, equations)
  # right
  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  v3_rr = rho_v3_rr / rho_rr
  v_mag_rr = sqrt(v1_rr^2 + v2_rr^2 + v3_rr^2)
  cf_rr = calc_fast_wavespeed(u_rr, orientation, equations)

  λ_max = max(v_mag_ll, v_mag_rr) + max(cf_ll, cf_rr)
end

@inline function max_abs_speeds(u, equations::IdealGlmMhdMulticomponentEquations1D)
  rho_v1, _ = u

  rho = density(u, equations)

  v1 = rho_v1 / rho

  cf_x_direction = calc_fast_wavespeed(u, 1, equations)

  return (abs(v1) + cf_x_direction, )
end

# Convert conservative variables to primitive
function cons2prim(u, equations::IdealGlmMhdMulticomponentEquations1D)
  rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u

  prim_rho = SVector{ncomponents(equations), real(equations)}(u[i+7] for i in eachcomponent(equations))
  rho = density(u, equations)

  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  v3 = rho_v3 / rho

  gamma = totalgamma(u, equations)

  p = (gamma - 1) * (rho_e - 0.5*rho*(v1^2 + v2^2 + v3^2) - 0.5*(B1^2 + B2^2 + B3^2))
  prim_other =  SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3)
  return vcat(prim_other, prim_rho)
end

# Convert conservative variables to entropy
@inline function cons2entropy(u, equations::IdealGlmMhdMulticomponentEquations1D)
  rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
  @unpack cv, gammas, gas_constants = equations

  rho = density(u, equations)

  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  v3 = rho_v3 / rho
  v_square = v1^2 + v2^2 + v3^2
  gamma = totalgamma(u, equations)
  p = (gamma - 1) * (rho_e - 0.5*rho*v_square - 0.5*(B1^2 + B2^2 + B3^2))
  s = log(p) - gamma*log(rho)
  rho_p = rho / p

  # Multicomponent stuff
  help1 = zero(v1)

  for i in eachcomponent(equations)
    help1 += u[i+7] * cv[i]
  end

  T         = (rho_e - 0.5 * rho * v_square - 0.5*(B1^2 + B2^2 + B3^2)) / (help1)

  entrop_rho  = SVector{ncomponents(equations), real(equations)}( -1.0 * (cv[i] * log(T) - gas_constants[i] * log(u[i+7])) + gas_constants[i] + cv[i] - (v_square / (2*T)) for i in eachcomponent(equations))

  w1 = v1 / T
  w2 = v2 / T
  w3 = v3 / T
  w4 = -1.0 / T
  w5 = B1 / T
  w6 = B2 / T
  w7 = B3 / T

  entrop_other = SVector{7, real(equations)}(w1, w2, w3, w4, w5, w6, w7)

  return vcat(entrop_other, entrop_rho)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equations::IdealGlmMhdMulticomponentEquations1D)
  v1, v2, v3, p, B1, B2, B3 = prim

  cons_rho = SVector{ncomponents(equations), real(equations)}(prim[i+7] for i in eachcomponent(equations))
  rho = density(prim, equations)

  rho_v1 = rho * v1
  rho_v2 = rho * v2
  rho_v3 = rho * v3

  gamma = totalgamma(prim, equations)
  rho_e = p/(gamma-1) + 0.5 * (rho_v1*v1 + rho_v2*v2 + rho_v3*v3) +
                                 0.5 * (B1^2 + B2^2 + B3^2)

  cons_other = SVector{7, real(equations)}(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3)

  return vcat(cons_other, cons_rho)
end

@inline function density_pressure(u, equations::IdealGlmMhdMulticomponentEquations1D)
  rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
  rho = density(u, equations)
  gamma = totalgamma(u, equations)
  p = (gamma - 1)*(rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
                                   - 0.5 * (B1^2 + B2^2 + B3^2)
                                   )
  return rho * p
end

# Compute the fastest wave speed for ideal MHD equations: c_f, the fast magnetoacoustic eigenvalue
@inline function calc_fast_wavespeed(cons, direction, equations::IdealGlmMhdMulticomponentEquations1D)
  rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = cons
  rho = density(cons, equations)
  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  v3 = rho_v3 / rho
  v_mag = sqrt(v1^2 + v2^2 + v3^2)
  gamma = totalgamma(cons, equations)
  p = (gamma - 1)*(rho_e - 0.5*rho*v_mag^2 - 0.5*(B1^2 + B2^2 + B3^2))
  a_square = gamma * p / rho
  sqrt_rho = sqrt(rho)
  b1 = B1 / sqrt_rho
  b2 = B2 / sqrt_rho
  b3 = B3 / sqrt_rho
  b_square = b1^2 + b2^2 + b3^2

  c_f = sqrt(0.5*(a_square + b_square) + 0.5*sqrt((a_square + b_square)^2 - 4.0*a_square*b1^2))

  return c_f
end

@inline function density(u, equations::IdealGlmMhdMulticomponentEquations1D)
  rho = zero(u[1])

  for i in eachcomponent(equations)
    rho += u[i+7]
  end

  return rho
 end

 @inline function totalgamma(u, equations::IdealGlmMhdMulticomponentEquations1D)
  @unpack cv, gammas = equations

  help1 = zero(u[1])
  help2 = zero(u[1])

  for i in eachcomponent(equations)
    help1 += u[i+7] * cv[i] * gammas[i]
    help2 += u[i+7] * cv[i]
  end

  return help1/help2
end

@inline function densities(u, v, equations::IdealGlmMhdMulticomponentEquations1D)

  return SVector{ncomponents(equations), real(equations)}(u[i+7]*v for i in eachcomponent(equations))
 end

end # @muladd


text:
@testset "chempotential.jl" begin
    C₁C₅ = CubicEoS.load(BrusilovskyEoSMixture; names=("methane", "n-pentane"))
    N = 1e-2 .+ rand(2)
    V = 1e-2 + rand()
    RT = 300 * CubicEoS.GAS_CONSTANT_SI

    @testset "func vs func! versions" begin
        vec = similar(N)
        log_c_activity!(vec, C₁C₅, N, V, RT)
        @test vec == log_c_activity(C₁C₅, N, V, RT)

        matr = vec * vec'
        log_c_activity_wj!(vec, matr, C₁C₅, N, V, RT)
        @test (vec, matr) == log_c_activity_wj(C₁C₅, N, V, RT)
    end

    @testset "Buffered versions" begin
        buf = thermo_buffer(C₁C₅)
        @test log_c_activity(C₁C₅, N, V, RT) == log_c_activity(C₁C₅, N, V, RT; buf=buf)
        @test log_c_activity_wj(C₁C₅, N, V, RT) == log_c_activity_wj(C₁C₅, N, V, RT; buf=buf)
    end

    @testset "Homogeneity" begin
        for α in 0.1:0.2:1.9
            @test log_c_activity(C₁C₅, N, V, RT) ≈ log_c_activity(C₁C₅, α .* N, α * V, RT)

            v, m = log_c_activity_wj(C₁C₅, N, V, RT)
            vα, mα = log_c_activity_wj(C₁C₅, α .* N, α * V, RT)
            @test v ≈ vα
            # m is jacobian
            # while activity has 0-order homogeneity,
            # its derivatives (jacobian) is -1-order homogeneous
            @test m ≈ α .* mα
        end
    end
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
using Plots

"Enumerate possible states of a single plant"
@enum InfectionStatus uninfected infected dead recovered immune

"Data structure containing the infection status of a plant"
mutable struct Plant
    status::InfectionStatus
    infection_time::Int8
end

"Parameters for a simulation"
mutable struct Parameters
    infection_rate::Float32
    death_probability::Float32
    recovery_time::Int8
    immunity_rate::Float32
end

"Create a map of immune and uninfected plants with 1 infected in the middle."
function make_plants(width::Integer=7, height::Integer=7, immunity_rate::Float64=0.0)
    plants = Matrix{Plant}(undef, width, height)
    for i in 1:size(plants)[1]
        for j in 1:size(plants)[2]
            if rand(1)[1] < immunity_rate
                plants[i,j] = Plant(immune, 0)
            else
                plants[i,j] = Plant(uninfected, 0)
            end
        end
    end
    plants[width÷2+1,height÷2+1] = Plant(infected, 0)
    return plants
end

"Show the plant as an empty or a filled in box for quick viewing"
function Base.show(io::IO, plant::Plant)
    if plant.status == infected
        print(io, "◼")
    else
        print(io, "◻")
    end
end

"Map the plants to colors for plotting"
function to_colors(plant::Plant)
    if plant.status == uninfected
        return RGB(0.0,0.8,0.0)
    end
    if plant.status == infected
        return RGB(0.8,0.0,0.0)
    end
    if plant.status == dead
        return RGB(0.1,0.1,0.1)
    end
    if plant.status == recovered
        return RGB(0.0,0.0,0.8)
    end
    if plant.status == immune
        return RGB(0.2,0.2,1.0)
    end
end

"""
Simulate the interaction between one plant and a neighbour.

If the neighbour is infected, it infect this plant with the propability parameters.infection_rate.
"""
function interact!(new_plant::Plant, other_plant::Plant, infection_rate)
    if new_plant.status == uninfected && other_plant.status == infected
        if rand(1)[1] < infection_rate
            new_plant.status = infected
            new_plant.infection_time = 0
        end
    end
end

"""
Update a single plant, not accounting for it's interactions with the neighbours.
"""
function update!(new_plant::Plant, recovery_time, death_rate)
    if new_plant.status == infected
        # Increase the number of time steps since infection
        new_plant.infection_time += 1

        # Check if the plant recovers
        if new_plant.infection_time > recovery_time
            new_plant.status = recovered
        end

        # Check if it dies
        if rand(1)[1] < death_rate
            new_plant.status = dead
        end
    end
end

"""
Update the plants in the 2D array of plants, using given parameters.
"""
function update!(plants::Matrix{Plant}, recovery_time, death_rate, infection_rate)
    old_plants = deepcopy(plants)
    for i in 1:size(plants)[1]
        for j in 1:size(plants)[2]
            update!(plants[i,j], recovery_time, death_rate)
        end
    end
    for i in 1:size(plants)[1]-1
        for j in 1:size(plants)[2]
            interact!(plants[i,j], old_plants[i+1,j], infection_rate)
            interact!(plants[i+1,j], old_plants[i,j], infection_rate)
        end
    end
    for i in 1:size(plants)[1]
        for j in 1:size(plants)[2]-1
            interact!(plants[i,j], old_plants[i,j+1], infection_rate)
            interact!(plants[i,j+1], old_plants[i,j], infection_rate)
        end
    end
end

"Count the current number of infections"
function count_infections(plants::Matrix{Plant})
    infections = 0
    for i in 1:size(plants)[1]
        for j in 1:size(plants)[2]
            if plants[i,j].status == infected
                infections += 1
            end
        end
    end
    return infections
end

"Count the number of dead plants"
function count_deaths(plants::Matrix{Plant})
    deaths = 0
    for i in 1:size(plants)[1]
        for j in 1:size(plants)[2]
            if plants[i,j].status == dead
                deaths += 1
            end
        end
    end
    return deaths
end

max_stars_repo_name: nnhjy/julia-introduction

max_issues_repo_name: nnhjy/julia-introduction

max_forks_repo_name: nnhjy/julia-introduction


text:
"This script computes the Laguerre-Kravchuk functions, who are the Kravchuk functions on a angular momentum space."

using PixelsGT # Call the parent module

export  ϕd, ϕu, LKd, LKu, LK

"`Φd(q₁, q₂, n, m, j)` rotates π/2 the Kravchuk function modes at (n, |m| ≤ n) for n ∈ [0, 2j]. It needs to be evaluted at (q₁, q₂) ∈ [0, 2j + 1]."
function Φd(q₁, q₂, n, m, j)
    v₀ = zeros(Integer(2 * j + 1)) # (Sep 9, 2021) Added Integer() to meet the half-spin feature in j
    for n₁ in 0:n
        v₀[Integer(n₁ + 1)] = Φ2(n₁, n - n₁, q₁, q₂, j) * dd(n, j)[n₁ + 1, Integer((m + n)/2) + 1]
    end
    return sum(v₀)
end

"`Φu(q₁, q₂, n, m, j)` rotates π/2 the Kravchuk function modes at (n, |m| ≤ 4j - n) for n ∈ [4j, 2j]. It needs to be evaluted at (q₁, q₂) ∈ [0, 2j + 1]."
function Φu(q₁, q₂, n, m, j)
    v₁ = zeros(Integer(2 * j + 1))
    for n₁ in 0:Integer((4* j - n))
        v₁[n₁ + 1] = Φ2(2 * j - n₁, n + n₁ - 2 * j, q₁, q₂, j) * du(n, j)[n₁ + 1, Integer((m + 4 * j - n)/2) + 1]
    end
    return sum(v₁)
end

"Laguerre-Kravchuk modes functions are the core of the rotation of images, they encodes the original position information in a mode-angular momentum space, ready to apply an angle phase who rotates the image"

"`LKd(q₁, q₂, n, m, j)` gives Laguerre-Kravchuk modes at (n, |m| ≤ n) for n ∈ [0, 2j]. It needs to be evaluted at (q₁, q₂) ∈ [0, 2j + 1]."
function LKd(q₁, q₂, n, m, j)
    vlkd = zeros(ComplexF64, Integer(2 * j + 1))
    for n₁ in 0:n
        vlkd[n₁ + 1] = (-1.0im)^(n - n₁) * Φ2(n₁, n - n₁, q₁, q₂, j) * dd(n, j)[n₁ + 1, Integer((m + n)/2) + 1]
    end
    return (-1.0)^((0.5) * (abs(m) - m)) * sum(vlkd)
end

"`LKu(q₁, q₂, n, m, j)` gives the Laguerre-Kravchuk modes at (n, |m| ≤ 4j - n) for n ∈ [4j, 2j]. It needs to be evaluted at (q₁, q₂) ∈ [0, 2j + 1]."
function LKu(q₁, q₂, n, m, j)
    vlku = zeros(ComplexF64, Integer(2 * j + 1))
    for n₁ in 0:Integer(4* j - n)
        vlku[n₁ + 1] = (-1.0im)^(n - n₁) * Φ2(2 * j - n₁, n + n₁ - 2 * j, q₁, q₂, j) * du(n, j)[n₁ + 1, Integer((m + 4 * j - n)/2) + 1]
    end
    return (-1.0)^((0.5) * (abs(m) - m)) * sum(vlku)
end

"`LK(q₁, q₂, n, m, j)` gives the complete Laguerre-Kravchuk modes. Their index ranges are n ∈ [0, 4j] and m ∈ [-n,n] in steps of 2."
function LK(q₁, q₂, n, m, j)
    if 0 ≤ n ≤ 2 * j
        return LKd(q₁, q₂, n, m, j)
    elseif 2 * j < n ≤ 4 * j
        return LKu(q₁, q₂, n, m, j)
    end
end


text:
abstract type AbstractMonotonizedSymmetricConnectivityModel end

@inline (mdl::AbstractMonotonizedSymmetricConnectivityModel)(n_vertices, n_over_pins, n_dia_nets, k) = mdl(n_vertices, n_over_pins, n_dia_nets)

struct AffineMonotonizedSymmetricConnectivityModel{Tv} <: AbstractMonotonizedSymmetricConnectivityModel
    α::Tv
    β_vertex::Tv
    β_over_pin::Tv
    β_dia_net::Tv
    Δ_pins::Tv
end

function AffineMonotonizedSymmetricConnectivityModel(; α = false, β_vertex = false, β_over_pin = false, β_dia_net = false, Δ_pins = false)
    AffineMonotonizedSymmetricConnectivityModel(promote(α, β_vertex, β_over_pin, β_dia_net, Δ_pins)...)
end

function AffineMonotonizedSymmetricConnectivityModel(mdl::AffineSymmetricConnectivityModel{Tv}) where {Tv}
    α = mdl.α
    β_dia_net = mdl.β_remote_net
    β_over_pin = mdl.β_pin
    if mdl.β_vertex < mdl.β_remote_net
        Δ_pins = cld(mdl.β_remote_net - mdl.β_vertex, mdl.β_pin)
        β_vertex = zero(Tv)
    else
        Δ_pins = 0
        β_vertex = mdl.β_vertex - mdl.β_remote_net
    end
    return AffineMonotonizedSymmetricConnectivityModel(α, β_vertex, β_over_pin, β_dia_net, Δ_pins)
end

@inline cost_type(::Type{AffineMonotonizedSymmetricConnectivityModel{Tv}}) where {Tv} = Tv

(mdl::AffineMonotonizedSymmetricConnectivityModel)(n_vertices, n_over_pins, n_dia_nets) = mdl.α + n_vertices * mdl.β_vertex + n_over_pins * mdl.β_over_pin + n_dia_nets * mdl.β_dia_net

function bound_stripe(A::SparseMatrixCSC, K, ocl::AbstractOracleCost{<:AffineMonotonizedSymmetricConnectivityModel})
    m, n = size(A)
    @assert m == n
    N = nnz(A)
    mdl = oracle_model(ocl)
    @assert mdl.β_vertex >= 0
    @assert mdl.β_over_pin >= 0
    @assert mdl.β_dia_net >= 0
    c_hi = ocl(1, n + 1)
    c_lo = mdl.α + fld(c_hi - mdl.α, K)
    return (c_lo, c_hi)
end

function bound_stripe(A::SparseMatrixCSC, K, mdl::AffineMonotonizedSymmetricConnectivityModel)
    @inbounds begin
        @assert mdl.β_vertex >= 0
        @assert mdl.β_over_pin >= 0
        @assert mdl.β_dia_net >= 0
        m, n = size(A)
        @assert m == n
        N = nnz(A)
        n_over_pins = 0
        for j = 1:n
            n_over_pins += max(A.colptr[j + 1] - A.colptr[j] - mdl.Δ_pins, 0)
        end
        c_hi = mdl.α + mdl.β_vertex * n + mdl.β_over_pin * n_over_pins + mdl.β_dia_net * m
        c_lo = mdl.α + fld(c_hi - mdl.α, K)
        return (c_lo, c_hi)
    end
end

struct MonotonizedSymmetricConnectivityOracle{Ti, DiaNet, Mdl} <: AbstractOracleCost{Mdl}
    overpos::Vector{Ti}
    dianet::DiaNet
    mdl::Mdl
end

oracle_model(ocl::MonotonizedSymmetricConnectivityOracle) = ocl.mdl

function oracle_stripe(hint::AbstractHint, mdl::AbstractMonotonizedSymmetricConnectivityModel, A::SparseMatrixCSC{Tv, Ti}; net=nothing, adj_A=nothing, kwargs...) where {Tv, Ti}
    @inbounds begin
        m, n = size(A)
        @assert m == n
        N = nnz(A)
        overpos = undefs(eltype(A.colptr), n + 1)
        overpos[1] = 1
        for j = 1:n
            overpos[j + 1] = overpos[j] + max(A.colptr[j + 1] - A.colptr[j] - mdl.Δ_pins, 0)
        end

        dianet = dianetcount(hint, A; kwargs...)

        return MonotonizedSymmetricConnectivityOracle(overpos, dianet, mdl)
    end
end

function bound_stripe(A::SparseMatrixCSC, K, ocl::MonotonizedSymmetricConnectivityOracle{<:Any, <:Any, <:AffineMonotonizedSymmetricConnectivityModel})
    m, n = size(A)
    @assert m == n
    N = nnz(A)
    mdl = oracle_model(ocl)
    @assert mdl.β_vertex >= 0
    @assert mdl.β_over_pin >= 0
    @assert mdl.β_dia_net >= 0
    c_hi = mdl.α + mdl.β_vertex * n + mdl.β_over_pin * (ocl.overpos[end] - ocl.overpos[1]) + mdl.β_dia_net * m
    c_lo = mdl.α + fld(c_hi - mdl.α, K)
    return (c_lo, c_hi)
end

@inline function (cst::MonotonizedSymmetricConnectivityOracle{Ti, Mdl})(j::Ti, j′::Ti, k...) where {Ti, Mdl}
    @inbounds begin
        w = cst.overpos[j′] - cst.overpos[j]
        d = cst.dianet[j, j′]
        return cst.mdl(j′ - j, w, d, k...)
    end
end

@inline function (stp::Step{Ocl})(_j, _j′, _k...) where {Ti, Mdl, Ocl <: MonotonizedSymmetricConnectivityOracle{Ti, Mdl}}
    @inbounds begin
        cst = stp.ocl
        j = destep(_j)
        j′ = destep(_j′)
        k = maptuple(destep, _k...)
        w = cst.overpos[j′] - cst.overpos[j]
        d = Step(cst.dianet)(_j, _j′)
        return cst.mdl(j′ - j, w, d, k...)
    end
end


text:
import Base: ==, ≈, copy, eltype, rationalize, extrema
import Random.rand

export LazySet,
       basetype,
       ρ, support_function,
       σ, support_vector,
       complement,
       dim,
       norm,
       radius,
       diameter,
       an_element,
       isbounded,
       isboundedtype,
       neutral,
       absorbing,
       tosimplehrep,
       isuniversal,
       translate,
       affine_map,
       exponential_map,
       reflect,
       is_interior_point,
       isoperation,
       isoperationtype,
       isequivalent,
       isconvextype,
       is_polyhedral,
       area,
       surface,
       singleton_list,
       concretize,
       constraints,
       vertices,
       project,
       rectify,
       permute

"""
    LazySet{N}

Abstract type for convex sets, i.e., sets characterized by a (possibly infinite)
intersection of halfspaces, or equivalently, sets ``S`` such that for any two
elements ``x, y ∈ S`` and ``0 ≤ λ ≤ 1`` it holds that ``λ·x + (1-λ)·y ∈ S``.

### Notes

`LazySet` types should be parameterized with a type `N`, typically `N<:Real`,
for using different numeric types.

Every concrete `LazySet` must define the following functions:
- `σ(d::AbstractVector, S::LazySet)` -- the support vector of `S` in a given
    direction `d`
- `dim(S::LazySet)` -- the ambient dimension of `S`

The function
- `ρ(d::AbstractVector, S::LazySet)` -- the support function of `S` in a given
    direction `d`
is optional because there is a fallback implementation relying on `σ`.
However, for unbounded sets (which includes most lazy set types) this fallback
cannot be used and an explicit method must be implemented.

The subtypes of `LazySet` (including abstract interfaces):

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(LazySet, false)
14-element Vector{Any}:
 AbstractAffineMap
 AbstractCentrallySymmetric
 AbstractPolyhedron
 Bloating
 CachedMinkowskiSumArray
 CartesianProduct
 CartesianProductArray
 ConvexHull
 ConvexHullArray
 EmptySet
 Intersection
 IntersectionArray
 MinkowskiSum
 MinkowskiSumArray
```

If we only consider *concrete* subtypes, then:

```jldoctest; setup = :(using LazySets: subtypes)
julia> concrete_subtypes = subtypes(LazySet, true);

julia> length(concrete_subtypes)
44

julia> println.(concrete_subtypes);
AffineMap
Ball1
Ball2
BallInf
Ballp
Bloating
CachedMinkowskiSumArray
CartesianProduct
CartesianProductArray
ConvexHull
ConvexHullArray
Ellipsoid
EmptySet
ExponentialMap
ExponentialProjectionMap
HParallelotope
HPolygon
HPolygonOpt
HPolyhedron
HPolytope
HalfSpace
Hyperplane
Hyperrectangle
Intersection
IntersectionArray
Interval
InverseLinearMap
Line
Line2D
LineSegment
LinearMap
MinkowskiSum
MinkowskiSumArray
ResetMap
RotatedHyperrectangle
Singleton
Star
SymmetricIntervalHull
Translation
Universe
VPolygon
VPolytope
ZeroSet
Zonotope
```
"""
abstract type LazySet{N} end

# --- common LazySet functions ---

"""
    eltype(::Type{<:LazySet{N}}) where {N}

Return the numeric type (`N`) of the given set type.

### Input

- `T` -- set type, used for dispatch

### Output

The numeric type of `T`.
"""
eltype(::Type{<:LazySet{N}}) where {N} = N

"""
    eltype(::LazySet{N}) where {N}

Return the numeric type (`N`) of the given set.

### Input

- `X` -- set instance, used for dispatch

### Output

The numeric type of `X`.
"""
eltype(::LazySet{N}) where {N} = N

"""
    basetype(T::Type{<:LazySet})

Return the base type of the given set type (i.e., without type parameters).

### Input

- `T` -- set type, used for dispatch

### Output

The base type of `T`.
"""
basetype(T::Type{<:LazySet}) = Base.typename(T).wrapper

"""
    basetype(S::LazySet)

Return the base type of the given set (i.e., without type parameters).

### Input

- `S` -- set instance, used for dispatch

### Output

The base type of `S`.

### Examples

```jldoctest
julia> z = rand(Zonotope);

julia> basetype(z)
Zonotope

julia> basetype(z + z)
MinkowskiSum

julia> basetype(LinearMap(rand(2, 2), z + z))
LinearMap
```
"""
basetype(S::LazySet) = Base.typename(typeof(S)).wrapper

"""
    ρ(d::AbstractVector, S::LazySet)

Evaluate the support function of a set in a given direction.

### Input

- `d` -- direction
- `S` -- convex set

### Output

The support function of the set `S` for the direction `d`.
"""
function ρ(d::AbstractVector, S::LazySet)
    return dot(d, σ(d, S))
end

"""
    support_function

Alias for the support function ρ.
"""
const support_function = ρ

"""
    σ

Function to compute the support vector σ.
"""
function σ end

"""
    support_vector

Alias for the support vector σ.
"""
const support_vector = σ

"""
    isboundedtype(::Type{<:LazySet})

Determine whether a set type only represents bounded sets.

### Input

- `LazySet` -- set type for dispatch

### Output

`true` if the set type only represents bounded sets.
Note that some sets may still represent an unbounded set even though their type
actually does not (example: [`HPolytope`](@ref), because the construction with
non-bounding linear constraints is allowed).

### Notes

By default this function returns `false`.
All set types that can determine boundedness should override this behavior.
"""
function isboundedtype(::Type{T}) where {T<:LazySet}
    return false
end

"""
    isbounded(S::LazySet)

Determine whether a set is bounded.

### Input

- `S`         -- set
- `algorithm` -- (optional, default: `"support_function"`) algorithm choice,
                 possible options are `"support_function"` and `"stiemke"`

### Output

`true` iff the set is bounded.

### Algorithm

See the documentation of `_isbounded_unit_dimensions` or `_isbounded_stiemke`
for details.
"""
function isbounded(S::LazySet; algorithm="support_function")
    if algorithm == "support_function"
        return _isbounded_unit_dimensions(S)
    elseif algorithm == "stiemke"
        return _isbounded_stiemke(S)
    else
        throw(ArgumentError("unknown algorithm $algorithm"))
    end
end

"""
    _isbounded_unit_dimensions(S::LazySet{N}) where {N}

Determine whether a set is bounded in each unit dimension.

### Input

- `S` -- set

### Output

`true` iff the set is bounded in each unit dimension.

### Algorithm

This function performs ``2n`` support function checks, where ``n`` is the
ambient dimension of `S`.
"""
function _isbounded_unit_dimensions(S::LazySet{N}) where {N}
    n = dim(S)
    @inbounds for i in 1:n
        for o in [one(N), -one(N)]
            d = SingleEntryVector(i, n, o)
            if ρ(d, S) == N(Inf)
                return false
            end
        end
    end
    return true
end

"""
    is_polyhedral(S::LazySet)

Trait for polyhedral sets.

### Input

- `S` -- set

### Output

`true` only if the set behaves like an [`AbstractPolyhedron`](@ref).
The answer is conservative, i.e., may sometimes be `false` even if the set is
polyhedral.
"""
function is_polyhedral(S::LazySet)
    return false
end

"""
    norm(S::LazySet, [p]::Real=Inf)

Return the norm of a convex set.
It is the norm of the enclosing ball (of the given ``p``-norm) of minimal volume
that is centered in the origin.

### Input

- `S` -- convex set
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the norm.
"""
function norm(S::LazySet, p::Real=Inf)
    if p == Inf
        return norm(box_approximation(S), p)
    elseif applicable(vertices_list, S)
        return maximum(norm(v, p) for v in vertices_list(S))
    else
        error("the norm for this value of p=$p is not implemented")
    end
end

"""
    radius(S::LazySet, [p]::Real=Inf)

Return the radius of a convex set.
It is the radius of the enclosing ball (of the given ``p``-norm) of minimal
volume with the same center.

### Input

- `S` -- convex set
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the radius.
"""
function radius(S::LazySet, p::Real=Inf)
    if p == Inf
        return radius(Approximations.ballinf_approximation(S), p)
    else
        error("the radius for this value of p=$p is not implemented")
    end
end

"""
    diameter(S::LazySet, [p]::Real=Inf)

Return the diameter of a convex set.
It is the maximum distance between any two elements of the set, or,
equivalently, the diameter of the enclosing ball (of the given ``p``-norm) of
minimal volume with the same center.

### Input

- `S` -- convex set
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the diameter.
"""
function diameter(S::LazySet, p::Real=Inf)
    return radius(S, p) * 2
end

"""
    affine_map(M::AbstractMatrix, X::LazySet, v::AbstractVector; kwargs...)

Compute a concrete affine map.

### Input

- `M` -- linear map
- `X` -- convex set
- `v` -- translation vector

### Output

A set representing the affine map of `X`.

### Algorithm

The implementation applies the functions `linear_map` and `translate`.
"""
function affine_map(M::AbstractMatrix, X::LazySet, v::AbstractVector; kwargs...)
    return translate(linear_map(M, X; kwargs...), v)
end

"""
    exponential_map(M::AbstractMatrix, X::LazySet)

Compute the concrete exponential map of `M` and `X`, i.e., `exp(M) * X`.

### Input

- `M` -- matrix
- `X` -- set

### Output

A set representing the exponential map of `M` and `X`.

### Algorithm

The implementation applies the functions `exp` and `linear_map`.
"""
function exponential_map(M::AbstractMatrix, X::LazySet)
    return linear_map(exp(M), X)
end

"""
    an_element(S::LazySet{N}) where {N}

Return some element of a convex set.

### Input

- `S` -- convex set

### Output

An element of a convex set.

### Algorithm

An element of the set is obtained by evaluating its support vector along
direction ``[1, 0, …, 0]``.
"""
function an_element(S::LazySet{N}) where {N}
    e₁ = SingleEntryVector(1, dim(S), one(N))
    return σ(e₁, S)
end

"""
    ==(X::LazySet, Y::LazySet)

Return whether two LazySets of the same type are exactly equal.

### Input

- `X` -- any `LazySet`
- `Y` -- another `LazySet` of the same type as `X`

### Output

- `true` iff `X` is equal to `Y`.

### Notes

The check is purely syntactic and the sets need to have the same base type.
For instance, `X::VPolytope == Y::HPolytope` returns `false` even if `X` and `Y`
represent the same polytope.
However `X::HPolytope{Int64} == Y::HPolytope{Float64}` is a valid comparison.

### Algorithm

We recursively compare the fields of `X` and `Y` until a mismatch is found.

### Examples

```jldoctest
julia> HalfSpace([1], 1) == HalfSpace([1], 1)
true

julia> HalfSpace([1], 1) == HalfSpace([1.0], 1.0)
true

julia> Ball1([0.], 1.) == Ball2([0.], 1.)
false
```
"""
function ==(X::LazySet, Y::LazySet)
    # if the common supertype of X and Y is abstract, they cannot be compared
    if isabstracttype(promote_type(typeof(X), typeof(Y)))
        return false
    end

    for f in fieldnames(typeof(X))
        if getfield(X, f) != getfield(Y, f)
            return false
        end
    end

    return true
end
"""
    ≈(X::LazySet, Y::LazySet)

Return whether two LazySets of the same type are approximately equal.

### Input

- `X` -- any `LazySet`
- `Y` -- another `LazySet` of the same type as `X`

### Output

- `true` iff `X` is equal to `Y`.

### Notes

The check is purely syntactic and the sets need to have the same base type.
For instance, `X::VPolytope ≈ Y::HPolytope` returns `false` even if `X` and `Y`
represent the same polytope.
However `X::HPolytope{Int64} ≈ Y::HPolytope{Float64}` is a valid comparison.

### Algorithm

We recursively compare the fields of `X` and `Y` until a mismatch is found.

### Examples

```jldoctest
julia> HalfSpace([1], 1) ≈ HalfSpace([1], 1)
true

julia> HalfSpace([1], 1) ≈ HalfSpace([1.00000001], 0.99999999)
true

julia> HalfSpace([1], 1) ≈ HalfSpace([1.0], 1.0)
true

julia> Ball1([0.], 1.) ≈ Ball2([0.], 1.)
false
```
"""
function ≈(X::LazySet, Y::LazySet)
    # if the common supertype of X and Y is abstract, they cannot be compared
    if isabstracttype(promote_type(typeof(X), typeof(Y)))
        return false
    end

    for f in fieldnames(typeof(X))
        if !_isapprox(getfield(X, f), getfield(Y, f))
            return false
        end
    end

    return true
end

# hook into random API
function rand(rng::AbstractRNG, ::SamplerType{T}) where T<:LazySet
    rand(T, rng=rng)
end

"""
    copy(S::LazySet)

Return a copy of the given set by copying its values recursively.

### Input

- `S` -- any `LazySet`

### Output

A copy of `S`.

### Notes

This function performs a `copy` of each field in `S`.
See the documentation of `?copy` for further details.
"""
function copy(S::T) where {T<:LazySet}
    args = [copy(getfield(S, f)) for f in fieldnames(T)]
    BT = basetype(S)
    return BT(args...)
end

"""
    tosimplehrep(S::LazySet)

Return the simple H-representation ``Ax ≤ b`` of a set from its list of linear
constraints.

### Input

- `S` -- set

### Output

The tuple `(A, b)` where `A` is the matrix of normal directions and `b` is the
vector of offsets.

### Notes

This function only works for sets that can be represented exactly by a finite
list of linear constraints.
This fallback implementation relies on `constraints_list(S)`.
"""
tosimplehrep(S::LazySet) = tosimplehrep(constraints_list(S))

"""
    reflect(P::LazySet)

Concrete reflection of a convex set `P`, resulting in the reflected set `-P`.

### Note

This function requires that the list of constraints of the set `P` is
available, i.e. such that it can be written as
``P = \\{z ∈ ℝⁿ: ⋂ sᵢᵀz ≤ rᵢ, i = 1, ..., N\\}.``

This function can be used to implement the alternative definition of the
Minkowski Difference, which writes as
```math
A ⊖ B = \\{a − b | a ∈ A, b ∈ B\\} = A ⊕ (-B)
```
by calling `minkowski_sum(A, reflect(B))`.
"""
function reflect(P::LazySet)
    @assert applicable(constraints_list, P)  "this function " *
        "requires that the list of constraints is available, but it is not; " *
        "if the set is bounded, try overapproximating with an `HPolytope` first"

    F, g = tosimplehrep(P)
    T = isbounded(P) ? HPolytope : HPolyhedron
    return T(-F, g)
end

"""
    isuniversal(X::LazySet{N}, [witness]::Bool=false) where {N}

Check whether a given convex set is universal, and otherwise optionally compute
a witness.

### Input

- `X`       -- convex set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X`` is universal
* If `witness` option is activated:
  * `(true, [])` iff ``X`` is universal
  * `(false, v)` iff ``X`` is not universal and ``v ∉ X``

### Notes

This is a naive fallback implementation.
"""
function isuniversal(X::LazySet{N}, witness::Bool=false) where {N}
    if isbounded(X)
        result = false
    else
        error("cannot determine universality of the set")
    end

    if result
        return witness ? (true, N[]) : true
    elseif witness
        error("witness production is currently not supported")
    else
        return false
    end
end

"""
    is_interior_point(d::AbstractVector{N}, P::LazySet{N};
                      p=N(Inf), ε=_rtol(N)) where {N<:Real}

Check if the point `d` is contained in the interior of the convex set `P`.

### Input

- `d`  -- point
- `P`  -- set
- `p`  -- (optional; default: `N(Inf)`) norm of the ball used to apply the error
          tolerance
- `ε`  -- (optional; default: `_rtol(N)`) error tolerance of check

### Output

Boolean which indicates if the point `d` is contained in `P`.

### Algorithm

The implementation checks if a `Ballp` of norm `p` with center `d` and radius
`ε` is contained in the set `P`.
This is a numerical check for `d ∈ interior(P)` with error tolerance `ε`.
"""
function is_interior_point(d::AbstractVector{N}, P::LazySet{N};
                           p=N(Inf), ε=_rtol(N)) where {N<:Real}
    return Ballp(p, d, ε) ⊆ P
end

"""
    plot_recipe(X::LazySet{N}, [ε]=N(PLOT_PRECISION)) where {N}

Convert a convex set to a pair `(x, y)` of points for plotting.

### Input

- `X` -- convex set
- `ε` -- (optional, default: `PLOT_PRECISION`) approximation error bound

### Output

A pair `(x, y)` of points that can be plotted.

### Algorithm

We first assert that `X` is bounded.

One-dimensional sets are converted to an `Interval`.
We do not support three-dimensional or higher-dimensional sets at the moment.

For two-dimensional sets, we first compute a polygonal overapproximation.
The second argument, `ε`, corresponds to the error in Hausdorff distance between
the overapproximating set and `X`.
The default value `PLOT_PRECISION` is chosen such that the unit ball in the
2-norm is approximated with reasonable accuracy.
On the other hand, if you only want to produce a fast box-overapproximation of
`X`, pass `ε=Inf`.
Finally, we use the plot recipe for polygons.
"""
function plot_recipe(X::LazySet{N}, ε=N(PLOT_PRECISION)) where {N}
    @assert dim(X) <= 2 "cannot plot a $(dim(X))-dimensional $(typeof(X))"
    @assert isbounded(X) "cannot plot an unbounded $(typeof(X))"

    if dim(X) == 1
        Y = convert(Interval, X)
    else
        Y = overapproximate(X, ε)
    end
    return plot_recipe(Y, ε)
end

"""
    isoperation(X::LazySet)

Check whether the given `LazySet` is an instance of a set operation or not.

### Input

- `X` -- a `LazySet`

### Output

`true` if `X` is an instance of a set-based operation and `false` otherwise.

### Notes

The fallback implementation returns whether the set type of the input is an
operation or not using `isoperationtype`.

See also [`isoperationtype(X::Type{<:LazySet})`](@ref).

### Examples

```jldoctest
julia> B = BallInf([0.0, 0.0], 1.0);

julia> isoperation(B)
false

julia> isoperation(B ⊕ B)
true
```
"""
function isoperation(X::LazySet)
    return isoperationtype(typeof(X))
end

isoperation(::Type{<:LazySet}) = error("`isoperation` cannot be applied to a set " *
                                       "type; use `isoperationtype` instead")

"""
    isoperationtype(X::Type{<:LazySet})

Check whether the given `LazySet` type is an operation or not.

### Input

- `X` -- subtype of `LazySet`

### Output

`true` if the given set type is a set-based operation and `false` otherwise.

### Notes

The fallback for this function returns an error that `isoperationtype` is not
implemented. Subtypes of `LazySet` should dispatch on this function as required.

See also [`isoperation(X<:LazySet)`](@ref).

### Examples

```jldoctest
julia> isoperationtype(BallInf)
false

julia> isoperationtype(LinearMap)
true
```
"""
function isoperationtype(X::Type{<:LazySet})
    error("`isoperationtype` is not implemented for type $X")
end

isoperationtype(::LazySet) = error("`isoperationtype` cannot be applied to " *
                                   "a set instance; use `isoperation` instead")

"""
    isequivalent(X::LazySet, Y::LazySet)

Return whether two LazySets are equal in the mathematical sense, i.e. equivalent.

### Input

- `X` -- any `LazySet`
- `Y` -- another `LazySet`

### Output

`true` iff `X` is equivalent to `Y`.

## Algorithm

First we check `X ≈ Y`, which returns `true` if and only if `X` and `Y` have the
same type and approximately the same values (checked with `LazySets._isapprox`).
If that fails, we check the double inclusion `X ⊆ Y && Y ⊆ X`.

### Examples

```jldoctest
julia> X = BallInf([0.1, 0.2], 0.3);

julia> Y = convert(HPolytope, X);

julia> X == Y
false

julia> isequivalent(X, Y)
true
```
"""
function isequivalent(X::LazySet, Y::LazySet)
    try  # TODO temporary try-catch construct until ≈ is fixed for all set types
        if X ≈ Y
            return true
        end
    catch e
    end
    return _isequivalent_inclusion(X, Y)
end

function _isequivalent_inclusion(X::LazySet, Y::LazySet)
    return X ⊆ Y && Y ⊆ X
end

"""
    isconvextype(X::Type{<:LazySet})

Check whether the given `LazySet` type is convex.

### Input

- `X` -- subtype of `LazySet`

### Output

`true` if the given set type is guaranteed to be convex by using only type
information, and `false` otherwise.

### Notes

Since this operation only acts on types (not on values), it can return false
negatives, i.e. there may be instances where the set is convex, even though the
answer of this function is `false`. The examples below illustrate this point.

### Examples

A ball in the infinity norm is always convex, hence we get:

```jldoctest convex_types
julia> isconvextype(BallInf)
true
```

For instance, the union (`UnionSet`) of two sets may in general be either convex
or not, since convexity cannot be decided by just using type information.
Hence, `isconvextype` returns `false` if `X` is `Type{<:UnionSet}`.

```jldoctest convex_types
julia> isconvextype(UnionSet)
false
```

However, the type parameters from the set operations allow to decide convexity
in some cases, by falling back to the convexity of the type of its arguments.
Consider for instance the lazy intersection. The intersection of two convex sets
is always convex, hence we can get:

```jldoctest convex_types
julia> isconvextype(Intersection{Float64, BallInf{Float64}, BallInf{Float64}})
true
```
"""
isconvextype(X::Type{<:LazySet}) = false

"""
    surface(X::LazySet{N}) where {N}

Compute the surface area of a set.

### Input

- `X` -- set

### Output

A number representing the surface area of `X`.
"""
function surface(X::LazySet{N}) where {N}
    if dim(X) == 2
        return area(X)
    else
        throw(ArgumentError("the surface function is only implemented for " *
                    "two-dimensional sets, but the given set is $(dim(X))-dimensional"))
    end
end

"""
    area(X::LazySet{N}) where {N}

Compute the area of a two-dimensional polytopic set using the Shoelace formula.

### Input

- `X` -- two-dimensional set

### Output

A number representing the area of `X`.

### Notes

This algorithm is applicable to any lazy set `X` such that its list of vertices,
`vertices_list`, can be computed.

### Algorithm

Let `m` be the number of vertices of `X`. The following instances are considered:

- `m = 0, 1, 2`: the output is zero.
- `m = 3`: the triangle case is computed using the Shoelace formula with 3 points.
- `m = 4`: the quadrilateral case is obtained by the factored version of the Shoelace
           formula with 4 points.

Otherwise, the general Shoelace formula is used; for detals see the wikipedia
article [Shoelace formula](https://en.wikipedia.org/wiki/Shoelace_formula).
"""
function area(X::LazySet{N}) where {N}
    @assert dim(X) == 2 "this function only applies to two-dimensional sets, " *
    "but the given set is $(dim(X))-dimensional"

    Xpoly = convert(VPolygon, X)
    vlist = vertices_list(Xpoly)
    m = length(vlist)

    if m <= 2
        return zero(N)
    end

    if m == 3 # triangle
        res = _area_triangle(vlist)

    elseif m == 4 # quadrilateral
        res = _area_quadrilateral(vlist)

    else # general case
        res = _area_polygon(vlist)
    end

    return res
end

function _area_triangle(v::Vector{VN}) where {N, VN<:AbstractVector{N}}
    A = v[1]
    B = v[2]
    C = v[3]
    res = A[1] * (B[2] - C[2]) + B[1] * (C[2] - A[2]) + C[1] * (A[2] - B[2])
    return abs(res/2)
end

function _area_quadrilateral(v::Vector{VN}) where {N, VN<:AbstractVector{N}}
    A = v[1]
    B = v[2]
    C = v[3]
    D = v[4]
    res = A[1] * (B[2] - D[2]) + B[1] * (C[2] - A[2]) + C[1] * (D[2] - B[2]) + D[1] * (A[2] - C[2])
    return abs(res/2)
end

function _area_polygon(v::Vector{VN}) where {N, VN<:AbstractVector{N}}
    m = length(v)
    @inbounds res = v[m][1] * v[1][2] - v[1][1] * v[m][2]
    for i in 1:m-1
        @inbounds res += v[i][1] * v[i+1][2] - v[i+1][1] * v[i][2]
    end
    return abs(res/2)
end

"""
    singleton_list(P::LazySet)

Return the vertices of a polytopic set as a list of singletons.

### Input

- `P`                 -- polytopic set

### Output

The list of vertices of `P`, as `Singleton`.

### Notes

This function relies on `vertices_list`, which raises an error if the set is
not polytopic (e.g., unbounded).
"""
function singleton_list(P::LazySet)
    return [Singleton(x) for x in vertices_list(P)]
end

"""
    concretize(X::LazySet)

Construct a concrete representation of a (possibly lazy) set.

### Input

- `X` -- set

### Output

A concrete representation of `X` (as far as possible).

### Notes

Since not every lazy set has a concrete set representation in this library, the
result may be partially lazy.
"""
function concretize(X::LazySet)
    return X
end

"""
    constraints(X::LazySet)

Construct an iterator over the constraints of a polyhedral set.

### Input

- `X` -- polyhedral set

### Output

An iterator over the constraints of `X`.
"""
function constraints(X::LazySet)
    return _constraints_fallback(X)
end

"""
    vertices(X::LazySet)

Construct an iterator over the vertices of a polyhedral set.

### Input

- `X` -- polyhedral set

### Output

An iterator over the vertices of `X`.
"""
function vertices(X::LazySet)
    return _vertices_fallback(X)
end

function _constraints_fallback(X::LazySet)
    return VectorIterator(constraints_list(X))
end

function _vertices_fallback(X::LazySet)
    return VectorIterator(vertices_list(X))
end

# =========================
# Code requiring MiniQhull
# =========================

function load_delaunay_MiniQhull()
return quote

import .MiniQhull: delaunay
export delaunay

"""
    delaunay(X::LazySet)

Compute the Delaunay triangulation of the given convex set.

### Input

- `X` -- set

### Output

A union of polytopes in vertex representation.

### Notes

This function requires that you have properly installed the package
[MiniQhull.jl](https://github.com/gridap/MiniQhull.jl), including the library
[Qhull](http://www.qhull.org/).

The method works in arbitrary dimension and the requirement is that the list of
vertices of `X` can be obtained.
"""
function delaunay(X::LazySet)
    n = dim(X)
    v = vertices_list(X)
    m = length(v)
    coordinates = vcat(v...)
    connectivity_matrix = delaunay(n, m, coordinates)
    nelements = size(connectivity_matrix, 2)
    elements = [VPolytope(v[connectivity_matrix[:, i]]) for i in 1:nelements]
    return UnionSetArray(elements)
end

end end  # load_delaunay_MiniQhull

"""
    complement(X::LazySet)

Return the complement of a set.

### Input

- `X` -- set

### Output

A `UnionSetArray` of half-spaces, i.e. the output is the union of the linear
constraints which are obtained by complementing each constraint of `X`.

### Algorithm

The principle used in this function is that if ``X`` and ``Y`` are any pair of sets,
then ``(X ∩ Y)^C = X^C ∪ Y^C``. In particular, we can apply this rule for each constraint
that defines a polyhedral set, hence the concrete complement can be represented as the set
union of the complement of each constraint.
"""
function complement(X::LazySet)
    return UnionSetArray(constraints_list(Complement(X)))
end

# -- concrete projection --

"""
    project(S::LazySet{N},
            block::AbstractVector{Int},
            [::Nothing=nothing],
            [n]::Int=dim(S);
            [kwargs...]
           ) where {N}

Project a high-dimensional set to a given block by using a concrete linear map.

### Input

- `S`       -- set
- `block`   -- block structure - a vector with the dimensions of interest
- `nothing` -- (default: `nothing`) used for dispatch
- `n`       -- (optional, default: `dim(S)`) ambient dimension of the set `S`

### Output

A set representing the projection of the set `S` to block `block`.

### Algorithm

We apply the function `linear_map`.
"""
@inline function project(S::LazySet{N},
                         block::AbstractVector{Int},
                         ::Nothing=nothing,
                         n::Int=dim(S);
                         kwargs...
                        ) where {N}
    return _project_linear_map(S, block, n; kwargs...)
end

@inline function _project_linear_map(S::LazySet{N},
                                     block::AbstractVector{Int},
                                     n::Int=dim(S);
                                     kwargs...
                                    ) where {N}
    M = projection_matrix(block, n, N)
    return linear_map(M, S)
end

"""
    project(S::LazySet,
            block::AbstractVector{Int},
            set_type::Type{TS},
            [n]::Int=dim(S);
            [kwargs...]
           ) where {TS<:LazySet}

Project a high-dimensional set to a given block and set type, possibly involving
an overapproximation.

### Input

- `S`        -- set
- `block`    -- block structure - a vector with the dimensions of interest
- `set_type` -- target set type
- `n`        -- (optional, default: `dim(S)`) ambient dimension of the set `S`

### Output

A set of type `set_type` representing an overapproximation of the projection of
`S`.

### Algorithm

1. Project the set `S` with `M⋅S`, where `M` is the identity matrix in the block
coordinates and zero otherwise.
2. Overapproximate the projected lazy set using `overapproximate` and
`set_type`.
"""
@inline function project(S::LazySet,
                         block::AbstractVector{Int},
                         set_type::Type{TS},
                         n::Int=dim(S);
                         kwargs...
                        ) where {TS<:LazySet}
    lm = project(S, block, LinearMap, n)
    return overapproximate(lm, set_type)
end

"""
    project(S::LazySet,
            block::AbstractVector{Int},
            set_type_and_precision::Pair{T, N},
            [n]::Int=dim(S);
            [kwargs...]
           ) where {T<:UnionAll, N<:Real}

Project a high-dimensional set to a given block and set type with a certified
error bound.

### Input

- `S`     -- set
- `block` -- block structure - a vector with the dimensions of interest
- `set_type_and_precision` -- pair `(T, ε)` of a target set type `T` and an
                              error bound `ε` for approximation
- `n`     -- (optional, default: `dim(S)`) ambient dimension of the set `S`

### Output

A set representing the epsilon-close approximation of the projection of `S`.

### Notes

Currently we only support `HPolygon` as set type, which implies that the set
must be two-dimensional.

### Algorithm

1. Project the set `S` with `M⋅S`, where `M` is the identity matrix in the block
coordinates and zero otherwise.
2. Overapproximate the projected lazy set with the given error bound `ε`.
"""
@inline function project(S::LazySet,
                         block::AbstractVector{Int},
                         set_type_and_precision::Pair{T, N},
                         n::Int=dim(S);
                         kwargs...
                        ) where {T<:UnionAll, N<:Real}
    set_type = set_type_and_precision[1]
    ε = set_type_and_precision[2]
    @assert length(block) == 2 && set_type == HPolygon "currently only 2D " *
        "HPolygon decomposition is supported"

    lm = project(S, block, LinearMap, n)
    return overapproximate(lm, set_type, ε)
end

"""
    project(S::LazySet,
            block::AbstractVector{Int},
            ε::Real,
            [n]::Int=dim(S);
            [kwargs...]
           )

Project a high-dimensional set to a given block and set type with a certified
error bound.

### Input

- `S`     -- set
- `block` -- block structure - a vector with the dimensions of interest
- `ε`     -- error bound for approximation
- `n`     -- (optional, default: `dim(S)`) ambient dimension of the set `S`

### Output

A set representing the epsilon-close approximation of the projection of `S`.

### Algorithm

1. Project the set `S` with `M⋅S`, where `M` is the identity matrix in the block
coordinates and zero otherwise.
2. Overapproximate the projected lazy set with the given error bound `ε`.
The target set type is chosen automatically.
"""
@inline function project(S::LazySet,
                         block::AbstractVector{Int},
                         ε::Real,
                         n::Int=dim(S);
                         kwargs...
                        )
    # currently we only support HPolygon
    if length(block) == 2
        set_type = HPolygon
    else
        throw(ArgumentError("ε-close approximation is only supported for 2D " *
                            "blocks"))
    end
    return project(S, block, set_type => ε, n)
end

"""
    rectify(X::LazySet, [concrete_intersection]::Bool=false)

Concrete rectification of a set.

### Input

- `X`                     -- set
- `concrete_intersection` -- (optional, default: `false`) flag to compute
                             concrete intersections for intermediate results

### Output

A set corresponding to the rectification of `X`, which is in general a union of
linear maps of intersections.

### Algorithm

For each dimension in which `X` is both positive and negative we split `X` into
these two parts. Additionally we project the negative part to zero.
"""
function rectify(X::LazySet, concrete_intersection::Bool=false)
    return to_union_of_projections(Rectification(X), concrete_intersection)
end

"""
    low(X::LazySet, i::Int)

Return the lower coordinate of a set in a given dimension.

### Input

- `X` -- set
- `i` -- dimension of interest

### Output

The lower coordinate of the set in the given dimension.
"""
function low(X::LazySet{N}, i::Int) where {N}
    n = dim(X)
    d = SingleEntryVector(i, n, -one(N))
    return -ρ(d, X)
end

"""
    low(X::LazySet)

Return a vector with the lowest coordinates of the set for each canonical direction.

### Input

- `X` -- set

### Output

A vector with the lower coordinate of the set for each dimension.

### Notes

See also [`low(X::LazySet, i::Int)`](@ref).
"""
function low(X::LazySet)
    n = dim(X)
    return [low(X, i) for i in 1:n]
end

"""
    high(X::LazySet, i::Int)

Return the higher coordinate of a set in a given dimension.

### Input

- `X` -- set
- `i` -- dimension of interest

### Output

The higher coordinate of the set in the given dimension.
"""
function high(X::LazySet{N}, i::Int) where {N}
    n = dim(X)
    d = SingleEntryVector(i, n, one(N))
    return ρ(d, X)
end

"""
    high(X::LazySet)

Return a vector with the highest coordinate of the set for each canonical direction.

### Input

- `X` -- set

### Output

A vector with the highest coordinate of the set for each dimension.

### Notes

See also [`high(X::LazySet, i::Int)`](@ref).
"""
function high(X::LazySet)
    n = dim(X)
    return [high(X, i) for i in 1:n]
end

"""
    extrema(X::LazySet)

Return two vectors with the lowest and highest coordinate of `X` for each
dimension.

### Input

- `X` -- set

### Output

Two vectors with the lowest and highest coordinates of `X` for each dimension.

### Notes

The result is equivalent to `(low(X), high(X))`, but sometimes it can be
computed more efficiently.
"""
function extrema(X::LazySet)
    l = low(X)
    h = high(X)
    return (l, h)
end

"""
    extrema(X::LazySet, i::Int)

Return the lower and higher coordinate of a set in a given dimension.

### Input

- `X` -- set
- `i` -- dimension of interest

### Output

The lower and higher coordinate of the set in the given dimension.

### Notes

The result is equivalent to `(low(X, i), high(X, i))`, but sometimes it can be
computed more efficiently.
"""
function extrema(X::LazySet, i::Int)
    l = low(X, i)
    h = high(X, i)
    return (l, h)
end

"""
    rationalize(::Type{T}, X::LazySet{N}, tol::Real) where {T<:Integer, N<:AbstractFloat}

Approximate a LazySet of floating point numbers as a set whose entries are
rationals of the given integer type.

### Input

- `T`   -- (optional, default: `Int`) integer type to represent the rationals
- `X`   -- set which has floating-point components
- `tol` -- (optional, default: `eps(N)`) tolerance of the result; each rationalized
           component will differ by no more than `tol` with respect to the floating-point value

### Output

A LazySet of the same base type of `X` where each numerical component is of
type `Rational{T}`.
"""
function rationalize(::Type{T}, X::LazySet{N}, tol::Real) where {T<:Integer, N<:AbstractFloat}
    m = length(fieldnames(typeof(X)))
    frat = ntuple(fi -> _rationalize(T, getfield(X, fi), tol), m)
    ST = basetype(X)
    return ST(frat...)
end

rationalize(X::LazySet{N}; kwargs...) where {N<:AbstractFloat} = rationalize(Int, X; kwargs...)
rationalize(::Type{T}, X::LazySet{N}; tol::Real=eps(N)) where {T<:Integer, N<:AbstractFloat} = rationalize(T, X, tol)

# method extension for lazy sets
_rationalize(::Type{T}, X::AbstractVector{<:LazySet{N}}, tol::Real) where {T<:Integer, N<:AbstractFloat} = rationalize.(Ref(T), X, Ref(tol))
_rationalize(::Type{T}, X::LazySet{N}, tol::Real) where {T<:Integer, N<:AbstractFloat} = rationalize(T, X, tol)

"""
    permute(X::LazySet, p::AbstractVector{Int})

Permute the dimensions of a set according to a given permutation vector.

### Input

- `X` -- set
- `p` -- permutation vector

### Output

A new set corresponding to `X` where the dimensions have been permuted according
to `p`.
"""
function permute end


text:
#   This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2016 Nathan Smith

const EOT = Char(4)					# end of transmission character
const synerr = "incorrect syntax"
const maxerr = "-- an error"

struct MaximaSession <: Base.AbstractPipe

    input::Pipe
    output::Pipe
    process::Base.Process

    function MaximaSession(;args::Cmd=``)
        # If windows, executable is .bat
        cmd = Sys.isunix() ? `maxima --very-quiet $args` :
            `maxima.bat --very-quiet $args`

        # Setup pipes and maxima process
        input = Pipe()
        output = Pipe()
        pipe = pipeline(cmd, stdin=input, stdout=output)
        process = run(pipe, wait=false)

        # Close the unneeded ends of Pipes
        close(input.out)
        close(output.in)

        # Set display to 1-dimensionsal
        write(input, "display2d: false\$")

        return new(input, output, process)
    end
end

Base.kill(ms::MaximaSession) = kill(ms.process)
Base.process_exited(ms::MaximaSession) = process_exited(ms.process)

function Base.write(ms::MaximaSession, input::String)
    # The line break right ..v.. there is apparently very important...
    write(ms.input, "$input;\n")
    write(ms.input, "print(ascii(4))\$")
end

function Base.read(ms::MaximaSession)
    readuntil(ms.output, EOT) |> String |> str -> rstrip(str, EOT)
end


text:
struct Interval{T<:Real,A<:Bound{T},B<:Bound{T}}
    a::A
    b::B
    function Interval{T}(a::A, b::B) where {T<:Real,A<:Bound{T},B<:Bound{T}}
        if !(A <: NullBound || B <: NullBound)
            va = a.value
            vb = b.value
            if A <: ClosedBound && B <: ClosedBound
                va <= vb || error("Invalid bounds: a=$va must be less than or equal to b=$vb")
            else
                va < vb || error("Invalid bounds: a=$va must be less than b=$vb")
            end
        end
        new{T,A,B}(a,b)
    end
end
Interval(a::Bound{T}, b::Bound{T}) where {T<:Real} = Interval{T}(a,b)

eltype(::Interval{T}) where {T} = T

interval(a::Nothing, b::Nothing) = Interval(NullBound{Float64}(), NullBound{Float64}())
interval(a::Bound{T}, b::Nothing) where {T<:Real} = Interval(a, NullBound{T}())
interval(a::Nothing, b::Bound{T}) where {T<:Real} = Interval(NullBound{T}(), b)
interval(a::Bound{T}, b::Bound{T}) where {T<:Real} = Interval(a,b)
interval(::Type{T}) where {T<:Real} = Interval(NullBound{T}(), NullBound{T}())

checkvalue(I::Interval, x::Real) = checkvalue(I.a, x) && checkvalue(x, I.b)

function theta(I::Interval{T,A,B}, x::T) where {T<:AbstractFloat,A,B}
    checkvalue(I,x) || throw(DomainError(x))
    if A <: OpenBound
        return B <: OpenBound ? log(x-I.a.value) - log(I.b.value-x) : log(x-I.a.value)
    else
        return B <: OpenBound ? log(I.b.value-x) : x
    end
end

function deriv_theta(I::Interval{T,A,B}, x::T) where {T<:AbstractFloat,A,B}
    checkvalue(I,x) || throw(DomainError(x))
    if A <: OpenBound
        return B <: OpenBound ? one(T)/(x-I.a.value) + one(T)/(I.b.value-x) : one(T)/(x-I.a.value)
    else
        return B <: OpenBound ? one(T)/(I.b.value-x) : one(T)
    end
end

function upperboundtheta(I::Interval{T,A,B}) where {T<:AbstractFloat,A,B}
    if B <: ClosedBound
        return A <: OpenBound ? log(I.b.value - I.a.value) : I.b.value
    elseif B <: OpenBound
        return A <: ClosedBound ? log(I.b.value - I.a.value) : convert(T,Inf)
    else
        return convert(T,Inf)
    end
end

function lowerboundtheta(I::Interval{T,A,B}) where {T<:AbstractFloat,A,B}
    A <: ClosedBound && !(B <: OpenBound) ? I.a.value : convert(T,-Inf)
end

function checktheta(I::Interval{T}, x::T) where {T<:AbstractFloat}
    lowerboundtheta(I) <= x <= upperboundtheta(I)
end

function eta(I::Interval{T,A,B}, x::T) where {T<:AbstractFloat,A,B}
    checktheta(I,x) || throw(DomainError(x))
    if A <: OpenBound
        if B <: OpenBound
            return (I.b.value*exp(x) + I.a.value)/(one(T) + exp(x))
        else
            return exp(x) + I.a.value
        end
    else
        return B <: OpenBound ? I.b.value - exp(x) : x
    end
end

function deriv_eta(I::Interval{T,A,B}, x::T) where {T<:AbstractFloat,A,B}
    checktheta(I,x) || throw(DomainError(x))
    if A <: OpenBound
        if B <: OpenBound
            return (I.a.value - I.b.value)*exp(x)/(one(T) + exp(x))^2
        else
            return exp(x)
        end
    else
        return B <: OpenBound ? -exp(x) : x
    end
end

function string(I::Interval{T1,T2,T3}) where {T1,T2,T3}
    if T2 <: NullBound
        if T3 <: NullBound
            string("interval(", T1, ")")
        else
            string("interval(nothing,", string(I.b), ")")
        end
    else
        string("interval(", string(I.a), ",", T3 <: NullBound ? "nothing" : string(I.b), ")")
    end
end

function show(io::IO, I::Interval)
    print(io, string(I))
end


text:
function exposurize_model(model; exposure = 1)
    result = Array{Union{Nothing, LRMoE.AnyExpert}}(nothing, size(model)[1], size(model)[2], length(exposure))
    for i in 1:length(exposure)
        result[:,:,i] = LRMoE.exposurize_expert.(model, exposure = exposure[i])
    end
    return result
end

function loglik_aggre_dim(ls)
    return Matrix(sum(ls, dims = 1)[1,:,:]')
end

function loglik_aggre_gate_dim(gate, ll_mat)
    return gate + ll_mat
end

function loglik_aggre_gate_dim_comp(ll_mat)
    return rowlogsumexp(ll_mat)
end

# Exact case
function expert_ll_ind_mat_exact(Y, model)
    result = fill(NaN, size(model)[1], size(model)[2])
    for j in 1:size(model)[1]
        for k in 1:size(model)[2]
            result[j, k] = LRMoE.expert_ll_exact.(model[j, k], Y[j])
        end
    end
    return result # cat([[LRMoE.expert_ll_exact.(model[j, k], Y[j]) for j in 1:size(model)[1]] for k in 1:size(model)[2]]..., dims = 2)
end

function expert_ll_list_exact(Y, model)
    result = fill(NaN, size(model)[1], size(model)[2], size(Y)[1])
    for i in 1:size(Y)[1]
        result[:,:,i] = expert_ll_ind_mat_exact(Y[i,:], model[:,:,i])
    end
    return result # cat([expert_ll_ind_mat_exact(Y[i,:], model[:,:,i]) for i in 1:size(Y)[1]]..., dims = 3)
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
# power analysis functions for difference in two proportions

"""
    powerTwoPTest(h::Real = 0, n::Real = 0, alpha::Float64 = 0.05, alternative = "two.sided")

Compute power of two equal size samples of `n` observations to test the effect size `h`
for a difference between two independent proportions (P) at type I error = `alpha` (default: 0.05).
The effect size `h` = ϕ₁ - ϕ₂ where ϕ = 2*asin(sqrt(P)). The option
`alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function powerTwoPTest(;
    h::Real = 0,
    n::Real = 0,
    alpha::Float64 = 0.05,
    alternative::String = "two"
    )

    check_args(h=h,n=n,alpha=alpha)

    if alternative == "less"
        tside = 1
    elseif alternative in ("two","two.sided","two-sided")
        tside = 2
        h = abs(h)
    elseif alternative == "greater"
        tside = 3
    else
        error("`alternative` must be `less`, `two`, or `greater`.")
    end

    if tside == 1
        return cdf(Normal(),quantile(Normal(),alpha) - h*sqrt(n/2))
    elseif tside == 2
        return ccdf(Normal(),cquantile(Normal(),alpha/2) - h*sqrt(n/2)) + cdf(Normal(),quantile(Normal(),alpha/2) - h*sqrt(n/2))
    elseif tside == 3
        return ccdf(Normal(),cquantile(Normal(),alpha) - h*sqrt(n/2))
    end
    error("internal error")
end

"""
    samplesizeTwoPTest(h::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the sample size `n` of two equal size samples to test the effect size `h`
for a difference between two independent proportions (P) with power > `power` (default: 0.8)
at type I error = `alpha` (default: 0.05).
The effect size `h` = ϕ₁ - ϕ₂ where ϕ = 2*asin(sqrt(P)). The option
`alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function samplesizeTwoPTest(;
    h::Real = 0,
    alpha::Float64 = 0.05,
    power::Float64 = 0.8,
    alternative::String = "two"
    )

    check_args(h=h,alpha=alpha,power=power)

    return ceil(Int64,fzero(x->powerTwoPTest(h = h, n = x, alpha = alpha, alternative = alternative) - power, 2 + 1e-10, 1e+09))
end

"""
    effectsizeTwoPTest(n::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the effect size that two equal size samples of `n` observations can detect for
a difference between two independent proportions (P) with power > `power` (default: 0.8)
at type I error = `alpha` (default: 0.05).
The effect size `h` = ϕ₁ - ϕ₂ where ϕ = 2*asin(sqrt(P)). The option
`alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function effectsizeTwoPTest(;
    n::Real = 0,
    alpha::Float64 = 0.05,
    power::Float64 = 0.8,
    alternative::String = "two"
    )

    check_args(n=n,alpha=alpha,power=power)

    if alternative == "less"
        return fzero(x->powerTwoPTest(h = x, n = n, alpha = alpha, alternative = alternative) - power, -10.0, 5.0)
    elseif alternative == "two"
        return fzero(x->powerTwoPTest(h = x, n = n, alpha = alpha, alternative = alternative) - power, 1e-10, 10.0)
    elseif side == "greater"
        return fzero(x->powerTwoPTest(h = x, n = n, alpha = alpha, alternative = alternative) - power, -5.0, 10.0)
    end
end

"""
    alphaTwoPTest(h::Real = 0, n::Real = 0, power::Float64 = 0.8, alternative = "two.sided")

Compute the probability of type I error that two equal size samples of `n` observations have in detecting
the effect size `h` for a difference between two independent proportions (P) with power > `power` (default: 0.8)
at type I error = `alpha` (default: 0.05).
The effect size `h` = ϕ₁ - ϕ₂ where ϕ = 2*asin(sqrt(P)). The option
`alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function alphaTwoPTest(;
    h::Real = 0,
    n::Real = 0,
    power::Float64 = 0.8,
    alternative::String = "two"
    )

    check_args(h=h,n=n,power=power)

    return fzero(x->powerTwoPTest(h = h, n = n, alpha = x, alternative = alternative) - power, 1e-10, 1 - 1e-10)
end

"""
    TwoPTest(h::Real = 0, n::Real = 0, alpha::Float64 = 0.0, power::Float64 = 0.0, alternative = "two.sided")

Compute one of the test parameters such as sample size (`n`),
effect size (`h`), type I error (`alpha`), and power (`power`).
The parameter to be estimated must be set to zero (default).
The effect size `h` = ϕ₁ - ϕ₂ where ϕ = 2*asin(sqrt(P)).
The option `alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function TwoPTest(;
    h::Real = 0,
    n::Real = 0,
    alpha::Float64 = 0.0,
    power::Float64 = 0.0,
    alternative::String = "two"
    )
    if sum([x == 0 for x in (h,n,alpha,power)]) != 1
        error("exactly one of `h`, `n`, `power`, and `alpha` must be zero")
    end

    if power == 0.0
        power = powerTwoPTest(h = h, n = n, alpha = alpha, alternative = alternative)
    elseif alpha == 0.0
        alpha = alphaTwoPTest(h = h, n = n, power = power, alternative = alternative)
    elseif h == 0
        h = effectsizeTwoPTest(n = n, alpha = alpha, power = power, alternative = alternative)
    elseif n == 0
        n = samplesizeTwoPTest(h = h, alpha = alpha, power = power, alternative = alternative)
    end

    alt = Dict("two" => "two-sided", "two.sided" => "two-sided","less" => "less", "greater" => "greater")

    return htest(
        "Difference of proportion power calculation for binomial distribution (same sample)",
        OrderedDict(
            "h" => h,
            "n" => n,
            "alpha" => alpha,
            "power" => power,
            "alternative" => alt[alternative])
        )
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
# ---
# title: 1053. Previous Permutation With One Swap
# id: problem1053
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Greedy
# link: <https://leetcode.com/problems/previous-permutation-with-one-swap/description/>
# hidden: true
# ---
#
# Given an array `A` of positive integers (not necessarily distinct), return the
# lexicographically largest permutation that is smaller than `A`, that can be
# **made with one swap** (A _swap_ exchanges the positions of two numbers `A[i]`
# and `A[j]`).  If it cannot be done, then return the same array.
#
#
#
# **Example 1:**
#
#
#
#     Input: [3,2,1]
#     Output: [3,1,2]
#     Explanation: Swapping 2 and 1.
#
#
# **Example 2:**
#
#
#
#     Input: [1,1,5]
#     Output: [1,1,5]
#     Explanation: This is already the smallest permutation.
#
#
# **Example 3:**
#
#
#
#     Input: [1,9,4,6,7]
#     Output: [1,7,4,6,9]
#     Explanation: Swapping 9 and 7.
#
#
# **Example 4:**
#
#
#
#     Input: [3,1,1,3]
#     Output: [1,3,1,3]
#     Explanation: Swapping 1 and 3.
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 10000`
#   2. `1 <= A[i] <= 10000`
#
#
## @lc code=start
using LeetCode

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


text:
module MLE

using ProgressMeter
using LinearAlgebra, Statistics, StatsBase
using GSL, Optim, NLSolversBase
using SpecialFunctions: erfinv, loggamma

incbeta(a,b,x) = GSL.sf_beta_inc(a,b,x)
incgamma(a,x)  = GSL.sf_gamma_inc_P(a,x)

const ∞ = Inf
function clamp(x,lo,hi)
    x < lo && return lo
    x > hi && return hi
    return x
end

# ---------------------------------------------------------------------
# generalized linear models
# Γ denotes parameters for prior distributions

# the common formulation i.e. NB2
function negative_binomial(x⃗, z⃗, Γ)
    β̄,δβ¯²,Γᵧ = Γ

    f = if isnothing(Γᵧ)
        (Θ) -> let
            α,β,γ = Θ
            L  = (loggamma(x+γ)
                - loggamma(x+1)
                - loggamma(γ)
                + x*(α+β*z)
                + γ*log(γ)
                - (x+γ)*log(exp(α+β*z)+γ) for (x,z) ∈ zip(x⃗,z⃗))

            return -sum(L) + 0.5*δβ¯²*(β-β̄)^2
        end
    else
        μ,σ,ν = Γᵧ
        (Θ) -> let
            α,β,γ = Θ
            L  = (loggamma(x+γ)
                - loggamma(x+1)
                - loggamma(γ)
                + x*(α+β*z)
                + γ*log(γ)
                - (x+γ)*log(exp(α+β*z)+γ) for (x,z) ∈ zip(x⃗,z⃗))

            return -sum(L) + 0.5*δβ¯²*(β-β̄)^2 + (abs(log(γ)-μ)/σ)^ν
        end
    end

    # TODO: 1st and 2nd derivatives

    μ  = mean(x⃗)
    Θ₀ = [
        log(μ),
        β̄,
        var(x⃗)/μ - 1,
    ]

    if Θ₀[end] < 0 || isinf(Θ₀[end])
        Θ₀[end] = 1
    end

    return (
        Θ₀         = Θ₀,
        loss       = TwiceDifferentiable(f, Θ₀; autodiff=:forward),
        constraint = TwiceDifferentiableConstraints([-∞,-∞,0],[+∞,+∞,+∞]),
        residual   = (Θ) -> let
            α, β, γ = Θ
            μ = @. exp(α + β*z⃗)
            σ = @. √(μ + μ^2/γ)
            return @. (x⃗ - μ) / σ
        end,
        cumulative = (Θ) -> let
            α, β, γ = Θ
            μ = @. exp(α + β*z⃗)
            p = @. μ / (μ + γ)
            return @. 1 - incbeta(x⃗+1, γ, p)
        end
    )
end

function gamma(x⃗, z⃗, Γ)
    β̄, δβ¯² = Γ

    function f(Θ)
        α, β, γ = Θ

        k⃗ = (exp(α+β*z)/γ for z ∈ z⃗)
        L = (-loggamma(k)
             -k*log(γ)
             +(k-1)*log(x)
             -x/γ for (k,x) ∈ zip(k⃗,x⃗))

        return -sum(L) + 0.5*δβ¯²*(β-β̄)^2
    end

    μ  = mean(x⃗)
    Θ₀ = [
        log(μ),
        β̄,
        μ^2 / (var(x⃗)-μ),
    ]

    if Θ₀[end] < 0 || isinf(Θ₀[end])
        Θ₀[end] = 1
    end

    return (
        Θ₀         = Θ₀,
        loss       = TwiceDifferentiable(f, Θ₀; autodiff=:forward),
        constraint = TwiceDifferentiableConstraints([-∞,-∞,0],[+∞,+∞,+∞]),
        residual   = (Θ) -> let
            α, β, γ = Θ
            μ = @. exp(α + β*z⃗)
            k = @. μ / γ
            Φ = @. incgamma(k, x⃗/γ)
            ρ = @. erfinv(clamp(2*Φ-1,-1,1))
            ρ[isinf.(ρ)] = 10*sign.(ρ[isinf.(ρ)])
            return ρ
        end,
        cumulative = (Θ) -> let
            α, β, γ = Θ
            μ = @. exp(α + β*z⃗)
            k = @. μ / γ
            return @. incgamma(k, x⃗/γ)
        end
    )
end

# ---------------------------------------------------------------------
# simple distributions (no latent variables)

# univariate gamma
function gamma(x)
    f = (Θ) -> let
        k, θ = Θ
        return -sum((k-1)*log.(x) .- (x./θ) .- loggamma(k) .- k*log(θ))
    end

    μ  = mean(x)
    σ² = var(x)
    Θ₀ = [
        μ^2/σ²,
        σ²/μ,
    ]

    return (
        Θ₀         = Θ₀,
        loss       = TwiceDifferentiable(f, Θ₀; autodiff=:forward),
        constraint = TwiceDifferentiableConstraints([0,0],[+∞,+∞]),
        cumulative = (Θ) -> let
            k, θ = Θ
            return @. incgamma(k, x/θ)
        end
    )
end

# univariate lognormal
function log_normal(x)
    f = (Θ) -> let
        μ, σ = Θ
        return sum(log(σ) .+ (log.(x) .- μ).^2 ./ (2*σ^2))
    end

    μ  = mean(log.(x))
    σ  = std(log.(x))
    Θ₀ = [
        μ,
        σ
    ]

    return (
        Θ₀         = Θ₀,
        loss       = TwiceDifferentiable(f, Θ₀; autodiff=:forward),
        constraint = TwiceDifferentiableConstraints([0,0],[+∞,+∞]),
        cumulative = (Θ) -> let
            μ, σ = Θ
            return @. .5*(1+erf((log(x) - μ)/(√2*σ)))
        end
    )
end

# univariate generalized normal (variant 1 wikipedia)
function generalized_normal(x)
    f = (Θ) -> let
        μ, σ, β = Θ
        return sum(loggamma(1/β) + log(σ) .+ (abs.(x .- μ)./ σ).^β .- log(β))
    end

    μ  = mean(x)
    σ  = std(x)
    Θ₀ = [
        μ,
        σ,
        2
    ]

    return (
        Θ₀         = Θ₀,
        loss       = TwiceDifferentiable(f, Θ₀; autodiff=:forward),
        constraint = TwiceDifferentiableConstraints([0,0,0],[+∞,+∞,+∞]),
        cumulative = (Θ) -> let
            μ, σ, β = Θ
            return @. .5*(1+sign(x-μ)*incgamma((1/β), abs(x-μ)/(σ^β)))
        end
   )
end

# ---------------------------------------------------------------------
# drivers

const FitType = NamedTuple{
    (
        :likelihood,
        :parameters,
        :uncertainty,
        :cumulative,
        :residuals,
    ),
    Tuple{
        Float64,
        Array{Float64,1},
        Array{Float64,1},
        Array{Float64,1},
        Array{Float64,1},
    }
}

# log link function is assumed
function fit_glm(model::Symbol, data; Γ=(β̄=1,δβ¯²=10,Γᵧ=nothing), run=(x)->true)
    Σ = vec(mean(data, dims=1))

    foundmodel = try
        getfield(MLE, model)
    catch
        error("model '$model' not implemented")
    end
    modelfor(x) = foundmodel(x, log.(Σ), Γ)

    function fit(row, i)::FitType
        model    = modelfor(vec(row))
        estimate = optimize(model.loss, model.constraint, model.Θ₀, IPNewton())

        Θ  = Optim.minimizer(estimate)
        E  = Optim.minimum(estimate)
        δΘ = diag(inv(hessian!(model.loss, Θ)))

        return (
            likelihood  = E,
            parameters  = Θ,
            uncertainty = δΘ,
            cumulative  = model.cumulative(Θ),
            residuals   = model.residual(Θ),
        )
    end

    ι = zeros(Int, size(data,1))
    j = 1
    for (i, row) ∈ enumerate(eachrow(data))
        if run(row)
            ι[i] = j
            j   += 1
        end
    end

    fits     = Array{FitType,1}(undef,j-1)
    progress = Progress(sum(ι .!= 0); desc="--> fitting:", output=stderr)
    for (i, row) ∈ enumerate(eachrow(data))
        ι[i] == 0 && continue

        fits[ι[i]] = fit(row,i)
        next!(progress)
    end

    return (
        likelihood  = map((f)->f.likelihood,  fits),

        # XXX: assumes particular form for parameters
        α  = map((f)->f.parameters[1],  fits),
        β  = map((f)->f.parameters[2],  fits),
        γ  = map((f)->f.parameters[3],  fits),

        δα = map((f)->f.uncertainty[1], fits),
        δβ = map((f)->f.uncertainty[2], fits),
        δγ = map((f)->f.uncertainty[3], fits),

        cdf = map((f)->f.cumulative, fits),
    )

end

function fit(model)
    estimate = optimize(model.loss, model.constraint, model.Θ₀, IPNewton())
    return Optim.minimizer(estimate)
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
using Documenter, Luna

makedocs(sitename = "Luna")

max_stars_repo_name: basp/luna

max_issues_repo_name: basp/luna

max_forks_repo_name: basp/luna


text:
@testset "789.escape-the-ghosts.jl" begin
    @test escape_ghosts([[1,0],[0,3]], [0, 1]) == true
    @test escape_ghosts([[1,0]], [2,0]) == false
    @test escape_ghosts([[2,0]], [1,0]) == false
    @test escape_ghosts([[5,0],[-10,-2],[0,-5],[-2,-2],[-7,1]], [7, 7]) == false
    @test escape_ghosts([[-1,0],[0,1],[-1,0],[0,1],[-1,0]], [0,0])
end


text:
"""
Load data from file using metadata
"""
column_loader(t::Type{T}, io, metadata) where {T} = begin
    buffer = Vector{UInt8}(undef, metadata.len)
    column_loader!(buffer, t, io, metadata)
end

# load bytes bytes from io decompress into type
column_loader!(buffer, ::Type{T}, io, metadata) where {T} = begin
    readbytes!(io, buffer, metadata.len)
    return Blosc.decompress(T, buffer)
end


text:
import TensorOperations.optimaltree

export netcon_contraction!, netcon

# *************************************************************************** #
#                              Netcon functions
# *************************************************************************** #

"""
    function netcon(network::TensorNetworkCircuit)

This function will call the netcon implementation in tensor operations.
"""
function netcon(network::TensorNetworkCircuit, return_type::Symbol=:plan)

    # Generate the arguments for the netcon method from the given network.
    tensor_indices, index_dims, tensor_labels = create_netcon_input(network)

    # Run the netcon method to find the optimal contraction plan.
    contraction_tree, cost = optimaltree(tensor_indices, index_dims)

    if return_type == :tree
        # replace integers in contraction tree with node labels before returning
        contraction_tree = replace_leaf_integers_with_labels(contraction_tree,
                                                             tensor_labels)
        return contraction_tree

    elseif return_type == :plan
        # Convert the tree to an array of tensor pairs to contract in sequence.
        last_node_label = network.counters["node"]
        a, b, contraction_plan = convert_tree_to_plan(contraction_tree,
                                                      last_node_label,
                                                      tensor_labels)
        return contraction_plan

    else
        # If neither return option is picked, return to the user the raw netcon
        # output and labels of the tensors separately.
        return contraction_tree, tensor_labels
    end
end

"""
    function create_netcon_input(network::TensorNetworkCircuit)

This function perpares arguments for netcon to find the optimal contraction
plan of the given tensor network.
"""
function create_netcon_input(network::TensorNetworkCircuit)

    # Collect the indices and corresponding dimensions of the network
    # into arrays.
    num_tensors = length(network.nodes)

    tensor_indices = Array{Array{Symbol, 1}, 1}(undef, num_tensors)
    tensor_dims = Array{Array{Int, 1}, 1}(undef, num_tensors)
    tensor_labels = Array{Symbol, 1}(undef, num_tensors)

    for (i, node) in enumerate(values(network.nodes))
        tensor_indices[i] = node.indices
        tensor_dims[i] = node.dims
        tensor_labels[i] = node.data_label
    end

    # Make a dictionary of index dimensions for netcon.
    labels = reduce(vcat, tensor_indices); dims = reduce(vcat, tensor_dims)
    index_dims = Dict{Symbol, Int}(labels .=> dims)

    # The netcon implementation only needs tensor_indices and index_dims as
    # input but it returns a nested array of integers which need to be replaced
    # by node labels for contraction. To this end, tensor_labels is also
    # returned.
    return tensor_indices, index_dims, tensor_labels
end

"""
    function convert_tree_to_plan(tree::Union{Array{<:Any, 1}, Int},
                                  last_node_label::Int,
                                  tensor_labels::Array{Symbol, 1},
                                  plan::Array{Array{Symbol, 1}, 1})

This function converts a contraction tree, found by the netcon method, to a
contraction plan (array of node-label pairs) for PicoQuant contract functions.

tree - The contraction tree found by netcon.

last_node_label - Initially, this should be the number used to create the last
                  node added to the network being contracted, it is subsequantly
                  incremented to account for nodes created during contraction.

tensor_labels - An array of all the node labels for the nodes in a tensor
                network.

converted_plan - Initially, this argument should be an empty array. It will hold
                 the sequence of pairs when the function completes.
"""
function convert_tree_to_plan(tree::Union{Array{<:Any, 1}, Int},
                              last_node_label::Int,
                              tensor_labels::Array{Symbol, 1},
                              plan::Array{Array{Symbol, 1}, 1}
                                    =Array{Array{Symbol, 1}, 1}())

    # Recursively convert the plan into a sequence of pairs.
    if typeof(tree) == Int
        # Base case: If the plan is an integer corresponding to a single tensor,
        # return that integer as the label for this tensor and don't change the
        # number of nodes or add anything to the converted plan.
        return tensor_labels[tree], last_node_label, plan

    else
        # First, get the label for the tensor created by contracting the tree
        # in plan[1] and append its contraction plan to converted plan. Then
        # do the same for plan[2]. Finally, append the contraction of these
        # trees to the converted plan.
        a, last_node_label = convert_tree_to_plan(tree[1], last_node_label,
                                                  tensor_labels, plan)
        b, last_node_label = convert_tree_to_plan(tree[2], last_node_label,
                                                  tensor_labels, plan)
        append!(plan, [[a, b]])
        last_node_label = last_node_label + 1
        return Symbol("node_$last_node_label"), last_node_label, plan
    end
end

"""
    function replace_leaf_integers_with_labels(tree, labels::Array{Symbol, 1})

Function to replace the integers in a contraction tree with the corresponding
node labels.
"""
function replace_leaf_integers_with_labels(tree, labels::Array{Symbol, 1})
    # Base Case: If the given sub tree is an integer, meaning it is a leaf
    # of the parent tree, then return the corresponding node label.
    if typeof(tree) == Int
        return labels[tree]

    else
        # Replace the leaves of each branch of the given sub tree with
        # node labels.
        branch1 = replace_leaf_integers_with_labels(tree[1], labels)
        branch2 = replace_leaf_integers_with_labels(tree[2], labels)
        return [branch1, branch2]
    end
end

"""
    function netcon_contraction!(network::TensorNetworkCircuit,
                            output_shape::Union{String, Array{Int, 1}}="")

Function to contraction a network according to the optimal contraction plan
found by the netcon method.
"""
function netcon_contraction!(network::TensorNetworkCircuit,
                             output_shape::Union{String, Array{Int, 1}}="")

    if length(network.nodes) > 36
        error("The number of nodes in the network exceed the maximum for netcon")
    end
    contraction_plan = netcon(network)
    contract_network!(network, contraction_plan, output_shape)
end


text:
mutable struct Network
    inputs::Vector{Node}
    outputs::Vector{Node}

    fully_developed::Bool
    num_development_steps::Int8

    Network() = new([], [], false, 0)
end

function Network(starting_node::Node)
    network = Network()

    apply_to_all(starting_node) do x
        if x.type == InputNode()
            push!(network.inputs, x)
        elseif x.type == OutputNode()
            push!(network.outputs, x)
        end
    end

    return network
end

function get_number_neurons(network)
    num_neurons = 0
    apply_to_all(_ -> num_neurons += 1, network)
    return num_neurons
end

function update_inputs_outputs!(network, starting_nodes)
    clear!(network.inputs)
    clear!(network.outputs)

    if length(starting_nodes) == 0 return end

    apply_to_all(starting_nodes) do x
        if x.type == InputNode()
            push!(network.inputs, x)
        elseif x.type == OutputNode()
            push!(network.outputs, x)
        end
    end
end

function remove_non_reachable_nodes!(network)
    apply_to_all(network) do node
        node.reachable_from_input = false
    end

    apply_to_all(network.inputs, only_to_outputs = true) do node
        node.reachable_from_input = true
    end

    apply_to_all(network) do node
        filter!(x -> x.reachable_from_input, node.in_inhibitory)
        filter!(x -> x.reachable_from_input, node.in_excitatory)
        filter!(x -> x.reachable_from_input, node.out_inhibitory)
        filter!(x -> x.reachable_from_input, node.out_excitatory)
    end

    filter!(x -> x.reachable_from_input, network.outputs)
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark